Added a routine for writing the problem to be solved to compute the BGP of a model.
- Only works with backward models. - Probably doesn't work if the model includes auxiliary variables. - Assumes that the trends are multiplicative.time-shift
Stéphane Adjemia (Scylla)
parent
f7b332efa6
commit
b7c60ddf59
GPG Key ID:
A6D44CB9C64CE77B
|
|
@ -0,0 +1,168 @@
|
|||
function write(DynareModel)
|
||||
|
||||
% Writes the nonlinear problem to be solved for computing the growth
|
||||
% rates and levels along the Balanced Growth Path. Note that for
|
||||
% the variables growing along the BGP, the identified levels are
|
||||
% meaningless, only the relative levels of the growing variables
|
||||
% sharing the same trend(s) are relevant.
|
||||
%
|
||||
% INPUTS
|
||||
% - DynareModel [struct] Dynare generated stucture describing the model (M_).
|
||||
%
|
||||
% OUTPUTS
|
||||
% None
|
||||
%
|
||||
% REMARKS
|
||||
% - The trends are assumed to be multiplicative.
|
||||
|
||||
% Copyright (C) 2019 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/>.
|
||||
|
||||
i0 = transpose(DynareModel.lead_lag_incidence(1,:)); % Indices of the lagged variables.
|
||||
i1 = transpose(DynareModel.lead_lag_incidence(2,:)); % Indices of the current variables.
|
||||
n0 = length(find(i0)); % Number of lagged variables.
|
||||
n1 = length(find(i1)); % Number of current variables.
|
||||
|
||||
% Create function in mod namespace.
|
||||
fid = fopen(sprintf('+%s/bgpfun.m', DynareModel.fname), 'w');
|
||||
|
||||
% Write header.
|
||||
fprintf(fid, 'function [F, JAC] = bgpfun(z)\n\n');
|
||||
fprintf(fid, '%% This file has been generated by dynare (%s).\n\n', datestr(now));
|
||||
|
||||
% The function admits a unique vector as input argument. The first
|
||||
% half of the elements are for the levels of the endogenous
|
||||
% variables, the second half for the growth factors.
|
||||
fprintf(fid, 'y = z(1:%u);\n\n', DynareModel.endo_nbr);
|
||||
fprintf(fid, 'g = z(%u:%u);\n', DynareModel.endo_nbr+1, 2*DynareModel.endo_nbr);
|
||||
|
||||
% Define the point where the dynamic model is to be evaluated.
|
||||
fprintf(fid, 'Y = zeros(%u, 1);\n', 2*(n0+n1));
|
||||
for i=1:length(i0) % period t equations, lagged variables.
|
||||
if i0(i)
|
||||
fprintf(fid, 'Y(%u) = y(%u);\n', i0(i), i);
|
||||
end
|
||||
end
|
||||
for i=1:length(i1) % period t equations, current variables.
|
||||
if i1(i)
|
||||
fprintf(fid, 'Y(%u) = y(%u)*g(%u);\n', i1(i), i, i);
|
||||
end
|
||||
end
|
||||
for i=1:length(i0) % period t+1 equaions lagged variables.
|
||||
if i0(i)
|
||||
fprintf(fid, 'Y(%u) = y(%u)*g(%u);\n', n0+n1+i0(i), i, i);
|
||||
end
|
||||
end
|
||||
for i=1:length(i1) % period t+1 equations current variables.
|
||||
if i1(i)
|
||||
fprintf(fid, 'Y(%u) = y(%u)*g(%u)*g(%u);\n', n0+n1+i1(i), i, i, i);
|
||||
end
|
||||
end
|
||||
fprintf(fid, '\n');
|
||||
|
||||
% Define the vector of parameters.
|
||||
fprintf(fid, 'p = zeros(%u, 1);\n', DynareModel.param_nbr);
|
||||
for i = 1:DynareModel.param_nbr
|
||||
fprintf(fid, 'p(%u) = %16.12f;\n', i, DynareModel.params(i));
|
||||
end
|
||||
fprintf(fid, '\n');
|
||||
|
||||
% Initialize the vector holding the residuals over the two periods.
|
||||
fprintf(fid, 'F = NaN(%u, 1);\n', 2*DynareModel.endo_nbr);
|
||||
|
||||
% Set vector of exogenous variables to 0.
|
||||
fprintf(fid, 'x = zeros(1, %u);\n\n', DynareModel.exo_nbr);
|
||||
|
||||
% Evaluate the residuals and jacobian of the dynamic model in periods t and t+1.
|
||||
fprintf(fid, 'if nargout>1\n');
|
||||
fprintf(fid, ' J = zeros(%u, %u);\n', 2*DynareModel.endo_nbr, n0+n1+DynareModel.endo_nbr);
|
||||
fprintf(fid, ' [F(1:%u), tmp] = %s.dynamic(Y(1:%u), x, p, y, 1);\n', DynareModel.endo_nbr, DynareModel.fname, n1+n0);
|
||||
fprintf(fid, ' J(1:%u,1:%u) = tmp(:,1:%u);\n', DynareModel.endo_nbr, n0+n1, n0+n1);
|
||||
fprintf(fid, ' [F(%u:%u), tmp] = %s.dynamic(Y(1+%u:%u), x, p, y, 1);\n', DynareModel.endo_nbr+1, 2*DynareModel.endo_nbr, DynareModel.fname, n1+n0, 2*(n1+n0));
|
||||
fprintf(fid, ' J(%u:%u,1:%u) = tmp(:,1:%u);\n', DynareModel.endo_nbr+1, 2*DynareModel.endo_nbr, n0+n1, n0+n1);
|
||||
fprintf(fid, 'else\n');
|
||||
fprintf(fid, ' F(1:%u) = %s.dynamic(Y(1:%u), x, p, y, 1);\n', DynareModel.endo_nbr, DynareModel.fname, n1+n0);
|
||||
fprintf(fid, ' F(%u:%u) = %s.dynamic(Y(1+%u:%u), x, p, y, 1);\n', DynareModel.endo_nbr+1, 2*DynareModel.endo_nbr, DynareModel.fname, n1+n0, 2*(n1+n0));
|
||||
fprintf(fid, 'end\n\n');
|
||||
|
||||
% Compute the jacobian if required.
|
||||
fprintf(fid, 'if nargout>1\n');
|
||||
fprintf(fid, ' JAC = zeros(%u,%u);\n', 2*DynareModel.endo_nbr, 2*DynareModel.endo_nbr);
|
||||
|
||||
% Compute the derivatives of the first block of equations (period t)
|
||||
% with respect to the endogenous variables.
|
||||
for i=1:DynareModel.eq_nbr
|
||||
for j=1:DynareModel.endo_nbr
|
||||
if i1(j)
|
||||
if i0(j)
|
||||
fprintf(fid, ' JAC(%u,%u) = J(%u,%u)+J(%u,%u)*g(%u);\n', i, j, i, i0(j), i, i1(j), j);
|
||||
else
|
||||
fprintf(fid, ' JAC(%u,%u) = J(%u,%u)*g(%u);\n', i, j, i, i1(j), j);
|
||||
end
|
||||
else
|
||||
if i0(j)
|
||||
fprintf(fid, ' JAC(%u,%u) = J(%u,%u);\n', i, j, i, i0(j));
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
% Compute the derivatives of the second block of equations (period t+1)
|
||||
% with respect to the endogenous variables.
|
||||
for i=DynareModel.eq_nbr+1:2*DynareModel.eq_nbr
|
||||
for j=1:DynareModel.endo_nbr
|
||||
if i1(j)
|
||||
if i0(j)
|
||||
fprintf(fid, ' JAC(%u,%u) = J(%u,%u)*g(%u)+J(%u,%u)*g(%u)*g(%u);\n', i, j, i, i0(j), j, i, i1(j), j, j);
|
||||
else
|
||||
fprintf(fid, ' JAC(%u,%u) = J(%u,%u)*g(%u)*g(%u);\n', i, j, i, i1(j), j, j);
|
||||
end
|
||||
else
|
||||
if i0(j)
|
||||
fprintf(fid, ' JAC(%u, %u) = J(%u, %u)*g(%u);\n', i, j, i, i0(j), j);
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
% Compute the derivatives of the first block of equations (period t)
|
||||
% with respect to the growth factors.
|
||||
for i=1:DynareModel.eq_nbr
|
||||
for j=1:DynareModel.endo_nbr
|
||||
if i1(j)
|
||||
fprintf(fid, ' J(%u,%u) = J(%u,%u)*y(%u);\n', i, n0+n1+j, i, i1(j), j);
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
% Compute the derivatives of the second block of equations (period t+1)
|
||||
% with respect to the endogenous variables.
|
||||
for i=DynareModel.eq_nbr+1:2*DynareModel.eq_nbr
|
||||
for j=1:DynareModel.endo_nbr
|
||||
if i0(j)
|
||||
fprintf(fid, ' J(%u,%u) = J(%u,%u)+J(%u,%u)*y(%u);\n', i, n0+n1+j, i, n0+n1+j, i, i0(j), j);
|
||||
end
|
||||
if i1(j)
|
||||
fprintf(fid, ' J(%u,%u) = J(%u,%u)+2*J(%u,%u)*y(%u)*g(%u);\n', i, n0+n1+j, i, n0+n1+j, i, i1(j), j, j);
|
||||
end
|
||||
end
|
||||
end
|
||||
fprintf(fid, ' JAC(:,%u:%u) = J(:,%u:%u);\n', DynareModel.endo_nbr+1, 2*DynareModel.endo_nbr, n0+n1+1, n0+n1+DynareModel.endo_nbr);
|
||||
fprintf(fid,'end\n');
|
||||
|
||||
% Close file.
|
||||
fclose(fid);
|
||||
|
|
@ -357,7 +357,8 @@ MODFILES = \
|
|||
observation_trends_and_prefiltering/calib_smoother/Tr_prefilter_loglin_calib_smoother.mod \
|
||||
observation_trends_and_prefiltering/calib_smoother/Tr_prefil_f_obs_loglin_cal_smoother.mod \
|
||||
observation_trends_and_prefiltering/ML/Trend_no_prefilter_selected_var.mod \
|
||||
dynare-command-options/ramst.mod
|
||||
bgp/solow-1/solow.mod \
|
||||
dynare-command-options/ramst.mod
|
||||
|
||||
PARTICLEFILES = \
|
||||
particle/dsge_base2.mod \
|
||||
|
|
|
|||
|
|
@ -0,0 +1,80 @@
|
|||
// This mod file test the bgp.write function by characterizing the Balanced Growth Path of the Solow model. This is done 10000 times in a Monte Carlo loop
|
||||
// randomizing the initial guess of the nonlinear solver.
|
||||
|
||||
var Efficiency
|
||||
EfficiencyGrowth
|
||||
Population
|
||||
PopulationGrowth
|
||||
Output
|
||||
PhysicalCapitalStock ;
|
||||
|
||||
varexo e_x
|
||||
e_n ;
|
||||
|
||||
parameters alpha
|
||||
delta
|
||||
s
|
||||
rho_x
|
||||
rho_n
|
||||
EfficiencyGrowth_ss
|
||||
PopulationGrowth_ss ;
|
||||
|
||||
alpha = .33;
|
||||
delta = .02;
|
||||
s = .20;
|
||||
rho_x = .90;
|
||||
rho_n = .95;
|
||||
EfficiencyGrowth_ss = 1.01;
|
||||
PopulationGrowth_ss = 1.01;
|
||||
|
||||
model;
|
||||
Efficiency = EfficiencyGrowth*Efficiency(-1);
|
||||
EfficiencyGrowth/EfficiencyGrowth_ss = (EfficiencyGrowth(-1)/EfficiencyGrowth_ss)^(rho_x)*exp(e_x);
|
||||
Population = PopulationGrowth*Population(-1);
|
||||
PopulationGrowth/PopulationGrowth_ss = (PopulationGrowth(-1)/PopulationGrowth_ss)^(rho_n)*exp(e_n);
|
||||
Output = PhysicalCapitalStock(-1)^alpha*(Efficiency*Population)^(1-alpha);
|
||||
PhysicalCapitalStock = (1-delta)*PhysicalCapitalStock(-1) + s*Output;
|
||||
end;
|
||||
|
||||
|
||||
verbatim;
|
||||
|
||||
bgp.write(M_);
|
||||
MC = 10000;
|
||||
KY = NaN(MC,1);
|
||||
GY = NaN(MC,1);
|
||||
GK = NaN(MC,1);
|
||||
EG = NaN(MC,1);
|
||||
options = optimoptions('fsolve','Display','off','Algorithm','levenberg-marquardt','MaxFunctionEvaluations',1000000,'MaxIterations',100000,'SpecifyObjectiveGradient',true,'FunctionTolerance',1e-8,'StepTolerance',1e-8);
|
||||
for i=1:MC
|
||||
y = 1+(rand(6,1)-.5)*.2;
|
||||
g = ones(6,1);
|
||||
[y, fval, exitflag] = fsolve(@solow.bgpfun, [y;g], options);
|
||||
if exitflag>0
|
||||
KY(i) = y(6)/y(5);
|
||||
GY(i) = y(11);
|
||||
GK(i) = y(12);
|
||||
EG(i) = y(2);
|
||||
end
|
||||
end
|
||||
% Compute the physical capital stock over output ratio along the BGP as
|
||||
% a function of the deep parameters...
|
||||
theoretical_long_run_ratio = s*EfficiencyGrowth_ss*PopulationGrowth_ss/(EfficiencyGrowth_ss*PopulationGrowth_ss-1+delta);
|
||||
% ... and compare to the average ratio in the Monte Carlo
|
||||
mu = mean(KY(~isnan(KY)));
|
||||
s2 = var(KY(~isnan(KY)));
|
||||
nn = sum(isnan(KY))/MC;
|
||||
disp(sprintf('Theoretical K/Y BGP ratio = %s.', num2str(theoretical_long_run_ratio)));
|
||||
disp(sprintf('mean(K/Y) = %s.', num2str(mu)));
|
||||
disp(sprintf('var(K/Y) = %s.', num2str(s2)));
|
||||
disp(sprintf('Number of failures: %u (over %u problems).', sum(isnan(KY)), MC));
|
||||
assert(abs(mu-theoretical_long_run_ratio)<1e-8);
|
||||
assert(s2<1e-16);
|
||||
assert(abs(mean(GY(~isnan(GY)))-EfficiencyGrowth_ss*PopulationGrowth_ss)<1e-8)
|
||||
assert(var(GY(~isnan(GY)))<1e-16);
|
||||
assert(abs(mean(GK(~isnan(GK)))-EfficiencyGrowth_ss*PopulationGrowth_ss)<1e-8)
|
||||
assert(var(GK(~isnan(GK)))<1e-16);
|
||||
assert(abs(mean(EG(~isnan(EG)))-EfficiencyGrowth_ss)<1e-8)
|
||||
assert(var(EG(~isnan(EG)))<1e-16);
|
||||
|
||||
end;
|
||||
Loading…
Reference in New Issue