stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'mcp_steady' into 'master' 

Provisions for solving steady state with MCP-tag

See merge request Dynare/dynare!1877
bgp-dev
Sébastien Villemot 2022-09-13 09:26:08 +00:00
parent b071e85214 3fbab4cbe0
commit 9b787d8417
5 changed files with 122 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
matlab/dynare_solve.m
View File

@ -354,12 +354,14 @@ elseif options.solve_algo==10
% LMMCP
olmmcp = options.lmmcp;
[x, fvec, errorcode, ~, fjac] = lmmcp(f, x, olmmcp.lb, olmmcp.ub, olmmcp, args{:});
eq_to_check=find(isfinite(olmmcp.lb) | isfinite(olmmcp.ub));
eq_to_ignore=eq_to_check(x(eq_to_check,:)<=olmmcp.lb(eq_to_check)+eps | x(eq_to_check,:)>=olmmcp.ub(eq_to_check)-eps);
fvec(eq_to_ignore)=0;
if errorcode==1
errorflag = false;
else
errorflag = true;
end
[fvec, fjac] = feval(f, x, args{:});
elseif options.solve_algo == 11
% PATH mixed complementary problem
% PATH linear mixed complementary problem
@ -378,7 +380,9 @@ elseif options.solve_algo == 11
errorflag = true;
end
errorcode = nan; % There is no error code for this algorithm, as PATH is closed source it is unlikely we can fix that.
[fvec, fjac] = feval(f, x, args{:});
eq_to_check=find(isfinite(omcppath.lb) | isfinite(omcppath.ub));
eq_to_ignore=eq_to_check(x(eq_to_check,:)<=omcppath.lb(eq_to_check)+eps | x(eq_to_check,:)>=omcppath.ub(eq_to_check)-eps);
fvec(eq_to_ignore)=0;
elseif ismember(options.solve_algo, [13, 14])
if ~jacobian_flag
error('DYNARE_SOLVE: option solve_algo=13|14 needs computed Jacobian')

30
matlab/evaluate_steady_state.m
View File

@ -280,11 +280,27 @@ elseif steadystate_flag
end
elseif ~options.bytecode && ~options.block
if ~options.linear
% non linear model
static_model = str2func(sprintf('%s.static', M.fname));
[ys, check] = dynare_solve(@static_problem, ys_init, ...
options.steady.maxit, options.solve_tolf, options.solve_tolx, ...
options, exo_ss, params, M.endo_nbr, static_model);
% non linear model
if ismember(options.solve_algo,[10,11])
[lb,ub,eq_index] = get_complementarity_conditions(M,options.ramsey_policy);
if options.solve_algo == 10
options.lmmcp.lb = lb;
options.lmmcp.ub = ub;
elseif options.solve_algo == 11
options.mcppath.lb = lb;
options.mcppath.ub = ub;
end
[ys,check,fvec] = dynare_solve(@static_mcp_problem,...
ys_init,...
options.steady.maxit, options.solve_tolf, options.solve_tolx, ...
options, exo_ss, params,...
M.endo_nbr,static_model,eq_index);
else
[ys, check] = dynare_solve(@static_problem, ys_init, ...
options.steady.maxit, options.solve_tolf, options.solve_tolx, ...
options, exo_ss, params, M.endo_nbr, static_model);
end
if check && options.debug
[ys, check, fvec, fjac, errorcode] = dynare_solve(@static_problem, ys_init, ...
options.steady.maxit, options.solve_tolf, options.solve_tolx, ...
@ -412,3 +428,9 @@ function [resids,jac] = static_problem(y,x,params,nvar,fh_static_model)
[r,j] = fh_static_model(y,x,params);
resids = r(1:nvar);
jac = j(1:nvar,1:nvar);
function [resids,jac] = static_mcp_problem(y,x,params,nvar,fh_static_model,eq_index)
[r,j] = fh_static_model(y,x,params);
resids = r(eq_index);
jac = j(eq_index,1:nvar);

19
matlab/resid.m
View File

@ -37,7 +37,7 @@ non_zero = nargin > 0 && isfield(options_resid_, 'non_zero') && options_resid_.n
tags = M_.equations_tags;
istag = 0;
if length(tags)
if ~isempty(tags)
istag = 1;
end
@ -70,14 +70,27 @@ z = evaluate_static_model(oo_.steady_state, [oo_.exo_steady_state; ...
oo_.exo_det_steady_state], ...
M_.params, M_, options_);
if ismember(options_.solve_algo,[10,11])
[lb,ub,eq_index] = get_complementarity_conditions(M_,options_.ramsey_policy);
eq_to_check=find(isfinite(lb) | isfinite(ub));
eq_to_ignore=eq_to_check(oo_.steady_state(eq_to_check,:)<=lb(eq_to_check)+eps | oo_.steady_state(eq_to_check,:)>=ub(eq_to_check)-eps);
z(eq_index(eq_to_ignore))=0;
disp_string=' (accounting for MCP tags)';
else
if istag && ~isempty(strmatch('mcp',M_.equations_tags(:,2),'exact'))
disp_string=' (ignoring MCP tags)';
else
disp_string='';
end
end
M_.Sigma_e = Sigma_e;
% Display the non-zero residuals if no return value
if nargout == 0
skipline(4)
skipline(2)
ind = [];
disp('Residuals of the static equations:')
fprintf('Residuals of the static equations%s:',disp_string)
skipline()
any_non_zero_residual = false;
for i=1:M_.orig_eq_nbr

1
tests/Makefile.am
View File

@ -14,6 +14,7 @@ MODFILES = \
lmmcp/purely_backward.mod \
lmmcp/purely_forward.mod \
lmmcp/MCP_ramsey.mod \
lmmcp/SMOPEC.mod \
ep/rbc_mc.mod \
estimation/TaRB/fs2000_tarb.mod \
observation_trends_and_prefiltering/MCMC/Trend_loglin_no_prefilt_first_obs_MC.mod \

73
tests/lmmcp/SMOPEC.mod Normal file
View File

@ -0,0 +1,73 @@
/*
* This file implements a small open economy with impatient agents, i.e. beta<1/(1+r),
* but subject to a borrowing constraint d<2. The problem is set up as a
* mixed complementarity problem (MCP), requiring a setup of the form:
* LB =X => F(X)>0,
* LB<=X<=UB => F(X)=0,
* X =UB => F(X)<0.
* Thus, if d=2, the associated complementary Euler equation needs to return a
* negative residual. Economically, if the borrowing constraint binds, consumption
* today is lower than the unrestricted Euler equation
*
* 1/c=beta*(1+r_bar)*1/c(+1)
*
* would require and marginal utility on the left would be higher. Dynare by
* default computes the residual of an equation LHS=RHS as residual=LHS-RHS. For the
* above equation this would imply
*
* residual=1/c-beta*(1+r_bar)*1/c(+1)
*
* and therefore would return a positive residual, while a negative one is required.
* The equation must therefore be entered as
* beta*(1+r_bar)*1/c(+1)-1/c=0;
* which returns a negative residual when the upper bound is binding.
*
* This implementation was written by Johannes Pfeifer. If you spot any mistakes,
* email me at jpfeifer@gmx.de.
* Please note that the following copyright notice only applies to this Dynare
* implementation of the model.
*/
/*
* Copyright (C) 2022 Johannes Pfeifer
*
* This 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.
*
* It 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.
*
* For a copy of the GNU General Public License,
* see <http://www.gnu.org/licenses/>.
*/
var c d y;
parameters r_bar ${\bar r}$
beta;
r_bar = 0.04; %world interest rate
beta=0.9;
d_bar=2;
model;
[name='Evolution of debt']
y + d= c + (1+r_bar)*d(-1);
[name='Endowment']
y = 1;
[name='Euler equation', mcp = 'd<2']
beta*(1+r_bar)*1/c(+1)-1/c=0;
end;
initval;
d=d_bar;
y=1;
c=y-r_bar*d;
end;
resid;
steady(solve_algo=10);
resid;