octave_packages/tsa-4.2.4/pacf.m

   1 function [PARCOR,sig,cil,ciu]= pacf(Z,KMAX);
   2 % Partial Autocorrelation function
   3 % [parcor,sig,cil,ciu] = pacf(Z,N);
   4 %
   5 % Input:
   6 %       Z    Signal, each row is analysed
   7 %       N    # of coefficients
   8
   9 % Output:
  10 %       parcor autocorrelation function
  11 %       sig     p-value for significance test
  12 %       cil     lower confidence interval
  13 %       ciu     upper confidence interval
  14 %
  15 % see also: DURLEV, LATTICE, AC2RC, AR2RC,
  16 %       FLAG_IMPLICIT_SIGNIFICANCE
  17
  18 %       $Id: pacf.m 5090 2008-06-05 08:12:04Z schloegl $
  19 %       Copyright (C) 1997-2002,2008 by Alois Schloegl <a.schloegl@ieee.org>
  20 %
  21 %    This program is free software: you can redistribute it and/or modify
  22 %    it under the terms of the GNU General Public License as published by
  23 %    the Free Software Foundation, either version 3 of the License, or
  24 %    (at your option) any later version.
  25 %
  26 %    This program is distributed in the hope that it will be useful,
  27 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  28 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  29 %    GNU General Public License for more details.
  30 %
  31 %    You should have received a copy of the GNU General Public License
  32 %    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  33
  34 [nr,nc] = size(Z);
  35 if nc<KMAX,
  36         warning('too less elements.\nmake sure the data is row order\n')
  37 end;
  38 [s,n] = sumskipnan(Z,2);
  39 Z = Z - repmat(s./n,1,nc);      % remove mean
  40
  41 if (nargin == 1), KMAX = N-1; end;
  42
  43 AutoCov = acovf(Z,KMAX);
  44 [AR,PARCOR,PE] = durlev(AutoCov); % PARCOR are the reflection coefficients
  45 %[AR,PARCOR,PE] = lattice(Z,KMAX); % PARCOR are the reflection coefficients
  46 PARCOR = -PARCOR;                       % the partial correlation coefficients are the negative reflection coefficient.
  47
  48 if nargout<2, return, end;
  49
  50
  51 % significance test
  52 s = 1./sqrt(repmat(n,1,KMAX)-1-ones(nr,1)*(1:KMAX));
  53 sig = normcdf(PARCOR,0,s);
  54 sig = min(sig,1-sig);
  55
  56
  57 if nargout<3, return, end;
  58 % calculate confidence interval
  59 if exist('flag_implicit_significance')==2;
  60         alpha = flag_implicit_significance;
  61 else
  62         alpha = 0.05;
  63 end;
  64
  65 fprintf(1,'PACF: confidence interval for alpha=%f\n', alpha);
  66 ciu = norminv(alpha/2).*s;
  67 cil = -ciu;
  68
  69