Commit 97de83c8 authored by Johannes Pfeifer 's avatar Johannes Pfeifer

Further changes to method of moments

- move scaling to objective functions instead of weighting matrix
- fix prefilter option
- Implement iterative GMM
parent 9e746fc5
Pipeline #3960 failed with stages
in 17 seconds
......@@ -105,7 +105,7 @@ for plt = 1:nbplt
z1 = l1:((xparam(kk)-l1)/(options_.mode_check.number_of_points/2)):xparam(kk);
z2 = xparam(kk):((l2-xparam(kk))/(options_.mode_check.number_of_points/2)):l2;
z = union(z1,z2);
if options_.penalized_estimator
if options_.mom.penalized_estimator
y = zeros(length(z),2);
dy=(xx-bayestopt_.p1)'/diag(bayestopt_.p2.^2)*(xx-bayestopt_.p1);
else
......@@ -122,7 +122,7 @@ for plt = 1:nbplt
fprintf('mode_check:: could not solve model for parameter %s at value %4.3f, error code: %u\n',name,z(i),info(1))
end
end
if options_.penalized_estimator
if options_.mom.penalized_estimator
prior=(xx-bayestopt_.p1)'/diag(bayestopt_.p2.^2)*(xx-bayestopt_.p1);
y(i,2) = (y(i,1)+prior-dy);
end
......@@ -152,7 +152,7 @@ for plt = 1:nbplt
hold off
drawnow
end
if options_.penalized_estimator
if options_.mom.penalized_estimator
if isoctave
axes('outerposition',[0.3 0.93 0.42 0.07],'box','on'),
else
......
......@@ -105,7 +105,7 @@ if info(1)
end
end
if strcmp(options_mom_.mom_method,'GMM')
if strcmp(options_mom_.mom.mom_method,'GMM')
%--------------------------------------------------------------------------
% 3. Set up pruned state-space system and compute model moments
%--------------------------------------------------------------------------
......@@ -134,7 +134,7 @@ if strcmp(options_mom_.mom_method,'GMM')
end
% Lead/lags covariance
if isfield(options_mom_.index,'E_yyt') && nnz(options_mom_.index.E_yyt) > 0
if options_mom_.centered_moments
if options_mom_.prefilter
E_yyt = pruned_state_space.Var_yi;
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)]);
......@@ -143,7 +143,7 @@ if strcmp(options_mom_.mom_method,'GMM')
oo_.mom.model_moments(offset+(1:E_yyt_nbr),1) = E_yyt(options_mom_.index.E_yyt);
end
elseif strcmp(options_mom_.mom_method,'SMM')
elseif strcmp(options_mom_.mom.mom_method,'SMM')
%------------------------------------------------------------------------------
% 3. Compute Moments of the model solution for normal innovations
%------------------------------------------------------------------------------
......@@ -195,16 +195,16 @@ end
% 4. Compute quadratic target function
%--------------------------------------------------------------------------
moments_difference = oo_.mom.data_moments - oo_.mom.model_moments;
residuals = oo_.mom.Sw*moments_difference;
residuals = sqrt(options_mom_.mom.weighting_matrix_scaling_factor)*oo_.mom.Sw*moments_difference;
oo_.mom.Q = residuals'*residuals;
if options_mom_.vector_output == 1 % lsqnonlin requires vector output
fval = residuals;
if options_mom_.penalized_estimator
if options_mom_.mom.penalized_estimator
fval=[fval;(xparam1-oo_.prior.mean)./sqrt(diag(oo_.prior.variance))];
end
else
fval = oo_.mom.Q;
if options_mom_.penalized_estimator
if options_mom_.mom.penalized_estimator
fval=fval+(xparam1-oo_.prior.mean)'/oo_.prior.variance*(xparam1-oo_.prior.mean);
end
end
......
......@@ -74,10 +74,6 @@ catch err
error('method_of_moments: The optimal weighting matrix is not positive definite. Check whether your model implies stochastic singularity.\n')
end
%normalize scale of W_opt to avoid numerical issues
normalization_factor=1000/max(diag(W_opt));
W_opt=W_opt*normalization_factor;
end
% The correlation matrix
......
......@@ -56,7 +56,7 @@ function [SE_values, Asympt_Var] = method_of_moments_standard_errors(xparam, obj
num_mom = size(oo_.mom.model_moments,1);
dim_params = size(xparam,1);
D = zeros(num_mom,dim_params);
eps_value = options_mom_.dynatol.x;
eps_value = options_mom_.mom.se_tolx;
for i=1:dim_params
%Positive step
......@@ -91,12 +91,12 @@ end
if Wopt_flag
% We have the optimal weighting matrix
WW = oo_.mom.Sw'*oo_.mom.Sw/oo_.mom.W_mat_normalization_factor;
WW = oo_.mom.Sw'*oo_.mom.Sw;
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/oo_.mom.W_mat_normalization_factor;
WWopt = method_of_moments_optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.bartlett_kernel_lag);
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;
......
......@@ -19,7 +19,7 @@
% =========================================================================
% Define testscenario
@#define orderApp = 3
@#define orderApp = 1
@#define estimParams = 1
% Note that we will set the numerical optimization tolerance levels very large to speed up the testsuite
......@@ -29,7 +29,7 @@
@#include "RBC_MoM_common.inc"
shocks;
var u_a = 0.0072;
var u_a; stderr 0.0072;
end;
varobs n c iv;
......@@ -37,13 +37,13 @@ varobs n c iv;
@#if estimParams == 0
estimated_params;
DELTA, 0.02;
BETTA, 0.9;
B, 0.4;
DELTA, 0.025;
BETTA, 0.98;
B, 0.45;
%ETAl, 1;
ETAc, 1.5;
ALFA, 0.6;
RHOA, 0.9;
ETAc, 1.8;
ALFA, 0.65;
RHOA, 0.95;
stderr u_a, 0.01;
%THETA, 3.48;
end;
......@@ -51,38 +51,40 @@ end;
@#if estimParams == 1
estimated_params;
DELTA, 0.02, 0, 1;
BETTA, 0.90, 0, 1;
B, 0.40, 0, 1;
DELTA, , 0, 1;
BETTA, , 0, 1;
B, , 0, 1;
%ETAl, 1, 0, 10;
ETAc, 1.80, 0, 10;
ALFA, 0.60, 0, 1;
RHOA, 0.90, 0, 1;
stderr u_a, 0.01, 0, 1;
ETAc, , 0, 10;
ALFA, , 0, 1;
RHOA, , 0, 1;
stderr u_a, , 0, 1;
%THETA, 3.48, 0, 10;
end;
@#endif
@#if estimParams == 2
estimated_params;
DELTA, 0.02, 0, 1, normal_pdf, 0.02, 0.5;
BETTA, 0.90, 0, 1, beta_pdf, 0.90, 0.25;
B, 0.40, 0, 1, normal_pdf, 0.40, 0.5;
DELTA, 0.025, 0, 1, normal_pdf, 0.02, 0.5;
BETTA, 0.98, 0, 1, beta_pdf, 0.90, 0.25;
B, 0.45, 0, 1, normal_pdf, 0.40, 0.5;
%ETAl, 1, 0, 10, normal_pdf, 0.25, 0.0.1;
ETAc, 1.80, 0, 10, normal_pdf, 1.80, 0.5;
ALFA, 0.60, 0, 1, normal_pdf, 0.60, 0.5;
RHOA, 0.90, 0, 1, normal_pdf, 0.90, 0.5;
stderr u_a, -0.99, 0, 1, normal_pdf, 0.01, 0.5;
ETAc, 1.8, 0, 10, normal_pdf, 1.80, 0.5;
ALFA, 0.65, 0, 1, normal_pdf, 0.60, 0.5;
RHOA, 0.95, 0, 1, normal_pdf, 0.90, 0.5;
stderr u_a, 0.01, 0, 1, normal_pdf, 0.01, 0.5;
%THETA, 3.48, 0, 10, normal_pdf, 0.25, 0.0.1;
end;
@#endif
% Simulate data
stoch_simul(order=@{orderApp},pruning,nodisplay,nomoments,periods=750,drop=500);
stoch_simul(order=@{orderApp},pruning,nodisplay,nomoments,periods=500);
save('RBC_MoM_data_@{orderApp}.mat', options_.varobs{:} );
pause(1);
estimated_params_init(use_calibration);
end;
%--------------------------------------------------------------------------
% Method of Moments Estimation
......@@ -130,8 +132,9 @@ matched_moments_ = {
% , penalized_estimator % use penalized optimization
, pruning % use pruned state space system at higher-order
% , verbose % display and store intermediate estimation results
, weighting_matrix = OPTIMAL % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename
, additional_optimizer_steps = [4] % vector of additional mode-finders run after mode_compute
, weighting_matrix = ['optimal','optimal'] % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename
, weighting_matrix_scaling_factor=1
%, additional_optimizer_steps = [4] % vector of additional mode-finders run after mode_compute
% , prefilter=0 % demean each data series by its empirical mean and use centered moments
%
% Options for SMM
......
......@@ -25,7 +25,7 @@
@#include "RBC_MoM_common.inc"
shocks;
var u_a = 0.0072;
var u_a; stderr 0.0072;
var n; stderr 0.01;
end;
......
% Tests SMM and GMM routines with prefilter, explicit initialization, and estimated_params_init(use_calibration);
%
% Copyright (C) 2020 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 <http://www.gnu.org/licenses/>.
% =========================================================================
% Define testscenario
@#define orderApp = 2
% Note that we will set the numerical optimization tolerance levels very large to speed up the testsuite
@#define optimizer = 13
@#include "RBC_MoM_common.inc"
shocks;
var u_a; stderr 0.0072;
end;
varobs n c iv;
% Simulate data
stoch_simul(order=@{orderApp},pruning,nodisplay,nomoments,periods=250,TeX);
save('RBC_MoM_data_@{orderApp}.mat', options_.varobs{:} );
pause(1);
set_param_value('DELTA',NaN);
estimated_params;
DELTA, 0.025, 0, 1;
BETTA, , 0, 1;
B, , 0, 1;
%ETAl, 1, 0, 10;
ETAc, , 0, 10;
ALFA, , 0, 1;
RHOA, , 0, 1;
stderr u_a, , 0, 1;
%THETA, 3.48, 0, 10;
end;
estimated_params_init(use_calibration);
end;
%--------------------------------------------------------------------------
% Method of Moments Estimation
%--------------------------------------------------------------------------
% matched_moments blocks : We don't have an interface yet
% get indices in declaration order
ic = strmatch('c', M_.endo_names,'exact');
iiv = strmatch('iv', M_.endo_names,'exact');
in = strmatch('n', M_.endo_names,'exact');
% first entry: number of variable in declaration order
% second entry: lag
% third entry: power
matched_moments_ = {
[ic ] [0 ], [1 ];
[in ] [0 ], [1 ];
[iiv ] [0 ], [1 ];
[ic ic ] [0 0], [1 1];
[ic iiv] [0 0], [1 1];
[ic in ] [0 0], [1 1];
[iiv ic ] [0 0], [1 1];
[iiv iiv] [0 0], [1 1];
[iiv in ] [0 0], [1 1];
% [in ic ] [0 0], [1 1];
% [in iiv] [0 0], [1 1];
[in in ] [0 0], [1 1];
[ic ic ] [0 -1], [1 1];
[in in ] [0 -1], [1 1];
[iiv iiv] [0 -1], [1 1];
% [iiv iiv] [0 -1], [1 1];
};
weighting_matrix=diag([1000;ones(6,1)]);
save('test_matrix.mat','weighting_matrix')
@#for mommethod in ["GMM", "SMM"]
method_of_moments(
% Necessery options
mom_method = @{mommethod} % method of moments method; possible values: GMM|SMM
, datafile = 'RBC_MoM_data_@{orderApp}.mat' % name of filename with data
% Options for both GMM and SMM
% , bartlett_kernel_lag = 20 % bandwith in optimal weighting matrix
, order = @{orderApp} % order of Taylor approximation in perturbation
% , penalized_estimator % use penalized optimization
, pruning % use pruned state space system at higher-order
% , verbose % display and store intermediate estimation results
% , weighting_matrix = 'test_matrix.mat' % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename
, weighting_matrix =['test_matrix.mat','optimal']
%, weighting_matrix = optimal % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename
%, additional_optimizer_steps = [4] % vector of additional mode-finders run after mode_compute
, prefilter=1 % demean each data series by its empirical mean and use centered moments
, se_tolx=1e-5
%
% Options for SMM
% , bounded_shock_support % trim shocks in simulation to +- 2 stdev
, burnin = 500 % number of periods dropped at beginning of simulation
% , seed = 24051986 % seed used in simulations
% , simulation_multiple = 5 % multiple of the data length used for simulation
%
% General options
%, dirname = 'MM' % directory in which to store estimation output
% , graph_format = EPS % specify the file format(s) for graphs saved to disk
% , nodisplay % do not display the graphs, but still save them to disk
% , nograph % do not create graphs (which implies that they are not saved to the disk nor displayed)
% , noprint % do not print stuff to console
% , plot_priors = 1 % control plotting of priors
% , prior_trunc = 1e-10 % probability of extreme values of the prior density that is ignored when computing bounds for the parameters
% , TeX % print TeX tables and graphics
%
% Data and model options
%, first_obs = 501 % number of first observation
% , logdata % if loglinear is set, this option is necessary if the user provides data already in logs, otherwise the log transformation will be applied twice (this may result in complex data)
% , loglinear % computes a log-linear approximation of the model instead of a linear approximation
%, nobs = 500 % number of observations
% , xls_sheet = willi % name of sheet with data in Excel
% , xls_range = B2:D200 % range of data in Excel sheet
%
% Optimization options that can be set by the user in the mod file, otherwise default values are provided
% , analytic_derivation % uses analytic derivatives to compute standard errors for GMM
%, huge_number=1D10 % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
, mode_compute = @{optimizer} % specifies the optimizer for minimization of moments distance, note that by default there is a new optimizer
%, optim = ('TolFun', 1e-3
% ,'TolX', 1e-5
% ) % 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
%
% % Numerical algorithms options
% , aim_solver % Use AIM algorithm to compute perturbation approximation
% , dr=default % method used to compute the decision rule; possible values are DEFAULT, CYCLE_REDUCTION, LOGARITHMIC_REDUCTION
% , dr_cycle_reduction_tol = 1e-7 % convergence criterion used in the cycle reduction algorithm
% , dr_logarithmic_reduction_maxiter = 100 % maximum number of iterations used in the logarithmic reduction algorithm
% , dr_logarithmic_reduction_tol = 1e-12 % convergence criterion used in the cycle reduction algorithm
% , k_order_solver % use k_order_solver in higher order perturbation approximations
% , lyapunov = DEFAULT % algorithm used to solve lyapunov equations; possible values are DEFAULT, FIXED_POINT, DOUBLING, SQUARE_ROOT_SOLVER
% , lyapunov_complex_threshold = 1e-15 % complex block threshold for the upper triangular matrix in symmetric Lyapunov equation solver
% , lyapunov_fixed_point_tol = 1e-10 % convergence criterion used in the fixed point Lyapunov solver
% , lyapunov_doubling_tol = 1e-16 % convergence criterion used in the doubling algorithm
% , sylvester = default % algorithm to solve Sylvester equation; possible values are DEFAULT, FIXED_POINT
% , 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
);
@#endfor
\ 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