diff --git a/matlab/dr1.m b/matlab/dr1.m
index d5fbf267b..6419c4421 100644
--- a/matlab/dr1.m
+++ b/matlab/dr1.m
@@ -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
diff --git a/matlab/dyn_ramsey_static_.m b/matlab/dyn_ramsey_static.m
similarity index 56%
rename from matlab/dyn_ramsey_static_.m
rename to matlab/dyn_ramsey_static.m
index c156f2a57..4750cd2c1 100644
--- a/matlab/dyn_ramsey_static_.m
+++ b/matlab/dyn_ramsey_static.m
@@ -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 .
-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);
diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index 464264918..db69c0919 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -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
diff --git a/matlab/evaluate_static_model.m b/matlab/evaluate_static_model.m
new file mode 100644
index 000000000..29b0e3e4e
--- /dev/null
+++ b/matlab/evaluate_static_model.m
@@ -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 .
+
+ 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
diff --git a/matlab/evaluate_steady_state.m b/matlab/evaluate_steady_state.m
new file mode 100644
index 000000000..f3dc5b2fa
--- /dev/null
+++ b/matlab/evaluate_steady_state.m
@@ -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 .
+
+ 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
+
diff --git a/matlab/evaluate_steady_state_file.m b/matlab/evaluate_steady_state_file.m
new file mode 100644
index 000000000..9671ff06d
--- /dev/null
+++ b/matlab/evaluate_steady_state_file.m
@@ -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 .
+
+
+ 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
+
\ No newline at end of file
diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index 053fdddc9..740a519e9 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -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';
diff --git a/matlab/initial_estimation_checks.m b/matlab/initial_estimation_checks.m
index 02638fc76..c5bb8f1bb 100644
--- a/matlab/initial_estimation_checks.m
+++ b/matlab/initial_estimation_checks.m
@@ -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)
diff --git a/matlab/osr1.m b/matlab/osr1.m
index 4effeac90..478dff673 100644
--- a/matlab/osr1.m
+++ b/matlab/osr1.m
@@ -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;
diff --git a/matlab/print_info.m b/matlab/print_info.m
index 3874c8301..0c143226d 100644
--- a/matlab/print_info.m
+++ b/matlab/print_info.m
@@ -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
diff --git a/matlab/resol.m b/matlab/resol.m
index 1d51b542a..86d071634 100644
--- a/matlab/resol.m
+++ b/matlab/resol.m
@@ -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 .
-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 = [];
\ No newline at end of file
diff --git a/matlab/steady.m b/matlab/steady.m
index 4fd8c4fe0..db2fdd4e3 100644
--- a/matlab/steady.m
+++ b/matlab/steady.m
@@ -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_);
diff --git a/matlab/steady_.m b/matlab/steady_.m
index f090816d7..822d4fbf4 100644
--- a/matlab/steady_.m
+++ b/matlab/steady_.m
@@ -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 .
-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);
\ No newline at end of file
diff --git a/preprocessor/ExprNode.cc b/preprocessor/ExprNode.cc
index c3aa019b3..6b94fd7bb 100644
--- a/preprocessor/ExprNode.cc
+++ b/preprocessor/ExprNode.cc
@@ -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);
diff --git a/preprocessor/SteadyStateModel.cc b/preprocessor/SteadyStateModel.cc
index d665428bb..c4b2cfd28 100644
--- a/preprocessor/SteadyStateModel.cc
+++ b/preprocessor/SteadyStateModel.cc
@@ -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++)
{
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 6af009226..3319038c7 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -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 \
diff --git a/tests/auxiliary_variables/test1.mod b/tests/auxiliary_variables/test1.mod
new file mode 100644
index 000000000..c4c01c61d
--- /dev/null
+++ b/tests/auxiliary_variables/test1.mod
@@ -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;
diff --git a/tests/expectations/expectation.mod b/tests/expectations/expectation.mod
index 696eba4da..972521456 100644
--- a/tests/expectations/expectation.mod
+++ b/tests/expectations/expectation.mod
@@ -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);
diff --git a/tests/expectations/expectation_ss.mod b/tests/expectations/expectation_ss.mod
new file mode 100644
index 000000000..592905af6
--- /dev/null
+++ b/tests/expectations/expectation_ss.mod
@@ -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);
diff --git a/tests/expectations/expectation_ss_old.mod b/tests/expectations/expectation_ss_old.mod
new file mode 100644
index 000000000..3d25ba835
--- /dev/null
+++ b/tests/expectations/expectation_ss_old.mod
@@ -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);
diff --git a/tests/steady_state/walsh1_initval.mod b/tests/steady_state/walsh1_initval.mod
new file mode 100644
index 000000000..b99016724
--- /dev/null
+++ b/tests/steady_state/walsh1_initval.mod
@@ -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;
+
diff --git a/tests/steady_state/walsh1_old_ss.mod b/tests/steady_state/walsh1_old_ss.mod
new file mode 100644
index 000000000..39eb8910d
--- /dev/null
+++ b/tests/steady_state/walsh1_old_ss.mod
@@ -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;
+
diff --git a/tests/steady_state/walsh1_ssm.mod b/tests/steady_state/walsh1_ssm.mod
new file mode 100644
index 000000000..8d1de749e
--- /dev/null
+++ b/tests/steady_state/walsh1_ssm.mod
@@ -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;
+
diff --git a/tests/steady_state/walsh1_ssm_block.mod b/tests/steady_state/walsh1_ssm_block.mod
new file mode 100644
index 000000000..05bd18125
--- /dev/null
+++ b/tests/steady_state/walsh1_ssm_block.mod
@@ -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;
+