octave_packages/tsa-4.2.4/aar.m

   1 function [a,e,REV,TOC,CPUTIME,ESU] = aar(y, Mode, arg3, arg4, arg5, arg6, arg7, arg8, arg9);
   2 % Calculates adaptive autoregressive (AAR) and adaptive autoregressive moving average estimates (AARMA)
   3 % of real-valued data series using Kalman filter algorithm.
   4 % [a,e,REV] = aar(y, mode, MOP, UC, a0, A, W, V);
   5 %
   6 % The AAR process is described as following
   7 %       y(k) - a(k,1)*y(t-1) -...- a(k,p)*y(t-p) = e(k);
   8 % The AARMA process is described as following
   9 %       y(k) - a(k,1)*y(t-1) -...- a(k,p)*y(t-p) = e(k) + b(k,1)*e(t-1) + ... + b(k,q)*e(t-q);
  10 %
  11 % Input:
  12 %       y       Signal (AR-Process)
  13 %       Mode    is a two-element vector [aMode, vMode],
  14 %               aMode determines 1 (out of 12) methods for updating the co-variance matrix (see also [1])
  15 %               vMode determines 1 (out of 7) methods for estimating the innovation variance (see also [1])
  16 %               aMode=1, vmode=2 is the RLS algorithm as used in [2]
  17 %               aMode=-1, LMS algorithm (signal normalized)
  18 %               aMode=-2, LMS algorithm with adaptive normalization
  19 %
  20 %       MOP     model order, default [10,0]
  21 %               MOP=[p]         AAR(p) model. p AR parameters
  22 %               MOP=[p,q]       AARMA(p,q) model, p AR parameters and q MA coefficients
  23 %       UC      Update Coefficient, default 0
  24 %       a0      Initial AAR parameters [a(0,1), a(0,2), ..., a(0,p),b(0,1),b(0,2), ..., b(0,q)]
  25 %                (row vector with p+q elements, default zeros(1,p) )
  26 %       A       Initial Covariance matrix (positive definite pxp-matrix, default eye(p))
  27 %       W       system noise (required for aMode==0)
  28 %       V       observation noise (required for vMode==0)
  29 %
  30 % Output:
  31 %       a       AAR(MA) estimates [a(k,1), a(k,2), ..., a(k,p),b(k,1),b(k,2), ..., b(k,q]
  32 %       e       error process (Adaptively filtered process)
  33 %       REV     relative error variance MSE/MSY
  34 %
  35 %
  36 % Hint:
  37 % The mean square (prediction) error of different variants is useful for determining the free parameters (Mode, MOP, UC)
  38 %
  39 % REFERENCE(S):
  40 % [1] A. Schloegl (2000), The electroencephalogram and the adaptive autoregressive model: theory and applications.
  41 %     ISBN 3-8265-7640-3 Shaker Verlag, Aachen, Germany.
  42 %
  43 % More references can be found at
  44 %     http://www.dpmi.tu-graz.ac.at/~schloegl/publications/
  45
  46 %
  47 %       $Id: aar.m 8383 2011-07-16 20:06:59Z schloegl $
  48 %       Copyright (C) 1998-2003 by Alois Schloegl <a.schloegl@ieee.org>
  49 %
  50 %    This program is free software: you can redistribute it and/or modify
  51 %    it under the terms of the GNU General Public License as published by
  52 %    the Free Software Foundation, either version 3 of the License, or
  53 %    (at your option) any later version.
  54 %
  55 %    This program is distributed in the hope that it will be useful,
  56 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  57 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  58 %    GNU General Public License for more details.
  59 %
  60 %    You should have received a copy of the GNU General Public License
  61 %    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  62
  63
  64 [nc,nr]=size(y);
  65 %if nc<nr y=y'; end; tmp=nr;nc=nr; nr=tmp;end;
  66
  67 if nargin<2 Mode=0; end;
  68 % check Mode (argument2)
  69 if prod(size(Mode))==2
  70         aMode=Mode(1);
  71         vMode=Mode(2);
  72 end;
  73 if any(aMode==(0:14)) && any(vMode==(0:7)),
  74         fprintf(1,['a' int2str(aMode) 'e' int2str(vMode) ' ']);
  75 else
  76         fprintf(2,'Error AAR.M: invalid Mode argument\n');
  77         return;
  78 end;
  79
  80 % check model order (argument3)
  81 if nargin<3 MOP=[10,0]; else MOP= arg3; end;
  82 if length(MOP)==0 p=10; q=0; MOP=p;
  83 elseif length(MOP)==1 p=MOP(1); q=0; MOP=p;
  84 elseif length(MOP)>=2 p=MOP(1); q=MOP(2); MOP=p+q;
  85 end;
  86
  87 if nargin<4 UC=0; else UC= arg4; end;
  88
  89 a0=zeros(1,MOP);
  90 A0=eye(MOP);
  91 if nargin>4,
  92         if all(size(arg5)==([1,1]*(MOP+1)));    % extended covariance matrix of AAR parameters
  93                 a0 = arg5(1,2:size(arg5,2));
  94                 A0 = arg5(2:size(arg5,1),2:size(arg5,2)) - a0'*a0;
  95         else
  96                 a0 = arg5;
  97                 if nargin>5
  98                         A0 = arg6;
  99                 end;
 100         end;
 101 end;
 102
 103 if nargin<7, W  = []; else W  = arg7; end;
 104
 105 if all(size(W)==MOP),
 106         if aMode ~= 0,
 107                 fprintf(1,'aMode should be 0, because W is given.\n');
 108         end;
 109 elseif isempty(W),
 110         if aMode == 0,
 111                 fprintf(1,'aMode must be non-zero, because W is not given.\n');
 112         end;
 113 elseif any(size(W)~=MOP),
 114         fprintf(1,'size of W does not fit. It must be %i x %i.\n',MOP,MOP);
 115         return;
 116 end;
 117
 118 if nargin<8, V0 = []; else V0 = arg8; end;
 119 if all(size(V0)==nr),
 120         if vMode ~= 0,
 121                 fprintf(1,'vMode should be 0, because V is given.\n');
 122         end;
 123 elseif isempty(V0),
 124         if aMode == 0,
 125                 fprintf(1,'vMode must be non-zero, because V is not given.\n');
 126         end;
 127 else
 128         fprintf(1,'size of V does not fit. It must be 1x1.\n');
 129         return;
 130 end;
 131
 132 % if nargin<7 TH=3; else TH = arg7;  end;
 133 %       TH=TH*var(y);
 134 %       TH=TH*mean(detrend(y,0).^2);
 135 MSY=mean(detrend(y,0).^2);
 136
 137 e=zeros(nc,1);
 138 Q=zeros(nc,1);
 139 V=zeros(nc,1);
 140 T=zeros(nc,1);
 141 %DET=zeros(nc,1);
 142 SPUR=zeros(nc,1);
 143 ESU=zeros(nc,1);
 144 a=a0(ones(nc,1),:);
 145 %a=zeros(nc,MOP);
 146 %b=zeros(nc,q);
 147
 148 mu=1-UC; % Patomaeki 1995
 149 lambda=(1-UC); % Schloegl 1996
 150 arc=poly((1-UC*2)*[1;1]);b0=sum(arc); % Whale forgettting factor for Mode=258,(Bianci et al. 1997)
 151
 152 dW=UC/MOP*eye(MOP);                % Schloegl
 153
 154
 155 %------------------------------------------------
 156 %       First Iteration
 157 %------------------------------------------------
 158 Y=zeros(MOP,1);
 159 C=zeros(MOP);
 160 %X=zeros(q,1);
 161 at=a0;
 162 A=A0;
 163 E=y(1);
 164 e(1)=E;
 165 if ~isempty(V0)
 166         V(1) = V0;
 167 else
 168         V(1) = (1-UC) + UC*E*E;
 169 end;
 170 ESU(1) = 1; %Y'*A*Y;
 171
 172 A1=zeros(MOP);A2=A1;
 173 tic;CPUTIME=cputime;
 174 %------------------------------------------------
 175 %       Update Equations
 176 %------------------------------------------------
 177 T0=2;
 178
 179 for t=T0:nc,
 180
 181         %Y=[y(t-1); Y(1:p-1); E ; Y(p+1:MOP-1)]
 182
 183         if t<=p Y(1:t-1)=y(t-1:-1:1);           % Autoregressive
 184         else    Y(1:p)=y(t-1:-1:t-p);
 185         end;
 186
 187         if t<=q Y(p+(1:t-1))=e(t-1:-1:1);       % Moving Average
 188         else    Y(p+1:MOP)=e(t-1:-1:t-q);
 189         end;
 190
 191         % Prediction Error
 192         E = y(t) - a(t-1,:)*Y;
 193         e(t) = E;
 194         E2=E*E;
 195
 196         AY=A*Y;
 197         esu=Y'*AY;
 198         ESU(t)=esu;
 199
 200         if isnan(E),
 201                 a(t,:)=a(t-1,:);
 202         else
 203                 V(t) = V(t-1)*(1-UC)+UC*E2;
 204                 if aMode == -1, % LMS
 205                         %       V(t) = V(t-1)*(1-UC)+UC*E2;
 206                         a(t,:)=a(t-1,:) + (UC/MSY)*E*Y';
 207                 elseif aMode == -2, % LMS with adaptive estimation of the variance
 208                         a(t,:)=a(t-1,:) + UC/V(t)*E*Y';
 209
 210                 else    % Kalman filtering (including RLS)
 211                         if vMode==0,            %eMode==4
 212                                 Q(t) = (esu + V0);
 213                         elseif vMode==1,            %eMode==4
 214                                 Q(t) = (esu + V(t));
 215                         elseif vMode==2,        %eMode==2
 216                                 Q(t) = (esu + 1);
 217                         elseif vMode==3,        %eMode==3
 218                                 Q(t) = (esu + lambda);
 219                         elseif vMode==4,        %eMode==1
 220                                 Q(t) = (esu + V(t-1));
 221                         elseif vMode==5,        %eMode==6
 222                                 if E2>esu
 223                                         V(t)=(1-UC)*V(t-1)+UC*(E2-esu);
 224                                 else
 225                                         V(t)=V(t-1);
 226                                 end;
 227                                 Q(t) = (esu + V(t));
 228                         elseif vMode==6,        %eMode==7
 229                                 if E2>esu
 230                                         V(t)=(1-UC)*V(t-1)+UC*(E2-esu);
 231                                 else
 232                                         V(t)=V(t-1);
 233                                 end;
 234                                 Q(t) = (esu + V(t-1));
 235                         elseif vMode==7,        %eMode==8
 236                                 Q(t) = esu;
 237                         end;
 238
 239                         k = AY / Q(t);          % Kalman Gain
 240                         a(t,:) = a(t-1,:) + k'*E;
 241
 242                         if aMode==0,                    %AMode=0
 243                                 A = A - k*AY' + W;                   % Schloegl et al. 2003
 244                         elseif aMode==1,                    %AMode=1
 245                                 A = (1+UC)*(A - k*AY');                   % Schloegl et al. 1997
 246                         elseif aMode==2,                %AMode=11
 247                                 A = A - k*AY';
 248                                 A = A + sum(diag(A))*dW;
 249                         elseif aMode==3,                %AMode=5
 250                                 A = A - k*AY' + sum(diag(A))*dW;
 251                         elseif aMode==4,                %AMode=6
 252                                 A = A - k*AY' + UC*eye(MOP);               % Schloegl 1998
 253                         elseif aMode==5,                %AMode=2
 254                                 A = A - k*AY' + UC*UC*eye(MOP);
 255                         elseif aMode==6,                %AMode=2
 256                                 T(t)=(1-UC)*T(t-1)+UC*(E2-Q(t))/(Y'*Y);
 257                                 A=A*V(t-1)/Q(t);
 258                                 if T(t)>0 A=A+T(t)*eye(MOP); end;
 259                         elseif aMode==7,                %AMode=6
 260                                 T(t)=(1-UC)*T(t-1)+UC*(E2-Q(t))/(Y'*Y);
 261                                 A=A*V(t)/Q(t);
 262                                 if T(t)>0 A=A+T(t)*eye(MOP); end;
 263                         elseif aMode==8,                %AMode=5
 264                                 Q_wo = (Y'*C*Y + V(t-1));
 265                                 C=A-k*AY';
 266                                 T(t)=(1-UC)*T(t-1)+UC*(E2-Q_wo)/(Y'*Y);
 267                                 if T(t)>0 A=C+T(t)*eye(MOP); else A=C; end;
 268                         elseif aMode==9,                %AMode=3
 269                                 A = A - (1+UC)*k*AY';
 270                                 A = A + sum(diag(A))*dW;
 271                         elseif aMode==10,               %AMode=7
 272                                 A = A - (1+UC)*k*AY' + sum(diag(A))*dW;
 273                         elseif aMode==11,               %AMode=8
 274
 275                                 A = A - (1+UC)*k*AY' + UC*eye(MOP);        % Schloegl 1998
 276                         elseif aMode==12,               %AMode=4
 277                                 A = A - (1+UC)*k*AY' + UC*UC*eye(MOP);
 278                         elseif aMode==13
 279                                 A = A - k*AY' + UC*diag(diag(A));
 280                         elseif aMode==14
 281                                 A = A - k*AY' + (UC*UC)*diag(diag(A));
 282                         end;
 283                 end;
 284         end;
 285 end;
 286
 287 %a=a(end,:);
 288 TOC = toc;
 289 CPUTIME = cputime - CPUTIME;
 290 %REV = (e'*e)/(y'*y);
 291
 292 REV = mean(e.*e)./mean(y.*y);