From efcf6bd9c09fa43bcf78fdb88b7b2d3a4dae5e05 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?=
<stephane.adjemian@univ-lemans.fr>
Date: Mon, 16 Jun 2014 17:41:59 +0200
Subject: [PATCH] Use dseries object in the estimation routines.
---
matlab/dsge_likelihood.m | 34 +--
matlab/dynare_estimation_1.m | 72 ++++---
matlab/dynare_estimation_init.m | 105 +--------
matlab/initial_estimation_checks.m | 13 +-
matlab/metropolis_hastings_initialization.m | 8 +-
matlab/random_walk_metropolis_hastings.m | 5 +-
matlab/random_walk_metropolis_hastings_core.m | 3 +-
matlab/utilities/dataset/makedataset.m | 203 ++++++++++++++++++
tests/fs2000/fs2000_dseries.mod | 2 +-
9 files changed, 279 insertions(+), 166 deletions(-)
create mode 100644 matlab/utilities/dataset/makedataset.m
diff --git a/matlab/dsge_likelihood.m b/matlab/dsge_likelihood.m
index e3ffc7ceac..1ef6952621 100644
--- a/matlab/dsge_likelihood.m
+++ b/matlab/dsge_likelihood.m
@@ -1,4 +1,4 @@
-function [fval,DLIK,Hess,exit_flag,ys,trend_coeff,info,Model,DynareOptions,BayesInfo,DynareResults] = dsge_likelihood(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults,derivatives_info)
+function [fval,DLIK,Hess,exit_flag,ys,trend_coeff,info,Model,DynareOptions,BayesInfo,DynareResults] = dsge_likelihood(xparam1,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults,derivatives_info)
% Evaluates the posterior kernel of a dsge model using the specified
% kalman_algo; the resulting posterior includes the 2*pi constant of the
% likelihood function
@@ -295,7 +295,7 @@ BayesInfo.mf = BayesInfo.mf1;
% Define the constant vector of the measurement equation.
if DynareOptions.noconstant
- constant = zeros(DynareDataset.info.nvobs,1);
+ constant = zeros(DynareDataset.vobs,1);
else
if DynareOptions.loglinear
constant = log(SteadyState(BayesInfo.mfys));
@@ -306,28 +306,28 @@ end
% Define the deterministic linear trend of the measurement equation.
if BayesInfo.with_trend
- trend_coeff = zeros(DynareDataset.info.nvobs,1);
+ trend_coeff = zeros(DynareDataset.vobs,1);
t = DynareOptions.trend_coeffs;
for i=1:length(t)
if ~isempty(t{i})
trend_coeff(i) = evalin('base',t{i});
end
end
- trend = repmat(constant,1,DynareDataset.info.ntobs)+trend_coeff*[1:DynareDataset.info.ntobs];
+ trend = repmat(constant,1,DynareDataset.nobs)+trend_coeff*[1:DynareDataset.nobs];
else
- trend = repmat(constant,1,DynareDataset.info.ntobs);
+ trend = repmat(constant,1,DynareDataset.nobs);
end
% Get needed informations for kalman filter routines.
start = DynareOptions.presample+1;
-Z = BayesInfo.mf; % old mf
-no_missing_data_flag = ~DynareDataset.missing.state;
-mm = length(T); % old np
-pp = DynareDataset.info.nvobs;
+Z = BayesInfo.mf;
+no_missing_data_flag = ~DatasetInfo.missing.state;
+mm = length(T);
+pp = DynareDataset.vobs;
rr = length(Q);
kalman_tol = DynareOptions.kalman_tol;
riccati_tol = DynareOptions.riccati_tol;
-Y = DynareDataset.data-trend;
+Y = transpose(DynareDataset.data)-trend;
%------------------------------------------------------------------------------
% 3. Initial condition of the Kalman filter
@@ -406,7 +406,7 @@ switch DynareOptions.lik_init
kalman_tol, riccati_tol, DynareOptions.presample, ...
T,R,Q,H,Z,mm,pp,rr);
else
- [dLIK,dlik,a,Pstar] = missing_observations_kalman_filter_d(DynareDataset.missing.aindex,DynareDataset.missing.number_of_observations,DynareDataset.missing.no_more_missing_observations, ...
+ [dLIK,dlik,a,Pstar] = missing_observations_kalman_filter_d(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations, ...
Y, 1, size(Y,2), ...
zeros(mm,1), Pinf, Pstar, ...
kalman_tol, riccati_tol, DynareOptions.presample, ...
@@ -441,9 +441,9 @@ switch DynareOptions.lik_init
% no need to test again for correlation elements
correlated_errors_have_been_checked = 1;
- [dLIK,dlik,a,Pstar] = univariate_kalman_filter_d(DynareDataset.missing.aindex,...
- DynareDataset.missing.number_of_observations,...
- DynareDataset.missing.no_more_missing_observations, ...
+ [dLIK,dlik,a,Pstar] = univariate_kalman_filter_d(DatasetInfo.missing.aindex,...
+ DatasetInfo.missing.number_of_observations,...
+ DatasetInfo.missing.no_more_missing_observations, ...
Y, 1, size(Y,2), ...
zeros(mmm,1), Pinf, Pstar, ...
kalman_tol, riccati_tol, DynareOptions.presample, ...
@@ -655,10 +655,10 @@ if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
end
else
if 0 %DynareOptions.block
- [err, LIK,lik] = block_kalman_filter(DynareDataset.missing.aindex,DynareDataset.missing.number_of_observations,DynareDataset.missing.no_more_missing_observations,...
+ [err, LIK,lik] = block_kalman_filter(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations,...
T,R,Q,H,Pstar,Y,start,Z,kalman_tol,riccati_tol, Model.nz_state_var, Model.n_diag, Model.nobs_non_statevar);
else
- [LIK,lik] = missing_observations_kalman_filter(DynareDataset.missing.aindex,DynareDataset.missing.number_of_observations,DynareDataset.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ...
+ [LIK,lik] = missing_observations_kalman_filter(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ...
a, Pstar, ...
kalman_tol, DynareOptions.riccati_tol, ...
DynareOptions.presample, ...
@@ -733,7 +733,7 @@ if (kalman_algo==2) || (kalman_algo==4)
end
end
- [LIK, lik] = univariate_kalman_filter(DynareDataset.missing.aindex,DynareDataset.missing.number_of_observations,DynareDataset.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ...
+ [LIK, lik] = univariate_kalman_filter(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ...
a,Pstar, ...
DynareOptions.kalman_tol, ...
DynareOptions.riccati_tol, ...
diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index d617bfbf56..861af4d58d 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -29,7 +29,7 @@ function dynare_estimation_1(var_list_,dname)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
-global M_ options_ oo_ estim_params_ bayestopt_ dataset_
+global M_ options_ oo_ estim_params_ bayestopt_ dataset_ dataset_info
% Set particle filter flag.
if options_.order > 1
@@ -78,7 +78,8 @@ else
objective_function = str2func('dsge_var_likelihood');
end
-[dataset_,xparam1, hh, M_, options_, oo_, estim_params_,bayestopt_] = dynare_estimation_init(var_list_, dname, [], M_, options_, oo_, estim_params_, bayestopt_);
+[dataset_,dataset_info, xparam1, hh, M_, options_, oo_, estim_params_, bayestopt_] = ...
+ dynare_estimation_init(var_list_, dname, [], M_, options_, oo_, estim_params_, bayestopt_);
if options_.dsge_var
check_dsge_var_model(M_, estim_params_, bayestopt_);
@@ -114,12 +115,12 @@ if estim_params_.ncx || any(nnz(tril(M_.Correlation_matrix,-1))) || isfield(esti
end
data = dataset_.data;
-rawdata = dataset_.rawdata;
-data_index = dataset_.missing.aindex;
-missing_value = dataset_.missing.state;
+rawdata = dataset_.data;
+data_index = dataset_info.missing.aindex;
+missing_value = dataset_info.missing.state;
% Set number of observations
-gend = options_.nobs;
+gend = dataset_.nobs;
% Set the number of observed variables.
n_varobs = length(options_.varobs);
% Get the number of parameters to be estimated.
@@ -156,7 +157,7 @@ end
% compute sample moments if needed (bvar-dsge)
if options_.dsge_var
- if dataset_.missing.state
+ if dataset_info.missing.state
error('I cannot estimate a DSGE-VAR model with missing observations!')
end
if options_.noconstant
@@ -176,7 +177,7 @@ end
%% perform initial estimation checks;
try
- oo_ = initial_estimation_checks(objective_function,xparam1,dataset_,M_,estim_params_,options_,bayestopt_,oo_);
+ oo_ = initial_estimation_checks(objective_function,xparam1,dataset_,dataset_info,M_,estim_params_,options_,bayestopt_,oo_);
catch % if check fails, provide info on using calibration if present
e = lasterror();
if full_calibration_detected %calibrated model present and no explicit starting values
@@ -246,7 +247,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
optim_options = optimset(optim_options,'GradObj','on','TolX',1e-7);
end
[xparam1,fval,exitflag,output,lamdba,grad,hessian_fmincon] = ...
- fmincon(objective_function,xparam1,[],[],[],[],lb,ub,[],optim_options,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+ fmincon(objective_function,xparam1,[],[],[],[],lb,ub,[],optim_options,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
case 2
error('ESTIMATION: mode_compute=2 option (Lester Ingber''s Adaptive Simulated Annealing) is no longer available')
case 3
@@ -264,10 +265,10 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
optim_options = optimset(optim_options,'GradObj','on');
end
if ~isoctave
- [xparam1,fval,exitflag] = fminunc(objective_function,xparam1,optim_options,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+ [xparam1,fval,exitflag] = fminunc(objective_function,xparam1,optim_options,dataset_,dataset_info_,options_,M_,estim_params_,bayestopt_,oo_);
else
% Under Octave, use a wrapper, since fminunc() does not have a 4th arg
- func = @(x) objective_function(x, dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+ func = @(x) objective_function(x, dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
[xparam1,fval,exitflag] = fminunc(func,xparam1,optim_options);
end
case 4
@@ -306,7 +307,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
end
% Call csminwell.
[fval,xparam1,grad,hessian_csminwel,itct,fcount,retcodehat] = ...
- csminwel1(objective_function, xparam1, H0, analytic_grad, crit, nit, numgrad, epsilon, dataset_, options_, M_, estim_params_, bayestopt_, oo_);
+ csminwel1(objective_function, xparam1, H0, analytic_grad, crit, nit, numgrad, epsilon, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_);
% Disp value at the mode.
disp(sprintf('Objective function at mode: %f',fval))
case 5
@@ -334,7 +335,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
else
nit=1000;
end
- [xparam1,hh,gg,fval,invhess] = newrat(objective_function,xparam1,analytic_grad,crit,nit,flag,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+ [xparam1,hh,gg,fval,invhess] = newrat(objective_function,xparam1,analytic_grad,crit,nit,flag,dataset_, dataset_info, options_,M_,estim_params_,bayestopt_,oo_);
if options_.analytic_derivation,
options_.analytic_derivation = ana_deriv;
end
@@ -384,7 +385,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
end
end
% Evaluate the objective function.
- fval = feval(objective_function,xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+ fval = feval(objective_function,xparam1,dataset_, dataset_info, options_,M_,estim_params_,bayestopt_,oo_);
OldMode = fval;
if ~exist('MeanPar','var')
MeanPar = xparam1;
@@ -416,8 +417,8 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
flag = 'LastCall';
end
[xparam1,PostVar,Scale,PostMean] = ...
- gmhmaxlik(objective_function,xparam1,[lb ub],gmhmaxlikOptions,Scale,flag,MeanPar,CovJump,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
- fval = feval(objective_function,xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+ gmhmaxlik(objective_function,xparam1,[lb ub],gmhmaxlikOptions,Scale,flag,MeanPar,CovJump,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
+ fval = feval(objective_function,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
options_.mh_jscale = Scale;
mouvement = max(max(abs(PostVar-OldPostVar)));
skipline()
@@ -436,11 +437,11 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
end
[xparam1,PostVar,Scale,PostMean] = ...
gmhmaxlik(objective_function,xparam1,[lb ub],...
- gmhmaxlikOptions,Scale,flag,PostMean,PostVar,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
- fval = feval(objective_function,xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+ gmhmaxlikOptions,Scale,flag,PostMean,PostVar,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
+ fval = feval(objective_function,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
options_.mh_jscale = Scale;
mouvement = max(max(abs(PostVar-OldPostVar)));
- fval = feval(objective_function,xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+ fval = feval(objective_function,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
skipline()
disp('========================================================== ')
disp([' Change in the covariance matrix = ' num2str(mouvement) '.'])
@@ -473,7 +474,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
if isfield(options_,'optim_opt')
eval(['optim_options = optimset(optim_options,' options_.optim_opt ');']);
end
- [xparam1,fval,exitflag] = fminsearch(objective_function,xparam1,optim_options,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+ [xparam1,fval,exitflag] = fminsearch(objective_function,xparam1,optim_options,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
case 8
% Dynare implementation of the simplex algorithm.
simplexOptions = options_.simplex;
@@ -498,7 +499,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
end
end
end
- [xparam1,fval,exitflag] = simplex_optimization_routine(objective_function,xparam1,simplexOptions,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+ [xparam1,fval,exitflag] = simplex_optimization_routine(objective_function,xparam1,simplexOptions,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
case 9
% Set defaults
H0 = 1e-4*ones(nx,1);
@@ -523,7 +524,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
end
warning('off','CMAES:NonfinitenessRange');
warning('off','CMAES:InitialSigma');
- [x, fval, COUNTEVAL, STOPFLAG, OUT, BESTEVER] = cmaes(func2str(objective_function),xparam1,H0,cmaesOptions,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+ [x, fval, COUNTEVAL, STOPFLAG, OUT, BESTEVER] = cmaes(func2str(objective_function),xparam1,H0,cmaesOptions,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
xparam1=BESTEVER.x;
disp(sprintf('\n Objective function at mode: %f',fval))
case 10
@@ -554,16 +555,16 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
end
simpsaOptionsList = options2cell(simpsaOptions);
simpsaOptions = simpsaset(simpsaOptionsList{:});
- [xparam1, fval, exitflag] = simpsa(func2str(objective_function),xparam1,lb,ub,simpsaOptions,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+ [xparam1, fval, exitflag] = simpsa(func2str(objective_function),xparam1,lb,ub,simpsaOptions,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
case 11
options_.cova_compute = 0 ;
- [xparam1,stdh,lb_95,ub_95,med_param] = online_auxiliary_filter(xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_) ;
+ [xparam1,stdh,lb_95,ub_95,med_param] = online_auxiliary_filter(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_) ;
case 101
myoptions=soptions;
myoptions(2)=1e-6; %accuracy of argument
myoptions(3)=1e-6; %accuracy of function (see Solvopt p.29)
myoptions(5)= 1.0;
- [xparam1,fval]=solvopt(xparam1,objective_function,[],myoptions,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+ [xparam1,fval]=solvopt(xparam1,objective_function,[],myoptions,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
case 102
%simulating annealing
% LB=zeros(size(xparam1))-20;
@@ -600,10 +601,10 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
disp(['c vector ' num2str(c')]);
[xparam1, fval, nacc, nfcnev, nobds, ier, t, vm] = sa(objective_function,xparam1,maxy,rt_,epsilon,ns,nt ...
- ,neps,maxevl,LB,UB,c,idisp ,t,vm,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+ ,neps,maxevl,LB,UB,c,idisp ,t,vm,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
otherwise
if ischar(options_.mode_compute)
- [xparam1, fval, retcode ] = feval(options_.mode_compute,objective_function,xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+ [xparam1, fval, retcode ] = feval(options_.mode_compute,objective_function,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
else
error(['dynare_estimation:: mode_compute = ' int2str(options_.mode_compute) ' option is unknown!'])
end
@@ -614,11 +615,11 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
ana_deriv = options_.analytic_derivation;
options_.analytic_derivation = 2;
[junk1, junk2, hh] = feval(objective_function,xparam1, ...
- dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+ dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
options_.analytic_derivation = ana_deriv;
else
hh = reshape(hessian(objective_function,xparam1, ...
- options_.gstep,dataset_,options_,M_,estim_params_,bayestopt_,oo_),nx,nx);
+ options_.gstep,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_),nx,nx);
end
end
end
@@ -698,7 +699,7 @@ end
if options_.mode_check.status == 1 && ~options_.mh_posterior_mode_estimation
ana_deriv = options_.analytic_derivation;
options_.analytic_derivation = 0;
- mode_check(objective_function,xparam1,hh,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+ mode_check(objective_function,xparam1,hh,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
options_.analytic_derivation = ana_deriv;
end
@@ -732,7 +733,7 @@ if any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation
estim_params_nbr = size(xparam1,1);
scale_factor = -sum(log10(diag(invhess)));
log_det_invhess = -estim_params_nbr*log(scale_factor)+log(det(scale_factor*invhess));
- likelihood = feval(objective_function,xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+ likelihood = feval(objective_function,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
oo_.MarginalDensity.LaplaceApproximation = .5*estim_params_nbr*log(2*pi) + .5*log_det_invhess - likelihood;
skipline()
disp(sprintf('Log data density [Laplace approximation] is %f.',oo_.MarginalDensity.LaplaceApproximation))
@@ -779,7 +780,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
ana_deriv = options_.analytic_derivation;
options_.analytic_derivation = 0;
if options_.cova_compute
- feval(options_.posterior_sampling_method,objective_function,options_.proposal_distribution,xparam1,invhess,bounds,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+ feval(options_.posterior_sampling_method,objective_function,options_.proposal_distribution,xparam1,invhess,bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
else
error('I Cannot start the MCMC because the Hessian of the posterior kernel at the mode was not computed.')
end
@@ -842,7 +843,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
> 0) && options_.load_mh_file)) ...
|| ~options_.smoother ) && options_.partial_information == 0 % to be fixed
%% ML estimation, or posterior mode without metropolis-hastings or metropolis without bayesian smooth variable
- [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp] = DsgeSmoother(xparam1,dataset_.info.ntobs,dataset_.data,dataset_.missing.aindex,dataset_.missing.state);
+ [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp] = DsgeSmoother(xparam1,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state);
oo_.Smoother.SteadyState = ys;
oo_.Smoother.TrendCoeffs = trend_coeff;
oo_.Smoother.Variance = P;
@@ -949,8 +950,11 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
%% Smooth observational errors...
%%
if options_.prefilter == 1
- yf = atT(bayestopt_.mf,:)+repmat(dataset_.descriptive.mean',1,gend);
+ yf = atT(bayestopt_.mf,:)+repmat(dataset_info.descriptive.mean',1,gend);
elseif options_.loglinear == 1
+ gend
+ bayestopt_.mfys
+ ys
yf = atT(bayestopt_.mf,:)+repmat(log(ys(bayestopt_.mfys)),1,gend)+...
trend_coeff*[1:gend];
else
diff --git a/matlab/dynare_estimation_init.m b/matlab/dynare_estimation_init.m
index a9e9e18589..90189e7c85 100644
--- a/matlab/dynare_estimation_init.m
+++ b/matlab/dynare_estimation_init.m
@@ -1,4 +1,4 @@
-function [dataset_, xparam1, hh, M_, options_, oo_, estim_params_,bayestopt_, fake] = dynare_estimation_init(var_list_, dname, gsa_flag, M_, options_, oo_, estim_params_, bayestopt_)
+function [dataset_, dataset_info, xparam1, hh, M_, options_, oo_, estim_params_,bayestopt_, fake] = dynare_estimation_init(var_list_, dname, gsa_flag, M_, options_, oo_, estim_params_, bayestopt_)
% function dynare_estimation_init(var_list_, gsa_flag)
% preforms initialization tasks before estimation or
@@ -447,106 +447,8 @@ end
k = find(isnan(bayestopt_.jscale));
bayestopt_.jscale(k) = options_.mh_jscale;
-% Test if the dataset is declared.
-if isempty(options_.datafile) && isempty(options_.dataset)
- if gsa_flag
- dataset_ = [];
- return
- else
- error('datafile option is missing')
- end
-end
-
-if ~isfield(options_,'nobs')
- options_.nobs = [];
-end
-
-if isempty(options_.datafile) && ~isempty(options_.dataset.file)
- datafile = options_.dataset.file;
-elseif ~isempty(options_.datafile) && isempty(options_.dataset.file)
- datafile = options_.datafile;
-elseif isempty(options_.datafile) && ~isempty(options_.dataset.file)
- error('You cannot use simultaneously the data command and the datafile option (in the estimation command)!')
-else
- error('You have to specify the datafile!')
-end
-
-
-% Check extension.
-allowed_extensions = {'m','mat','csv','xls','xlsx'};
-datafile_extension = get_file_extension(datafile);
-if isempty(datafile_extension)
- available_extensions = {}; j = 1;
- for i=1:length(allowed_extensions)
- if exist([datafile '.' allowed_extensions{i}])
- available_extensions(j) = {allowed_extensions{i}};
- j = j+1;
- end
- end
- if isempty(available_extensions)
- error(['I can''t find a datafile (with allowed extension)!'])
- end
- if length(available_extensions)>1
- error(sprintf(['You did not specify an extension for the datafile, but more than one candidate ' ...
- 'are available in the designed folder!\nPlease, add an extension to the datafile ' ...
- '(m, mat, csv, xls or xlsx are legal extensions).']));
- end
- datafile = [datafile '.' available_extensions{1}];
-end
-
-% Load the data in a dseries object.
-dataset = dseries(datafile);
-
-% Select a subset of the variables.
-dataset = dataset{options_.varobs{:}};
-
-% Apply log function if needed.
-if options_.loglinear && ~options_.logdata
- dataset = dataset.log();
-end
-
-% Test if an initial period (different from its default value) is explicitely defined in the datafile.
-if isequal(dataset.init, dates(1))
- dataset_default_initial_period = 1;
-else
- dataset_default_initial_period = 0;
-end
-
-% Test if an initial period (different from its default value) is explicitely defined in the mod file with the set_time command.
-if isequal(options_.initial_period, dates(1))
- set_time_default_initial_period = 1;
-else
- set_time_default_initial_period = 0;
-end
-
-if ~set_time_default_initial_period && dataset_default_initial_period
- % Overwrite the initial period in dataset (it was set to default).
- % Note that the update of freq and time members is auto-magically
- % done by dseries::subsasgn overload method.
- dataset.init = options_.initial_period;
-end
-
-if set_time_default_initial_period && ~dataset_default_initial_period
- % Overwrite the global initial period defined by set_time (it was set to default).
- options_.initial_period = dataset.init;
-end
-
-if ~set_time_default_initial_period && ~dataset_default_initial_period
- % Check if dataset.init and options_.initial_period are identical.
- if ~isequal(options_.initial_period, dataset.init)
- error('dynare_estimation_init:: The date as defined by the set_time command is not consistent with the initial period in the database!')
- end
-end
-
-
-if options_.initial_period>options_.dataset.first_obs
- skipline()
- error(['First observation (' format(options_.dataset.first_obs) ') must be greater than the initial period (' format(options_.initial_period) ') as set by the set_time command']);
-end
-
-dataset_ = initialize_dataset(options_.datafile,options_.varobs,options_.first_obs,options_.nobs,logged_data_flag,options_.prefilter,xls);
-
-options_.nobs = dataset_.info.ntobs;
+% Build the dataset
+[dataset_, dataset_info] = makedataset(options_);
% setting steadystate_check_flag option
if options_.diffuse_filter
@@ -554,6 +456,7 @@ if options_.diffuse_filter
else
steadystate_check_flag = 1;
end
+
% If the steady state of the observed variables is non zero, set noconstant equal 0 ()
M = M_;
nvx = estim_params_.nvx;
diff --git a/matlab/initial_estimation_checks.m b/matlab/initial_estimation_checks.m
index 9dee5461e9..c4cbf5ab96 100644
--- a/matlab/initial_estimation_checks.m
+++ b/matlab/initial_estimation_checks.m
@@ -1,11 +1,12 @@
-function DynareResults = initial_estimation_checks(objective_function,xparam1,DynareDataset,Model,EstimatedParameters,DynareOptions,BayesInfo,DynareResults)
+function DynareResults = initial_estimation_checks(objective_function,xparam1,DynareDataset,DatasetInfo,Model,EstimatedParameters,DynareOptions,BayesInfo,DynareResults)
% function initial_estimation_checks(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations)
% Checks data (complex values, ML evaluation, initial values, BK conditions,..)
%
% INPUTS
% objective_function [function handle] of the objective function
% xparam1: [vector] of parameters to be estimated
-% DynareDataset: [structure] storing the dataset information
+% 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
@@ -35,11 +36,11 @@ function DynareResults = initial_estimation_checks(objective_function,xparam1,Dy
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
-if DynareDataset.info.nvobs>Model.exo_nbr+EstimatedParameters.nvn
+if DynareDataset.vobs>Model.exo_nbr+EstimatedParameters.nvn
error(['initial_estimation_checks:: Estimation can''t take place because there are less declared shocks than observed variables!'])
end
-
-if DynareDataset.info.nvobs>length(find(diag(Model.Sigma_e)))+EstimatedParameters.nvn
+
+if DynareDataset.vobs>length(find(diag(Model.Sigma_e)))+EstimatedParameters.nvn
error(['initial_estimation_checks:: Estimation can''t take place because too many shocks have been calibrated with a zero variance!'])
end
@@ -72,7 +73,7 @@ end
% Evaluate the likelihood.
ana_deriv = DynareOptions.analytic_derivation;
DynareOptions.analytic_derivation=0;
-[fval,junk1,junk2,a,b,c,d] = feval(objective_function,xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
+[fval,junk1,junk2,a,b,c,d] = feval(objective_function,xparam1,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
DynareOptions.analytic_derivation=ana_deriv;
if DynareOptions.dsge_var || strcmp(func2str(objective_function),'non_linear_dsge_likelihood')
diff --git a/matlab/metropolis_hastings_initialization.m b/matlab/metropolis_hastings_initialization.m
index 8c6ed5a9f6..b475bfcc5e 100644
--- a/matlab/metropolis_hastings_initialization.m
+++ b/matlab/metropolis_hastings_initialization.m
@@ -1,7 +1,7 @@
function [ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = ...
- metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,options_,M_,estim_params_,bayestopt_,oo_)
+ metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
%function [ ix2, ilogpo2, ModelName, MhDirectoryName, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] =
-% metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds, dataset_,options_,M_,estim_params_,bayestopt_,oo_)
+% metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds, dataset_,dataset_info,,options_,M_,estim_params_,bayestopt_,oo_)
% Metropolis-Hastings initialization.
%
% INPUTS
@@ -110,7 +110,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
candidate = rand_multivariate_normal( transpose(xparam1), d * options_.mh_init_scale, npar);
if all(candidate(:) > mh_bounds(:,1)) && all(candidate(:) < mh_bounds(:,2))
ix2(j,:) = candidate;
- ilogpo2(j) = - feval(TargetFun,ix2(j,:)',dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+ ilogpo2(j) = - feval(TargetFun,ix2(j,:)',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
if ~isfinite(ilogpo2(j)) % if returned log-density is
% Inf or Nan (penalized value)
validate = 0;
@@ -152,7 +152,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
candidate = transpose(xparam1(:));%
if all(candidate(:) > mh_bounds(:,1)) && all(candidate(:) < mh_bounds(:,2))
ix2 = candidate;
- ilogpo2 = - feval(TargetFun,ix2',dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+ ilogpo2 = - feval(TargetFun,ix2',dataset_,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
disp('Estimation::mcmc: Initialization at the posterior mode.')
skipline()
fprintf(fidlog,[' Blck ' int2str(1) 'params:\n']);
diff --git a/matlab/random_walk_metropolis_hastings.m b/matlab/random_walk_metropolis_hastings.m
index f524b44c13..c2a49e7a29 100644
--- a/matlab/random_walk_metropolis_hastings.m
+++ b/matlab/random_walk_metropolis_hastings.m
@@ -1,4 +1,4 @@
-function random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,dataset_,options_,M_,estim_params_,bayestopt_,oo_)
+function random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
%function record=random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,dataset_,options_,M_,estim_params_,bayestopt_,oo_)
% Random walk Metropolis-Hastings algorithm.
%
@@ -60,7 +60,7 @@ objective_function_penalty_base = Inf;
% Initialization of the random walk metropolis-hastings chains.
[ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = ...
- metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+ metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
InitSizeArray = min([repmat(MAX_nruns,nblck,1) fline+nruns-1],[],2);
@@ -101,6 +101,7 @@ localVars = struct('TargetFun', TargetFun, ...
'InitSizeArray',InitSizeArray, ...
'record', record, ...
'dataset_', dataset_, ...
+ 'dataset_info', dataset_info, ...
'options_', options_, ...
'M_',M_, ...
'bayestopt_', bayestopt_, ...
diff --git a/matlab/random_walk_metropolis_hastings_core.m b/matlab/random_walk_metropolis_hastings_core.m
index 481e92b77a..62e1728eb9 100644
--- a/matlab/random_walk_metropolis_hastings_core.m
+++ b/matlab/random_walk_metropolis_hastings_core.m
@@ -88,6 +88,7 @@ d=myinputs.d;
InitSizeArray=myinputs.InitSizeArray;
record=myinputs.record;
dataset_ = myinputs.dataset_;
+dataset_info = myinputs.dataset_info;
bayestopt_ = myinputs.bayestopt_;
estim_params_ = myinputs.estim_params_;
options_ = myinputs.options_;
@@ -164,7 +165,7 @@ for b = fblck:nblck,
par = feval(ProposalFun, ix2(b,:), proposal_covariance_Cholesky_decomposition, n);
if all( par(:) > mh_bounds(:,1) ) && all( par(:) < mh_bounds(:,2) )
try
- logpost = - feval(TargetFun, par(:),dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+ logpost = - feval(TargetFun, par(:),dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
catch
logpost = -inf;
end
diff --git a/matlab/utilities/dataset/makedataset.m b/matlab/utilities/dataset/makedataset.m
new file mode 100644
index 0000000000..52064d441b
--- /dev/null
+++ b/matlab/utilities/dataset/makedataset.m
@@ -0,0 +1,203 @@
+function [DynareDataset, DatasetInfo] = makedataset(DynareOptions)
+
+% Initialize a dataset as a dseries object.
+%
+%
+% INPUTS
+% ======
+%
+% DynareOptions [struct] Structure of options built by Dynare's preprocessor.
+%
+%
+% OUTPUTS
+% =======
+%
+% DynareDataset [dseries] The dataset.
+% DatasetInfo [struct] Various informations about the dataset (descriptive statistics and missing observations).
+%
+% EXAMPLE
+% =======
+%
+% [dataset_, dataset_info] = makedataset(options_) ;
+%
+%
+% See also dynare_estimation_init
+
+if isempty(DynareOptions.datafile) && isempty(DynareOptions.dataset.file)
+ if gsa_flag
+ DynareDataset = dseries();
+ return
+ else
+ error('datafile option is missing')
+ end
+end
+
+if isempty(DynareOptions.datafile) && ~isempty(DynareOptions.dataset.file)
+ datafile = DynareOptions.dataset.file;
+ newdatainterface = 1;
+elseif ~isempty(DynareOptions.datafile) && isempty(DynareOptions.dataset.file)
+ datafile = DynareOptions.datafile;
+ newdatainterface = 0;
+elseif isempty(DynareOptions.datafile) && ~isempty(DynareOptions.dataset.file)
+ error('You cannot use simultaneously the data command and the datafile option (in the estimation command)!')
+else
+ error('You have to specify the datafile!')
+end
+
+% Check extension.
+allowed_extensions = {'m','mat','csv','xls','xlsx'};
+datafile_extension = get_file_extension(datafile);
+if isempty(datafile_extension)
+ available_extensions = {}; j = 1;
+ for i=1:length(allowed_extensions)
+ if exist([datafile '.' allowed_extensions{i}])
+ available_extensions(j) = {allowed_extensions{i}};
+ j = j+1;
+ end
+ end
+ if isempty(available_extensions)
+ error(['I can''t find a datafile (with allowed extension)!'])
+ end
+ if length(available_extensions)>1
+ error(sprintf(['You did not specify an extension for the datafile, but more than one candidate ' ...
+ 'are available in the designed folder!\nPlease, add an extension to the datafile ' ...
+ '(m, mat, csv, xls or xlsx are legal extensions).']));
+ end
+ datafile = [datafile '.' available_extensions{1}];
+end
+
+% Load the data in a dseries object.
+DynareDataset = dseries(datafile);
+
+% Select a subset of the variables.
+DynareDataset = DynareDataset{DynareOptions.varobs{:}};
+
+% Apply log function if needed.
+if DynareOptions.loglinear && ~DynareOptions.logdata
+ DynareDataset = DynareDataset.log();
+end
+
+% Test if an initial period (different from its default value) is explicitely defined in the datafile.
+if isequal(DynareDataset.init, dates(1,1))
+ dataset_default_initial_period = 1;
+else
+ dataset_default_initial_period = 0;
+end
+
+% Test if an initial period (different from its default value) is explicitely defined in the mod file with the set_time command.
+if isequal(DynareOptions.initial_period, dates(1,1))
+ set_time_default_initial_period = 1;
+else
+ set_time_default_initial_period = 0;
+end
+
+if ~set_time_default_initial_period && dataset_default_initial_period
+ % Overwrite the initial period in dataset (it was set to default).
+ % Note that the updates of freq and time members are auto-magically
+ % done by dseries::subsasgn overloaded method.
+ DynareDataset.init = DynareOptions.initial_period;
+end
+
+if set_time_default_initial_period && ~dataset_default_initial_period
+ % Overwrite the global initial period defined by set_time (it was set to default).
+ DynareOptions.initial_period = DynareDataset.init;
+end
+
+if ~set_time_default_initial_period && ~dataset_default_initial_period
+ % Check if dataset.init and options_.initial_period are identical.
+ if ~isequal(DynareOptions.initial_period, DynareDataset.init)
+ error('The date as defined by the set_time command is not consistent with the initial period in the database!')
+ end
+end
+
+% Set firstobs, lastobs and nobs
+if newdatainterface
+ if isempty(DynareOptions.dataset.firstobs)
+ % first_obs option was not used in the data command.
+ firstobs = DynareDataset.init;
+ else
+ firstobs = DynareOptions.dataset.firstobs;
+ end
+ if isnan(DynareOptions.dataset.nobs)
+ % nobs option was not used in the data command.
+ if isempty(DynareOptions.dataset.lastobs)
+ % last_obs option was not used in the data command.
+ nobs = DynareDataset.nobs;
+ lastobs = DynareDataset.dates(end);
+ else
+ lastobs = DynareOptions.dataset.lastobs;
+ nobs = lastobs-firstobs+1;
+ end
+ else
+ nobs = DynareOptions.dataset.nobs;
+ if isempty(DynareOptions.dataset.lastobs)
+ % last_obs option was not used in the data command.
+ lastobs = firstobs+(nobs-1);
+ else
+ % last_obs and nobs were used in the data command. Check that they are consistent (with firstobs).
+ if ~isequal(lastobs,firstobs+(nobs-1))
+ error(sprintf('Options last_obs (%s), first_obs (%s) and nobs (%s) are not consistent!',char(lastobs),char(firstobs),num2str(nobs)));
+ end
+ end
+ end
+else
+ if isnan(DynareOptions.first_obs)
+ firstobs = DynareDataset.init;
+ else
+ firstobs = DynareDataset.dates(DynareOptions.first_obs);
+ end
+ if isnan(DynareOptions.nobs)
+ lastobs = DynareDataset.dates(end);
+ nobs = lastobs-firstobs+1;
+ else
+ nobs = DynareOptions.nobs;
+ lastobs = firstobs+(nobs-1);
+ end
+end
+
+% Check that firstobs belongs to DynareDataset.dates
+if firstobs<DynareDataset.init
+ error(sprintf('first_obs (%s) cannot be less than the first date in the dataset (%s)!',char(firstobs),char(DynareDataset.init)))
+end
+
+% Check that lastobs belongs to DynareDataset.dates...
+if newdatainterface
+ if lastobs>DynareDataset.dates(end)
+ error(sprintf('last_obs (%s) cannot be greater than the last date in the dataset (%s)!',char(lastobs),char(DynareDataset.dates(end))))
+ end
+else
+ % ... or check that nobs is smaller than the number of observations in DynareDataset.
+ if nobs>DynareDataset.nobs
+ error(sprintf('nobs (%s) cannot be greater than the last date in the dataset (%s)!', num2str(nobs), num2str(nobs)))
+ end
+end
+
+% Select a subsample.
+DynareDataset = DynareDataset(firstobs:lastobs);
+
+% Initialize DatasetInfo structure.
+DatasetInfo = struct('missing', struct('state', NaN, 'aindex', [], 'vindex', [], 'number_of_observations', NaN, 'no_more_missing_observations', NaN), ...
+ 'descriptive', struct('mean', [], 'covariance', [], 'correlation', [], 'autocovariance', []));
+
+% Fill DatasetInfo.missing if some observations are missing
+DatasetInfo.missing.state = isanynan(DynareDataset.data);
+if DatasetInfo.missing.state
+ [DatasetInfo.missing.aindex, DatasetInfo.missing.number_of_observations, DatasetInfo.missing.no_more_missing_observations, DatasetInfo.missing.vindex] = ...
+ describe_missing_data(DynareDataset.data);
+else
+ DatasetInfo.missing.aindex = num2cell(transpose(repmat(1:DynareDataset.vobs,DynareDataset.nobs,1)),1);
+ DatasetInfo.missing.no_more_missing_observations = 1;
+end
+
+% Compute the empirical mean of the observed variables.
+DatasetInfo.descriptive.mean = nanmean(DynareDataset.data);
+
+% Compute the empirical covariance matrix of the observed variables.
+DatasetInfo.descriptive.covariance = nancovariance(DynareDataset.data);
+
+% Compute the empirical correlation matrix of the observed variables.
+normalization_matrix = diag(1./sqrt(diag(DatasetInfo.descriptive.covariance)));
+DatasetInfo.descriptive.correlation = normalization_matrix*DatasetInfo.descriptive.covariance*normalization_matrix;
+
+% Compute autocorrelation function.
+DatasetInfo.descriptive.autocovariance = nanautocovariance(DynareDataset.data, DynareOptions.ar);
\ No newline at end of file
diff --git a/tests/fs2000/fs2000_dseries.mod b/tests/fs2000/fs2000_dseries.mod
index c0d851721f..1210e7425c 100644
--- a/tests/fs2000/fs2000_dseries.mod
+++ b/tests/fs2000/fs2000_dseries.mod
@@ -72,6 +72,6 @@ varobs gp_obs gy_obs;
options_.solve_tolf = 1e-12;
-data(file=fsdat_simul,first_obs=1950Q3, nobs=192);
+data(file=fsdat_simul_dseries,first_obs=1950Q3, nobs=192);
estimation(order=1,loglinear,mh_replic=2000,mh_nblocks=2,mh_jscale=0.8);
--
GitLab