カルマン、Matlab、および状態空間モデル

clear all; clc; close all; format long; %------------------- 1. Load and prepare data ------------------------------ load otgruzka.mat; % Structure ssmopt contains several important records used throughout the code ssmopt.frcst=70; % forecast length yend=y(end-ssmopt.frcst+1:end); % saved last observation for the forecast comparison y=y(1:end-ssmopt.frcst); ssmopt.N=length(y); ssmopt.trans=1; % transform estimated parameters to preserve positiveness of variances ssmopt.sv=[5e+8;500;5e+8;5e+8;0.05;0.9]; % starting values ssmopt.mle='f_obj'; % name of the objective function for the maximization ssmopt.sv=transform(ssmopt.sv,'in',ssmopt); ssmopt.filter='kfmy2'; % name of the function computing Kalman Filter ssmopt.smooth='ksmy2'; % name of the function computing Kalman Smoother %------------------- 2. Estimate the model ------------------------------ result=mle_my(y,ssmopt); % call Maximum Likelihood maximization function b=transform(result.b,'out',ssmopt); % transform parameters back % recompute data based on the correct non-transformed parameters ssmopt.trans=0; ssmopt.sv=b; [LH,KF_out,Ksm_out] = feval(ssmopt.mle,b,y, ssmopt); % Recover filtered states series - trend, cycle, and seasonal a_trend=KF_out.Afilt(1,:); a_cycle=KF_out.Afilt(3,:); a_seas=KF_out.Afilt(5,:)+KF_out.Afilt(7,:)+KF_out.Afilt(9,:); y_filt=a_trend+a_cycle+a_seas; % build the estimated filtered series Y % Recover smoothed states series - trend, cycle, and seasonal a_trendsm=Ksm_out.Asm(1,:); a_cyclesm=Ksm_out.Asm(3,:); a_seassm=Ksm_out.Asm(5,:)+Ksm_out.Asm(7,:)+Ksm_out.Asm(9,:); y_smooth=a_trendsm+a_cyclesm+a_seassm; % build the estimated smoothed series Y result=mle_my(y,ssmopt); % find correct Hessian for non-transformed parameters %------------------- 3. Compute estimation statistics ------------------------------ %Find standard errors, and p-values se=sqrt(abs(diag(inv(result.hessian)/ssmopt.N))); % se(b) tstat=b./se; % t-statistics pval=2*(1-tcdf(abs(tstat),ssmopt.N-length(ssmopt.sv))); % p-value % Display output fprintf('Estimated parameters and p-values:\n'); out=[b pval] period=2*pi/b(end-1) % Compute R-squared resid=y-y_filt; % estimation errors SSE=resid*resid'; % Sum of Squared Errors SST=(y-mean(y))*(y-mean(y))'; % Sum of Squares Total R2=1-SSE/SST % R-squared' %------------------- 4. Make Forecast ------------------------------ af0=KF_out.Afilt(:,end-ssmopt.frcst); [yf,af]=frcst(b,y,ssmopt, af0); %------------------- 5. Plot results ------------------------------ %p=ssmopt.N; p=600; t=[1:1:p]; figure(1) plot(t,y(1:p),'k',t,y_filt(1:p),'b',t,y_smooth(1:p),'r--') title('Original, Filtered, and Smoothed data') legend('y(t)','y filtered','y smoothed'); figure(2) plot(t,y(1:p),'k',t,a_trend(1:p),'b',t,a_trendsm(1:p),'r--') title('Original data, Filtered and Smoothed trend') legend('y(t)','Filtered trend','Smoothed trend'); figure(3) plot(t,a_cycle(1:p),'b',t,a_cyclesm(1:p),'r--') title('Filtered and Smoothed cycle') legend('Filtered cycle','Smoothed cycle'); figure(4) % Filtered + smoothed seasonal plot(t,a_seas(1:p),'b',t,a_seassm(1:p),'r--') title('Filtered and Smoothed weekly seasonal') legend('Filtered seasonal','Smoothed seasonal'); t=[1:1:ssmopt.frcst]; figure(5) plot(t,yend,'b',t,yf,'r--') title('Original data and Forecast') legend('Original data','Forecast'); RMSE=sqrt(sum((yend - yf).^2)/ssmopt.frcst) % Root Mean Squared Error
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

 function result_mle=mle_my(y,mleopt); warning off; %---------------- 1. Set-up minimization options ---------------- options=optimset('fminunc'); options=optimset('LargeScale', 'off' , ... 'HessUpdate', 'bfgs' , ... 'LineSearchType', 'quadcubic' , ... 'GradObj' , 'off' , ... 'Display','off' , ... 'MaxIter' , 1000 , ... 'TolX', 1e-12 , ... 'TolFun' , 1e-12, ... 'DerivativeCheck' , 'off' , ... 'Diagnostics' , 'off' , ... 'MaxFunEvals', 1000); %---------------- 2. Run minimization ---------------- [b,fval,exitflag,output,grad,hessian]=fminunc(mleopt.mle,mleopt.sv,options,y,mleopt); warning on; result_mle.b=b; result_mle.fval=fval; result_mle.output=output; result_mle.hessian=hessian;
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

 function [obj,KF_out,Ksm_out]=f_obj_tr(b,y,ssmopt); %---------------- 1. Recover parameters ------------------------------------ b=transform(b,'out',ssmopt); s2_irr=b(1); s2_tr=b(2); s2_cyc=b(3); s2_seas=b(4); freq=b(5); rho=b(6); %---------------- 2. Build the model ------------------------------------ ssmopt.ssmodel.states=10; ssmopt.ssmodel.Z=[1 0 1 0 1 0 1 0 1 0]; T1 = [1 1 0 0; 0 1 0 0; 0 0 rho*cos(freq) rho*sin(freq); 0 0 -rho*sin(freq) rho*cos(freq)]; T2=[cos(2*pi/7) sin(2*pi/7) 0 0 0 0;... -sin(2*pi/7) cos(2*pi/7) 0 0 0 0;... 0 0 cos(4*pi/7) sin(4*pi/7) 0 0;... 0 0 -sin(4*pi/7) cos(4*pi/7) 0 0;... 0 0 0 0 cos(6*pi/7) sin(6*pi/7);... 0 0 0 0 -sin(6*pi/7) cos(6*pi/7)]; ssmopt.ssmodel.T = [T1 zeros(rows(T1),cols(T2));zeros(rows(T2),cols(T1)) T2]; ssmopt.ssmodel.R=eye(10); ssmopt.ssmodel.R(1,1)=0; H=s2_irr; Q=zeros(10,10); Q(2,2)=s2_tr; Q(3,3)=s2_cyc; Q(4,4)=s2_cyc; Q(5,5)=s2_seas; Q(6,6)=s2_seas; Q(7,7)=s2_seas; Q(8,8)=s2_seas; Q(9,9)=s2_seas; Q(10,10)=s2_seas; %---------------- 3. Suggest starting conditions for the states ------------------------ a0=[y(1);0;0;0;0;0;0;0;0;0]; P0=eye(ssmopt.ssmodel.states)*1e+10; %---------------- 4. Run Kalman filter ------------------------ KF_out = feval(ssmopt.filter,y, ssmopt, Q, H, a0, P0); obj=KF_out.LH; %---------------- 5. Run Kalman smoother ------------------------ if nargout>2 ssmopt.ssmodel.num_etas=3; % number of the state variances (required for Kalman smoother) Ksm_out = feval(ssmopt.smooth,KF_out, ssmopt); end
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

 % Kalman filter % y[t] = Z*alpha[t] + eps, eps ~ N(0,H). % alpha[t] = T*alpha[t-1] + R*eta, eta ~ N(0,Q). % v[t]=y[t]-E(y[t]) = y[t]-Z*a[t] % F[t]=var(v[t]) function KF_out = kfmy_koop(y, ssmopt, Q, H, a, P); N=ssmopt.N; m=ssmopt.ssmodel.states; %---------------- 1. Recover parameters and prepare matrices ---------------- T=ssmopt.ssmodel.T; Z=ssmopt.ssmodel.Z; R=ssmopt.ssmodel.R; KF_out.Vmat=zeros(1,N); KF_out.Fmat=zeros(1,N); KF_out.Afilt=zeros(m,N); KF_out.Pfilt=zeros(m,m,N); KF_out.Kmat=zeros(m,N); KF_out.Lmat=zeros(m,m,N); LHmat=zeros(1,N); if ~isfield(ssmopt,'exactcheck'); ssmopt.exactcheck=1; % set exact filter initialization by default end; %---------------- 2. Set default starting values for a and P , if none provided ---------------- Pinf=eye(m); if nargin < 6 if ssmopt.exactcheck==1 P=zeros(m,m); else P=eye(m)*1000000000; end end if nargin < 5 a=[y(1); zeros(m-1,1)]; end KF_out.Afilt(:,1)= a; KF_out.Pfilt(:,:,1) = P; d=0; %---------------- 3. Exact Filtering ---------------- if ssmopt.exactcheck==1 evals=10; % number of time steps to evaluate until Pinf converges to zero KF_out.exact.F1=zeros(1,evals); KF_out.exact.F2=zeros(1,evals); KF_out.exact.L1=zeros(m,m,evals); KF_out.exact.Pinf=zeros(m,m,evals); KF_out.exact.Pinf(:,:,1)=Pinf; for i=1:evals if sum(sum(Pinf))<1e-20; d=i-1; % time point after which Pinf-->0, and after which we may start regular Kalman filter break; else if sum(Pinf*Z')>0 % Pinf is not singular F1=inv(Z*Pinf*Z'); F2=-F1*(Z*P*Z'+H)*F1; K=T*Pinf*Z'*F1; K1=T*(P*Z'*F1+Pinf*Z'*F2); L=TK*Z; L1=-K1*Z; P=T*Pinf*L1' + T*P*L' + R*Q*R'; Pinf=T*Pinf*L'; else F1=Z*P*Z'+H; F2=F1; K=T*P*Z'/F1; L=TK*Z; L1=L; P=T*P*L' + R*Q*R'; Pinf=T*Pinf*T'; end v=y(i) - Z*a; a=T*a+K*v; %save filtered estimates KF_out.Afilt(:,i+1)=a; KF_out.Pfilt(:,:,i+1)=P; KF_out.Vmat(i)=v; KF_out.Fmat(i)=F1; KF_out.Kmat(:,i)=K; KF_out.Lmat(:,:,i)=L; LHmat(i) = -0.5*(log(2*pi*F1) + v^2/F1); %save exact values for smoother KF_out.exact.F1(i)=F1; KF_out.exact.F2(i)=F2; KF_out.exact.L1(:,:,i)=L1; KF_out.exact.Pinf(:,:,i+1)=Pinf; end end end %---------------- 4. Regular Filtering ---------------- for i=d+1:N v=y(i) - Z*a; f=Z*P*Z' + H; K=T*P*Z'/f; L=TK*Z; a=T*a+K*v; P=T*P*L'+R*Q*R'; if i<N KF_out.Afilt(:,i+1)=a; KF_out.Pfilt(:,:,i+1)=P; end KF_out.Vmat(i)=v; KF_out.Fmat(i)=f; KF_out.Kmat(:,i)=K; KF_out.Lmat(:,:,i)=L; LHmat(i) = -0.5*(log(2*pi*f) + v^2/f); end %---------------- 5. Prepare output ---------------- KF_out.LH=-sum(LHmat); KF_out.Q=Q; KF_out.H=H; KF_out.exact.d=d; end
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

 function [Ksm_out, Kdism_out] = ksmy2(KF_out, ssmopt); [m,N]=size(KF_out.Afilt); meta=ssmopt.ssmodel.num_etas; %---------------- 1. Recover filtered matrices ---------------- Fmat=KF_out.Fmat; Vmat=KF_out.Vmat; Pfilt=KF_out.Pfilt; Afilt=KF_out.Afilt; Q=KF_out.Q; H=KF_out.H; %---------------- 2. Recover Model structure ---------------- Z=ssmopt.ssmodel.Z; T=ssmopt.ssmodel.T; R=ssmopt.ssmodel.R; Asm=zeros(m,N); Psm=zeros(m,m,N); rmat=zeros(m,N); Nmat=zeros(m,m,N); Eps=zeros(1,N); Eta=zeros(meta,N); Kmat=KF_out.Kmat; Lmat=KF_out.Lmat; if ~isfield(KF_out,'exact'); KF_out.exact.d=0; end d=KF_out.exact.d; if KF_out.exact.d>0 L1=KF_out.exact.L1; F1=KF_out.exact.F1; F2=KF_out.exact.F2; Pinf=KF_out.exact.Pinf; end %---------------- 3. Regular Smoothing for t=N..d+1 observations ---------------- for i=N:-1:d+1 r=Z'/Fmat(i)*Vmat(i) + Lmat(:,:,i)'*rmat(:,i); N=Z'/Fmat(i)*Z + Lmat(:,:,i)'*Nmat(:,:,i)*Lmat(:,:,i); Asm(:,i)=Afilt(:,i) + Pfilt(:,:,i)*r; Psm(:,:,i)=Pfilt(:,:,i)-Pfilt(:,:,i)*N*Pfilt(:,:,i); if i>1 rmat(:,i-1)=r; Nmat(:,:,i-1)=N; end if nargout>1 Eps(i)=H*(1/(Fmat(i))*Vmat(i)-Kmat(:,i)'*rmat(:,i)); Eta(:,i)=Q*R'*rmat(:,i); end end %---------------- 4. Exact Smoothing for t=d..1 observations ---------------- if KF_out.exact.d>0 r1=zeros(m,1); N1=zeros(m,m); N2=zeros(m,m); for i=d:-1:1 if sum(Pinf(:,:,i)*Z')>0 %cond(Pinf)<1e+12 % Pinf is not singular r1=Z'*F1(i)*Vmat(i) + Lmat(:,:,i)'*r1 + L1(:,:,i)'*rmat(:,i); N2=Z'*F2(i)*Z + Lmat(:,:,i)'*N2*Lmat(:,:,i) + Lmat(:,:,i)'*N1*L1(:,:,i) + L1(:,:,i)'*N1*Lmat(:,:,i) + L1(:,:,i)'*Nmat(:,:,i)*L1(:,:,i); N1=Z'*F1(i)*Z + Lmat(:,:,i)'*N1*Lmat(:,:,i) + L1(:,:,i)'*Nmat(:,:,i)*Lmat(:,:,i); r=Lmat(:,:,i)'*r1; N=Lmat(:,:,i)'*Nmat(:,:,i)*Lmat(:,:,i); if nargout>1 Eps(i)=-H*Kmat(:,i)'*rmat(:,i); Eta(:,i)=Q*R'*rmat(:,i); end else % Pinf is singular r1=T'*rmat(:,i); N2=T'*N2*T; N1=T'*N1*Lmat(:,:,i); r=Z'/(Fmat(i))*Vmat(i) + Lmat(:,:,i)'*rmat(:,i); N=Z'/(Fmat(i))*Z + Lmat(:,:,i)'*Nmat(:,:,i)*Lmat(:,:,i); if nargout>1 Eps(i)=H*(1/Fmat(i)*Vmat(i) - Kmat(:,i)'*rmat(:,i)); Eta(:,i)=Q*R'*rmat(:,i); end end if i>1 rmat(:,i-1)=r; Nmat(:,:,i-1)=N; end Asm(:,i)=Afilt(:,i) + Pfilt(:,:,i)*r + Pinf(:,:,i)*r1; Psm(:,:,i)=Pfilt(:,:,i)-Pfilt(:,:,i)*N*Pfilt(:,:,i) - (Pinf(:,:,i)*N1*Pfilt(:,:,i))' - Pinf(:,:,i)*N1*Pfilt(:,:,i) - Pinf(:,:,i)*N2*Pinf(:,:,i); end end %---------------- 5. Prepare output ---------------- Ksm_out.Asm=Asm; Ksm_out.Psm=Psm; Ksm_out.Kmat=Kmat; Ksm_out.Lmat=Lmat; Ksm_out.Nmat=Nmat; Ksm_out.rmat=rmat; Kdism_out.Eps=Eps; Kdism_out.Eta=Eta;
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

 function b=transform(b,howto,ssmopt); k=length(b); if strcmp(howto,'in') % in-transformation if ssmopt.trans==0 % no transformation b=b; end; if ssmopt.trans==1 % transformation to preserve the positiveness of variances b(1:k-1,:)=log(b(1:k-1,:)); b(k)=log(1/b(k)-1); end; else % out-transformation if ssmopt.trans==0 % no transformation b=b; end; if ssmopt.trans==1 b(1:k-1,:)=exp(b(1:k-1,:)); b(k)=1/(1+exp(b(k))); end; end
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

• Durbin、J。、およびKoopman、状態空間法によるSJ時系列分析。 オックスフォード：オックスフォード大学出版局、2001。
• Durbin、J。、およびKoopman、SJ「状態空間時系列分析のためのよりシンプルで効率的なシミュレーション」Biometrika vol。 89、第3号、2002年。
• Harvey、AC、およびJaeger、A。「トレンド除去、定型化された事実およびビジネスサイクル。」Journal of Applied Econometrics（8）、1993。
• Harvey、AC「予測、構造時系列モデル、およびカルマンフィルター」 ケンブリッジ：ケンブリッジ大学出版局、1989年。