
%
% Included are 5 Matlab functions:
% 
% ac      - calculates autocovariance estimates
% ssaeig  - calculates the T-EOFs and eigenvalues
% pc      - calculates the PCs
% rc      - calculates the RCs
% ssa     - this just calls the other functions, 
%             and returns T-EOFs, eigenvalues, and RCs.
% 
% Note
% Although the RCs sum correctly to give the original series, the spectra,
% as computed by Maximum Entropy, don't sum to the MEM spectrum 
% of the original series.
% The source of this error is not well understood.
% Many people use the T-PCs as opposed to RCs for stack spectra, and
% in this case the sum of T-PCs' power equals that of
% the original time series.
% 

----------------------------------------------------------------------------

function  c=ac(x,k,bias)
%    Syntax c=ac(x,k);  c=ac(x,k,'biased'); 
%    Ac calculates the auto-covariance for series 
%  x out to k lags. The result is output in c, which 
%  has k+1 elements. The first element is the covariance at lag zero;
%  succeeding elements 2:k+1 are the covariances at lags 1 to k.
%  Bias can be either:   'biased'       (Yule-Walker)  
%                        'unbiased'     The default is 'unbiased'.
%  If you want the autocorrelation rather than the autocovariance, just
%  call c=ac(x/std(x), lag).
%
%  Written by Eric Breitenberger.    Version date 5/17/95
%  Please send comments and suggestions to eric@gi.alaska.edu   
 
[N,col]=size(x);
if col>N, x=x'; [N,col]=size(x); end
if k>=N, error('Hey! Too big a lag!'), end
if col>=2, error('Hey! Vectors only!'), end
x=x-mean(x);
c=zeros(1,k+1);

if nargin==2, 
  for i=1:k+1
    c(i)=sum(x(1:N-i+1).*x(i:N));
    c(i)=c(i)/(N-i+1);
  end
elseif nargin==3,
  for i=1:k+1
    c(i)=sum(x(1:N-i+1).*x(i:N));
    if strcmp(bias,'unbiased'), c(i)=c(i)./(N-i+1); 
    elseif strcmp(bias,'biased'), c(i)=c(i)./N;
    else, error('Bias incorrectly specified!'), end 
  end
end

------------------------------------------------------------------------------

  function [E,V]=ssaeig(x, M)
%  Syntax: [E,V]=ssaeig(x, M);
%  This function starts an SSA of series 'x', for embedding dimension 'M'.
%  Returns:    E - eigenfunction matrix in standard form
%                  (columns are the eigenvectors, or T-EOFs)
%              V - vector containing variances (unnormalized eigenvalues)
%  E and V are ordered from large to small.
%  See section 2 of Vautard, Yiou, and Ghil, Physica D 58, 95-126, 1992.
%
%  Written by Eric Breitenberger.    Version date 5/17/95
%  Please send comments and suggestions to eric@gi.alaska.edu   

[N,col]=size(x);                
if col>N, x=x'; [N,col]=size(x); end   % change x to column vector if needed
if M-1>=N, error('Hey! Too big a lag!'), end
if col>=2, error('Hey! Vectors only!'), end

c=ac(x, M-1);          % calculate autocovariance estimates
T=toeplitz(c);         % create Toeplitz matrix (trajectory matrix)
[E,L]=eig(T);          % calculate eigenvectors, values of T
V=sum(L);              % create eigenvalue vector
V=V(M:-1:1);   
[r,col]=size(E);
E=E(:,col:-1:1);       % create standard eigenvector matrix

-------------------------------------------------------------------------------

   function [A]=pc(x, E)
%        Syntax: [A]=pc(x, E); 
%  This function calculates the principal components of the series x
%  from the eigenfunction matrix E, which is output from ssaeig.m
%  Returns:      A - principal components matrix (N-M+1 x M)
%  See section 2.4 of Vautard, Yiou, and Ghil, Physica D 58, 95-126, 1992.
%
%  Written by Eric Breitenberger.    Version date 5/17/95
%  Please send comments and suggestions to eric@gi.alaska.edu   

[N,col]=size(x);
if min(N,col)>1, error('x must be a vector.'), end
if col>1, x=x'; N=col; end     % convert x to column if necessary.
[M,c]=size(E);                
if M~=c, error('E is improperly dimensioned'), end
A=zeros(N-M+1,M);
for i=1:N-M+1;                 
  w=x(i:i+M-1);          
  A(i,:)=w'*E;
end

-------------------------------------------------------------------------------

  function [R]=rc(A,E)
% Syntax: [R]=rc(A,E);
% This function calculates the 'reconstructed components' using the 
% eigenvectors (E, from ssaeig.m) and principal components (A, from pc.m).
% R is N x M, where M is the embedding dimension used in ssaeig.m.
%
% See section 2.5 of Vautard, Yiou, and Ghil, Physica D 58, 95-126, 1992.
%
%  Written by Eric Breitenberger.   Version date 5/18/95
%  Please send comments and suggestions to eric@gi.alaska.edu   

[M,c]=size(E);
[ra, ca]=size(A);
if M~=c, error('E is improperly dimensioned.'),end
if ca~=M, error('A is improperly dimensioned.'),end
N=ra+M-1;  % Assumes A has N-M+1 rows.

R=zeros(N,M);
Z=zeros(M-1,M);
A=[A' Z'];
A=A';

% Calculate RCs
for k=1:M
  R(:,k)=filter(E(:,k),M,A(:,k));
end

% Adjust first M-1 rows and last M-1 rows
for i=1:M-1
  R(i,:)=R(i,:)*(M/i);
  R(N-i+1,:)=R(N-i+1,:)*(M/i);
end

-------------------------------------------------------------------------------

  function [E,V,R]=ssa(x, M) 
%  Syntax: [E,V,R]=ssa(x, M);
%  This function performs an SSA of series 'x', for embedding dimension 'M'.
%  Returns:  E - eigenfunction (T-EOF) matrix in standard form
%            V - vector containing variances (unnormalized eigenvalues)
%            R - Matrix of reconstructed components
%
%  See Vautard, Yiou, and Ghil, Physica D 58, 95-126, 1992.
%
%  Written by Eric Breitenberger.     Version date 5/17/95
%  Please send comments and suggestions to eric@gi.alaska.edu   

[E,V]=ssaeig(x,M);
[A]=pc(x,E);
[R]=rc(A,E);

% That's all for now!


