1function [residuals,JJacobian] = perfect_foresight_mcp_problem(y, dynamic_function, Y0, YT, ... 2 exo_simul, params, steady_state, ... 3 maximum_lag, T, ny, i_cols, ... 4 i_cols_J1, i_cols_1, i_cols_T, ... 5 i_cols_j, i_cols_0,i_cols_J0, eq_index) 6% function [residuals,JJacobian] = perfect_foresight_mcp_problem(y, dynamic_function, Y0, YT, ... 7% exo_simul, params, steady_state, ... 8% maximum_lag, T, ny, i_cols, ... 9% i_cols_J1, i_cols_1, i_cols_T, ... 10% i_cols_j,eq_index) 11% Computes the residuals and the Jacobian matrix for a perfect foresight problem over T periods 12% in a mixed complementarity problem context 13% 14% INPUTS 15% y [double] N*1 array, terminal conditions for the endogenous variables 16% dynamic_function [handle] function handle to _dynamic-file 17% Y0 [double] N*1 array, initial conditions for the endogenous variables 18% YT [double] N*1 array, terminal conditions for the endogenous variables 19% exo_simul [double] nperiods*M_.exo_nbr matrix of exogenous variables (in declaration order) 20% for all simulation periods 21% params [double] nparams*1 array, parameter values 22% steady_state [double] endo_nbr*1 vector of steady state values 23% maximum_lag [scalar] maximum lag present in the model 24% T [scalar] number of simulation periods 25% ny [scalar] number of endogenous variables 26% i_cols [double] indices of variables appearing in M.lead_lag_incidence 27% and that need to be passed to _dynamic-file 28% i_cols_J1 [double] indices of contemporaneous and forward looking variables 29% appearing in M.lead_lag_incidence 30% i_cols_1 [double] indices of contemporaneous and forward looking variables in 31% M.lead_lag_incidence in dynamic Jacobian (relevant in first period) 32% i_cols_T [double] columns of dynamic Jacobian related to contemporaneous and backward-looking 33% variables (relevant in last period) 34% i_cols_j [double] indices of variables in M.lead_lag_incidence 35% in dynamic Jacobian (relevant in intermediate periods) 36% eq_index [double] N*1 array, index vector describing residual mapping resulting 37% from complementarity setup 38% OUTPUTS 39% residuals [double] (N*T)*1 array, residuals of the stacked problem 40% JJacobian [double] (N*T)*(N*T) array, Jacobian of the stacked problem 41% ALGORITHM 42% None 43% 44% SPECIAL REQUIREMENTS 45% None. 46 47% Copyright (C) 1996-2020 Dynare Team 48% 49% This file is part of Dynare. 50% 51% Dynare is free software: you can redistribute it and/or modify 52% it under the terms of the GNU General Public License as published by 53% the Free Software Foundation, either version 3 of the License, or 54% (at your option) any later version. 55% 56% Dynare is distributed in the hope that it will be useful, 57% but WITHOUT ANY WARRANTY; without even the implied warranty of 58% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 59% GNU General Public License for more details. 60% 61% You should have received a copy of the GNU General Public License 62% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 63 64 65YY = [Y0; y; YT]; 66 67residuals = zeros(T*ny,1); 68if nargout == 2 69 iJacobian = cell(T,1); 70end 71 72i_rows = 1:ny; 73offset = 0; 74i_cols_J = i_cols; 75 76for it = maximum_lag+(1:T) 77 if nargout == 1 78 res = dynamic_function(YY(i_cols),exo_simul, params, ... 79 steady_state,it); 80 residuals(i_rows) = res(eq_index); 81 elseif nargout == 2 82 [res,jacobian] = dynamic_function(YY(i_cols),exo_simul, params, steady_state,it); 83 residuals(i_rows) = res(eq_index); 84 if T==1 && it==maximum_lag+1 85 [rows, cols, vals] = find(jacobian(eq_index,i_cols_0)); 86 if size(jacobian, 1) == 1 % find() will return row vectors in this case 87 rows = rows'; 88 cols = cols'; 89 vals = vals'; 90 end 91 iJacobian{1} = [rows, i_cols_J0(cols), vals]; 92 elseif it == maximum_lag+1 93 [rows,cols,vals] = find(jacobian(eq_index,i_cols_1)); 94 if numel(eq_index) == 1 % find() will return row vectors in this case 95 rows = rows'; 96 cols = cols'; 97 vals = vals'; 98 end 99 iJacobian{1} = [offset+rows, i_cols_J1(cols), vals]; 100 elseif it == maximum_lag+T 101 [rows,cols,vals] = find(jacobian(eq_index,i_cols_T)); 102 if numel(eq_index) == 1 % find() will return row vectors in this case 103 rows = rows'; 104 cols = cols'; 105 vals = vals'; 106 end 107 iJacobian{T} = [offset+rows, i_cols_J(i_cols_T(cols)), vals]; 108 else 109 [rows,cols,vals] = find(jacobian(eq_index,i_cols_j)); 110 if numel(eq_index) == 1 % find() will return row vectors in this case 111 rows = rows'; 112 cols = cols'; 113 vals = vals'; 114 end 115 iJacobian{it-maximum_lag} = [offset+rows, i_cols_J(cols), vals]; 116 i_cols_J = i_cols_J + ny; 117 end 118 offset = offset + ny; 119 end 120 121 i_rows = i_rows + ny; 122 i_cols = i_cols + ny; 123end 124 125if nargout == 2 126 iJacobian = cat(1,iJacobian{:}); 127 JJacobian = sparse(iJacobian(:,1),iJacobian(:,2),iJacobian(:,3),T*ny,T*ny); 128end 129