method_of_moments: use underscore instead of camel case for variables
Willi Mutschler
parent
937ee0ef77
commit
c59daa6139
|
|
@ -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);
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
irf_index = find(~isnan(data_irfs));
|
||||
data_irfs = data_irfs(irf_index);
|
||||
weight_mat = weight_mat(irf_index,irf_index);
|
||||
|
|
@ -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
|
||||
|
|
@ -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));
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
gamma_corr = gamma_corr/T;
|
||||
end % corr_matrix end
|
||||
|
|
@ -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');
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -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));
|
||||
stderr_values = sqrt(diag(asympt_cov_mat));
|
||||
Loading…
Reference in New Issue