SSA data pre-processing for rainfall prediction

function [y,r,vr]=SCALE_DATA(p,upper)

% p Original time series (column vector form)
% upper  Window length
% y  Reconstructed time series
% r  Residual time series r=p-y
% vr Relative value of the norm of the approximated trajectory matrix with respect
%	  to the original trajectory matrix

% The program output is the Singular Spectrum of p (must be a column vector),
% using a window length upper. You must choose the components be used to reconstruct 
%the series in the form [i1,i2:ik,...,iupper], based on the Singular Spectrum appearance.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Step1 : Build trajectory matrix

   N=length(p); 
   if upper>N/2;upper=N-upper;end
	K=N-upper+1; 
   X=zeros(upper,K);  
	for i=1:K
	  X(1:upper,i)=p(i:upper+i-1); 
	end
    
% Step 2: SVD

   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,...,iupper]  ')
   Vt=V';
   rca=U(:,I)*Vt(I,:);

% Step 4: Reconstruction

   y=zeros(N,1);  
   upperp=min(upper,K);
   Kp=max(upper,K);

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

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

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

   figure;subplot(2,1,1);hold on;xlabel('Data poit');ylabel('Original and reconstructed series')
   plot(p);grid on;plot(y,'r')

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