aarmam.m

时间序列分析的matlab程序

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 <[e-mail]>
%
% 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;