From 44093cbbb996783fce279a81e973ef6287563806 Mon Sep 17 00:00:00 2001 From: adjemian Date: Fri, 26 Jan 2007 14:50:30 +0000 Subject: [PATCH] BVAR-DSGE with constant (not yet ready) git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1167 ac1d8469-bf42-47a9-8791-bf33cf982152 --- matlab/DsgeVarLikelihood.m | 134 ++++++++++++++++++------------------ matlab/dynare_estimation.m | 28 +++++--- matlab/var_sample_moments.m | 4 +- 3 files changed, 87 insertions(+), 79 deletions(-) diff --git a/matlab/DsgeVarLikelihood.m b/matlab/DsgeVarLikelihood.m index f47bb36a6..b10f44cae 100644 --- a/matlab/DsgeVarLikelihood.m +++ b/matlab/DsgeVarLikelihood.m @@ -1,7 +1,6 @@ function [fval,cost_flag,ys,trend_coeff,info,PHI,SIGMAu,iXX] = DsgeVarLikelihood(xparam1,gend) % stephane.adjemian@ens.fr -global bayestopt_ estim_params_ M_ options_ xparam1_test -global xparam_ +global bayestopt_ estim_params_ M_ options_ nvx = estim_params_.nvx; nvn = estim_params_.nvn; @@ -29,71 +28,75 @@ cost_flag = 1; nobs = size(options_.varobs,1); if options_.mode_compute ~= 1 & any(xparam1 < bayestopt_.lb) - k = find(xparam1 < bayestopt_.lb); - fval = bayestopt_.penalty*min(1e3,exp(sum(bayestopt_.lb(k)-xparam1(k)))); - cost_flag = 0; - return; + k = find(xparam1 < bayestopt_.lb); + fval = bayestopt_.penalty*min(1e3,exp(sum(bayestopt_.lb(k)-xparam1(k)))); + cost_flag = 0; + return; end if options_.mode_compute ~= 1 & any(xparam1 > bayestopt_.ub) - k = find(xparam1 > bayestopt_.ub); - fval = bayestopt_.penalty*min(1e3,exp(sum(xparam1(k)-bayestopt_.ub(k)))); - cost_flag = 0; - return; + k = find(xparam1 > bayestopt_.ub); + fval = bayestopt_.penalty*min(1e3,exp(sum(xparam1(k)-bayestopt_.ub(k)))); + cost_flag = 0; + return; end Q = M_.Sigma_e; for i=1:estim_params_.nvx - k = estim_params_.var_exo(i,1); - Q(k,k) = xparam1(i)*xparam1(i); + k = estim_params_.var_exo(i,1); + Q(k,k) = xparam1(i)*xparam1(i); end offset = estim_params_.nvx; if estim_params_.nvn - H = zeros(nobs,nobs); - for i=1:estim_params_.nvn - k = estim_params_.var_endo(i,1); - H(k,k) = xparam1(i+offset)*xparam1(i+offset); - end - offset = offset+estim_params_.nvn; + disp('DsgeVarLikelihood :: Measurement errors are not implemented!') + return +% $$$ H = zeros(nobs,nobs); +% $$$ for i=1:estim_params_.nvn +% $$$ k = estim_params_.var_endo(i,1); +% $$$ H(k,k) = xparam1(i+offset)*xparam1(i+offset); +% $$$ end +% $$$ offset = offset+estim_params_.nvn; end if estim_params_.ncx - for i=1:estim_params_.ncx - k1 =estim_params_.corrx(i,1); - k2 =estim_params_.corrx(i,2); - Q(k1,k2) = xparam1(i+offset)*sqrt(Q(k1,k1)*Q(k2,k2)); - Q(k2,k1) = Q(k1,k2); - end - [CholQ,testQ] = chol(Q); - if testQ%% The variance-covariance matrix of the structural innovations is not definite positive. - %% We have to compute the eigenvalues of this matrix in order to build the penalty. - a = eig(Q); - k = a<0; - if k > 0 - fval = bayestopt_.penalty*min(1e3,exp(sum(-a(k)))); - cost_flag = 0; - return + for i=1:estim_params_.ncx + k1 =estim_params_.corrx(i,1); + k2 =estim_params_.corrx(i,2); + Q(k1,k2) = xparam1(i+offset)*sqrt(Q(k1,k1)*Q(k2,k2)); + Q(k2,k1) = Q(k1,k2); end - end - offset = offset+estim_params_.ncx; + [CholQ,testQ] = chol(Q); + if testQ%% The variance-covariance matrix of the structural innovations is not definite positive. + %% We have to compute the eigenvalues of this matrix in order to build the penalty. + a = eig(Q); + k = a<0; + if k > 0 + fval = bayestopt_.penalty*min(1e3,exp(sum(-a(k)))); + cost_flag = 0; + return + end + end + offset = offset+estim_params_.ncx; end if estim_params_.nvn & estim_params_.ncn - for i=1:estim_params_.ncn - k1 = options_.lgyidx2varobs(estim_params_.corrn(i,1)); - k2 = options_.lgyidx2varobs(estim_params_.corrn(i,2)); - H(k1,k2) = xparam1(i+offset)*sqrt(H(k1,k1)*H(k2,k2)); - H(k2,k1) = H(k1,k2); - end - [CholH,testH] = chol(H); - if testH - a = eig(H); - k = a<0; - if k > 0 - fval = bayestopt_.penalty*min(1e3,exp(sum(-a(k)))); - cost_flag = 0; - return - end - end - offset = offset+estim_params_.ncn; + disp('DsgeVarLikelihood :: Measurement errors are not implemented!') + return +% $$$ for i=1:estim_params_.ncn +% $$$ k1 = options_.lgyidx2varobs(estim_params_.corrn(i,1)); +% $$$ k2 = options_.lgyidx2varobs(estim_params_.corrn(i,2)); +% $$$ H(k1,k2) = xparam1(i+offset)*sqrt(H(k1,k1)*H(k2,k2)); +% $$$ H(k2,k1) = H(k1,k2); +% $$$ end +% $$$ [CholH,testH] = chol(H); +% $$$ if testH +% $$$ a = eig(H); +% $$$ k = a<0; +% $$$ if k > 0 +% $$$ fval = bayestopt_.penalty*min(1e3,exp(sum(-a(k)))); +% $$$ cost_flag = 0; +% $$$ return +% $$$ end +% $$$ end +% $$$ offset = offset+estim_params_.ncn; end M_.params(estim_params_.param_vals(:,1)) = xparam1(offset+1:end); M_.Sigma_e = Q; @@ -107,27 +110,24 @@ dsge_prior_weight = M_.params(strmatch('dsge_prior_weight',M_.param_names)); bayestopt_.restrict_columns,... bayestopt_.restrict_aux); if info(1) == 1 | info(1) == 2 | info(1) == 5 - fval = bayestopt_.penalty; - cost_flag = 0; - return + fval = bayestopt_.penalty; + cost_flag = 0; + return elseif info(1) == 3 | info(1) == 4 | info(1) == 20 - fval = bayestopt_.penalty*min(1e3,exp(info(2))); - cost_flag = 0; - return + fval = bayestopt_.penalty*min(1e3,exp(info(2))); + cost_flag = 0; + return end if options_.loglinear == 1 - constant = log(SteadyState(bayestopt_.mfys)); + constant = log(SteadyState(bayestopt_.mfys)); else - constant = SteadyState(bayestopt_.mfys); + constant = SteadyState(bayestopt_.mfys); end if bayestopt_.with_trend == 1 - trend_coeff = zeros(nobs,1); - for i=1:nobs - trend_coeff(i) = evalin('base',bayestopt_.trend_coeff{i}); - end - trend = constant*ones(1,gend)+trend_coeff*(1:gend); -else - trend = constant*ones(1,gend); + disp('DsgeVarLikelihood :: Linear trend is not yet implemented!') + return +% $$$ else +% $$$ trend = constant*ones(1,gend); end %------------------------------------------------------------------------------ diff --git a/matlab/dynare_estimation.m b/matlab/dynare_estimation.m index 7d7725c96..0ee505449 100644 --- a/matlab/dynare_estimation.m +++ b/matlab/dynare_estimation.m @@ -140,7 +140,6 @@ if ~isempty(options_.unit_root_vars) end i_ur(i) = i1; end - [junk,bayestopt_.var_list_stationary] = ... setdiff((1:M_.endo_nbr)',i_ur); [junk,bayestopt_.restrict_var_list_stationary] = ... @@ -194,22 +193,31 @@ if length(options_.mode_file) > 0 & options_.posterior_mode_estimation eval(['load ' options_.mode_file ';']'); end -%% compute sample moments if needed (bvar-dsge) -if ~isempty(strmatch('dsge_prior_weight',M_.param_names)) - evalin('base',['[mYY,mXY,mYX,mXX,Ydata,Xdata] = ' ... - 'var_sample_moments(options_.first_obs,options_.first_obs+options_.nobs-1,options_.varlag,-1);']) + +%% Compute the steadyn state: +if options_.steadystate_flag% if the _steadystate.m file is provided. + [oo_.steady_state,tchek] = feval([M_.fname '_steadystate'],[],[]); +else% if the steady state file is not provided. + [dd,info] = resol(oo_.steady_state,0) + oo_.steady_state = dd.ys; clear('dd'); end -%% Compute the steadyn state if the _steadystate.m file is provided -if options_.steadystate_flag - [oo_.steady_state,tchek] = feval([M_.fname '_steadystate'],[],[]); -end initial_estimation_checks(xparam1,gend,data); if options_.mode_compute == 0 & length(options_.mode_file) == 0 return; end +%% compute sample moments if needed (bvar-dsge) +if ~isempty(strmatch('dsge_prior_weight',M_.param_names)) + if all(abs(oo_.steady_state)<10e-9) + evalin('base',['[mYY,mXY,mYX,mXX,Ydata,Xdata] = ' ... + 'var_sample_moments(options_.first_obs,options_.first_obs+options_.nobs-1,options_.varlag,-1);']) + 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_.varlag,0);']) + end +end %% Estimation of the posterior mode or likelihood mode if options_.mode_compute > 0 & options_.posterior_mode_estimation @@ -827,7 +835,7 @@ end if ~((any(bayestopt_.pshape > 0) & options_.mh_replic) | (any(bayestopt_.pshape ... > 0) & options_.load_mh_file)) | ~options_.smoother - %% ML estimation, or posterior mode without metropolis-hastings or metropolis without bayesian smooth variables + %% ML estimation, or posterior mode without metropolis-hastings or metropolis without bayesian smooth variable options_.lik_algo = 2; [atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff] = DsgeSmoother(xparam1,gend,data); for i=1:M_.endo_nbr diff --git a/matlab/var_sample_moments.m b/matlab/var_sample_moments.m index 62e0e8512..e8256f011 100644 --- a/matlab/var_sample_moments.m +++ b/matlab/var_sample_moments.m @@ -29,7 +29,7 @@ if trend == -1% No constant X = zeros(NumberOfObservations,NumberOfVariables*qlag); elseif trend == 0% Constant X = zeros(NumberOfObservations,NumberOfVariables*qlag+1); - indx = NumberOfVariables*qlag+1; + indx = NumberOfVariables*qlag+1:NumberOfVariables*qlag+NumberOfVariables; elseif trend == 1;% Constant + Trend X = zeros(NumberOfObservations,NumberOfVariables*qlag+2); indx = NumberOfVariables*qlag+1:NumberOfVariables*qlag+2; @@ -45,7 +45,7 @@ for t=1:NumberOfObservations X(t,(lag-1)*NumberOfVariables+1:lag*NumberOfVariables) = data(line-lag,:); end if trend == 0 - X(t,indx) = 1; + X(t,indx) = ones(1,NumberOfVariables); elseif trend == 1 X(t,indx) = [ 1 , t ]; end