Commit 3615962a authored by Willi Mutschler's avatar Willi Mutschler Committed by Johannes Pfeifer

First draft of method of moments toolbox with GMM and SMM

parent 77dbeab9
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
function [dataMoments, m_data] = method_of_moments_datamoments(data, DynareResults, MatchedMoments, OptionsMoM)
% [dataMoments, m_data] = method_of_moments_datamoments(data, DynareResults, MatchedMoments, OptionsMoM)
% This function computes the user-selected empirical moments from data
% =========================================================================
% INPUTS
% o data [T x varobs_nbr] data set
% o DynareResults: [structure] storage for results (oo_)
% o MatchedMoments: [structure] information about selected moments to match in estimation (matched_moments_)
% o OptionsMoM: [structure] information about all settings (specified by the user, preprocessor, and taken from global options_)
% -------------------------------------------------------------------------
% OUTPUTS
% o dataMoments [numMom x 1] mean of selected empirical moments
% o m_data [T x numMom] selected empirical moments at each point in time
% -------------------------------------------------------------------------
% This function is called by
% o method_of_moments.m
% o method_of_moments_SMM.m
% =========================================================================
% 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/>.
% -------------------------------------------------------------------------
% Author(s):
% o Willi Mutschler (willi@mutschler.eu)
% o Johannes Pfeifer (jpfeifer@uni-koeln.de)
% =========================================================================
% Initialization
T = size(data,1); % Number of observations (T) and number of observables (ny)
mom_nbr = OptionsMoM.mom_nbr;
dataMoments = nan(mom_nbr,1);
m_data = nan(T,mom_nbr);
% Product moment for each time period, i.e. each row t contains yt1(l1)^p1*yt2(l2)^p2*...
% note that here we already are able to treat leads and lags and any power product moments
for jm = 1:mom_nbr
vars = DynareResults.dr.inv_order_var(MatchedMoments{jm,1})';
leadlags = MatchedMoments{jm,2}; % note that lags are negative numbers and leads are positive numbers
powers = MatchedMoments{jm,3};
for jv = 1:length(vars)
jvar = DynareResults.dr.obs_var == vars(jv);
y = nan(T,1);
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;
else
m_data_tmp = m_data_tmp.*y;
end
end
dataMoments(jm,1) = sum(m_data_tmp,'omitnan')/(T-sum(abs(leadlags)));
% We replace nan (due to leads and lags) with the corresponding mean
% @wmutschl: this should also work for missing values, right?
m_data_tmp(isnan(m_data_tmp)) = dataMoments(jm,1);
m_data(:,jm) = m_data_tmp;
end
end %function end
function Wopt = method_of_moments_optimal_weighting_matrix(m_data, moments, qLag)
% Wopt = method_of_moments_optimal_weighting_matrix(m_data, moments, qLag)
% -------------------------------------------------------------------------
% This function computes the optimal weigthing matrix by a Bartlett kernel with maximum lag qlag
% Adapted from replication codes of
% 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 qlag [integer] Bartlett kernel maximum lag order
% -------------------------------------------------------------------------
% OUTPUTS
% o Wopt [numMom x numMom] optimal weighting matrix
% -------------------------------------------------------------------------
% This function is called by
% o method_of_moments.m
% -------------------------------------------------------------------------
% This function calls:
% o CorrMatrix (embedded)
% =========================================================================
% 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/>.
% -------------------------------------------------------------------------
% Author(s):
% o Willi Mutschler (willi@mutschler.eu)
% o Johannes Pfeifer (jpfeifer@uni-koeln.de)
% =========================================================================
% Initialize
[T,numMom] = size(m_data); %note that in m_data nan values (due to leads or lags in matchedmoments) are removed so T is the effective sample size
% center around moments (could be either datamoments or modelmoments)
hFunc = m_data - repmat(moments',T,1);
% The required correlation matrices
GAMA_array = zeros(numMom,numMom,qLag);
GAMA0 = CorrMatrix(hFunc,T,numMom,0);
if qLag > 0
for ii=1:qLag
GAMA_array(:,:,ii) = CorrMatrix(hFunc,T,numMom,ii);
end
end
% The estimate of S
S = GAMA0;
if qLag > 0
for ii=1:qLag
S = S + (1-ii/(qLag+1))*(GAMA_array(:,:,ii) + GAMA_array(:,:,ii)');
end
end
% The estimate of W
Wopt = S\eye(size(S,1));
% Check positive definite W
try
chol(Wopt);
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
function GAMAcorr = CorrMatrix(hFunc,T,numMom,v)
GAMAcorr = zeros(numMom,numMom);
for t = 1+v:T
GAMAcorr = GAMAcorr + hFunc(t-v,:)'*hFunc(t,:);
end
GAMAcorr = GAMAcorr/T;
end
\ No newline at end of file
function [SEvalues, AVar] = method_of_moments_standard_errors(xparam, objective_function, Bounds, DynareResults, EstimatedParameters, MatchedMoments, Model, OptionsMoM, Wopt_flag)
% [SEvalues, AVar] = method_of_moments_standard_errors(xparam, objective_function, Bounds, DynareResults, EstimatedParameters, MatchedMoments, Model, OptionsMoM, Wopt_flag)
% -------------------------------------------------------------------------
% This function computes standard errors to the method of moments estimates
% Adapted from replication codes of
% 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 xparam: value of estimated parameters as returned by set_prior()
% o objective_function string of objective function, either method_of_moments_GMM.m or method_of_moments_SMM.m
% o Bounds: structure containing parameter bounds
% o DynareResults: structure for results (oo_)
% o EstimatedParameters: structure describing the estimated_parameters (estim_params_)
% o MatchedMoments: structure containing information about selected moments to match in estimation (matched_moments_)
% o Model structure describing the Model
% o OptionsMoM: structure information about all settings (specified by the user, preprocessor, and taken from global options_)
% o Wopt_flag: indicator whether the optimal weighting is actually used
% -------------------------------------------------------------------------
% OUTPUTS
% o SEvalues [nparam x 1] vector of standard errors
% o AVar [nparam x nparam] asymptotic covariance matrix
% -------------------------------------------------------------------------
% This function is called by
% o method_of_moments.m
% -------------------------------------------------------------------------
% This function calls:
% o get_the_name
% o get_error_message
% o method_of_moments_GMM.m (objective function)
% o method_of_moments_SMM.m (objective function)
% o method_of_moments_optimal_weighting_matrix
% =========================================================================
% 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/>.
% -------------------------------------------------------------------------
% Author(s):
% o Willi Mutschler (willi@mutschler.eu)
% o Johannes Pfeifer (jpfeifer@uni-koeln.de)
% =========================================================================
% Some dimensions
numMom = size(DynareResults.mom.modelMoments,1);
dimParams = size(xparam,1);
D = zeros(numMom,dimParams);
epsValue = OptionsMoM.dynatol.x;
for i=1:dimParams
%Positive step
xparam_eps_p = xparam;
xparam_eps_p(i,1) = xparam_eps_p(i) + epsValue;
[~, info_p, exit_flag_p, DynareResults_p, ~, ~] = feval(objective_function, xparam_eps_p, Bounds, DynareResults, EstimatedParameters, MatchedMoments, Model, OptionsMoM);
% Negative step
xparam_eps_m = xparam;
xparam_eps_m(i,1) = xparam_eps_m(i) - epsValue;
[~, info_m, exit_flag_m, DynareResults_m, ~, ~] = feval(objective_function, xparam_eps_m, Bounds, DynareResults, EstimatedParameters, MatchedMoments, Model, OptionsMoM);
% The Jacobian:
if nnz(info_p)==0 && nnz(info_m)==0
D(:,i) = (DynareResults_p.mom.modelMoments - DynareResults_m.mom.modelMoments)/(2*epsValue);
else
problpar = get_the_name(i,OptionsMoM.TeX, Model, EstimatedParameters, OptionsMoM);
message_p = get_error_message(info_p, OptionsMoM);
message_m = get_error_message(info_m, OptionsMoM);
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)
AVar = NaN(length(xparam),length(xparam));
SEvalues = NaN(length(xparam),1);
return
end
end
T = OptionsMoM.nobs; %Number of observations
if isfield(OptionsMoM,'variance_correction_factor')
T = T*OptionsMoM.variance_correction_factor;
end
if Wopt_flag
% We have the optimal weighting matrix
WW = DynareResults.mom.Sw'*DynareResults.mom.Sw;
AVar = 1/T*((D'*WW*D)\eye(dimParams));
else
% We do not have the optimal weighting matrix yet
WWused = DynareResults.mom.Sw'*DynareResults.mom.Sw;
WWopt = method_of_moments_optimal_weighting_matrix(DynareResults.mom.m_data, DynareResults.mom.modelMoments, OptionsMoM.bartlett_kernel_lag);
S = WWopt\eye(size(WWopt,1));
AA = (D'*WWused*D)\eye(dimParams);
AVar = 1/T*AA*D'*WWused*S*WWused*D*AA;
end
SEvalues = sqrt(diag(AVar));
\ No newline at end of file
......@@ -428,7 +428,7 @@ switch minimizer_algorithm
case 12
if isoctave
error('Option mode_compute=12 is not available under Octave')
elseif ~user_has_matlab_license('global_optimization_toolbox')
elseif ~user_has_matlab_license('GADS_Toolbox')
error('Option mode_compute=12 requires the Global Optimization Toolbox')
end
[LB, UB] = set_bounds_to_finite_values(bounds, options_.huge_number);
......@@ -523,6 +523,21 @@ switch minimizer_algorithm
end
func = @(x)objective_function(x,varargin{:});
[opt_par_values,fval,exitflag,output] = simulannealbnd(func,start_par_value,bounds(:,1),bounds(:,2),optim_options);
case 13
% Matlab's lsqnonlin (Optimization toolbox needed).
if isoctave && ~user_has_octave_forge_package('optim')
error('Option mode_compute=13 requires the optim package')
elseif ~isoctave && ~user_has_matlab_license('optimization_toolbox')
error('Option mode_compute=13 requires the Optimization Toolbox')
end
optim_options = optimset('display','iter','MaxFunEvals',5000,'MaxIter',5000,'TolFun',1e-6,'TolX',1e-6);
if ~isempty(options_.optim_opt)
eval(['optim_options = optimset(optim_options,' options_.optim_opt ');']);
end
if options_.silent_optimizer
optim_options = optimset(optim_options,'display','off');
end
[opt_par_values,Resnorm,fval,exitflag,OUTPUT,LAMBDA,JACOB] = lsqnonlin(objective_function,start_par_value,bounds(:,1),bounds(:,2),optim_options,varargin{:});
otherwise
if ischar(minimizer_algorithm)
if exist(minimizer_algorithm)
......
......@@ -47,6 +47,8 @@ MODFILES = \
estimation/MH_recover/fs2000_recover_3.mod \
estimation/t_proposal/fs2000_student.mod \
estimation/tune_mh_jscale/fs2000.mod \
estimation/method_of_moments/AnScho_MoM.mod \
estimation/method_of_moments/RBCmodel_MoM.mod \
moments/example1_var_decomp.mod \
moments/example1_bp_test.mod \
moments/test_AR1_spectral_density.mod \
......
This diff is collapsed.
This diff is collapsed.
% By Willi Mutschler, September 26, 2016. Email: willi@mutschler.eu
function [ys,params,check] = RBCmodel_steadystate(ys,exo,M_,options_)
%% Step 0: initialize indicator and set options for numerical solver
check = 0;
options = optimset('Display','off','TolX',1e-12,'TolFun',1e-12);
params = M_.params;
%% Step 1: read out parameters to access them with their name
for ii = 1:M_.param_nbr
eval([ M_.param_names{ii} ' = M_.params(' int2str(ii) ');']);
end
%% Step 2: Check parameter restrictions
if ETAc*ETAl<1 % parameter violates restriction (here it is artifical)
check=1; %set failure indicator
return; %return without updating steady states
end
%% Step 3: Enter model equations here
A = 1;
RK = 1/BETTA - (1-DELTA);
K_O_N = (RK/(A*(1-ALFA)))^(-1/ALFA);
if K_O_N <= 0
check = 1; % set failure indicator
return; % return without updating steady states
end
W = A*ALFA*(K_O_N)^(1-ALFA);
IV_O_N = DELTA*K_O_N;
Y_O_N = A*K_O_N^(1-ALFA);
C_O_N = Y_O_N - IV_O_N;
if C_O_N <= 0
check = 1; % set failure indicator
return; % return without updating steady states
end
% The labor level
if ETAc == 1 && ETAl == 1
N = (1-BETTA*B)*(C_O_N*(1-B))^-1*W/THETA/(1+(1-BETTA*B)*(C_O_N*(1-B))^-1*W/THETA);
else
% No closed-form solution use a fixed-point algorithm
N0 = 1/3;
[N,~,exitflag] = fsolve(@(N) THETA*(1-N)^(-ETAl)*N^ETAc - (1-BETTA*B)*(C_O_N*(1-B))^(-ETAc)*W, N0,options);
if exitflag <= 0
check = 1; % set failure indicator
return % return without updating steady states
end
end
C=C_O_N*N;
Y=Y_O_N*N;
IV=IV_O_N*N;
K=K_O_N*N;
LA = (C-B*C)^(-ETAc)-BETTA*B*(C-B*C)^(-ETAc);
k=log(K);
c=log(C);
a=log(A);
iv=log(IV);
y=log(Y);
la=log(LA);
n=log(N);
rk=log(RK);
w=log(W);
%% Step 4: Update parameters and variables
params=NaN(M_.param_nbr,1);
for iter = 1:M_.param_nbr %update parameters set in the file
eval([ 'params(' num2str(iter) ') = ' M_.param_names{iter} ';' ])
end
for ii = 1:M_.orig_endo_nbr %auxiliary variables are set automatically
eval(['ys(' int2str(ii) ') = ' M_.endo_names{ii} ';']);
end
end
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