ar1sig.m

蒙托卡罗模拟奇异谱分析

 function [g,a,mu2,c0,QUEUE,W,K]=ar1sig(x,C,E,R,k)
% AR1SIG - estimate AR1 parameters for noise + signal model
% Syntax: [g,a,mu2,c0,Q,W,K]=ar1sig(x,C,E,R,k);
%
% Estimates the parameters for the noise component of a composite 
% null-hypothesis (signal + AR(1) noise), using the method of 
% Allen and Smith 1995. 
%
% Input:  x - the data vector.
%         C - the covariance matrix of x, from SSAEIG.
%         E - EOFs of x.
%         R - the matrix of RCs, from SSA.
%         k - vector containing the indices of EOFs considered to be
%             signal: e. g. k=[8 9]. If k=0, or k is not specified,
%             parameters for the pure-noise null-hypothesis are computed.
%
% Output: g   - estimate of the lag-one autocorrelation.
%         a   - estimate of the noise innovation variance.
%         mu2 - estimated square on the mean.
%         c0  - estimated lag-zero covariance of the noise. 
%         Q   - noise projection matrix.
%         W   - estimated lag-covariance matrix of the red-noise.
%         K   - identity matrix, but with K(i,i)=0 if i is in k.
%         
% This function uses the algorithm described by Allen and Smith 1995,
% except that 'fzero' is used rather than Newton-Raphson root-finding.
% Fzero can sometimes be picky - keep an eye on the displayed
% convergence numbers to see if they make sense.
%
% Written by Eric Breitenberger.      Version 2/12/96
% Please send comments and suggestions to [e-mail]       
%

global QUEUE D_ZERO D_ONE NPOINTS

if nargin==3, k=0; end

x=x(:);
N=length(x);
NPOINTS=N;
m=mean(x);
x=x-m;

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

% Create signal and noise matrices:
[M,c]=size(E);
K=eye(M);
if k(1)~=0
  for i=1:length(k)
    K(k(i),k(i))=0;
  end
end

QUEUE=E*K*E';
QC=QUEUE*C;
% Sum the main diag. and first off-diag. 
D_ZERO=sum(diag(QC));
D_ONE=sum(diag(QC,1));

% Get initial estimates for gamma:
[g0,a0,c0]=ar1signv(x,sig);
disp(['Initial estimates: g=' num2str(g0) ', a=' num2str(a0)])

% Find g by getting zero of 'gammaest':
g=fzero('gammest2', g0, .0000001);

% Recompute everything:
gk=1:N-1;
gk=g.^gk;
mu2=(1/N)+(2/N^2)*sum((N-1:-1:1).*gk);
W=toeplitz(gk(1:M))/g-mu2;
QWQ=QUEUE*W*QUEUE;
c0=D_ZERO/sum(diag(QWQ));
a=sqrt((1-g^2)*c0);
disp(['  Final estimates: g=' num2str(g) ', a=' num2str(a)])