diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m index 929079add12b208366a19ca4451c5d887a9cb333..a655489ad0fe3367acba8919b42023bd259a4407 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 e82e4389ddcb1e9669c851d9068a8734f33bb00d..2856abd508a5cd29d58dd4a9d4723ecd6280d637 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 a1cc385e3b7611f9d4ebb34231904d99450a41a1..a59fe9b3ae95b5a1401b4801dc7442c2726c24fb 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 0000000000000000000000000000000000000000..cc28dbc632c579cbccf1f863753ccb910e2780a7 --- /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);