mcssa2.m

蒙托卡罗模拟奇异谱分析

  function [c,c2,L,s,vs,g,a,mu2]=mcssa2(x,C,E,R,k,N,p,method)
% MXCSSA2 - MCSSA for red noise + signal hypothesis. 
% Syntax: [c,c2,L,s,vs,g,a,mu2]=mcssa2(x,C,E,R,k,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 a deterministic signal (composed of the k RCs), contaminated 
% by an AR(1) process with unknown mean.
%
% Parameters for the AR(1) process are estimated from the data, with the 
% signal removed, and then used to generate realizations of AR(1) noise.
% After adding the signal to these surrogate realizations, the covariance 
% matrix is computed. The eigenspectra are then computed by projection 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.
%
%         C - the covariance matrix of x, from SSAEIG.
%
%         E - An eigenbasis, for example the SSA eigenvectors of x (from
%             SSAEIG) or the red-noise + signal basis (from AR1EOF2).
%
%             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).
%
%         R - the matrix of RCs (from SSA).
%
%         k - vector containing the indices of EOFs considered to be signal:
%             e. g. k=[1 8 9]. If k=0, parameters for the pure-noise 
%             null-hypothesis are computed.
%
%         N - (optional) number of surrogate realizations. The bigger N is,
%             the better, but this can be very computationally expensive. 
%             The default is N=1000.
%
%         p - (optional) 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 - A matrix containing the upper and lower confidence 
%             limits; 2*length(p) by M.
%
%         c2- A matrix with length(p) columns, 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 squared Frobenius norm of 
%             each red-noise covariance matrix. 
%
%         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 3/17/97
% Please send comments and suggestions to [e-mail]       
%

if nargin==5, N=1000; p=.95; method='unbiased';
elseif nargin==6
  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 fifth argument.')
  end
elseif nargin==7
  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);

% Signal model : sum of RCs
if k(1)~=0
  sig=R(:,k)*ones(length(k),1);
else
  sig=zeros(n,1);
end

[g,a,mu2]=ar1sig(x,C,ssa(x,M,method),R,k); % AR(1) parameters using Allen and Smith GLS method 

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
  X=X+sig*ones(1,Ni(k));     % add the signal to the surrogates
  X=X-ones(n,1)*mean(X);     % center the surrogates

  % 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 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