diff --git a/matlab/+mom/Jtest.m b/matlab/+mom/Jtest.m
index 5c5001a3dbb0604c7f20508f7e25f476fab1a740..d455b3754b0c2fdb9949a5d52bce22bdba1b22df 100644
--- a/matlab/+mom/Jtest.m
+++ b/matlab/+mom/Jtest.m
@@ -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
\ No newline at end of file
diff --git a/matlab/+mom/mode_compute_gmm_smm.m b/matlab/+mom/mode_compute_gmm_smm.m
index 5514b0bd6a25bca52a8f1d824e3c3c2efa6173a8..c4d28f53844b91239b3637c4cf14d87492c1a6db 100644
--- a/matlab/+mom/mode_compute_gmm_smm.m
+++ b/matlab/+mom/mode_compute_gmm_smm.m
@@ -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
diff --git a/matlab/+mom/objective_function.m b/matlab/+mom/objective_function.m
index c6b06948ebe7cfb4336bbbcd656091af0390f88f..2514a295603f851fd00b568549c7371e5d117809 100644
--- a/matlab/+mom/objective_function.m
+++ b/matlab/+mom/objective_function.m
@@ -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
\ No newline at end of file
diff --git a/matlab/+mom/run.m b/matlab/+mom/run.m
index 57464756a10c40cae337482a123fe978c3faf01e..3109e7d720ee695aa700deed0b1eae5afa30ddbd 100644
--- a/matlab/+mom/run.m
+++ b/matlab/+mom/run.m
@@ -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
diff --git a/matlab/+mom/standard_errors.m b/matlab/+mom/standard_errors.m
index 0fac91d94133cd3f3937a4df0bfc68fbf907c253..1eeed07c62be002dfb7a46935e3e91c00e364812 100644
--- a/matlab/+mom/standard_errors.m
+++ b/matlab/+mom/standard_errors.m
@@ -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));
\ No newline at end of file