From 0c08b0fdf42a64ffedd49e2536f0e63a2dccf098 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Fr=C3=A9d=C3=A9ric=20Karam=C3=A9?=
<frederic.karame@univ-lemans.fr>
Date: Thu, 13 Mar 2025 11:54:41 +0100
Subject: [PATCH] Add the One Step estimation method as mode_compute=99.
---
matlab/default_option_values.m | 10 +-
matlab/estimation/dynare_estimation_1.m | 57 ++++--
.../optimization/dynare_minimize_objective.m | 21 +++
tests/optimizers/fs2000_99.mod | 170 ++++++++++++++++++
4 files changed, 241 insertions(+), 17 deletions(-)
create mode 100644 tests/optimizers/fs2000_99.mod
diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index 929079add..a655489ad 100644
--- a/matlab/default_option_values.m
+++ b/matlab/default_option_values.m
@@ -662,10 +662,12 @@ particleswarm.UseVectorized = false;
options_.particleswarm = particleswarm;
% One step
-one_step.gradient = [] ;
-one_step.hessian = [] ;
-one_step.outer_product = 0 ;
-options_.one_step= one_step ;
+one_step.prelim = [];
+one_step.delta = 0.8;
+one_step.gradient = [];
+one_step.hessien = [];
+one_step.first_mode_compute = 3;
+options_.one_step = one_step;
% prior analysis
options_.prior_mc = 20000;
diff --git a/matlab/estimation/dynare_estimation_1.m b/matlab/estimation/dynare_estimation_1.m
index e82e4389d..2856abd50 100644
--- a/matlab/estimation/dynare_estimation_1.m
+++ b/matlab/estimation/dynare_estimation_1.m
@@ -123,6 +123,28 @@ end
[dataset_, dataset_info, xparam1, hh, M_, options_, oo_, estim_params_, bayestopt_, bounds] = ...
dynare_estimation_init(var_list_, dname, [], M_, options_, oo_, estim_params_, bayestopt_);
+if options_.mode_compute==99
+ if ~isempty(options_.optim_opt)
+ options_list = read_key_value_string(options_.optim_opt);
+ for i=1:rows(options_list)
+ switch options_list{i,1}
+ case 'delta'
+ options_.one_step.delta = options_list{i,2};
+ case 'first_mode_compute'
+ options_.one_step.first_mode_compute = options_list{i,2};
+ otherwise
+ warning(['OneStep: Unknown option (' options_list{i,1} ')!'])
+ end
+ end
+ options_.optim_opt = [];
+ end
+ old_dataset_ = dataset_;
+ old_dataset_info = dataset_info;
+ Tprelim = ceil(power(options_.nobs,options_.one_step.delta));
+ dataset_info.rawdata = dataset_info.rawdata(1:Tprelim,:);
+ dataset_ = dataset_(dataset_.dates(options_.first_obs:options_.first_obs+Tprelim-1));
+end
+
if options_.dsge_var
check_dsge_var_model(M_, estim_params_, bayestopt_);
if dataset_info.missing.state
@@ -234,9 +256,17 @@ end
%
if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation && ~issmc(options_)
- optimizer_vec = [options_.mode_compute;num2cell(options_.additional_optimizer_steps)];
+ if options_.mode_compute==99
+ optimizer_vec = [num2cell(options_.one_step.first_mode_compute);options_.mode_compute];
+ else
+ optimizer_vec = [options_.mode_compute;num2cell(options_.additional_optimizer_steps)];
+ end
for optim_iter = 1:length(optimizer_vec)
current_optimizer = optimizer_vec{optim_iter};
+ if current_optimizer==99
+ dataset_ = old_dataset_;
+ dataset_info.rawdata = old_dataset_info.rawdata;
+ end
if isnumeric(current_optimizer)
if current_optimizer==5
if options_.analytic_derivation
@@ -578,18 +608,19 @@ if options_.particle.status
end
%Run and store classical smoother if needed
-if options_.smoother && ... %Bayesian smoother requested before
- (any(bayestopt_.pshape > 0) && options_.mh_replic || ... % Bayesian with MCMC run
- any(bayestopt_.pshape > 0) && options_.load_mh_file) % Bayesian with loaded MCMC
- % nothing to do
-elseif options_.partial_information ||...
- options_.order>1 %no particle smoother
- % smoothing not yet supported
-else
- %% ML estimation, or posterior mode without Metropolis-Hastings or Metropolis without Bayesian smoothed variables
- oo_=save_display_classical_smoother_results(xparam1,M_,oo_,options_,bayestopt_,dataset_,dataset_info,estim_params_);
-end
-
+if options_.mode_compute ~= 99
+ if options_.smoother && ... %Bayesian smoother requested before
+ (any(bayestopt_.pshape > 0) && options_.mh_replic || ... % Bayesian with MCMC run
+ any(bayestopt_.pshape > 0) && options_.load_mh_file) % Bayesian with loaded MCMC
+ % nothing to do
+ elseif options_.partial_information ||...
+ options_.order>1 %no particle smoother
+ % smoothing not yet supported
+ else
+ %% ML estimation, or posterior mode without Metropolis-Hastings or Metropolis without Bayesian smoothed variables
+ oo_=save_display_classical_smoother_results(xparam1,M_,oo_,options_,bayestopt_,dataset_,dataset_info,estim_params_);
+ end
+end
if options_.forecast == 0 || options_.mh_replic > 0 || options_.load_mh_file
% nothing to do
elseif options_.order>1 && M_.exo_det_nbr == 0 || ...
diff --git a/matlab/optimization/dynare_minimize_objective.m b/matlab/optimization/dynare_minimize_objective.m
index a1cc385e3..a59fe9b3a 100644
--- a/matlab/optimization/dynare_minimize_objective.m
+++ b/matlab/optimization/dynare_minimize_objective.m
@@ -496,6 +496,27 @@ switch minimizer_algorithm
% subvarargin = [varargin(1), varargin(3:6), varargin(8)];
% opt_par_values = online_auxiliary_filter(start_par_value, subvarargin{:});
warning('Online particle filter is no more available with mode_compute=11; use posterior_sampler with option online')
+ case 99
+ options_.one_step.prelim = start_par_value;
+ optim_options = optimoptions('fminunc','display','off','MaxIterations',0,'TolFun',1e9,'TolX',1e-6);
+ [~,fval,exitflag,~,options_.one_step.gradient,options_.one_step.hessian] = fminunc(func2str(objective_function),start_par_value,optim_options,varargin{:});
+ opt_par_values = start_par_value - options_.one_step.hessian\options_.one_step.gradient;
+ outside_bound_pars=find(opt_par_values < varargin{7}.lb | opt_par_values > varargin{7}.ub);
+ if ~isempty(outside_bound_pars)
+ for ii=1:length(outside_bound_pars)
+ outside_bound_par_names{ii,1}=get_the_name(outside_bound_pars(ii),0,varargin{4},varargin{5},options_.varobs);
+ end
+ disp_string=[outside_bound_par_names{1,:}];
+ for ii=2:size(outside_bound_par_names,1)
+ disp_string=[disp_string,', ',outside_bound_par_names{ii,:}];
+ end
+ fprintf('\n')
+ if size(outside_bound_pars,1)==1
+ disp(['WARNING: ', disp_string ,' is outside parameter bounds. You should correct the prior or consider to increase delta.'])
+ elseif size(outside_bound_pars,1)>1
+ disp(['WARNING: ', disp_string ,' are outside parameter bounds. You should correct the prior or consider to increase delta.'])
+ end
+ end
case 12
if isoctave
error('Option mode_compute=12 is not available under Octave')
diff --git a/tests/optimizers/fs2000_99.mod b/tests/optimizers/fs2000_99.mod
new file mode 100644
index 000000000..cc28dbc63
--- /dev/null
+++ b/tests/optimizers/fs2000_99.mod
@@ -0,0 +1,170 @@
+/*
+ * This file replicates the estimation of the cash in advance model (termed M1
+ * in the paper) described in Frank Schorfheide (2000): "Loss function-based
+ * evaluation of DSGE models", Journal of Applied Econometrics, 15(6), 645-670.
+ *
+ * The data are taken from the replication package at
+ * http://dx.doi.org/10.15456/jae.2022314.0708799949
+ *
+ * The prior distribution follows the one originally specified in Schorfheide's
+ * paper. Note that the elicited beta prior for rho in the paper
+ * implies an asymptote and corresponding prior mode at 0. It is generally
+ * recommended to avoid this extreme type of prior.
+ *
+ * Because the data are already logged and we use the loglinear option to conduct
+ * a full log-linearization, we need to use the logdata option.
+ *
+ * The equations are taken from J. Nason and T. Cogley (1994): "Testing the
+ * implications of long-run neutrality for monetary business cycle models",
+ * Journal of Applied Econometrics, 9, S37-S70, NC in the following.
+ * Note that there is an initial minus sign missing in equation (A1), p. S63.
+ *
+ * This implementation was originally written by Michel Juillard. Please note that the
+ * following copyright notice only applies to this Dynare implementation of the
+ * model.
+ */
+
+/*
+ * Copyright © 2004-2023 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 <https://www.gnu.org/licenses/>.
+ */
+
+var m ${m}$ (long_name='money growth')
+ P ${P}$ (long_name='Price level')
+ c ${c}$ (long_name='consumption')
+ e ${e}$ (long_name='capital stock')
+ W ${W}$ (long_name='Wage rate')
+ R ${R}$ (long_name='interest rate')
+ k ${k}$ (long_name='capital stock')
+ d ${d}$ (long_name='dividends')
+ n ${n}$ (long_name='labor')
+ l ${l}$ (long_name='loans')
+ gy_obs ${\Delta \ln GDP}$ (long_name='detrended capital stock')
+ gp_obs ${\Delta \ln P}$ (long_name='detrended capital stock')
+ y ${y}$ (long_name='detrended output')
+ dA ${\Delta A}$ (long_name='TFP growth')
+ ;
+varexo e_a ${\epsilon_A}$ (long_name='TFP shock')
+ e_m ${\epsilon_M}$ (long_name='Money growth shock')
+ ;
+
+parameters alp ${\alpha}$ (long_name='capital share')
+ bet ${\beta}$ (long_name='discount factor')
+ gam ${\gamma}$ (long_name='long-run TFP growth')
+ logmst ${\log(m^*)}$ (long_name='long-run money growth')
+ rho ${\rho}$ (long_name='autocorrelation money growth')
+ phi ${\phi}$ (long_name='labor weight in consumption')
+ del ${\delta}$ (long_name='depreciation rate')
+ ;
+
+% roughly picked values to allow simulating the model before estimation
+alp = 0.33;
+bet = 0.99;
+gam = 0.003;
+logmst = log(1.011);
+rho = 0.7;
+phi = 0.787;
+del = 0.02;
+
+model;
+[name='NC before eq. (1), TFP growth equation']
+dA = exp(gam+e_a);
+[name='NC eq. (2), money growth rate']
+log(m) = (1-rho)*logmst + rho*log(m(-1))+e_m;
+[name='NC eq. (A1), Euler equation']
+-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
+[name='NC below eq. (A1), firm borrowing constraint']
+W = l/n;
+[name='NC eq. (A2), intratemporal labour market condition']
+-(phi/(1-phi))*(c*P/(1-n))+l/n = 0;
+[name='NC below eq. (A2), credit market clearing']
+R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
+[name='NC eq. (A3), credit market optimality']
+1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
+[name='NC eq. (18), aggregate resource constraint']
+c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
+[name='NC eq. (19), money market condition']
+P*c = m;
+[name='NC eq. (20), credit market equilibrium condition']
+m-1+d = l;
+[name='Definition TFP shock']
+e = exp(e_a);
+[name='Implied by NC eq. (18), production function']
+y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
+[name='Observation equation GDP growth']
+gy_obs = dA*y/y(-1);
+[name='Observation equation price level']
+gp_obs = (P/P(-1))*m(-1)/dA;
+end;
+
+shocks;
+var e_a; stderr 0.014;
+var e_m; stderr 0.005;
+end;
+
+steady_state_model;
+ dA = exp(gam);
+ gst = 1/dA;
+ m = exp(logmst);
+ khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1));
+ xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/m )^(-1);
+ nust = phi*m^2/( (1-alp)*(1-phi)*bet*gst^alp*khst^alp );
+ n = xist/(nust+xist);
+ P = xist + nust;
+ k = khst*n;
+
+ l = phi*m*n/( (1-phi)*(1-n) );
+ c = m/P;
+ d = l - m + 1;
+ y = k^alp*n^(1-alp)*gst^alp;
+ R = m/bet;
+ W = l/n;
+ ist = y-c;
+ q = 1 - d;
+
+ e = 1;
+
+ gp_obs = m/dA;
+ gy_obs = dA;
+end;
+
+steady;
+check;
+
+% Table 1 of Schorfheide (2000)
+estimated_params;
+alp, beta_pdf, 0.356, 0.02;
+bet, beta_pdf, 0.993, 0.002;
+gam, normal_pdf, 0.0085, 0.003;
+logmst, normal_pdf, 0.0002, 0.007;
+rho, beta_pdf, 0.129, 0.223;
+phi, beta_pdf, 0.65, 0.05;
+del, beta_pdf, 0.01, 0.005;
+stderr e_a, inv_gamma_pdf, 0.035449, inf;
+stderr e_m, inv_gamma_pdf, 0.008862, inf;
+end;
+
+varobs gp_obs gy_obs;
+
+@#define T = 200
+
+set_dynare_seed(123456);
+stoch_simul(order=1,periods=1000,loglinear,noprint,nograph);
+ds = dseries(oo_.endo_simul','1Y',M_.endo_names);
+data(series=ds,first_obs=801Y,nobs=@{T});
+
+estimation(order=1, loglinear, logdata, mode_compute=99, optim=('delta', 0.9, 'first_mode_compute', 4), mh_replic=0, mode_check, cova_compute=0, silent_optimizer);
--
GitLab