octave_packages/tsa-4.2.4/selmo2.m

   1 function X = selmo2(y,Pmax);
   2 % SELMO2 - model order selection for univariate and multivariate
   3 %   autoregressive models
   4 %
   5 %  X = selmo(y,Pmax);
   6 %
   7 %  y    data series
   8 %  Pmax maximum model order
   9 %  X.A, X.B, X.C parameters of AR model
  10 %  X.OPT... various optimization criteria
  11 %
  12 % see also: SELMO, MVAR,
  13
  14 %       $Id: selmo2.m 5090 2008-06-05 08:12:04Z schloegl $
  15 %       Copyright (C) 2007 by Alois Schloegl <a.schloegl@ieee.org>
  16 %       This is part of the TSA-toolbox. See also
  17 %       http://hci.tugraz.at/schloegl/matlab/tsa/
  18 %       http://octave.sourceforge.net/
  19 %       http://biosig.sourceforge.net/
  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
  35 [M,N]=size(y);
  36 if M>N,
  37         y=y';
  38 end;
  39 [M,N]=size(y);
  40
  41 % Univariate AR
  42 [AutoCov, AutoCorr, ARPMX, E, NC] = invest0(y,Pmax);
  43 [FPE,AIC,BIC,SBC,MDL,CATcrit,PHI,optFPE,optAIC,optBIC,optSBC,optMDL,optCAT,optPHI,s,C] = selmo(E,NC);
  44
  45
  46 % AIC for MVAR model
  47 [AR,RC,PE] = mvar(y',Pmax);
  48 E2 = repmat(NaN,1,Pmax);
  49 for k = 0:Pmax,
  50         S = PE(:,k*M+(1:M));
  51         E2(k+1) = trace(S);
  52
  53         %%%% FIX ME %%%%% this does not seem right because it depends on the scaling of y
  54         AIC_MV(k+1) = 2*log(det(S))+2*k*M*M/N; % Ding et al. 2000 refering to Akaike 1974
  55 end;
  56
  57
  58 X.A = [eye(M),-AR];
  59 X.B = eye(M);
  60 X.C = S;
  61 X.datatype = 'MVAR';
  62 X.MV.AIC = AIC_MV;
  63 X.UV.FPE = FPE;
  64 X.UV.AIC = AIC;
  65 X.UV.BIC = BIC;
  66 X.UV.MDL = MDL;
  67 X.UV.CAT = CATcrit;
  68 X.UV.PHI = PHI;
  69
  70 X.OPT.FPE = optFPE;
  71 X.OPT.AIC = optAIC;
  72 X.OPT.BIC = optBIC;
  73 X.OPT.MDL = optMDL;
  74 X.OPT.CAT = optCAT;
  75 X.OPT.PHI = optPHI;
  76
  77 [tmp,X.OPT.MVAIC] = min(AIC_MV);