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

method_of_moments: cosmetical changes to mom.objective function

kalman-mex
Willi Mutschler 2023-09-04 16:24:32 +02:00
parent 0cd65df72a
commit c0cae0ebaa
No known key found for this signature in database
GPG Key ID: 91E724BF17A73F6D
1 changed files with 115 additions and 118 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

141
matlab/+mom/objective_function.m
View File

@ -1,42 +1,43 @@
function [fval, info, exit_flag, df, junk1, oo_, M_, options_mom_] = objective_function(xparam1, Bounds, oo_, estim_params_, M_, options_mom_)
% [fval, info, exit_flag, df, junk1, oo_, M_, options_mom_] = objective_function(xparam1, Bounds, oo_, estim_params_, M_, options_mom_)
function [fval, info, exit_flag, df, junkHessian, oo_, M_] = objective_function(xparam, Bounds, oo_, estim_params_, M_, options_mom_, bayestopt_)
% [fval, info, exit_flag, df, junk1, oo_, M_] = objective_function(xparam, Bounds, oo_, estim_params_, M_, options_mom_, bayestopt_)
% -------------------------------------------------------------------------
% This function evaluates the objective function for GMM/SMM estimation
% This function evaluates the objective function for method of moments estimation
% =========================================================================
% INPUTS
% o xparam1: current value of estimated parameters as returned by set_prior()
% o Bounds: 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 xparam: [vector] current value of estimated parameters as returned by set_prior()
% o Bounds: [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 bayestopt_: [structure] information about the prior
% -------------------------------------------------------------------------
% OUTPUTS
% o fval: value of the quadratic form of the moment difference (except for lsqnonlin, where this is done implicitly)
% o info: vector storing error code and penalty
% o exit_flag: 0 if error, 1 if no error
% o df: analytical parameter Jacobian of the quadratic form of the moment difference (for GMM only)
% o junk1: empty matrix required for optimizer interface (Hessian would go here)
% o oo_: structure containing the results with the following updated fields:
% - mom.model_moments [numMom x 1] vector with model moments
% - mom.Q value of the quadratic form of the moment difference
% - mom.model_moments_params_derivs
% [numMom x numParams] Jacobian matrix of derivatives of model_moments with respect to estimated parameters
% (only for GMM with analytical derivatives)
% o M_: Matlab's structure describing the model
% o options_mom_: structure information about all settings (specified by the user, preprocessor, and taken from global options_)
% 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
% -------------------------------------------------------------------------
% This function is called by
% o mom.run.m
% o dynare_minimize_objective.m
% o mom.run
% o dynare_minimize_objective
% -------------------------------------------------------------------------
% This function calls
% o check_bounds_and_definiteness_estimation
% o get_perturbation_params_derivs
% o mom.get_data_moments
% o pruned_state_space_system
% o resol
% o set_all_parameters
% o simult_
% =========================================================================
% Copyright © 2020-2021 Dynare Team
% Copyright © 2020-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -52,41 +53,41 @@ function [fval, info, exit_flag, df, junk1, oo_, M_, options_mom_] = objective_f
%
% 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 (jpfeifer@uni-koeln.de)
% =========================================================================
%% TO DO
% check the info values and make use of meaningful penalties
% how do we do the penalty for the prior??
%------------------------------------------------------------------------------
% 0. Initialization of the returned variables and others...
% Initialization of the returned variables and others...
%------------------------------------------------------------------------------
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(xparam1),length(xparam1));
df = nan(size(oo_.mom.data_moments,1)+length(xparam),length(xparam));
else
df = nan(size(oo_.mom.data_moments,1),length(xparam1));
df = nan(size(oo_.mom.data_moments,1),length(xparam));
end
else
df = nan(1,length(xparam1));
df = nan(1,length(xparam));
end
else
df=[]; %required to be empty by e.g. newrat
end
junk1 = [];
junk2 = [];
end
%--------------------------------------------------------------------------
% 1. Get the structural parameters & define penalties
% Get the structural parameters and define penalties
%--------------------------------------------------------------------------
% Ensure that xparam1 is a column vector; particleswarm.m requires this.
xparam1 = xparam1(:);
M_ = set_all_parameters(xparam1, estim_params_, M_);
[fval,info,exit_flag]=check_bounds_and_definiteness_estimation(xparam1, M_, estim_params_, Bounds);
xparam = xparam(:);
M_ = set_all_parameters(xparam, estim_params_, M_);
[fval,info,exit_flag] = check_bounds_and_definiteness_estimation(xparam, M_, estim_params_, Bounds);
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;
@ -94,13 +95,13 @@ if info(1)
return
end
%--------------------------------------------------------------------------
% 2. call resol to compute steady state and model solution
% Call resol to compute steady state and model solution
%--------------------------------------------------------------------------
% Compute linear approximation around the deterministic steady state
[dr, info, M_, oo_] = resol(0, M_, options_mom_, oo_);
% 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 ||...
@ -125,10 +126,11 @@ if info(1)
end
end
%--------------------------------------------------------------------------
% GMM: Set up pruned state-space system and compute model moments
%--------------------------------------------------------------------------
if strcmp(options_mom_.mom.mom_method,'GMM')
%--------------------------------------------------------------------------
% 3. Set up pruned state-space system and compute model moments
%--------------------------------------------------------------------------
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
indpmodel = []; % initialize index for model parameters
if ~isempty(estim_params_.param_vals)
@ -155,7 +157,6 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
else
pruned_state_space = pruned_state_space_system(M_, options_mom_, dr, oo_.mom.obs_var, options_mom_.ar, 0, 0);
end
oo_.mom.model_moments = NaN(options_mom_.mom.mom_nbr,1);
for jm = 1:size(M_.matched_moments,1)
% First moments
@ -166,12 +167,12 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
oo_.mom.model_moments_params_derivs(jm,:) = pruned_state_space.dE_y(idx1,:);
end
end
% Second moments
% 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)) );
if nnz(M_.matched_moments{jm,2}) == 0
% Covariance
% covariance
if options_mom_.prefilter
oo_.mom.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 )
@ -186,8 +187,8 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
end
end
else
% Autocovariance
lag = -M_.matched_moments{jm,2}(2); %note that leads/lags in matched_moments are transformed such that first entry is always 0 and the second is a lag
% 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);
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
@ -204,19 +205,18 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
end
end
end
end
elseif strcmp(options_mom_.mom.mom_method,'SMM')
%------------------------------------------------------------------------------
% 3. Compute Moments of the model solution for normal innovations
% SMM: Compute Moments of the model solution for Gaussian innovations
%------------------------------------------------------------------------------
if strcmp(options_mom_.mom.mom_method,'SMM')
% create shock series with correct covariance matrix from iid standard normal shocks
i_exo_var = setdiff(1:M_.exo_nbr, find(diag(M_.Sigma_e) == 0 )); % find singular entries in covariance
chol_S = chol(M_.Sigma_e(i_exo_var,i_exo_var));
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_, dr.ys, dr, scaled_shock_series, options_mom_.order);
% provide meaningful penalty if data is nan or inf
@ -234,10 +234,8 @@ elseif strcmp(options_mom_.mom.mom_method,'SMM')
end
return
end
% Remove burn-in and focus on observables (note that y_sim is in declaration order)
% 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)';
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
@ -245,42 +243,40 @@ elseif strcmp(options_mom_.mom.mom_method,'SMM')
shock_mat(:,i_ME)=options_mom_.mom.ME_shock_series(:,i_ME)*chol_S;
y_sim = y_sim+shock_mat;
end
% Remove mean if centered moments
% remove mean if centered moments
if options_mom_.prefilter
y_sim = bsxfun(@minus, y_sim, mean(y_sim,1));
end
oo_.mom.model_moments = mom.data_moments(y_sim, oo_, M_.matched_moments, options_mom_);
end
%--------------------------------------------------------------------------
% 4. Compute quadratic target function
% Compute quadratic target function
%--------------------------------------------------------------------------
moments_difference = oo_.mom.data_moments - oo_.mom.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;
if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
fval = residuals;
if options_mom_.mom.penalized_estimator
fval=[fval;(xparam1-oo_.mom.prior.mean)./sqrt(diag(oo_.mom.prior.variance))];
fval=[fval;(xparam-oo_.mom.prior.mean)./sqrt(diag(oo_.mom.prior.variance))];
end
else
fval = oo_.mom.Q;
if options_mom_.mom.penalized_estimator
fval=fval+(xparam1-oo_.mom.prior.mean)'/oo_.mom.prior.variance*(xparam1-oo_.mom.prior.mean);
fval=fval+(xparam-oo_.mom.prior.mean)'/oo_.mom.prior.variance*(xparam-oo_.mom.prior.mean);
end
end
if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_jacobian
if options_mom_.mom.penalized_estimator
dxparam1 = eye(length(xparam1));
dxparam1 = eye(length(xparam));
end
for jp=1:length(xparam1)
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;
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))];
@ -290,7 +286,8 @@ if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_jacobian
else
df(:,jp) = dresiduals'*residuals + residuals'*dresiduals;
if options_mom_.mom.penalized_estimator
df(:,jp)=df(:,jp)+(dxparam1(:,jp))'/oo_.mom.prior.variance*(xparam1-oo_.mom.prior.mean)+(xparam1-oo_.mom.prior.mean)'/oo_.mom.prior.variance*(dxparam1(:,jp));
df(:,jp)=df(:,jp)+(dxparam1(:,jp))'/oo_.mom.prior.variance*(xparam-oo_.mom.prior.mean)+(xparam-oo_.mom.prior.mean)'/oo_.mom.prior.variance*(dxparam1(:,jp));
end
end
end
end