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

method_of_moments: refactor iterated GMM/SMM estimation

kalman-mex
Willi Mutschler 2023-09-04 20:21:29 +02:00
parent 9dc4878f9f
commit 07b62fe554
No known key found for this signature in database
GPG Key ID: 91E724BF17A73F6D
2 changed files with 150 additions and 92 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

129
matlab/+mom/mode_compute_gmm_smm.m Normal file
View File

@ -0,0 +1,129 @@
function [xparam1, oo_, Woptflag] = mode_compute_gmm_smm(xparam0, objective_function, oo_, M_, options_mom_, estim_params_, bayestopt_, Bounds)
% function [xparam1, oo_, Woptflag] = mode_compute_gmm_smm(xparam0, objective_function, oo_, M_, options_mom_, estim_params_, bayestopt_, Bounds)
% -------------------------------------------------------------------------
% 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
% Bounds: [structure] bounds for optimization
% -------------------------------------------------------------------------
% OUTPUT
% xparam1: [vector] mode of objective function
% oo_: [structure] updated results
% Woptflag: [logical] true if optimal weighting matrix was computed
% -------------------------------------------------------------------------
% This function is called by
% o mom.run
% -------------------------------------------------------------------------
% This function calls
% o mom.optimal_weighting_matrix
% o mom.display_estimation_results_table
% o dynare_minimize_objective
% o mom.objective_function
% =========================================================================
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
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;
switch lower(options_mom_.mom.weighting_matrix{stage_iter})
case 'identity_matrix'
fprintf(' - identity weighting matrix\n');
weighting_matrix = eye(options_mom_.mom.mom_nbr);
case 'diagonal'
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) ));
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) ));
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);
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;
end
otherwise % user specified matrix in file
fprintf(' - user-specified weighting matrix\n');
try
load(options_mom_.mom.weighting_matrix{stage_iter},'weighting_matrix')
catch
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!'])
end
end
try % check for positive definiteness of weighting_matrix
oo_.mom.Sw = chol(weighting_matrix);
catch
error('method_of_moments: Specified ''weighting_matrix'' is not positive definite. Check whether your model implies stochastic singularity!')
end
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, Bounds, oo_, estim_params_, M_, options_mom_);
else
if options_mom_.optimizer_vec{optim_iter} == 13
options_mom_.mom.vector_output = true;
else
options_mom_.mom.vector_output = false;
end
if strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_jacobian && ismember(options_mom_.optimizer_vec{optim_iter},options_mom_.mom.analytic_jacobian_optimizers) %do this only for gradient-based optimizers
options_mom_.mom.compute_derivs = true;
else
options_mom_.mom.compute_derivs = false;
end
[xparam1, fval, exitflag] = dynare_minimize_objective(objective_function, xparam0, options_mom_.optimizer_vec{optim_iter}, options_mom_, [Bounds.lb Bounds.ub], bayestopt_.name, bayestopt_, [],...
Bounds, oo_, estim_params_, M_, options_mom_);
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));
end
xparam0 = xparam1;
end
options_mom_.vector_output = false;
[~, ~, ~,~,~, oo_] = feval(objective_function, xparam1, Bounds, oo_, estim_params_, M_, options_mom_); % get oo_.mom.model_moments for iterated GMM/SMM to compute optimal weighting matrix
end

113
matlab/+mom/run.m
View File

@ -547,100 +547,28 @@ mom.print_info_on_estimation_settings(options_mom_, number_of_estimated_paramete
% -------------------------------------------------------------------------
% Step 7b: Iterated method of moments estimation
% GMM/SMM: iterated estimation
% -------------------------------------------------------------------------
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;
switch lower(options_mom_.mom.weighting_matrix{stage_iter})
case 'identity_matrix'
fprintf(' - identity weighting matrix\n');
weighting_matrix = eye(options_mom_.mom.mom_nbr);
case 'diagonal'
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) ));
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) ));
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);
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;
end
otherwise %user specified matrix in file
fprintf(' - user-specified weighting matrix\n');
try
load(options_mom_.mom.weighting_matrix{stage_iter},'weighting_matrix')
catch
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'])
end
end
try %check for positive definiteness of weighting_matrix
oo_.mom.Sw = chol(weighting_matrix);
catch
error('method_of_moments: Specified weighting_matrix is not positive definite. Check whether your model implies stochastic singularity.')
end
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, Bounds, oo_, estim_params_, M_, options_mom_);
else
if options_mom_.optimizer_vec{optim_iter}==13
options_mom_.mom.vector_output = true;
else
options_mom_.mom.vector_output = false;
end
if strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_jacobian && ismember(options_mom_.optimizer_vec{optim_iter},options_mom_.mom.analytic_jacobian_optimizers) %do this only for gradient-based optimizers
options_mom_.mom.compute_derivs = true;
else
options_mom_.mom.compute_derivs = false;
end
[xparam1, fval, exitflag] = dynare_minimize_objective(objective_function, xparam0, options_mom_.optimizer_vec{optim_iter}, options_mom_, [Bounds.lb Bounds.ub], bayestopt_.name, bayestopt_, [],...
Bounds, oo_, estim_params_, M_, options_mom_);
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));
end
xparam0=xparam1;
end
options_mom_.mom.vector_output = false;
% Update M_ and DynareResults (in particular to get oo_.mom.model_moments)
M_ = set_all_parameters(xparam1,estim_params_,M_);
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_, Bounds);
% 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 DynareResults (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
options_mom_.mom.compute_derivs = true; % for GMM we compute derivatives analytically in the objective function with this flag
end
[fval, ~, ~,~,~, oo_] = feval(objective_function, xparam1, Bounds, oo_, estim_params_, M_, options_mom_);
[~, ~, ~,~,~, oo_] = feval(objective_function, xparam1, Bounds, oo_, estim_params_, M_, options_mom_); % compute model moments and oo_.mom.model_moments_params_derivs
options_mom_.mom.compute_derivs = false; % reset to not compute derivatives in objective function during optimization
SE = mom.standard_errors(xparam1, objective_function, Bounds, oo_, estim_params_, M_, options_mom_, Woptflag);
% Store results in output structure
oo_.mom = display_estimation_results_table(xparam1,SE,M_,options_mom_,estim_params_,bayestopt_,oo_.mom,prior_dist_names,sprintf('%s (STAGE %u)',options_mom_.mom.mom_method,stage_iter),sprintf('%s_stage_%u',lower(options_mom_.mom.mom_method),stage_iter));
% mode check plots
if options_mom_.mode_check.status
mom.check_plot(objective_function,xparam1,SE,options_mom_,M_,estim_params_,Bounds,bayestopt_,...
Bounds, oo_, estim_params_, M_, options_mom_);
end
end
% -------------------------------------------------------------------------
% Step 8: J test
% -------------------------------------------------------------------------
@ -666,9 +594,14 @@ if options_mom_.mom.mom_nbr > length(xparam1)
fprintf('p-value of J-test statistic: %f\n',oo_.mom.J_test.p_val)
end
% -------------------------------------------------------------------------
% Step 9: Display estimation results
% display final estimation results
% -------------------------------------------------------------------------
if strcmp(options_mom_.mom.mom_method,'SMM') || strcmp(options_mom_.mom.mom_method,'GMM')
% Store results in output structure
oo_.mom = display_estimation_results_table(xparam1,SE,M_,options_mom_,estim_params_,bayestopt_,oo_.mom,prior_dist_names,options_mom_.mom.mom_method,lower(options_mom_.mom.mom_method));
title = ['Comparison of data moments and model moments (',options_mom_.mom.mom_method,')'];
headers = {'Moment','Data','Model'};
for jm = 1:size(M_.matched_moments,1)
@ -703,10 +636,6 @@ dyntable(options_mom_, title, headers, labels, data_mat, cellofchararraymaxlengt
if options_mom_.TeX
dyn_latex_table(M_, options_mom_, title, ['comparison_moments_', options_mom_.mom.mom_method], headers, labels_TeX, data_mat, cellofchararraymaxlength(labels)+2, 10, 7);
end
if options_mom_.mode_check.status
mom.check_plot(objective_function,xparam1,SE,options_mom_,M_,estim_params_,Bounds,bayestopt_,...
Bounds, oo_, estim_params_, M_, options_mom_)
end
fprintf('\n==== Method of Moments Estimation (%s) Completed ====\n\n',options_mom_.mom.mom_method)