stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

method_of_moments: don't use oo_ as input

covariance-quadratic-approximation
Willi Mutschler 2023-12-13 10:20:22 +01:00
parent dedbb3be57
commit 37c7ca2d97
No known key found for this signature in database
GPG Key ID: 91E724BF17A73F6D
5 changed files with 203 additions and 184 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

62
matlab/+mom/Jtest.m
View File

@ -1,22 +1,28 @@
function oo_ = Jtest(xparam, objective_function, Woptflag, oo_, options_mom_, bayestopt_, BoundsInfo, estim_params_, M_, nobs)
% function oo_ = Jtest(xparam, objective_function, Woptflag, oo_, options_mom_, bayestopt_, BoundsInfo, estim_params_, M_, nobs)
function J_test = Jtest(xparam, objective_function, Q, model_moments, m_data, data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% J_test = Jtest(xparam, objective_function, Q, model_moments, m_data, data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% -------------------------------------------------------------------------
% Computes the J-test statistic and p-value given a GMM/SMM estimation
% -------------------------------------------------------------------------
% Computes the J-test statistic and p-value for a GMM/SMM estimation
% =========================================================================
% INPUTS
% xparam: [vector] estimated parameter vector
% objective_function: [function handle] objective function
% Woptflag: [logical] flag if optimal weighting matrix has already been computed
% oo_: [struct] results
% options_mom_: [struct] options
% bayestopt_: [struct] information on priors
% BoundsInfo: [struct] bounds on parameters
% estim_params_: [struct] information on estimated parameters
% M_: [struct] information on the model
% nobs: [scalar] number of observations
% xparam: [vector] estimated parameter vector
% objective_function: [function handle] objective function
% Q: [scalar] value of moments distance criterion
% model_moments: [vector] model moments
% m_data: [matrix] selected empirical moments at each point in time
% data_moments: [vector] empirical moments
% weighting_info: [struct] information on weighting matrix
% options_mom_: [struct] options
% M_: [struct] model information
% estim_params_: [struct] estimated parameters
% bayestopt_: [struct] info on prior distributions
% BoundsInfo: [struct] info bounds on parameters
% dr: [struct] reduced form model
% endo_steady_state: [vector] steady state of endogenous variables (initval)
% exo_steady_state: [vector] steady state of exogenous variables (initval)
% exo_det_steady_state: [vector] steady state of deterministic exogenous variables (initval)
% -------------------------------------------------------------------------
% OUTPUT
% oo_: [struct] updated results
% J_test: [struct] results of J test
% -------------------------------------------------------------------------
% This function is called by
% o mom.run
@ -24,7 +30,8 @@ function oo_ = Jtest(xparam, objective_function, Woptflag, oo_, options_mom_, ba
% This function calls
% o mom.objective_function
% o mom.optimal_weighting_matrix
% =========================================================================
% -------------------------------------------------------------------------
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
@ -41,17 +48,16 @@ function oo_ = Jtest(xparam, objective_function, Woptflag, oo_, options_mom_, ba
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
if options_mom_.mom.mom_nbr > length(xparam)
% Get optimal weighting matrix for J test, if necessary
if ~Woptflag
W_opt = mom.optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag);
oo_J = oo_;
oo_J.mom.Sw = chol(W_opt);
fval = feval(objective_function, xparam, BoundsInfo, oo_J, estim_params_, M_, options_mom_);
if ~weighting_info.Woptflag
W_opt = mom.optimal_weighting_matrix(m_data, model_moments, options_mom_.mom.bartlett_kernel_lag);
weighting_info.Sw = chol(W_opt);
fval = feval(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);
else
fval = oo_.mom.Q;
fval = Q;
end
% Compute J statistic
if strcmp(options_mom_.mom.mom_method,'SMM')
@ -59,9 +65,9 @@ if options_mom_.mom.mom_nbr > length(xparam)
elseif strcmp(options_mom_.mom.mom_method,'GMM')
Variance_correction_factor = 1;
end
oo_.mom.J_test.j_stat = nobs*Variance_correction_factor*fval/options_mom_.mom.weighting_matrix_scaling_factor;
oo_.mom.J_test.degrees_freedom = length(oo_.mom.model_moments)-length(xparam);
oo_.mom.J_test.p_val = 1-chi2cdf(oo_.mom.J_test.j_stat, oo_.mom.J_test.degrees_freedom);
fprintf('\nValue of J-test statistic: %f\n',oo_.mom.J_test.j_stat);
fprintf('p-value of J-test statistic: %f\n',oo_.mom.J_test.p_val);
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);
fprintf('p-value of J-test statistic: %f\n',J_test.p_val);
end

72
matlab/+mom/mode_compute_gmm_smm.m
View File

@ -1,25 +1,30 @@
function [xparam1, oo_, Woptflag] = mode_compute_gmm_smm(xparam0, objective_function, oo_, M_, options_mom_, estim_params_, bayestopt_, BoundsInfo)
% function [xparam1, oo_, Woptflag] = mode_compute_gmm_smm(xparam0, objective_function, oo_, M_, options_mom_, estim_params_, bayestopt_, BoundsInfo)
function [xparam1, weighting_info, mom_verbose] = mode_compute_gmm_smm(xparam0, objective_function, m_data, data_moments, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% [xparam1, weighting_info, mom_verbose] = mode_compute_gmm_smm(xparam0, objective_function, m_data, data_moments, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% -------------------------------------------------------------------------
% Iterated method of moments for GMM and SMM, computes the minimum of the
% objective function (distance between data moments and model moments)
% for a sequence of optimizers and GMM/SMM iterations with different
% weighting matrices.
% =========================================================================
% -------------------------------------------------------------------------
% INPUTS
% xparam0: [vector] vector of initialized parameters
% objective_function: [func handle] name of the objective function
% oo_: [structure] results
% M_: [structure] model information
% options_mom_: [structure] options
% estim_params_: [structure] information on estimated parameters
% bayestopt_: [structure] information on priors
% BoundsInfo: [structure] bounds for optimization
% xparam0: [vector] vector of initialized parameters
% objective_function: [func handle] name of the objective function
% m_data: [matrix] selected data moments at each point in time
% data_moments: [vector] vector of data moments
% options_mom_: [structure] options
% M_: [structure] model information
% estim_params_: [structure] information on estimated parameters
% bayestopt_: [structure] information on priors
% BoundsInfo: [structure] bounds for optimization
% dr: [structure] reduced form model
% endo_steady_state: [vector] steady state for endogenous variables (initval)
% exo_steady_state: [vector] steady state for exogenous variables (initval)
% exo_det_steady_state: [vector] steady state for exogenous deterministic variables (initval)
% -------------------------------------------------------------------------
% OUTPUT
% xparam1: [vector] mode of objective function
% oo_: [structure] updated results
% Woptflag: [logical] true if optimal weighting matrix was computed
% xparam1: [vector] mode of objective function
% weighting_info: [structure] information on weighting matrix
% mom_verbose: [structure] information on intermediate estimation results
% -------------------------------------------------------------------------
% This function is called by
% o mom.run
@ -29,7 +34,9 @@ function [xparam1, oo_, Woptflag] = mode_compute_gmm_smm(xparam0, objective_func
% o mom.display_estimation_results_table
% o dynare_minimize_objective
% o mom.objective_function
% =========================================================================
% o prior_dist_names
% -------------------------------------------------------------------------
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
@ -46,15 +53,16 @@ function [xparam1, oo_, Woptflag] = mode_compute_gmm_smm(xparam0, objective_func
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
mom_verbose = [];
if size(options_mom_.mom.weighting_matrix,1)>1 && ~(any(strcmpi('diagonal',options_mom_.mom.weighting_matrix)) || any(strcmpi('optimal',options_mom_.mom.weighting_matrix)))
fprintf('\nYou did not specify the use of an optimal or diagonal weighting matrix. There is no point in running an iterated method of moments.\n')
end
for stage_iter = 1:size(options_mom_.mom.weighting_matrix,1)
fprintf('Estimation stage %u\n',stage_iter);
Woptflag = false;
weighting_info.Woptflag = false;
switch lower(options_mom_.mom.weighting_matrix{stage_iter})
case 'identity_matrix'
fprintf(' - identity weighting matrix\n');
@ -63,20 +71,20 @@ for stage_iter = 1:size(options_mom_.mom.weighting_matrix,1)
fprintf(' - diagonal of optimal weighting matrix (Bartlett kernel with %d lags)\n', options_mom_.mom.bartlett_kernel_lag);
if stage_iter == 1
fprintf(' and using data-moments as initial estimate of model-moments\n');
weighting_matrix = diag(diag( mom.optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.data_moments, options_mom_.mom.bartlett_kernel_lag) ));
weighting_matrix = diag(diag( mom.optimal_weighting_matrix(m_data, data_moments, options_mom_.mom.bartlett_kernel_lag) ));
else
fprintf(' and using previous stage estimate of model-moments\n');
weighting_matrix = diag(diag( mom.optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag) ));
weighting_matrix = diag(diag( mom.optimal_weighting_matrix(m_data, model_moments, options_mom_.mom.bartlett_kernel_lag) ));
end
case 'optimal'
fprintf(' - optimal weighting matrix (Bartlett kernel with %d lags)\n', options_mom_.mom.bartlett_kernel_lag);
if stage_iter == 1
fprintf(' and using data-moments as initial estimate of model-moments\n');
weighting_matrix = mom.optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.data_moments, options_mom_.mom.bartlett_kernel_lag);
weighting_matrix = mom.optimal_weighting_matrix(m_data, data_moments, options_mom_.mom.bartlett_kernel_lag);
else
fprintf(' and using previous stage estimate of model-moments\n');
weighting_matrix = mom.optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag);
Woptflag = true;
weighting_matrix = mom.optimal_weighting_matrix(m_data, model_moments, options_mom_.mom.bartlett_kernel_lag);
weighting_info.Woptflag = true;
end
otherwise % user specified matrix in file
fprintf(' - user-specified weighting matrix\n');
@ -86,12 +94,12 @@ for stage_iter = 1:size(options_mom_.mom.weighting_matrix,1)
error(['method_of_moments: No matrix named ''weighting_matrix'' could be found in ',options_mom_.mom.weighting_matrix{stage_iter},'.mat !'])
end
[nrow, ncol] = size(weighting_matrix);
if ~isequal(nrow,ncol) || ~isequal(nrow,length(oo_.mom.data_moments)) %check if square and right size
error(['method_of_moments: ''weighting_matrix'' must be square and have ',num2str(length(oo_.mom.data_moments)),' rows and columns!'])
if ~isequal(nrow,ncol) || ~isequal(nrow,length(data_moments)) %check if square and right size
error(['method_of_moments: ''weighting_matrix'' must be square and have ',num2str(length(data_moments)),' rows and columns!'])
end
end
try % check for positive definiteness of weighting_matrix
oo_.mom.Sw = chol(weighting_matrix);
weighting_info.Sw = chol(weighting_matrix);
catch
error('method_of_moments: Specified ''weighting_matrix'' is not positive definite. Check whether your model implies stochastic singularity!')
end
@ -99,8 +107,8 @@ for stage_iter = 1:size(options_mom_.mom.weighting_matrix,1)
for optim_iter = 1:length(options_mom_.optimizer_vec)
options_mom_.current_optimizer = options_mom_.optimizer_vec{optim_iter};
if options_mom_.optimizer_vec{optim_iter} == 0
xparam1 = xparam0; % no minimization, evaluate objective at current values
fval = feval(objective_function, xparam1, BoundsInfo, oo_, estim_params_, M_, options_mom_);
xparam1 = xparam0; % no minimization, evaluate objective at current values
fval = feval(objective_function, xparam1, data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
else
if options_mom_.optimizer_vec{optim_iter} == 13
options_mom_.mom.vector_output = true;
@ -113,17 +121,21 @@ for stage_iter = 1:size(options_mom_.mom.weighting_matrix,1)
options_mom_.mom.compute_derivs = false;
end
[xparam1, fval] = dynare_minimize_objective(objective_function, xparam0, options_mom_.optimizer_vec{optim_iter}, options_mom_, [BoundsInfo.lb BoundsInfo.ub], bayestopt_.name, bayestopt_, [],...
BoundsInfo, oo_, estim_params_, M_, options_mom_);
data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
if options_mom_.mom.vector_output
fval = fval'*fval;
end
end
fprintf('\nStage %d Iteration %d: Value of minimized moment distance objective function: %12.10f.\n',stage_iter,optim_iter,fval)
if options_mom_.mom.verbose
oo_.mom = display_estimation_results_table(xparam1,NaN(size(xparam1)),M_,options_mom_,estim_params_,bayestopt_,oo_.mom,prior_dist_names,sprintf('%s (STAGE %d ITERATION %d) VERBOSE',options_mom_.mom.mom_method,stage_iter,optim_iter),sprintf('verbose_%s_stage_%d_iter_%d',lower(options_mom_.mom.mom_method),stage_iter,optim_iter));
fprintf('\n''verbose'' option: ')
std_via_invhessian_xparam1_iter = NaN(size(xparam1));
tbl_title_iter = sprintf('FREQUENTIST %s (STAGE %d ITERATION %d) VERBOSE',options_mom_.mom.mom_method,stage_iter,optim_iter);
field_name_iter = sprintf('%s_stage_%d_iter_%d',lower(options_mom_.mom.mom_method),stage_iter,optim_iter);
mom_verbose.(field_name_iter) = display_estimation_results_table(xparam1,std_via_invhessian_xparam1_iter,M_,options_mom_,estim_params_,bayestopt_,[],prior_dist_names,tbl_title_iter,field_name_iter);
end
xparam0 = xparam1;
end
options_mom_.vector_output = false;
[~, ~, ~,~,~, oo_] = feval(objective_function, xparam1, BoundsInfo, oo_, estim_params_, M_, options_mom_); % get oo_.mom.model_moments for iterated GMM/SMM to compute optimal weighting matrix
[~, ~, ~, ~, ~, ~, model_moments] = feval(objective_function, xparam1, data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); % get model_moments for iterated GMM/SMM to compute optimal weighting matrix
end

140
matlab/+mom/objective_function.m
View File

@ -1,27 +1,31 @@
function [fval, info, exit_flag, df, junkHessian, oo_, M_] = objective_function(xparam, BoundsInfo, oo_, estim_params_, M_, options_mom_)
% [fval, info, exit_flag, df, junk1, oo_, M_] = objective_function(xparam, BoundsInfo, oo_, estim_params_, M_, options_mom_)
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);
% -------------------------------------------------------------------------
% This function evaluates the objective function for method of moments estimation
% =========================================================================
% INPUTS
% o xparam: [vector] current value of estimated parameters as returned by set_prior()
% o BoundsInfo: [structure] containing parameter bounds
% o oo_: [structure] for results
% o estim_params_: [structure] describing the estimated_parameters
% o M_ [structure] describing the model
% o options_mom_: [structure] information about all settings (specified by the user, preprocessor, and taken from global options_)
% This function evaluates the objective function for a method of moments estimation
% -------------------------------------------------------------------------
% 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()
% - data_moments: [vector] data with moments/irfs to match (corresponds to dataset_ in likelihood-based estimation)
% - weighting_info: [structure] storing information on weighting matrices (corresponds to dataset_info in likelihood-based estimation)
% - options_mom_: [structure] information about all settings (specified by the user, preprocessor, and taken from global options_)
% - M_ [structure] model information
% - estim_params_: [structure] information from estimated_params block
% - bayestopt_: [structure] information on the prior distributions
% - BoundsInfo: [structure] parameter bounds
% - dr: [structure] reduced form model
% - endo_steady_state: [vector] steady state value for endogenous variables (initval)
% - exo_steady_state: [vector] steady state value for exogenous variables (initval)
% - exo_det_steady_state: [vector] steady state value for exogenous deterministic variables (initval)
% -------------------------------------------------------------------------
% OUTPUTS
% o fval: [double] value of the quadratic form of the moment difference (except for lsqnonlin, where this is done implicitly)
% o info: [vector] information on error codes and penalties
% o exit_flag: [double] flag for exit status (0 if error, 1 if no error)
% o df: [matrix] analytical jacobian of the moment difference (wrt paramters), currently for GMM only
% o junkHessian: [matrix] empty matrix required for optimizer interface (Hessian would typically go here)
% o oo_: [structure] results with the following updated fields:
% - oo_.mom.model_moments: [vector] model moments
% - oo_.mom.Q: [double] value of the quadratic form of the moment difference
% - oo_.mom.model_moments_params_derivs: [matrix] analytical jacobian of the model moments wrt estimated parameters (currently for GMM only)
% o M_: [structure] updated model structure
% - fval: [double] value of the quadratic form of the moment difference (except for lsqnonlin, where this is done implicitly)
% - 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)
% - 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)
% -------------------------------------------------------------------------
% This function is called by
% o mom.run
@ -35,7 +39,8 @@ function [fval, info, exit_flag, df, junkHessian, oo_, M_] = objective_function(
% o resol
% o set_all_parameters
% o simult_
% =========================================================================
% -------------------------------------------------------------------------
% Copyright © 2020-2023 Dynare Team
%
% This file is part of Dynare.
@ -52,25 +57,23 @@ function [fval, info, exit_flag, df, junkHessian, oo_, M_] = objective_function(
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
%% TO DO
% check the info values and make use of meaningful penalties
% how do we do the penalty for the prior??
%------------------------------------------------------------------------------
% Initialization of the returned variables and others...
%------------------------------------------------------------------------------
model_moments_params_derivs = [];
model_moments = [];
Q = [];
junkHessian = [];
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
if options_mom_.mom.vector_output == 1
if options_mom_.mom.penalized_estimator
df = nan(size(oo_.mom.data_moments,1)+length(xparam),length(xparam));
df = nan(options_mom_.mom.mom_nbr+length(xparam),length(xparam));
else
df = nan(size(oo_.mom.data_moments,1),length(xparam));
df = nan(options_mom_.mom.mom_nbr,length(xparam));
end
else
df = nan(length(xparam),1);
@ -82,14 +85,13 @@ end
%--------------------------------------------------------------------------
% Get the structural parameters and define penalties
%--------------------------------------------------------------------------
% Ensure that xparam1 is a column vector; particleswarm.m requires this.
xparam = xparam(:);
M_ = set_all_parameters(xparam, estim_params_, M_);
[fval,info,exit_flag] = check_bounds_and_definiteness_estimation(xparam, M_, estim_params_, BoundsInfo);
if info(1)
if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
fval = ones(size(oo_.mom.data_moments,1),1)*options_mom_.huge_number;
fval = ones(options_mom_.mom.mom_nbr,1)*options_mom_.huge_number;
end
return
end
@ -98,9 +100,8 @@ end
%--------------------------------------------------------------------------
% Call resol to compute steady state and model solution
%--------------------------------------------------------------------------
% Compute linear approximation around the deterministic steady state
[oo_.dr, info, M_.params] = resol(0, M_, options_mom_, oo_.dr ,oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
[dr, info, M_.params] = resol(0, M_, options_mom_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
% Return, with endogenous penalty when possible, if resol issues an error code
if info(1)
if info(1) == 3 || info(1) == 4 || info(1) == 5 || info(1)==6 ||info(1) == 19 ||...
@ -111,7 +112,7 @@ if info(1)
info(4) = info(2);
exit_flag = 0;
if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
fval = ones(size(oo_.mom.data_moments,1),1)*options_mom_.huge_number;
fval = ones(options_mom_.mom.mom_nbr,1)*options_mom_.huge_number;
end
return
else
@ -119,7 +120,7 @@ if info(1)
info(4) = 0.1;
exit_flag = 0;
if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
fval = ones(size(oo_.mom.data_moments,1),1)*options_mom_.huge_number;
fval = ones(options_mom_.mom.mom_nbr,1)*options_mom_.huge_number;
end
return
end
@ -150,38 +151,38 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
stderrparam_nbr = estim_params_.nvx; % number of stderr parameters
corrparam_nbr = estim_params_.ncx; % number of corr parameters
totparam_nbr = stderrparam_nbr+corrparam_nbr+modparam_nbr;
oo_.dr.derivs = identification.get_perturbation_params_derivs(M_, options_mom_, estim_params_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state, indpmodel, indpstderr, indpcorr, 0); %analytic derivatives of perturbation matrices
oo_.mom.model_moments_params_derivs = NaN(options_mom_.mom.mom_nbr,totparam_nbr);
pruned_state_space = pruned_SS.pruned_state_space_system(M_, options_mom_, oo_.dr, oo_.mom.obs_var, options_mom_.ar, 0, 1);
dr.derivs = identification.get_perturbation_params_derivs(M_, options_mom_, estim_params_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, indpmodel, indpstderr, indpcorr, 0); %analytic derivatives of perturbation matrices
model_moments_params_derivs = NaN(options_mom_.mom.mom_nbr,totparam_nbr);
pruned_state_space = pruned_SS.pruned_state_space_system(M_, options_mom_, dr, options_mom_.mom.obs_var, options_mom_.ar, 0, 1);
else
pruned_state_space = pruned_SS.pruned_state_space_system(M_, options_mom_, oo_.dr, oo_.mom.obs_var, options_mom_.ar, 0, 0);
pruned_state_space = pruned_SS.pruned_state_space_system(M_, options_mom_, dr, options_mom_.mom.obs_var, options_mom_.ar, 0, 0);
end
oo_.mom.model_moments = NaN(options_mom_.mom.mom_nbr,1);
model_moments = NaN(options_mom_.mom.mom_nbr,1);
for jm = 1:size(M_.matched_moments,1)
% First moments
if ~options_mom_.prefilter && (sum(M_.matched_moments{jm,3}) == 1)
idx1 = (oo_.mom.obs_var == find(oo_.dr.order_var==M_.matched_moments{jm,1}) );
oo_.mom.model_moments(jm,1) = pruned_state_space.E_y(idx1);
idx1 = (options_mom_.mom.obs_var == find(dr.order_var==M_.matched_moments{jm,1}) );
model_moments(jm,1) = pruned_state_space.E_y(idx1);
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
oo_.mom.model_moments_params_derivs(jm,:) = pruned_state_space.dE_y(idx1,:);
model_moments_params_derivs(jm,:) = pruned_state_space.dE_y(idx1,:);
end
end
% second moments
if (sum(M_.matched_moments{jm,3}) == 2)
idx1 = (oo_.mom.obs_var == find(oo_.dr.order_var==M_.matched_moments{jm,1}(1)) );
idx2 = (oo_.mom.obs_var == find(oo_.dr.order_var==M_.matched_moments{jm,1}(2)) );
idx1 = (options_mom_.mom.obs_var == find(dr.order_var==M_.matched_moments{jm,1}(1)) );
idx2 = (options_mom_.mom.obs_var == find(dr.order_var==M_.matched_moments{jm,1}(2)) );
if nnz(M_.matched_moments{jm,2}) == 0
% covariance
if options_mom_.prefilter
oo_.mom.model_moments(jm,1) = pruned_state_space.Var_y(idx1,idx2);
model_moments(jm,1) = pruned_state_space.Var_y(idx1,idx2);
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
oo_.mom.model_moments_params_derivs(jm,:) = pruned_state_space.dVar_y(idx1,idx2,:);
model_moments_params_derivs(jm,:) = pruned_state_space.dVar_y(idx1,idx2,:);
end
else
oo_.mom.model_moments(jm,1) = pruned_state_space.Var_y(idx1,idx2) + pruned_state_space.E_y(idx1)*pruned_state_space.E_y(idx2)';
model_moments(jm,1) = pruned_state_space.Var_y(idx1,idx2) + pruned_state_space.E_y(idx1)*pruned_state_space.E_y(idx2)';
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
for jp=1:totparam_nbr
oo_.mom.model_moments_params_derivs(jm,jp) = pruned_state_space.dVar_y(idx1,idx2,jp) + pruned_state_space.dE_y(idx1,jp)*pruned_state_space.E_y(idx2)' + pruned_state_space.E_y(idx1)*pruned_state_space.dE_y(idx2,jp)';
model_moments_params_derivs(jm,jp) = pruned_state_space.dVar_y(idx1,idx2,jp) + pruned_state_space.dE_y(idx1,jp)*pruned_state_space.E_y(idx2)' + pruned_state_space.E_y(idx1)*pruned_state_space.dE_y(idx2,jp)';
end
end
end
@ -189,15 +190,15 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
% autocovariance
lag = -M_.matched_moments{jm,2}(2); %note that leads/lags in M_.matched_moments are transformed such that first entry is always 0 and the second is a lag
if options_mom_.prefilter
oo_.mom.model_moments(jm,1) = pruned_state_space.Var_yi(idx1,idx2,lag);
model_moments(jm,1) = pruned_state_space.Var_yi(idx1,idx2,lag);
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
oo_.mom.model_moments_params_derivs(jm,:) = pruned_state_space.dVar_yi(idx1,idx2,lag,:);
model_moments_params_derivs(jm,:) = pruned_state_space.dVar_yi(idx1,idx2,lag,:);
end
else
oo_.mom.model_moments(jm,1) = pruned_state_space.Var_yi(idx1,idx2,lag) + pruned_state_space.E_y(idx1)*pruned_state_space.E_y(idx2)';
model_moments(jm,1) = pruned_state_space.Var_yi(idx1,idx2,lag) + pruned_state_space.E_y(idx1)*pruned_state_space.E_y(idx2)';
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
for jp=1:totparam_nbr
oo_.mom.model_moments_params_derivs(jm,jp) = vec( pruned_state_space.dVar_yi(idx1,idx2,lag,jp) + pruned_state_space.dE_y(idx1,jp)*pruned_state_space.E_y(idx2)' + pruned_state_space.E_y(idx1)*pruned_state_space.dE_y(idx2,jp)');
model_moments_params_derivs(jm,jp) = vec( pruned_state_space.dVar_yi(idx1,idx2,lag,jp) + pruned_state_space.dE_y(idx1,jp)*pruned_state_space.E_y(idx2)' + pruned_state_space.E_y(idx1)*pruned_state_space.dE_y(idx2,jp)');
end
end
end
@ -217,11 +218,11 @@ if strcmp(options_mom_.mom.mom_method,'SMM')
scaled_shock_series = zeros(size(options_mom_.mom.shock_series)); % initialize
scaled_shock_series(:,i_exo_var) = options_mom_.mom.shock_series(:,i_exo_var)*chol_S; % set non-zero entries
% simulate series
y_sim = simult_(M_, options_mom_, oo_.dr.ys, oo_.dr, scaled_shock_series, options_mom_.order);
y_sim = simult_(M_, options_mom_, dr.ys, dr, scaled_shock_series, options_mom_.order);
% provide meaningful penalty if data is nan or inf
if any(any(isnan(y_sim))) || any(any(isinf(y_sim)))
if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
fval = Inf(size(oo_.mom.Sw,1),1);
fval = Inf(size(weighting_info.Sw,1),1);
else
fval = Inf;
end
@ -229,12 +230,12 @@ if strcmp(options_mom_.mom.mom_method,'SMM')
info(4) = 0.1;
exit_flag = 0;
if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
fval = ones(size(oo_.mom.data_moments,1),1)*options_mom_.huge_number;
fval = ones(options_mom_.mom.mom_nbr,1)*options_mom_.huge_number;
end
return
end
% remove burn-in and focus on observables (note that y_sim is in declaration order)
y_sim = y_sim(oo_.dr.order_var(oo_.mom.obs_var) , end-options_mom_.mom.long+1:end)';
y_sim = y_sim(dr.order_var(options_mom_.mom.obs_var) , end-options_mom_.mom.long+1:end)';
if ~all(diag(M_.H)==0)
i_ME = setdiff(1:size(M_.H,1),find(diag(M_.H) == 0)); % find ME with 0 variance
chol_S = chol(M_.H(i_ME,i_ME)); % decompose rest
@ -246,27 +247,27 @@ if strcmp(options_mom_.mom.mom_method,'SMM')
if options_mom_.prefilter
y_sim = bsxfun(@minus, y_sim, mean(y_sim,1));
end
oo_.mom.model_moments = mom.get_data_moments(y_sim, oo_.mom.obs_var, oo_.dr.inv_order_var, M_.matched_moments, options_mom_);
model_moments = mom.get_data_moments(y_sim, options_mom_.mom.obs_var, dr.inv_order_var, M_.matched_moments, options_mom_);
end
%--------------------------------------------------------------------------
% Compute quadratic target function
%--------------------------------------------------------------------------
moments_difference = oo_.mom.data_moments - oo_.mom.model_moments;
moments_difference = data_moments - model_moments;
if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')
residuals = sqrt(options_mom_.mom.weighting_matrix_scaling_factor)*oo_.mom.Sw*moments_difference;
oo_.mom.Q = residuals'*residuals;
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
fval = residuals;
if options_mom_.mom.penalized_estimator
fval=[fval;(xparam-oo_.mom.prior.mean)./sqrt(diag(oo_.mom.prior.variance))];
fval=[fval;(xparam-bayestopt_.p1)./sqrt(diag(diag(bayestopt_.p2.^2)))];
end
else
fval = oo_.mom.Q;
fval = Q;
if options_mom_.mom.penalized_estimator
fval=fval+(xparam-oo_.mom.prior.mean)'/oo_.mom.prior.variance*(xparam-oo_.mom.prior.mean);
fval=fval+(xparam-bayestopt_.p1)'/(diag(bayestopt_.p2.^2))*(xparam-bayestopt_.p1);
end
end
if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_jacobian
@ -274,18 +275,18 @@ if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_meth
dxparam1 = eye(length(xparam));
end
for jp=1:length(xparam)
dmoments_difference = - oo_.mom.model_moments_params_derivs(:,jp);
dresiduals = sqrt(options_mom_.mom.weighting_matrix_scaling_factor)*oo_.mom.Sw*dmoments_difference;
dmoments_difference = - model_moments_params_derivs(:,jp);
dresiduals = sqrt(options_mom_.mom.weighting_matrix_scaling_factor)*weighting_info.Sw*dmoments_difference;
if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
if options_mom_.mom.penalized_estimator
df(:,jp)=[dresiduals;dxparam1(:,jp)./sqrt(diag(oo_.mom.prior.variance))];
df(:,jp)=[dresiduals;dxparam1(:,jp)./sqrt(diag(diag(bayestopt_.p2.^2)))];
else
df(:,jp) = dresiduals;
end
else
df(jp,1) = dresiduals'*residuals + residuals'*dresiduals;
if options_mom_.mom.penalized_estimator
df(jp,1)=df(jp,1)+(dxparam1(:,jp))'/oo_.mom.prior.variance*(xparam-oo_.mom.prior.mean)+(xparam-oo_.mom.prior.mean)'/oo_.mom.prior.variance*(dxparam1(:,jp));
df(jp,1)=df(jp,1)+(dxparam1(:,jp))'/(diag(bayestopt_.p2.^2))*(xparam-bayestopt_.p1)+(xparam-bayestopt_.p1)'/(diag(bayestopt_.p2.^2))*(dxparam1(:,jp));
end
end
end
@ -293,5 +294,4 @@ if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_meth
end
end % main function end
end % main function end

23
matlab/+mom/run.m
View File

@ -234,9 +234,9 @@ options_mom_.mom.compute_derivs = false; % flag to compute derivs in objective f
options_mom_.mom.vector_output = false; % specifies whether the objective function returns a vector
% decision rule
oo_.dr = set_state_space(oo_.dr,M_); % get state-space representation
oo_.mom.obs_var = []; % create index of observed variables in DR order
options_mom_.mom.obs_var = []; % create index of observed variables in DR order
for i = 1:options_mom_.obs_nbr
oo_.mom.obs_var = [oo_.mom.obs_var; find(strcmp(options_mom_.varobs{i}, M_.endo_names(oo_.dr.order_var)))];
options_mom_.mom.obs_var = [options_mom_.mom.obs_var; find(strcmp(options_mom_.varobs{i}, M_.endo_names(oo_.dr.order_var)))];
end
@ -255,7 +255,7 @@ if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_meth
% Get maximum lag number for autocovariances/autocorrelations
options_mom_.ar = max(cellfun(@max,M_.matched_moments(:,2))) - min(cellfun(@min,M_.matched_moments(:,2)));
% Check that only observed variables are involved in moments
not_observed_variables=setdiff(oo_.dr.inv_order_var([M_.matched_moments{:,1}]),oo_.mom.obs_var);
not_observed_variables=setdiff(oo_.dr.inv_order_var([M_.matched_moments{:,1}]),options_mom_.mom.obs_var);
if ~isempty(not_observed_variables)
skipline;
error('method_of_moments: You specified moments involving %s, but it is not a varobs!',M_.endo_names{oo_.dr.order_var(not_observed_variables)})
@ -419,7 +419,7 @@ if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_meth
% Provide info on data moments handling
fprintf('Computing data moments. Note that NaN values in the moments (due to leads and lags or missing data) are replaced by the mean of the corresponding moment.\n');
% Get data moments for the method of moments
[oo_.mom.data_moments, oo_.mom.m_data] = mom.get_data_moments(dataset_.data, oo_.mom.obs_var, oo_.dr.inv_order_var, M_.matched_moments, options_mom_);
[oo_.mom.data_moments, oo_.mom.m_data] = mom.get_data_moments(dataset_.data, options_mom_.mom.obs_var, oo_.dr.inv_order_var, M_.matched_moments, options_mom_);
if ~isreal(dataset_.data)
error('method_of_moments: The data moments contain complex values!')
end
@ -490,10 +490,10 @@ objective_function = str2func('mom.objective_function');
try
% Check for NaN or complex values of moment-distance-funtion evaluated at initial parameters
if strcmp(options_mom_.mom.mom_method,'SMM') || strcmp(options_mom_.mom.mom_method,'GMM')
oo_.mom.Sw = eye(options_mom_.mom.mom_nbr); % initialize with identity weighting matrix
weighting_info.Sw = eye(options_mom_.mom.mom_nbr); % initialize with identity weighting matrix
end
tic_id = tic;
[fval, info, ~, ~, ~, oo_, M_] = feval(objective_function, xparam0, BoundsInfo, oo_, estim_params_, M_, options_mom_);
[fval, info] = feval(objective_function, xparam0, oo_.mom.data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
elapsed_time = toc(tic_id);
if isnan(fval)
error('method_of_moments: The initial value of the objective function with identity weighting matrix is NaN!')
@ -535,16 +535,17 @@ mom.print_info_on_estimation_settings(options_mom_, number_of_estimated_paramete
% -------------------------------------------------------------------------
if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')
% compute mode
[xparam1, oo_, Woptflag] = mom.mode_compute_gmm_smm(xparam0, objective_function, oo_, M_, options_mom_, estim_params_, bayestopt_, BoundsInfo);
[xparam1, oo_.mom.weighting_info, oo_.mom.verbose] = mom.mode_compute_gmm_smm(xparam0, objective_function, oo_.mom.m_data, 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);
% compute standard errors at mode
options_mom_.mom.vector_output = false; % make sure flag is reset
M_ = set_all_parameters(xparam1,estim_params_,M_); % update M_ and oo_ (in particular to get oo_.mom.model_moments)
if strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_standard_errors
options_mom_.mom.compute_derivs = true; % for GMM we compute derivatives analytically in the objective function with this flag
end
[~, ~, ~,~,~, oo_] = feval(objective_function, xparam1, BoundsInfo, oo_, estim_params_, M_, options_mom_); % compute model moments and oo_.mom.model_moments_params_derivs
[~, ~, ~, ~, ~, oo_.mom.Q, oo_.mom.model_moments, oo_.mom.model_moments_params_derivs] = feval(objective_function, xparam1, 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);
options_mom_.mom.compute_derivs = false; % reset to not compute derivatives in objective function during optimization
[stdh,hessian_xparam1] = mom.standard_errors(xparam1, objective_function, BoundsInfo, oo_, estim_params_, M_, options_mom_, Woptflag);
[stdh, hessian_xparam1] = mom.standard_errors(xparam1, objective_function, oo_.mom.model_moments, oo_.mom.model_moments_params_derivs, oo_.mom.m_data, 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);
hessian_xparam1 = inv(hessian_xparam1); % mom.standard_errors returns the asymptotic covariance matrix
end
@ -553,7 +554,7 @@ end
% -------------------------------------------------------------------------
if options_mom_.mode_check.status
mode_check(objective_function, xparam1, hessian_xparam1, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, true,...
BoundsInfo, oo_, estim_params_, M_, options_mom_);
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
@ -564,7 +565,7 @@ if strcmp(options_mom_.mom.mom_method,'SMM') || strcmp(options_mom_.mom.mom_meth
% Store results in output structure
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));
% J test
oo_ = mom.Jtest(xparam1, objective_function, Woptflag, oo_, options_mom_, bayestopt_, BoundsInfo, estim_params_, M_, dataset_.nobs);
oo_.mom.J_test = mom.Jtest(xparam1, objective_function, oo_.mom.Q, oo_.mom.model_moments, oo_.mom.m_data, 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);
% display comparison of model moments and data moments
mom.display_comparison_moments(M_, options_mom_, oo_.mom.data_moments, oo_.mom.model_moments);
end

90
matlab/+mom/standard_errors.m
View File

@ -1,19 +1,28 @@
function [SE_values, Asympt_Var] = standard_errors(xparam, objective_function, BoundsInfo, oo_, estim_params_, M_, options_mom_, Wopt_flag)
% [SE_values, Asympt_Var] = standard_errors(xparam, objective_function, BoundsInfo, oo_, estim_params_, M_, options_mom_, Wopt_flag)
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)
% -------------------------------------------------------------------------
% This function computes standard errors to the method of moments estimates
% Adapted from replication codes of
% o Andreasen, Fernández-Villaverde, Rubio-Ramírez (2018): "The Pruned State-Space System for Non-Linear DSGE Models: Theory and Empirical Applications", Review of Economic Studies, 85(1):1-49.
% =========================================================================
% Adapted from replication codes of Andreasen, Fernández-Villaverde, Rubio-Ramírez (2018):
% "The Pruned State-Space System for Non-Linear DSGE Models: Theory and Empirical Applications",
% Review of Economic Studies, 85(1):1-49.
% -------------------------------------------------------------------------
% INPUTS
% o xparam: value of estimated parameters as returned by set_prior()
% o objective_function string of objective function
% o BoundsInfo: structure containing parameter bounds
% o oo_: structure for results
% o estim_params_: structure describing the estimated_parameters
% o M_ structure describing the model
% o options_mom_: structure information about all settings (specified by the user, preprocessor, and taken from global options_)
% o Wopt_flag: indicator whether the optimal weighting is actually used
% - xparam: [vector] value of estimated parameters as returned by set_prior()
% - objective_function [func] function handle with string of objective function
% - model_moments: [vector] model moments
% - model_moments_params_derivs: [matrix] analytical jacobian of the model moments wrt estimated parameters (currently for GMM only)
% - m_data [matrix] selected empirical moments at each point in time
% - data_moments: [vector] data with moments/irfs to match
% - weighting_info: [structure] storing information on weighting matrices
% - options_mom_: [structure] information about all settings (specified by the user, preprocessor, and taken from global options_)
% - M_ [structure] model information
% - estim_params_: [structure] information from estimated_params block
% - bayestopt_: [structure] information on the prior distributions
% - BoundsInfo: [structure] parameter bounds
% - dr: [structure] reduced form model
% - endo_steady_state: [vector] steady state value for endogenous variables (initval)
% - exo_steady_state: [vector] steady state value for exogenous variables (initval)
% - exo_det_steady_state: [vector] steady state value for exogenous deterministic variables (initval)
% -------------------------------------------------------------------------
% OUTPUTS
% o SE_values [nparam x 1] vector of standard errors
@ -26,9 +35,10 @@ function [SE_values, Asympt_Var] = standard_errors(xparam, objective_function, B
% o get_the_name
% o get_error_message
% o mom.objective_function
% o mom.optimal_weighting_matrix
% =========================================================================
% Copyright © 2020-2021 Dynare Team
% o mom.optimal_weighting_matrix
% -------------------------------------------------------------------------
% Copyright © 2020-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -44,21 +54,16 @@ function [SE_values, Asympt_Var] = standard_errors(xparam, objective_function, B
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------
% Author(s):
% o Willi Mutschler (willi@mutschler.eu)
% o Johannes Pfeifer (johannes.pfeifer@unibw.de)
% =========================================================================
% Some dimensions
num_mom = size(oo_.mom.model_moments,1);
num_mom = size(model_moments,1);
dim_params = size(xparam,1);
D = zeros(num_mom,dim_params);
eps_value = options_mom_.mom.se_tolx;
if strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_standard_errors
fprintf('\nComputing standard errors using analytical derivatives of moments\n');
D = oo_.mom.model_moments_params_derivs; %already computed in objective function via get_perturbation_params.m
D = model_moments_params_derivs; % already computed in objective function via get_perturbation_params.m
idx_nan = find(any(isnan(D)));
if any(idx_nan)
for i = idx_nan
@ -72,19 +77,17 @@ if strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_standa
else
fprintf('\nComputing standard errors using numerical derivatives of moments\n');
for i=1:dim_params
%Positive step
% positive step
xparam_eps_p = xparam;
xparam_eps_p(i,1) = xparam_eps_p(i) + eps_value;
[~, info_p, ~, ~,~, oo__p] = feval(objective_function, xparam_eps_p, BoundsInfo, oo_, estim_params_, M_, options_mom_);
% Negative step
xparam_eps_p(i,1) = xparam_eps_p(i) + eps_value;
[~, info_p, ~, ~, ~, ~, model_moments_p] = feval(objective_function, xparam_eps_p, data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
% negative step
xparam_eps_m = xparam;
xparam_eps_m(i,1) = xparam_eps_m(i) - eps_value;
[~, info_m, ~, ~,~, oo__m] = feval(objective_function, xparam_eps_m, BoundsInfo, oo_, estim_params_, M_, options_mom_);
% The Jacobian:
xparam_eps_m(i,1) = xparam_eps_m(i) - eps_value;
[~, info_m, ~, ~, ~, ~, model_moments_m] = feval(objective_function, xparam_eps_m, data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
% the Jacobian
if nnz(info_p)==0 && nnz(info_m)==0
D(:,i) = (oo__p.mom.model_moments - oo__m.mom.model_moments)/(2*eps_value);
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);
if info_p(1)==42
@ -95,7 +98,7 @@ else
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)
else
message_m = get_error_message(info_m, options_mom_);
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)
@ -106,22 +109,19 @@ else
end
end
end
T = options_mom_.nobs; %Number of observations
T = options_mom_.nobs;
if isfield(options_mom_,'variance_correction_factor')
T = T*options_mom_.variance_correction_factor;
end
WW = oo_.mom.Sw'*oo_.mom.Sw;
if Wopt_flag
% We have the optimal weighting matrix
Asympt_Var = 1/T*((D'*WW*D)\eye(dim_params));
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));
else
% We do not have the optimal weighting matrix yet
WWopt = mom.optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag);
% 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;
end
SE_values = sqrt(diag(Asympt_Var));
SE_values = sqrt(diag(Asympt_Var));