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

fixing bug in mcp model setup 

makes perfect_foresight_problem more efficient
added test case for stack_solve_algo == 7
time-shift
Michel Juillard 2016-07-02 22:04:21 +02:00
parent 2f2413e64a
commit 22f49971bc
7 changed files with 295 additions and 15 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
matlab/dynare_solve.m
View File

@ -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

82
matlab/perfect-foresight-models/perfect_foresight_mcp_problem.m Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
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

25
matlab/perfect-foresight-models/perfect_foresight_problem.m
View File

@ -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

34
matlab/perfect-foresight-models/solve_stacked_problem.m
View File

@ -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

5
tests/Makefile.am
View File

@ -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

78
tests/deterministic_simulations/rbc_det.mod Normal file
View File

@ -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;

85
tests/deterministic_simulations/rbc_det_stack_solve_algo_7.mod Normal file
View File

@ -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;