octave_packages/tsa-4.2.4/arfit2.m

   1 function [w, MAR, C, sbc, fpe, th] = arfit2(Y, pmin, pmax, selector, no_const)
   2 % ARFIT2 estimates multivariate autoregressive parameters
   3 % of the MVAR process Y
   4 %
   5 %   Y(t,:)' = w' + A1*Y(t-1,:)' + ... + Ap*Y(t-p,:)' + x(t,:)'
   6 %
   7 % ARFIT2 uses the Nutall-Strand method (multivariate Burg algorithm)
   8 % which provides better estimates the ARFIT [1], and uses the
   9 % same arguments. Moreover, ARFIT2 is faster and can deal with
  10 % missing values encoded as NaNs.
  11 %
  12 % [w, A, C, sbc, fpe] = arfit2(v, pmin, pmax, selector, no_const)
  13 %
  14 % INPUT:
  15 %  v            data - each channel in a column
  16 %  pmin, pmax   minimum and maximum model order
  17 %  selector     'fpe' or 'sbc' [default]
  18 %  no_const     'zero' indicates no bias/offset need to be estimated
  19 %               in this case is w = [0, 0, ..., 0]';
  20 %
  21 % OUTPUT:
  22 %  w            mean of innovation noise
  23 %  A            [A1,A2,...,Ap] MVAR estimates
  24 %  C            covariance matrix of innovation noise
  25 %  sbc, fpe     criteria for model order selection
  26 %
  27 % see also: ARFIT, MVAR
  28 %
  29 % REFERENCES:
  30 %  [1] A. Schloegl, 2006, Comparison of Multivariate Autoregressive Estimators.
  31 %       Signal processing, p. 2426-9.
  32 %  [2] T. Schneider and A. Neumaier, 2001.
  33 %       Algorithm 808: ARFIT-a Matlab package for the estimation of parameters and eigenmodes
  34 %       of multivariate autoregressive models. ACM-Transactions on Mathematical Software. 27, (Mar.), 58-65.
  35
  36 %       $Id: arfit2.m 9609 2012-02-10 10:18:00Z schloegl $
  37 %       Copyright (C) 1996-2005,2008,2012 by Alois Schloegl <alois.schloegl@ist.ac.at>
  38 %       This is part of the TSA-toolbox. See also
  39 %       http://pub.ist.ac.at/~schloegl/matlab/tsa/
  40 %       http://octave.sourceforge.net/
  41 %       http://biosig.sourceforge.net/
  42
  43 %    This program is free software: you can redistribute it and/or modify
  44 %    it under the terms of the GNU General Public License as published by
  45 %    the Free Software Foundation, either version 3 of the License, or
  46 %    (at your option) any later version.
  47 %
  48 %    This program is distributed in the hope that it will be useful,
  49 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  50 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  51 %    GNU General Public License for more details.
  52 %
  53 %    You should have received a copy of the GNU General Public License
  54 %    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  55
  56
  57 %%%%% checking of the input arguments was done the same way as ARFIT
  58 if (pmin ~= round(pmin) || pmax ~= round(pmax))
  59         error('Order must be integer.');
  60 end
  61 if (pmax < pmin)
  62         error('PMAX must be greater than or equal to PMIN.')
  63 end
  64
  65 % set defaults and check for optional arguments
  66 if (nargin == 3)                % no optional arguments => set default values
  67         mcor       = 1;         % fit intercept vector
  68         selector   = 'sbc';     % use SBC as order selection criterion
  69 elseif (nargin == 4)            % one optional argument
  70         if strcmp(selector, 'zero')
  71                 mcor     = 0;   % no intercept vector to be fitted
  72                 selector = 'sbc';       % default order selection
  73         else
  74                 mcor     = 1;   % fit intercept vector
  75         end
  76 elseif (nargin == 5)            % two optional arguments
  77         if strcmp(no_const, 'zero')
  78                 mcor     = 0;           % no intercept vector to be fitted
  79         else
  80                 error(['Bad argument. Usage: ', '[w,A,C,SBC,FPE,th]=AR(v,pmin,pmax,SELECTOR,''zero'')'])
  81         end
  82 end
  83
  84
  85
  86 %%%%% Implementation of the MVAR estimation
  87 [N,M]=size(Y);
  88
  89 if mcor,
  90         [m,N] = sumskipnan(Y,1);                    % calculate mean
  91         m = m./N;
  92         Y = Y - repmat(m,size(Y)./size(m));    % remove mean
  93 end;
  94
  95 [MAR,RCF,PE] = mvar(Y, pmax, 2);   % estimate MVAR(pmax) model
  96
  97 N = min(N);
  98
  99 %if 1;nargout>3;
 100 ne = N-mcor-(pmin:pmax);
 101 for p=pmin:pmax,
 102         % Get downdated logarithm of determinant
 103         logdp(p-pmin+1) = log(det(PE(:,p*M+(1:M))*(N-p-mcor)));
 104 end;
 105
 106 % Schwarz's Bayesian Criterion
 107 sbc = logdp/M - log(ne) .* (1-(M*(pmin:pmax)+mcor)./ne);
 108
 109 % logarithm of Akaike's Final Prediction Error
 110 fpe = logdp/M - log(ne.*(ne-M*(pmin:pmax)-mcor)./(ne+M*(pmin:pmax)+mcor));
 111
 112 % Modified Schwarz criterion (MSC):
 113 % msc(i) = logdp(i)/m - (log(ne) - 2.5) * (1 - 2.5*np(i)/(ne-np(i)));
 114
 115 % get index iopt of order that minimizes the order selection
 116 % criterion specified by the variable selector
 117 if strcmpi(selector,'fpe');
 118     [val, iopt]  = min(fpe);
 119 else %if strcmpi(selector,'sbc');
 120     [val, iopt]  = min(sbc);
 121 end;
 122
 123 % select order of model
 124 popt = pmin + iopt-1; % estimated optimum order
 125
 126 if popt<pmax,
 127         [MAR, RCF, PE] = mvar(Y, popt, 2);
 128 end;
 129 %end
 130
 131 C = PE(:,size(PE,2)+(1-M:0));
 132
 133 if mcor,
 134         I = eye(M);
 135         for k = 1:popt,
 136                 I = I - MAR(:,k*M+(1-M:0));
 137         end;
 138         w = -I*m';
 139 else
 140         w = zeros(M,1);
 141 end;
 142
 143 if nargout>5,   th=[];  fprintf(2,'Warning ARFIT2: output TH not defined\n'); end;