diff --git a/doc/dynare.texi b/doc/dynare.texi
index cf674053e3d17d89b3e6b71f16f3633067d69271..6cdd6852b02c1ecd27d82dc11563b7cd47ee6c27 100644
--- a/doc/dynare.texi
+++ b/doc/dynare.texi
@@ -4969,6 +4969,39 @@ end;
@end deffn
+@anchor{filter_initial_state}
+@deffn Block filter_initial_state ;
+
+@descriptionhead
+
+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 @code{end;}. Each line inside of the block should be of the form:
+
+@example
+@var{VARIABLE_NAME(INTEGER)}=(@var{EXPRESSION});
+@end example
+@code{EXPRESSION} is any valid expression returning a numerical value and can contain parameter values. This
+allows specifying relationships that will be honored during estimation. @code{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 @code{filter_initial_state} block does not take non-state variables.
+
+@examplehead
+
+@example
+filter_initial_state;
+k(0)= ((1/bet-(1-del))/alp)^(1/(alp-1))*l_ss;
+P(0)=2.5258;
+m(0)= mst;
+end;
+@end example
+
+@end deffn
@anchor{estimated_params}
@deffn Block estimated_params ;
diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m
index 6466ad2f1ef6fea5023645365987432abd775b6b..ef9c03c31f19f781e206ac79122e71689ca689ae 100644
--- a/matlab/DsgeSmoother.m
+++ b/matlab/DsgeSmoother.m
@@ -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
@@ -183,7 +182,7 @@ elseif options_.lik_init == 3 % Diffuse Kalman filter
Z = [Z, eye(vobs)];
end
end
- [Pstar,Pinf] = compute_Pinf_Pstar(mf,T,R,Q,options_.qz_criterium,oo_.dr.restrict_var_list);
+ [Pstar,Pinf] = compute_Pinf_Pstar(mf,T,R,Q,options_.qz_criterium);
elseif options_.lik_init == 4 % Start from the solution of the Riccati equation.
[err, Pstar] = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(mf,np,vobs)),H);
mexErrCheck('kalman_steady_state',err);
@@ -228,7 +227,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);
@@ -270,7 +272,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, ...
@@ -300,3 +305,35 @@ 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 ~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
\ No newline at end of file
diff --git a/matlab/dsge_likelihood.m b/matlab/dsge_likelihood.m
index 0b62e516753e59b6fd4e6c8332f1d08d1e7416ee..7ab6fd9571b25d4b33965cf70846727007d2ad89 100644
--- a/matlab/dsge_likelihood.m
+++ b/matlab/dsge_likelihood.m
@@ -354,6 +354,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
@@ -361,8 +363,10 @@ switch DynareOptions.lik_init
kalman_algo = 1;
end
Pstar = DynareOptions.Harvey_scale_factor*eye(mm);
- Pinf = [];
- a = zeros(mm,1);
+ 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.
@@ -381,16 +385,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
@@ -430,15 +437,17 @@ switch DynareOptions.lik_init
Pinf = blkdiag(Pinf,zeros(pp));
H1 = zeros(pp,1);
mmm = mm+pp;
-
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);
@@ -466,7 +475,9 @@ switch DynareOptions.lik_init
Pstar=lyapunov_solver(T,R,Q,DynareOptions);
end
Pinf = [];
- a = zeros(mm,1);
+ 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);
@@ -487,6 +498,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!')
@@ -644,14 +657,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, ...
@@ -664,7 +677,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, ...
@@ -744,9 +757,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
@@ -754,7 +769,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, ...
@@ -881,3 +896,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 ~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
diff --git a/matlab/initial_estimation_checks.m b/matlab/initial_estimation_checks.m
index 7dfc4d9a2a74e2bd2bdd7800b4a887c16234b944..3396a3c5afefcbbb5126a828d5f7f21dcf693e0a 100644
--- a/matlab/initial_estimation_checks.m
+++ b/matlab/initial_estimation_checks.m
@@ -111,6 +111,19 @@ if info
error('The initial value of the prior is -Inf')
end
+if ~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);
@@ -187,3 +200,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
diff --git a/matlab/kalman/likelihood/kalman_filter_d.m b/matlab/kalman/likelihood/kalman_filter_d.m
index 79b3607fb21de53cfc895071501657715dc9c332..bbdbfef31869acecc17d348e02d44d9c5bc9514a 100644
--- a/matlab/kalman/likelihood/kalman_filter_d.m
+++ b/matlab/kalman/likelihood/kalman_filter_d.m
@@ -5,7 +5,7 @@ function [dLIK,dlik,a,Pstar] = kalman_filter_d(Y, start, last, a, Pinf, Pstar, k
% Y [double] pp*smpl matrix of (detrended) data, where pp is the number of observed variables.
% start [integer] scalar, first observation.
% last [integer] scalar, last observation.
-% a [double] mm*1 vector, levels of the state variables.
+% a [double] mm*1 vector, levels of the predicted initial state variables (E_{0}(alpha_1)).
% Pinf [double] mm*mm matrix used to initialize the covariance matrix of the state vector.
% Pstar [double] mm*mm matrix used to initialize the covariance matrix of the state vector.
% kalman_tol [double] scalar, tolerance parameter (rcond) of F_star.
@@ -25,7 +25,7 @@ function [dLIK,dlik,a,Pstar] = kalman_filter_d(Y, start, last, a, Pinf, Pstar, k
% OUTPUTS
% LIK [double] scalar, minus loglikelihood
% lik [double] smpl*1 vector, log density of each vector of observations.
-% a [double] mm*1 vector, current estimate of the state vector.
+% a [double] mm*1 vector, current estimate of the state vector tomorrow (E_{T}(alpha_{T+1})).
% Pstar [double] mm*mm matrix, covariance matrix of the state vector.
%
% REFERENCES
diff --git a/matlab/kalman/likelihood/kalman_filter_fast.m b/matlab/kalman/likelihood/kalman_filter_fast.m
index db3dec3b755d406d55bf89005c80b6946b3a49e6..331019d95cd725a15435b99dec5d7cf4966f3441 100644
--- a/matlab/kalman/likelihood/kalman_filter_fast.m
+++ b/matlab/kalman/likelihood/kalman_filter_fast.m
@@ -20,7 +20,7 @@ function [LIK, LIKK, a, P] = kalman_filter_fast(Y,start,last,a,P,kalman_tol,ricc
%! @item last
%! Integer scalar, last period (@var{last}-@var{first} has to be inferior to T).
%! @item a
-%! Vector (@var{mm}*1) of doubles, initial mean of the state vector.
+%! Vector (@var{mm}*1) of doubles, levels of the predicted initial state variables (E_{0}(alpha_1)).
%! @item P
%! Matrix (@var{mm}*@var{mm}) of doubles, initial covariance matrix of the state vector.
%! @item kalman_tol
@@ -59,7 +59,7 @@ function [LIK, LIKK, a, P] = kalman_filter_fast(Y,start,last,a,P,kalman_tol,ricc
%! @item likk
%! Column vector of doubles, values of the density of each observation.
%! @item a
-%! Vector (@var{mm}*1) of doubles, mean of the state vector at the end of the (sub)sample.
+%! Vector (@var{mm}*1) of doubles, mean of the state vector at the end of the (sub)sample (E_{T}(alpha_{T+1})).
%! @item P
%! Matrix (@var{mm}*@var{mm}) of doubles, covariance of the state vector at the end of the (sub)sample.
%! @end table
diff --git a/matlab/kalman/likelihood/kalman_filter_ss.m b/matlab/kalman/likelihood/kalman_filter_ss.m
index 5510de9f36b9d262f470551878c6742b9e1e71fe..e94f7c3194d45bd932767172f66d099dc2761927 100644
--- a/matlab/kalman/likelihood/kalman_filter_ss.m
+++ b/matlab/kalman/likelihood/kalman_filter_ss.m
@@ -17,7 +17,7 @@ function [LIK, likk, a] = kalman_filter_ss(Y,start,last,a,T,K,iF,log_dF,Z,pp,Zfl
%! @item last
%! Integer scalar, last period (@var{last}-@var{first} has to be inferior to T).
%! @item a
-%! Vector (mm*1) of doubles, initial mean of the state vector.
+%! Vector (mm*1) of doubles, levels of the predicted initial state variables (E_{0}(alpha_1)).
%! @item T
%! Matrix (mm*mm) of doubles, transition matrix of the state equation.
%! @item K
@@ -42,7 +42,8 @@ function [LIK, likk, a] = kalman_filter_ss(Y,start,last,a,T,K,iF,log_dF,Z,pp,Zfl
%! @item likk
%! Column vector of doubles, values of the density of each observation.
%! @item a
-%! Vector (mm*1) of doubles, mean of the state vector at the end of the (sub)sample.
+%! Vector (mm*1) of doubles, current estimate of the state vector tomorrow
+%(E_{T}(alpha_{T+1})).
%! @end table
%! @sp 2
%! @strong{This function is called by:}
diff --git a/matlab/kalman/likelihood/missing_observations_kalman_filter.m b/matlab/kalman/likelihood/missing_observations_kalman_filter.m
index f1bd40cb59e5e0ce3198019d3e83c599228635e1..a2b33df30048bec5676b93996945b0434f16aed9 100644
--- a/matlab/kalman/likelihood/missing_observations_kalman_filter.m
+++ b/matlab/kalman/likelihood/missing_observations_kalman_filter.m
@@ -8,7 +8,7 @@ function [LIK, lik, a, P] = missing_observations_kalman_filter(data_index,numbe
% Y [double] pp*smpl matrix of data.
% start [integer] scalar, index of the first observation.
% last [integer] scalar, index of the last observation.
-% a [double] pp*1 vector, initial level of the state vector.
+% a [double] pp*1 vector, levels of the predicted initial state variables (E_{0}(alpha_1)).
% P [double] pp*pp matrix, covariance matrix of the initial state vector.
% kalman_tol [double] scalar, tolerance parameter (rcond).
% riccati_tol [double] scalar, tolerance parameter (riccati iteration).
@@ -25,7 +25,8 @@ function [LIK, lik, a, P] = missing_observations_kalman_filter(data_index,numbe
% OUTPUTS
% LIK [double] scalar, MINUS loglikelihood
% lik [double] vector, density of observations in each period.
-% a [double] mm*1 vector, estimated level of the states.
+% a [double] mm*1 vector, current estimate of the state
+% vector tomorrow (E_{T}(alpha_{T+1})).
% P [double] mm*mm matrix, covariance matrix of the states.
%
%
diff --git a/matlab/kalman/likelihood/missing_observations_kalman_filter_d.m b/matlab/kalman/likelihood/missing_observations_kalman_filter_d.m
index ec8c286a4cf67f650724ee1ddf43ce0d54d5169f..bc15d7a81a849faef9322d07f8dc130fd65e670d 100644
--- a/matlab/kalman/likelihood/missing_observations_kalman_filter_d.m
+++ b/matlab/kalman/likelihood/missing_observations_kalman_filter_d.m
@@ -12,7 +12,7 @@ function [dLIK,dlik,a,Pstar] = missing_observations_kalman_filter_d(data_index,n
% Y [double] pp*smpl matrix of (detrended) data, where pp is the number of observed variables.
% start [integer] scalar, first observation.
% last [integer] scalar, last observation.
-% a [double] mm*1 vector, levels of the state variables.
+% a [double] mm*1 vector, levels of the predicted initial state variables (E_{0}(alpha_1)).
% Pinf [double] mm*mm matrix used to initialize the covariance matrix of the state vector.
% Pstar [double] mm*mm matrix used to initialize the covariance matrix of the state vector.
% kalman_tol [double] scalar, tolerance parameter (rcond).
@@ -30,7 +30,8 @@ function [dLIK,dlik,a,Pstar] = missing_observations_kalman_filter_d(data_index,n
% OUTPUTS
% dLIK [double] scalar, MINUS loglikelihood
% dlik [double] vector, density of observations in each period.
-% a [double] mm*1 vector, estimated level of the states.
+% a [double] mm*1 vector, current estimate of the state
+% vector tomorrow (E_{T}(alpha_{T+1})).
% Pstar [double] mm*mm matrix, covariance matrix of the states.
%
% REFERENCES
diff --git a/matlab/missing_DiffuseKalmanSmootherH1_Z.m b/matlab/missing_DiffuseKalmanSmootherH1_Z.m
index 5ba614f945c059768ff3204063b8708da84e2ae0..bd3af2fa11ebd17a0fc75bb30b929160e8b96fd3 100644
--- a/matlab/missing_DiffuseKalmanSmootherH1_Z.m
+++ b/matlab/missing_DiffuseKalmanSmootherH1_Z.m
@@ -1,9 +1,10 @@
-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);
diff --git a/matlab/missing_DiffuseKalmanSmootherH3_Z.m b/matlab/missing_DiffuseKalmanSmootherH3_Z.m
index de40741fcf29f4d0c203c1d61b6c12437700e9a6..39126d57fcc0bbc020cca8a9f9739fd445f4f50a 100644
--- a/matlab/missing_DiffuseKalmanSmootherH3_Z.m
+++ b/matlab/missing_DiffuseKalmanSmootherH3_Z.m
@@ -1,9 +1,10 @@
-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);
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 671e0e1ed579487903d2275126e067609d7f36b7..a561fd54292f3423518c3bcc0b32d7e93d0c854c 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -235,6 +235,9 @@ MODFILES = \
kalman/likelihood_from_dynare/fs2000ns_corr_ME_missing.mod \
kalman/likelihood_from_dynare/fs2000ns_uncorr_ME.mod \
kalman/likelihood_from_dynare/fs2000ns_uncorr_ME_missing.mod \
+ kalman/kalman_initial_state/fs2000_smoother_only_initial_state.mod \
+ kalman/kalman_initial_state/fs2000_ns_smoother_only_initial_state.mod \
+ kalman/kalman_initial_state/fs2000_kalman_initial.mod \
second_order/burnside_1.mod \
kalman_filter_smoother/compare_results_simulation/fs2000_ML.mod \
kalman_filter_smoother/compare_results_simulation/fs2000_ML_loglinear.mod \
@@ -357,6 +360,7 @@ XFAIL_MODFILES = ramst_xfail.mod \
identification/ident_unit_root/ident_unit_root_xfail.mod \
steady_state/Linear_steady_state_xfail.mod \
optimal_policy/Ramsey/ramsey_histval_xfail.mod \
+ kalman/kalman_initial_state/fs2000_kalman_initial_xfail.mod \
example1_extra_exo_xfail.mod
MFILES = initval_file/ramst_initval_file_data.m
diff --git a/tests/kalman/kalman_initial_state/fs2000_kalman_initial.mod b/tests/kalman/kalman_initial_state/fs2000_kalman_initial.mod
new file mode 100644
index 0000000000000000000000000000000000000000..31211cb186b90a943bf2fbe8b080a73a401aba36
--- /dev/null
+++ b/tests/kalman/kalman_initial_state/fs2000_kalman_initial.mod
@@ -0,0 +1,92 @@
+// 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;
+
+steady;
+
+check;
+
+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=2000,mh_nblocks=2,mh_jscale=0.8,moments_varendo,consider_only_observed);
diff --git a/tests/kalman/kalman_initial_state/fs2000_kalman_initial_xfail.mod b/tests/kalman/kalman_initial_state/fs2000_kalman_initial_xfail.mod
new file mode 100644
index 0000000000000000000000000000000000000000..5b0f68532ba380dbc72a8db30fc5b6c9692f9054
--- /dev/null
+++ b/tests/kalman/kalman_initial_state/fs2000_kalman_initial_xfail.mod
@@ -0,0 +1,91 @@
+// 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;
+
+steady;
+
+check;
+
+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=2,mh_jscale=0.8,moments_varendo,consider_only_observed);
diff --git a/tests/kalman/kalman_initial_state/fs2000_ns_smoother_only_initial_state.mod b/tests/kalman/kalman_initial_state/fs2000_ns_smoother_only_initial_state.mod
new file mode 100644
index 0000000000000000000000000000000000000000..175dc79ea8e2f58faf1d6c9803ca87783596ffe5
--- /dev/null
+++ b/tests/kalman/kalman_initial_state/fs2000_ns_smoother_only_initial_state.mod
@@ -0,0 +1,130 @@
+/*
+ * 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
+ * following copyright notice only applies to this Dynare implementation of the
+ * model.
+ */
+
+/*
+ * Copyright (C) 2004-2010 Dynare Team
+ *
+ * This file is part of Dynare.
+ *
+ * Dynare is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Dynare is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Dynare. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+var m P c e W R k d n l gy_obs gp_obs Y_obs P_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;
+Y_obs/Y_obs(-1) = gy_obs;
+P_obs/P_obs(-1) = gp_obs;
+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;
+ Y_obs = gy_obs;
+ P_obs = gp_obs;
+end;
+
+filter_initial_state;
+P(0)=2.24;
+k(0)=5.8;
+y(0)=0.58;
+m(0)=1.01;
+Y_obs(0)=1;
+P_obs(0)=1;
+end;
+
+shocks;
+var e_a; stderr 0.014;
+var e_m; stderr 0.005;
+end;
+
+varobs P_obs Y_obs;
+
+observation_trends;
+P_obs (log(mst)-gam);
+Y_obs (gam);
+end;
+
+estimation(order=1, datafile='../../fs2000/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='../../fs2000/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='../../fs2000/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
+ * generate another random data file, comment the "estimation" line and uncomment
+ * the following lines.
+ */
+
+//stoch_simul(periods=200, order=1);
+//datatomfile('fsdat_simul', char('gy_obs', 'gp_obs'));
\ No newline at end of file
diff --git a/tests/kalman/kalman_initial_state/fs2000_smoother_only_initial_state.mod b/tests/kalman/kalman_initial_state/fs2000_smoother_only_initial_state.mod
new file mode 100644
index 0000000000000000000000000000000000000000..a8409f9d17a926b2c5733f1c033ea152af592a17
--- /dev/null
+++ b/tests/kalman/kalman_initial_state/fs2000_smoother_only_initial_state.mod
@@ -0,0 +1,123 @@
+/*
+ * 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
+ * following copyright notice only applies to this Dynare implementation of the
+ * model.
+ */
+
+/*
+ * Copyright (C) 2004-2010 Dynare Team
+ *
+ * This file is part of Dynare.
+ *
+ * Dynare is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Dynare is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Dynare. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+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;
+
+check;
+
+filter_initial_state;
+P(0)=2.24;
+k(0)=5.8;
+y(0)=0.58;
+m(0)=1.01;
+end;
+
+varobs gp_obs gy_obs;
+
+estimation(order=1, datafile='../../fs2000/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='../../fs2000/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='../../fs2000/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='../../fs2000/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='../../fs2000/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
+ * generate another random data file, comment the "estimation" line and uncomment
+ * the following lines.
+ */
+
+//stoch_simul(periods=200, order=1);
+//datatomfile('fsdat_simul', char('gy_obs', 'gp_obs'));
\ No newline at end of file