lmmcp: display norm of residuals in homotopy

matlab/perfect-foresight-models/perfect_foresight_solver_core.m

@@ -114,7 +114,7 @@ else
                     [oo_.endo_simul, oo_.deterministic_simulation] = ...
                         solve_stacked_linear_problem(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, M_, options_);
                 else
-                    [oo_.endo_simul, oo_.deterministic_simulation] = ...
+                    [oo_.endo_simul, oo_.deterministic_simulation, residuals] = ...
                         solve_stacked_problem(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
                 end
               otherwise
@@ -126,7 +126,7 @@ end
 
 if nargout>1
     if options_.lmmcp.status
-        maxerror = NaN; % Could be improved
+        maxerror = max(max(abs(residuals)));
     elseif options_.block && ~options_.bytecode
         maxerror = oo_.deterministic_simulation.error;
     else

matlab/perfect-foresight-models/solve_stacked_problem.m

@@ -1,4 +1,4 @@
-function [endogenousvariables, info] = solve_stacked_problem(endogenousvariables, exogenousvariables, steadystate, M, options)
+function [endogenousvariables, info, residuals] = solve_stacked_problem(endogenousvariables, exogenousvariables, steadystate, M, options)
 
 % Solves the perfect foresight model using dynare_solve
 %
@@ -12,6 +12,7 @@ function [endogenousvariables, info] = solve_stacked_problem(endogenousvariables
 % OUTPUTS
 % - endogenousvariables [double] N*T array, paths for the endogenous variables (solution of the perfect foresight model).
 % - info                [struct] contains informations about the results.
+% - residuals           [double] N*T array, residuals of the equations (with 0 for initial condition)
 
 % Copyright © 2015-2022 Dynare Team
 %
@@ -46,7 +47,7 @@ if (options.solve_algo == 10 || options.solve_algo == 11)% mixed complementarity
         options.mcppath.lb = repmat(lb,options.periods,1);
         options.mcppath.ub = repmat(ub,options.periods,1);
     end
-    [y, check, ~, ~, errorcode] = dynare_solve(@perfect_foresight_mcp_problem, z(:), ...
+    [y, check, res, ~, errorcode] = dynare_solve(@perfect_foresight_mcp_problem, z(:), ...
                                                options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
                                                options, ...
                                                dynamicmodel, y0, yT, ...
@@ -55,7 +56,7 @@ if (options.solve_algo == 10 || options.solve_algo == 11)% mixed complementarity
                                                i_cols_J1, i_cols_1, i_cols_T, i_cols_j, i_cols_0, i_cols_J0, ...
                                                eq_index);
 else
-    [y, check, ~, ~, errorcode] = dynare_solve(@perfect_foresight_problem, z(:), ...
+    [y, check, res, ~, errorcode] = dynare_solve(@perfect_foresight_problem, z(:), ...
                                                options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
                                                options, y0, yT, exogenousvariables, M.params, steadystate, options.periods, M, options);
 end
@@ -69,6 +70,8 @@ else
 end
 
 endogenousvariables(:, M.maximum_lag+(1:options.periods)) = reshape(y, M.endo_nbr, options.periods);
+residuals=zeros(size(endogenousvariables));
+residuals(:, M.maximum_lag+(1:options.periods)) = reshape(res, M.endo_nbr, options.periods);
 
 if check
     info.status = false;