www.gusucode.com > 时间序列分析工具箱 - tsa源码程序 > tsa/amarma.m
function [z,e,REV,ESU,V,Z,SPUR] = amarma(y, Mode, MOP, UC, z0, Z0, V0, W);
% Adaptive Mean-AutoRegressive-Moving-Average model estimation
% [z,E,ESU,REV,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)) + mu(t) + E(t)
%
% 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 = [?(t)/(1-sum_i(a(i,t))),a_1(t-1),..,a_p(t-p),b_1(t-1),...,b_q(t-q)];
% H = [M0/(1-sum_i(a(i,t))),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
% [0,0] uses V0 and W
%
% MOP Model order [m,p,q], default [0,10,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/
% $Id: amarma.m,v 1.1 2005/11/07 14:42:26 schloegl Exp $
% Copyright (c) 1998-2002,2005 by Alois Schloegl <a.schloegl@ieee.org>
%
% 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
%#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 nargin<8, W = nan ; 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,
%M0 = z0(1)/(1-sum(z0(2:p+1))); %transformierter Mittelwert
H(1) = 1;%M0;
%z0(1)= 1;
end;
Z = Z0;
zt= z0;
A1 = zeros(MOP); A2 = A1;
%------------------------------------------------
% Update Equations
%------------------------------------------------
for t=1:nc,
%H=[y(t-1); H(1:p-1); E ; H(p+1:MOP-1)]
if t<=p, H(m+(1:t-1)) = y(t-1:-1:1); %H(p)=mu0; % Autoregressive
else H(m+(1:p)) = y(t-1:-1:t-p); %mu0];
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;
% Prediction Error
E = y(t) - zt*H;
e(t) = E;
if ~isnan(E),
E2 = E*E;
AY = Z*H;
% [zt, t, y(t), E,ESU(t),V(t),H,Z],pause,
ESU(t) = H'*AY;
if eMode==0
V(t) = V0;
elseif 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;
%[t,size(H),size(Z)]
k = AY / (ESU(t) + V0); % Kalman Gain
zt = zt + k'*E;
%z(t,:) = zt;
if aMode==0
%W = W; %nop % Schloegl et al. 2003
elseif 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;
if 0,m,
z(:,1)=M0*z(:,1)./(1-sum(z(:,2:p),2));
end;
REV = mean(e.*e)/mean(y.*y);
if any(~isfinite(Z(:))), REV=inf; end;