factoring out steady-state computations; steady_state_model now

generates <fname>_steadystate2.m returning parameters as well in case
they have been modified by the user. Added several test cases.
time-shift
Michel Juillard 2011-10-12 21:46:50 +02:00
parent 8f350dc744
commit ee7078e56c
24 changed files with 749 additions and 466 deletions

View File

@ -73,81 +73,6 @@ if M_.exo_nbr == 0
oo_.exo_steady_state = [] ;
end
% expanding system for Optimal Linear Regulator
if options_.ramsey_policy
if isfield(M_,'orig_model')
orig_model = M_.orig_model;
M_.endo_nbr = orig_model.endo_nbr;
M_.orig_endo_nbr = orig_model.orig_endo_nbr;
M_.aux_vars = orig_model.aux_vars;
M_.endo_names = orig_model.endo_names;
M_.lead_lag_incidence = orig_model.lead_lag_incidence;
M_.maximum_lead = orig_model.maximum_lead;
M_.maximum_endo_lead = orig_model.maximum_endo_lead;
M_.maximum_lag = orig_model.maximum_lag;
M_.maximum_endo_lag = orig_model.maximum_endo_lag;
oo_.steady_state = oo_.steady_state(1:M_.endo_nbr);
end
if options_.steadystate_flag
k_inst = [];
instruments = options_.instruments;
inst_nbr = size(options_.instruments);
for i = 1:inst_nbr
k_inst = [k_inst; strmatch(options_.instruments(i,:), ...
M_.endo_names,'exact')];
end
ys = oo_.steady_state;
if inst_nbr == 1
nl_func = @(x) dyn_ramsey_static_(x,M_,options_,oo_,it_);
% inst_val = fzero(nl_func,oo_.steady_state(k_inst));
inst_val = csolve(nl_func,oo_.steady_state(k_inst),'',options_.solve_tolf,100);
else
[inst_val,info1] = dynare_solve('dyn_ramsey_static_', ...
oo_.steady_state(k_inst),0, ...
M_,options_,oo_,it_);
end
M_.params = evalin('base','M_.params;');
ys(k_inst) = inst_val;
[x,check] = feval([M_.fname '_steadystate'],...
ys,[oo_.exo_steady_state; ...
oo_.exo_det_steady_state]);
M_.params = evalin('base','M_.params;');
if size(x,1) < M_.endo_nbr
if length(M_.aux_vars) > 0
x = add_auxiliary_variables_to_steadystate(x,M_.aux_vars,...
M_.fname,...
oo_.exo_steady_state,...
oo_.exo_det_steady_state,...
M_.params,...
options_.bytecode);
else
error([M_.fname '_steadystate.m doesn''t match the model']);
end
end
oo_.steady_state = x;
[junk,junk,multbar] = dyn_ramsey_static_(oo_.steady_state(k_inst),M_,options_,oo_,it_);
oo_.steady_state = [x(1:M_.orig_endo_nbr); multbar];
else
% xx = oo_.steady_state([1:M_.orig_endo_nbr (M_.orig_endo_nbr+M_.orig_eq_nbr+1):end]);
xx = oo_.steady_state(1:M_.orig_endo_nbr);
[xx,info1] = dynare_solve('dyn_ramsey_static_', ...
xx,0,M_,options_,oo_,it_);
[junk,junk,multbar] = dyn_ramsey_static_(xx,M_,options_,oo_,it_);
oo_.steady_state = [xx; multbar];
end
check1 = max(abs(feval([M_.fname '_static'],...
oo_.steady_state,...
[oo_.exo_steady_state; ...
oo_.exo_det_steady_state], M_.params))) > options_.dynatol ;
if check1
info(1) = 20;
info(2) = check1'*check1;
return
end
dr.ys = oo_.steady_state;
end
klen = M_.maximum_lag + M_.maximum_lead + 1;
iyv = M_.lead_lag_incidence';
iyv = iyv(:);
@ -163,24 +88,25 @@ z = repmat(dr.ys,1,klen);
if ~options_.bytecode
z = z(iyr0) ;
end;
exo_ss = [oo_.exo_steady_state' oo_.exo_det_steady_state'];
if options_.order == 1
if (options_.bytecode)
[chck, junk, loc_dr] = bytecode('dynamic','evaluate', z,[oo_.exo_simul ...
oo_.exo_det_simul], M_.params, dr.ys, 1);
[chck, junk, loc_dr] = bytecode('dynamic','evaluate', z,exo_ss, ...
M_.params, dr.ys, 1);
jacobia_ = [loc_dr.g1 loc_dr.g1_x loc_dr.g1_xd];
else
[junk,jacobia_] = feval([M_.fname '_dynamic'],z,[oo_.exo_simul ...
oo_.exo_det_simul], M_.params, dr.ys, it_);
[junk,jacobia_] = feval([M_.fname '_dynamic'],z,exo_ss, ...
M_.params, dr.ys, 1);
end;
elseif options_.order == 2
if (options_.bytecode)
[chck, junk, loc_dr] = bytecode('dynamic','evaluate', z,[oo_.exo_simul ...
oo_.exo_det_simul], M_.params, dr.ys, 1);
[chck, junk, loc_dr] = bytecode('dynamic','evaluate', z,exo_ss, ...
M_.params, dr.ys, 1);
jacobia_ = [loc_dr.g1 loc_dr.g1_x];
else
[junk,jacobia_,hessian1] = feval([M_.fname '_dynamic'],z,...
[oo_.exo_simul ...
oo_.exo_det_simul], M_.params, dr.ys, it_);
exo_ss, ...
M_.params, dr.ys, 1);
end;
if options_.use_dll
% In USE_DLL mode, the hessian is in the 3-column sparse representation

View File

@ -1,10 +1,10 @@
function [resids, rJ,mult] = dyn_ramsey_static_(x,M,options_,oo,it_)
function [steady_state,params,check] = dyn_ramsey_static(x,M,options_,oo)
% function [resids, rJ,mult] = dyn_ramsey_static_(x)
% function [steady_state,params,check] = dyn_ramsey_static_(x)
% Computes the static first order conditions for optimal policy
%
% INPUTS
% x: vector of endogenous variables
% x: vector of endogenous variables or instruments
%
% OUTPUTS
% resids: residuals of non linear equations
@ -30,13 +30,60 @@ function [resids, rJ,mult] = dyn_ramsey_static_(x,M,options_,oo,it_)
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global oo_ M_
steady_state = [];
params = M.params;
check = 0;
nl_func = @(x) dyn_ramsey_static_1(x,M,options_,oo);
if options_.steadystate_flag
k_inst = [];
instruments = options_.instruments;
inst_nbr = size(options_.instruments);
for i = 1:inst_nbr
k_inst = [k_inst; strmatch(options_.instruments(i,:), ...
M.endo_names,'exact')];
end
ys = oo.steady_state;
if inst_nbr == 1
inst_val = csolve(nl_func,oo_.steady_state(k_inst),'',options_.solve_tolf,100);
else
[inst_val,info1] = dynare_solve(nl_func,ys(k_inst),0);
end
ys(k_inst) = inst_val;
[x,params,check] = evaluate_steadystate_file(ys,exo_ss,params,M.fname,options_.steadystate_flag);
if size(x,1) < M.endo_nbr
if length(M.aux_vars) > 0
x = feval([M.fname '_set_auxiliary_variables'],xx,...
[oo.exo_steady_state,...
oo.exo_det_steady_state],...
M.params);
else
error([M.fname '_steadystate.m doesn''t match the model']);
end
end
[junk,junk,multbar] = dyn_ramsey_static_1(x(k_inst),M,options_,oo_);
steady_state = [x(1:M.orig_endo_nbr); multbar];
else
xx = oo.steady_state(1:M.orig_endo_nbr);
[xx,info1] = dynare_solve(nl_func,xx,0);
[junk,junk,multbar] = nl_func(xx);
steady_state = [xx; multbar];
end
function [resids, rJ,mult] = dyn_ramsey_static_1(x,M,options_,oo)
resids = [];
rJ = [];
mult = [];
% recovering usefull fields
endo_nbr = M.endo_nbr;
exo_nbr = M.exo_nbr;
orig_endo_nbr = M_.orig_endo_nbr;
orig_eq_nbr = M_.orig_eq_nbr;
orig_endo_nbr = M.orig_endo_nbr;
orig_eq_nbr = M.orig_eq_nbr;
inst_nbr = orig_endo_nbr - orig_eq_nbr;
% indices of Lagrange multipliers
i_mult = [orig_endo_nbr+(1:orig_eq_nbr)]';
@ -63,41 +110,26 @@ if options_.steadystate_flag
oo.steady_state,...
[oo.exo_steady_state; ...
oo.exo_det_steady_state]);
if size(x,1) < M.endo_nbr
if length(M.aux_vars) > 0
x = add_auxiliary_variables_to_steadystate(x,M.aux_vars,...
M.fname,...
oo.exo_steady_state,...
oo.exo_det_steady_state,...
M_.params,...
options_.bytecode);
else
error([M.fname '_steadystate.m doesn''t match the model']);
end
end
else
xx = zeros(endo_nbr,1);
xx(1:orig_endo_nbr) = x;
xx = feval([M_.fname '_set_auxiliary_variables'],xx,...
[oo.exo_steady_state,...
oo.exo_det_steady_state],...
M_.params);
% x = [x(1:orig_endo_nbr); zeros(orig_eq_nbr,1); x(orig_endo_nbr+1:end)];
end
% setting steady state of auxiliary variables
xx = zeros(endo_nbr,1);
xx(1:orig_endo_nbr) = x(1:orig_endo_nbr);
x = feval([M.fname '_set_auxiliary_variables'],xx,...
[oo.exo_steady_state,...
oo.exo_det_steady_state],...
M.params);
% value and Jacobian of objective function
ex = zeros(1,M.exo_nbr);
[U,Uy,Uyy] = feval([fname '_objective_static'],xx,ex, M_.params);
[U,Uy,Uyy] = feval([fname '_objective_static'],x,ex, M.params);
Uy = Uy';
Uyy = reshape(Uyy,endo_nbr,endo_nbr);
% set multipliers to 0 to compute residuals
it_ = 1;
[f,fJ] = feval([fname '_static'],xx,[oo.exo_simul oo.exo_det_simul], ...
M_.params);
[f,fJ] = feval([fname '_static'],x,[oo.exo_simul oo.exo_det_simul], ...
M.params);
aux_eq = [1:orig_endo_nbr orig_endo_nbr+orig_eq_nbr+1:size(fJ,1)];
A = fJ(aux_eq,orig_endo_nbr+1:end);

View File

@ -66,31 +66,10 @@ if ~isempty(options_.mode_file) && ~options_.mh_posterior_mode_estimation
load(options_.mode_file);
end
%% Compute the steady state:
if ~isempty(estim_params_)
set_parameters(xparam1);
end
if options_.steadystate_flag% if the *_steadystate.m file is provided.
[ys,tchek] = feval([M_.fname '_steadystate'],...
[zeros(M_.exo_nbr,1);...
oo_.exo_det_steady_state]);
if size(ys,1) < M_.endo_nbr
if length(M_.aux_vars) > 0
ys = add_auxiliary_variables_to_steadystate(ys,M_.aux_vars,...
M_.fname,...
zeros(M_.exo_nbr,1),...
oo_.exo_det_steady_state,...
M_.params,...
options_.bytecode);
else
error([M_.fname '_steadystate.m doesn''t match the model']);
end
end
oo_.steady_state = ys;
else% if the steady state file is not provided.
[dd,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
oo_.steady_state = dd.ys; clear('dd');
end
if all(abs(oo_.steady_state(bayestopt_.mfys))<1e-9)
options_.noconstant = 1;
else

View File

@ -0,0 +1,59 @@
function [residuals,check1,jacob] = evaluate_static_model(ys,exo_ss,params,M,options)
% function [ys,info] = evaluate_steady_state(M,options,oo)
% Computes the steady state
%
% INPUTS
% ys vector initial values used to compute the steady
% state
% exo_ss vector exogenous steady state
% params vector parameters
% M struct model structure
% options struct options
%
% OUTPUTS
% residuals vector residuals when ys is not
% the steady state
% check1 scalar error flag
% jacob matrix Jacobian of static model
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2001-2011 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/>.
check1 = 0;
fh_static = str2func([M.fname '_static']);
if options.bytecode
[check1, residuals] = bytecode('evaluate','static',ys,...
exo_ss, params, ys, 1);
mexErrCheck('bytecode', check1);
elseif options.block
residuals = zeros(M.endo_nbr,1);
mfs = M.blocksMFS;
for b = 1:size(mfs,1)
mfsb = mfs{b};
% blocks that can be directly evaluated (mfsb is empty)
% have zero residuals by construction
if ~isempty(mfsb)
residuals(mfsb) = feval(fh_static,b,ys,exo_ss,params);
end
end
else
residuals = feval(fh_static,ys,exo_ss,params);
end

View File

@ -0,0 +1,154 @@
function [ys,params,info] = evaluate_steady_state(ys_init,M,options,oo,steadystate_check_flag)
% function [ys,info] = evaluate_steady_state(M,options,oo)
% Computes the steady state
%
% INPUTS
% ys_init vector initial values used to compute the steady
% state
% M struct model structure
% options struct options
% oo struct output results
% steadystate_check_flag boolean if true, check that the
% steadystate verifies the
% static model
%
% OUTPUTS
% ys vector steady state
% params vector model parameters possibly
% modified by user steadystate
% function
% info 2x1 vector error codes
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2001-2011 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/>.
info = 0;
check = 0;
steadystate_flag = options.steadystate_flag;
params = M.params;
exo_ss = [oo.exo_steady_state; oo.exo_det_steady_state];
updated_params_flag = 0;
fh_static = str2func([M.fname '_static']);
if length(M.aux_vars) > 0
h_set_auxiliary_variables = str2func([M.fname '_set_auxiliary_variables']);
ys_init = h_set_auxiliary_variables(ys_init,exo_ss,M.params);
end
if options.ramsey_policy
[ys,params] = dyn_ramsey_static(ys_init,M,options,oo);
elseif steadystate_flag
% explicit steady state file
[ys,params1,check] = evaluate_steady_state_file(ys_init,exo_ss,params,M.fname,steadystate_flag);
updated_params_flag = max(abs(params1-params)) > 1e-12;
if updated_params_flag
params = params1;
end
% adding values for auxiliary variables
if length(M.aux_vars) > 0
ys = h_set_auxiliary_variables(ys,exo_ss,M.params);
end
check1 = 0;
if steadystate_check_flag
% Check whether the steady state obtained from the _steadystate file is a steady state.
[residuals,check] = evaluate_static_model(ys,exo_ss,params,M,options);
if check
info(1) = 19;
info(2) = check; % to be improved
return;
end
if max(abs(residuals)) > options.dynatol
info(1) = 19;
info(2) = residuals'*residuals;
return
end
elseif ~isempty(options.steadystate_partial)
ssvar = options.steadystate_partial.ssvar;
nov = length(ssvar);
indv = zeros(nov,1);
for i = 1:nov
indv(i) = strmatch(ssvar(i),M.endo_names,'exact');
end
[ys,check] = dynare_solve('restricted_steadystate',...
ys(indv),...
options.jacobian_flag, ...
exo_ss,indv);
end
elseif (options.bytecode == 0 && options.block == 0)
if options.linear == 0
% non linear model
[ys,check] = dynare_solve(fh_static,...
ys_init,...
options.jacobian_flag, ...
exo_ss, params);
else
% linear model
[fvec,jacob] = feval(fh_static,ys_init,exo_ss, ...
params);
if max(abs(fvec)) > 1e-12
ys = ys_init-jacob\fvec;
else
ys = ys_init;
end
end
else
% block or bytecode
[ys,check] = dynare_solve_block_or_bytecode(ys_init,exo_ss, params);
end
if check
if options.steadystate_flag
info(1)= 19;
resid = check1 ;
else
info(1)= 20;
resid = evaluate_static_model(ys_init,exo_ss,params,M,options);
end
info(2) = resid'*resid ;
return
end
if ~isreal(ys)
info(1) = 21;
info(2) = sum(imag(ys).^2);
ys = real(ys);
return
end
if ~isempty(find(isnan(ys)))
info(1) = 22;
info(2) = NaN;
return
end
if options.steadystate_flag && updated_params_flag && ~isreal(M.params)
info(1) = 23;
info(2) = sum(imag(M.params).^2);
return
end
if options.steadystate_flag && updated_params_flag && ~isempty(find(isnan(M.params)))
info(1) = 24;
info(2) = NaN;
return
end

View File

@ -0,0 +1,52 @@
function [ys,params1,check] = evaluate_steady_state_file(ys_init,exo_ss,params,fname,steadystate_flag)
% function [ys,params1,check] = evaluate_steady_state_file(ys_init,exo_ss,params,steadystate_flag)
% Evaluates steady state files
%
% INPUTS
% ys_init vector initial values used to compute the steady
% state
% exo_ss vector exogenous steady state
% params vector parameters
% steadystate_check_flag boolean if true, check that the
% steadystate verifies the
% static model
%
% OUTPUTS
% ys vector steady state
% params1 vector model parameters possibly
% modified by user steadystate
% function
% check 2x1 vector error codes
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2001-2011 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/>.
if steadystate_flag == 1
% old format
h_steadystate = str2func([fname '_steadystate']);
[ys,check] = h_steadystate(ys_init, exo_ss);
params1 = evalin('base','M_.params');
else % steadystate_flag == 2
% new format
h_steadystate = str2func([fname '_steadystate2']);
[ys,params1,check] = h_steadystate(ys_init, exo_ss, params);
end

View File

@ -60,12 +60,15 @@ options_.threads.kronecker.A_times_B_kronecker_C = 1;
options_.threads.kronecker.sparse_hessian_times_B_kronecker_C = 1;
% steady state file
if exist([M_.fname '_steadystate.m'],'file')
if exist([M_.fname '_steadystate2.m'],'file')
options_.steadystate_flag = 2;
elseif exist([M_.fname '_steadystate.m'],'file')
options_.steadystate_flag = 1;
else
options_.steadystate_flag = 0;
end
options_.steadystate_partial = [];
options_.steadystate.nocheck = 0;
% subset of the estimated deep parameters
options_.ParamSubSet = 'None';

View File

@ -3,12 +3,12 @@ function DynareResults = initial_estimation_checks(xparam1,DynareDataset,Model,E
% Checks data (complex values, ML evaluation, initial values, BK conditions,..)
%
% INPUTS
% xparam1: vector of parameters to be estimated
% gend: scalar specifying the number of observations
% data: matrix of data
% xparam1: vector of parameters to be estimated
% gend: scalar specifying the number of observations
% data: matrix of data
%
% OUTPUTS
% none
% DynareResults structure of temporary results
%
% SPECIAL REQUIREMENTS
% none
@ -34,56 +34,16 @@ if DynareDataset.info.nvobs>Model.exo_nbr+EstimatedParameters.nvn
error(['initial_estimation_checks:: Estimation can''t take place because there are less shocks than observed variables!'])
end
% check if steady state solves static model (except if diffuse_filter == 1)
[DynareResults.steady_state] = ...
evaluate_steady_state(DynareResults.steady_state,Model,DynareOptions,DynareResults,DynareOptions.diffuse_filter==0);
if DynareOptions.dsge_var
[fval,cost_flag,info] = DsgeVarLikelihood(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
else
[fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
end
% when their is an analytical steadystate, check that the values
% returned by *_steadystate match with the static model
if DynareOptions.steadystate_flag
[ys,check] = feval([Model.fname '_steadystate'],...
DynareResults.steady_state,...
[DynareResults.exo_steady_state; ...
DynareResults.exo_det_steady_state]);
if size(ys,1) < Model.endo_nbr
if length(Model.aux_vars) > 0
ys = add_auxiliary_variables_to_steadystate(ys,Model.aux_vars,...
Model.fname,...
DynareResults.exo_steady_state,...
DynareResults.exo_det_steady_state,...
Model.params,...
DynareOptions.bytecode);
else
error([Model.fname '_steadystate.m doesn''t match the model']);
end
end
DynareResults.steady_state = ys;
% Check if the steady state obtained from the _steadystate file is a
% steady state.
check1 = 0;
if isfield(DynareOptions,'unit_root_vars') && DynareOptions.diffuse_filter == 0
if isempty(DynareOptions.unit_root_vars)
if ~DynareOptions.bytecode
check1 = max(abs(feval([Model.fname '_static'],...
DynareResults.steady_state,...
[DynareResults.exo_steady_state; ...
DynareResults.exo_det_steady_state], Model.params))) > DynareOptions.dynatol ;
else
[info, res] = bytecode('static','evaluate',DynareResults.steady_state,...
[DynareResults.exo_steady_state; ...
DynareResults.exo_det_steady_state], Model.params);
check1 = max(abs(res)) > DynareOptions.dynatol;
end
if check1
error(['The seadystate values returned by ' Model.fname ...
'_steadystate.m don''t solve the static model!' ])
end
end
end
end
if info(1) > 0
disp('Error in computing likelihood for initial parameter values')
print_info(info, DynareOptions.noprint)

View File

@ -42,54 +42,10 @@ exe =zeros(M_.exo_nbr,1);
dr = set_state_space(oo_.dr,M_);
% check if ys is steady state
if exist([M_.fname '_steadystate'])
[ys,check1] = feval([M_.fname '_steadystate'],oo_.steady_state,...
[oo_.exo_steady_state; oo_.exo_det_steady_state]);
if size(ys,1) < M_.endo_nbr
if length(M_.aux_vars) > 0
ys = add_auxiliary_variables_to_steadystate(ys,M_.aux_vars,...
M_.fname,...
oo_.exo_steady_state,...
oo_.exo_det_steady_state,M_.params,...
options_.bytecode);
else
error([M_.fname '_steadystate.m doesn''t match the model']);
end
end
dr.ys = ys;
else
% testing if ys isn't a steady state or if we aren't computing Ramsey policy
fh = str2func([M_.fname '_static']);
if max(abs(feval(fh,oo_.steady_state,[oo_.exo_steady_state; oo_.exo_det_steady_state], M_.params))) ...
> options_.dynatol && options_.ramsey_policy == 0
if options_.linear == 0
% nonlinear models
[dr.ys,check1] = dynare_solve(fh,dr.ys,options_.jacobian_flag,...
[oo_.exo_steady_state; ...
oo_.exo_det_steady_state], M_.params);
else
% linear models
[fvec,jacob] = feval(fh,oo_.steady_state,[oo_.exo_steady_state;...
oo_.exo_det_steady_state], M_.params);
dr.ys = oo_.steady_state-jacob\fvec;
end
end
end
oo_.dr = dr;
% $$$ if max(abs(feval(fh, oo_.steady_state, exe, M_.params))) > options_.dynatol
% $$$ [oo_.dr.ys, check] = dynare_solve([M_.fname '_static'], oo_.steady_state, 1, exe, M_.params);
% $$$ if check
% $$$ error('OLR: convergence problem in DYNARE_SOLVE')
% $$$ end
% $$$ else
% $$$ oo_.dr.ys = oo_.steady_state;
% $$$ end
np = size(i_params,1);
t0 = M_.params(i_params);
inv_order_var = oo_.dr.inv_order_var;
inv_order_var = dr.inv_order_var;
H0 = 1e-4*eye(np);
crit = 1e-7;

View File

@ -60,11 +60,11 @@ if ~noprint
case 21
error('The steady state is complex.')
case 22
error('The steady state contains NaN.')
error('The steady state contains NaN or Inf.')
case 23
error('Some updated params are complex.')
case 24
error('Some updated params contain NaN.')
error('Some updated params contain NaN or Inf.')
case 30
error('Variance can''t be computed')
case 41

View File

@ -93,10 +93,6 @@ function [dr,info,M,options,oo] = resol(check_flag,M,options,oo)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global it_
jacobian_flag = 0;
if isfield(oo,'dr');
dr = oo.dr;
end
@ -110,136 +106,19 @@ if M.exo_nbr == 0
oo.exo_steady_state = [] ;
end
params0 = M.params;
[dr.ys,M.params,info] = evaluate_steady_state(oo.steady_state,M,options,oo,0);
% check if steady_state_0 (-> oo.steady_state) is steady state
tempex = oo.exo_simul;
oo.exo_simul = repmat(oo.exo_steady_state',M.maximum_lag+M.maximum_lead+1,1);
if M.exo_det_nbr > 0
tempexdet = oo.exo_det_simul;
oo.exo_det_simul = repmat(oo.exo_det_steady_state',M.maximum_lag+M.maximum_lead+1,1);
end
steady_state = oo.steady_state;
check1 = 0;
% testing for steadystate file
if (~options.bytecode)
fh = [M.fname '_static'];
end
if options.steadystate_flag
[steady_state,check1] = feval([M.fname '_steadystate'],steady_state,...
[oo.exo_steady_state; ...
oo.exo_det_steady_state]);
if size(steady_state,1) < M.endo_nbr
if length(M.aux_vars) > 0
steady_state = add_auxiliary_variables_to_steadystate(steady_state,M.aux_vars,...
M.fname,...
oo.exo_steady_state,...
oo.exo_det_steady_state,...
M.params,...
options.bytecode);
else
error([M.fname '_steadystate.m doesn''t match the model']);
end
end
else
% testing if steady_state_0 (-> oo.steady_state) isn't a steady state or if we aren't computing Ramsey policy
if options.ramsey_policy == 0
if options.linear == 0
% nonlinear models
if (options.block == 0 && options.bytecode == 0)
if max(abs(feval(fh,steady_state,[oo.exo_steady_state; ...
oo.exo_det_steady_state], M.params))) > options.dynatol
[steady_state,check1] = dynare_solve(fh,steady_state,options.jacobian_flag,...
[oo.exo_steady_state; ...
oo.exo_det_steady_state], M.params);
end
else
[steady_state,check1] = dynare_solve_block_or_bytecode(steady_state,...
[oo.exo_steady_state; ...
oo.exo_det_steady_state], M.params);
end;
else
if (options.block == 0 && options.bytecode == 0)
% linear models
[fvec,jacob] = feval(fh,steady_state,[oo.exo_steady_state;...
oo.exo_det_steady_state], M.params);
if max(abs(fvec)) > 1e-12
steady_state = steady_state-jacob\fvec;
end
else
[steady_state,check1] = dynare_solve_block_or_bytecode(steady_state,...
[oo.exo_steady_state; ...
oo.exo_det_steady_state], M.params);
end;
end
end
end
% test if M.params_has changed.
if options.steadystate_flag
updated_params_flag = max(abs(M.params-params0))>1e-12;
else
updated_params_flag = 0;
end
% testing for problem.
dr.ys = steady_state;
if check1
if options.steadystate_flag
info(1)= 19;
resid = check1 ;
else
info(1)= 20;
resid = feval(fh,oo.steady_state,oo.exo_steady_state, M.params);
end
info(2) = resid'*resid ;
if info(1)
return
end
if ~isreal(steady_state)
info(1) = 21;
info(2) = sum(imag(steady_state).^2);
steady_state = real(steady_state);
dr.ys = steady_state;
return
end
if ~isempty(find(isnan(steady_state)))
info(1) = 22;
info(2) = NaN;
dr.ys = steady_state;
return
end
if options.steadystate_flag && updated_params_flag && ~isreal(M.params)
info(1) = 23;
info(2) = sum(imag(M.params).^2);
dr.ys = steady_state;
return
end
if options.steadystate_flag && updated_params_flag && ~isempty(find(isnan(M.params)))
info(1) = 24;
info(2) = NaN;
dr.ys = steady_state;
return
end
if options.block
[dr,info,M,options,oo] = dr_block(dr,check_flag,M,options,oo);
else
[dr,info,M,options,oo] = dr1(dr,check_flag,M,options,oo);
end
if info(1)
return
end
if M.exo_det_nbr > 0
oo.exo_det_simul = tempexdet;
end
oo.exo_simul = tempex;
tempex = [];

View File

@ -34,9 +34,6 @@ test_for_deep_parameters_calibration(M_);
options_ = set_default_option(options_,'jacobian_flag',1);
options_ = set_default_option(options_,'steadystate_flag',0);
if exist([M_.fname '_steadystate.m'])
options_.steadystate_flag = 1;
end
if options_.steadystate_flag && options_.homotopy_mode
error('STEADY: Can''t use homotopy when providing a steady state external file');
@ -51,7 +48,11 @@ switch options_.homotopy_mode
homotopy3(options_.homotopy_values, options_.homotopy_steps);
end
steady_;
[oo_.steady_state,M_.params,info] = steady_(M_,options_,oo_);
if info(1)
print_info(info,options_.noprint);
end
disp_steady_state(M_,oo_);

View File

@ -1,12 +1,19 @@
function steady_()
% function steady_()
function [steady_state,params,info] = steady_(M_,options_,oo_)
% function [steady_state,params,info] = steady_(M_,options_,oo_)
% Computes the steady state
%
% INPUTS
% none
% M struct model structure
% options struct options
% oo struct output results
%
% OUTPUTS
% none
% steady_state vector steady state
% params vector parameters (may have been
% modified by user in
% explicit computation of
% the steady state)
% info 2x1 vector error codes
%
% SPECIAL REQUIREMENTS
% none
@ -28,8 +35,6 @@ function steady_()
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global M_ oo_ options_
if options_.solve_algo < 0 || options_.solve_algo > 8
error('STEADY: solve_algo must be between 0 and 8')
end
@ -46,88 +51,4 @@ if exist('OCTAVE_VERSION') && options_.solve_algo == 7
error('SIMUL: you can''t use solve_algo = 7 under Octave')
end
if options_.steadystate_flag
[ys,check] = feval([M_.fname '_steadystate'],...
oo_.steady_state,...
[oo_.exo_steady_state; ...
oo_.exo_det_steady_state]);
if size(ys,1) < M_.endo_nbr
if length(M_.aux_vars) > 0
ys = add_auxiliary_variables_to_steadystate(ys,M_.aux_vars,...
M_.fname,...
oo_.exo_steady_state,...
oo_.exo_det_steady_state,...
M_.params,...
options_.bytecode);
else
error([M_.fname '_steadystate.m doesn''t match the model']);
end
end
oo_.steady_state = ys;
% Check if the steady state obtained from the _steadystate file is a steady state.
check1 = 0;
if isempty(options_.unit_root_vars)
if options_.block && ~options_.bytecode
check2 = zeros(size(M_.blocksMFS,1),1);
for b = 1:size(M_.blocksMFS,1)
n = size(M_.blocksMFS{b}, 1);
ss = oo_.steady_state;
if n
check2(b) = max(abs(feval([M_.fname '_static'], b, ss, ...
[oo_.exo_steady_state; oo_.exo_det_steady_state], M_.params))) > options_.dynatol;
else
[r, g1, ssss] = feval([M_.fname '_static'], b, ss, ...
[oo_.exo_steady_state; oo_.exo_det_steady_state], M_.params);
end
end
check1 = any(check2);
idx = find(abs(ssss-oo_.steady_state)>10*options_.dynatol);
if ~isempty(idx)
check1 = 1;
end
elseif options_.bytecode
[check1, residuals] = bytecode('evaluate','static',oo_.steady_state,...
[oo_.exo_steady_state; ...
oo_.exo_det_steady_state], M_.params, oo_.steady_state, 1);
mexErrCheck('bytecode', check1);
check1 = max(abs(residuals)) > options_.dynatol ;
else
check1 = 0;
check1 = max(abs(feval([M_.fname '_static'],...
oo_.steady_state,...
[oo_.exo_steady_state; ...
oo_.exo_det_steady_state], M_.params))) > options_.dynatol ;
end
end
if check1
resid;
error(['The steadystate values returned by ' M_.fname ...
'_steadystate.m don''t solve the static model!' ])
end
if ~isempty(options_.steadystate_partial)
ssvar = options_.steadystate_partial.ssvar;
nov = length(ssvar);
indv = zeros(nov,1);
for i = 1:nov
indv(i) = strmatch(ssvar(i),M_.endo_names,'exact');
end
[oo_.steady_state,check] = dynare_solve('restricted_steadystate',...
oo_.steady_state(indv),...
options_.jacobian_flag, ...
[oo_.exo_steady_state;oo_.exo_det_steady_state],indv);
end
elseif (options_.bytecode == 0 && options_.block == 0)
[oo_.steady_state,check] = dynare_solve([M_.fname '_static'],...
oo_.steady_state,...
options_.jacobian_flag, ...
[oo_.exo_steady_state; ...
oo_.exo_det_steady_state], M_.params);
else
[oo_.steady_state,check] = dynare_solve_block_or_bytecode(oo_.steady_state,...
[oo_.exo_steady_state; ...
oo_.exo_det_steady_state], M_.params);
end
if check ~= 0
error('STEADY: convergence problems')
end
[steady_state,params,info] = evaluate_steady_state(oo_.steady_state,M_,options_,oo_,~options_.steadystate.nocheck);

View File

@ -584,7 +584,7 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
switch (type)
{
case eParameter:
if (output_type == oMatlabOutsideModel || output_type == oSteadyStateFile)
if (output_type == oMatlabOutsideModel)
output << "M_.params" << "(" << tsid + 1 << ")";
else
output << "params" << LEFT_ARRAY_SUBSCRIPT(output_type) << tsid + ARRAY_SUBSCRIPT_OFFSET(output_type) << RIGHT_ARRAY_SUBSCRIPT(output_type);

View File

@ -103,7 +103,7 @@ SteadyStateModel::writeSteadyStateFile(const string &basename, bool ramsey_polic
if (recursive_order.size() == 0)
return;
string filename = basename + "_steadystate.m";
string filename = basename + "_steadystate2.m";
ofstream output;
output.open(filename.c_str(), ios::out | ios::binary);
@ -113,16 +113,11 @@ SteadyStateModel::writeSteadyStateFile(const string &basename, bool ramsey_polic
exit(EXIT_FAILURE);
}
output << "function [ys_, check_] = " << basename << "_steadystate(";
if (ramsey_policy)
output << "ys_";
else
output << "ys_orig_";
output << ", exo_)" << endl
<< "% Steady state generated by Dynare preprocessor" << endl;
if (!ramsey_policy)
output << " ys_=zeros(" << symbol_table.endo_nbr() << ",1);" << endl;
output << " global M_" << endl;
output << "function [ys_, params, info] = " << basename << "_steadystate2("
<< "ys_, exo_, params)" << endl
<< "% Steady state generated by Dynare preprocessor" << endl
<< " info = 0;" << endl
<< " ys_ = zeros(" << symbol_table.endo_nbr() << ",1);" << endl;
for (size_t i = 0; i < recursive_order.size(); i++)
{

View File

@ -21,7 +21,14 @@ MODFILES = \
comments.mod \
histval_sto.mod \
histval_det.mod \
auxiliary_variables/test1.mod \
expectations/expectation.mod \
expectations/expectation_ss.mod \
expectations/expectation_ss_old.mod \
steady_state/walsh1_initval.mod \
steady_state/walsh1_old_ss.mod \
steady_state/walsh1_ssm.mod \
steady_state/walsh1_ssm_block.mod \
steady_state_operator/standard.mod \
steady_state_operator/use_dll.mod \
steady_state_operator/block.mod \

View File

@ -0,0 +1,14 @@
var x, y;
varexo e;
model;
expectation(-2)(x(+4)*y(+1)) = 1;
y = exp(e);
end;
initval;
y = 1;
x = 1;
end;
steady;

View File

@ -41,4 +41,4 @@ var u; stderr 0.009;
var e, u = phi*0.009*0.009;
end;
stoch_simul;
stoch_simul(order=1,irf=0);

View File

@ -0,0 +1,46 @@
// Example 1 from Collard's guide to Dynare
var y, c, k, a, h, b;
varexo e, u;
parameters beta, rho, alpha, delta, theta, psi, tau;
alpha = 0.36;
rho = 0.95;
tau = 0.025;
beta = 0.99;
delta = 0.025;
psi = 0;
theta = 2.95;
phi = 0.1;
model;
c*theta*h^(1+psi)=(expectation(1)((1-alpha)*y)+expectation(-2)((1-alpha)*y))/2;
k = beta*(((exp(b)*c)/(exp(b(+1))*c(+1)))
*(exp(b(+1))*alpha*y(+1)+(1-delta)*k));
y = exp(a)*(k(-1)^alpha)*(h^(1-alpha));
k = exp(b)*(y-c)+(1-delta)*k(-1);
a = rho*a(-1)+tau*b(-1) + e;
b = tau*a(-1)+rho*b(-1) + u;
end;
steady_state_model;
a = 0;
b = 0;
h = 1/3;
k = ((1/beta-(1-delta))/(alpha*h^(1-alpha)))^(1/(alpha-1));
y = k^alpha*h^(1-alpha);
c = y - delta*k;
theta = (1-alpha)*y/(c*h^(1+psi));
end;
steady;
shocks;
var e; stderr 0.009;
var u; stderr 0.009;
var e, u = phi*0.009*0.009;
end;
stoch_simul(order=1,irf=0);

View File

@ -0,0 +1,36 @@
// Example 1 from Collard's guide to Dynare
var y, c, k, a, h, b;
varexo e, u;
parameters beta, rho, alpha, delta, theta, psi, tau;
alpha = 0.36;
rho = 0.95;
tau = 0.025;
beta = 0.99;
delta = 0.025;
psi = 0;
theta = 2.95;
phi = 0.1;
model;
c*theta*h^(1+psi)=(expectation(1)((1-alpha)*y)+expectation(-2)((1-alpha)*y))/2;
k = beta*(((exp(b)*c)/(exp(b(+1))*c(+1)))
*(exp(b(+1))*alpha*y(+1)+(1-delta)*k));
y = exp(a)*(k(-1)^alpha)*(h^(1-alpha));
k = exp(b)*(y-c)+(1-delta)*k(-1);
a = rho*a(-1)+tau*b(-1) + e;
b = tau*a(-1)+rho*b(-1) + u;
end;
steady;
shocks;
var e; stderr 0.009;
var u; stderr 0.009;
var e, u = phi*0.009*0.009;
end;
stoch_simul(order=1,irf=0);

View File

@ -0,0 +1,68 @@
var y c k m n R pi z u;
varexo e sigma;
// sigma stands for phi in the eq 2.37 p.69
parameters alpha beta delta gamm phi1 eta a b rho phi2 Psi thetass;
//phi1 stands for capital phi in eq.2.68 and 2.69
//phi2 stands for lowercase phi in eq. 2.66
change_type(var) Psi;
change_type(parameters) n;
alpha = 0.36;
beta = 0.989;
gamm = 0.5;
delta = 0.019;
phi1 = 2;
phi2 = 0;
eta = 1;
a = 0.95;
b = 2.56;
rho = 0.95;
//Psi = 1.47630583;
thetass = 1.0125;
model;
(a*exp(c)^(1-b)+(1-a)*exp(m)^(1-b))^((b-phi1)/(1-b))*a*exp(c)^(-b) = (a*exp(c)^(1-b)+(1-a)*exp(m)^(1-b))^((b-phi1)/(1-b))*(1-a)*exp(m)^(-b)+beta*(a*exp(c(+1))^(1-b)+(1-a)*exp(m(+1))^(1-b))^((b-phi1)/(1-b))*a*exp(c(+1))^(-b)/(1+pi(+1));
Psi*(1-exp(n))^(-eta)/(a*exp(c)^(-b)*(a*exp(c)^(1-b) + (1-a)*exp(m)^(1-b))^((b-phi1)/(1-b))) = (1-alpha)*exp(y)/exp(n);
(a*exp(c)^(1-b)+(1-a)*exp(m)^(1-b))^((b-phi1)/(1-b))*a*exp(c)^(-b) = beta*exp(R(+1))*(a*exp(c(+1))^(1-b)+(1-a)*exp(m(+1))^(1-b))^((b-phi1)/(1-b))*a*exp(c(+1))^(-b);
exp(R) = alpha*exp(y)/exp(k(-1)) + 1-delta;
exp(k) = (1-delta)*exp(k(-1))+exp(y)-exp(c);
exp(y) = exp(z)*exp(k(-1))^alpha*exp(n)^(1-alpha);
exp(m) = exp(m(-1))*(u+thetass)/(1+pi);
z = rho*z(-1) + e;
u = gamm*u(-1) + phi2*z(-1) + sigma;
end;
shocks;
var e; stderr 0.007;
var sigma;stderr 0.0089;
end;
n=log(1/3);
initval;
y = 0.2;
c = 0.03;
k = 2.7;
m = 0.3;
Psi = 1.5;
R = 0.01;
pi = 0.01;
z = 0;
u = 0;
end;
steady;

View File

@ -0,0 +1,51 @@
var y c k m n R pi z u;
varexo e sigma;
// sigma stands for phi in the eq 2.37 p.69
parameters alpha beta delta gamm phi1 eta a b rho phi2 Psi thetass;
//phi1 stands for capital phi in eq.2.68 and 2.69
//phi2 stands for lowercase phi in eq. 2.66
alpha = 0.36;
beta = 0.989;
gamm = 0.5;
delta = 0.019;
phi1 = 2;
phi2 = 0;
eta = 1;
a = 0.95;
b = 2.56;
rho = 0.95;
Psi = 1.47630583;
thetass = 1.0125;
model;
(a*exp(c)^(1-b)+(1-a)*exp(m)^(1-b))^((b-phi1)/(1-b))*a*exp(c)^(-b) = (a*exp(c)^(1-b)+(1-a)*exp(m)^(1-b))^((b-phi1)/(1-b))*(1-a)*exp(m)^(-b)+beta*(a*exp(c(+1))^(1-b)+(1-a)*exp(m(+1))^(1-b))^((b-phi1)/(1-b))*a*exp(c(+1))^(-b)/(1+pi(+1));
Psi*(1-exp(n))^(-eta)/(a*exp(c)^(-b)*(a*exp(c)^(1-b) + (1-a)*exp(m)^(1-b))^((b-phi1)/(1-b))) = (1-alpha)*exp(y)/exp(n);
(a*exp(c)^(1-b)+(1-a)*exp(m)^(1-b))^((b-phi1)/(1-b))*a*exp(c)^(-b) = beta*exp(R(+1))*(a*exp(c(+1))^(1-b)+(1-a)*exp(m(+1))^(1-b))^((b-phi1)/(1-b))*a*exp(c(+1))^(-b);
exp(R) = alpha*exp(y)/exp(k(-1)) + 1-delta;
exp(k) = (1-delta)*exp(k(-1))+exp(y)-exp(c);
exp(y) = exp(z)*exp(k(-1))^alpha*exp(n)^(1-alpha);
exp(m) = exp(m(-1))*(u+thetass)/(1+pi);
z = rho*z(-1) + e;
u = gamm*u(-1) + phi2*z(-1) + sigma;
end;
shocks;
var e; stderr 0.007;
var sigma;stderr 0.0089;
end;
steady;

View File

@ -0,0 +1,72 @@
var y c k m n R pi z u;
varexo e sigma;
// sigma stands for phi in the eq 2.37 p.69
parameters alpha beta delta gamm phi1 eta a b rho phi2 Psi thetass;
//phi1 stands for capital phi in eq.2.68 and 2.69
//phi2 stands for lowercase phi in eq. 2.66
alpha = 0.36;
beta = 0.989;
gamm = 0.5;
delta = 0.019;
phi1 = 2;
phi2 = 0;
eta = 1;
a = 0.95;
b = 2.56;
rho = 0.95;
Psi = 1.47630583;
thetass = 1.0125;
model;
(a*exp(c)^(1-b)+(1-a)*exp(m)^(1-b))^((b-phi1)/(1-b))*a*exp(c)^(-b) = (a*exp(c)^(1-b)+(1-a)*exp(m)^(1-b))^((b-phi1)/(1-b))*(1-a)*exp(m)^(-b)+beta*(a*exp(c(+1))^(1-b)+(1-a)*exp(m(+1))^(1-b))^((b-phi1)/(1-b))*a*exp(c(+1))^(-b)/(1+pi(+1));
Psi*(1-exp(n))^(-eta)/(a*exp(c)^(-b)*(a*exp(c)^(1-b) + (1-a)*exp(m)^(1-b))^((b-phi1)/(1-b))) = (1-alpha)*exp(y)/exp(n);
(a*exp(c)^(1-b)+(1-a)*exp(m)^(1-b))^((b-phi1)/(1-b))*a*exp(c)^(-b) = beta*exp(R(+1))*(a*exp(c(+1))^(1-b)+(1-a)*exp(m(+1))^(1-b))^((b-phi1)/(1-b))*a*exp(c(+1))^(-b);
exp(R) = alpha*exp(y)/exp(k(-1)) + 1-delta;
exp(k) = (1-delta)*exp(k(-1))+exp(y)-exp(c);
exp(y) = exp(z)*exp(k(-1))^alpha*exp(n)^(1-alpha);
exp(m) = exp(m(-1))*(u+thetass)/(1+pi);
z = rho*z(-1) + e;
u = gamm*u(-1) + phi2*z(-1) + sigma;
end;
shocks;
var e; stderr 0.007;
var sigma;stderr 0.0089;
end;
steady_state_model;
// solving in levels
// calibrating n = 1/3 and recovering the value of Psi
// adapting solution Walsh (2003) p. 84
pi = thetass-1;
en = 1/3;
eR = 1/beta;
y_k = (1/alpha)*(1/beta-1+delta);
ek = en*y_k^(-1/(1-alpha));
ec = ek*(y_k-delta);
em = ec*(a/(1-a))^(-1/b)*((thetass-beta)/thetass)^(-1/b);
ey = ek*y_k;
Xss = a*ec^(1-b)*(1+(a/(1-a))^(-1/b)*((thetass-beta)/thetass)^((b-1)/b));
Psi = (1-alpha)*(ey/en)*Xss^((b-phi1)/(1-b))*a*ec^(-b)*(1-en)^eta;
n = log(en);
k = log(ek);
m = log(em);
c = log(ec);
y = log(ey);
R = log(eR);
end;
steady;

View File

@ -0,0 +1,72 @@
var y c k m n R pi z u;
varexo e sigma;
// sigma stands for phi in the eq 2.37 p.69
parameters alpha beta delta gamm phi1 eta a b rho phi2 Psi thetass;
//phi1 stands for capital phi in eq.2.68 and 2.69
//phi2 stands for lowercase phi in eq. 2.66
alpha = 0.36;
beta = 0.989;
gamm = 0.5;
delta = 0.019;
phi1 = 2;
phi2 = 0;
eta = 1;
a = 0.95;
b = 2.56;
rho = 0.95;
Psi = 1.47630583;
thetass = 1.0125;
model(block);
(a*exp(c)^(1-b)+(1-a)*exp(m)^(1-b))^((b-phi1)/(1-b))*a*exp(c)^(-b) = (a*exp(c)^(1-b)+(1-a)*exp(m)^(1-b))^((b-phi1)/(1-b))*(1-a)*exp(m)^(-b)+beta*(a*exp(c(+1))^(1-b)+(1-a)*exp(m(+1))^(1-b))^((b-phi1)/(1-b))*a*exp(c(+1))^(-b)/(1+pi(+1));
Psi*(1-exp(n))^(-eta)/(a*exp(c)^(-b)*(a*exp(c)^(1-b) + (1-a)*exp(m)^(1-b))^((b-phi1)/(1-b))) = (1-alpha)*exp(y)/exp(n);
(a*exp(c)^(1-b)+(1-a)*exp(m)^(1-b))^((b-phi1)/(1-b))*a*exp(c)^(-b) = beta*exp(R(+1))*(a*exp(c(+1))^(1-b)+(1-a)*exp(m(+1))^(1-b))^((b-phi1)/(1-b))*a*exp(c(+1))^(-b);
exp(R) = alpha*exp(y)/exp(k(-1)) + 1-delta;
exp(k) = (1-delta)*exp(k(-1))+exp(y)-exp(c);
exp(y) = exp(z)*exp(k(-1))^alpha*exp(n)^(1-alpha);
exp(m) = exp(m(-1))*(u+thetass)/(1+pi);
z = rho*z(-1) + e;
u = gamm*u(-1) + phi2*z(-1) + sigma;
end;
shocks;
var e; stderr 0.007;
var sigma;stderr 0.0089;
end;
steady_state_model;
// solving in levels
// calibrating n = 1/3 and recovering the value of Psi
// adapting solution Walsh (2003) p. 84
pi = thetass-1;
en = 1/3;
eR = 1/beta;
y_k = (1/alpha)*(1/beta-1+delta);
ek = en*y_k^(-1/(1-alpha));
ec = ek*(y_k-delta);
em = ec*(a/(1-a))^(-1/b)*((thetass-beta)/thetass)^(-1/b);
ey = ek*y_k;
Xss = a*ec^(1-b)*(1+(a/(1-a))^(-1/b)*((thetass-beta)/thetass)^((b-1)/b));
Psi = (1-alpha)*(ey/en)*Xss^((b-phi1)/(1-b))*a*ec^(-b)*(1-en)^eta;
n = log(en);
k = log(ek);
m = log(em);
c = log(ec);
y = log(ey);
R = log(eR);
end;
steady;