diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m
index 06e2b266dc3d2beeebe13deb28050bb6a39f6443..17934f30ee560ff7966432599c6fec6d943b40bd 100644
--- a/matlab/dynare_solve.m
+++ b/matlab/dynare_solve.m
@@ -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')
diff --git a/matlab/evaluate_steady_state.m b/matlab/evaluate_steady_state.m
index 0dfd84d95f7bf4129c3944ee1a38aad64363bed2..fb1a3c4a904e958a6a69e155d9f7d6741908aec8 100644
--- a/matlab/evaluate_steady_state.m
+++ b/matlab/evaluate_steady_state.m
@@ -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);
diff --git a/matlab/resid.m b/matlab/resid.m
index 2c6d544edbd38c628a3de847488db3cebd02f40a..b7586ccea1490c72468cc91716cdd9e1edbef726 100644
--- a/matlab/resid.m
+++ b/matlab/resid.m
@@ -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
diff --git a/tests/Makefile.am b/tests/Makefile.am
index d537e3bdd91966f341de38333ddd9ad51ecf883e..e162bc767511d87c24d24f6aceade63ebf7e9761 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -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 \
diff --git a/tests/lmmcp/SMOPEC.mod b/tests/lmmcp/SMOPEC.mod
new file mode 100644
index 0000000000000000000000000000000000000000..9aa97bf9d7fe4b7d9f31cffaa2155a60f8eabac5
--- /dev/null
+++ b/tests/lmmcp/SMOPEC.mod
@@ -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;
\ No newline at end of file