diff --git a/matlab/+bgp/write.m b/matlab/+bgp/write.m new file mode 100644 index 000000000..35909154f --- /dev/null +++ b/matlab/+bgp/write.m @@ -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 . + +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); \ No newline at end of file diff --git a/tests/Makefile.am b/tests/Makefile.am index ca4c30634..fceb5fc65 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -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 \ diff --git a/tests/bgp/solow-1/solow.mod b/tests/bgp/solow-1/solow.mod new file mode 100644 index 000000000..74e979296 --- /dev/null +++ b/tests/bgp/solow-1/solow.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; \ No newline at end of file