|
2414
|
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 |
|
4907
|
35 % $Id$ |
|
12580
|
36 % Copyright (C) 1998-2002,2008 by Alois Schloegl <alois.schloegl@gmail.com> |
|
4907
|
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. |
|
2414
|
42 % |
|
4907
|
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 |
|
2414
|
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; |
