octave_packages/tsa-4.2.4/mvaar.m

   1 function [x,e,Kalman,Q2] = mvaar(y,p,UC,mode,Kalman)
   2 % Multivariate (Vector) adaptive AR estimation base on a multidimensional
   3 % Kalman filer algorithm. A standard VAR model (A0=I) is implemented. The
   4 % state vector is defined as X=(A1|A2...|Ap) and x=vec(X')
   5 %
   6 % [x,e,Kalman,Q2] = mvaar(y,p,UC,mode,Kalman)
   7 %
   8 % The standard MVAR model is defined as:
   9 %
  10 %               y(n)-A1(n)*y(n-1)-...-Ap(n)*y(n-p)=e(n)
  11 %
  12 %       The dimension of y(n) equals s
  13 %
  14 %       Input Parameters:
  15 %
  16 %               y                       Observed data or signal
  17 %               p                       prescribed maximum model order (default 1)
  18 %               UC                      update coefficient      (default 0.001)
  19 %               mode            update method of the process noise covariance matrix 0...4 ^
  20 %                                       correspond to S0...S4 (default 0)
  21 %
  22 %       Output Parameters
  23 %
  24 %               e                       prediction error of dimension s
  25 %               x                       state vector of dimension s*s*p
  26 %               Q2                      measurement noise covariance matrix of dimension s x s
  27 %
  28
  29 %       $Id: mvaar.m 5090 2008-06-05 08:12:04Z schloegl $
  30 %       Copyright (C) 2001-2002 Christian Kasess
  31 %       Copyright (C) 2003, 2008 Alois Schloegl
  32 %
  33 %    This program is free software: you can redistribute it and/or modify
  34 %    it under the terms of the GNU General Public License as published by
  35 %    the Free Software Foundation, either version 3 of the License, or
  36 %    (at your option) any later version.
  37 %
  38 %    This program is distributed in the hope that it will be useful,
  39 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  40 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  41 %    GNU General Public License for more details.
  42 %
  43 %    You should have received a copy of the GNU General Public License
  44 %    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  45
  46
  47 if nargin<4,
  48         mode=0;
  49 end;
  50 if nargin<3,
  51         UC=0.001
  52 end;
  53 if nargin<2,
  54         p=1;
  55 end
  56 if nargin<1,
  57         fprintf(2,'No arguments supplied\n');
  58         return
  59 end;
  60
  61 if ~any(mode==(0:4))
  62         fprintf(2,'Invalid mode (0...4)\n');
  63         return
  64 end;
  65
  66
  67 [M,LEN] = size(y');             %number of channels, total signal length
  68 L = M*M*p;
  69
  70 if LEN<(p+1),
  71         fprintf(2,'Not enough observed data supplied for given model order\n');
  72         return
  73 end
  74
  75 ye = zeros(size(y));    %prediction of y
  76
  77 if nargout>1,
  78         x=zeros(L,LEN);
  79 end;
  80 if nargout>3,
  81         Q2=zeros(M,M,LEN);
  82 end
  83
  84
  85
  86 if nargin<5,
  87         %Kalman Filter initialsiation (Kp (K predicted or a-priori) equals K(n+1,n) )
  88         Kalman=struct('F',eye(L),'H',zeros(M,L),'G',zeros(L,M),'x',zeros(L,1),'Kp',eye(L),'Q1',eye(L)*UC,'Q2',eye(M),'ye',zeros(M,1));
  89         end;
  90
  91 upd = eye(L)/L*UC;              %diagonal matrix containing UC
  92
  93
  94 if(mode==3)
  95         Block=kron(eye(M),ones(M*p));
  96
  97 elseif(mode==4)
  98         index=[];
  99         Block1=[];
 100         Block0=[];
 101         for i=1:M,
 102                 index=[index ((i-1)*M*p+i:M:i*M*p)];
 103                 mone=eye(M);
 104                 mone(i,i)=0;
 105                 mzero=eye(M)-mone;
 106                 Block1=Blkdiag(Block1,kron(eye(p),mone));
 107                 Block0=Blkdiag(Block0,kron(eye(p),mzero));
 108         end;
 109 end;
 110
 111
 112 for n = 2:LEN,
 113         if(n<=p)
 114                 Yr=[y(n-1:-1:1,:)' zeros(M,p-n+1)];     %vector of past observations
 115                 Yr=Yr(:)';
 116         else
 117                 Yr=y(n-1:-1:n-p,:)';                                            %vector of past observations
 118                 Yr=Yr(:)';
 119         end
 120
 121         %Update of measurement matrix
 122         Kalman.H=kron(eye(M),Yr);
 123
 124
 125         %calculate prediction error
 126         ye(n,:)=(Kalman.H*Kalman.x)';
 127         err=y(n,:)-ye(n,:);
 128
 129         if ~any(isnan(err(:))),
 130                 %update of Q2 using the prediction error of the previous step
 131                 Kalman.Q2=(1-UC)*Kalman.Q2+UC*err'*err;
 132
 133
 134                 KpH=Kalman.Kp*Kalman.H';
 135                 HKp=Kalman.H*Kalman.Kp;
 136
 137                 %Kalman gain
 138                 Kalman.G=KpH*inv(Kalman.H*KpH+Kalman.Q2);
 139
 140                 %calculation of the a-posteriori state error covariance matrix
 141                 %K=Kalman.Kp-Kalman.G*KpH'; Althouh PK is supposed to be symmetric, this operation makes the filter unstable
 142                 K=Kalman.Kp-Kalman.G*HKp;
 143
 144                 %mode==0 no update of Q1
 145                 %update of Q1 using the predicted state error cov matrix
 146                 if(mode==1)
 147                         Kalman.Q1=diag(diag(K)).*UC;
 148                 elseif(mode==2)
 149                         Kalman.Q1=upd*trace(K);
 150                 elseif(mode==3)
 151                         Kalman.Q1=diag(sum((Block*diag(diag(K)))'))/(p*M)*UC;
 152                 elseif(mode==4)
 153                         avg=trace(K(index,index))/(p*M)*UC;
 154                         Kalman.Q1=Block1*UC+Block0*avg;
 155                 end
 156
 157                 %a-priori state error covariance matrix for the next time step
 158                 Kalman.Kp=K+Kalman.Q1;
 159
 160                 %current estimation of state x
 161                 Kalman.x=Kalman.x+Kalman.G*(err)';
 162         end; % isnan>(err)
 163
 164         if nargout>1,
 165                 x(:,n) = Kalman.x;
 166         end;
 167         if nargout>3,
 168                 Q2(:,:,n)=Kalman.Q2;
 169         end;
 170 end;
 171
 172 e = y - ye;
 173 x = x';
 174