reconstruct time series from SSA

问题:

let us consider following code

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;

The fragment of code was used from this site:

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

When I am running this code 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);

What is problem in my code? Please help me.

回答1:

As the error message suggests, you are trying to access the 16th column of rca and rca only has 15 columns.

Running your code in Octave, I get the following values:

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

The code errors when k=15 and m1=1, which gives k-m1+2 = 16, hence the error message.

EDIT

Rather than modify the function, why not simply call it with your data?

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);

回答2:

Have a look at this working GNU Octave compatible Basic SSA implementation: ssa-octave.m. This is an adaptation of my Scilab code which works fine for me. Pay attention to the comments inside, they may contain some important notes.

If you're interested I can also provide an M-SSA (multivariate) implementation. It's pretty close to Basic SSA and differs only in time delay embedding and diagonal averaging phases.

I want to say that SSA looks to me like (and I believe it is) PCA applied to time delay embedding of a time series.