diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m index 820144f46..1b2ad753d 100644 --- a/matlab/DsgeSmoother.m +++ b/matlab/DsgeSmoother.m @@ -155,7 +155,8 @@ end if kalman_algo == 1 || kalman_algo == 3 [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp] = missing_DiffuseKalmanSmootherH1_Z(ST, ... Z,R1,Q,H,Pinf,Pstar, ... - data1,nobs,np,smpl,kalman_tol,riccati_tol,data_index); + data1,nobs,np,smpl,data_index, ... + options_.nk,kalman_tol,options_.filter_decomposition); if isequal(alphahat,0) if kalman_algo == 1 kalman_algo = 2; @@ -180,7 +181,10 @@ if kalman_algo == 2 || kalman_algo == 4 end [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp] = missing_DiffuseKalmanSmootherH3_Z(ST, ... - Z,R1,Q,diag(H),Pinf,Pstar,data1,nobs,np,smpl,data_index); + Z,R1,Q,diag(H), ... + Pinf,Pstar,data1,nobs,np,smpl,data_index, ... + options_.nk,kalman_tol,... + options_.filter_decomposition); end if kalman_algo == 3 || kalman_algo == 4 @@ -192,8 +196,10 @@ if kalman_algo == 3 || kalman_algo == 4 for i=1:size(PK,4) PK(jnk,:,:,i) = QT*dynare_squeeze(PK(jnk,:,:,i))*QT'; end - for i=1:size(decomp,4) - decomp(jnk,:,:,i) = QT*dynare_squeeze(decomp(jnk,:,:,i)); + if options_.filter_decomposition + for i=1:size(decomp,4) + decomp(jnk,:,:,i) = QT*dynare_squeeze(decomp(jnk,:,:,i)); + end end end for i=1:size(P,4) diff --git a/matlab/missing_DiffuseKalmanSmootherH1_Z.m b/matlab/missing_DiffuseKalmanSmootherH1_Z.m index b2c644396..852a7f9b8 100644 --- a/matlab/missing_DiffuseKalmanSmootherH1_Z.m +++ b/matlab/missing_DiffuseKalmanSmootherH1_Z.m @@ -1,6 +1,6 @@ -function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp] = missing_DiffuseKalmanSmootherH1_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,kalman_tol,riccati_tol,data_index) +function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp] = missing_DiffuseKalmanSmootherH1_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,decomp_flag) -% function [alphahat,epsilonhat,etahat,a, aK] = DiffuseKalmanSmoother1(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl) +% function [alphahat,epsilonhat,etahat,a,aK,PK,decomp] = DiffuseKalmanSmoother1(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,decomp_flag) % Computes the diffuse kalman smoother without measurement error, in the case of a non-singular var-cov matrix % % INPUTS @@ -15,9 +15,10 @@ function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp] = missing_DiffuseKal % pp: number of observed variables % mm: number of state variables % smpl: sample size -% kalman_tol tolerance for singular matrix -% riccati_tol % data_index [cell] 1*smpl cell of column vectors of indices. +% nk number of forecasting periods +% kalman_tol tolerance for reciprocal condition number +% decomp_flag if true, compute filter decomposition % % OUTPUTS % alphahat: smoothed variables (a_{t|T}) @@ -60,12 +61,8 @@ function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp] = missing_DiffuseKal % new crit1 value for rank of Pinf % it is assured that P is symmetric - -global options_ - d = 0; decomp = []; -nk = options_.nk; spinf = size(Pinf1); spstar = size(Pstar1); v = zeros(pp,smpl); @@ -83,7 +80,6 @@ Kstar = zeros(mm,pp,smpl); P = zeros(mm,mm,smpl+1); Pstar = zeros(spstar(1),spstar(2),smpl+1); Pstar(:,:,1) = Pstar1; Pinf = zeros(spinf(1),spinf(2),smpl+1); Pinf(:,:,1) = Pinf1; -crit = options_.kalman_tol; crit1 = 1.e-8; steady = smpl; rr = size(Q,1); @@ -170,7 +166,7 @@ while notsteady && t 8 +if decomp_flag decomp = zeros(nk,mm,rr,smpl+nk); ZRQinv = inv(Z*QQ*Z'); for t = max(d,1):smpl diff --git a/matlab/missing_DiffuseKalmanSmootherH3_Z.m b/matlab/missing_DiffuseKalmanSmootherH3_Z.m index f836e3cb2..4d1668adf 100644 --- a/matlab/missing_DiffuseKalmanSmootherH3_Z.m +++ b/matlab/missing_DiffuseKalmanSmootherH3_Z.m @@ -1,5 +1,5 @@ -function [alphahat,epsilonhat,etahat,a,P,aK,PK,decomp] = missing_DiffuseKalmanSmootherH3_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index) -% function [alphahat,epsilonhat,etahat,a1,P,aK,PK,d,decomp_filt] = missing_DiffuseKalmanSmoother3_Z(T,Z,R,Q,Pinf1,Pstar1,Y,pp,mm,smpl) +function [alphahat,epsilonhat,etahat,a,P,aK,PK,decomp] = missing_DiffuseKalmanSmootherH3_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,decomp_flag) +% function [alphahat,epsilonhat,etahat,a1,P,aK,PK,d,decomp] = missing_DiffuseKalmanSmootherH3_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,decomp_flag) % Computes the diffuse kalman smoother without measurement error, in the case of a singular var-cov matrix. % Univariate treatment of multivariate time series. % @@ -16,7 +16,10 @@ function [alphahat,epsilonhat,etahat,a,P,aK,PK,decomp] = missing_DiffuseKalmanSm % mm: number of state variables % smpl: sample size % data_index [cell] 1*smpl cell of column vectors of indices. -% +% nk number of forecasting periods +% kalman_tol tolerance for zero divider +% decomp_flag if true, compute filter decomposition +% % OUTPUTS % alphahat: smoothed state variables (a_{t|T}) % epsilonhat: measurement errors @@ -55,20 +58,9 @@ function [alphahat,epsilonhat,etahat,a,P,aK,PK,decomp] = missing_DiffuseKalmanSm % Modified by M. Ratto % New output argument aK: 1-step to nk-stpe ahed predictions) % New input argument nk: max order of predictions in aK -% New option options_.diffuse_d where the DKF stops (common with -% diffuselikelihood3) -% New icc variable to count number of iterations for Finf steps -% Pstar % Pinf simmetric -% New termination of DKF iterations based on options_.diffuse_d -% Li also stored during DKF iterations !! -% some bugs corrected in the DKF part of the smoother (Z matrix and -% alphahat) - -global options_ d = 0; decomp = []; -nk = options_.nk; spinf = size(Pinf1); spstar = size(Pstar1); v = zeros(pp,smpl); @@ -76,28 +68,18 @@ a = zeros(mm,smpl); a1 = zeros(mm,smpl+1); aK = zeros(nk,mm,smpl+nk); -if isempty(options_.diffuse_d), - smpl_diff = 1; -else - smpl_diff=rank(Pinf1); -end - -Fstar = zeros(pp,smpl_diff); -Finf = zeros(pp,smpl_diff); +Fstar = zeros(pp,smpl); +Finf = zeros(pp,smpl); Fi = zeros(pp,smpl); Ki = zeros(mm,pp,smpl); -%Li = zeros(mm,mm,pp,smpl); -%Linf = zeros(mm,mm,pp,smpl_diff); -%L0 = zeros(mm,mm,pp,smpl_diff); -Kstar = zeros(mm,pp,smpl_diff); +Kstar = zeros(mm,pp,smpl); P = zeros(mm,mm,smpl+1); P1 = P; PK = zeros(nk,mm,mm,smpl+nk); -Pstar = zeros(spstar(1),spstar(2),smpl_diff+1); Pstar(:,:,1) = Pstar1; -Pinf = zeros(spinf(1),spinf(2),smpl_diff+1); Pinf(:,:,1) = Pinf1; +Pstar = zeros(spstar(1),spstar(2),smpl); Pstar(:,:,1) = Pstar1; +Pinf = zeros(spinf(1),spinf(2),smpl); Pinf(:,:,1) = Pinf1; Pstar1 = Pstar; Pinf1 = Pinf; -crit = options_.kalman_tol; crit1 = 1.e-6; steady = smpl; rr = size(Q,1); % number of structural shocks @@ -123,7 +105,7 @@ while newRank && t < smpl Fstar(i,t) = Zi*Pstar(:,:,t)*Zi' +H(i); Finf(i,t) = Zi*Pinf(:,:,t)*Zi'; Kstar(:,i,t) = Pstar(:,:,t)*Zi'; - if Finf(i,t) > crit && newRank + if Finf(i,t) > kalman_tol && newRank icc=icc+1; Kinf(:,i,t) = Pinf(:,:,t)*Zi'; Kinf_Finf = Kinf(:,i,t)/Finf(i,t); @@ -133,7 +115,7 @@ while newRank && t < smpl Kstar(:,i,t)*Kinf_Finf' - ... Kinf_Finf*Kstar(:,i,t)'; Pinf(:,:,t) = Pinf(:,:,t) - Kinf(:,i,t)*Kinf(:,i,t)'/Finf(i,t); - elseif Fstar(i,t) > crit + elseif Fstar(i,t) > kalman_tol a(:,t) = a(:,t) + Kstar(:,i,t)*v(i,t)/Fstar(i,t); Pstar(:,:,t) = Pstar(:,:,t) - Kstar(:,i,t)*Kstar(:,i,t)'/Fstar(i,t); end @@ -180,8 +162,7 @@ while notsteady && t crit - %Li(:,:,i,t) = eye(mm)-Ki(:,i,t)*Z(i,:)/Fi(i,t); + if Fi(i,t) > kalman_tol a(:,t) = a(:,t) + Ki(:,i,t)*v(i,t)/Fi(i,t); P(:,:,t) = P(:,:,t) - Ki(:,i,t)*Ki(:,i,t)'/Fi(i,t); end @@ -197,7 +178,7 @@ while notsteady && t crit +% $$$ if Fi_s(i) > kalman_tol % $$$ a(:,t) = a(:,t) + Ki_s(:,i)*v(i,t)/Fi_s(i); % $$$ end % $$$ end @@ -236,7 +217,7 @@ while t > d+1 t = t-1; di = flipud(data_index{t})'; for i = di - if Fi(i,t) > crit + if Fi(i,t) > kalman_tol ri = Z(i,:)'/Fi(i,t)*v(i,t)+ri-Ki(:,i,t)'*ri/Fi(i,t)*Z(i,:)'; end end @@ -252,13 +233,12 @@ if d for t = d:-1:1 di = flipud(data_index{t})'; for i = di - % if Finf(i,t) > crit && ~(t==d && i>options_.diffuse_d), % use of options_.diffuse_d to be sure of DKF termination - if Finf(i,t) > crit + if Finf(i,t) > kalman_tol r1(:,t) = Z(i,:)'*v(i,t)/Finf(i,t) + ... (Kinf(:,i,t)'*Fstar(i,t)/Finf(i,t)-Kstar(:,i,t)')*r0(:,t)/Finf(i,t)*Z(i,:)' + ... r1(:,t)-Kinf(:,i,t)'*r1(:,t)/Finf(i,t)*Z(i,:)'; r0(:,t) = r0(:,t)-Kinf(:,i,t)'*r0(:,t)/Finf(i,t)*Z(i,:)'; - elseif Fstar(i,t) > crit % step needed whe Finf == 0 + elseif Fstar(i,t) > kalman_tol % step needed whe Finf == 0 r0(:,t) = Z(i,:)'/Fstar(i,t)*v(i,t)+r0(:,t)-(Kstar(:,i,t)'*r0(:,t))/Fstar(i,t)*Z(i,:)'; end end @@ -272,14 +252,14 @@ if d end end -if nargout > 8 +if decomp_flag decomp = zeros(nk,mm,rr,smpl+nk); ZRQinv = inv(Z*QQ*Z'); for t = max(d,1):smpl ri_d = zeros(mm,1); di = flipud(data_index{t})'; for i = di - if Fi(i,t) > crit + if Fi(i,t) > kalman_tol ri_d = Z(i,:)'/Fi(i,t)*v(i,t)+ri_d-Ki(:,i,t)'*ri_d/Fi(i,t)*Z(i,:)'; end end diff --git a/tests/kalman_filter_smoother/gen_data.mod b/tests/kalman_filter_smoother/gen_data.mod index a6fb2f4f1..32c12a7b2 100644 --- a/tests/kalman_filter_smoother/gen_data.mod +++ b/tests/kalman_filter_smoother/gen_data.mod @@ -31,7 +31,7 @@ var e_y; stderr 0.05; var e_z; stderr 0.05; end; -stoch_simul(periods=2000,irf=0,simul_seed=7); +stoch_simul(periods=2000,irf=0); plot([w x y z]);