mcssa1.m

蒙托卡罗模拟奇异谱分析

  function [c,c2,L,s,vs,g,a,mu2]=mcssa1(x,E,N,p,method)
% MCSSA1 - Monte Carlo SSA test against pure red noise.
% Syntax: [c,c2,L,s,vs,g,a,mu2]=mcssa1(x,E,N,p,method);
%
% Generates confidence limits for an SSA eigenspectrum, given the time 
% series x and eigenvectors E. The null hypothesis tested is that the data
% have been generated by an AR(1) process with unknown mean.
%
% Parameters for an AR(1) process are estimated from the data, and used
% to generate realizations of AR(1) noise. From these surrogate realizations,
% covariance matrices are computed. Eigenspectra are then computed by 
% projection of the covariance matrices on the basis E. These surrogate 
% eigenspectra are then used to construct the confidence interval.
%
%  See Allen 1992, Allen and Smith 1994, 1996.
%
% Input:  x - the time series.
%
%         E - a fixed M by M eigenbasis, for example the SSA eigenvectors 
%             of x (from SSAEIG), or a red-noise basis (from AR1EOF).
%
%             Alternatively, if E is an integer (the embedding dimension M),
%             the SSA eigenspectrum for each surrogate will be computed (the
%             test proposed by Elsner and Tsonis, also by Dettinger et al.,
%             which Allen and Smith call the FMU test).
%
%         N - number of surrogate realizations (default 1000).
%             For a good test, the bigger N the better, but
%             this can be very computationally expensive.
%
%         p - confidence limit, in a two-tailed form, e. g. for p=.95
%             the confidence limits will be at 2.5% and 97.5% of the
%             Monte Carlo distribution. The default is p=.95.
%             p may be a vector, e. g. p=[.9 .95 .99].
%
%    method - (optional) method of calculating the covariance matrix:
%                 'unbiased' (N-k weighted) (default)
%                 'biased'   (N-weighted or Yule-Walker)  
%                 'BK'       (Broomhead/King type estimate)
%
% Output: c - An M by 2*length(p) matrix containing upper and 
%             lower confidence limits.
%
%        c2 - A matrix where the i-th row contains the probabilities
%             of realizing i 'significant' eigenvalues in a single 
%             surrogate eigenspectrum. 
%
%         L - The M by N matrix of surrogate eigenspectra.
%
%         s - A vector containing the variance (N-weighted) of each surrogate. 
%
%        vs - A vector containing the squared Frobenius norm of each  
%             surrogate eigenvalue matrix. 
%
%         g - estimate of the lag-one autocorrelation of the AR(1) process.
%
%         a - estimate of the noise variance of the AR(1) process.
%
%       mu2 - estimated square on the mean of the AR(1) process.
%
% Written by Eric Breitenberger.      Version 5/22/96
% Please send comments and suggestions to [e-mail]       
%

if nargin==2, N=1000; p=.95; method='unbiased';
elseif nargin==3
  if isstr(N), method=N; p=.95; N=1000;
  elseif N<=1, p=N; N=1000; method='unbiased';
  elseif  N>1, p=.95; method='unbiased';
  else error('Improper specification of third argument.')
  end
elseif nargin==4
  if N<=1, method=p;  p=N; N=1000;
  elseif isstr(p), method=p; p=.95;
  else, method='unbiased';
  end
end 

disp(['MCSSA: N=' num2str(N) ', p=' num2str(p) ', using ' method ' covariance estimation.'])

[M,r]=size(E);
if M~=r, error('E is not square.'), end
if M==1, M=E; end
n=length(x);

[g,a,mu2]=ar1(x);    % AR(1) parameters using Allen and Smith GLS method 
disp(['AR(1) estimates: g=' num2str(g) ', a=' num2str(a)])

L=zeros(M,N);
vs=zeros(1,N);
s=vs;

% Computational considerations:
% The matrix L is of size M by N - this is the principal limitation
% on memory usage. As L must be sorted, two copies of it are needed,
% consuming 2*M*N*8 bytes. To reduce looping, it is desireable to 
% compute many surrogate processes at once - storage of all of these
% would require N*length(x)*8 bytes, which can easily get very large.
% For this reason, the surrogates are computed and stored B at a time.
% B can be varied depending on the available memory. 

B=1000;
if N>B  % set up to do surrogates in K blocks of size B
  if rem(N,B)==0
    K=floor(N/B);
    Ni=B*ones(1,K);
  else
    K=floor(N/B)+1;
    Ni=B*ones(1,K);
    Ni(K)=rem(N,B);
  end
else
  Ni=N; 
  K=1;
  B=N;
end

for k=1:K
  X=ar1noise(n,Ni(k),g,a);   % generate the AR(1) surrogate data

  % Calculate the surrogate eigenspectra:
  if length(E)>1     % fixed basis
    for i=1:Ni(k)
      if strcmp(method,'BK')==0
        c=ac(X(:,i), M-1,method);     % calculate acv estimates
        T=toeplitz(c);                % create covariance matrix
      else
        T=bk(X(:,i), M);              % BK autocovariance matrix
      end
      V=E'*T*E;                       % project T on the data eigenbasis
      s(i+(k-1)*B)=X(:,i)'*X(:,i);    % calculate variance of each surrogate
      vs(i+(k-1)*B)=sum(sum(V.^2));   % Calculate the squared Frobenius norm of V
      L(:,i+(k-1)*B)=diag(V);
    end
  else % variable basis
    for i=1:Ni(k)
      s(i+(k-1)*B)=X(:,i)'*X(:,i);    % calculate variance of each surrogate
      [Essa,V]=ssaeig(X(:,i),M,method);
      L(:,i+(k-1)*B)=V';
    end
  end
end

clear X

% Finally, extract the confidence limits for each element of p:
n=length(p);
len=zeros(1,n);
c=zeros(M,2*n);

for i=1:n
  ind=2*(i-1);
  [c(:,ind+1:ind+2),c2]=conflim(L,p(i));  % Extract the confidence limits from L
  eval(['q' int2str(i) '=c2; len(i)=length(q' int2str(i) ');']);
end

c2=zeros(n,max(len));
for i=1:n
  eval(['c2(i,1:len(i))=q' int2str(i) ';']);
end