//wip GMM first stage working

Added moment_estimation to preprocesser
Started rebase of Andreasen branch (a7cbc637)
parent ea4d3f4d
Pipeline #3686 failed with stages
in 5 seconds
This diff is collapsed.
function [fval,info,exit_flag,moments_difference,modelMoments,junk1,junk2,M,options_mom,mom_info,oo]...
= method_of_moments_GMM(xparam1, DynareDataset, options_mom, M, estim_params, mom_info, bounds, oo)
% [fval,info,exit_flag,moments_difference,modelMoments,junk1,junk2,Model,DynareOptions,GMMinfo,DynareResults]...
% = GMM_Objective_Function(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,GMMinfo,BoundsInfo,DynareResults)
% This function evaluates the objective function for GMM estimation
%
% INPUTS
% o xparam1: initial value of estimated parameters as returned by set_prior()
% o DynareDataset: data after required transformation
% o DynareOptions: Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
% o Model Matlab's structure describing the Model (initialized by dynare, see @ref{M_}).
% o EstimatedParameters: Matlab's structure describing the estimated_parameters (initialized by dynare, see @ref{estim_params_}).
% o GMMInfo Matlab's structure describing the GMM settings (initialized by dynare, see @ref{bayesopt_}).
% o BoundsInfo Matlab's structure containing prior bounds
% o DynareResults Matlab's structure gathering the results (initialized by dynare, see @ref{oo_}).
%
% OUTPUTS
% o fval: value of the quadratic form of the moment difference
% o info: vector storing error code and penalty
% o exit_flag: 0 if no error, 1 of error
% o moments_difference: [numMom x 1] vector with difference of empirical and model moments
% o modelMoments: [numMom x 1] vector with model moments
% o Model: Matlab's structure describing the Model (initialized by dynare, see @ref{M_}).
% o DynareOptions: Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
% o GMMinfo: Matlab's structure describing the GMM parameter options (initialized by dynare, see @ref{GMMinfo_}).
% o DynareResults: Matlab's structure gathering the results (initialized by dynare, see @ref{oo_}).
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2013-17 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/>.
% Initialization of the returned variables and others...
fval = NaN;
exit_flag = 1;
info = 0;
junk2 = [];
junk1 = [];
moments_difference = NaN(mom_info.numMom,1);
modelMoments = NaN(mom_info.numMom,1);
%------------------------------------------------------------------------------
% 1. Get the structural parameters & define penalties
%------------------------------------------------------------------------------
% Return, with endogenous penalty, if some parameters are smaller than the lower bound of the parameters.
% Why is mode compute == 1 excluded??? @wmutschl
if ~isequal(options_mom.mode_compute,1) && any(xparam1<bounds.lb)
k = find(xparam1<bounds.lb);
fval = Inf;
exit_flag = 0;
info(1) = 41;
info(4)= sum((bounds.lb(k)-xparam1(k)).^2);
return
end
% Return, with endogenous penalty, if some parameters are greater than the upper bound of the parameters.
if ~isequal(options_mom.mode_compute,1) && any(xparam1>bounds.ub)
k = find(xparam1>bounds.ub);
fval = Inf;
exit_flag = 0;
info(1) = 42;
info(4)= sum((xparam1(k)-bounds.ub(k)).^2);
return
end
% Set all parameters
M = set_all_parameters(xparam1,estim_params,M);
% Test if Q is positive definite.
if ~issquare(M.Sigma_e) || estim_params.ncx || isfield(estim_params,'calibrated_covariances')
[Q_is_positive_definite, penalty] = ispd(M.Sigma_e(estim_params.Sigma_e_entries_to_check_for_positive_definiteness,estim_params.Sigma_e_entries_to_check_for_positive_definiteness));
if ~Q_is_positive_definite
fval = Inf;
exit_flag = 0;
info(1) = 43;
info(4) = penalty;
return
end
if isfield(estim_params,'calibrated_covariances')
correct_flag=check_consistency_covariances(M.Sigma_e);
if ~correct_flag
penalty = sum(M.Sigma_e(estim_params.calibrated_covariances.position).^2);
fval = Inf;
exit_flag = 0;
info(1) = 71;
info(4) = penalty;
return
end
end
end
%------------------------------------------------------------------------------
% 2. call resol to compute steady state and model solution
%------------------------------------------------------------------------------
% Compute linear approximation around the deterministic steady state and fill dr structure
[oo.dr, info, M, options_mom, oo] = resol(0, M, options_mom, oo);
% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol).
if info(1)
if info(1) == 3 || info(1) == 4 || info(1) == 5 || info(1)==6 ||info(1) == 19 ||...
info(1) == 20 || info(1) == 21 || info(1) == 23 || info(1) == 26 || ...
info(1) == 81 || info(1) == 84 || info(1) == 85 || info(1) == 86
%meaningful second entry of output that can be used
fval = Inf;
info(4) = info(2);
exit_flag = 0;
return
else
fval = Inf;
info(4) = 0.1;
exit_flag = 0;
return
end
end
% % check endogenous prior restrictions
% info=endogenous_prior_restrictions(T,R,Model,DynareOptions,DynareResults);
% if info(1)
% fval = Inf;
% info(4)=info(2);
% exit_flag = 0;
% return
% end
%------------------------------------------------------------------------------
% 3. Set up pruned state-space system and compute model moments
%------------------------------------------------------------------------------
useautocorr = 0; %@wmutschl_ need to add autocorr yet
pruned_state_space = pruned_state_space_system(M, options_mom, oo.dr, mom_info.varobs_index_dr, options_mom.ar, useautocorr, 0);
% Get the moments implied by the model solution that are matched
% Add centered moments @wmutschl
modelMoments = [];
% First moments
if nnz(mom_info.index_moms_dr.E_y) > 0
E_y = pruned_state_space.E_y;
modelMoments = [modelMoments; E_y(mom_info.index_moms_dr.E_y)];
end
% Second moments
if nnz(mom_info.index_moms_dr.E_yy) > 0
E_yy = pruned_state_space.Var_y + pruned_state_space.E_y*pruned_state_space.E_y';
modelMoments = [modelMoments; E_yy(tril(mom_info.index_moms_dr.E_yy))];
end
if nnz(mom_info.index_moms_dr.E_yyt) > 0
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)]);
modelMoments = [modelMoments; E_yyt(mom_info.index_moms_dr.E_yyt)];
end
%------------------------------------------------------------------------------
% 4. Compute quadratic target function using weighting matrix W
%------------------------------------------------------------------------------
moments_difference = mom_info.datamoments-modelMoments;
fval = moments_difference'*mom_info.W*moments_difference;
if options_mom.penalized_estimator
fval=fval+(xparam1-mom_info.p1)'/diag(mom_info.p2)*(xparam1-mom_info.p1);
end
end
function [Wopt] = method_of_moments_Optimal_Weighting_Matrix(xparam1, DynareDataset, options_mom, M, estim_params, mom_info, oo, bounds)
% [Wopt] = get_Optimal_Weighting_Matrix_GMM_SMM(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,GMM_SMM_info,DynareResults,BoundsInfo,GMM_SMM_indicator)
% This function computes the optimal weigthing matrix by a Bartlett kernel with maximum lag qlag
% INPUTS
% o xparam1: initial value of estimated parameters as returned by set_prior()
% o DynareDataset: data after required transformation
% o DynareOptions Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
% o Model Matlab's structure describing the Model (initialized by dynare, see @ref{M_}).
% o EstimatedParameters: Matlab's structure describing the estimated_parameters (initialized by dynare, see @ref{estim_params_}).
% o GMM_SMM_Info Matlab's structure describing the GMM settings (initialized by dynare, see @ref{bayesopt_}).
% o DynareResults Matlab's structure gathering the results (initialized by dynare, see @ref{oo_}).
% o BoundsInfo Matlab's structure containing prior bounds
% o GMM_SMM_indicator string indicating SMM or GMM
%
% OUTPUTS
% o Wopt [numMom x numMom] optimal weighting matrix
%
% SPECIAL REQUIREMENTS
% None.
% Copyright (C) 2013-17 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/>.
qLag=options_mom.bartlett_kernel_lag;
% We compute the h-function for all observations
T = DynareDataset.nobs;
% Evaluating the objective function to get modelMoments
if strcmp(options_mom.mom_method,'GMM')
[fval,info,exit_flag,moments_difference,modelMoments] = method_of_moments_GMM(xparam1, DynareDataset, options_mom, M, estim_params, mom_info, bounds, oo);
% centered around theoretical moments
hFunc = oo.(lower(GMM_SMM_indicator)).datamoments.m_data - repmat(modelMoments',T,1);
elseif strcmp(GMM_SMM_indicator,'SMM')
[fval,info,exit_flag,moments_difference,modelMoments] = SMM_Objective_Function(xparam1,DynareDataset,options_mom,M,estim_params,mom_info,bounds,oo);
% centered around data moments
hFunc = oo.(lower(GMM_SMM_indicator)).datamoments.m_data - repmat(mean(oo.(lower(GMM_SMM_indicator)).datamoments.m_data),T,1);
end
% The required correlation matrices
numMom = mom_info.numMom;
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
Wopt = S\eye(size(S,1));
try
chol(Wopt);
catch err
error('%s Error: The optimal weighting matrix is not positive definite. Check whether your model implies stochastic singularity\n',GMM_SMM_indicator)
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 [fval,info,exit_flag,moments_difference,modelMoments,junk1,junk2,Model,DynareOptions,SMMinfo,DynareResults]...
= SMM_Objective_Function(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,SMMinfo,BoundsInfo,DynareResults)
% [fval,info,exit_flag,moments_difference,modelMoments,junk1,junk2,Model,DynareOptions,SMMinfo,DynareResults]...
% = SMM_Objective_Function(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,SMMinfo,BoundsInfo,DynareResults)
% This function evaluates the objective function for SMM estimation
%
% INPUTS
% o xparam1: initial value of estimated parameters as returned by set_prior()
% o DynareDataset: data after required transformation
% o DynareOptions: Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
% o Model Matlab's structure describing the Model (initialized by dynare, see @ref{M_}).
% o EstimatedParameters: Matlab's structure describing the estimated_parameters (initialized by dynare, see @ref{estim_params_}).
% o SMMInfo Matlab's structure describing the SMM settings (initialized by dynare, see @ref{bayesopt_}).
% o BoundsInfo Matlab's structure containing prior bounds
% o DynareResults Matlab's structure gathering the results (initialized by dynare, see @ref{oo_}).
%
% OUTPUTS
% o fval: value of the quadratic form of the moment difference
% o moments_difference: [numMom x 1] vector with difference of empirical and model moments
% o modelMoments: [numMom x 1] vector with model moments
% o exit_flag: 0 if no error, 1 of error
% o info: vector storing error code and penalty
% o Model: Matlab's structure describing the Model (initialized by dynare, see @ref{M_}).
% o DynareOptions: Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
% o SMMinfo: Matlab's structure describing the SMM parameter options (initialized by dynare, see @ref{SMMinfo_}).
% o DynareResults: Matlab's structure gathering the results (initialized by dynare, see @ref{oo_}).
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2013-17 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/>.
global objective_function_penalty_base
% Initialization of the returned variables and others...
fval = NaN;
exit_flag = 1;
info = 0;
junk2 = [];
junk1 = [];
moments_difference=NaN(SMMinfo.numMom,1);
modelMoments=NaN(SMMinfo.numMom,1);
%------------------------------------------------------------------------------
% 1. Get the structural parameters & define penalties
%------------------------------------------------------------------------------
% Return, with endogenous penalty, if some parameters are smaller than the lower bound of the parameters.
if ~isequal(DynareOptions.mode_compute,1) && any(xparam1<BoundsInfo.lb)
k = find(xparam1<BoundsInfo.lb);
fval = Inf;
exit_flag = 0;
info(1) = 41;
info(4)= sum((BoundsInfo.lb(k)-xparam1(k)).^2);
return
end
% Return, with endogenous penalty, if some parameters are greater than the upper bound of the parameters.
if ~isequal(DynareOptions.mode_compute,1) && any(xparam1>BoundsInfo.ub)
k = find(xparam1>BoundsInfo.ub);
fval = Inf;
exit_flag = 0;
info(1) = 42;
info(4)= sum((xparam1(k)-BoundsInfo.ub(k)).^2);
return
end
% Set all parameters
Model = set_all_parameters(xparam1,EstimatedParameters,Model);
% Test if Q is positive definite.
if ~issquare(Model.Sigma_e) || EstimatedParameters.ncx || isfield(EstimatedParameters,'calibrated_covariances')
[Q_is_positive_definite, penalty] = ispd(Model.Sigma_e(EstimatedParameters.Sigma_e_entries_to_check_for_positive_definiteness,EstimatedParameters.Sigma_e_entries_to_check_for_positive_definiteness));
if ~Q_is_positive_definite
fval = Inf;
exit_flag = 0;
info(1) = 43;
info(4) = penalty;
return
end
if isfield(EstimatedParameters,'calibrated_covariances')
correct_flag=check_consistency_covariances(Model.Sigma_e);
if ~correct_flag
penalty = sum(Model.Sigma_e(EstimatedParameters.calibrated_covariances.position).^2);
fval = Inf;
exit_flag = 0;
info(1) = 71;
info(4) = penalty;
return
end
end
end
%------------------------------------------------------------------------------
% 2. call resol to compute steady state and model solution
%------------------------------------------------------------------------------
[dr_dynare_state_space,info,Model,DynareOptions,DynareResults] = resol(0,Model,DynareOptions,DynareResults);
% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol).
if info(1)
if info(1) == 3 || info(1) == 4 || info(1) == 5 || info(1)==6 ||info(1) == 19 ||...
info(1) == 20 || info(1) == 21 || info(1) == 23 || info(1) == 26 || ...
info(1) == 81 || info(1) == 84 || info(1) == 85 || info(1) == 86
%meaningful second entry of output that can be used
fval = Inf;
info(4) = info(2);
exit_flag = 0;
return
else
fval = Inf;
info(4) = 0.1;
exit_flag = 0;
return
end
end
% % check endogenous prior restrictions
% info=endogenous_prior_restrictions(T,R,Model,DynareOptions,DynareResults);
% if info(1)
% fval = Inf;
% info(4)=info(2);
% exit_flag = 0;
% return
% end
%------------------------------------------------------------------------------
% 3. Compute Moments of the model solution for normal innovations
%------------------------------------------------------------------------------
% create shock series with correct covariance matrix from iid standard
% normal shocks
i_exo_var = setdiff([1:Model.exo_nbr],find(diag(Model.Sigma_e) == 0 )); %find singular entries in covariance
chol_S = chol(Model.Sigma_e(i_exo_var,i_exo_var));
scaled_shock_series=zeros(size(DynareResults.smm.shock_series)); %initialize
scaled_shock_series(:,i_exo_var) = DynareResults.smm.shock_series(:,i_exo_var)*chol_S; %set non-zero entries
%% simulate series
y_sim = simult_(dr_dynare_state_space.ys,dr_dynare_state_space,scaled_shock_series,DynareOptions.order);
if any(any(isnan(y_sim))) || any(any(isinf(y_sim)))
fval = Inf;
info(1)=180;
info(4) = 0.1;
exit_flag = 0;
return
end
y_sim_after_burnin = y_sim(SMMinfo.varsindex,end-DynareOptions.smm.long:end)';
autolag=max(DynareOptions.smm.autolag);
if DynareOptions.smm.centeredmoments
y_sim_after_burnin=bsxfun(@minus,y_sim_after_burnin,mean(y_sim_after_burnin,1));
end
[modelMoments, E_y, E_yy, autoE_yy] = moments_GMM_SMM_Data(y_sim_after_burnin,DynareOptions);
% write centered and uncentered simulated moments to results
DynareResults.smm.unconditionalmoments.E_y=E_y;
DynareResults.smm.unconditionalmoments.E_yy=E_yy;
DynareResults.smm.unconditionalmoments.autoE_yy=autoE_yy;
DynareResults.smm.unconditionalmoments.Var_y=E_yy-E_y*E_y';
DynareResults.smm.unconditionalmoments.Cov_y=autoE_yy-repmat(E_y*E_y',[1 1 autolag]);
%------------------------------------------------------------------------------
% 4. Compute quadratic target function using weighting matrix W
%------------------------------------------------------------------------------
moments_difference = DynareResults.smm.datamoments.momentstomatch-modelMoments;
fval = moments_difference'*DynareResults.smm.W*moments_difference;
if DynareOptions.smm.penalized_estimator
fval=fval+(xparam1-SMMinfo.p1)'/diag(SMMinfo.p2)*(xparam1-SMMinfo.p1);
end
end
\ No newline at end of file
function [dataMoments,m_data] = method_of_moments_datamoments(data, matched_moments, index_moms_dr)
% [dataMoments E_y E_yy autoE_yy m_data]= moments_GMM_SMM_Data(data,DynareOptions)
% This function computes the following empirical moments from data
% - E[y]
% - E[y(i)*y(j)] for i=1:ny and j=i,ny
% - E[y(i)_t*y(i)_t-k] for i=1:ny and k=1,2,...autoLoag
%
% INPUTS
% data [T x ny] matrix data set
% numMom scalar number of selected moments
% DynareOptions Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
%
% OUTPUTS
% dataMoments [numMom x 1] vector selected moments
% E_y [ny x 1] vector (un)centered first moments
% E_yy [ny x ny] matrix (un)centered contemporaneous second moments
% autoE_yy [ny x ny x nLags] matrix (un)centered lagged second moments
% m_data [T x numMom] matrix empirical moments at each point in time
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2013-17 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/>.
% Number of observations (T) and number of observables (ny)
[T,ny] = size(data);
dataMoments = [];
%% First moments
if nnz(index_moms_dr.E_y) > 0
E_y = sum(data,1)'/T;
dataMoments = [dataMoments; E_y(index_moms_dr.E_y)];
end
%% Second moments
if nnz(index_moms_dr.E_yy) > 0
E_yy = data'*data/T;
dataMoments = [dataMoments; E_yy(tril(index_moms_dr.E_yy))];
end
if nnz(index_moms_dr.E_yyt) > 0
E_yyt = zeros(ny,ny,size(index_moms_dr.E_yyt,3));
for t=1:size(index_moms_dr.E_yyt,3)
E_yyt(:,:,t) = data(1+t:T,:)'*data(1:T-t,:)/T;
end
dataMoments = [dataMoments; E_yyt(index_moms_dr.E_yyt)];
end
if nargout>1
m_data = zeros(T,length(dataMoments));
for jm = 1:length(dataMoments)
if all(matched_moments{jm,2}==0) %no lags
m_data(:,jm) = prod(data(:,matched_moments{jm,1}).^matched_moments{jm,3},2);
else
m_data(:,jm) = lagmatrix(data(:,matched_moments{jm,1}(1)),abs(matched_moments{jm,2}(1))).^matched_moments{jm,3}(1); %first entry
for jv = 2:length(matched_moments{jm,1})
m_data(:,jm) = m_data(:,jm).*lagmatrix(data(:,matched_moments{jm,1}(jv)),abs(matched_moments{jm,2}(jv))).^matched_moments{jm,3}(jv);
end
end
end
end
% % For Ey
% m_data(:,1:nnz(index_moms_dr.E_y))= data(:,index_moms_dr.E_y');
%
% % % For Eyy
% % for i=1:ny
% % for j=i:ny
% % idxMom = idxMom + 1;
% % m_data(:,idxMom) = data(:,i).*data(:,j);
% % end
% % end
%
% % For Eyy
% for idxMom = 1:numMom
% hFunc(:,idxMom) = prod(data.^idx_moments(idxMom,:),2);
% end
% hFunc = hFunc - repmat(modelMoments',T,1);
%
%
% % For autoEyy
% maxLag = size(index_moms_dr.E_yyt,3)
% autoEyy = zeros(ny,ny,maxLag);
% for i=1:maxLag
% for j=1:ny
% autoEyy(1+i:T,j,i) = data(1+i:T,j).*data(1:T-i,j);
% end
% end
% for i=1:numLags
% for j=1:ny
% idxMom = idxMom + 1;
% m_data(:,idxMom) = autoEyy(:,j,autoLagsIdx(1,i));
% end
% end
end %function end
function oo = method_of_moments_initial_checks(objective_function, xparam1, DynareDataset, M, estim_params, options_mom, matched_moments, mom_info, bounds, oo)
% function DynareResults = initial_moment estimation_checks(objective_function,xparam1,DynareDataset,DatasetInfo,Model,EstimatedParameters,DynareOptions,BoundsInfo,DynareResults)
% Checks data (complex values, moment evaluation, initial values, BK conditions,..)
%
% INPUTS
% objective_function [function handle] of the objective function
% xparam1: [vector] of parameters to be estimated
% DynareDataset: [dseries] object storing the dataset
% DataSetInfo: [structure] storing informations about the sample.
% Model: [structure] decribing the model
% EstimatedParameters [structure] characterizing parameters to be estimated
% DynareOptions [structure] describing the options
% BoundsInfo [structure] containing bounds
% DynareResults [structure] storing the results
%
% OUTPUTS
% DynareResults structure of temporary results
%
% SPECIAL REQUIREMENTS
% none
% 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/>.
%get maximum number of simultaneously observed variables for stochastic
%singularity check
if options_mom.order>1
if any(any(isnan(DynareDataset.data)))
error('initial_moment_estimation_checks: moment estimation does not support missing observations')
end
end
if estim_params.nvn || estim_params.ncn
error('initial_moment_estimation_checks: moment estimation does not support measurement error(s). Please specifiy them as a structural shock')
end
if ~isempty(options_mom.endogenous_prior_restrictions.irf) && ~isempty(options_mom.endogenous_prior_restrictions.moment)
error('initial_moment_estimation_checks: Endogenous prior restrictions are not supported.')
end
old_steady_params=M.params; %save initial parameters for check if steady state changes param values
% % check if steady state solves static model
[oo.steady_state, new_steady_params] = evaluate_steady_state(oo.steady_state,M,options_mom,oo,1);
if isfield(estim_params,'param_vals') && ~isempty(estim_params.param_vals)
%check whether steady state file changes estimated parameters
Model_par_varied=M; %store Model structure
Model_par_varied.params(estim_params.param_vals(:,1))=Model_par_varied.params(estim_params.param_vals(:,1))*1.01; %vary parameters
[~, new_steady_params_2] = evaluate_steady_state(oo.steady_state,Model_par_varied,options_mom,oo,1);
changed_par_indices=find((old_steady_params(estim_params.param_vals(:,1))-new_steady_params(estim_params.param_vals(:,1))) ...
| (Model_par_varied.params(estim_params.param_vals(:,1))-new_steady_params_2(estim_params.param_vals(:,1))));
if ~isempty(changed_par_indices)
fprintf('\nThe steady state file internally changed the values of the following estimated parameters:\n')
disp(char(M.param_names(estim_params.param_vals(changed_par_indices,1))))
fprintf('This will override the parameter values drawn from the optimizer and may lead to wrong results.\n')
fprintf('Check whether this is really intended.\n')
warning('initial_moment_estimation_checks: The steady state file internally changes the values of the estimated parameters.')
end
end
% display warning if some parameters are still NaN
test_for_deep_parameters_calibration(M);
% Evaluate the moment-function.
tic_id=tic;
[fval,info,exit_flag,moments_difference,modelMoments] = feval(objective_function, xparam1, DynareDataset, options_mom, M, estim_params, mom_info, bounds, oo);
elapsed_time=toc(tic_id);
if isnan(fval)
error('initial_moment_estimation_checks: The initial value of the target function is NaN')
elseif imag(fval)
error('initial_moment_estimation_checks: The initial value of the target function is complex')
end
if info(1) >