diff --git a/doc/dynare.texi b/doc/dynare.texi index d884995b2..330f920da 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -5919,7 +5919,7 @@ Use the Univariate Diffuse Kalman Filter Default value is @code{0}. In case of missing observations of single or all series, Dynare treats those missing values as unobserved states and uses the Kalman filter to infer their value (see e.g. @cite{Durbin and Koopman (2012), Ch. 4.10}) This procedure has the advantage of being capable of dealing with observations where the forecast error variance matrix becomes singular for some variable(s). If this happens, the respective observation enters with a weight of zero in the log-likelihood, i.e. this observation for the respective variable(s) is dropped -from the likelihood computations (for details see @cite{Durbin and Koopman (2012), Ch. 6.4 and 7.2.5}). If the use of a multivariate Kalman filter is specified and a +from the likelihood computations (for details see @cite{Durbin and Koopman (2012), Ch. 6.4 and 7.2.5} and @cite{Koopman and Durbin (2000)}). If the use of a multivariate Kalman filter is specified and a singularity is encountered, Dynare by default automatically switches to the univariate Kalman filter for this parameter draw. This behavior can be changed via the @ref{use_univariate_filters_if_singularity_is_detected} option. @@ -5954,12 +5954,16 @@ Triggers the computation k-step ahead filtered values, i.e. @math{E_{t}{y_{t+k}} @anchor{filter_decomposition} Triggers the computation of the shock decomposition of the above k-step ahead filtered values. Stores results in @code{oo_.FilteredVariablesShockDecomposition}. +@item smoothed_state_uncertainty +@anchor{smoothed_state_uncertainty} Triggers the computation of the variance of smoothed estimates, i.e. +@code{Var_T(y_t)}. Stores results in @code{oo_.Smoother.State_uncertainty}. @item diffuse_filter @anchor{diffuse_filter} Uses the diffuse Kalman filter (as described in @cite{Durbin and Koopman (2012)} and @cite{Koopman and Durbin -(2003)}) to estimate models with non-stationary observed variables. +(2003)} for the multivariate and @cite{Koopman and Durbin +(2000)} for the univariate filter) to estimate models with non-stationary observed variables. When @code{diffuse_filter} is used the @code{lik_init} option of @code{estimation} has no effect. @@ -6512,6 +6516,19 @@ two dimensions. The third dimension of the array provides the observation for which the forecast has been made. @end defvr +@defvr {MATLAB/Octave variable} oo_.Smoother.State_uncertainty +@anchor{oo_.Smoother.State_uncertainty} +Three-dimensional array set by the @code{estimation} command (if used with the +@code{smoother}) without Metropolis, +or by the @code{calib_smoother} command, if the @code{o_smoothed_state_uncertainty} option +has been requested. +Contains the series of covariance matrices for the state estimate given the full data +from the Kalman smoother. The @code{M_.endo_nbr} times @code{M_.endo_nbr} times +@code{T} array contains the variables in declaration order along the first +two dimensions. The third dimension of the array provides the +observation for which the smoothed estimate has been made. +@end defvr + @defvr {MATLAB/Octave variable} oo_.Smoother.SteadyState @anchor{oo_.Smoother.SteadyState} Variable set by the @code{estimation} command (if used with the @@ -14128,6 +14145,11 @@ Dynamics and Control}, 32(11), 3397--3414 @item Koop, Gary (2003), @i{Bayesian Econometrics}, John Wiley & Sons +@item +Koopman, S. J. and J. Durbin (2000): ``Fast Filtering and Smoothing for +Multivariate State Space Models,'' @i{Journal of Time +Series Analysis}, 21(3), 281--296 + @item Koopman, S. J. and J. Durbin (2003): ``Filtering and Smoothing of State Vector for Diffuse State Space Models,'' @i{Journal of Time diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m index 065ddd68c..743454590 100644 --- a/matlab/DsgeSmoother.m +++ b/matlab/DsgeSmoother.m @@ -1,4 +1,4 @@ -function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,decomp,trend_addition] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value) +function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,decomp,trend_addition,state_uncertainty] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value) % Estimation of the smoothed variables and innovations. % % INPUTS @@ -25,6 +25,8 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,de % o decomp (K*m*r*(T+K)) 4D array of shock decomposition of k-step ahead % filtered variables (decision-rule order) % o trend_addition [double] (n*T) pure trend component; stored in options_.varobs order +% o state_uncertainty [double] (K,K,T) array, storing the uncertainty +% about the smoothed state (decision-rule order) % % Notes: % m: number of endogenous variables (M_.endo_nbr) @@ -233,10 +235,10 @@ ST = T; R1 = R; if kalman_algo == 1 || kalman_algo == 3 - [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp] = missing_DiffuseKalmanSmootherH1_Z(ST, ... + [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty] = missing_DiffuseKalmanSmootherH1_Z(ST, ... Z,R1,Q,H,Pinf,Pstar, ... data1,vobs,np,smpl,data_index, ... - options_.nk,kalman_tol,diffuse_kalman_tol,options_.filter_decomposition); + options_.nk,kalman_tol,diffuse_kalman_tol,options_.filter_decomposition,options_.smoothed_state_uncertainty); if isinf(alphahat) if kalman_algo == 1 kalman_algo = 2; @@ -273,11 +275,11 @@ if kalman_algo == 2 || kalman_algo == 4 end end - [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp] = missing_DiffuseKalmanSmootherH3_Z(ST, ... + [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty] = missing_DiffuseKalmanSmootherH3_Z(ST, ... Z,R1,Q,diag(H), ... Pinf,Pstar,data1,vobs,np,smpl,data_index, ... options_.nk,kalman_tol,diffuse_kalman_tol, ... - options_.filter_decomposition); + options_.filter_decomposition,options_.smoothed_state_uncertainty); end @@ -299,4 +301,7 @@ if expanded_state_vector_for_univariate_filter && (kalman_algo == 2 || kalman_al if ~isempty(P) P = P(k,k,:); end + if ~isempty(state_uncertainty) + state_uncertainty = state_uncertainty(k,k,:); + end end diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index db0e1c650..ff9e86836 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -163,8 +163,8 @@ end if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && options_.mh_posterior_mode_estimation==0 if options_.smoother == 1 - [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend] = DsgeSmoother(xparam1,gend,transpose(data),data_index,missing_value); - [oo_]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend); + [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty] = DsgeSmoother(xparam1,gend,transpose(data),data_index,missing_value); + [oo_]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty); end %reset qz_criterium options_.qz_criterium=qz_criterium_old; @@ -483,8 +483,8 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha > 0) && options_.load_mh_file)) ... || ~options_.smoother ) && options_.partial_information == 0 % to be fixed %% ML estimation, or posterior mode without Metropolis-Hastings or Metropolis without Bayesian smoothes variables - [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend] = DsgeSmoother(xparam1,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state); - [oo_,yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend); + [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty] = DsgeSmoother(xparam1,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state); + [oo_,yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty); if ~options_.nograph, [nbplt,nr,nc,lr,lc,nstar] = pltorg(M_.exo_nbr); diff --git a/matlab/evaluate_smoother.m b/matlab/evaluate_smoother.m index 3c2a81f00..6e420d95a 100644 --- a/matlab/evaluate_smoother.m +++ b/matlab/evaluate_smoother.m @@ -97,9 +97,9 @@ if ischar(parameters) end end -[atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend] = ... +[atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty] = ... DsgeSmoother(parameters,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state); -[oo_]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend); +[oo_]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty); if nargout==2 Smoothed_variables_declaration_order_deviation_form=atT(oo_.dr.inv_order_var(bayestopt_.smoother_var_list),:); diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m index 37b581266..12273dd17 100644 --- a/matlab/global_initialization.m +++ b/matlab/global_initialization.m @@ -413,6 +413,7 @@ options_.bayesian_th_moments = 0; options_.diffuse_filter = 0; options_.filter_step_ahead = []; options_.filtered_vars = 0; +options_.smoothed_state_uncertainty = 0; options_.first_obs = NaN; options_.nobs = NaN; options_.kalman_algo = 0; diff --git a/matlab/imcforecast.m b/matlab/imcforecast.m index 99ee00219..1a9d1607e 100644 --- a/matlab/imcforecast.m +++ b/matlab/imcforecast.m @@ -123,7 +123,9 @@ if estimated_model %store qz_criterium qz_criterium_old=options_.qz_criterium; options_=select_qz_criterium_value(options_); + options_state_uncertainty_old=options_.state_uncertainty; [atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff,aK,T,R,P,PK,decomp,trend_addition] = DsgeSmoother(xparam,gend,data,data_index,missing_value); + options_.state_uncertainty=options_state_uncertainty_old; %get constant part if options_.noconstant constant = zeros(size(ys,1),options_cond_fcst.periods+1); diff --git a/matlab/kalman/likelihood/kalman_filter_d.m b/matlab/kalman/likelihood/kalman_filter_d.m index a1e3281d7..a10cf6f4d 100644 --- a/matlab/kalman/likelihood/kalman_filter_d.m +++ b/matlab/kalman/likelihood/kalman_filter_d.m @@ -87,7 +87,7 @@ while rank(Pinf,diffuse_kalman_tol) && (t<=last) else iFstar = inv(Fstar); dFstar = det(Fstar); - Kstar = Pstar*Z'*iFstar; %(5.15) of DK (2012) with Kstar=K^(0)*T^{-1} + Kstar = Pstar*Z'*iFstar; %(5.15) of DK (2012) with Kstar=T^{-1}*K^(0) dlik(s)= log(dFstar) + v'*iFstar*v; %set w_t to bottom case in bottom equation page 172, DK (2012) Pinf = T*Pinf*transpose(T); % (5.16) DK (2012) Pstar = T*(Pstar-Pstar*Z'*Kstar')*T'+QQ; % (5.17) DK (2012) @@ -96,8 +96,8 @@ while rank(Pinf,diffuse_kalman_tol) && (t<=last) end else %F_{\infty,t} positive definite %To compare to DK (2012), this block makes use of the following transformation - %Kstar=K^{(1)}*T^{-1}=M_{*}*F^{(1)}+M_{\infty}*F^{(2)} - % =P_{*}*Z'*F^{(1)}+P_{\infty}*Z'*((-1)*(F_{\infty}^{-1})*F_{*}*(F_{\infty}^{-1})) + %Kstar=T^{-1}*K^{(1)}=M_{*}*F^{(1)}+M_{\infty}*F^{(2)} + % =P_{*}*Z'*F^{(1)}+P_{\infty}*Z'*((-1)*(-F_{\infty}^{-1})*F_{*}*(F_{\infty}^{-1})) % =[P_{*}*Z'-Kinf*F_{*})]*F^{(1)} %Make use of L^{0}'=(T-K^{(0)}*Z)'=(T-T*M_{\infty}*F^{(1)}*Z)' % =(T-T*P_{\infty*Z'*F^{(1)}*Z)'=(T-T*Kinf*Z)' @@ -106,9 +106,9 @@ while rank(Pinf,diffuse_kalman_tol) && (t<=last) % =T*[(P_{\infty}*(-K^{(1)*Z}))+P_{*}*(I-Z'*Kinf')*T'+RQR] dlik(s)= log(det(Finf)); %set w_t to top case in bottom equation page 172, DK (2012) iFinf = inv(Finf); - Kinf = Pinf*Z'*iFinf; %define Kinf=K_0*T^{-1} with M_{\infty}=Pinf*Z' + Kinf = Pinf*Z'*iFinf; %define Kinf=T^{-1}*K_0 with M_{\infty}=Pinf*Z' Fstar = Z*Pstar*Z' + H; %(5.7) DK(2012) - Kstar = (Pstar*Z'-Kinf*Fstar)*iFinf; %(5.12) DK(2012) + Kstar = (Pstar*Z'-Kinf*Fstar)*iFinf; %(5.12) DK(2012); note that there is a typo in DK (2003) with "+ Kinf" instead of "- Kinf", but it is correct in their appendix Pstar = T*(Pstar-Pstar*Z'*Kinf'-Pinf*Z'*Kstar')*T'+QQ; %(5.14) DK(2012) Pinf = T*(Pinf-Pinf*Z'*Kinf')*T'; %(5.14) DK(2012) a = T*(a+Kinf*v); %(5.13) DK(2012) diff --git a/matlab/kalman/likelihood/missing_observations_kalman_filter_d.m b/matlab/kalman/likelihood/missing_observations_kalman_filter_d.m index acfa5563c..a497d6e53 100644 --- a/matlab/kalman/likelihood/missing_observations_kalman_filter_d.m +++ b/matlab/kalman/likelihood/missing_observations_kalman_filter_d.m @@ -106,10 +106,10 @@ while rank(Pinf,diffuse_kalman_tol) && (t<=last) else iFstar = inv(Fstar); dFstar = det(Fstar); - Kstar = Pstar*ZZ'*iFstar; %(5.15) of DK (2012) with Kstar=K^(0)*T^{-1} + Kstar = Pstar*ZZ'*iFstar; %(5.15) of DK (2012) with Kstar=T^{-1}*K^(0) dlik(s) = log(dFstar) + v'*iFstar*v + length(d_index)*log(2*pi); %set w_t to bottom case in bottom equation page 172, DK (2012) Pinf = T*Pinf*transpose(T); % (5.16) DK (2012) - Pstar = T*(Pstar-Pstar*ZZ'*Kstar')*T'+QQ; % (5.17) DK (2012) + Pstar = T*(Pstar-Pstar*ZZ'*Kstar')*T'+QQ; % (5.17) DK (2012) with L_0 plugged in a = T*(a+Kstar*v); % (5.13) DK (2012) end end @@ -119,7 +119,7 @@ while rank(Pinf,diffuse_kalman_tol) && (t<=last) Kinf = Pinf*ZZ'*iFinf; %see notes in kalman_filter_d.m for details of computations Fstar = ZZ*Pstar*ZZ' + H(d_index,d_index); %(5.7) DK(2012) - Kstar = (Pstar*ZZ'-Kinf*Fstar)*iFinf; %(5.12) DK(2012) + Kstar = (Pstar*ZZ'-Kinf*Fstar)*iFinf; %(5.12) DK(2012); note that there is a typo in DK (2003) with "+ Kinf" instead of "- Kinf", but it is correct in their appendix Pstar = T*(Pstar-Pstar*ZZ'*Kinf'-Pinf*ZZ'*Kstar')*T'+QQ; %(5.14) DK(2012) Pinf = T*(Pinf-Pinf*ZZ'*Kinf')*T'; %(5.14) DK(2012) a = T*(a+Kinf*v); %(5.13) DK(2012) diff --git a/matlab/missing_DiffuseKalmanSmootherH1_Z.m b/matlab/missing_DiffuseKalmanSmootherH1_Z.m index 80b5c6efa..a19d3fff6 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,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag) +function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp,V] = missing_DiffuseKalmanSmootherH1_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag,state_uncertainty_flag) -% 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,diffuse_kalman_tol,decomp_flag) +% 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,diffuse_kalman_tol,decomp_flag,state_uncertainty_flag) % Computes the diffuse kalman smoother without measurement error, in the case of a non-singular var-cov matrix. % % INPUTS @@ -20,6 +20,8 @@ function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp] = missing_DiffuseKal % kalman_tol tolerance for reciprocal condition number % diffuse_kalman_tol tolerance for reciprocal condition number (for Finf) and the rank of Pinf % decomp_flag if true, compute filter decomposition +% state_uncertainty_flag if true, compute uncertainty about smoothed +% state estimate % % OUTPUTS % alphahat: smoothed variables (a_{t|T}) @@ -33,6 +35,7 @@ function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp] = missing_DiffuseKal % PK: 4D array of k-step ahead forecast error variance % matrices (meaningless for periods 1:d) % decomp: decomposition of the effect of shocks on filtered values +% V: 3D array of state uncertainty matrices % % Notes: % Outputs are stored in decision-rule order, i.e. to get variables in order of declaration @@ -43,6 +46,8 @@ function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp] = missing_DiffuseKal % See "Filtering and Smoothing of State Vector for Diffuse State Space % Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series % Analysis, vol. 24(1), pp. 85-98). +% Durbin/Koopman (2012): "Time Series Analysis by State Space Methods", Oxford University Press, +% Second Edition, Ch. 5 % Copyright (C) 2004-2016 Dynare Team % @@ -78,76 +83,92 @@ aK = zeros(nk,mm,smpl+nk); PK = zeros(nk,mm,mm,smpl+nk); iF = zeros(pp,pp,smpl); Fstar = zeros(pp,pp,smpl); +iFstar = zeros(pp,pp,smpl); iFinf = zeros(pp,pp,smpl); K = zeros(mm,pp,smpl); L = zeros(mm,mm,smpl); Linf = zeros(mm,mm,smpl); +Lstar = zeros(mm,mm,smpl); Kstar = zeros(mm,pp,smpl); +Kinf = 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; -steady = smpl; +Pstar = zeros(spstar(1),spstar(2),smpl+1); +Pstar(:,:,1) = Pstar1; +Pinf = zeros(spinf(1),spinf(2),smpl+1); +Pinf(:,:,1) = Pinf1; rr = size(Q,1); QQ = R*Q*transpose(R); QRt = Q*transpose(R); alphahat = zeros(mm,smpl); etahat = zeros(rr,smpl); -epsilonhat = zeros(rr,smpl); +epsilonhat = zeros(rr,smpl); r = zeros(mm,smpl+1); +Finf_singular = zeros(1,smpl); +if state_uncertainty_flag + V = zeros(mm,mm,smpl); + N = zeros(mm,mm,smpl+1); +else + V=[]; +end t = 0; while rank(Pinf(:,:,t+1),diffuse_kalman_tol) && td+1 t = t-1; di = data_index{t}; if isempty(di) - r(:,t) = L(:,:,t)'*r(:,t+1); + % in this case, L is simply T due to Z=0, so that DK (2012), eq. 4.93 obtains + r(:,t) = L(:,:,t)'*r(:,t+1); %compute r_{t-1}, DK (2012), eq. 4.38 with Z=0 + if state_uncertainty_flag + N(:,:,t)=L(:,:,t)'*N(:,:,t+1)*L(:,:,t); %compute N_{t-1}, DK (2012), eq. 4.42 with Z=0 + end else ZZ = Z(di,:); - r(:,t) = ZZ'*iF(di,di,t)*v(di,t) + L(:,:,t)'*r(:,t+1); + r(:,t) = ZZ'*iF(di,di,t)*v(di,t) + L(:,:,t)'*r(:,t+1); %compute r_{t-1}, DK (2012), eq. 4.38 + if state_uncertainty_flag + N(:,:,t)=ZZ'*iF(di,di,t)*ZZ+L(:,:,t)'*N(:,:,t+1)*L(:,:,t); %compute N_{t-1}, DK (2012), eq. 4.42 + end + end + alphahat(:,t) = a(:,t) + P(:,:,t)*r(:,t); %DK (2012), eq. 4.35 + etahat(:,t) = QRt*r(:,t); %DK (2012), eq. 4.63 + if state_uncertainty_flag + V(:,:,t) = P(:,:,t)-P(:,:,t)*N(:,:,t)*P(:,:,t); %DK (2012), eq. 4.43 end - alphahat(:,t) = a(:,t) + P(:,:,t)*r(:,t); - etahat(:,t) = QRt*r(:,t); end -if d +if d %diffuse periods + % initialize r_d^(0) and r_d^(1) as below DK (2012), eq. 5.23 r0 = zeros(mm,d+1); - r0(:,d+1) = r(:,d+1); - r1 = zeros(mm,d+1); + r0(:,d+1) = r(:,d+1); %set r0_{d}, i.e. shifted by one period + r1 = zeros(mm,d+1); %set r1_{d}, i.e. shifted by one period + if state_uncertainty_flag + %N_0 at (d+1) is N(d+1), so we can use N for continuing and storing N_0-recursion + N_1=zeros(mm,mm,d+1); %set N_1_{d}=0, i.e. shifted by one period, below DK (2012), eq. 5.26 + N_2=zeros(mm,mm,d+1); %set N_2_{d}=0, i.e. shifted by one period, below DK (2012), eq. 5.26 + end for t = d:-1:1 - r0(:,t) = Linf(:,:,t)'*r0(:,t+1); di = data_index{t}; if isempty(di) r1(:,t) = Linf(:,:,t)'*r1(:,t+1); else - r1(:,t) = Z(di,:)'*(iFinf(di,di,t)*v(di,t)-Kstar(:,di,t)'*r0(:,t+1)) ... - + Linf(:,:,t)'*r1(:,t+1); + if ~Finf_singular(1,t) + r0(:,t) = Linf(:,:,t)'*r0(:,t+1); % DK (2012), eq. 5.21 where L^(0) is named Linf + r1(:,t) = Z(di,:)'*(iFinf(di,di,t)*v(di,t)-Kstar(:,di,t)'*T'*r0(:,t+1)) ... + + Linf(:,:,t)'*r1(:,t+1); % DK (2012), eq. 5.21, noting that i) F^(1)=(F^Inf)^(-1)(see 5.10), ii) where L^(0) is named Linf, and iii) Kstar=T^{-1}*K^(1) + if state_uncertainty_flag + L_1=(-T*Kstar(:,di,t)*Z(di,:)); % noting that Kstar=T^{-1}*K^(1) + N(:,:,t)=Linf(:,:,t)'*N(:,:,t+1)*Linf(:,:,t); % DK (2012), eq. 5.19, noting that L^(0) is named Linf + N_1(:,:,t)=Z(di,:)'*iFinf(di,di,t)*Z(di,:)+Linf(:,:,t)'*N_1(:,:,t+1)*Linf(:,:,t)... + +L_1'*N(:,:,t+1)*Linf(:,:,t); % DK (2012), eq. 5.29; note that, compared to DK (2003) this drops the term (L_1'*N(:,:,t+1)*Linf(:,:,t))' in the recursion due to it entering premultiplied by Pinf when computing V, and Pinf*Linf'*N=0 + N_2(:,:,t)=Z(di,:)'*(-iFinf(di,di,t)*Fstar(di,di,t)*iFinf(di,di,t))*Z(di,:) ... + + Linf(:,:,t)'*N_2(:,:,t+1)*Linf(:,:,t)... + + Linf(:,:,t)'*N_1(:,:,t+1)*L_1... + + L_1'*N_1(:,:,t+1)'*Linf(:,:,t)... + + L_1'*N(:,:,t+1)*L_1; % DK (2012), eq. 5.29 + end + else + r0(:,t) = Z(di,:)'*iFstar(di,di,t)*v(di,t)-Lstar(:,di,t)'*r0(:,t+1); % DK (2003), eq. (14) + r1(:,t) = T'*r1(:,t+1); % DK (2003), eq. (14) + if state_uncertainty_flag + N(:,:,t)=Z(di,:)'*iFstar(di,di,t)*Z(di,:)... + +Lstar(:,:,t)'*N(:,:,t+1)*Lstar(:,:,t); % DK (2003), eq. (14) + N_1(:,:,t)=T'*N_1(:,:,t+1)*Lstar(:,:,t); % DK (2003), eq. (14) + N_2(:,:,t)=T'*N_2(:,:,t+1)*T'; % DK (2003), eq. (14) + end + end + end + alphahat(:,t) = a(:,t) + Pstar(:,:,t)*r0(:,t) + Pinf(:,:,t)*r1(:,t); % DK (2012), eq. 5.23 + etahat(:,t) = QRt*r0(:,t); % DK (2012), p. 135 + if state_uncertainty_flag + V(:,:,t)=Pstar(:,:,t)-Pstar(:,:,t)*N(:,:,t)*Pstar(:,:,t)... + -(Pinf(:,:,t)*N_1(:,:,t)*Pstar(:,:,t))'... + - Pinf(:,:,t)*N_1(:,:,t)*Pstar(:,:,t)... + - Pinf(:,:,t)*N_2(:,:,t)*Pinf(:,:,t); % DK (2012), eq. 5.30 end - alphahat(:,t) = a(:,t) + Pstar(:,:,t)*r0(:,t) + Pinf(:,:,t)*r1(:,t); - etahat(:,t) = QRt*r0(:,t); end end @@ -260,7 +327,6 @@ if decomp_flag eta_tm1t = QRt*Z(di,:)'*iF(di,di,t)*v(di,t); AAA = P(:,:,t)*Z(di,:)'*ZRQinv(di,di)*bsxfun(@times,Z(di,:)*R,eta_tm1t'); % calculate decomposition - Ttok = eye(mm,mm); decomp(1,:,:,t+1) = AAA; for h = 2:nk AAA = T*AAA; diff --git a/matlab/missing_DiffuseKalmanSmootherH3_Z.m b/matlab/missing_DiffuseKalmanSmootherH3_Z.m index 489d2be5f..598f1c05a 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,nk,kalman_tol,diffuse_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) +function [alphahat,epsilonhat,etahat,a,P,aK,PK,decomp,V] = missing_DiffuseKalmanSmootherH3_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag,state_uncertainty_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,state_uncertainty_flag) % Computes the diffuse kalman smoother in the case of a singular var-cov matrix. % Univariate treatment of multivariate time series. % @@ -20,6 +20,8 @@ function [alphahat,epsilonhat,etahat,a,P,aK,PK,decomp] = missing_DiffuseKalmanSm % kalman_tol tolerance for zero divider % diffuse_kalman_tol tolerance for zero divider % decomp_flag if true, compute filter decomposition +% state_uncertainty_flag if true, compute uncertainty about smoothed +% state estimate % % OUTPUTS % alphahat: smoothed state variables (a_{t|T}) @@ -33,12 +35,22 @@ function [alphahat,epsilonhat,etahat,a,P,aK,PK,decomp] = missing_DiffuseKalmanSm % PK: 4D array of k-step ahead forecast error variance % matrices (meaningless for periods 1:d) % decomp: decomposition of the effect of shocks on filtered values +% V: 3D array of state uncertainty matrices % % Notes: % Outputs are stored in decision-rule order, i.e. to get variables in order of declaration % as in M_.endo_names, ones needs code along the lines of: % variables_declaration_order(dr.order_var,:) = alphahat +% +% Algorithm: % +% Uses the univariate filter as described in Durbin/Koopman (2012): "Time +% Series Analysis by State Space Methods", Oxford University Press, +% Second Edition, Ch. 6.4 + 7.2.5 +% and +% Koopman/Durbin (2000): "Fast Filtering and Smoothing for Multivariatze State Space +% Models", in Journal of Time Series Analysis, vol. 21(3), pp. 281-296. +% % 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 @@ -84,21 +96,29 @@ Finf = zeros(pp,smpl); Fi = zeros(pp,smpl); Ki = zeros(mm,pp,smpl); Kstar = zeros(mm,pp,smpl); +Kinf = zeros(spstar(1),pp,smpl); P = zeros(mm,mm,smpl+1); P1 = P; PK = zeros(nk,mm,mm,smpl+nk); -Pstar = zeros(spstar(1),spstar(2),smpl); Pstar(:,:,1) = Pstar1; -Pinf = zeros(spinf(1),spinf(2),smpl); 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; -steady = smpl; rr = size(Q,1); % number of structural shocks QQ = R*Q*transpose(R); -QRt = Q*transpose(R); +QRt = Q*transpose(R); alphahat = zeros(mm,smpl); etahat = zeros(rr,smpl); epsilonhat = zeros(rr,smpl); r = zeros(mm,smpl); +if state_uncertainty_flag + V = zeros(mm,mm,smpl); + N = zeros(mm,mm,smpl); +else + V=[]; +end t = 0; icc=0; @@ -111,23 +131,28 @@ while newRank && t < smpl di = data_index{t}'; for i=di Zi = Z(i,:); - v(i,t) = Y(i,t)-Zi*a(:,t); - 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) > diffuse_kalman_tol && newRank + v(i,t) = Y(i,t)-Zi*a(:,t); % nu_{t,i} in 6.13 in DK (2012) + Fstar(i,t) = Zi*Pstar(:,:,t)*Zi' +H(i); % F_{*,t} in 5.7 in DK (2012), relies on H being diagonal + Finf(i,t) = Zi*Pinf(:,:,t)*Zi'; % F_{\infty,t} in 5.7 in DK (2012) + Kstar(:,i,t) = Pstar(:,:,t)*Zi'; % KD (2000), eq. (15) + if Finf(i,t) > diffuse_kalman_tol && newRank % F_{\infty,t,i} = 0, use upper part of bracket on p. 175 DK (2012) for w_{t,i} icc=icc+1; - Kinf(:,i,t) = Pinf(:,:,t)*Zi'; + Kinf(:,i,t) = Pinf(:,:,t)*Zi'; % KD (2000), eq. (15) Kinf_Finf = Kinf(:,i,t)/Finf(i,t); - a(:,t) = a(:,t) + Kinf_Finf*v(i,t); + a(:,t) = a(:,t) + Kinf_Finf*v(i,t); % KD (2000), eq. (16) Pstar(:,:,t) = Pstar(:,:,t) + ... Kinf(:,i,t)*Kinf_Finf'*(Fstar(i,t)/Finf(i,t)) - ... Kstar(:,i,t)*Kinf_Finf' - ... - Kinf_Finf*Kstar(:,i,t)'; - Pinf(:,:,t) = Pinf(:,:,t) - Kinf(:,i,t)*Kinf(:,i,t)'/Finf(i,t); + Kinf_Finf*Kstar(:,i,t)'; % KD (2000), eq. (16) + Pinf(:,:,t) = Pinf(:,:,t) - Kinf(:,i,t)*Kinf(:,i,t)'/Finf(i,t); % KD (2000), eq. (16) 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); + a(:,t) = a(:,t) + Kstar(:,i,t)*v(i,t)/Fstar(i,t); % KD (2000), eq. (17) + Pstar(:,:,t) = Pstar(:,:,t) - Kstar(:,i,t)*Kstar(:,i,t)'/Fstar(i,t); % KD (2000), eq. (17) + % Pinf is passed through unaltered, see eq. (17) of + % Koopman/Durbin (2000) + else + % do nothing as a_{t,i+1}=a_{t,i} and P_{t,i+1}=P_{t,i}, see + % p. 157, DK (2012) end end if newRank @@ -142,8 +167,7 @@ while newRank && t < smpl end Pstar(:,:,t+1) = T*Pstar(:,:,t)*T'+ QQ; Pinf(:,:,t+1) = T*Pinf(:,:,t)*T'; - P0=Pinf(:,:,t+1); - if newRank, + if newRank newRank = rank(Pinf(:,:,t+1),diffuse_kalman_tol); end if oldRank ~= newRank @@ -169,15 +193,18 @@ while notsteady && 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); + a(:,t) = a(:,t) + Ki(:,i,t)*v(i,t)/Fi(i,t); %filtering according to (6.13) in DK (2012) + P(:,:,t) = P(:,:,t) - Ki(:,i,t)*Ki(:,i,t)'/Fi(i,t); %filtering according to (6.13) in DK (2012) + else + % do nothing as a_{t,i+1}=a_{t,i} and P_{t,i+1}=P_{t,i}, see + % p. 157, DK (2012) end end - a1(:,t+1) = T*a(:,t); + a1(:,t+1) = T*a(:,t); %transition according to (6.14) in DK (2012) Pf = P(:,:,t); aK(1,:,t+1) = a1(:,t+1); for jnk=1:nk @@ -187,7 +214,7 @@ while notsteady && t d+1 t = t-1; di = flipud(data_index{t})'; for i = di 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,:)'; + Li = eye(mm)-Ki(:,i,t)*Z(i,:)/Fi(i,t); + ri = Z(i,:)'/Fi(i,t)*v(i,t)+Li'*ri; % DK (2012), 6.15, equation for r_{t,i-1} + if state_uncertainty_flag + Ni = Z(i,:)'/Fi(i,t)*Z(i,:)+Li'*Ni*Li; % KD (2000), eq. (23) + end end end - r(:,t) = ri; + r(:,t) = ri; % DK (2012), below 6.15, r_{t-1}=r_{t,0} alphahat(:,t) = a1(:,t) + P1(:,:,t)*r(:,t); etahat(:,t) = QRt*r(:,t); - ri = T'*ri; + ri = T'*ri; % KD (2003), eq. (23), equation for r_{t-1,p_{t-1}} + if state_uncertainty_flag + N(:,:,t) = Ni; % DK (2012), below 6.15, N_{t-1}=N_{t,0} + V(:,:,t) = P1(:,:,t)-P1(:,:,t)*N(:,:,t)*P1(:,:,t); % KD (2000), eq. (7) with N_{t-1} stored in N(:,:,t) + Ni = T'*Ni*T; % KD (2000), eq. (23), equation for N_{t-1,p_{t-1}} + end end if d r0 = zeros(mm,d); r0(:,d) = ri; r1 = zeros(mm,d); + if state_uncertainty_flag + %N_0 at (d+1) is N(d+1), so we can use N for continuing and storing N_0-recursion + N_0=zeros(mm,mm,d); %set N_1_{d}=0, below KD (2000), eq. (24) + N_0(:,:,d) = Ni; + N_1=zeros(mm,mm,d); %set N_1_{d}=0, below KD (2000), eq. (24) + N_2=zeros(mm,mm,d); %set N_2_{d}=0, below KD (2000), eq. (24) + end for t = d:-1:1 di = flipud(data_index{t})'; for i = di - if Finf(i,t) > diffuse_kalman_tol + if Finf(i,t) > diffuse_kalman_tol + % recursions need to be from highest to lowest term in order to not + % overwrite lower terms still needed in thisstep 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,:)'; + r1(:,t)-Kinf(:,i,t)'*r1(:,t)/Finf(i,t)*Z(i,:)'; % KD (2000), eq. (25) for r_1 + r0(:,t) = r0(:,t)-Kinf(:,i,t)'*r0(:,t)/Finf(i,t)*Z(i,:)'; % KD (2000), eq. (25) for r_0 + if state_uncertainty_flag + Linf = eye(mm) - Kinf(:,i,t)'/Finf(i,t); + L0 = (Kinf(:,i,t)*(Fstar(i,t)/Finf(i,t))-Kstar(:,i,t))*Z(i,:)/Finf(i,t); + N_2(:,:,t)=Z(i,:)'/Finf(i,t)^2*Z(i,:)*Fstar(i,t) ... + + Linf'*N_2(:,:,t)*Linf... + + Linf'*N_1(:,:,t)*L0... + + L0'*N_1(:,:,t)'*Linf... + + L0'*N_0(:,:,t)*L0; % DK (2012), eq. 5.29 + N_1(:,:,t)=Z(i,:)'/Finf(i,t)*Z(i,:)+Linf'*N_1(:,:,t)*Linf... + +L0'*N_0(:,:,t)*Linf; % DK (2012), eq. 5.29; note that, compared to DK (2003) this drops the term (L_1'*N(:,:,t+1)*Linf(:,:,t))' in the recursion due to it entering premultiplied by Pinf when computing V, and Pinf*Linf'*N=0 + N_0(:,:,t)=Linf'*N_0(:,:,t)*Linf; % DK (2012), eq. 5.19, noting that L^(0) is named Linf + end 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,:)'; + L_i=eye(mm) - Kstar(:,i,t)*Z(i,:)*Fstar(i,t); + r0(:,t) = Z(i,:)'/Fstar(i,t)*v(i,t)+L_i'*r0(:,t); % propagate r0 and keep r1 fixed + if state_uncertainty_flag + N_0(:,:,t)=Z(i,:)'/Fstar(i,t)*Z(i,:)+L_i'*N_0(:,:,t)*L_i; % propagate N_0 and keep N_1 and N_2 fixed + end end end - alphahat(:,t) = a1(:,t) + Pstar1(:,:,t)*r0(:,t) + Pinf1(:,:,t)*r1(:,t); + alphahat(:,t) = a1(:,t) + Pstar1(:,:,t)*r0(:,t) + Pinf1(:,:,t)*r1(:,t); % KD (2000), eq. (26) r(:,t) = r0(:,t); - etahat(:,t) = QRt*r(:,t); + etahat(:,t) = QRt*r(:,t); % KD (2000), eq. (27) + if state_uncertainty_flag + V(:,:,t)=Pstar(:,:,t)-Pstar(:,:,t)*N_0(:,:,t)*Pstar(:,:,t)... + -(Pinf(:,:,t)*N_1(:,:,t)*Pstar(:,:,t))'... + - Pinf(:,:,t)*N_1(:,:,t)*Pstar(:,:,t)... + - Pinf(:,:,t)*N_2(:,:,t)*Pinf(:,:,t); % DK (2012), eq. 5.30 + end if t > 1 - r0(:,t-1) = T'*r0(:,t); - r1(:,t-1) = T'*r1(:,t); + r0(:,t-1) = T'*r0(:,t); % KD (2000), below eq. (25) r_{t-1,p_{t-1}}=T'*r_{t,0} + r1(:,t-1) = T'*r1(:,t); % KD (2000), below eq. (25) r_{t-1,p_{t-1}}=T'*r_{t,0} + if state_uncertainty_flag + N_0(:,:,t-1)= T'*N_0(:,t)*T; % KD (2000), below eq. (25) N_{t-1,p_{t-1}}=T'*N_{t,0}*T + N_1(:,:,t-1)= T'*N_1(:,t)*T; % KD (2000), below eq. (25) N^1_{t-1,p_{t-1}}=T'*N^1_{t,0}*T + N_2(:,:,t-1)= T'*N_2(:,t)*T; % KD (2000), below eq. (25) N^2_{t-1,p_{t-1}}=T'*N^2_{t,0}*T + end end end end diff --git a/matlab/pm3.m b/matlab/pm3.m index 3c6857c62..21a486162 100644 --- a/matlab/pm3.m +++ b/matlab/pm3.m @@ -85,6 +85,7 @@ fprintf(['Estimation::mcmc: ' tit1 '\n']); k = 0; filter_step_ahead_indicator=0; filter_covar_indicator=0; +state_uncert_indicator=0; for file = 1:ifil loaded_file=load([DirectoryName '/' M_.fname var_type int2str(file)]); @@ -118,6 +119,13 @@ for file = 1:ifil end k = k(end)+(1:size(stock,2)); stock1(:,k) = stock; + elseif strcmp(var_type,'_state_uncert') + if file==1 %on first run, initialize variable for storing filter_step_ahead + stock1_state_uncert=NaN(n1,n2,size(stock,3),B); + end + state_uncert_indicator=1; + k = k(end)+(1:size(stock,4)); + stock1_state_uncert(:,:,:,k) = stock; else if file==1 %on first run, initialize variable for storing filter_step_ahead stock1 = zeros(n1,n2,B); @@ -138,8 +146,7 @@ if filter_step_ahead_indicator if options_.estimation.moments_posterior_density.indicator Density_filter_step_ahead = zeros(options_.estimation.moments_posterior_density.gridpoints,2,filter_steps,nvar,n2); end -end -if filter_covar_indicator +elseif filter_covar_indicator draw_dimension=4; oo_.FilterCovariance.Mean = squeeze(mean(stock1_filter_covar,draw_dimension)); oo_.FilterCovariance.Median = squeeze(median(stock1_filter_covar,draw_dimension)); @@ -161,6 +168,28 @@ if filter_covar_indicator oo_.FilterCovariance.HPDsup=squeeze(hpd_interval(:,:,:,2)); fprintf(['Estimation::mcmc: ' tit1 ', done!\n']); return +elseif state_uncert_indicator + draw_dimension=4; + oo_.Smoother.State_uncertainty.Mean = squeeze(mean(stock1_state_uncert,draw_dimension)); + oo_.Smoother.State_uncertainty.Median = squeeze(median(stock1_state_uncert,draw_dimension)); + oo_.Smoother.State_uncertainty.var = squeeze(var(stock1_state_uncert,0,draw_dimension)); + if size(stock1_state_uncert,draw_dimension)>2 + hpd_interval = quantile(stock1_state_uncert,[(1-options_.mh_conf_sig)/2 (1-options_.mh_conf_sig)/2+options_.mh_conf_sig],draw_dimension); + else + size_matrix=size(stock1_state_uncert); + hpd_interval=NaN([size_matrix(1:3),2]); + end + if size(stock1_state_uncert,draw_dimension)>9 + post_deciles =quantile(stock1_state_uncert,[0.1:0.1:0.9],draw_dimension); + else + size_matrix=size(stock1_state_uncert); + post_deciles=NaN([size_matrix(1:3),9]); + end + oo_.Smoother.State_uncertainty.post_deciles=post_deciles; + oo_.Smoother.State_uncertainty.HPDinf=squeeze(hpd_interval(:,:,:,1)); + oo_.Smoother.State_uncertainty.HPDsup=squeeze(hpd_interval(:,:,:,2)); + fprintf(['Estimation::mcmc: ' tit1 ', done!\n']); + return end if strcmp(var_type,'_trend_coeff') %two dimensional arrays diff --git a/matlab/prior_posterior_statistics.m b/matlab/prior_posterior_statistics.m index f4832ef47..5d29d01a9 100644 --- a/matlab/prior_posterior_statistics.m +++ b/matlab/prior_posterior_statistics.m @@ -117,6 +117,10 @@ if options_.filter_covariance MAX_filter_covariance = min(B,ceil(MaxNumberOfBytes/(endo_nbr^2*(gend+1))/8)); end +if options_.smoothed_state_uncertainty + MAX_n_smoothed_state_uncertainty = min(B,ceil(MaxNumberOfBytes/((endo_nbr*endo_nbr)*gend)/8)); +end + varlist = options_.varlist; if isempty(varlist) varlist = M_.endo_names(1:M_.orig_endo_nbr, :); @@ -129,13 +133,13 @@ for i=1:nvar end end -n_variables_to_fill=12; +n_variables_to_fill=13; irun = ones(n_variables_to_fill,1); ifil = zeros(n_variables_to_fill,1); run_smoother = 0; -if options_.smoother || options_.forecast || options_.filter_step_ahead +if options_.smoother || options_.forecast || options_.filter_step_ahead || options_.smoothed_state_uncertainty run_smoother = 1; end @@ -144,11 +148,17 @@ if options_.filter_covariance filter_covariance=1; end +smoothed_state_uncertainty=0; +if options_.smoothed_state_uncertainty + smoothed_state_uncertainty=1; +end + % Store the variable mandatory for local/remote parallel computing. localVars.type=type; localVars.run_smoother=run_smoother; localVars.filter_covariance=filter_covariance; +localVars.smoothed_state_uncertainty=smoothed_state_uncertainty; localVars.gend=gend; localVars.Y=Y; localVars.data_index=data_index; @@ -181,6 +191,9 @@ end if options_.filter_covariance localVars.MAX_filter_covariance = MAX_filter_covariance; end +if options_.smoothed_state_uncertainty + localVars.MAX_n_smoothed_state_uncertainty = MAX_n_smoothed_state_uncertainty ; +end localVars.MAX_n_smoothed_constant=MAX_n_smoothed_constant; localVars.MAX_n_smoothed_trend=MAX_n_smoothed_trend; localVars.MAX_n_trend_coeff=MAX_n_trend_coeff; @@ -246,6 +259,10 @@ else ifil(10,j+1) =ifil(10,j)+nfiles; nfiles = ceil(nBlockPerCPU(j)/MAX_n_smoothed_trend); ifil(11,j+1) =ifil(11,j)+nfiles; + if smoothed_state_uncertainty + nfiles = ceil(nBlockPerCPU(j)/MAX_n_smoothed_state_uncertainty); + ifil(13,j+1) =ifil(13,j)+nfiles; + end end end localVars.ifil = ifil; @@ -297,7 +314,12 @@ if options_.smoother pm3(endo_nbr,gend,ifil(11),B,'Smoothed trend',... '',M_.endo_names(1:M_.orig_endo_nbr, :),M_.endo_names_tex,M_.endo_names,... varlist,'Trend',DirectoryName,'_smoothed_trend'); - + if smoothed_state_uncertainty + pm3(endo_nbr,endo_nbr,ifil(13),B,'State Uncertainty',... + '',varlist,M_.endo_names_tex,M_.endo_names,... + varlist,'StateUncertainty',DirectoryName,'_state_uncert'); + end + if nvn for obs_iter=1:length(options_.varobs) meas_error_names{obs_iter,1}=['SE_EOBS_' M_.endo_names(strmatch(options_.varobs{obs_iter},M_.endo_names,'exact'),:)]; diff --git a/matlab/prior_posterior_statistics_core.m b/matlab/prior_posterior_statistics_core.m index 5b2a8e1bc..df250d79a 100644 --- a/matlab/prior_posterior_statistics_core.m +++ b/matlab/prior_posterior_statistics_core.m @@ -21,6 +21,7 @@ function myoutput=prior_posterior_statistics_core(myinputs,fpar,B,whoiam, ThisMa % _trend_coeff; % _smoothed_trend; % _smoothed_constant; +% _state_uncert; % % ALGORITHM % Portion of prior_posterior.m function. @@ -58,6 +59,7 @@ end type=myinputs.type; run_smoother=myinputs.run_smoother; filter_covariance=myinputs.filter_covariance; +smoothed_state_uncertainty=myinputs.smoothed_state_uncertainty; gend=myinputs.gend; Y=myinputs.Y; data_index=myinputs.data_index; @@ -83,8 +85,12 @@ if naK MAX_naK=myinputs.MAX_naK; end if filter_covariance - MAX_filter_covariance=myinputs.MAX_filter_covariance; + MAX_filter_covariance=myinputs.MAX_filter_covariance; end +if smoothed_state_uncertainty + MAX_n_smoothed_state_uncertainty=myinputs.MAX_n_smoothed_state_uncertainty; +end + exo_nbr=myinputs.exo_nbr; maxlag=myinputs.maxlag; MAX_nsmoo=myinputs.MAX_nsmoo; @@ -177,6 +183,9 @@ stock_ys = NaN(MAX_nruns,endo_nbr); if filter_covariance stock_filter_covariance = zeros(endo_nbr,endo_nbr,gend+1,MAX_filter_covariance); end +if smoothed_state_uncertainty + stock_smoothed_uncert = zeros(endo_nbr,endo_nbr,gend,MAX_n_smoothed_state_uncertainty); +end for b=fpar:B if strcmpi(type,'prior') @@ -195,7 +204,7 @@ for b=fpar:B if run_smoother [dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_); - [alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK,junk1,junk2,P,junk4,junk5,trend_addition] = ... + [alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK,junk1,junk2,P,junk4,junk5,trend_addition,state_uncertainty] = ... DsgeSmoother(deep,gend,Y,data_index,missing_value); stock_trend_coeff(options_.varobs_id,irun(9))=trend_coeff; @@ -322,6 +331,9 @@ for b=fpar:B if filter_covariance stock_filter_covariance(dr.order_var,dr.order_var,:,irun(8)) = P; end + if smoothed_state_uncertainty + stock_smoothed_uncert(dr.order_var,dr.order_var,:,irun(13)) = state_uncertainty; + end else [T,R,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_); end @@ -330,7 +342,7 @@ for b=fpar:B stock_ys(irun(5),:) = SteadyState'; - irun = irun + ones(12,1); + irun = irun + ones(13,1); if run_smoother && (irun(1) > MAX_nsmoo || b == B), @@ -463,6 +475,17 @@ for b=fpar:B irun(irun_index) = 1; end + irun_index=13; + if run_smoother && smoothed_state_uncertainty && (irun(irun_index) > MAX_n_smoothed_state_uncertainty || b == B) + stock = stock_smoothed_uncert(:,:,:,1:irun(irun_index)-1); + ifil(irun_index) = ifil(irun_index) + 1; + save([DirectoryName '/' M_.fname '_state_uncert' int2str(ifil(irun_index)) '.mat'],'stock'); + if RemoteFlag==1, + OutputFileName_state_uncert = [OutputFileName_state_uncert; {[DirectoryName filesep], [M_.fname '_state_uncert' int2str(ifil(irun_index)) '.mat']}]; + end + irun(irun_index) = 1; + end + dyn_waitbar((b-fpar+1)/(B-fpar+1),h); end @@ -480,7 +503,8 @@ if RemoteFlag==1, OutputFileName_filter_covar; OutputFileName_trend_coeff; OutputFileName_smoothed_trend; - OutputFileName_smoothed_constant]; + OutputFileName_smoothed_constant; + OutputFileName_state_uncert]; end dyn_waitbar_close(h); diff --git a/matlab/store_smoother_results.m b/matlab/store_smoother_results.m index 44782ed28..61af7b6a6 100644 --- a/matlab/store_smoother_results.m +++ b/matlab/store_smoother_results.m @@ -1,4 +1,4 @@ -function [oo_, yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend) +function [oo_, yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty) % oo_=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend) % Writes the smoother results into respective fields in oo_ % @@ -24,6 +24,8 @@ function [oo_, yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,da % decomp [4D array] (K*m*r*(T+K)) 4D array of shock decomposition of k-step ahead % filtered variables (decision-rule order) % Trend [double] [nvarobs*T] matrix of trends in observables +% state_uncertainty [double] (K,K,T) array, storing the uncertainty +% about the smoothed state (decision-rule order) % % Outputs: % oo_ [structure] storing the results: @@ -38,7 +40,8 @@ function [oo_, yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,da % oo_.UpdatedVariables: structure storing the updated variables % oo_.SmoothedShocks: structure storing the smoothed shocks % oo_.SmoothedMeasurementErrors: structure storing the smoothed measurement errors -% +% oo_.Smoother.State_uncertainty: smoothed state uncertainty (declaration order) + % yf [double] (nvarobs*T) matrix storing the smoothed observed variables (order of options_.varobs) % % Notes: @@ -117,6 +120,9 @@ if options_.filter_covariance oo_.Smoother.Variance = P; end +if options_.smoothed_state_uncertainty + oo_.Smoother.State_uncertainty=state_uncertainty; +end %get indices of smoothed variables i_endo_in_bayestopt_smoother_varlist = bayestopt_.smoother_saved_var_list; i_endo_in_dr_matrices=bayestopt_.smoother_var_list(i_endo_in_bayestopt_smoother_varlist); @@ -207,6 +213,9 @@ end if options_.filter_covariance oo_.Smoother.Variance(oo_.dr.order_var,oo_.dr.order_var,:)=oo_.Smoother.Variance; end +if options_.smoothed_state_uncertainty + oo_.Smoother.State_uncertainty(oo_.dr.order_var,oo_.dr.order_var,:)=state_uncertainty; +end %% get smoothed shocks for exo_iter=1:M_.exo_nbr diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy index a2e55f055..48aacafbc 100644 --- a/preprocessor/DynareBison.yy +++ b/preprocessor/DynareBison.yy @@ -153,7 +153,7 @@ class ParsingDriver; %token ALPHA BETA ABAND NINV CMS NCMS CNUM GAMMA INV_GAMMA INV_GAMMA1 INV_GAMMA2 NORMAL UNIFORM EPS PDF FIG DR NONE PRIOR PRIOR_VARIANCE HESSIAN IDENTITY_MATRIX DIRICHLET %token GSIG2_LMDM Q_DIAG FLAT_PRIOR NCSK NSTD WEIBULL WEIBULL_PDF %token INDXPARR INDXOVR INDXAP APBAND INDXIMF IMFBAND INDXFORE FOREBAND INDXGFOREHAT INDXGIMFHAT -%token INDXESTIMA INDXGDLS EQ_MS FILTER_COVARIANCE FILTER_DECOMPOSITION +%token INDXESTIMA INDXGDLS EQ_MS FILTER_COVARIANCE FILTER_DECOMPOSITION SMOOTHED_STATE_UNCERTAINTY %token EQ_CMS TLINDX TLNUMBER BANACT RESTRICTIONS POSTERIOR_SAMPLER_OPTIONS %token OUTPUT_FILE_TAG DRAWS_NBR_BURN_IN_1 DRAWS_NBR_BURN_IN_2 HORIZON %token SBVAR TREND_VAR DEFLATOR GROWTH_FACTOR MS_IRF MS_VARIANCE_DECOMPOSITION @@ -1777,6 +1777,7 @@ estimation_options : o_datafile | o_partial_information | o_filter_covariance | o_filter_decomposition + | o_smoothed_state_uncertainty | o_selected_variables_only | o_conditional_variance_decomposition | o_cova_compute @@ -2581,6 +2582,7 @@ calib_smoother_option : o_filtered_vars | o_filter_decomposition | o_diffuse_kalman_tol | o_diffuse_filter + | o_smoothed_state_uncertainty ; extended_path : EXTENDED_PATH ';' @@ -3129,6 +3131,9 @@ o_filter_covariance : FILTER_COVARIANCE o_filter_decomposition : FILTER_DECOMPOSITION { driver.option_num("filter_decomposition","1");} ; +o_smoothed_state_uncertainty : SMOOTHED_STATE_UNCERTAINTY + { driver.option_num("smoothed_state_uncertainty","1");} + ; o_selected_variables_only : SELECTED_VARIABLES_ONLY { driver.option_num("selected_variables_only","1");} ; diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll index e9213a45d..936bb7b96 100644 --- a/preprocessor/DynareFlex.ll +++ b/preprocessor/DynareFlex.ll @@ -573,6 +573,7 @@ DATE -?[0-9]+([YyAa]|[Mm]([1-9]|1[0-2])|[Qq][1-4]|[Ww]([1-9]{1}|[1-4][0-9]|5[0-2 k_order_solver {return token::K_ORDER_SOLVER; } filter_covariance {return token::FILTER_COVARIANCE; } filter_decomposition {return token::FILTER_DECOMPOSITION; } +smoothed_state_uncertainty {return token::SMOOTHED_STATE_UNCERTAINTY; } selected_variables_only {return token::SELECTED_VARIABLES_ONLY; } pruning {return token::PRUNING; } save_draws {return token::SAVE_DRAWS; } diff --git a/tests/TeX/fs2000_corr_ME.mod b/tests/TeX/fs2000_corr_ME.mod index 419851fb5..fe432c944 100644 --- a/tests/TeX/fs2000_corr_ME.mod +++ b/tests/TeX/fs2000_corr_ME.mod @@ -168,7 +168,8 @@ end; write_latex_prior_table; -estimation(mode_compute=8,order=1,datafile='../fs2000/fsdat_simul',mode_check,smoother,filter_decomposition,mh_replic=4000, mh_nblocks=1, mh_jscale=0.8,forecast = 8,bayesian_irf,filtered_vars,filter_step_ahead=[1,3],irf=20,moments_varendo,contemporaneous_correlation,conditional_variance_decomposition=[1 2 4]) m P c e W R k d y gy_obs; +estimation(mode_compute=8,order=1,datafile='../fs2000/fsdat_simul',mode_check,smoother,filter_decomposition,mh_replic=4000, mh_nblocks=1, mh_jscale=0.8,forecast = 8,bayesian_irf,filtered_vars,filter_step_ahead=[1,3],irf=20, + moments_varendo,contemporaneous_correlation,conditional_variance_decomposition=[1 2 4],smoothed_state_uncertainty) m P c e W R k d y gy_obs; trace_plot(options_,M_,estim_params_,'PosteriorDensity',1); trace_plot(options_,M_,estim_params_,'StructuralShock',1,'e_a') diff --git a/tests/kalman_filter_smoother/algo1.mod b/tests/kalman_filter_smoother/algo1.mod index 29e5b0203..95e34daf1 100644 --- a/tests/kalman_filter_smoother/algo1.mod +++ b/tests/kalman_filter_smoother/algo1.mod @@ -31,8 +31,8 @@ stderr e_z, uniform_pdf,,, 0.01, 0.1; end; varobs dw dx dy z; - -estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,filtered_vars); + +estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,filtered_vars,smoothed_state_uncertainty); //checking smoother consistency X = oo_.SmoothedVariables; diff --git a/tests/kalman_filter_smoother/algo2.mod b/tests/kalman_filter_smoother/algo2.mod index 0c54b4cee..b5dde589c 100644 --- a/tests/kalman_filter_smoother/algo2.mod +++ b/tests/kalman_filter_smoother/algo2.mod @@ -31,8 +31,8 @@ stderr e_z, uniform_pdf,,, 0.01, 0.1; end; varobs dw dx dy z; - -estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,mode_compute=0,mode_file=algo1_mode,kalman_algo=2,filtered_vars); + +estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,mode_compute=0,mode_file=algo1_mode,kalman_algo=2,filtered_vars,smoothed_state_uncertainty); //checking smoother consistency X = oo_.SmoothedVariables; diff --git a/tests/kalman_filter_smoother/algo3.mod b/tests/kalman_filter_smoother/algo3.mod index 22f1ed5b3..4a7b0c668 100644 --- a/tests/kalman_filter_smoother/algo3.mod +++ b/tests/kalman_filter_smoother/algo3.mod @@ -34,8 +34,8 @@ stderr e_z, inv_gamma_pdf,0.01, inf; end; varobs w x y; - -estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,diffuse_filter,filtered_vars); + +estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,diffuse_filter,filtered_vars,smoothed_state_uncertainty); //checking smoother consistency X = oo_.SmoothedVariables; diff --git a/tests/kalman_filter_smoother/algo4.mod b/tests/kalman_filter_smoother/algo4.mod index 7ad9e208e..e46b9b1b6 100644 --- a/tests/kalman_filter_smoother/algo4.mod +++ b/tests/kalman_filter_smoother/algo4.mod @@ -34,8 +34,8 @@ stderr e_z, inv_gamma_pdf,0.01, inf; end; varobs w x y; - -estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,mode_compute=0,mode_file=algo3_mode,diffuse_filter,kalman_algo=4,filtered_vars); + +estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,mode_compute=0,mode_file=algo3_mode,diffuse_filter,kalman_algo=4,filtered_vars,smoothed_state_uncertainty); //checking smoother consistency X = oo_.SmoothedVariables; diff --git a/tests/kalman_filter_smoother/algo4a.mod b/tests/kalman_filter_smoother/algo4a.mod index 9ca467501..37d743bbb 100644 --- a/tests/kalman_filter_smoother/algo4a.mod +++ b/tests/kalman_filter_smoother/algo4a.mod @@ -33,7 +33,7 @@ end; varobs dw dx y z; -estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,diffuse_filter); +estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,diffuse_filter,smoothed_state_uncertainty); //estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,mode_compute=0,mode_file=algo3_mode,diffuse_filter); //checking smoother consistency diff --git a/tests/kalman_filter_smoother/algo4b.mod b/tests/kalman_filter_smoother/algo4b.mod index 9ca467501..37d743bbb 100644 --- a/tests/kalman_filter_smoother/algo4b.mod +++ b/tests/kalman_filter_smoother/algo4b.mod @@ -33,7 +33,7 @@ end; varobs dw dx y z; -estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,diffuse_filter); +estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,diffuse_filter,smoothed_state_uncertainty); //estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,mode_compute=0,mode_file=algo3_mode,diffuse_filter); //checking smoother consistency diff --git a/tests/kalman_filter_smoother/algoH1.mod b/tests/kalman_filter_smoother/algoH1.mod index b09c3ef24..200dec37f 100644 --- a/tests/kalman_filter_smoother/algoH1.mod +++ b/tests/kalman_filter_smoother/algoH1.mod @@ -34,7 +34,7 @@ end; varobs dw dx dy z; -estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,filtered_vars); +estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,filtered_vars,smoothed_state_uncertainty); //checking smoother consistency X = oo_.SmoothedVariables; diff --git a/tests/kalman_filter_smoother/algoH2.mod b/tests/kalman_filter_smoother/algoH2.mod index 70a7b1de0..bd54deabe 100644 --- a/tests/kalman_filter_smoother/algoH2.mod +++ b/tests/kalman_filter_smoother/algoH2.mod @@ -35,7 +35,7 @@ end; varobs dw dx dy z; //estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,kalman_algo=2); -estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,mode_compute=0,mode_file=algoH1_mode,kalman_algo=2,filtered_vars); +estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,mode_compute=0,mode_file=algoH1_mode,kalman_algo=2,filtered_vars,smoothed_state_uncertainty); //checking smoother consistency X = oo_.SmoothedVariables; diff --git a/tests/kalman_filter_smoother/algoH3.mod b/tests/kalman_filter_smoother/algoH3.mod index d3f923d34..bc4d77436 100644 --- a/tests/kalman_filter_smoother/algoH3.mod +++ b/tests/kalman_filter_smoother/algoH3.mod @@ -37,7 +37,7 @@ end; varobs w x y; -estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,diffuse_filter); +estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,diffuse_filter,smoothed_state_uncertainty); stoch_simul(irf=0); diff --git a/tests/kalman_filter_smoother/fs2000_smoother_only.mod b/tests/kalman_filter_smoother/fs2000_smoother_only.mod index 696f68c05..b24c8b94b 100644 --- a/tests/kalman_filter_smoother/fs2000_smoother_only.mod +++ b/tests/kalman_filter_smoother/fs2000_smoother_only.mod @@ -100,11 +100,11 @@ check; varobs gp_obs gy_obs; -estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear, smoother) m P c e W R k d n l gy_obs gp_obs y dA; -estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear, smoother,kalman_algo=1) m P c e W R k d n l gy_obs gp_obs y dA; -estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear, smoother,kalman_algo=2) m P c e W R k d n l gy_obs gp_obs y dA; -estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear, smoother,kalman_algo=3) m P c e W R k d n l gy_obs gp_obs y dA; -estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear, smoother,kalman_algo=4) m P c e W R k d n l gy_obs gp_obs y dA; +estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear, smoother, smoothed_state_uncertainty) m P c e W R k d n l gy_obs gp_obs y dA; +estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear, smoother,kalman_algo=1, smoothed_state_uncertainty) m P c e W R k d n l gy_obs gp_obs y dA; +estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear, smoother,kalman_algo=2, smoothed_state_uncertainty) m P c e W R k d n l gy_obs gp_obs y dA; +estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear, smoother,kalman_algo=3, smoothed_state_uncertainty) m P c e W R k d n l gy_obs gp_obs y dA; +estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear, smoother,kalman_algo=4, smoothed_state_uncertainty) m P c e W R k d n l gy_obs gp_obs y dA; /* * The following lines were used to generate the data file. If you want to diff --git a/tests/kalman_filter_smoother/fs2000_smoother_only_ns.mod b/tests/kalman_filter_smoother/fs2000_smoother_only_ns.mod index 5688611a3..557486b79 100644 --- a/tests/kalman_filter_smoother/fs2000_smoother_only_ns.mod +++ b/tests/kalman_filter_smoother/fs2000_smoother_only_ns.mod @@ -107,9 +107,9 @@ P_obs (log(mst)-gam); Y_obs (gam); end; -estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear,diffuse_filter, smoother) m P c e W R k d n l gy_obs gp_obs y dA; -estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear,diffuse_filter, smoother,kalman_algo=3) m P c e W R k d n l gy_obs gp_obs y dA; -estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear,diffuse_filter, smoother,kalman_algo=4) m P c e W R k d n l gy_obs gp_obs y dA; +estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear,diffuse_filter, smoother,smoothed_state_uncertainty) m P c e W R k d n l gy_obs gp_obs y dA; +estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear,diffuse_filter, smoother,kalman_algo=3,smoothed_state_uncertainty) m P c e W R k d n l gy_obs gp_obs y dA; +estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear,diffuse_filter, smoother,kalman_algo=4,smoothed_state_uncertainty) m P c e W R k d n l gy_obs gp_obs y dA; /* * The following lines were used to generate the data file. If you want to