diff --git a/matlab/+mom/Jtest.m b/matlab/+mom/Jtest.m index d455b3754b0c2fdb9949a5d52bce22bdba1b22df..955756fb1dc78dff0f9389f3972dc9d6bf34c1b9 100644 --- a/matlab/+mom/Jtest.m +++ b/matlab/+mom/Jtest.m @@ -61,11 +61,11 @@ if options_mom_.mom.mom_nbr > length(xparam) end % Compute J statistic if strcmp(options_mom_.mom.mom_method,'SMM') - Variance_correction_factor = options_mom_.mom.variance_correction_factor; + variance_correction_factor = options_mom_.mom.variance_correction_factor; elseif strcmp(options_mom_.mom.mom_method,'GMM') - Variance_correction_factor = 1; + variance_correction_factor = 1; end - J_test.j_stat = options_mom_.nobs*Variance_correction_factor*fval/options_mom_.mom.weighting_matrix_scaling_factor; + J_test.j_stat = options_mom_.nobs*variance_correction_factor*fval/options_mom_.mom.weighting_matrix_scaling_factor; J_test.degrees_freedom = length(model_moments)-length(xparam); J_test.p_val = 1-chi2cdf(J_test.j_stat, J_test.degrees_freedom); fprintf('\nValue of J-test statistic: %f\n',J_test.j_stat); diff --git a/matlab/+mom/default_option_mom_values.m b/matlab/+mom/default_option_mom_values.m index f465f76a08831e9ff7121b35eecea6adeede31e4..9b390ecb321c06ed306d9cac988ee3784c1468a4 100644 --- a/matlab/+mom/default_option_mom_values.m +++ b/matlab/+mom/default_option_mom_values.m @@ -1,5 +1,5 @@ -function options_mom_ = default_option_mom_values(options_mom_, options_, dname, doBayesianEstimation) -% options_mom_ = default_option_mom_values(options_mom_, options_, dname, doBayesianEstimation) +function options_mom_ = default_option_mom_values(options_mom_, options_, dname, do_bayesian_estimation) +% options_mom_ = default_option_mom_values(options_mom_, options_, dname, do_bayesian_estimation) % ------------------------------------------------------------------------- % Returns structure containing the options for method_of_moments command. % Note 1: options_mom_ is local and contains default and user-specified @@ -13,10 +13,10 @@ function options_mom_ = default_option_mom_values(options_mom_, options_, dname, % computed from the model and the moments/irfs computed from the data. % ------------------------------------------------------------------------- % INPUTS -% o options_mom_: [structure] all user-specified settings (from the method_of_moments command) -% o options_: [structure] global options -% o dname: [string] default name of directory to store results -% o doBayesianEstimation [boolean] indicator whether we do Bayesian estimation +% o options_mom_: [structure] all user-specified settings (from the method_of_moments command) +% o options_: [structure] global options +% o dname: [string] default name of directory to store results +% o do_bayesian_estimation [boolean] indicator whether we do Bayesian estimation % ------------------------------------------------------------------------- % OUTPUTS % o options_mom_: [structure] all user-specified and updated settings required for method_of_moments estimation @@ -96,7 +96,7 @@ options_mom_ = set_default_option(options_mom_,'nograph',false); % do no options_mom_ = set_default_option(options_mom_,'noprint',false); % do not print output to console options_mom_ = set_default_option(options_mom_,'TeX',false); % print TeX tables and graphics options_mom_.mom = set_default_option(options_mom_.mom,'verbose',false); % display and store intermediate estimation results -if doBayesianEstimation +if do_bayesian_estimation options_mom_ = set_default_option(options_mom_,'plot_priors',true); % control plotting of priors options_mom_ = set_default_option(options_mom_,'prior_trunc',1e-10); % probability of extreme values of the prior density that is ignored when computing bounds for the parameters end @@ -193,7 +193,7 @@ options_mom_ = set_default_option(options_mom_,'lyapunov_fixed_point_tol',1e-10) options_mom_ = set_default_option(options_mom_,'lyapunov_doubling_tol',1e-16); % convergence criterion used in the doubling algorithm % Bayesian MCMC related -if doBayesianEstimation +if do_bayesian_estimation options_mom_ = set_default_option(options_mom_,'mh_replic',0); % number of draws in Metropolis-Hastings and slice samplers options_mom_ = set_default_option(options_mom_,'mh_posterior_mode_estimation',false); % skip optimizer-based mode-finding and instead compute the mode based on a run of a MCMC options_mom_ = set_default_option(options_mom_,'load_mh_file',false); % add to previous Metropolis-Hastings or slice simulations instead of starting from scratch @@ -448,7 +448,7 @@ if strcmp(mom_method,'GMM') error('method_of_moments: Perturbation orders higher than 3 are not implemented for GMM estimation, try using SMM!'); end end -if strcmp(mom_method,'IRF_MATCHING') && doBayesianEstimation +if strcmp(mom_method,'IRF_MATCHING') && do_bayesian_estimation if isfield(options_mom_,'mh_tune_jscale') && options_mom_.mh_tune_jscale.status && (options_mom_.mh_tune_jscale.maxiter<options_mom_.mh_tune_jscale.stepsize) warning('method_of_moments: You specified mh_tune_jscale, but the maximum number of iterations is smaller than the step size. No update will take place.') end diff --git a/matlab/+mom/get_data_moments.m b/matlab/+mom/get_data_moments.m index 16b128813ad414bcade16aca42b0aea8fe41b3b7..e00f6de2ed5922dd342357c68c497a3d04336b89 100644 --- a/matlab/+mom/get_data_moments.m +++ b/matlab/+mom/get_data_moments.m @@ -1,5 +1,5 @@ -function [dataMoments, m_data] = get_data_moments(data, obs_var, inv_order_var, matched_moments_, options_mom_) -% [dataMoments, m_data] = get_data_moments(data, obs_var, inv_order_var, matched_moments_, options_mom_) +function [data_moments, m_data] = get_data_moments(data, obs_var, inv_order_var, matched_moments_, options_mom_) +% [data_moments, m_data] = get_data_moments(data, obs_var, inv_order_var, matched_moments_, options_mom_) % ------------------------------------------------------------------------- % Computes the user-selected empirical moments from data % ------------------------------------------------------------------------- @@ -11,7 +11,7 @@ function [dataMoments, m_data] = get_data_moments(data, obs_var, inv_order_var, % o options_mom_: [structure] information about all settings (specified by the user, preprocessor, and taken from global options_) % ------------------------------------------------------------------------- % OUTPUTS -% o dataMoments [numMom x 1] mean of selected empirical moments +% o data_moments [numMom x 1] mean of selected empirical moments % o m_data [T x numMom] selected empirical moments at each point in time % ------------------------------------------------------------------------- % This function is called by @@ -39,7 +39,7 @@ function [dataMoments, m_data] = get_data_moments(data, obs_var, inv_order_var, % Initialization T = size(data,1); % Number of observations (T) -dataMoments = NaN(options_mom_.mom.mom_nbr,1); +data_moments = NaN(options_mom_.mom.mom_nbr,1); m_data = NaN(T,options_mom_.mom.mom_nbr); % Product moment for each time period, i.e. each row t contains y_t1(l1)^p1*y_t2(l2)^p2*... % note that here we already are able to treat leads and lags and any power product moments @@ -59,10 +59,10 @@ for jm = 1:options_mom_.mom.mom_nbr end % We replace NaN (due to leads and lags and missing values) with the corresponding mean if isoctave && octave_ver_less_than('8') - dataMoments(jm,1) = nanmean(m_data_tmp); + data_moments(jm,1) = nanmean(m_data_tmp); else - dataMoments(jm,1) = mean(m_data_tmp,'omitnan'); + data_moments(jm,1) = mean(m_data_tmp,'omitnan'); end - m_data_tmp(isnan(m_data_tmp)) = dataMoments(jm,1); + m_data_tmp(isnan(m_data_tmp)) = data_moments(jm,1); m_data(:,jm) = m_data_tmp; end diff --git a/matlab/+mom/matched_irfs_blocks.m b/matlab/+mom/matched_irfs_blocks.m index 1b422478864e709918f426c8880f254839983609..810b6de110bedd91773b03178d5a1fbb0db2b47b 100644 --- a/matlab/+mom/matched_irfs_blocks.m +++ b/matlab/+mom/matched_irfs_blocks.m @@ -1,5 +1,5 @@ -function [data_irfs, weightMat, irfIndex, maxIrfHorizon] = matched_irfs_blocks(matched_irfs, matched_irfs_weight, varobs_id, obs_nbr, exo_nbr, endo_names) -% [data_irfs, weightMat, irfIndex, maxIrfHorizon] = matched_irfs_blocks(matched_irfs, matched_irfs_weight, varobs_id, obs_nbr, exo_nbr, endo_names) +function [data_irfs, weight_mat, irf_index, max_irf_horizon] = matched_irfs_blocks(matched_irfs, matched_irfs_weight, varobs_id, obs_nbr, exo_nbr, endo_names) +% [data_irfs, weight_mat, irf_index, max_irf_horizon] = matched_irfs_blocks(matched_irfs, matched_irfs_weight, varobs_id, obs_nbr, exo_nbr, endo_names) % ------------------------------------------------------------------------- % Checks and transforms matched_irfs and matched_irfs_weight blocks % for further use in the estimation. @@ -14,9 +14,9 @@ function [data_irfs, weightMat, irfIndex, maxIrfHorizon] = matched_irfs_blocks(m % ------------------------------------------------------------------------- % OUTPUT % data_irfs: [matrix] irfs for VAROBS as declared in matched_irfs block -% weightMat: [matrix] weighting matrix for irfs as declared in matched_irfs_weight block -% irfIndex: [vector] index for selecting specific irfs from full irf matrix of observables -% maxIrfHorizon: [scalar] maximum irf horizon as declared in matched_irfs block +% weight_mat: [matrix] weighting matrix for irfs as declared in matched_irfs_weight block +% irf_index: [vector] index for selecting specific irfs from full irf matrix of observables +% max_irf_horizon: [scalar] maximum irf horizon as declared in matched_irfs block % ------------------------------------------------------------------------- % This function is called by % o mom.run @@ -40,42 +40,42 @@ function [data_irfs, weightMat, irfIndex, maxIrfHorizon] = matched_irfs_blocks(m % along with Dynare. If not, see <https://www.gnu.org/licenses/>. -maxIrfHorizon = max(cellfun(@(x) x(end), matched_irfs(:,1))); % get maximum irf horizon +max_irf_horizon = max(cellfun(@(x) x(end), matched_irfs(:,1))); % get maximum irf horizon % create full matrix where 1st dimension are irf periods, 2nd dimension are variables as declared in VAROBS, 3rd dimension are shocks. -data_irfs = nan(maxIrfHorizon,obs_nbr,exo_nbr); +data_irfs = nan(max_irf_horizon,obs_nbr,exo_nbr); % overwrite nan values if they are declared in matched_irfs block; remaining nan values will be later ignored in the matching for jj = 1:size(matched_irfs,1) - idVar = matched_irfs{jj,1}(1); - idVarobs = find(varobs_id==idVar,1); - idShock = matched_irfs{jj,1}(2); - idIrfPeriod = matched_irfs{jj,1}(3); - irfValue = matched_irfs{jj,2}; - if isempty(idVarobs) + id_var = matched_irfs{jj,1}(1); + id_varobs = find(varobs_id==id_var,1); + id_shock = matched_irfs{jj,1}(2); + id_irf_period = matched_irfs{jj,1}(3); + irf_value = matched_irfs{jj,2}; + if isempty(id_varobs) skipline; - error('method_of_moments: You specified an irf matching involving variable %s, but it is not declared as a varobs!',endo_names{idVar}) + error('method_of_moments: You specified an irf matching involving variable %s, but it is not declared as a varobs!',endo_names{id_var}) end - data_irfs(idIrfPeriod,idVarobs,idShock) = irfValue; + data_irfs(id_irf_period,id_varobs,id_shock) = irf_value; end % create (full) empirical weighting matrix -weightMat = eye(maxIrfHorizon*obs_nbr*exo_nbr); % identity matrix by default: all irfs are equally important +weight_mat = eye(max_irf_horizon*obs_nbr*exo_nbr); % identity matrix by default: all irfs are equally important for jj = 1:size(matched_irfs_weight,1) - idVar1 = matched_irfs_weight{jj,1}(1); idVarobs1 = find(varobs_id==idVar1,1); idShock1 = matched_irfs_weight{jj,1}(2); idIrfPeriod1 = matched_irfs_weight{jj,1}(3); - idVar2 = matched_irfs_weight{jj,2}(1); idVarobs2 = find(varobs_id==idVar2,1); idShock2 = matched_irfs_weight{jj,2}(2); idIrfPeriod2 = matched_irfs_weight{jj,2}(3); - weightMatValue = matched_irfs_weight{jj,3}; - if isempty(idVarobs1) + id_var1 = matched_irfs_weight{jj,1}(1); id_varobs1 = find(varobs_id==id_var1,1); id_shock1 = matched_irfs_weight{jj,1}(2); id_irf_period1 = matched_irfs_weight{jj,1}(3); + id_var2 = matched_irfs_weight{jj,2}(1); id_varobs2 = find(varobs_id==id_var2,1); id_shock2 = matched_irfs_weight{jj,2}(2); id_irf_period2 = matched_irfs_weight{jj,2}(3); + weight_mat_value = matched_irfs_weight{jj,3}; + if isempty(id_varobs1) skipline; - error('method_of_moments: You specified a weight for an irf matching involving variable %s, but it is not a varobs!',endo_names{idVar1}) + error('method_of_moments: You specified a weight for an irf matching involving variable %s, but it is not a varobs!',endo_names{id_var1}) end - if isempty(idVarobs2)s + if isempty(id_varobs2) skipline; - error('method_of_moments: You specified a weight for an irf matching involving variable %s, but it is not a varobs!',endo_names{idVar2}) + error('method_of_moments: You specified a weight for an irf matching involving variable %s, but it is not a varobs!',endo_names{id_var2}) end - idweightMat1 = sub2ind(size(data_irfs),idIrfPeriod1,idVarobs1,idShock1); - idweightMat2 = sub2ind(size(data_irfs),idIrfPeriod2,idVarobs2,idShock2); - weightMat(idweightMat1,idweightMat2) = weightMatValue; - weightMat(idweightMat2,idweightMat1) = weightMatValue; % symmetry + idweight_mat1 = sub2ind(size(data_irfs),id_irf_period1,id_varobs1,id_shock1); + idweight_mat2 = sub2ind(size(data_irfs),id_irf_period2,id_varobs2,id_shock2); + weight_mat(idweight_mat1,idweight_mat2) = weight_mat_value; + weight_mat(idweight_mat2,idweight_mat1) = weight_mat_value; % symmetry end % focus only on specified irfs -irfIndex = find(~isnan(data_irfs)); -data_irfs = data_irfs(irfIndex); -weightMat = weightMat(irfIndex,irfIndex); \ No newline at end of file +irf_index = find(~isnan(data_irfs)); +data_irfs = data_irfs(irf_index); +weight_mat = weight_mat(irf_index,irf_index); \ No newline at end of file diff --git a/matlab/+mom/matched_moments_block.m b/matlab/+mom/matched_moments_block.m index c25e2740fdcbceb08de51093be7269fff367cd38..5d7a5154552461fe27426eec0c1ec232a8dceaf2 100644 --- a/matlab/+mom/matched_moments_block.m +++ b/matlab/+mom/matched_moments_block.m @@ -60,22 +60,22 @@ for jm = 1:size(matched_moments,1) end % find duplicate rows in cell array by making groups according to powers as we can then use cell2mat for the unique function powers = cellfun(@sum,matched_moments(:,3))'; -UniqueMomIdx = []; +unique_mom_idx = []; for jpow = unique(powers) idx1 = find(powers==jpow); [~,idx2] = unique(cell2mat(matched_moments(idx1,:)),'rows'); - UniqueMomIdx = [UniqueMomIdx idx1(idx2)]; + unique_mom_idx = [unique_mom_idx idx1(idx2)]; end % remove duplicate elements -DuplicateMoms = setdiff(1:size(matched_moments_orig,1),UniqueMomIdx); -if ~isempty(DuplicateMoms) - fprintf('Duplicate declared moments found and removed in ''matched_moments'' block in rows:\n %s.\n',num2str(DuplicateMoms)); +duplicate_moms = setdiff(1:size(matched_moments_orig,1),unique_mom_idx); +if ~isempty(duplicate_moms) + fprintf('Duplicate declared moments found and removed in ''matched_moments'' block in rows:\n %s.\n',num2str(duplicate_moms)); fprintf('Dynare will continue with remaining moment conditions\n'); end if strcmp(mom_method, 'SMM') % for SMM: keep the original structure, but get rid of duplicate moments - matched_moments = matched_moments_orig(sort(UniqueMomIdx),:); + matched_moments = matched_moments_orig(sort(unique_mom_idx),:); elseif strcmp(mom_method, 'GMM') % for GMM we use the transformed matched_moments structure - matched_moments = matched_moments(sort(UniqueMomIdx),:); + matched_moments = matched_moments(sort(unique_mom_idx),:); end \ No newline at end of file diff --git a/matlab/+mom/mode_compute_irf_matching.m b/matlab/+mom/mode_compute_irf_matching.m index 71bb8f9dfc8ed847e4954ab001dccf1ae05e7480..4ae7de813c2c3ff49694967831e672baee85bae6 100644 --- a/matlab/+mom/mode_compute_irf_matching.m +++ b/matlab/+mom/mode_compute_irf_matching.m @@ -78,7 +78,7 @@ for optim_iter = 1:length(options_mom_.optimizer_vec) if options_mom_.optimizer_vec{optim_iter}==0 hessian_xparam1_iter = hessian_xparam1; else - fprintf('computing hessian'); + fprintf('computing Hessian'); hessian_xparam1_iter = hessian(objective_function, xparam1, options_mom_.gstep,... data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); hessian_xparam1_iter = reshape(hessian_xparam1_iter, length(xparam1), length(xparam1)); diff --git a/matlab/+mom/objective_function.m b/matlab/+mom/objective_function.m index 53916fe7509a74d33c67d2e830a3cb5006a8de7b..e18d775016645366b3ce59cb5efd10d440cffedf 100644 --- a/matlab/+mom/objective_function.m +++ b/matlab/+mom/objective_function.m @@ -1,5 +1,5 @@ -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) +function [fval, info, exit_flag, df, junk_hessian, 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, junk_hessian, 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 method of moments estimation, % i.e. distance between model and data moments/irfs, possibly augmented with priors. @@ -23,7 +23,7 @@ function [fval, info, exit_flag, df, junkHessian, Q, model_moments, model_moment % - info: [vector] information on error codes and penalties % - exit_flag: [double] flag for exit status (0 if error, 1 if no error) % - df: [matrix] analytical jacobian of the moment difference (wrt paramters), currently for GMM only -% - junkHessian: [matrix] empty matrix required for optimizer interface (Hessian would typically go here) +% - junk_hessian: [matrix] empty matrix required for optimizer interface (Hessian would typically go here) % - 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) @@ -72,7 +72,7 @@ irf_model_varobs = []; model_moments_params_derivs = []; model_moments = []; Q = []; -junkHessian = []; +junk_hessian = []; df = []; % required to be empty by e.g. newrat if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM') if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_jacobian @@ -262,7 +262,7 @@ end 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); + model_irf = 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); @@ -285,11 +285,11 @@ if strcmp(options_mom_.mom.mom_method,'IRF_MATCHING') && strcmp(options_mom_.mom y_irf = 100*y_irf/cs(i,i); end end - modelIrf(:,:,i) = transpose(y_irf); + model_irf(:,:,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_, dr.ys); + [model_irf, check] = feval(str2func(options_mom_.mom.irf_matching_file.name), model_irf, M_, options_mom_, dr.ys); if check fval = Inf; info(1) = 180; @@ -301,7 +301,7 @@ if strcmp(options_mom_.mom.mom_method,'IRF_MATCHING') && strcmp(options_mom_.mom return end end - irf_model_varobs = modelIrf(:,options_mom_.varobs_id,:); % focus only on observables (this will be used later for plotting) + irf_model_varobs = model_irf(:,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 diff --git a/matlab/+mom/optimal_weighting_matrix.m b/matlab/+mom/optimal_weighting_matrix.m index 11459f08c9f26b1dc64edf973d14294c9d415645..271983c48c4958c5c2f7d4325049ff0031569521 100644 --- a/matlab/+mom/optimal_weighting_matrix.m +++ b/matlab/+mom/optimal_weighting_matrix.m @@ -18,7 +18,7 @@ function W_opt = optimal_weighting_matrix(m_data, moments, q_lag) % o mom.run.m % ------------------------------------------------------------------------- % This function calls: -% o CorrMatrix (embedded) +% o corr_matrix (embedded) % ------------------------------------------------------------------------- % Copyright © 2020-2023 Dynare Team @@ -43,36 +43,36 @@ function W_opt = optimal_weighting_matrix(m_data, moments, q_lag) [T,num_Mom] = size(m_data); % note that in m_data NaN values (due to leads or lags in matched_moments and missing data) were replaced by the mean % center around moments (could be either data_moments or model_moments) -h_Func = m_data - repmat(moments',T,1); +h_func = m_data - repmat(moments',T,1); % the required correlation matrices -GAMA_array = zeros(num_Mom,num_Mom,q_lag); -GAMA0 = Corr_Matrix(h_Func,T,num_Mom,0); +gamma_array = zeros(num_Mom,num_Mom,q_lag); +gamma0 = corr_matrix(h_func,T,num_Mom,0); if q_lag > 0 for ii=1:q_lag - GAMA_array(:,:,ii) = Corr_Matrix(h_Func,T,num_Mom,ii); + gamma_array(:,:,ii) = corr_matrix(h_func,T,num_Mom,ii); end end % the estimate of S -S = GAMA0; +S = gamma0; if q_lag > 0 for ii=1:q_lag - S = S + (1-ii/(q_lag+1))*(GAMA_array(:,:,ii) + GAMA_array(:,:,ii)'); + S = S + (1-ii/(q_lag+1))*(gamma_array(:,:,ii) + gamma_array(:,:,ii)'); end end % the estimate of W W_opt = S\eye(size(S,1)); -W_opt=(W_opt+W_opt')/2; % ensure symmetry +W_opt = (W_opt+W_opt')/2; % ensure symmetry end % main function end % The correlation matrix -function GAMA_corr = Corr_Matrix(h_Func,T,num_Mom,v) - GAMA_corr = zeros(num_Mom,num_Mom); +function gamma_corr = corr_matrix(h_func,T,num_Mom,v) + gamma_corr = zeros(num_Mom,num_Mom); for t = 1+v:T - GAMA_corr = GAMA_corr + h_Func(t-v,:)'*h_Func(t,:); + gamma_corr = gamma_corr + h_func(t-v,:)'*h_func(t,:); end - GAMA_corr = GAMA_corr/T; -end % Corr_Matrix end \ No newline at end of file + gamma_corr = gamma_corr/T; +end % corr_matrix end \ No newline at end of file diff --git a/matlab/+mom/print_info_on_estimation_settings.m b/matlab/+mom/print_info_on_estimation_settings.m index 7aa843a7a00e695bef15a0b5dbaa7d9abcc06fe5..5d6ad6bfff8596d5cdc83acf3e10806e80742591 100644 --- a/matlab/+mom/print_info_on_estimation_settings.m +++ b/matlab/+mom/print_info_on_estimation_settings.m @@ -1,12 +1,12 @@ -function print_info_on_estimation_settings(options_mom_, number_of_estimated_parameters, doBayesianEstimation) -% print_info_on_estimation_settings(options_mom_, number_of_estimated_parameters, doBayesianEstimation) +function print_info_on_estimation_settings(options_mom_, number_of_estimated_parameters, do_bayesian_estimation) +% print_info_on_estimation_settings(options_mom_, number_of_estimated_parameters, do_bayesian_estimation) % ------------------------------------------------------------------------- % Print information on the method of moments estimation settings to the console % ------------------------------------------------------------------------- % INPUTS % options_mom_ [struct] options for the method of moments estimation % number_of_estimated_parameters [integer] number of estimated parameters -% doBayesianEstimation [boolean] true if the estimation is Bayesian +% do_bayesian_estimation [boolean] true if the estimation is Bayesian % ------------------------------------------------------------------------- % OUTPUT % No output, just displays the chosen settings @@ -53,7 +53,7 @@ if strcmp(options_mom_.mom.mom_method,'SMM') || strcmp(options_mom_.mom.mom_meth end end if strcmp(options_mom_.mom.mom_method,'IRF_MATCHING') - if doBayesianEstimation + if do_bayesian_estimation fprintf('Bayesian Impulse Response Function Matching with'); else fprintf('Frequentist Impulse Response Function Matching with'); diff --git a/matlab/+mom/run.m b/matlab/+mom/run.m index e762e63c487e5d0b65eef9eeb27156b345ebf7aa..14e98e375f01c96f736d59e37e5ea1c60cbf3f7f 100644 --- a/matlab/+mom/run.m +++ b/matlab/+mom/run.m @@ -185,11 +185,11 @@ end if (~isempty(estim_params_.var_endo) || ~isempty(estim_params_.corrn)) && strcmp(options_mom_.mom.mom_method, 'GMM') error('method_of_moments: GMM estimation does not support measurement error(s) yet. Please specify them as a structural shock!'); end -doBayesianEstimation = [estim_params_.var_exo(:,5); estim_params_.var_endo(:,5); estim_params_.corrx(:,6); estim_params_.corrn(:,6); estim_params_.param_vals(:,5)]; -if all(doBayesianEstimation~=0) - doBayesianEstimation = true; -elseif all(doBayesianEstimation==0) - doBayesianEstimation = false; +do_bayesian_estimation = [estim_params_.var_exo(:,5); estim_params_.var_endo(:,5); estim_params_.corrx(:,6); estim_params_.corrn(:,6); estim_params_.param_vals(:,5)]; +if all(do_bayesian_estimation~=0) + do_bayesian_estimation = true; +elseif all(do_bayesian_estimation==0) + do_bayesian_estimation = false; else error('method_of_moments: Estimation must be either fully Frequentist or fully Bayesian. Maybe you forgot to specify a prior distribution!'); end @@ -208,7 +208,7 @@ check_varobs_are_endo_and_declared_once(options_.varobs,M_.endo_names); % The idea is to be independent of options_ and have full control of the % estimation instead of possibly having to deal with options chosen somewhere % else in the mod file. -options_mom_ = mom.default_option_mom_values(options_mom_, options_, M_.dname, doBayesianEstimation); +options_mom_ = mom.default_option_mom_values(options_mom_, options_, M_.dname, do_bayesian_estimation); % ------------------------------------------------------------------------- @@ -234,12 +234,12 @@ CheckPath(M_.dname,'.'); CheckPath('method_of_moments',M_.dname); CheckPath('graphs',M_.dname); -if doBayesianEstimation +if do_bayesian_estimation oo_.mom.posterior.optimization.mode = []; oo_.mom.posterior.optimization.Variance = []; oo_.mom.posterior.optimization.log_density=[]; end -doBayesianEstimationMCMC = doBayesianEstimation && ( (options_mom_.mh_replic>0) || options_mom_.load_mh_file ); +do_bayesian_estimation_mcmc = do_bayesian_estimation && ( (options_mom_.mh_replic>0) || options_mom_.load_mh_file ); invhess = []; % decision rule oo_.dr = set_state_space(oo_.dr,M_); % get state-space representation @@ -314,7 +314,7 @@ if exist([M_.fname '_prior_restrictions.m'],'file') options_mom_.prior_restrictions.routine = str2func([M_.fname '_prior_restrictions']); end % check that the provided mode_file is compatible with the current estimation settings -if ~isempty(options_mom_.mode_file) && ( ~doBayesianEstimation || (doBayesianEstimation && ~options_mom_.mh_posterior_mode_estimation) ) +if ~isempty(options_mom_.mode_file) && ( ~do_bayesian_estimation || (do_bayesian_estimation && ~options_mom_.mh_posterior_mode_estimation) ) [xparam0, hessian_xparam0] = check_mode_file(xparam0, hessian_xparam0, options_mom_, bayestopt_); end % check on specified priors and penalized estimation (which uses Laplace approximated priors) @@ -340,18 +340,18 @@ else estim_params_.full_calibration_detected = false; end if options_mom_.use_calibration_initialization % set calibration as starting values - if ~isempty(bayestopt_) && ~doBayesianEstimation && any(all(isnan([xparam_calib xparam0]),2)) + if ~isempty(bayestopt_) && ~do_bayesian_estimation && any(all(isnan([xparam_calib xparam0]),2)) error('method_of_moments: When using the use_calibration option with %s without prior, the parameters must be explicitly initialized!',options_mom_.mom.mom_method); else [xparam0,estim_params_] = do_parameter_initialization(estim_params_,xparam_calib,xparam0); % get explicitly initialized parameters that have precedence over calibrated values end end % check initialization -if ~isempty(bayestopt_) && ~doBayesianEstimation && any(isnan(xparam0)) +if ~isempty(bayestopt_) && ~do_bayesian_estimation && any(isnan(xparam0)) error('method_of_moments: Frequentist %s requires all estimated parameters to be initialized, either in an estimated_params or estimated_params_init-block!',options_mom_.mom.mom_method); end % set and check parameter bounds -if ~isempty(bayestopt_) && doBayesianEstimation +if ~isempty(bayestopt_) && do_bayesian_estimation % plot prior densities if ~options_mom_.nograph && options_mom_.plot_priors if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM') @@ -402,7 +402,7 @@ if estim_params_.ncx || any(nnz(tril(M_.Correlation_matrix,-1))) || isfield(esti M_.sigma_e_is_diagonal = false; end % storing prior parameters in results structure -if doBayesianEstimation || ( (strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')) && options_mom_.mom.penalized_estimator) +if do_bayesian_estimation || ( (strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')) && options_mom_.mom.penalized_estimator) oo_.mom.prior.mean = bayestopt_.p1; oo_.mom.prior.mode = bayestopt_.p5; oo_.mom.prior.variance = diag(bayestopt_.p2.^2); @@ -414,7 +414,7 @@ M_ = set_all_parameters(xparam0,estim_params_,M_); % provide warning if there is NaN in parameters test_for_deep_parameters_calibration(M_); % set jscale -if doBayesianEstimationMCMC +if do_bayesian_estimation_mcmc if ~strcmp(options_mom_.posterior_sampler_options.posterior_sampling_method,'slice') if isempty(options_mom_.mh_jscale) options_mom_.mh_jscale = 2.38/sqrt(number_of_estimated_parameters); % use optimal value for univariate normal distribution (check_posterior_sampler_options and mode_compute=6 may overwrite this setting) @@ -423,7 +423,7 @@ if doBayesianEstimationMCMC end end % initialization of posterior sampler options -if doBayesianEstimationMCMC +if do_bayesian_estimation_mcmc [current_options, options_mom_, bayestopt_] = check_posterior_sampler_options([], M_.fname, M_.dname, options_mom_, BoundsInfo, bayestopt_); options_mom_.posterior_sampler_options.current_options = current_options; if strcmp(current_options.posterior_sampling_method,'slice') && current_options.use_mh_covariance_matrix && ~current_options.rotated @@ -431,7 +431,7 @@ if doBayesianEstimationMCMC end end % warning if prior allows that stderr parameters are negative or corr parameters are outside the unit circle -if doBayesianEstimation +if do_bayesian_estimation check_prior_stderr_corr(estim_params_,bayestopt_); % check value of prior density [~,~,~,info] = priordens(xparam0,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4); @@ -578,7 +578,7 @@ end % ------------------------------------------------------------------------- % print some info to console % ------------------------------------------------------------------------- -mom.print_info_on_estimation_settings(options_mom_, number_of_estimated_parameters, doBayesianEstimation); +mom.print_info_on_estimation_settings(options_mom_, number_of_estimated_parameters, do_bayesian_estimation); % ------------------------------------------------------------------------- @@ -592,8 +592,8 @@ end % compute mode for IRF matching % ------------------------------------------------------------------------- if strcmp(options_mom_.mom.mom_method,'IRF_MATCHING') - if ~doBayesianEstimation || (doBayesianEstimation && ~options_mom_.mh_posterior_mode_estimation) - [xparam1, hessian_xparam1, fval, oo_.mom.verbose] = mom.mode_compute_irf_matching(xparam0, hessian_xparam0, objective_function, doBayesianEstimation, oo_.mom.weighting_info, oo_.mom.data_moments, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state); + if ~do_bayesian_estimation || (do_bayesian_estimation && ~options_mom_.mh_posterior_mode_estimation) + [xparam1, hessian_xparam1, fval, oo_.mom.verbose] = mom.mode_compute_irf_matching(xparam0, hessian_xparam0, objective_function, do_bayesian_estimation, oo_.mom.weighting_info, oo_.mom.data_moments, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state); else xparam1 = xparam0; hessian_xparam1 = hessian_xparam0; @@ -618,8 +618,8 @@ if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_meth hessian_xparam1 = inv(invhess); end else - if ~doBayesianEstimation || ~options_mom_.mh_posterior_mode_estimation - if doBayesianEstimation + if ~do_bayesian_estimation || ~options_mom_.mh_posterior_mode_estimation + if do_bayesian_estimation oo_.mom.posterior.optimization.mode = xparam1; if exist('fval','var') oo_.mom.posterior.optimization.log_density = -fval; @@ -629,18 +629,18 @@ else hsd = sqrt(diag(hessian_xparam1)); % represent the curvature (or second derivatives) of the likelihood with respect to each parameter being estimated. invhess = inv(hessian_xparam1./(hsd*hsd'))./(hsd*hsd'); % before taking the inverse scale the Hessian matrix by dividing each of its elements by the outer product of hsd such that the diagonal of the resulting matrix is approximately 1. This kind of scaling can help in regularizing the matrix and potentially improves its condition number, which in turn can make the matrix inversion more stable. stdh = sqrt(diag(invhess)); - if doBayesianEstimation + if do_bayesian_estimation oo_.mom.posterior.optimization.Variance = invhess; end end else variances = bayestopt_.p2.*bayestopt_.p2; - idInf = isinf(variances); - variances(idInf) = 1; + id_inf = isinf(variances); + variances(id_inf) = 1; invhess = options_mom_.mh_posterior_mode_estimation*diag(variances); xparam1 = bayestopt_.p5; - idNaN = isnan(xparam1); - xparam1(idNaN) = bayestopt_.p1(idNaN); + id_nan = isnan(xparam1); + xparam1(id_nan) = bayestopt_.p1(id_nan); outside_bound_pars=find(xparam1 < BoundsInfo.lb | xparam1 > BoundsInfo.ub); xparam1(outside_bound_pars) = bayestopt_.p1(outside_bound_pars); end @@ -653,7 +653,7 @@ end % ------------------------------------------------------------------------- % display estimation results at mode % ------------------------------------------------------------------------- -if doBayesianEstimation && ~options_mom_.mom.penalized_estimator && ~options_mom_.mh_posterior_mode_estimation +if do_bayesian_estimation && ~options_mom_.mom.penalized_estimator && ~options_mom_.mh_posterior_mode_estimation % display table with Bayesian mode estimation results and store parameter estimates and standard errors in oo_ oo_.mom = display_estimation_results_table(xparam1, stdh, M_, options_mom_, estim_params_, bayestopt_, oo_.mom, prior_dist_names, 'Posterior', 'posterior'); % Laplace approximation to the marginal log density @@ -668,7 +668,7 @@ if doBayesianEstimation && ~options_mom_.mom.penalized_estimator && ~options_mom end fprintf('\nLog data density [Laplace approximation] is %f.\n',oo_.mom.MarginalDensity.LaplaceApproximation); end -elseif ~doBayesianEstimation || (doBayesianEstimation && options_mom_.mom.penalized_estimator) +elseif ~do_bayesian_estimation || (do_bayesian_estimation && options_mom_.mom.penalized_estimator) % display table with Frequentist estimation results and store parameter estimates and standard errors in oo_ oo_.mom = display_estimation_results_table(xparam1, stdh, M_, options_mom_, estim_params_, bayestopt_, oo_.mom, prior_dist_names, options_mom_.mom.mom_method, lower(options_mom_.mom.mom_method)); end @@ -677,11 +677,11 @@ end % ------------------------------------------------------------------------- % checks for mode and hessian at the mode % ------------------------------------------------------------------------- -if (~doBayesianEstimation && options_mom_.cova_compute) || (doBayesianEstimation && ~options_mom_.mh_posterior_mode_estimation && options_mom_.cova_compute) +if (~do_bayesian_estimation && options_mom_.cova_compute) || (do_bayesian_estimation && ~options_mom_.mh_posterior_mode_estimation && options_mom_.cova_compute) check_hessian_at_the_mode(hessian_xparam1, xparam1, M_, estim_params_, options_, BoundsInfo); end if options_mom_.mode_check.status - if ~doBayesianEstimation || (doBayesianEstimation && ~options_mom_.mh_posterior_mode_estimation) + if ~do_bayesian_estimation || (do_bayesian_estimation && ~options_mom_.mh_posterior_mode_estimation) mode_check(objective_function, xparam1, diag(stdh), options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, true,... % use diag(stdh) instead of hessian_xparam1 as mode_check uses diagonal elements oo_.mom.data_moments, oo_.mom.weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state); end @@ -690,7 +690,7 @@ end % ------------------------------------------------------------------------- % Bayesian MCMC estimation % ------------------------------------------------------------------------- -if doBayesianEstimationMCMC +if do_bayesian_estimation_mcmc invhess = set_mcmc_jumping_covariance(invhess, length(xparam1), options_mom_.MCMC_jumping_covariance, bayestopt_, 'method_of_moments'); % reset bounds as lb and ub must only be operational during mode-finding BoundsInfo = set_mcmc_prior_bounds(xparam1, bayestopt_, options_mom_, 'method_of_moments'); @@ -801,20 +801,20 @@ if doBayesianEstimationMCMC end % MAYBE USEFUL???? % % Posterior correlations - % ExtremeCorrBound=0.7; - % if ~isnan(ExtremeCorrBound) + % extreme_corr_bound = 0.7; + % if ~isnan(extreme_corr_bound) % tril_para_correlation_matrix=tril(para_correlation_matrix,-1); - % [RowIndex,ColIndex]=find(abs(tril_para_correlation_matrix)>ExtremeCorrBound); - % ExtremeCorrParams=cell(length(RowIndex),3); - % for i=1:length(RowIndex) - % ExtremeCorrParams{i,1}=char(parameter_names(RowIndex(i),:)); - % ExtremeCorrParams{i,2}=char(parameter_names(ColIndex(i),:)); - % ExtremeCorrParams{i,3}=tril_para_correlation_matrix(RowIndex(i),ColIndex(i)); + % [row_index,col_index]=find(abs(tril_para_correlation_matrix)>extreme_corr_bound); + % extreme_corr_params=cell(length(row_index),3); + % for i=1:length(row_index) + % extreme_corr_params{i,1}=char(parameter_names(row_index(i),:)); + % extreme_corr_params{i,2}=char(parameter_names(col_index(i),:)); + % extreme_corr_params{i,3}=tril_para_correlation_matrix(row_index(i),col_index(i)); % end % end % disp(' '); - % disp(['Correlations of Parameters (at Posterior Mode) > ',num2str(ExtremeCorrBound)]); - % disp(ExtremeCorrParams) + % disp(['Correlations of Parameters (at Posterior Mode) > ',num2str(extreme_corr_bound)]); + % disp(extreme_corr_params) end diff --git a/matlab/+mom/standard_errors.m b/matlab/+mom/standard_errors.m index 1eeed07c62be002dfb7a46935e3e91c00e364812..55c64d0d28c3c2b2c22551f533202e5fcd0ba06e 100644 --- a/matlab/+mom/standard_errors.m +++ b/matlab/+mom/standard_errors.m @@ -1,5 +1,5 @@ -function [SE_values, Asympt_Var] = standard_errors(xparam, objective_function, model_moments, model_moments_params_derivs, m_data, data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state) -% [SE_values, Asympt_Var] = standard_errors(xparam, objective_function, model_moments, model_moments_params_derivs, m_data, data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state) +function [stderr_values, asympt_cov_mat] = standard_errors(xparam, objective_function, model_moments, model_moments_params_derivs, m_data, data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state) +% [stderr_values, asympt_cov_mat] = standard_errors(xparam, objective_function, model_moments, model_moments_params_derivs, m_data, data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state) % ------------------------------------------------------------------------- % This function computes standard errors to the method of moments estimates % Adapted from replication codes of Andreasen, Fernández-Villaverde, Rubio-Ramírez (2018): @@ -25,8 +25,8 @@ function [SE_values, Asympt_Var] = standard_errors(xparam, objective_function, m % - exo_det_steady_state: [vector] steady state value for exogenous deterministic variables (initval) % ------------------------------------------------------------------------- % OUTPUTS -% o SE_values [nparam x 1] vector of standard errors -% o Asympt_Var [nparam x nparam] asymptotic covariance matrix +% o stderr_values [nparam x 1] vector of standard errors +% o asympt_cov_mat [nparam x nparam] asymptotic covariance matrix % ------------------------------------------------------------------------- % This function is called by % o mom.run.m @@ -70,8 +70,8 @@ if strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_standa fprintf('No standard errors available for parameter %s\n',get_the_name(i,options_mom_.TeX, M_, estim_params_, options_mom_.varobs)) end warning('There are NaN in the analytical Jacobian of Moments. Check your bounds and/or priors, or use a different optimizer.') - Asympt_Var = NaN(length(xparam),length(xparam)); - SE_values = NaN(length(xparam),1); + asympt_cov_mat = NaN(length(xparam),length(xparam)); + stderr_values = NaN(length(xparam),1); return end else @@ -89,22 +89,22 @@ else if nnz(info_p)==0 && nnz(info_m)==0 D(:,i) = (model_moments_p - model_moments_m)/(2*eps_value); else - problpar = get_the_name(i,options_mom_.TeX, M_, estim_params_, options_mom_.varobs); + problematic_parameter = get_the_name(i,options_mom_.TeX, M_, estim_params_, options_mom_.varobs); if info_p(1)==42 - warning('method_of_moments:info','Cannot compute the Jacobian using finite differences for parameter %s due to hitting the upper bound - no standard errors available.\n',problpar) + warning('method_of_moments:info','Cannot compute the Jacobian using finite differences for parameter %s due to hitting the upper bound - no standard errors available.\n',problematic_parameter) else message_p = get_error_message(info_p, options_mom_); end if info_m(1)==41 - warning('method_of_moments:info','Cannot compute the Jacobian using finite differences for parameter %s due to hitting the lower bound - no standard errors available.\n',problpar) + warning('method_of_moments:info','Cannot compute the Jacobian using finite differences for parameter %s due to hitting the lower bound - no standard errors available.\n',problematic_parameter) else message_m = get_error_message(info_m, options_mom_); end if info_m(1)~=41 && info_p(1)~=42 - warning('method_of_moments:info','Cannot compute the Jacobian using finite differences for parameter %s - no standard errors available\n %s %s\nCheck your priors or use a different optimizer.\n',problpar, message_p, message_m) + warning('method_of_moments:info','Cannot compute the Jacobian using finite differences for parameter %s - no standard errors available\n %s %s\nCheck your priors or use a different optimizer.\n',problematic_parameter, message_p, message_m) end - Asympt_Var = NaN(length(xparam),length(xparam)); - SE_values = NaN(length(xparam),1); + asympt_cov_mat = NaN(length(xparam),length(xparam)); + stderr_values = NaN(length(xparam),1); return end end @@ -116,12 +116,12 @@ end WW = weighting_info.Sw'*weighting_info.Sw; if weighting_info.Woptflag % we already have the optimal weighting matrix - Asympt_Var = 1/T*((D'*WW*D)\eye(dim_params)); + asympt_cov_mat = 1/T*((D'*WW*D)\eye(dim_params)); else % we do not have the optimal weighting matrix yet WWopt = mom.optimal_weighting_matrix(m_data, model_moments, options_mom_.mom.bartlett_kernel_lag); S = WWopt\eye(size(WWopt,1)); AA = (D'*WW*D)\eye(dim_params); - Asympt_Var = 1/T*AA*D'*WW*S*WW*D*AA; + asympt_cov_mat = 1/T*AA*D'*WW*S*WW*D*AA; end -SE_values = sqrt(diag(Asympt_Var)); \ No newline at end of file +stderr_values = sqrt(diag(asympt_cov_mat)); \ No newline at end of file