xc.m

多通道奇异谱分析程序

function  q=xc(x,y,k,method) 
%  Syntax: q=xc(x,y,k); q=xc(x,y,k,'BK');
%  XC calculates the cross-covariance between series 
%  x and y out to k lags. The result is output in q, which 
%  has 2k+1 elements. q(k+1) corresponds to zero lag.
% 
%  NOTE THAT X LEADS Y FOR POSITIVE LAGS.
%
% Optionally, different methods may be used to 
% estimate the cross-covariances:
%       'unbiased' - N-i weighting. (default)
%       'biased'   - N weighting.
%       'BK' - Broomhead/King type estimation.
%
%  Note that the BK method is of limited use unless you are computing
%  the full covariance matrix - see BK.M.                    
% 
%  Written by Eric Breitenberger.   Version date 1/11/96
%  Please send comments and suggestions to [e-mail]   
%
if nargin==3, method='unbiased'; end

[N,cx]=size(x);
[Ny,cy]=size(y);
if cx>N, tmp=cx; cx=N; N=tmp; end
if cy>Ny, tmp=cy; cy=Ny; Ny=tmp; end
if cx>1, error('X must be a vector.'), end
if cy>1, error('Y must be a vector.'), end
if N~=Ny, error('Lengths of input series must be the same'),end
if k>=N, error('Too big a lag!'), end

x=x(:);
y=y(:);
x=x-mean(x);
y=y-mean(y);
q=ones(1,2*k+1);

if strcmp(method, 'unbiased')
  q(k+1)=x'*y/N;          % zero lag correlation
  for i=1:k               % positive lag correlations
    q(k+1+i)= x(1:N-i)' * y(i+1:N) /(N-i);
  end
  for i=1:k               % negative lag correlations
    q(k+1-i)= y(1:N-i)' * x(i+1:N) /(N-i);
  end

elseif strcmp(method,'biased')
  q(k+1)=x'*y;          % zero lag correlation
  for i=1:k               % positive lag correlations
    q(k+1+i)= x(1:N-i)' * y(i+1:N);
  end
  for i=1:k               % negative lag correlations
    q(k+1-i)= y(1:N-i)' * x(i+1:N);
  end
  q=q/N;

elseif strcmp(method,'BK')
  ind=1:N-k;
  q(k+1)=x(ind)'*y(ind);    % zero lag correlation
  for i=1:k                 % positive lag correlations
    q(k+1+i)= x(ind)' * y(ind+i);
  end
  for i=1:k                 % negative lag correlations
    q(k+1-i)= y(ind)' * x(ind+i);
  end
  q=q/(N-k);

else
  error('Improper specification for method.')
end