octave_packages/tsa-4.2.4/acorf.m

   1 function [AUTOCOV,stderr,lpq,qpval] = acorf(Z,N);
   2 %  Calculates autocorrelations for multiple data series.
   3 %  Missing values in Z (NaN) are considered.
   4 %  Also calculates Ljung-Box Q stats and p-values.
   5 %
   6 %    [AutoCorr,stderr,lpq,qpval] = acorf(Z,N);
   7 %  If mean should be removed use
   8 %    [AutoCorr,stderr,lpq,qpval] = acorf(detrend(Z',0)',N);
   9 %  If trend should be removed use
  10 %    [AutoCorr,stderr,lpq,qpval] = acorf(detrend(Z')',N);
  11 %
  12 % INPUT
  13 %  Z    is data series for which autocorrelations are required
  14 %       each in a row
  15 %  N    maximum lag
  16 %
  17 % OUTPUT
  18 %  AutoCorr nr x N matrix of autocorrelations
  19 %  stderr   nr x N matrix of (approx) std errors
  20 %  lpq      nr x M matrix of Ljung-Box Q stats
  21 %  qpval    nr x N matrix of p-values on Q stats
  22 %
  23 % All input and output parameters are organized in rows, one row
  24 % corresponds to one series
  25 %
  26 % REFERENCES:
  27 %  S. Haykin "Adaptive Filter Theory" 3ed. Prentice Hall, 1996.
  28 %  M.B. Priestley "Spectral Analysis and Time Series" Academic Press, 1981.
  29 %  W.S. Wei "Time Series Analysis" Addison Wesley, 1990.
  30 %  J.S. Bendat and A.G.Persol "Random Data: Analysis and Measurement procedures", Wiley, 1986.
  31
  32 % calculating lpq, stderr, qpval from
  33 % suggested by Philip Gray, University of New South Wales,
  34
  35 %       $Id: acorf.m 5090 2008-06-05 08:12:04Z schloegl $
  36 %       Copyright (C) 1998-2002,2008 by Alois Schloegl <a.schloegl@ieee.org>
  37 %
  38 %    This program is free software: you can redistribute it and/or modify
  39 %    it under the terms of the GNU General Public License as published by
  40 %    the Free Software Foundation, either version 3 of the License, or
  41 %    (at your option) any later version.
  42 %
  43 %    This program is distributed in the hope that it will be useful,
  44 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  45 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  46 %    GNU General Public License for more details.
  47 %
  48 %    You should have received a copy of the GNU General Public License
  49 %    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  50
  51
  52 [nr,nc] = size(Z);
  53 NC = sum(~isnan(Z),2); % missing values
  54
  55 AUTOCOV = acovf(Z,N);
  56 AUTOCOV = AUTOCOV(:,2:N+1) ./ AUTOCOV(:,ones(1,N));
  57
  58 if nargout > 1
  59         stderr = sqrt(1./NC)*ones(1,N);
  60         lpq = zeros(nr,N);
  61         qpval = zeros(nr,N);
  62
  63         cum=zeros(nr,1);
  64         for k=1:N,
  65                 cum = cum+AUTOCOV(:,k).*conj(AUTOCOV(:,k))./(NC-k);
  66                 lpq(:,k) = NC.*(NC+2).*cum;                 % Ljung box Q for k lags
  67                 %qpval(:,k) = 1 - chi2cdf(lpq(:,k),k);  % p-value of Q stat
  68                 qpval(:,k) = 1 - gammainc(lpq(:,k)/2,k/2);      % replace chi2cdf by gammainc
  69         end;
  70 end;