mcssadem.m

蒙托卡罗模拟奇异谱分析

% MCSSADEM - Demonstration of the MCSSA toolkit.
% Reproduces some of the figures in Allen and Smith 1996.
% Note that the reproductions may not be exact. This is because
% fewer Monte Carlo simulations are used in each example here,
% because Allen and Smith's AR(1) parameter estimates are slightly
% different in some cases than the ones used here (current versions
% of their code and mine give nearly identical results), and also 
% because in some cases I estimate AR(1) parameters where the known
% values are used in the paper.
%
% I suggest that you start the demo and go get a beverage - it will 
% take a while to run. 

echo off
N=100;             % this will be the number of Monte-Carlo simulations.
load mratest.dat    % load the test data.
x=mratest(:,1);     % x will be used as the example time series throughout.

%---------------------------------------------------------------
% Recreate Figure 1: the test signal alone and with noise added.
%---------------------------------------------------------------

figure
plot(x,'r')
hold on
plot(mratest(:,3),'y')
title('Figure 1: Test signal alone, and with added noise.')
drawnow

%---------------------------------------------------------------
% Recreate Figure 3:
% This differs slightly from Allen and Smith's because I use the
% estimated AR(1) parameters rather than the known parameters.
%---------------------------------------------------------------

[E,V,A,R,C]=ssa(x,40);     % do the SSA of x

[c]=mcssa1(x,E,N);         % generate the confidence limits

figure
semilogy(V,'+r')
hold on
semilogy(c,'c')
title('Figure 3: Eigenspectrum of test series and surrogate projections.')
drawnow

%disp('Press any key to continue'), pause

%---------------------------------------------------------------
% Recreate Figure 4:
% This differs slightly from Allen and Smith's because I use the
% estimated AR(1) parameters rather than the known parameters.
%---------------------------------------------------------------

[mE,fE]=eoffreq(E);  % calculate the dominant frequencies of the EOFs

figure
cm=mean(c')';
semilogy(fE,V,'+r')
hold on
errorbar(fE,cm,c(:,2)-cm,cm-c(:,1),'.y')
axis([0 .5 .05 15]);
title('Figure 4: ')
drawnow

%---------------------------------------------------------------
% Recreate Figure 5a:
% This differs slightly from Allen and Smith's because I use the
% estimated AR(1) parameters rather than the known parameters.
%---------------------------------------------------------------

[E,V,A,R,C]=ssa(x,30);   % do the SSA of x

[mE,fE]=eoffreq(E);      % calculate the dominant frequencies of the EOFs

[c]=mcssa1(x,E,N);       % generate the confidence limits

figure
cm=mean(c')';
semilogy(fE,V,'+r')
hold on
errorbar(fE,cm,c(:,2)-cm,cm-c(:,1),'.y')
axis([0 .5 .05 15]);
title('Figure 5a: M=30')
drawnow

%---------------------------------------------------------------
% Recreate Figure 5b:
% This differs slightly from Allen and Smith's because I use the
% estimated AR(1) parameters rather than the known parameters.
%---------------------------------------------------------------
[E,V,A,R,C]=ssa(x,60);    % do the SSA of x

[mE,fE]=eoffreq(E);       % calculate the dominant frequencies of the EOFs

[c]=mcssa1(x,E,N);        % generate the confidence limits

figure
cm=mean(c')';
semilogy(fE,V,'+r')
hold on
errorbar(fE,cm,c(:,2)-cm,cm-c(:,1),'.y')
axis([0 .5 .05 15]);
title('Figure 5b: M=60')
drawnow

%---------------------------------------------------------------
% Recreate Figure 6
% This differs slightly from Allen and Smith's because I use the
% estimated AR(1) parameters rather than the known parameters.
%---------------------------------------------------------------

[E,V,A,R,C]=ssa(x,40);   % do the SSA of x

[c]=mcssa1(x,40,N);      % generate confidence limits 

figure
semilogy(V,'+r')
hold on
semilogy(c,'c')
hold off
axis([1 40 .03 11]);
title('Figure 6: Eigenspectrum shape test.')
drawnow

%---------------------------------------------------------------
% Recreate Figure 7:
% This differs slightly from Allen and Smith's because my
% estimated AR(1) parameters are slightly different than those
% used in the preprint.
%---------------------------------------------------------------

[g,a,mu2]=ar1(x);        % get the AR(1) parameters

Er=ar1eof(g,a,mu2,40);   % generate the red-noise basis

Lr=eofproj(C,Er);        % project the covariance matrix on the red-noise eigenbasis 

[c,c2]=mcssa1(x,Er,N);   % generate confidence limits

[mEr,fEr]=eoffreq(Er);   % calculate the dominant frequencies for the red-noise EOFs.

cm=mean(c')';
figure
semilogy(fEr,Lr,'+r')
hold on
errorbar(fEr,cm,c(:,2)-cm,cm-c(:,1),'.y')
hold off
axis([0 .5 .04 10]);
title('Figure 7: Test against AR(1) eigenbasis.')
drawnow

%---------------------------------------------------------------
% Recreate Figure 8:
% This differs slightly from Allen and Smith's because my
% estimated AR(1) parameters are slightly different than those
% used in the preprint.
%---------------------------------------------------------------

[Er]=ar1eof2(x,C,E,R,[8 9],V);  % Compute the red-noise + signal eigenbasis

[mEr,fEr]=eoffreq(Er);  % compute the dominant frequencies of the EOFs

Lr=eofproj(C,Er);       % project the covariance matrix on the red-noise eigenbasis 

[c,c2]=mcssa2(x,C,Er,R,[8 9],N); % generate the confidence limits

cm=mean(c')';
figure
semilogy(fEr,Lr,'+r')
hold on
errorbar(fEr,cm,c(:,2)-cm,cm-c(:,1),'.y')
hold off
axis([0 .5 .04 10]);
title('Figure 8: Test against AR(1) + signal eigenbasis.')
drawnow