www.gusucode.com > 时间序列分析工具箱 - tsa源码程序 > tsa/aarmam.m
function [z,e,REV,ESU,V,Z,SPUR] = aarmam(y, Mode, MOP, UC, z0, Z0, V0, W);
% Estimating Adaptive AutoRegressive-Moving-Average-and-mean model (includes mean term)
% [z,E,REV,ESU,V,Z,SPUR] = amarma(y, mode, MOP, UC, z0, Z0, V0, W);
% Estimates AAR parameters with Kalman filter algorithm
% y(t) = sum_i(a_i(t)*y(t-i)) + m(t) + e(t) + sum_i(b_i(t)*e(t-i))
%
% State space model
% z(t) = G*z(t-1) + w(t) w(t)=N(0,W)
% y(t) = H*z(t) + v(t) v(t)=N(0,V)
%
% G = I,
% z = [m(t),a_1(t-1),..,a_p(t-p),b_1(t-1),...,b_q(t-q)];
% H = [1,y(t-1),..,y(t-p),e(t-1),...,e(t-q)];
% W = E{(z(t)-G*z(t-1))*(z(t)-G*z(t-1))'}
% V = E{(y(t)-H*z(t-1))*(y(t)-H*z(t-1))'}
%
%
% Input:
% y Signal (AR-Process)
% Mode determines the type of algorithm
%
% MOP Model order [m,p,q], default [0,10,0]
% m=1 includes the mean term, m=0 does not.
% p and q must be positive integers
% it is recommended to set q=0.
% UC Update Coefficient, default 0
% z0 Initial state vector
% Z0 Initial Covariance matrix
%
% Output:
% z AR-Parameter
% E error process (Adaptively filtered process)
% REV relative error variance MSE/MSY
%
% REFERENCE(S):
% [1] A. Schloegl (2000), The electroencephalogram and the adaptive autoregressive model: theory and applications.
% ISBN 3-8265-7640-3 Shaker Verlag, Aachen, Germany.
%
% More references can be found at
% http://www.dpmi.tu-graz.ac.at/~schloegl/publications/
% Version 2.99 07 Aug 2002
% Copyright (C) 1998-2002 by Alois Schloegl <a.schloegl@ieee.org>
%
% This library is free software; you can redistribute it and/or
% modify it under the terms of the GNU Library General Public
% License as published by the Free Software Foundation; either
% version 2 of the License, or (at your option) any later version.
%
% This library is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
% Library General Public License for more details.
%
% You should have received a copy of the GNU Library General Public
% License along with this library; if not, write to the
% Free Software Foundation, Inc., 59 Temple Place - Suite 330,
% Boston, MA 02111-1307, USA.
% 12.04.1999 ESU included
% 19.10.2000 aMode 13,14 included
% NaN handling
% 06.07.2002 Docu improved, included into TSA
% 11.07.2002 nanmean replaced by mean
% 07.08.2002 cosmetic changes
%#realonly
%#inbounds
[nc,nr]=size(y);
if nargin<2 Mode=0;
elseif ischar(Mode) Mode=bin2dec(Mode);
elseif isnan(Mode) return; end;
if nargin<3 MOP=[0,10,0]; end;
if length(MOP)==0, m=0;p=10; q=0; MOP=p;
elseif length(MOP)==1, m=0;p=MOP(1); q=0; MOP=p;
elseif length(MOP)==2, fprintf(1,'Error AMARMA: MOP is ambiguos\n');
elseif length(MOP)>2, m=MOP(1); p=MOP(2); q=MOP(3);MOP=m+p+q;
end;
if prod(size(Mode))>1
aMode=Mode(1);
eMode=Mode(2);
end;
%fprintf(1,['a' int2str(aMode) 'e' int2str(eMode) ' ']);
e=zeros(nc,1);
V=zeros(nc,1);V(1)=V0;
T=zeros(nc,1);
ESU=zeros(nc,1)+nan;
SPUR=zeros(nc,1)+nan;
z=z0(ones(nc,1),:);
arc=poly((1-UC*2)*[1;1]);b0=sum(arc); % Whale forgetting factor for Mode=258,(Bianci et al. 1997)
dW=UC/MOP*eye(MOP); % Schloegl
%------------------------------------------------
% First Iteration
%------------------------------------------------
H = zeros(MOP,1);
if m,
H(1) = 1;%M0;
if m~=1,
fprintf(2,'Warning AARMAM: m must be 0 or 1\n');
return;
end;
end;
if (p<0) | (q<0) | (round(p)~=p) | (round(q)~=q),
fprintf(2,'Error AARMAM: p and q must be positive integers\n');
return;
end;
E = 0;
Z = Z0;
zt= z0;
A1 = zeros(MOP); A2 = A1;
y_1=0;
%------------------------------------------------
% Update Equations
%------------------------------------------------
for t=1:nc,
% make measurement matrix
if 0,
if t>1,
y_1 = y(t-1);
end;
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
else % this seem to be slightly faster
if t<=p, H(m+(1:t-1)) = y(t-1:-1:1); % Autoregressive
else H(m+(1:p)) = y(t-1:-1:t-p);
end;
if t<=q, H(m+p+(1:t-1)) = e(t-1:-1:1); % Moving Average
else H(m+p+(1:q)) = e(t-1:-1:t-q);
end;
end;
% Prediction Error
E = y(t) - zt*H;
e(t) = E;
if ~isnan(E),
E2 = E*E;
AY = Z*H;
ESU(t) = H'*AY;
if eMode==1
V0 = V(t-1);
V(t) = V0*(1-UC)+UC*E2;
elseif eMode==2
V0 = 1;
V(t) = V0; %V(t-1)*(1-UC)+UC*E2;
elseif eMode==3
V0 = 1-UC;
V(t) = V0; %(t-1)*(1-UC)+UC*E2;
elseif eMode==4
V0 = V0*(1-UC)+UC*E2;
V(t) = V0;
elseif eMode==5
V(t)=V0;
%V0 = V0;
elseif eMode==6
if E2>ESU(t),
V0=(1-UC)*V0+UC*(E2-ESU(t));
end;
V(t)=V0;
elseif eMode==7
V0=V(t);
if E2>ESU(t)
V(t) = (1-UC)*V0+UC*(E2-ESU(t));
else
V(t) = V0;
end;
elseif eMode==8
V0=0;
V(t) = V0; % (t-1)*(1-UC)+UC*E2;
end;
k = AY / (ESU(t) + V0); % Kalman Gain
zt = zt + k'*E;
%z(t,:) = zt;
if aMode==2
T(t)=(1-UC)*T(t-1)+UC*(E2-Q(t))/(H'*H); % Roberts I 1998
Z=Z*V(t-1)/Q(t);
if T(t)>0 W=T(t)*eye(MOP); else W=zeros(MOP);end;
elseif aMode==5
Q_wo = (H'*C*H + V(t-1)); % Roberts II 1998
T(t)=(1-UC)*T(t-1)+UC*(E2-Q_wo)/(H'*H);
if T(t)>0 W=T(t)*eye(MOP); else W=zeros(MOP); end;
elseif aMode==6
T(t)=(1-UC)*T(t-1)+UC*(E2-Q(t))/(H'*H);
Z=Z*V(t)/Q(t);
if T(t)>0 W=T(t)*eye(MOP); else W=zeros(MOP); end;
elseif aMode==11
%Z = Z - k*AY';
W = sum(diag(Z))*dW;
elseif aMode==12
W = UC*UC*eye(MOP);
elseif aMode==13
W = UC*diag(diag(Z));
elseif aMode==14
W = (UC*UC)*diag(diag(Z));
elseif aMode==15
W = sum(diag(Z))*dW;
elseif aMode==16
W = UC*eye(MOP); % Schloegl 1998
%elseif aMode==17
%W=W;
end;
Z = Z - k*AY'; % Schloegl 1998
else
V(t) = V0;
end;
if any(any(isnan(W))), W=UC*Z;end;
z(t,:) = zt;
Z = Z + W; % Schloegl 1998
SPUR(t) = trace(Z);
end;
REV = mean(e.*e)/mean(y.*y);
if any(~isfinite(Z(:))), REV=inf; end;