Commit 7c8f4867 authored by Johannes Pfeifer's avatar Johannes Pfeifer
Browse files

Allow setting initial state for Kalman filter

Supersedes Dynare/dynare!1522
parent 8f47b276
......@@ -8646,6 +8646,37 @@ the :comm:`bvar_forecast` command.
necessary to specify the logarithm of the controlled variables.
.. block:: filter_initial_state ;
|br| This block specifies the initial values of the endogenous states at the beginning
of the Kalman filter recursions. That is, if the Kalman filter recursion starts with
time t=1 being the first observation, this block provides the state estimate at time 0
given information at time 0, :math:`E_0(x_0)`. If nothing is specified, the initial condition is assumed to be at
the steady state (which is the unconditional mean for a stationary model).
This block is terminated by ``end;``.
Each line inside of the block should be of the form::
VARIABLE_NAME(INTEGER)=EXPRESSION;
``EXPRESSION`` is any valid expression returning a numerical value and can contain parameter values. This
allows specifying relationships that will be honored during estimation. ``INTEGER`` refers to the lag
with which a variable appears. By convention in Dynare, period 1 is the first period. Going backwards in time,
the first period before the start of the simulation is period 0, then period -1, and so on.
Note that the ``filter_initial_state`` block does not take non-state variables.
*Example*
::
filter_initial_state;
k(0)= ((1/bet-(1-del))/alp)^(1/(alp-1))*l_ss;
P(0)=2.5258;
m(0)= mst;
end;
.. command:: plot_conditional_forecast [VARIABLE_NAME...];
plot_conditional_forecast (periods = INTEGER) [VARIABLE_NAME...];
......
......@@ -97,6 +97,7 @@ end
%------------------------------------------------------------------------------
% 2. call model setup & reduction program
%------------------------------------------------------------------------------
%store old setting of restricted var_list
oldoo.restrict_var_list = oo_.dr.restrict_var_list;
oldoo.restrict_columns = oo_.dr.restrict_columns;
oo_.dr.restrict_var_list = bayestopt_.smoother_var_list;
......@@ -108,8 +109,6 @@ if info~=0
print_info(info,options_.noprint, options_);
return
end
oo_.dr.restrict_var_list = oldoo.restrict_var_list;
oo_.dr.restrict_columns = oldoo.restrict_columns;
%get location of observed variables and requested smoothed variables in
%decision rules
......@@ -227,7 +226,10 @@ ST = T;
R1 = R;
if kalman_algo == 1 || kalman_algo == 3
[alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty] = missing_DiffuseKalmanSmootherH1_Z(ST, ...
a_initial = zeros(np,1);
a_initial=set_Kalman_smoother_starting_values(a_initial,M_,oo_,options_);
a_initial=T*a_initial; %set state prediction for first Kalman step;
[alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty] = missing_DiffuseKalmanSmootherH1_Z(a_initial,ST, ...
Z,R1,Q,H,Pinf,Pstar, ...
data1,vobs,np,smpl,data_index, ...
options_.nk,kalman_tol,diffuse_kalman_tol,options_.filter_decomposition,options_.smoothed_state_uncertainty);
......@@ -271,7 +273,10 @@ if kalman_algo == 2 || kalman_algo == 4
end
end
[alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty] = missing_DiffuseKalmanSmootherH3_Z(ST, ...
a_initial = zeros(np,1);
a_initial=set_Kalman_smoother_starting_values(a_initial,M_,oo_,options_);
a_initial=ST*a_initial; %set state prediction for first Kalman step;
[alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty] = missing_DiffuseKalmanSmootherH3_Z(a_initial,ST, ...
Z,R1,Q,diag(H), ...
Pinf,Pstar,data1,vobs,np,smpl,data_index, ...
options_.nk,kalman_tol,diffuse_kalman_tol, ...
......@@ -301,3 +306,36 @@ if expanded_state_vector_for_univariate_filter && (kalman_algo == 2 || kalman_al
state_uncertainty = state_uncertainty(k,k,:);
end
end
%reset old setting of restricted var_list
oo_.dr.restrict_var_list = oldoo.restrict_var_list;
oo_.dr.restrict_columns = oldoo.restrict_columns;
function a=set_Kalman_smoother_starting_values(a,M_,oo_,options_)
% function a=set_Kalman_smoother_starting_values(a,M_,oo_,options_)
% Sets initial states guess for Kalman filter/smoother based on M_.filter_initial_state
%
% INPUTS
% o a [double] (p*1) vector of states
% o M_ [structure] decribing the model
% o oo_ [structure] storing the results
% o options_ [structure] describing the options
%
% OUTPUTS
% o a [double] (p*1) vector of set initial states
if isfield(M_,'filter_initial_state') && ~isempty(M_.filter_initial_state)
state_indices=oo_.dr.order_var(oo_.dr.restrict_columns);
for ii=1:size(state_indices,1)
if ~isempty(M_.filter_initial_state{state_indices(ii),1})
if options_.loglinear && ~options_.logged_steady_state
a(oo_.dr.restrict_columns(ii)) = log(eval(M_.filter_initial_state{state_indices(ii),2})) - log(oo_.dr.ys(state_indices(ii)));
elseif ~options_.loglinear && ~options_.logged_steady_state
a(oo_.dr.restrict_columns(ii)) = eval(M_.filter_initial_state{state_indices(ii),2}) - oo_.dr.ys(state_indices(ii));
else
error(['The steady state is logged. This should not happen. Please contact the developers'])
end
end
end
end
......@@ -287,6 +287,8 @@ switch DynareOptions.lik_init
Pstar=lyapunov_solver(T,R,Q,DynareOptions);
Pinf = [];
a = zeros(mm,1);
a=set_Kalman_starting_values(a,Model,DynareResults,DynareOptions,BayesInfo);
a_0_given_tm1=T*a; %set state prediction for first Kalman step;
Zflag = 0;
case 2% Initialization with large numbers on the diagonal of the covariance matrix if the states (for non stationary models).
if kalman_algo ~= 2
......@@ -296,6 +298,8 @@ switch DynareOptions.lik_init
Pstar = DynareOptions.Harvey_scale_factor*eye(mm);
Pinf = [];
a = zeros(mm,1);
a = set_Kalman_starting_values(a,Model,DynareResults,DynareOptions,BayesInfo);
a_0_given_tm1 = T*a; %set state prediction for first Kalman step;
Zflag = 0;
case 3% Diffuse Kalman filter (Durbin and Koopman)
% Use standard kalman filter except if the univariate filter is explicitely choosen.
......@@ -314,16 +318,19 @@ switch DynareOptions.lik_init
% Run diffuse kalman filter on first periods.
if (kalman_algo==3)
% Multivariate Diffuse Kalman Filter
a = zeros(mm,1);
a = set_Kalman_starting_values(a,Model,DynareResults,DynareOptions,BayesInfo);
a_0_given_tm1 = T*a; %set state prediction for first Kalman step;
Pstar0 = Pstar; % store Pstar
if no_missing_data_flag
[dLIK,dlik,a,Pstar] = kalman_filter_d(Y, 1, size(Y,2), ...
zeros(mm,1), Pinf, Pstar, ...
[dLIK,dlik,a_0_given_tm1,Pstar] = kalman_filter_d(Y, 1, size(Y,2), ...
a_0_given_tm1, Pinf, Pstar, ...
kalman_tol, diffuse_kalman_tol, riccati_tol, DynareOptions.presample, ...
T,R,Q,H,Z,mm,pp,rr);
else
[dLIK,dlik,a,Pstar] = missing_observations_kalman_filter_d(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations, ...
[dLIK,dlik,a_0_given_tm1,Pstar] = missing_observations_kalman_filter_d(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations, ...
Y, 1, size(Y,2), ...
zeros(mm,1), Pinf, Pstar, ...
a_0_given_tm1, Pinf, Pstar, ...
kalman_tol, diffuse_kalman_tol, riccati_tol, DynareOptions.presample, ...
T,R,Q,H,Z,mm,pp,rr);
end
......@@ -367,11 +374,14 @@ switch DynareOptions.lik_init
end
end
[dLIK,dlik,a,Pstar] = univariate_kalman_filter_d(DatasetInfo.missing.aindex,...
a = zeros(mmm,1);
a = set_Kalman_starting_values(a,Model,DynareResults,DynareOptions,BayesInfo);
a_0_given_tm1 = T*a;
[dLIK,dlik,a_0_given_tm1,Pstar] = univariate_kalman_filter_d(DatasetInfo.missing.aindex,...
DatasetInfo.missing.number_of_observations,...
DatasetInfo.missing.no_more_missing_observations, ...
Y, 1, size(Y,2), ...
zeros(mmm,1), Pinf, Pstar, ...
a_0_given_tm1, Pinf, Pstar, ...
kalman_tol, diffuse_kalman_tol, riccati_tol, DynareOptions.presample, ...
T,R,Q,H1,Z,mmm,pp,rr);
diffuse_periods = size(dlik,1);
......@@ -402,6 +412,8 @@ switch DynareOptions.lik_init
end
Pinf = [];
a = zeros(mm,1);
a = set_Kalman_starting_values(a,Model,DynareResults,DynareOptions,BayesInfo);
a_0_given_tm1 = T*a;
Zflag = 0;
case 5 % Old diffuse Kalman filter only for the non stationary variables
[eigenvect, eigenv] = eig(T);
......@@ -422,6 +434,8 @@ switch DynareOptions.lik_init
Pstar(stable, stable) = Pstar_tmp;
Pinf = [];
a = zeros(mm,1);
a = set_Kalman_starting_values(a,Model,DynareResults,DynareOptions,BayesInfo);
a_0_given_tm1 = T*a;
Zflag = 0;
otherwise
error('dsge_likelihood:: Unknown initialization approach for the Kalman filter!')
......@@ -595,14 +609,14 @@ if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
return
end
[LIK,lik] = kalman_filter_fast(Y,diffuse_periods+1,size(Y,2), ...
a,Pstar, ...
a_0_given_tm1,Pstar, ...
kalman_tol, riccati_tol, ...
DynareOptions.presample, ...
T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods, ...
analytic_deriv_info{:});
else
[LIK,lik] = kalman_filter(Y,diffuse_periods+1,size(Y,2), ...
a,Pstar, ...
a_0_given_tm1,Pstar, ...
kalman_tol, riccati_tol, ...
DynareOptions.rescale_prediction_error_covariance, ...
DynareOptions.presample, ...
......@@ -615,7 +629,7 @@ if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
T,R,Q,H,Pstar,Y,start,Z,kalman_tol,riccati_tol, Model.nz_state_var, Model.n_diag, Model.nobs_non_statevar);
else
[LIK,lik] = missing_observations_kalman_filter(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ...
a, Pstar, ...
a_0_given_tm1, Pstar, ...
kalman_tol, DynareOptions.riccati_tol, ...
DynareOptions.rescale_prediction_error_covariance, ...
DynareOptions.presample, ...
......@@ -695,9 +709,11 @@ if (kalman_algo==2) || (kalman_algo==4)
end
mmm = mm+pp;
if singularity_has_been_detected
a = zeros(mmm,1);
a_tmp = zeros(mmm,1);
a_tmp(1:length(a_0_given_tm1)) = a_0_given_tm1;
a_0_given_tm1 = a_tmp;
elseif ~expanded_state_vector_for_univariate_filter
a = [a; zeros(pp,1)];
a_0_given_tm1 = [a_0_given_tm1; zeros(pp,1)];
end
end
end
......@@ -705,7 +721,7 @@ if (kalman_algo==2) || (kalman_algo==4)
analytic_deriv_info{5}=DH;
end
[LIK, lik] = univariate_kalman_filter(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ...
a,Pstar, ...
a_0_given_tm1,Pstar, ...
DynareOptions.kalman_tol, ...
DynareOptions.riccati_tol, ...
DynareOptions.presample, ...
......@@ -832,3 +848,32 @@ if analytic_derivation==0 && nargout>3
lik=lik(start:end,:);
DLIK=[-lnprior; lik(:)];
end
function a=set_Kalman_starting_values(a,M_,oo_,options_,bayestopt_)
% function a=set_Kalman_starting_values(a,M_,oo_,options_,bayestopt_)
% Sets initial states guess for Kalman filter/smoother based on M_.filter_initial_state
%
% INPUTS
% o a [double] (p*1) vector of states
% o M_ [structure] decribing the model
% o oo_ [structure] storing the results
% o options_ [structure] describing the options
% o bayestopt_ [structure] describing the priors
%
% OUTPUTS
% o a [double] (p*1) vector of set initial states
if isfield(M_,'filter_initial_state') && ~isempty(M_.filter_initial_state)
state_indices=oo_.dr.order_var(oo_.dr.restrict_var_list(bayestopt_.mf0));
for ii=1:size(state_indices,1)
if ~isempty(M_.filter_initial_state{state_indices(ii),1})
if options_.loglinear && ~options_.logged_steady_state
a(bayestopt_.mf0(ii)) = log(eval(M_.filter_initial_state{state_indices(ii),2})) - log(oo_.dr.ys(state_indices(ii)));
elseif ~options_.loglinear && ~options_.logged_steady_state
a(bayestopt_.mf0(ii)) = eval(M_.filter_initial_state{state_indices(ii),2}) - oo_.dr.ys(state_indices(ii));
else
error(['The steady state is logged. This should not happen. Please contact the developers'])
end
end
end
end
\ No newline at end of file
......@@ -148,6 +148,19 @@ if any(info)
error('The initial value of the prior is -Inf')
end
if isfield(Model,'filter_initial_state') && ~isempty(Model.filter_initial_state)
state_indices=DynareResults.dr.order_var(DynareResults.dr.restrict_var_list(BayesInfo.mf0));
for ii=1:size(state_indices,1)
if ~isempty(Model.filter_initial_state{state_indices(ii),1})
try
evaluate_expression(Model.filter_initial_state{state_indices(ii),2},Model,DynareResults)
catch
fprintf('Unable to evaluate the expression\n %s \nfor the filter_initial_state of variable %s\n',Model.filter_initial_state{state_indices(ii),2},Model.endo_names(state_indices(ii),:))
end
end
end
end
if DynareOptions.ramsey_policy
%test whether specification matches
inst_nbr = size(DynareOptions.instruments,1);
......@@ -233,3 +246,8 @@ end
if ~isequal(DynareOptions.mode_compute,11)
disp(['Initial value of the log posterior (or likelihood): ' num2str(-fval)]);
end
function evaluate_expression(expression,M_,oo_)
% function evaluate_expression(expression,M_,oo_)
%evaluates expressions relying on M_ and oo_ having their original names
eval(expression);
\ No newline at end of file
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,atilde,P,aK,PK,decomp,V] = missing_DiffuseKalmanSmootherH1_Z(a_initial,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,state_uncertainty_flag)
% Computes the diffuse kalman smoother without measurement error, in the case of a non-singular var-cov matrix.
%
% INPUTS
% a_initial:mm*1 vector of initial (predicted) states
% T: mm*mm matrix
% Z: pp*mm matrix
% R: mm*rr matrix
......@@ -78,6 +79,7 @@ spinf = size(Pinf1);
spstar = size(Pstar1);
v = zeros(pp,smpl);
a = zeros(mm,smpl+1);
a(:,1) = a_initial;
atilde = zeros(mm,smpl);
aK = zeros(nk,mm,smpl+nk);
PK = zeros(nk,mm,mm,smpl+nk);
......
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,a,P,aK,PK,decomp,V] = missing_DiffuseKalmanSmootherH3_Z(a_initial,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.
%
% INPUTS
% a_initial:mm*1 vector of initial states
% T: mm*mm matrix state transition matrix
% Z: pp*mm matrix selector matrix for observables in augmented state vector
% R: mm*rr matrix second matrix of the state equation relating the structural innovations to the state variables
......@@ -89,6 +90,8 @@ spstar = size(Pstar1);
v = zeros(pp,smpl);
a = zeros(mm,smpl);
a1 = zeros(mm,smpl+1);
a(:,1) = a_initial;
a1(:,1) = a_initial;
aK = zeros(nk,mm,smpl+nk);
Fstar = zeros(pp,smpl);
......
Subproject commit 7468b4d8d7457cb11beae557d2fd8a9f2ce3e246
Subproject commit ad22e628749560d9c59b5f7002700a505102a234
......@@ -265,6 +265,10 @@ MODFILES = \
kalman/lik_init/fs2000_lik_init_3.mod \
kalman/lik_init/fs2000_lik_init_4.mod \
kalman/lik_init/fs2000_lik_init_5.mod \
kalman_initial_state/fs2000_smoother_only_initial_state.mod \
kalman_initial_state/fs2000_ns_smoother_only_initial_state.mod \
kalman_initial_state/fs2000_kalman_initial.mod \
kalman_initial_state/fs2000_kalman_initial_2_lag.mod \
kalman_filter_smoother/gen_data.mod \
kalman_filter_smoother/algo1.mod \
kalman_filter_smoother/algo2.mod \
......@@ -561,6 +565,7 @@ XFAIL_MODFILES = ramst_xfail.mod \
optimal_policy/Ramsey/ramsey_histval_xfail.mod \
particle/first_spec_xfail_0.mod \
particle/first_spec_xfail_1.mod \
kalman_initial_state/fs2000_kalman_initial_xfail.mod \
example1_extra_exo_xfail.mod \
estimation/tune_mh_jscale/fs2000_1_xfail.mod \
estimation/tune_mh_jscale/fs2000_2_xfail.mod
......@@ -1139,6 +1144,7 @@ EXTRA_DIST = \
kalman_filter_smoother/fsdat_simul.m \
kalman/lik_init/fs2000_common.inc \
kalman/lik_init/fs2000_ns_common.inc \
kalman_initial_state/fs2000_common.inc \
kalman_filter_smoother/compare_results_simulation/fsdat_simul_logged.m \
kalman_filter_smoother/SOE_data_file.m \
kalman/likelihood_from_dynare/fs2000_model.inc \
......
// See fs2000.mod in the examples/ directory for details on the model
var m P c e W R k d n l gy_obs gp_obs y dA;
varexo e_a e_m;
parameters alp bet gam mst rho psi del;
alp = 0.33;
bet = 0.99;
gam = 0.003;
mst = 1.011;
rho = 0.7;
psi = 0.787;
del = 0.02;
model;
dA = exp(gam+e_a);
log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
W = l/n;
-(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
P*c = m;
m-1+d = l;
e = exp(e_a);
y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
gy_obs = dA*y/y(-1);
gp_obs = (P/P(-1))*m(-1)/dA;
end;
steady_state_model;
dA = exp(gam);
gst = 1/dA;
m = mst;
khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1));
xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1);
nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp );
n = xist/(nust+xist);
P = xist + nust;
k = khst*n;
l = psi*mst*n/( (1-psi)*(1-n) );
c = mst/P;
d = l - mst + 1;
y = k^alp*n^(1-alp)*gst^alp;
R = mst/bet;
W = l/n;
ist = y-c;
q = 1 - d;
e = 1;
gp_obs = m/dA;
gy_obs = dA;
end;
shocks;
var e_a; stderr 0.014;
var e_m; stderr 0.005;
end;
@#include "fs2000_common.inc"
estimated_params;
alp, beta_pdf, 0.356, 0.02;
bet, beta_pdf, 0.993, 0.002;
gam, normal_pdf, 0.0085, 0.003;
mst, normal_pdf, 1.0002, 0.007;
rho, beta_pdf, 0.129, 0.223;
psi, beta_pdf, 0.65, 0.05;
del, beta_pdf, 0.01, 0.005;
stderr e_a, inv_gamma_pdf, 0.035449, inf;
stderr e_m, inv_gamma_pdf, 0.008862, inf;
end;
varobs gp_obs gy_obs;
options_.solve_tolf = 1e-12;
filter_initial_state;
k(0)= (((1-1/exp(gam)*bet*(1-del)) / (alp*1/exp(gam)^alp*bet) )^(1/(alp-1)))*(( (((((1-1/exp(gam)*bet*(1-del)) / (alp*1/exp(gam)^alp*bet) )^(1/(alp-1)))*1/exp(gam))^alp - (1-1/exp(gam)*(1-del))*(((1-1/exp(gam)*bet*(1-del)) / (alp*1/exp(gam)^alp*bet) )^(1/(alp-1))))/mst )^(-1))/(psi*mst^2/( (1-alp)*(1-psi)*bet*1/exp(gam)^alp*(((1-1/exp(gam)*bet*(1-del)) / (alp*1/exp(gam)^alp*bet) )^(1/(alp-1)))^alp )+(( (((((1-1/exp(gam)*bet*(1-del)) / (alp*1/exp(gam)^alp*bet) )^(1/(alp-1)))*1/exp(gam))^alp - (1-1/exp(gam)*(1-del))*(((1-1/exp(gam)*bet*(1-del)) / (alp*1/exp(gam)^alp*bet) )^(1/(alp-1))))/mst)^(-1)));
P(0)=2.5258;
m(0) = mst;
end;
estimation(order=1,datafile='../fs2000/fsdat_simul',nobs=192,loglinear,mh_replic=2001,mh_nblocks=1,mh_jscale=0.8,moments_varendo,consider_only_observed,smoother);
// See fs2000.mod in the examples/ directory for details on the model
var m P c e W R k d n l gy_obs gp_obs y dA;
varexo e_a e_m;
parameters alp bet gam mst rho psi del;
alp = 0.33;
bet = 0.99;
gam = 0.003;
mst = 1.011;
rho = 0.7;
psi = 0.787;
del = 0.02;
model;
dA = exp(gam+e_a);
log(m) = (1-rho)*log(mst) + rho*log(m(-2))+e_m;
-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
W = l/n;
-(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
P*c = m;
m-1+d = l;
e = exp(e_a);
y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
gy_obs = dA*y/y(-1);
gp_obs = (P/P(-1))*m(-1)/dA;
end;
steady_state_model;
dA = exp(gam);
gst = 1/dA;
m = mst;
khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1));
xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1);
nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp );
n = xist/(nust+xist);
P = xist + nust;
k = khst*n;
l = psi*mst*n/( (1-psi)*(1-n) );
c = mst/P;
d = l - mst + 1;
y = k^alp*n^(1-alp)*gst^alp;
R = mst/bet;
W = l/n;
ist = y-c;
q = 1 - d;
e = 1;
gp_obs = m/dA;
gy_obs = dA;
end;
shocks;
var e_a; stderr 0.014;
var e_m; stderr 0.005;
end;
estimated_params;
alp, beta_pdf, 0.356, 0.02;
bet, beta_pdf, 0.993, 0.002;
gam, normal_pdf, 0.0085, 0.003;
mst, normal_pdf, 1.0002, 0.007;
rho, beta_pdf, 0.129, 0.223;
psi, beta_pdf, 0.65, 0.05;
del, beta_pdf, 0.01, 0.005;
stderr e_a, inv_gamma_pdf, 0.035449, inf;
stderr e_m, inv_gamma_pdf, 0.008862, inf;
end;
varobs gp_obs gy_obs;
options_.solve_tolf = 1e-12;
filter_initial_state;
k(0)= (((1-1/exp(gam)*bet*(1-del)) / (alp*1/exp(gam)^alp*bet) )^(1/(alp-1)))*(( (((((1-1/exp(gam)*bet*(1-del)) / (alp*1/exp(gam)^alp*bet) )^(1/(alp-1)))*1/exp(gam))^alp - (1-1/exp(gam)*(1-del))*(((1-1/exp(gam)*bet*(1-del)) / (alp*1/exp(gam)^alp*bet) )^(1/(alp-1))))/mst )^(-1))/(psi*mst^2/( (1-alp)*(1-psi)*bet*1/exp(gam)^alp*(((1-1/exp(gam)*bet*(1-del)) / (alp*1/exp(gam)^alp*bet) )^(1/(alp-1)))^alp )+(( (((((1-1/exp(gam)*bet*(1-del)) / (alp*1/exp(gam)^alp*bet) )^(1/(alp-1)))*1/exp(gam))^alp - (1-1/exp(gam)*(1-del))*(((1-1/exp(gam)*bet*(1-del)) / (alp*1/exp(gam)^alp*bet) )^(1/(alp-1))))/mst)^(-1)));
P(0)=2.5258;
m(-1) = mst;
end;
estimation(order=1,datafile='../fs2000/fsdat_simul',nobs=192,loglinear,mh_replic=2001,mh_nblocks=1,mh_jscale=0.8,moments_varendo,consider_only_observed,smoother);
@#include "fs2000_common.inc"
estimated_params;
alp, beta_pdf, 0.356, 0.02;
bet, beta_pdf, 0.993, 0.002;
gam, normal_pdf, 0.0085, 0.003;
mst, normal_pdf, 1.0002, 0.007;
rho, beta_pdf, 0.129, 0.223;
psi, beta_pdf, 0.65, 0.05;
del, beta_pdf, 0.01, 0.005;
stderr e_a, inv_gamma_pdf, 0.035449, inf;
stderr e_m, inv_gamma_pdf, 0.008862, inf;
end;
varobs gp_obs gy_obs;
options_.solve_tolf = 1e-12;
filter_initial_state;
k(0)= (((1-1/exp(gam)*bet*(1-del)) / (alp*1/exp(gam)^alp*bet) )^(1/(alp-1)))*(( (((((1-1/exp(gam)*bet*(1-del)) / (alp*1/exp(gam)^alp*bet) )^(1/(alp-1)))*1/exp(gam))^alp - (1-1/exp(gam)*(1-del))*(((1-1/exp(gam)*bet*(1-del)) / (alp*1/exp(gam)^alp*bet) )^(1/(alp-1))))/mst )^(-1))/(psi*mst^2/( (1-alp)*(1-psi)*bet*1/exp(gam)^alp*(((1-1/exp(gam)*bet*(1-del)) / (alp*1/exp(gam)^alp*bet) )^(1/(alp-1)))^alp )+(( (((((1-1/exp(gam)*bet*(1-del)) / (alp*1/exp(gam)^alp*bet) )^(1/(alp-1)))*1/exp(gam))^alp - (1-1/exp(gam)*(1-del))*(((1-1/exp(gam)*bet*(1-del)) / (alp*1/exp(gam)^alp*bet) )^(1/(alp-1))))/mst)^(-1)));
c(0)=1;
end;
estimation(order=1,datafile='../fs2000/fsdat_simul',nobs=192,loglinear,mh_replic=2000,mh_nblocks=1,mh_jscale=0.8,moments_varendo,consider_only_observed);
/*
* This file replicates the estimation of the cash in advance model described
* Frank Schorfheide (2000): "Loss function-based evaluation of DSGE models",
* Journal of Applied Econometrics, 15(6), 645-670.
*
* The data are in file "fsdat_simul.m", and have been artificially generated.
* They are therefore different from the original dataset used by Schorfheide.
*
* The equations are taken from J. Nason and T. Cogley (1994): "Testing the
* implications of long-run neutrality for monetary business cycle models",
* Journal of Applied Econometrics, 9, S37-S70.
* Note that there is an initial minus sign missing in equation (A1), p. S63.
*
* This implementation was written by Michel Juillard. Please note that the