diff --git a/matlab/dsge_var_likelihood.m b/matlab/dsge_var_likelihood.m index 6c37fd1a4..f4cf0a8b2 100644 --- a/matlab/dsge_var_likelihood.m +++ b/matlab/dsge_var_likelihood.m @@ -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 \ No newline at end of file diff --git a/matlab/dsgevar_posterior_density.m b/matlab/dsgevar_posterior_density.m index 16f593f06..8e359f3ce 100644 --- a/matlab/dsgevar_posterior_density.m +++ b/matlab/dsgevar_posterior_density.m @@ -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. diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index 861af4d58..7dbca0d8d 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -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 diff --git a/matlab/utilities/dataset/makedataset.m b/matlab/utilities/dataset/makedataset.m index d6ff48d4d..6af59b065 100644 --- a/matlab/utilities/dataset/makedataset.m +++ b/matlab/utilities/dataset/makedataset.m @@ -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 firstobsDynareDataset.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), ... diff --git a/matlab/var_sample_moments.m b/matlab/var_sample_moments.m index 94cd63de6..f9d4f8b02 100644 --- a/matlab/var_sample_moments.m +++ b/matlab/var_sample_moments.m @@ -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 . -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; \ No newline at end of file +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); \ No newline at end of file