diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m index 792fe213c..0e61f07b9 100644 --- a/matlab/dynare_solve.m +++ b/matlab/dynare_solve.m @@ -205,6 +205,7 @@ elseif options.solve_algo == 11 global mcp_data mcp_data.func = func; mcp_data.args = varargin; + [x,fval,jac,mu,status] = pathmcp(x,omcppath.lb,omcppath.ub,'mcp_func',omcppath.A,omcppath.b,omcppath.t,omcppath.mu0); info = ~status; else diff --git a/matlab/perfect-foresight-models/perfect_foresight_mcp_problem.m b/matlab/perfect-foresight-models/perfect_foresight_mcp_problem.m new file mode 100644 index 000000000..22156a1b7 --- /dev/null +++ b/matlab/perfect-foresight-models/perfect_foresight_mcp_problem.m @@ -0,0 +1,82 @@ +function [residuals,JJacobian] = perfect_foresight_problem(y, dynamic_function, Y0, YT, ... + exo_simul, params, steady_state, ... + maximum_lag, T, ny, i_cols, ... + i_cols_J1, i_cols_1, i_cols_T, ... + i_cols_j,nnzJ,eq_index) +% function [residuals,JJacobian] = perfect_foresight_problem(x, model_dynamic, Y0, YT,exo_simul, +% params, steady_state, maximum_lag, periods, ny, i_cols,i_cols_J1, i_cols_1, +% i_cols_T, i_cols_j, nnzA) +% computes the residuals and th Jacobian matrix +% for a perfect foresight problem over T periods. +% +% INPUTS +% ... +% OUTPUTS +% ... +% ALGORITHM +% ... +% +% SPECIAL REQUIREMENTS +% None. + +% Copyright (C) 1996-2015 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 . + + + YY = [Y0; y; YT]; + + residuals = zeros(T*ny,1); + if nargout == 2 + iJacobian = cell(T,1); + end + + i_rows = 1:ny; + offset = 0; + i_cols_J = i_cols; + + for it = 2:(T+1) + if nargout == 1 + res = dynamic_function(YY(i_cols),exo_simul, params, ... + steady_state,it); + residuals(i_rows) = res(eq_index); + elseif nargout == 2 + [res,jacobian] = dynamic_function(YY(i_cols),exo_simul, params, ... + steady_state,it); + residuals(i_rows) = res(eq_index); + if it == 2 + [rows,cols,vals] = find(jacobian(eq_index,i_cols_1)); + iJacobian{1} = [offset+rows, i_cols_J1(cols), vals]; + elseif it == T + 1 + [rows,cols,vals] = find(jacobian(eq_index,i_cols_T)); + iJacobian{T} = [offset+rows, i_cols_J(i_cols_T(cols)), vals]; + else + [rows,cols,vals] = find(jacobian(eq_index,i_cols_j)); + iJacobian{it-1} = [offset+rows, i_cols_J(cols), vals]; + i_cols_J = i_cols_J + ny; + end + offset = offset + ny; + end + + i_rows = i_rows + ny; + i_cols = i_cols + ny; + end + + if nargout == 2 + iJacobian = cat(1,iJacobian{:}); + JJacobian = sparse(iJacobian(:,1),iJacobian(:,2),iJacobian(:,3),T* ... + ny,T*ny); + end \ No newline at end of file diff --git a/matlab/perfect-foresight-models/perfect_foresight_problem.m b/matlab/perfect-foresight-models/perfect_foresight_problem.m index 9030a5ae5..adce18df5 100644 --- a/matlab/perfect-foresight-models/perfect_foresight_problem.m +++ b/matlab/perfect-foresight-models/perfect_foresight_problem.m @@ -41,31 +41,40 @@ function [residuals,JJacobian] = perfect_foresight_problem(y, dynamic_function, residuals = zeros(T*ny,1); if nargout == 2 - JJacobian = sparse([],[],[],T*ny,T*ny,T*nnzJ); + iJacobian = cell(T,1); end i_rows = 1:ny; i_cols_J = i_cols; - - for it = maximum_lag+(1:T) + offset = 0; + + for it = 2:(T+1) if nargout == 1 - residuals(i_rows) = dynamic_function(YY(i_cols),exo_simul, params, ... + residuals(i_rows) = dynamic_function(YY(i_cols),exo_simul, params, ... steady_state,it); elseif nargout == 2 [residuals(i_rows),jacobian] = dynamic_function(YY(i_cols),exo_simul, params, ... steady_state,it); - if it == 2 - JJacobian(i_rows,i_cols_J1) = jacobian(:,i_cols_1); + [rows,cols,vals] = find(jacobian(:,i_cols_1)); + iJacobian{1} = [offset+rows, i_cols_J1(cols), vals]; elseif it == T + 1 - JJacobian(i_rows,i_cols_J(i_cols_T)) = jacobian(:,i_cols_T); + [rows,cols,vals] = find(jacobian(:,i_cols_T)); + iJacobian{T} = [offset+rows, i_cols_J(i_cols_T(cols)), vals]; else - JJacobian(i_rows,i_cols_J) = jacobian(:,i_cols_j); + [rows,cols,vals] = find(jacobian(:,i_cols_j)); + iJacobian{it-1} = [offset+rows, i_cols_J(cols), vals]; i_cols_J = i_cols_J + ny; end + offset = offset + ny; end i_rows = i_rows + ny; i_cols = i_cols + ny; end + if nargout == 2 + iJacobian = cat(1,iJacobian{:}); + JJacobian = sparse(iJacobian(:,1),iJacobian(:,2),iJacobian(:,3),T* ... + ny,T*ny); + end \ No newline at end of file diff --git a/matlab/perfect-foresight-models/solve_stacked_problem.m b/matlab/perfect-foresight-models/solve_stacked_problem.m index 3211fbfb1..7b1ab780e 100644 --- a/matlab/perfect-foresight-models/solve_stacked_problem.m +++ b/matlab/perfect-foresight-models/solve_stacked_problem.m @@ -20,12 +20,33 @@ function [endogenousvariables, info] = solve_stacked_problem(endogenousvariables [options, y0, yT, z, i_cols, i_cols_J1, i_cols_T, i_cols_j, i_cols_1, dynamicmodel] = ... initialize_stacked_problem(endogenousvariables, options, M, steadystate); -[y, check] = dynare_solve(@perfect_foresight_problem,z(:),options, ... - dynamicmodel, y0, yT, ... - exogenousvariables, M.params, steadystate, ... - M.maximum_lag, options.periods, M.endo_nbr, i_cols, ... - i_cols_J1, i_cols_1, i_cols_T, i_cols_j, ... - M.NNZDerivatives(1)); +if (options.solve_algo == 10 || options.solve_algo == 11) + [lb,ub,eq_index] = get_complementarity_conditions(M,options.ramsey_policy); + if options.linear_approximation + lb = lb - steadystate_y; + ub = ub - steadystate_y; + end + if options.solve_algo == 10 + options.lmmcp.lb = repmat(lb,options.periods,1); + options.lmmcp.ub = repmat(ub,options.periods,1); + elseif options.solve_algo == 11 + options.mcppath.lb = repmat(lb,options.periods,1); + options.mcppath.ub = repmat(ub,options.periods,1); + end + [y, check] = dynare_solve(@perfect_foresight_mcp_problem,z(:),options, ... + dynamicmodel, y0, yT, ... + exogenousvariables, M.params, steadystate, ... + M.maximum_lag, options.periods, M.endo_nbr, i_cols, ... + i_cols_J1, i_cols_1, i_cols_T, i_cols_j, ... + M.NNZDerivatives(1),eq_index); +else + [y, check] = dynare_solve(@perfect_foresight_problem,z(:),options, ... + dynamicmodel, y0, yT, ... + exogenousvariables, M.params, steadystate, ... + M.maximum_lag, options.periods, M.endo_nbr, i_cols, ... + i_cols_J1, i_cols_1, i_cols_T, i_cols_j, ... + M.NNZDerivatives(1)); +end if all(imag(y)<.1*options.dynatol.x) if ~isreal(y) @@ -40,5 +61,6 @@ endogenousvariables = [y0 reshape(y, M.endo_nbr, options.periods) yT]; if check info.status = false; else + info.status = true; end \ No newline at end of file diff --git a/tests/Makefile.am b/tests/Makefile.am index 485d18946..458d8cd3b 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -242,7 +242,8 @@ MODFILES = \ deterministic_simulations/multiple_lead_lags/AR2.mod \ deterministic_simulations/multiple_lead_lags/AR2_forward.mod \ deterministic_simulations/multiple_lead_lags/ramst_augmented_histval.mod \ - deterministic_simulations/lbj/rbc.mod \ + deterministic_simulations/rbc_det.mod \ + deterministic_simulations/rbc_det_stack_solve_algo_7.mod \ lmmcp/rbc.mod \ lmmcp/rbcii.mod \ lmmcp/sw_lmmcp.mod \ @@ -402,6 +403,8 @@ deterministic_simulations/rbc_det_exo_lag_2b.m.trs: deterministic_simulations/rb deterministic_simulations/rbc_det_exo_lag_2c.m.trs: deterministic_simulations/rbc_det_exo_lag_2a.m.trs deterministic_simulations/rbc_det_exo_lag_2b.o.trs: deterministic_simulations/rbc_det_exo_lag_2a.o.trs deterministic_simulations/rbc_det_exo_lag_2c.o.trs: deterministic_simulations/rbc_det_exo_lag_2a.o.trs +deterministic_simulations/rbc_det_stack_solve_algo_7.m.trs: deterministic_simulations/rbc_det.m.trs +deterministic_simulations/rbc_det_stack_solve_algo_7.o.trs: deterministic_simulations/rbc_det.o.trs initval_file/ramst_initval_file.m.trs: initval_file/ramst_initval_file_data.m.tls initval_file/ramst_initval_file.o.trs: initval_file/ramst_initval_file_data.o.tls diff --git a/tests/deterministic_simulations/rbc_det.mod b/tests/deterministic_simulations/rbc_det.mod new file mode 100644 index 000000000..77388ac6f --- /dev/null +++ b/tests/deterministic_simulations/rbc_det.mod @@ -0,0 +1,78 @@ +var Capital, Output, Labour, Consumption, Efficiency, efficiency, ExpectedTerm; + +varexo EfficiencyInnovation; + +parameters beta, theta, tau, alpha, psi, delta, rho, effstar, sigma2; + +beta = 0.9900; +theta = 0.3570; +tau = 2.0000; +alpha = 0.4500; +psi = -0.1000; +delta = 0.0200; +rho = 0.8000; +effstar = 1.0000; +sigma2 = 0; + +model; + + // Eq. n°1: + efficiency = rho*efficiency(-1) + EfficiencyInnovation; + + // Eq. n°2: + Efficiency = effstar*exp(efficiency); + + // Eq. n°3: + Output = Efficiency*(alpha*(Capital(-1)^psi)+(1-alpha)*(Labour^psi))^(1/psi); + + // Eq. n°4: + Capital = Output-Consumption + (1-delta)*Capital(-1); + + // Eq. n°5: + ((1-theta)/theta)*(Consumption/(1-Labour)) - (1-alpha)*(Output/Labour)^(1-psi); + + // Eq. n°6: + (((Consumption^theta)*((1-Labour)^(1-theta)))^(1-tau))/Consumption = ExpectedTerm(1); + + // Eq. n°7: + ExpectedTerm = beta*((((Consumption^theta)*((1-Labour)^(1-theta)))^(1-tau))/Consumption)*(alpha*((Output/Capital(-1))^(1-psi))+(1-delta)); + +end; + +steady_state_model; +efficiency = EfficiencyInnovation/(1-rho); +Efficiency = effstar*exp(efficiency); +Output_per_unit_of_Capital=((1/beta-1+delta)/alpha)^(1/(1-psi)); +Consumption_per_unit_of_Capital=Output_per_unit_of_Capital-delta; +Labour_per_unit_of_Capital=(((Output_per_unit_of_Capital/Efficiency)^psi-alpha)/(1-alpha))^(1/psi); +Output_per_unit_of_Labour=Output_per_unit_of_Capital/Labour_per_unit_of_Capital; +Consumption_per_unit_of_Labour=Consumption_per_unit_of_Capital/Labour_per_unit_of_Capital; + +% Compute steady state share of capital. +ShareOfCapital=alpha/(alpha+(1-alpha)*Labour_per_unit_of_Capital^psi); + +% Compute steady state of the endogenous variables. +Labour=1/(1+Consumption_per_unit_of_Labour/((1-alpha)*theta/(1-theta)*Output_per_unit_of_Labour^(1-psi))); +Consumption=Consumption_per_unit_of_Labour*Labour; +Capital=Labour/Labour_per_unit_of_Capital; +Output=Output_per_unit_of_Capital*Capital; +ExpectedTerm=beta*((((Consumption^theta)*((1-Labour)^(1-theta)))^(1-tau))/Consumption) + *(alpha*((Output/Capital)^(1-psi))+1-delta); +end; + +steady; + +ik = varlist_indices('Capital',M_.endo_names); +CapitalSS = oo_.steady_state(ik); + +histval; +Capital(0) = CapitalSS/2; +end; + + +perfect_foresight_setup(periods=200); + +perfect_foresight_solver; + +rplot Consumption; +rplot Capital; \ No newline at end of file diff --git a/tests/deterministic_simulations/rbc_det_stack_solve_algo_7.mod b/tests/deterministic_simulations/rbc_det_stack_solve_algo_7.mod new file mode 100644 index 000000000..094da7933 --- /dev/null +++ b/tests/deterministic_simulations/rbc_det_stack_solve_algo_7.mod @@ -0,0 +1,85 @@ +var Capital, Output, Labour, Consumption, Efficiency, efficiency, ExpectedTerm; + +varexo EfficiencyInnovation; + +parameters beta, theta, tau, alpha, psi, delta, rho, effstar, sigma2; + +beta = 0.9900; +theta = 0.3570; +tau = 2.0000; +alpha = 0.4500; +psi = -0.1000; +delta = 0.0200; +rho = 0.8000; +effstar = 1.0000; +sigma2 = 0; + +model; + + // Eq. n°1: + efficiency = rho*efficiency(-1) + EfficiencyInnovation; + + // Eq. n°2: + Efficiency = effstar*exp(efficiency); + + // Eq. n°3: + Output = Efficiency*(alpha*(Capital(-1)^psi)+(1-alpha)*(Labour^psi))^(1/psi); + + // Eq. n°4: + Capital = Output-Consumption + (1-delta)*Capital(-1); + + // Eq. n°5: + ((1-theta)/theta)*(Consumption/(1-Labour)) - (1-alpha)*(Output/Labour)^(1-psi); + + // Eq. n°6: + (((Consumption^theta)*((1-Labour)^(1-theta)))^(1-tau))/Consumption = ExpectedTerm(1); + + // Eq. n°7: + ExpectedTerm = beta*((((Consumption^theta)*((1-Labour)^(1-theta)))^(1-tau))/Consumption)*(alpha*((Output/Capital(-1))^(1-psi))+(1-delta)); + +end; + +steady_state_model; +efficiency = EfficiencyInnovation/(1-rho); +Efficiency = effstar*exp(efficiency); +Output_per_unit_of_Capital=((1/beta-1+delta)/alpha)^(1/(1-psi)); +Consumption_per_unit_of_Capital=Output_per_unit_of_Capital-delta; +Labour_per_unit_of_Capital=(((Output_per_unit_of_Capital/Efficiency)^psi-alpha)/(1-alpha))^(1/psi); +Output_per_unit_of_Labour=Output_per_unit_of_Capital/Labour_per_unit_of_Capital; +Consumption_per_unit_of_Labour=Consumption_per_unit_of_Capital/Labour_per_unit_of_Capital; + +% Compute steady state share of capital. +ShareOfCapital=alpha/(alpha+(1-alpha)*Labour_per_unit_of_Capital^psi); + +% Compute steady state of the endogenous variables. +Labour=1/(1+Consumption_per_unit_of_Labour/((1-alpha)*theta/(1-theta)*Output_per_unit_of_Labour^(1-psi))); +Consumption=Consumption_per_unit_of_Labour*Labour; +Capital=Labour/Labour_per_unit_of_Capital; +Output=Output_per_unit_of_Capital*Capital; +ExpectedTerm=beta*((((Consumption^theta)*((1-Labour)^(1-theta)))^(1-tau))/Consumption) + *(alpha*((Output/Capital)^(1-psi))+1-delta); +end; + +steady; + +ik = varlist_indices('Capital',M_.endo_names); +CapitalSS = oo_.steady_state(ik); + +histval; +Capital(0) = CapitalSS/2; +end; + + +perfect_foresight_setup(periods=200); + +perfect_foresight_solver(stack_solve_algo=7,solve_algo=1); + +rplot Consumption; +rplot Capital; + +D = load('rbc_det_results'); + +if norm(D.oo_.endo_simul - oo_.endo_simul) > 1e-30; + disp(norm(D.oo_.endo_simul - oo_.endo_simul)); + error('rbc_det_stack_solve_algo_7 failed'); +end; \ No newline at end of file