从SSA重构时间序列

Question

让我们考虑一下下面的代码

clear all;
  B=xlsread('data_generations1','A1','g8:g301');
  n=length(B);
  L =input('Give the size of the interval: ' );% Number of columns in the Data matrix
   m=n-L+1;%number of rows in the Data matrix
 X = zeros(m,L);
  for i=1:m
        X(i,:)=B(i:i+L-1);
  end;
   S=X*X'; 
    [U,autoval]=eig(S);
    [d,i]=sort(-diag(autoval));  
   d=-d;
   U=U(:,i);sev=sum(d); 
    plot((d./sev)*100),hold on,plot((d./sev)*100,'rx');
    title('Singular Spectrum');xlabel('Eigenvalue Number');ylabel('Eigenvalue (% Norm of trajectory matrix retained)')
   V=(X')*U; 
   rc=U*V';
% Step 3: Grouping
   I=input('Choose the agrupation of components to reconstruct the series in the form I=[i1,i2:ik,...,iL]  ')
   Vt=V';
   rca=U(:,I)*Vt(I,:);
% Step 4: Reconstruction

   y=zeros(n,1);  
   Lp=min(L,m);
   Kp=max(L,m);

   for k=0:Lp-2
     for m1=1:k+1;
      y(k+1)=y(k+1)+(1/(k+1))*rca(m1,k-m1+2);
     end
   end

   for k=Lp-1:Kp-1
     for m1=1:Lp;
      y(k+1)=y(k+1)+(1/(Lp))*rca(m1,k-m1+2);
     end
   end

   for k=Kp:n
      for m1=k-Kp+2:n-Kp+1;
       y(k+1)=y(k+1)+(1/(n-k))*rca(m1,k-m1+2);
     end
   end

   figure;subplot(2,1,1);hold on;xlabel('Data poit');ylabel('Original and reconstructed series')
   plot(x1);grid on;plot(y,'r')
   r=x1-y;
   subplot(2,1,2);plot(r,'g');xlabel('Data poit');ylabel('Residual series');grid on
   vr=(sum(d(I))/sev)*100;

代码片段来自此站点：

http://www.mathworks.com/matlabcentral/fileexchange/8115-singular-spectrum-analysis-smoother/content/ssa.m

当我运行这段代码时，averaging

Give the size of the interval: 15
Choose the agrupation of components to reconstruct the series in the form I=[i1,i2:ik,...,iL]  4

I =

     4

Attempted to access rca(1,16); index out of bounds because size(rca)=[280,15].

Error in averaging (line 37)
      y(k+1)=y(k+1)+(1/(Lp))*rca(m1,k-m1+2);

我的代码中有什么问题？请帮帮我。

Answer 1

如错误消息所示，您正在尝试访问rca的第16列，而rca只有15列。

在Octave中运行您的代码，我得到以下值：

>> I
I =  4
>> Kp
Kp =  280
>> Lp
Lp =  15
>> n
n =  294
>> m
m =  280

当k=15和m1=1提供k-m1+2 = 16时，代码出错，因此出现错误消息。

编辑

与其修改函数，为什么不简单地用数据调用它呢？

clear all;
B=xlsread('data_generations1','A1','g8:g301');
L =input('Give the size of the interval: ' );% Number of columns in the Data matrix
[y,r,vr]=ssa(B,L);

Answer 2

看一下这个工作的GNU Octave兼容的Basic SSA实现：ssa-octave.m。这是我的Scilab代码的改编，它对我来说工作得很好。注意里面的注释，它们可能包含一些重要的注释。

如果您感兴趣，我还可以提供一个M-SSA (多变量)实现。它非常接近基本的SSA，只是在时间延迟嵌入和对角平均相位上有所不同。

我想说，在我看来，SSA就像(我相信是) PCA应用于时间序列的时延嵌入。