From 20ec0a6c9765f6d25ea4eee5066e1834406a54b5 Mon Sep 17 00:00:00 2001 From: Willi Mutschler Date: Mon, 9 Oct 2023 22:33:43 +0200 Subject: [PATCH] method_of_moments: add irf_matching to objective_function --- matlab/+mom/objective_function.m | 87 ++++++++++++++++++++++++++++++-- 1 file changed, 83 insertions(+), 4 deletions(-) diff --git a/matlab/+mom/objective_function.m b/matlab/+mom/objective_function.m index 1821cc566..2d82560a4 100644 --- a/matlab/+mom/objective_function.m +++ b/matlab/+mom/objective_function.m @@ -1,7 +1,8 @@ -function [fval, info, exit_flag, df, junkHessian, Q, model_moments, model_moments_params_derivs] = objective_function(xparam, data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state) -% [fval, info, exit_flag, df, junkHessian, Q, model_moments, model_moments_params_derivs] = objective_function(xparam, data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); +function [fval, info, exit_flag, df, junkHessian, Q, model_moments, model_moments_params_derivs, irf_model_varobs] = objective_function(xparam, data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state) +% [fval, info, exit_flag, df, junkHessian, Q, model_moments, model_moments_params_derivs, irf_model_varobs] = objective_function(xparam, data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state) % ------------------------------------------------------------------------- -% This function evaluates the objective function for a method of moments estimation +% This function evaluates the objective function for method of moments estimation, +% i.e. distance between model and data moments/irfs, possibly augmented with priors. % ------------------------------------------------------------------------- % INPUTS (same ones as in dsge_likelihood.m and dsge_var_likelihood.m) % - xparam: [vector] current value of estimated parameters as returned by set_prior() @@ -26,6 +27,7 @@ function [fval, info, exit_flag, df, junkHessian, Q, model_moments, model_moment % - Q: [double] value of the quadratic form of the moment difference % - model_moments: [vector] model moments % - model_moments_params_derivs: [matrix] analytical jacobian of the model moments wrt estimated parameters (currently for GMM only) +% - irf_model_varobs: [matrix] model irfs for observable variables (used for plotting matched irfs in mom.run) % ------------------------------------------------------------------------- % This function is called by % o mom.run @@ -33,8 +35,12 @@ function [fval, info, exit_flag, df, junkHessian, Q, model_moments, model_moment % ------------------------------------------------------------------------- % This function calls % o check_bounds_and_definiteness_estimation +% o get_lower_cholesky_covariance % o get_perturbation_params_derivs +% o getIrfShocksIndx +% o irf % o mom.get_data_moments +% o priordens % o pruned_state_space_system % o resol % o set_all_parameters @@ -62,6 +68,7 @@ function [fval, info, exit_flag, df, junkHessian, Q, model_moments, model_moment %------------------------------------------------------------------------------ % Initialization of the returned variables and others... %------------------------------------------------------------------------------ +irf_model_varobs = []; model_moments_params_derivs = []; model_moments = []; Q = []; @@ -247,12 +254,84 @@ if strcmp(options_mom_.mom.mom_method,'SMM') end +%------------------------------------------------------------------------------ +% IRF_MATCHING using STOCH_SIMUL: Compute irfs given model solution and Cholesky +% decomposition on shock covariance matrix; this resembles the core codes in +% stoch_simul +%------------------------------------------------------------------------------ +if strcmp(options_mom_.mom.mom_method,'IRF_MATCHING') && strcmp(options_mom_.mom.simulation_method,'STOCH_SIMUL') + cs = get_lower_cholesky_covariance(M_.Sigma_e,options_mom_.add_tiny_number_to_cholesky); + irf_shocks_indx = getIrfShocksIndx(M_, options_mom_); + modelIrf = nan(options_mom_.irf,M_.endo_nbr,M_.exo_nbr); + for i = irf_shocks_indx + if options_mom_.order>1 && options_mom_.relative_irf % normalize shock to 0.01 before IRF generation for GIRFs; multiply with 100 later + y_irf = irf(M_, options_mom_, dr, cs(:,i)./cs(i,i)/100, options_mom_.irf, options_mom_.drop, options_mom_.replic, options_mom_.order); + else % for linear model, rescaling is done later + y_irf = irf(M_, options_mom_, dr, cs(:,i), options_mom_.irf, options_mom_.drop, options_mom_.replic, options_mom_.order); + end + if any(any(isnan(y_irf))) && ~options_mom_.pruning && ~(options_mom_.order==1) + fprintf('\nirf_matching: The simulations conducted for generating IRFs to %s were explosive: Either reduce the shock size, use pruning, or set the approximation order to 1.\n', M_.exo_names{i}) + fval = Inf; + info(1) = 180; + info(4) = 0.1; + exit_flag = 0; + if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output + fval = ones(options_mom_.mom.mom_nbr,1)*options_mom_.huge_number; + end + return + end + if options_mom_.relative_irf + if options_mom_.order==1 % multiply with 100 for backward compatibility + y_irf = 100*y_irf/cs(i,i); + end + end + modelIrf(:,:,i) = transpose(y_irf); + end + % do transformations on model irfs if irf_matching_file is provided + if ~isempty(options_mom_.mom.irf_matching_file.name) + [modelIrf,check] = feval(str2func(options_mom_.mom.irf_matching_file.name),modelIrf, M_, options_mom_); + if check + fval = Inf; + info(1) = 180; + info(4) = 0.1; + exit_flag = 0; + if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output + fval = ones(options_mom_.mom.mom_nbr,1)*options_mom_.huge_number; + end + return + end + end + irf_model_varobs = modelIrf(:,options_mom_.varobs_id,:); % focus only on observables (this will be used later for plotting) + model_moments = irf_model_varobs(options_mom_.mom.irfIndex); % focus only on selected irf periods +end + + %-------------------------------------------------------------------------- % Compute quadratic target function %-------------------------------------------------------------------------- moments_difference = data_moments - model_moments; -if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM') +if strcmp(options_mom_.mom.mom_method,'IRF_MATCHING') + Q = transpose(moments_difference)*weighting_info.W*moments_difference; + % log-likelihood + lnlik = options_mom_.mom.mom_nbr/2*log(1/2/pi) - 1/2*weighting_info.Winv_logdet - 1/2*Q; + if isinf(lnlik) + fval = Inf; info(1) = 50; info(4) = 0.1; exit_flag = 0; + return + end + if isnan(lnlik) + fval = Inf; info(1) = 45; info(4) = 0.1; exit_flag = 0; + return + end + if imag(lnlik)~=0 + fval = Inf; info(1) = 46; info(4) = 0.1; exit_flag = 0; + return + end + % add log prior if necessary + lnprior = priordens(xparam,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4); + fval = - (lnlik + lnprior); + +elseif strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM') residuals = sqrt(options_mom_.mom.weighting_matrix_scaling_factor)*weighting_info.Sw*moments_difference; Q = residuals'*residuals; if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output