From 2f3767d31750cb19447e84eda32213e463157535 Mon Sep 17 00:00:00 2001 From: michel Date: Mon, 4 Feb 2008 17:52:16 +0000 Subject: [PATCH] v4: more changes related to new filter more variables returned by the smoother saved in oo_ DiffuseKalmanSmoother1_Z.m isn't finished yet git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1689 ac1d8469-bf42-47a9-8791-bf33cf982152 --- matlab/DiffuseKalmanSmoother1_Z.m | 39 +++++++++++++------ matlab/DiffuseKalmanSmoother3_Z.m | 62 +++++++++++++++++++++++++++---- matlab/DsgeSmoother.m | 56 ++++++++++++++++++++++------ matlab/dynare_estimation.m | 17 +++++++-- 4 files changed, 141 insertions(+), 33 deletions(-) diff --git a/matlab/DiffuseKalmanSmoother1_Z.m b/matlab/DiffuseKalmanSmoother1_Z.m index 1beeba5a2..1cda2fe21 100644 --- a/matlab/DiffuseKalmanSmoother1_Z.m +++ b/matlab/DiffuseKalmanSmoother1_Z.m @@ -1,4 +1,4 @@ -function [alphahat,etahat,a, aK] = DiffuseKalmanSmoother1_Z(T,Z,R,Q,Pinf1,Pstar1,Y,pp,mm,smpl) +function [alphahat,etahat,a,aK,P,PK,d] = DiffuseKalmanSmoother1_Z(T,Z,R,Q,Pinf1,Pstar1,Y,pp,mm,smpl) % function [alphahat,etahat,a, aK] = DiffuseKalmanSmoother1(T,Z,R,Q,Pinf1,Pstar1,Y,pp,mm,smpl) % Computes the diffuse kalman smoother without measurement error, in the case of a non-singular var-cov matrix @@ -20,7 +20,14 @@ function [alphahat,etahat,a, aK] = DiffuseKalmanSmoother1_Z(T,Z,R,Q,Pinf1,Pstar1 % etahat: smoothed shocks % a: matrix of one step ahead filtered state variables % aK: 3D array of k step ahead filtered state variables - +% (meaningless for periods 1:d) +% P: 3D array of one-step ahead forecast error variance +% matrices +% PK: 4D array of k-step ahead forecast error variance +% matrices (meaningless for periods 1:d) +% d: number of periods where filter remains in diffuse part +% (should be equal to the order of integration of the model) +% % SPECIAL REQUIREMENTS % See "Filtering and Smoothing of State Vector for Diffuse State Space % Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series @@ -44,7 +51,8 @@ spinf = size(Pinf1); spstar = size(Pstar1); v = zeros(pp,smpl); a = zeros(mm,smpl+1); -aK = zeros(nk,mm,smpl+1); +aK = zeros(nk,mm,smpl+nk); +PK = zeros(nk,mm,mm,smpl+nk); iF = zeros(pp,pp,smpl); Fstar = zeros(pp,pp,smpl); iFinf = zeros(pp,pp,smpl); @@ -77,6 +85,7 @@ while rank(Pinf(:,:,t+1),crit1) & t 7 + decomp = zeros(nk,mm,rr,smpl+nk); + ZRQinv = inv(Z*QQ*Z'); + for t = d:smpl + ri_d = zeros(mm,1); + for i=pp:-1:1 + if Fi(i,t) > crit + ri_d = Z(i,:)'/Fi(i,t)*v(i,t)+Li(:,:,i,t)'*ri_d; + end + end + + % calculate eta_tm1t + eta_tm1t = QRt*ri_d; + % calculate decomposition + Ttok = eye(mm,mm); + for h = 1:nk + for j=1:rr + eta=zeros(rr,1); + eta(j) = eta_tm1t(j); + decomp(h,:,j,t+h) = Ttok*P1(:,:,t)*Z'*ZRQinv*Z*R*eta; + end + Ttok = T*Ttok; + end + end +end \ No newline at end of file diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m index 45d86b91e..4bad064a4 100644 --- a/matlab/DsgeSmoother.m +++ b/matlab/DsgeSmoother.m @@ -1,4 +1,4 @@ -function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R] = DsgeSmoother(xparam1,gend,Y) +function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,d,decomp] = DsgeSmoother(xparam1,gend,Y) % Estimation of the smoothed variables and innovations. % % INPUTS @@ -14,9 +14,17 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R] = Dsge % o SteadyState [double] (m*1) vector specifying the steady state level of each endogenous variable. % o trend_coeff [double] (n*1) vector, parameters specifying the slope of the trend associated to each observed variable. % o aK [double] (K,n,T+K) array, k (k=1,...,K) steps ahead filtered (endogenous) variables. -% o T and R [double] Matrices defining the state equation (T is the (m*m) transition matrix). +% o T and R [double] Matrices defining the state equation (T is +% the (m*m) transition matrix). +% P: 3D array of one-step ahead forecast error variance +% matrices +% PK: 4D array of k-step ahead forecast error variance +% matrices (meaningless for periods 1:d) +% d: number of periods where filter remains in diffuse part +% (should be equal to the order of integration of the model) +% % ALGORITHM -% Metropolis-Hastings. +% Diffuse Kalman filter (Durbin and Koopman) % % SPECIAL REQUIREMENTS % None. @@ -27,8 +35,18 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R] = Dsge global bayestopt_ M_ oo_ estim_params_ options_ alphahat = []; - epsilonhat = []; etahat = []; + epsilonhat = []; + ahat = []; + SteadyState = []; + trend_coeff = []; + aK = []; + T = []; + R = []; + P = []; + PK = []; + d = []; + decomp = []; nobs = size(options_.varobs,1); smpl = size(Y,2); @@ -168,15 +186,31 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R] = Dsge end elseif options_.kalman_algo == 3 [alphahat,etahat,ahat,aK] = DiffuseKalmanSmoother3(T,R,Q,Pinf,Pstar,Y,trend,nobs,np,smpl,mf); - elseif options_.kalman_algo == 4 + elseif options_.kalman_algo == 4 | options_.kalman_algo == 5 data1 = Y - trend; - [alphahat,etahat,ahat,aK] = DiffuseKalmanSmoother1_Z(ST,Z,R1,Q,Pinf,Pstar,data1,nobs,np,smpl); - alphahat = QT*alphahat; - ahat = QT*ahat; - elseif options_.kalman_algo == 5 - data1 = Y - trend; - [alphahat,etahat,ahat,aK] = DiffuseKalmanSmoother3_Z(ST,Z,R1,Q,Pinf,Pstar,data1,nobs,np,smpl); + if options_.kalman_algo == 4 + [alphahat,etahat,ahat,P,aK,PK,d] = DiffuseKalmanSmoother1_Z(ST, ... + Z,R1,Q,Pinf,Pstar,data1,nobs,np,smpl); + else + [alphahat,etahat,ahat,P,aK,PK,d,decomp] = DiffuseKalmanSmoother3_Z(ST, ... + Z,R1,Q,Pinf,Pstar,data1,nobs,np,smpl); + end alphahat = QT*alphahat; ahat = QT*ahat; + if options_.nk > 0 + nk = options_.nk; + for jnk=1:nk + aK(jnk,:,:) = QT*squeeze(aK(jnk,:,:)); + for i=1:size(PK,4) + PK(jnk,:,:,i) = QT*squeeze(PK(jnk,:,:,i))*QT'; + end + for i=1:size(decomp,4) + decomp(jnk,:,:,i) = QT*squeeze(decomp(jnk,:,:,i)); + end + end + for i=1:size(P,4) + P(:,:,i) = QT*squeeze(P(:,:,i))*QT'; + end + end end end diff --git a/matlab/dynare_estimation.m b/matlab/dynare_estimation.m index 3df400182..5a8274ada 100644 --- a/matlab/dynare_estimation.m +++ b/matlab/dynare_estimation.m @@ -880,11 +880,20 @@ if (any(bayestopt_.pshape >0 ) & options_.mh_replic) | ... % return end -if ~((any(bayestopt_.pshape > 0) & options_.mh_replic) | (any(bayestopt_.pshape ... - > 0) & options_.load_mh_file)) | ~options_.smoother +if (~((any(bayestopt_.pshape > 0) & options_.mh_replic) | (any(bayestopt_.pshape ... + > 0) & options_.load_mh_file)) ... + | ~options_.smoother ) & M_.endo_nbr^2*gend < 1e6 % to be fixed %% ML estimation, or posterior mode without metropolis-hastings or metropolis without bayesian smooth variable - options_.lik_algo = 2; - [atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff] = DsgeSmoother(xparam1,gend,data); + [atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff,aK,T,R,P,PK,d,decomp] = DsgeSmoother(xparam1,gend,data); + oo_.Smoother.SteadyState = ys; + oo_.Smoother.TrendCoeffs = trend_coeff; + oo_.Smoother.integration_order = d; + oo_.Smoother.variance = P; + if options_.nk ~= 0 + oo_.FilteredVariablesKStepAhead = aK(options_.filter_step_ahead,:,:); + oo_.FilteredVariablesKStepAheadVariances = PK(options_.filter_step_ahead,:,:,:); + oo_.FilteredVariablesShockDecomposition = decomp(options_.filter_step_ahead,:,:,:); + end for i=1:M_.endo_nbr eval(['oo_.SmoothedVariables.' deblank(M_.endo_names(dr.order_var(i),:)) ' = atT(i,:)'';']); eval(['oo_.FilteredVariables.' deblank(M_.endo_names(dr.order_var(i),:)) ' = filtered_state_vector(i,:)'';']);