modwt_wcor.m

时间序列分析中很用的源码,书的原名为时间序列分析的小波方法.

function [wcor, CI_wcor] = modwt_wcor(WJtX, WJtY, p)
% modwt_wcor -- Calculate the wavelet correlation of MODWT wavelet coefficients.
%
%****f* wmtsa.dwt/modwt_wcor
%
% NAME
%   modwt_wcor -- Calculate the wavelet correlation of two sets of MODWT
%                 wavelet coefficients.
%
% SYNOPSIS
%   [wcor, CI_wcor] = modwt_wcor(WJtX, WJtY, [p])
%
% INPUTS
%   WJtX         - NxJ array containing MODWT-computed wavelet coefficients
%                  for X dataset
%                  where N  = number of time intervals,
%                        J = number of levels.
%   WJtY         - NxJ matrix containing MODWT-computed wavelet coefficients
%                  for Y dataset.
%   p            - (optional) percentage point for confidence interval.
%                  default: 0.025 ==> 95% confidence interval
%
% OUTPUTS
%   wcor         - Jx1 vector containing wavelet correlations.
%   CI_wcor      - Jx2 array containing confidence interval of wcor,
%                  lower bound (column 1) and upper bound (column 2).
%
% SIDE EFFECTS
%
%
% DESCRIPTION
%
%
% EXAMPLE
%
%
% NOTES
%
%
% ALGORITHM
%
%
% REFERENCES
%   Whitcher, B., P. Guttorp and D. B. Percival (2000)
%      Wavelet Analysis of Covariance with Application to Atmospheric Time
%      Series, \emph{Journal of Geophysical Research}, \bold{105}, No. D11,
%      14,941-14,962.
%
% SEE ALSO
%   modwt
%

% AUTHOR
%   Charlie Cornish
%   Brandon Whitcher
%
% CREATION DATE
%   2003-04-23
%
% Credits:
%   Based on original function wave_cov.m by Brandon Whitcher.
%
% COPYRIGHT
%
%
% REVISION
%   $Revision: 612 $
%
%***
%

% $Id: modwt_wcor.m 612 2005-10-28 21:42:24Z ccornish $


%%
%% Compute wavelet correlation with approximate 95% confidence interval
%% -----------------------------------------------------------------------
%% Input: X  Matrix containing wavelet coefficients with appropriate 
%%           boundary condition
%%        Y  Matrix containing wavelet coefficients with appropriate 
%%           boundary condition
%%
%% Output: C  Matrix containing the wavelet correlation (column 1), lower 
%%            95% quantile for confidence interval, upper 95% quantile 
%%            for confidence interval
%%

  usage_str = ['Usage:  [wcor, CI_wcor] = ', mfilename, ...
               '(WJtX, WJtY, [p])'];

  %%  Check input arguments and set defaults.
  error(nargerr(mfilename, nargin, [2:3], nargout, [0:2], 1, usage_str, 'struct'));

 set_default('p', 0.025);

  [N, J] = size(WJtX);
  
  WJtXY = WJtX .* WJtY;
  
  SSX  = NaN * zeros(J, 1);
  SSY  = NaN * zeros(J, 1);
  SSXY = NaN * zeros(J, 1);
  
  for (j = 1:J)
    WJtXNaN = WJtX(~isnan(WJtX(:,j)),j)';
    SSX(j) = sum(WJtXNaN.^2);
    WJtYNaN = WJtY(~isnan(WJtY(:,j)),j)';
    SSY(j) = sum(WJtYNaN.^2);
    WJtXYNaN = WJtXY(~isnan(WJtXY(:,j)),j)';
    SSXY(j) = sum(WJtXYNaN);
  end
  NX = sum(~isnan(WJtX),1)';
  COR = (SSXY./NX) ./ (sqrt(SSX./NX) .* sqrt(SSY./NX));
  
  % Return as column vector
  wcor = COR(:);
  
  if (nargout > 1)
    NDWT = floor(N ./ 2.^(1:J))';
  
    CI_wcor = zeros(J, 2);
    CI_wcor = [tanh(atanh(COR) - norminv(1-p) ./ sqrt(NDWT-3)), ...
               tanh(atanh(COR) + norminv(1-p) ./ sqrt(NDWT-3))];
  end
  

return