Adapted code for dsge-var models.
Stéphane Adjemian (Scylla)
parent
91d74fabb3
commit
b11f6e2505
|
|
@ -1,4 +1,4 @@
|
|||
function [fval,grad,hess,exit_flag,info,PHI,SIGMAu,iXX,prior] = dsge_var_likelihood(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
|
||||
function [fval,grad,hess,exit_flag,info,PHI,SIGMAu,iXX,prior] = dsge_var_likelihood(xparam1,DynareDataset,DynareInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
|
||||
% Evaluates the posterior kernel of the bvar-dsge model.
|
||||
%
|
||||
% INPUTS
|
||||
|
|
@ -55,7 +55,11 @@ end
|
|||
nx = EstimatedParameters.nvx + EstimatedParameters.np;
|
||||
|
||||
% Get the number of observed variables in the VAR model.
|
||||
NumberOfObservedVariables = DynareDataset.info.nvobs;
|
||||
NumberOfObservedVariables = DynareDataset.vobs;
|
||||
|
||||
% Get the number of observations.
|
||||
NumberOfObservations = DynareDataset.nobs;
|
||||
|
||||
|
||||
% Get the number of lags in the VAR model.
|
||||
NumberOfLags = DynareOptions.dsge_varlag;
|
||||
|
|
@ -110,8 +114,8 @@ Model.Sigma_e = Q;
|
|||
dsge_prior_weight = Model.params(dsge_prior_weight_idx);
|
||||
|
||||
% Is the dsge prior proper?
|
||||
if dsge_prior_weight<(NumberOfParameters+NumberOfObservedVariables)/DynareDataset.info.ntobs;
|
||||
fval = objective_function_penalty_base+abs(DynareDataset.info.ntobs*dsge_prior_weight-(NumberOfParameters+NumberOfObservedVariables));
|
||||
if dsge_prior_weight<(NumberOfParameters+NumberOfObservedVariables)/NumberOfObservations;
|
||||
fval = objective_function_penalty_base+abs(NumberOfObservations*dsge_prior_weight-(NumberOfParameters+NumberOfObservedVariables));
|
||||
exit_flag = 0;
|
||||
info = 51;
|
||||
return
|
||||
|
|
@ -197,9 +201,9 @@ assignin('base','GXX',GXX);
|
|||
assignin('base','GYX',GYX);
|
||||
|
||||
if ~isinf(dsge_prior_weight)% Evaluation of the likelihood of the dsge-var model when the dsge prior weight is finite.
|
||||
tmp0 = dsge_prior_weight*DynareDataset.info.ntobs*TheoreticalAutoCovarianceOfTheObservedVariables(:,:,1) + mYY ;
|
||||
tmp1 = dsge_prior_weight*DynareDataset.info.ntobs*GYX + mYX;
|
||||
tmp2 = inv(dsge_prior_weight*DynareDataset.info.ntobs*GXX+mXX);
|
||||
tmp0 = dsge_prior_weight*NumberOfObservations*TheoreticalAutoCovarianceOfTheObservedVariables(:,:,1) + mYY ;
|
||||
tmp1 = dsge_prior_weight*NumberOfObservations*GYX + mYX;
|
||||
tmp2 = inv(dsge_prior_weight*NumberOfObservations*GXX+mXX);
|
||||
SIGMAu = tmp0 - tmp1*tmp2*tmp1'; clear('tmp0');
|
||||
[SIGMAu_is_positive_definite, penalty] = ispd(SIGMAu);
|
||||
if ~SIGMAu_is_positive_definite
|
||||
|
|
@ -208,28 +212,28 @@ if ~isinf(dsge_prior_weight)% Evaluation of the likelihood of the dsge-var model
|
|||
exit_flag = 0;
|
||||
return;
|
||||
end
|
||||
SIGMAu = SIGMAu / (DynareDataset.info.ntobs*(1+dsge_prior_weight));
|
||||
SIGMAu = SIGMAu / (NumberOfObservations*(1+dsge_prior_weight));
|
||||
PHI = tmp2*tmp1'; clear('tmp1');
|
||||
prodlng1 = sum(gammaln(.5*((1+dsge_prior_weight)*DynareDataset.info.ntobs- ...
|
||||
prodlng1 = sum(gammaln(.5*((1+dsge_prior_weight)*NumberOfObservations- ...
|
||||
NumberOfObservedVariables*NumberOfLags ...
|
||||
+1-(1:NumberOfObservedVariables)')));
|
||||
prodlng2 = sum(gammaln(.5*(dsge_prior_weight*DynareDataset.info.ntobs- ...
|
||||
prodlng2 = sum(gammaln(.5*(dsge_prior_weight*NumberOfObservations- ...
|
||||
NumberOfObservedVariables*NumberOfLags ...
|
||||
+1-(1:NumberOfObservedVariables)')));
|
||||
lik = .5*NumberOfObservedVariables*log(det(dsge_prior_weight*DynareDataset.info.ntobs*GXX+mXX)) ...
|
||||
+ .5*((dsge_prior_weight+1)*DynareDataset.info.ntobs-NumberOfParameters)*log(det((dsge_prior_weight+1)*DynareDataset.info.ntobs*SIGMAu)) ...
|
||||
- .5*NumberOfObservedVariables*log(det(dsge_prior_weight*DynareDataset.info.ntobs*GXX)) ...
|
||||
- .5*(dsge_prior_weight*DynareDataset.info.ntobs-NumberOfParameters)*log(det(dsge_prior_weight*DynareDataset.info.ntobs*(GYY-GYX*inv(GXX)*GYX'))) ...
|
||||
+ .5*NumberOfObservedVariables*DynareDataset.info.ntobs*log(2*pi) ...
|
||||
- .5*log(2)*NumberOfObservedVariables*((dsge_prior_weight+1)*DynareDataset.info.ntobs-NumberOfParameters) ...
|
||||
+ .5*log(2)*NumberOfObservedVariables*(dsge_prior_weight*DynareDataset.info.ntobs-NumberOfParameters) ...
|
||||
lik = .5*NumberOfObservedVariables*log(det(dsge_prior_weight*NumberOfObservations*GXX+mXX)) ...
|
||||
+ .5*((dsge_prior_weight+1)*NumberOfObservations-NumberOfParameters)*log(det((dsge_prior_weight+1)*NumberOfObservations*SIGMAu)) ...
|
||||
- .5*NumberOfObservedVariables*log(det(dsge_prior_weight*NumberOfObservations*GXX)) ...
|
||||
- .5*(dsge_prior_weight*NumberOfObservations-NumberOfParameters)*log(det(dsge_prior_weight*NumberOfObservations*(GYY-GYX*inv(GXX)*GYX'))) ...
|
||||
+ .5*NumberOfObservedVariables*NumberOfObservations*log(2*pi) ...
|
||||
- .5*log(2)*NumberOfObservedVariables*((dsge_prior_weight+1)*NumberOfObservations-NumberOfParameters) ...
|
||||
+ .5*log(2)*NumberOfObservedVariables*(dsge_prior_weight*NumberOfObservations-NumberOfParameters) ...
|
||||
- prodlng1 + prodlng2;
|
||||
else% Evaluation of the likelihood of the dsge-var model when the dsge prior weight is infinite.
|
||||
iGXX = inv(GXX);
|
||||
SIGMAu = GYY - GYX*iGXX*transpose(GYX);
|
||||
PHI = iGXX*transpose(GYX);
|
||||
lik = DynareDataset.info.ntobs * ( log(det(SIGMAu)) + NumberOfObservedVariables*log(2*pi) + ...
|
||||
trace(inv(SIGMAu)*(mYY - transpose(mYX*PHI) - mYX*PHI + transpose(PHI)*mXX*PHI)/DynareDataset.info.ntobs));
|
||||
lik = NumberOfObservations * ( log(det(SIGMAu)) + NumberOfObservedVariables*log(2*pi) + ...
|
||||
trace(inv(SIGMAu)*(mYY - transpose(mYX*PHI) - mYX*PHI + transpose(PHI)*mXX*PHI)/NumberOfObservations));
|
||||
lik = .5*lik;% Minus likelihood
|
||||
end
|
||||
|
||||
|
|
@ -254,7 +258,7 @@ if (nargout==9)
|
|||
iGXX = inv(GXX);
|
||||
prior.SIGMAstar = GYY - GYX*iGXX*GYX';
|
||||
prior.PHIstar = iGXX*transpose(GYX);
|
||||
prior.ArtificialSampleSize = fix(dsge_prior_weight*DynareDataset.info.ntobs);
|
||||
prior.ArtificialSampleSize = fix(dsge_prior_weight*NumberOfObservations);
|
||||
prior.DF = prior.ArtificialSampleSize - NumberOfParameters - NumberOfObservedVariables;
|
||||
prior.iGXX = iGXX;
|
||||
end
|
||||
|
|
@ -1,4 +1,4 @@
|
|||
function bvar = dsgevar_posterior_density(deep,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
|
||||
function bvar = dsgevar_posterior_density(deep,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
|
||||
% This function characterizes the posterior distribution of a bvar with
|
||||
% a dsge prior (as in Del Negro and Schorfheide 2003) for a given value
|
||||
% of the deep parameters (structural parameters + the size of the
|
||||
|
|
@ -47,7 +47,7 @@ if ~options_.noconstant
|
|||
bvar.NumberOfVariables;
|
||||
end
|
||||
|
||||
[fval,cost_flag,info,PHI,SIGMAu,iXX,prior] = dsge_var_likelihood(deep',DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
|
||||
[fval,cost_flag,info,PHI,SIGMAu,iXX,prior] = dsge_var_likelihood(deep',DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
|
||||
|
||||
% Conditionnal posterior density of the lagged matrices (given Sigma) ->
|
||||
% Matric-variate normal distribution.
|
||||
|
|
|
|||
|
|
@ -83,6 +83,15 @@ end
|
|||
|
||||
if options_.dsge_var
|
||||
check_dsge_var_model(M_, estim_params_, bayestopt_);
|
||||
if dataset_info.missing.state
|
||||
error('Estimation::DsgeVarLikelihood: I cannot estimate a DSGE-VAR model with missing observations!')
|
||||
end
|
||||
if options_.noconstant
|
||||
var_sample_moments(options_.dsge_varlag, -1, dataset_);
|
||||
else
|
||||
% The steady state is non zero ==> a constant in the VAR is needed!
|
||||
var_sample_moments(options_.dsge_varlag, 0, dataset_);
|
||||
end
|
||||
end
|
||||
|
||||
%check for calibrated covariances before updating parameters
|
||||
|
|
@ -155,25 +164,6 @@ if ~isempty(estim_params_)
|
|||
M_ = set_all_parameters(xparam1,estim_params_,M_);
|
||||
end
|
||||
|
||||
% compute sample moments if needed (bvar-dsge)
|
||||
if options_.dsge_var
|
||||
if dataset_info.missing.state
|
||||
error('I cannot estimate a DSGE-VAR model with missing observations!')
|
||||
end
|
||||
if options_.noconstant
|
||||
evalin('base',...
|
||||
['[mYY,mXY,mYX,mXX,Ydata,Xdata] = ' ...
|
||||
'var_sample_moments(options_.first_obs,' ...
|
||||
'options_.first_obs+options_.nobs-1,options_.dsge_varlag,-1,' ...
|
||||
'options_.datafile, options_.varobs,options_.xls_sheet,options_.xls_range);'])
|
||||
else% The steady state is non zero ==> a constant in the VAR is needed!
|
||||
evalin('base',['[mYY,mXY,mYX,mXX,Ydata,Xdata] = ' ...
|
||||
'var_sample_moments(options_.first_obs,' ...
|
||||
'options_.first_obs+options_.nobs-1,options_.dsge_varlag,0,' ...
|
||||
'options_.datafile, options_.varobs,options_.xls_sheet,options_.xls_range);'])
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
%% perform initial estimation checks;
|
||||
try
|
||||
|
|
|
|||
|
|
@ -1,4 +1,4 @@
|
|||
function [DynareDataset, DatasetInfo] = makedataset(DynareOptions)
|
||||
function [DynareDataset, DatasetInfo] = makedataset(DynareOptions,initialconditions)
|
||||
|
||||
% Initialize a dataset as a dseries object.
|
||||
%
|
||||
|
|
@ -23,6 +23,12 @@ function [DynareDataset, DatasetInfo] = makedataset(DynareOptions)
|
|||
%
|
||||
% See also dynare_estimation_init
|
||||
|
||||
if nargin<2
|
||||
% If a the sample is to be used for the estimation of a VAR or DSGE-VAR model
|
||||
% the second argument must be a strictly positive integer (the number of lags).
|
||||
initialconditions = 0;
|
||||
end
|
||||
|
||||
if isempty(DynareOptions.datafile) && isempty(DynareOptions.dataset.file) && isempty(DynareOptions.dataset.series)
|
||||
if gsa_flag
|
||||
DynareDataset = dseries();
|
||||
|
|
@ -174,11 +180,19 @@ else
|
|||
end
|
||||
end
|
||||
|
||||
% Add initial conditions if needed
|
||||
FIRSTOBS = firstobs-initialconditions;
|
||||
|
||||
% 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 FIRSTOBS belongs to DynareDataset.dates
|
||||
if initialconditions && FIRSTOBS<DynareDataset.init
|
||||
error(sprintf('first_obs (%s) - %i cannot be less than the first date in the dataset (%s)!\nReduce the number of lags in the VAR model or increase the value of first_obs.', char(firstobs), initialconditions, char(DynareDataset.init)));
|
||||
end
|
||||
|
||||
% Check that lastobs belongs to DynareDataset.dates...
|
||||
if newdatainterface
|
||||
if lastobs>DynareDataset.dates(end)
|
||||
|
|
@ -192,7 +206,7 @@ else
|
|||
end
|
||||
|
||||
% Select a subsample.
|
||||
DynareDataset = DynareDataset(firstobs:lastobs);
|
||||
DynareDataset = DynareDataset(FIRSTOBS:lastobs);
|
||||
|
||||
% Initialize DatasetInfo structure.
|
||||
DatasetInfo = struct('missing', struct('state', NaN, 'aindex', [], 'vindex', [], 'number_of_observations', NaN, 'no_more_missing_observations', NaN), ...
|
||||
|
|
|
|||
|
|
@ -1,29 +1,28 @@
|
|||
function [YtY,XtY,YtX,XtX,Y,X] = ...
|
||||
var_sample_moments(FirstObservation,LastObservation,qlag,var_trend_order,datafile,varobs,xls_sheet,xls_range)
|
||||
function var_sample_moments(nlag, var_trend_order, dataset_)%datafile,varobs,xls_sheet,xls_range)
|
||||
% Computes the sample moments of a VAR model.
|
||||
%
|
||||
% The VAR(p) model is defined by:
|
||||
%
|
||||
% y_t = \sum_{k=1}^p y_{t-k} A_k + z_t C + e_t for t = 1,...,T
|
||||
% y_t = \sum_{k=1}^p y_{t-k} A_k + z_t C + e_t for t = 1,...,T
|
||||
%
|
||||
% where y_t is a 1*m vector of observed endogenous variables, p is the
|
||||
% number of lags, A_k is an m*m real matrix, z_t is a 1*q vector of
|
||||
% exogenous (deterministic) variables, C is a q*m real matrix and
|
||||
% e_t is a vector of exogenous stochastic shocks. T is the number
|
||||
% of observations. The deterministic exogenous variables are assumed to
|
||||
% be a polynomial trend of order q = "var_trend_order".
|
||||
% of observations. The deterministic exogenous variables are assumed to
|
||||
% be a polynomial trend of order q = "var_trend_order".
|
||||
%
|
||||
% We define:
|
||||
% We define:
|
||||
%
|
||||
% <> Y = (y_1',y_2',...,y_T')' a T*m matrix,
|
||||
%
|
||||
% <> x_t = (y_{t-1},y_{t-2},...,y_{t-p},z_t) a 1*(mp+q) row vector,
|
||||
% <> x_t = (y_{t-1},y_{t-2},...,y_{t-p},z_t) a 1*(mp+q) row vector,
|
||||
%
|
||||
% <> X = (x_1',x_2',...,x_T')' a T*(mp+q) matrix,
|
||||
% <> X = (x_1',x_2',...,x_T')' a T*(mp+q) matrix,
|
||||
%
|
||||
% <> E = (e_1',e_2',...,e_T')' a T*m matrix and
|
||||
%
|
||||
% <> A = (A_1',A_2',...,A_p',C')' an (mp+q)*m matrix of coefficients.
|
||||
% <> A = (A_1',A_2',...,A_p',C')' an (mp+q)*m matrix of coefficients.
|
||||
%
|
||||
% So that we can equivalently write the VAR(p) model using the following
|
||||
% matrix representation:
|
||||
|
|
@ -31,18 +30,17 @@ function [YtY,XtY,YtX,XtX,Y,X] = ...
|
|||
% Y = X * A +E
|
||||
%
|
||||
%
|
||||
% INPUTS
|
||||
% o FirstObservation [integer] First observation.
|
||||
% o LastObservation [integer] Last observation.
|
||||
% o qlag [integer] Number of lags in the VAR model.
|
||||
% o var_trend_order [integer] Order of the polynomial exogenous trend:
|
||||
% INPUTS
|
||||
% o nlag [integer] Number of lags in the VAR model.
|
||||
% o var_trend_order [integer] Order of the polynomial exogenous trend:
|
||||
% = -1 no constant and no linear trend,
|
||||
% = 0 constant and no linear trend,
|
||||
% = 1 constant and linear trend.
|
||||
% o dataset_ [dseries] The sample.
|
||||
%
|
||||
% OUTPUTS
|
||||
% OUTPUTS
|
||||
% o YtY [double] Y'*Y an m*m matrix.
|
||||
% o XtY [double] X'*Y an (mp+q)*m matrix.
|
||||
% o XtY [double] X'*Y an (mp+q)*m matrix.
|
||||
% o YtX [double] Y'*X an m*(mp+q) matrix.
|
||||
% o XtX [double] X'*X an (mp+q)*(mp+q) matrix.
|
||||
% o Y [double] Y a T*m matrix.
|
||||
|
|
@ -50,8 +48,11 @@ function [YtY,XtY,YtX,XtX,Y,X] = ...
|
|||
%
|
||||
% SPECIAL REQUIREMENTS
|
||||
% None.
|
||||
%
|
||||
% REMARKS
|
||||
% Outputs are saved in the base workspace (not returned by the function).
|
||||
|
||||
% Copyright (C) 2007-2009 Dynare Team
|
||||
% Copyright (C) 2007-2014 Dynare Team
|
||||
%
|
||||
% This file is part of Dynare.
|
||||
%
|
||||
|
|
@ -68,42 +69,35 @@ function [YtY,XtY,YtX,XtX,Y,X] = ...
|
|||
% You should have received a copy of the GNU General Public License
|
||||
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
X = [];
|
||||
Y = [];
|
||||
YtY = [];
|
||||
YtX = [];
|
||||
XtY = [];
|
||||
XtX = [];
|
||||
LastObservation = dataset_.dates(end);
|
||||
FirstObservation = dataset_.dates(1)+nlag;
|
||||
|
||||
data = read_variables(datafile,char(varobs),[],xls_sheet,xls_range);
|
||||
NumberOfObservations = LastObservation-FirstObservation+1;
|
||||
NumberOfVariables = dataset_.vobs;
|
||||
|
||||
if qlag > FirstObservation
|
||||
error('VarSampleMoments :: not enough data to initialize! Try to increase FirstObservation.')
|
||||
return
|
||||
end
|
||||
|
||||
NumberOfObservations = LastObservation-FirstObservation+1;% This is T.
|
||||
NumberOfVariables = length(varobs);% This is m.
|
||||
if var_trend_order == -1% No constant no linear trend case.
|
||||
X = zeros(NumberOfObservations,NumberOfVariables*qlag);
|
||||
elseif var_trend_order == 0% Constant and no linear trend case.
|
||||
X = ones(NumberOfObservations,NumberOfVariables*qlag+1);
|
||||
indx = NumberOfVariables*qlag+1;
|
||||
elseif var_trend_order == 1;% Constant and linear trend case.
|
||||
X = ones(NumberOfObservations,NumberOfVariables*qlag+2);
|
||||
indx = NumberOfVariables*qlag+1:NumberOfVariables*qlag+2;
|
||||
if isequal(var_trend_order,-1)
|
||||
% No constant no linear trend case.
|
||||
X = zeros(NumberOfObservations,NumberOfVariables*nlag);
|
||||
elseif isequal(var_trend_order,0)
|
||||
% Constant and no linear trend case.
|
||||
X = ones(NumberOfObservations,NumberOfVariables*nlag+1);
|
||||
indx = NumberOfVariables*nlag+1;
|
||||
elseif isequal(var_trend_order,1)
|
||||
% Constant and linear trend case.
|
||||
X = ones(NumberOfObservations,NumberOfVariables*nlag+2);
|
||||
indx = NumberOfVariables*nlag+1:NumberOfVariables*nlag+2;
|
||||
else
|
||||
error('var_sample_moments :: trend must be equal to -1,0 or 1!')
|
||||
error('Estimation::var_sample_moments: trend must be equal to -1,0 or 1!')
|
||||
return
|
||||
end
|
||||
|
||||
% I build matrices Y and X
|
||||
Y = data(FirstObservation:LastObservation,:);
|
||||
% I build matrices Y and X
|
||||
Y = dataset_(FirstObservation:LastObservation).data;
|
||||
|
||||
for t=1:NumberOfObservations
|
||||
line = t + FirstObservation-1;
|
||||
for lag = 1:qlag
|
||||
X(t,(lag-1)*NumberOfVariables+1:lag*NumberOfVariables) = data(line-lag,:);
|
||||
currentdate = FirstObservation+(t-1);
|
||||
for lag = 1:nlag
|
||||
X(t,(lag-1)*NumberOfVariables+1:lag*NumberOfVariables) = dataset_(currentdate-lag).data;
|
||||
end
|
||||
end
|
||||
|
||||
|
|
@ -111,7 +105,9 @@ if (var_trend_order == 1)
|
|||
X(:,end) = transpose(1:NumberOfObservations)
|
||||
end
|
||||
|
||||
YtY = Y'*Y;
|
||||
YtX = Y'*X;
|
||||
XtY = X'*Y;
|
||||
XtX = X'*X;
|
||||
assignin('base', 'mYY', Y'*Y);
|
||||
assignin('base', 'mYX', Y'*X);
|
||||
assignin('base', 'mXY', X'*Y);
|
||||
assignin('base', 'mXX', X'*X);
|
||||
assignin('base', 'Ydata', Y);
|
||||
assignin('base', 'Xdata', X);
|
||||
Loading…
Reference in New Issue