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

MoM: Add analytical standard errors for GMM 

Reset analytic_derivation_mode for steadystate file parameter changes


MoM: Fix GMM analytical standard errors wrong dimensions in autocovs
time-shift
Willi Mutschler 2020-12-16 12:03:21 +01:00
parent ec8ea32b3e
commit eae5e2f029
No known key found for this signature in database
GPG Key ID: 91E724BF17A73F6D
4 changed files with 107 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
matlab/method_of_moments/method_of_moments.m
View File

@ -169,7 +169,9 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
fprintf('GMM at higher order only works with pruning, so we set pruning option to 1.\n');
options_mom_.pruning = true;
end
options_mom_.mom = set_default_option(options_mom_.mom,'analytic_standard_errors',false); % compute standard errors numerically (0) or analytically (1). Analytical derivatives are only available for GMM.
end
options_mom_.mom.compute_derivs = false;% flag to compute derivs in objective function (needed for analytic standard errors with GMM)
% General options that can be set by the user in the mod file, otherwise default values are provided
@ -326,6 +328,7 @@ options_mom_.solveopt = options_.solveopt;
options_mom_.gradient_method = options_.gradient_method;
options_mom_.gradient_epsilon = options_.gradient_epsilon;
options_mom_.analytic_derivation = 0;
options_mom_.analytic_derivation_mode = 0; % needed by get_perturbation_params_derivs.m, ie use efficient sylvester equation method to compute analytical derivatives as in Ratto & Iskrev (2012)
options_mom_.vector_output= false; % specifies whether the objective function returns a vector
@ -688,6 +691,11 @@ if isfield(estim_params_,'param_vals') && ~isempty(estim_params_.param_vals)
fprintf('This will override parameter values and may lead to wrong results.\n')
fprintf('Check whether this is really intended.\n')
warning('The steady state file internally changes the values of the estimated parameters.')
if strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_standard_errors
fprintf('For analytical standard errors, the parameter-Jacobians of the dynamic model and of the steady-state will be computed numerically\n'),
fprintf('(re-set options_mom_.analytic_derivation_mode= -2)'),
options_mom_.analytic_derivation_mode= -2;
end
end
end
@ -785,8 +793,14 @@ fprintf('\n - perturbation order: %d', options_mom_.order)
if options_mom_.order > 1 && options_mom_.pruning
fprintf(' (with pruning)')
end
if strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_standard_errors
fprintf('\n - standard errors: analytic derivatives');
else
fprintf('\n - standard errors: numerical derivatives');
end
fprintf('\n - number of matched moments: %d', options_mom_.mom.mom_nbr);
fprintf('\n - number of parameters: %d\n\n', length(xparam0));
fprintf('\n - number of parameters: %d', length(xparam0));
fprintf('\n\n');
% -------------------------------------------------------------------------
% Step 7b: Iterated method of moments estimation
@ -861,8 +875,12 @@ for stage_iter=1:size(options_mom_.mom.weighting_matrix,1)
options_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') && 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
[fval, ~, ~,~,~, oo_] = feval(objective_function, xparam1, Bounds, oo_, estim_params_, M_, options_mom_);
% Compute Standard errors
options_mom_.mom.compute_derivs = false; % reset to not compute derivatives in objective function during optimization
SE = method_of_moments_standard_errors(xparam1, objective_function, Bounds, oo_, estim_params_, M_, options_mom_, Woptflag);
% Store results in output structure

55
matlab/method_of_moments/method_of_moments_objective_function.m
View File

@ -109,7 +109,32 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
%--------------------------------------------------------------------------
% 3. Set up pruned state-space system and compute model moments
%--------------------------------------------------------------------------
pruned_state_space = pruned_state_space_system(M_, options_mom_, dr, oo_.dr.obs_var, options_mom_.ar, 0, 0);
if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_standard_errors
indpmodel = []; %initialize index for model parameters
if ~isempty(estim_params_.param_vals)
indpmodel = estim_params_.param_vals(:,1); %values correspond to parameters declaration order, row number corresponds to order in estimated_params
end
indpstderr=[]; %initialize index for stderr parameters
if ~isempty(estim_params_.var_exo)
indpstderr = estim_params_.var_exo(:,1); %values correspond to varexo declaration order, row number corresponds to order in estimated_params
end
indpcorr=[]; %initialize matrix for corr paramters
if ~isempty(estim_params_.corrx)
indpcorr = estim_params_.corrx(:,1:2); %values correspond to varexo declaration order, row number corresponds to order in estimated_params
end
if estim_params_.nvn || estim_params_.ncn %nvn is number of stderr parameters and ncn is number of corr parameters of measurement innovations as declared in estimated_params
error('Analytic computation of standard errrors does not (yet) support measurement errors.\nInstead, define them explicitly as varexo and provide measurement equations in the model definition.\nAlternatively, use numerical standard errors.')
end
modparam_nbr = estim_params_.np; % number of model parameters as declared in estimated_params
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;
dr.derivs = get_perturbation_params_derivs(M_, options_mom_, estim_params_, oo_, 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_state_space_system(M_, options_mom_, dr, oo_.dr.obs_var, options_mom_.ar, 0, 1);
else
pruned_state_space = pruned_state_space_system(M_, options_mom_, dr, oo_.dr.obs_var, options_mom_.ar, 0, 0);
end
oo_.mom.model_moments = NaN(options_mom_.mom.mom_nbr,1);
offset = 0;
@ -118,6 +143,9 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
E_y = pruned_state_space.E_y;
E_y_nbr = nnz(options_mom_.mom.index.E_y);
oo_.mom.model_moments(offset+1:E_y_nbr,1) = E_y(options_mom_.mom.index.E_y);
if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_standard_errors
oo_.mom.model_moments_params_derivs(offset+1:E_y_nbr,:) = pruned_state_space.dE_y(options_mom_.mom.index.E_y,:);
end
offset = offset + E_y_nbr;
end
% Second moments
@ -125,22 +153,47 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
if isfield(options_mom_.mom.index,'E_yy') && nnz(options_mom_.mom.index.E_yy) > 0
if options_mom_.prefilter
E_yy = pruned_state_space.Var_y;
if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_standard_errors
dE_yy = pruned_state_space.dVar_y;
end
else
E_yy = pruned_state_space.Var_y + pruned_state_space.E_y*pruned_state_space.E_y';
if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_standard_errors
dE_yy = pruned_state_space.dVar_y;
for jp=1:totparam_nbr
dE_yy(:,:,jp) = dE_yy(:,:,jp) + pruned_state_space.dE_y(:,jp)*pruned_state_space.E_y' + pruned_state_space.E_y*pruned_state_space.dE_y(:,jp)';
end
end
end
E_yy_nbr = nnz(tril(options_mom_.mom.index.E_yy));
oo_.mom.model_moments(offset+(1:E_yy_nbr),1) = E_yy(tril(options_mom_.mom.index.E_yy));
if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_standard_errors
oo_.mom.model_moments_params_derivs(offset+(1:E_yy_nbr),:) = reshape(dE_yy(repmat(tril(options_mom_.mom.index.E_yy),[1 1 totparam_nbr])),E_yy_nbr,totparam_nbr);
end
offset = offset + E_yy_nbr;
end
% Lead/lags covariance
if isfield(options_mom_.mom.index,'E_yyt') && nnz(options_mom_.mom.index.E_yyt) > 0
if options_mom_.prefilter
E_yyt = pruned_state_space.Var_yi;
if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_standard_errors
dE_yyt = pruned_state_space.dVar_yi;
end
else
E_yyt = pruned_state_space.Var_yi + repmat(pruned_state_space.E_y*pruned_state_space.E_y',[1 1 size(pruned_state_space.Var_yi,3)]);
if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_standard_errors
dE_yyt = pruned_state_space.dVar_yi;
for jp=1:totparam_nbr
dE_yyt(:,:,:,jp) = dE_yyt(:,:,:,jp) + repmat(pruned_state_space.dE_y(:,jp)*pruned_state_space.E_y',[1 1 size(pruned_state_space.Var_yi,3)])...
+ repmat(pruned_state_space.E_y*pruned_state_space.dE_y(:,jp)',[1 1 size(pruned_state_space.Var_yi,3)]);
end
end
end
E_yyt_nbr = nnz(options_mom_.mom.index.E_yyt);
oo_.mom.model_moments(offset+(1:E_yyt_nbr),1) = E_yyt(options_mom_.mom.index.E_yyt);
if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_standard_errors
oo_.mom.model_moments_params_derivs(offset+(1:E_yyt_nbr),:) = reshape(dE_yyt(repmat(options_mom_.mom.index.E_yyt,[1 1 1 totparam_nbr])),E_yyt_nbr,totparam_nbr);
end
end
elseif strcmp(options_mom_.mom.mom_method,'SMM')

50
matlab/method_of_moments/method_of_moments_standard_errors.m
View File

@ -57,29 +57,35 @@ dim_params = size(xparam,1);
D = zeros(num_mom,dim_params);
eps_value = options_mom_.mom.se_tolx;
for i=1:dim_params
%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, Bounds, oo_, estim_params_, M_, options_mom_);
% 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, Bounds, oo_, estim_params_, M_, options_mom_);
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
else
fprintf('\nComputing standard errors using numerical derivatives of moments\n');
for i=1:dim_params
%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, Bounds, oo_, estim_params_, M_, options_mom_);
% 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);
else
problpar = get_the_name(i,options_mom_.TeX, M_, estim_params_, options_mom_);
message_p = get_error_message(info_p, options_mom_);
message_m = get_error_message(info_m, options_mom_);
warning('method_of_moments:info','Cannot compute the Jacobian for parameter %s - no standard errors available\n %s %s\nCheck your bounds and/or priors, or use a different optimizer.\n',problpar, message_p, message_m)
Asympt_Var = NaN(length(xparam),length(xparam));
SE_values = NaN(length(xparam),1);
return
% 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, Bounds, oo_, estim_params_, M_, options_mom_);
% 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);
else
problpar = get_the_name(i,options_mom_.TeX, M_, estim_params_, options_mom_);
message_p = get_error_message(info_p, options_mom_);
message_m = get_error_message(info_m, options_mom_);
warning('method_of_moments:info','Cannot compute the Jacobian for parameter %s - no standard errors available\n %s %s\nCheck your bounds and/or priors, or use a different optimizer.\n',problpar, message_p, message_m)
Asympt_Var = NaN(length(xparam),length(xparam));
SE_values = NaN(length(xparam),1);
return
end
end
end

7
tests/estimation/method_of_moments/AnScho_MoM.mod
View File

@ -25,7 +25,7 @@
@#define estimParams = 1
% Note that we set the numerical optimization tolerance levels very large to speed up the testsuite
@#define optimizer = 4
@#define optimizer = 13
var c p R g y z INFL INT YGR;
varexo e_r e_g e_z;
@ -248,7 +248,7 @@ end
%, optim = ('TolFun', 1e-5
% ,'TolX', 1e-6
% ) % a list of NAME and VALUE pairs to set options for the optimization routines. Available options depend on mode_compute
%, silent_optimizer % run minimization of moments distance silently without displaying results or saving files in between
, silent_optimizer % run minimization of moments distance silently without displaying results or saving files in between
% , tolf = 1e-5 % convergence criterion on function value for numerical differentiation
% , tolx = 1e-6 % convergence criterion on funciton input for numerical differentiation
%
@ -267,6 +267,9 @@ end
% , sylvester_fixed_point_tol = 1e-12 % convergence criterion used in the fixed point Sylvester solver
% , qz_criterium = 0.999999 % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems [IS THIS CORRET @wmutschl]
% , qz_zero_threshold = 1e-6 % value used to test if a generalized eigenvalue is 0/0 in the generalized Schur decomposition
@#if mommethod == "GMM"
, analytic_standard_errors
@#endif
);
@#endfor