octave_packages/tsa-4.2.4/amarma.m

   1 function [z,e,REV,ESU,V,Z,SPUR] = amarma(y, Mode, MOP, UC, z0, Z0, V0, W);
   2 % Adaptive Mean-AutoRegressive-Moving-Average model estimation
   3 % [z,E,ESU,REV,V,Z,SPUR] = amarma(y, mode, MOP, UC, z0, Z0, V0, W);
   4 % Estimates AAR parameters with Kalman filter algorithm
   5 %       y(t) = sum_i(a(i,t)*y(t-i)) + mu(t) + E(t)
   6 %
   7 % State space model:
   8 %       z(t)=G*z(t-1) + w(t)      w(t)=N(0,W)
   9 %       y(t)=H*z(t)   + v(t)      v(t)=N(0,V)
  10 %
  11 % G = I,
  12 % z = [µ(t)/(1-sum_i(a(i,t))),a_1(t-1),..,a_p(t-p),b_1(t-1),...,b_q(t-q)];
  13 % H = [1,y(t-1),..,y(t-p),e(t-1),...,e(t-q)];
  14 % W = E{(z(t)-G*z(t-1))*(z(t)-G*z(t-1))'}
  15 % V = E{(y(t)-H*z(t-1))*(y(t)-H*z(t-1))'}
  16 %
  17 % Input:
  18 %       y       Signal (AR-Process)
  19 %       Mode
  20 %           [0,0] uses V0 and W
  21 %
  22 %       MOP     Model order [m,p,q], default [0,10,0]
  23 %                       m=1 includes the mean term, m=0 does not.
  24 %                       p and q must be positive integers
  25 %                       it is recommended to set q=0.
  26 %       UC      Update Coefficient, default 0
  27 %       z0      Initial state vector
  28 %       Z0      Initial Covariance matrix
  29 %
  30 % Output:
  31 %       z       AR-Parameter
  32 %       E       error process (Adaptively filtered process)
  33 %       REV     relative error variance MSE/MSY
  34 %
  35 %
  36 % see also: AAR
  37 %
  38 % REFERENCE(S):
  39 % [1] A. Schloegl (2000), The electroencephalogram and the adaptive autoregressive model: theory and applications.
  40 %     ISBN 3-8265-7640-3 Shaker Verlag, Aachen, Germany.
  41 % [2] Schlögl A, Lee FY, Bischof H, Pfurtscheller G
  42 %     Characterization of Four-Class Motor Imagery EEG Data for the BCI-Competition 2005.
  43 %     Journal of neural engineering 2 (2005) 4, S. L14-L22
  44 %
  45 % More references can be found at
  46 %     http://www.dpmi.tu-graz.ac.at/~schloegl/publications/
  47
  48 %       $Id: amarma.m 5376 2008-10-13 15:53:47Z schloegl $
  49 %       Copyright (C) 1998-2002,2005,2006,2007,2008 by  Alois Schloegl <a.schloegl@ieee.org>
  50 %
  51 %    This program is free software: you can redistribute it and/or modify
  52 %    it under the terms of the GNU General Public License as published by
  53 %    the Free Software Foundation, either version 3 of the License, or
  54 %    (at your option) any later version.
  55 %
  56 %    This program is distributed in the hope that it will be useful,
  57 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  58 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  59 %    GNU General Public License for more details.
  60 %
  61 %    You should have received a copy of the GNU General Public License
  62 %    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  63
  64
  65 [nc,nr]=size(y);
  66
  67 if nargin<2 Mode=0;
  68 elseif isnan(Mode) return; end;
  69 if nargin<3, MOP=[0,10,0]; end;
  70 if nargin<8, W  = nan ; end;
  71 if length(MOP)==0,      m=0;p=10; q=0; MOP=p;
  72 elseif length(MOP)==1,  m=0;p=MOP(1); q=0; MOP=p;
  73 elseif length(MOP)==2,  fprintf(1,'Error AMARMA: MOP is ambiguos\n');
  74 elseif length(MOP)>2,   m=MOP(1); p=MOP(2); q=MOP(3);MOP=m+p+q;
  75 end;
  76
  77 if prod(size(Mode))>1
  78         aMode=Mode(1);
  79         eMode=Mode(2);
  80 end;
  81 %fprintf(1,['a' int2str(aMode) 'e' int2str(eMode) ' ']);
  82
  83
  84 e = zeros(nc,1);
  85 V = zeros(nc,1);V(1)=V0;
  86 T = zeros(nc,1);
  87 ESU = zeros(nc,1)+nan;
  88 SPUR = zeros(nc,1)+nan;
  89 z = z0(ones(nc,1),:);
  90
  91 dW = UC/MOP*eye(MOP);                % Schloegl
  92
  93 %------------------------------------------------
  94 %       First Iteration
  95 %------------------------------------------------
  96
  97 H = zeros(MOP,1);
  98 if m,
  99     %M0   = z0(1)/(1-sum(z0(2:p+1))); %transformierter Mittelwert
 100     H(1) = 1;%M0;
 101     %z0(1)= 1;
 102 end;
 103
 104 Z = Z0;
 105 zt= z0;
 106
 107 A1 = zeros(MOP); A2 = A1;
 108
 109 %------------------------------------------------
 110 %       Update Equations
 111 %------------------------------------------------
 112
 113 for t=1:nc,
 114         %H=[y(t-1); H(1:p-1); E ; H(p+1:MOP-1)]
 115
 116
 117         if t<=p, H(m+(1:t-1)) = y(t-1:-1:1);     %H(p)=mu0;          % Autoregressive
 118         else     H(m+(1:p)) = y(t-1:-1:t-p); %mu0];
 119         end;
 120
 121         if t<=q, H(m+p+(1:t-1)) = e(t-1:-1:1);       % Moving Average
 122         else     H(m+p+(1:q)) = e(t-1:-1:t-q);
 123         end;
 124
 125         % Prediction Error
 126         E = y(t) - zt*H;
 127
 128         e(t) = E;
 129
 130         if ~isnan(E),
 131                 E2 = E*E;
 132                 AY = Z*H;
 133
 134 %                [zt, t, y(t), E,ESU(t),V(t),H,Z],pause,
 135
 136                 ESU(t) = H'*AY;
 137
 138                 if eMode==0
 139                         V(t) = V0;
 140                 elseif eMode==1
 141                         V0 = V(t-1);
 142                         V(t) = V0*(1-UC)+UC*E2;
 143                 elseif eMode==2
 144                         V0 = 1;
 145                         V(t) = V0; %V(t-1)*(1-UC)+UC*E2;
 146                 elseif eMode==3
 147                         V0 = 1-UC;
 148                         V(t) = V0; %(t-1)*(1-UC)+UC*E2;
 149                 elseif eMode==4
 150                         V0 = V0*(1-UC)+UC*E2;
 151                         V(t) = V0;
 152                 elseif eMode==5
 153                         V(t)=V0;
 154                         %V0 = V0;
 155                 elseif eMode==6
 156                         if E2>ESU(t)
 157                                 V0=(1-UC)*V0+UC*(E2-ESU(t));
 158                         end;
 159                         V(t)=V0;
 160                 elseif eMode==7
 161                         V0=V(t);
 162                         if E2>ESU(t)
 163                                 V(t) = (1-UC)*V0+UC*(E2-ESU(t));
 164                         else
 165                                 V(t) = V0;
 166                         end;
 167                 elseif eMode==8
 168                         V0=0;
 169                         V(t) = V0; % (t-1)*(1-UC)+UC*E2;
 170                 end;
 171
 172 %[t,size(H),size(Z)]
 173
 174                 k = AY / (ESU(t) + V0);         % Kalman Gain
 175                 zt = zt + k'*E;
 176                 %z(t,:) = zt;
 177
 178                 if aMode==0
 179                         %W = W; %nop                  % Schloegl et al. 2003
 180                 elseif aMode==2
 181                         T(t)=(1-UC)*T(t-1)+UC*(E2-Q(t))/(H'*H);   % Roberts I 1998
 182                         Z=Z*V(t-1)/Q(t);
 183                         if T(t)>0 W=T(t)*eye(MOP); else W=zeros(MOP);end;
 184                 elseif aMode==5
 185                         Q_wo = (H'*C*H + V(t-1));                 % Roberts II 1998
 186                         T(t)=(1-UC)*T(t-1)+UC*(E2-Q_wo)/(H'*H);
 187                         if T(t)>0 W=T(t)*eye(MOP); else W=zeros(MOP); end;
 188                 elseif aMode==6
 189                         T(t)=(1-UC)*T(t-1)+UC*(E2-Q(t))/(H'*H);
 190                         Z=Z*V(t)/Q(t);
 191                         if T(t)>0 W=T(t)*eye(MOP); else W=zeros(MOP); end;
 192                 elseif aMode==11
 193                         %Z = Z - k*AY';
 194                         W = sum(diag(Z))*dW;
 195                 elseif aMode==12
 196                         W = UC*UC*eye(MOP);
 197                 elseif aMode==13
 198                         W = UC*diag(diag(Z));
 199                 elseif aMode==14
 200                         W = (UC*UC)*diag(diag(Z));
 201                 elseif aMode==15
 202                         W = sum(diag(Z))*dW;
 203                 elseif aMode==16
 204                         W = UC*eye(MOP);               % Schloegl 1998
 205                 elseif aMode==17
 206                         Z = 0.5*(Z+Z');
 207                         W = UC*Z;
 208                 elseif aMode==18
 209                         W = 0.5*UC*(Z+Z');
 210                         %W=W;
 211                 end;
 212
 213                 Z = Z - k*AY';               % Schloegl 1998
 214         else
 215
 216                 V(t) = V0;
 217
 218         end;
 219         if any(any(isnan(W))), W=UC*Z; end;
 220
 221         z(t,:) = zt;
 222         Z   = Z + W;               % Schloegl 1998
 223         SPUR(t)=trace(Z);
 224 end;
 225
 226
 227 if 0,m,
 228     z(:,1)=M0*z(:,1)./(1-sum(z(:,2:p),2));
 229 end;
 230
 231 REV = mean(e.*e)/mean(y.*y);
 232 if any(~isfinite(Z(:))), REV=inf; end;
 233