Commit b273a279 authored by Willi Mutschler's avatar Willi Mutschler Committed by Johannes Pfeifer

Fix iterated method of moments

Also improve fprintf descriptions and comments and other cosmetical changes
parent 9c99eac5
......@@ -14,7 +14,7 @@ function [dataMoments, m_data] = method_of_moments_data_moments(data, oo_, match
% -------------------------------------------------------------------------
% This function is called by
% o method_of_moments.m
% o method_of_moments_SMM.m
% o method_of_moments_objective_function.m
% =========================================================================
% Copyright (C) 2020 Dynare Team
%
......@@ -39,7 +39,7 @@ function [dataMoments, m_data] = method_of_moments_data_moments(data, oo_, match
% =========================================================================
% Initialization
T = size(data,1); % Number of observations (T) and number of observables (ny)
T = size(data,1); % Number of observations (T)
dataMoments = NaN(options_mom_.mom.mom_nbr,1);
m_data = NaN(T,options_mom_.mom.mom_nbr);
% Product moment for each time period, i.e. each row t contains y_t1(l1)^p1*y_t2(l2)^p2*...
......@@ -50,7 +50,7 @@ for jm = 1:options_mom_.mom.mom_nbr
powers = matched_moments_{jm,3};
for jv = 1:length(vars)
jvar = (oo_.dr.obs_var == vars(jv));
y = NaN(T,1); %Take care of T_eff instead of T for lags and NaN via nanmean below
y = NaN(T,1); %Take care of T_eff instead of T for lags and NaN via mean with 'omitnan' option below
y( (1-min(leadlags(jv),0)) : (T-max(leadlags(jv),0)), 1) = data( (1+max(leadlags(jv),0)) : (T+min(leadlags(jv),0)), jvar).^powers(jv);
if jv==1
m_data_tmp = y;
......
function [fval, info, exit_flag, junk1, junk2, oo_, M_, options_mom_] = method_of_moments_objective_function(xparam1, Bounds, oo_, estim_params_, matched_moments_, M_, options_mom_)
% [fval, info, exit_flag, oo_, M_, options_mom_] = method_of_moments_objective_function(xparam1, Bounds, oo_, estim_params_, matched_moments_, M_, options_mom_)
% [fval, info, exit_flag, junk1, junk2, oo_, M_, options_mom_] = method_of_moments_objective_function(xparam1, Bounds, oo_, estim_params_, matched_moments_, M_, options_mom_)
% -------------------------------------------------------------------------
% This function evaluates the objective function for GMM/SMM estimation
% =========================================================================
......@@ -8,19 +8,19 @@ function [fval, info, exit_flag, junk1, junk2, oo_, M_, options_mom_] = method_o
% o Bounds: structure containing parameter bounds
% o oo_: structure for results
% o estim_params_: structure describing the estimated_parameters
% o matched_moments_: structure containing information about selected moments to match in estimation (matched_moments_)
% o matched_moments_: structure containing information about selected moments to match in estimation
% o M_ structure describing the model
% o options_mom_: structure information about all settings (specified by the user, preprocessor, and taken from global options_)
% -------------------------------------------------------------------------
% 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 no error, 1 of error
% o exit_flag: 0 if error, 1 if no error
% o junk1: empty matrix required for optimizer interface
% o junk2: empty matrix required for optimizer interface
% 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.Q value of the quadratic form of the moment difference
% o M_: Matlab's structure describing the model
% -------------------------------------------------------------------------
% This function is called by
......@@ -53,7 +53,6 @@ function [fval, info, exit_flag, junk1, junk2, oo_, M_, options_mom_] = method_o
% o Willi Mutschler (willi@mutschler.eu)
% o Johannes Pfeifer (jpfeifer@uni-koeln.de)
% =========================================================================
% To Do: check penalized estimation for different optimizers, what is special about mode_compute=1 [@wmutschl]
%------------------------------------------------------------------------------
% 0. Initialization of the returned variables and others...
......
function [W_opt, normalization_factor]= method_of_moments_optimal_weighting_matrix(m_data, moments, q_lag)
function W_opt = method_of_moments_optimal_weighting_matrix(m_data, moments, q_lag)
% W_opt = method_of_moments_optimal_weighting_matrix(m_data, moments, q_lag)
% -------------------------------------------------------------------------
% This function computes the optimal weigthing matrix by a Bartlett kernel with maximum lag q_lag
......@@ -6,8 +6,8 @@ function [W_opt, normalization_factor]= method_of_moments_optimal_weighting_matr
% 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.
% =========================================================================
% INPUTS
% o m_data [T x numMom] selected empirical or theoretical moments at each point in time
% o moments [numMom x 1] mean of selected empirical or theoretical moments
% o m_data [T x numMom] selected data moments at each point in time
% o moments [numMom x 1] selected estimated moments (either data_moments or estimated model_moments)
% o q_lag [integer] Bartlett kernel maximum lag order
% -------------------------------------------------------------------------
% OUTPUTS
......@@ -67,13 +67,6 @@ end
% The estimate of W
W_opt = S\eye(size(S,1));
% Check positive definite W
try
chol(W_opt);
catch err
error('method_of_moments: The optimal weighting matrix is not positive definite. Check whether your model implies stochastic singularity.\n')
end
end
% The correlation matrix
......
......@@ -89,17 +89,16 @@ 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
WW = oo_.mom.Sw'*oo_.mom.Sw;
Asympt_Var = 1/T*((D'*WW*D)\eye(dim_params));
Asympt_Var = 1/T*((D'*WW*D)\eye(dim_params));
else
% We do not have the optimal weighting matrix yet
WWused = oo_.mom.Sw'*oo_.mom.Sw;
WWopt = method_of_moments_optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag);
S = WWopt\eye(size(WWopt,1));
AA = (D'*WWused*D)\eye(dim_params);
Asympt_Var = 1/T*AA*D'*WWused*S*WWused*D*AA;
% We do not have the optimal weighting matrix yet
WWopt = method_of_moments_optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.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));
\ No newline at end of file
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment