octave_packages/tsa-4.2.4/aarmam.m

   1 function [z,e,REV,ESU,V,Z,SPUR] = aarmam(y, Mode, MOP, UC, z0, Z0, V0, W);
   2 % Estimating Adaptive AutoRegressive-Moving-Average-and-mean model (includes mean term)
   3 %
   4 % !! This function is obsolete and is replaced by AMARMA
   5 %
   6 % [z,E,REV,ESU,V,Z,SPUR] = aarmam(y, mode, MOP, UC, z0, Z0, V0, W);
   7 % Estimates AAR parameters with Kalman filter algorithm
   8 %       y(t) = sum_i(a_i(t)*y(t-i)) + m(t) + e(t) + sum_i(b_i(t)*e(t-i))
   9 %
  10 % State space model
  11 %       z(t) = G*z(t-1) + w(t)    w(t)=N(0,W)
  12 %       y(t) = H*z(t)   + v(t)    v(t)=N(0,V)
  13 %
  14 % G = I,
  15 % z = [m(t),a_1(t-1),..,a_p(t-p),b_1(t-1),...,b_q(t-q)];
  16 % H = [1,y(t-1),..,y(t-p),e(t-1),...,e(t-q)];
  17 % W = E{(z(t)-G*z(t-1))*(z(t)-G*z(t-1))'}
  18 % V = E{(y(t)-H*z(t-1))*(y(t)-H*z(t-1))'}
  19 %
  20 %
  21 % Input:
  22 %       y       Signal (AR-Process)
  23 %       Mode    determines the type of algorithm
  24 %
  25 %       MOP     Model order [m,p,q], default [0,10,0]
  26 %                       m=1 includes the mean term, m=0 does not.
  27 %                       p and q must be positive integers
  28 %                       it is recommended to set q=0.
  29 %       UC      Update Coefficient, default 0
  30 %       z0      Initial state vector
  31 %       Z0      Initial Covariance matrix
  32 %
  33 % Output:
  34 %       z       AR-Parameter
  35 %       E       error process (Adaptively filtered process)
  36 %       REV     relative error variance MSE/MSY
  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 %
  42 % More references can be found at
  43 %     http://pub.ist.ac.at/~schloegl/publications/
  44
  45 %       $Id: aarmam.m 9609 2012-02-10 10:18:00Z schloegl $
  46 %       Copyright (C) 1998-2002,2008,2012 by Alois Schloegl <alois.schloegl@ist.ac.at>
  47 %       This is part of the TSA-toolbox. See also
  48 %       http://pub.ist.ac.at/~schloegl/matlab/tsa/
  49 %       http://octave.sourceforge.net/
  50 %       http://biosig.sourceforge.net/
  51 %
  52 %    This program is free software: you can redistribute it and/or modify
  53 %    it under the terms of the GNU General Public License as published by
  54 %    the Free Software Foundation, either version 3 of the License, or
  55 %    (at your option) any later version.
  56 %
  57 %    This program is distributed in the hope that it will be useful,
  58 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  59 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  60 %    GNU General Public License for more details.
  61 %
  62 %    You should have received a copy of the GNU General Public License
  63 %    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  64
  65 %#realonly
  66 %#inbounds
  67
  68 warning('AARMAM is obsolete. Use AMARMA instead!')
  69
  70 [nc,nr]=size(y);
  71
  72 if nargin<2 Mode=0;
  73 elseif ischar(Mode) Mode=bin2dec(Mode);
  74 elseif isnan(Mode) return; end;
  75 if nargin<3 MOP=[0,10,0]; end;
  76 if length(MOP)==0,      m=0;p=10; q=0; MOP=p;
  77 elseif length(MOP)==1,  m=0;p=MOP(1); q=0; MOP=p;
  78 elseif length(MOP)==2,  fprintf(1,'Error AMARMA: MOP is ambiguos\n');
  79 elseif length(MOP)>2,   m=MOP(1); p=MOP(2); q=MOP(3);MOP=m+p+q;
  80 end;
  81
  82 if prod(size(Mode))>1
  83         aMode=Mode(1);
  84         eMode=Mode(2);
  85 end;
  86 %fprintf(1,['a' int2str(aMode) 'e' int2str(eMode) ' ']);
  87
  88
  89 e=zeros(nc,1);
  90 V=zeros(nc,1);V(1)=V0;
  91 T=zeros(nc,1);
  92 ESU=zeros(nc,1)+nan;
  93 SPUR=zeros(nc,1)+nan;
  94 z=z0(ones(nc,1),:);
  95
  96 arc=poly((1-UC*2)*[1;1]);b0=sum(arc); % Whale forgetting factor for Mode=258,(Bianci et al. 1997)
  97
  98 dW=UC/MOP*eye(MOP);                % Schloegl
  99
 100
 101 %------------------------------------------------
 102 %       First Iteration
 103 %------------------------------------------------
 104
 105 H = zeros(MOP,1);
 106 if m,
 107         H(1) = 1;%M0;
 108         if m~=1,
 109                 fprintf(2,'Warning AARMAM: m must be 0 or 1\n');
 110                 return;
 111         end;
 112 end;
 113 if (p<0) || (q<0) || (round(p)~=p) || (round(q)~=q),
 114         fprintf(2,'Error AARMAM: p and q must be positive integers\n');
 115         return;
 116 end;
 117
 118 E = 0;
 119 Z = Z0;
 120 zt= z0;
 121
 122 A1 = zeros(MOP); A2 = A1;
 123
 124 y_1=0;
 125
 126 %------------------------------------------------
 127 %       Update Equations
 128 %------------------------------------------------
 129
 130 for t=1:nc,
 131
 132         % make measurement matrix
 133         if 0,
 134                 if t>1,
 135                         y_1 = y(t-1);
 136                 end;
 137                 H=[1; y_1; H(m+(1:p-1)'); E(1:min(1,q-1)) ; H(p+m+(1:q-1)')];  % shift y and e
 138
 139         else    % this seem to be slightly faster
 140                 if t<=p, H(m+(1:t-1)) = y(t-1:-1:1);    % Autoregressive
 141                 else     H(m+(1:p))   = y(t-1:-1:t-p);
 142                 end;
 143
 144                 if t<=q, H(m+p+(1:t-1)) = e(t-1:-1:1);  % Moving Average
 145                 else     H(m+p+(1:q)) = e(t-1:-1:t-q);
 146                 end;
 147         end;
 148
 149         % Prediction Error
 150         E = y(t) - zt*H;
 151         e(t) = E;
 152
 153         if ~isnan(E),
 154                 E2 = E*E;
 155                 AY = Z*H;
 156
 157                 ESU(t) = H'*AY;
 158
 159                 if eMode==1
 160                         V0 = V(t-1);
 161                         V(t) = V0*(1-UC)+UC*E2;
 162                 elseif eMode==2
 163                         V0 = 1;
 164                         V(t) = V0; %V(t-1)*(1-UC)+UC*E2;
 165                 elseif eMode==3
 166                         V0 = 1-UC;
 167                         V(t) = V0; %(t-1)*(1-UC)+UC*E2;
 168                 elseif eMode==4
 169                         V0 = V0*(1-UC)+UC*E2;
 170                         V(t) = V0;
 171                 elseif eMode==5
 172                         V(t)=V0;
 173                         %V0 = V0;
 174                 elseif eMode==6
 175                         if E2>ESU(t),
 176                                 V0=(1-UC)*V0+UC*(E2-ESU(t));
 177                         end;
 178                         V(t)=V0;
 179                 elseif eMode==7
 180                         V0=V(t);
 181                         if E2>ESU(t)
 182                                 V(t) = (1-UC)*V0+UC*(E2-ESU(t));
 183                         else
 184                                 V(t) = V0;
 185                         end;
 186                 elseif eMode==8
 187                         V0=0;
 188                         V(t) = V0; % (t-1)*(1-UC)+UC*E2;
 189                 end;
 190
 191                 k = AY / (ESU(t) + V0);         % Kalman Gain
 192                 zt = zt + k'*E;
 193                 %z(t,:) = zt;
 194
 195                 if aMode==2
 196                         T(t)=(1-UC)*T(t-1)+UC*(E2-Q(t))/(H'*H);   % Roberts I 1998
 197                         Z=Z*V(t-1)/Q(t);
 198                         if T(t)>0 W=T(t)*eye(MOP); else W=zeros(MOP);end;
 199                 elseif aMode==5
 200                         Q_wo = (H'*C*H + V(t-1));                 % Roberts II 1998
 201                         T(t)=(1-UC)*T(t-1)+UC*(E2-Q_wo)/(H'*H);
 202                         if T(t)>0 W=T(t)*eye(MOP); else W=zeros(MOP); end;
 203                 elseif aMode==6
 204                         T(t)=(1-UC)*T(t-1)+UC*(E2-Q(t))/(H'*H);
 205                         Z=Z*V(t)/Q(t);
 206                         if T(t)>0 W=T(t)*eye(MOP); else W=zeros(MOP); end;
 207                 elseif aMode==11
 208                         %Z = Z - k*AY';
 209                         W = sum(diag(Z))*dW;
 210                 elseif aMode==12
 211                         W = UC*UC*eye(MOP);
 212                 elseif aMode==13
 213                         W = UC*diag(diag(Z));
 214                 elseif aMode==14
 215                         W = (UC*UC)*diag(diag(Z));
 216                 elseif aMode==15
 217                         W = sum(diag(Z))*dW;
 218                 elseif aMode==16
 219                         W = UC*eye(MOP);               % Schloegl 1998
 220                         %elseif aMode==17
 221                         %W=W;
 222                 end;
 223
 224                 Z = Z - k*AY';               % Schloegl 1998
 225         else
 226
 227                 V(t) = V0;
 228
 229         end;
 230
 231         if any(any(isnan(W))), W=UC*Z;end;
 232
 233         z(t,:) = zt;
 234         Z   = Z + W;               % Schloegl 1998
 235         SPUR(t) = trace(Z);
 236 end;
 237
 238 REV = mean(e.*e)/mean(y.*y);
 239 if any(~isfinite(Z(:))), REV=inf; end;
 240