BVARSELECTION - Korobilis (2011)
Posted: Wed Nov 09, 2011 7:04 am
I coded the approach for Bayesian variable selection in VAR models as proposed by Korobilis, D. (2011), VAR Forecasting Using Bayesian Variable Selection, Journal of Applied Econometrics, forthcoming. This is what I would call a beta version. If you have any hints on bugs or ways to improve the efficiency of the code, please feel free to publish them here. In particular, I wonder whether the gibbs sampler could be speed up a little bit as it becomes very slow for VARs of more than 4 or 5 variables. (I think the key line that slows down is the one where beta is computed during each iteration (comp beta_var = inv(D0_inv+tr(x_star)*V*x_star) which involves inverting a huge matrix.)
Code: Select all
****************************************************************************************************************************
* This procedure allows the estimation of Bayesian VAR models with variable selection as outlined in Korobilis, D. (2011),
* VAR Forecasting Using Bayesian Variable Selection, Journal of Applied Econometrics, forthcoming.
*
* @BVARSELECTION( options ) start end forecasts
* # varlist
* # detlist
*
* start end - first and last period of estimation period
* forecasts - (optional) vector of series with point forecasts basedn on the restricted B-VAR
* varlist - list of variables of the VAR
* detlist - list of exogenous variables in the VAR
*
* OPTIONS:
* lags [1] - lag order of VAR
* nburns [1000] - number of initial convergence iterations for Gibbs sampler
* ndraws [1000] - number of iterations for Gibbs sampler
* ftime - forecasting time
* fsteps [12] - number of forecast steps
* [report]/noreport - determines if results are displayed
* graph/[nograph] - determines if forecasts are plotted
* bprior [not used] - provides values for prior means for the model coefficients (ord. like with RATS-VAR)
* gprior [.5] - prior for probability of inclusion of parameters
* rlprior [9.] - prior for lambda corresponding to coefficients of VAR lags (with ridge prior)
* rldetprior [100.] - prior for lambda corresponding to coefficients of exogenous variables (with ridge prior)
* mlprior [.1] - prior for lambda (with Minnesota prior)
*
* Written by Jonas Dovern (based on Matlab codes by D. Korobilis),
* Kiel Economics Research & Forecasting GmbH & Co. KG (www.kieleconomics.de),
* v0.9 - 11/2011.
*
****************************************************************************************************************************
source(noecho) mcmcpostproc.src
procedure BVARSELECTION start end forecasts
*
type integer start end
type vec[ser] *forecasts
*
option integer lags 1
option integer nburns 1000
option integer ndraws 1000
option integer ftime
option integer fsteps 12
option switch report 1
option switch graph 0
option choice prior 1 ridge minnesota
option vec bprior
option real gprior .5
option real rlprior 9.
option real rldetprior 100.
option real mlprior .1
*
local integer nvar ndet i j l draw periods hfi vfi lftime
local real p_j_tilde c d
local index varlist detlist
local rec x y x2 y2 D0 D0_inv eta eta_inv beta_OLS2 sse sigma V x_star beta_var beta2 theta2
local rec[int] dettable
local vec b0 p_j gamma beta_OLS beta ux theta theta_star theta_star_star beta_mean beta_stderrs gamma_mean theta_mean theta_cd uni_var
local vec[ser] fc_ols simfc fc_gibbs
local ser[vec] beta_draws gamma_draws theta_draws
local ser[rec] sigma_draws
local ser[int] permutation
local model VARmodel
local vec[str] rowlabels
*** Read varialbes and deterministic components from cards
enter(varying,entries=nvar) varlist
linreg(noprint) varlist(1)
comp detlist = %reglist(), dettable = %tablefromrl(detlist), ndet=%nreg
*** Check data availability and set forecasting time
inquire(regressorlist) i j
# varlist
if i>start-lags.or.j<end {
disp "There are missing data points in the sample you specified for the estimation!"
disp "(Pay attention to lags!)"
return
}
if %defined(bprior).and.%rows(bprior)<>lags*nvar*nvar+nvar*ndet {
disp "You provided a prior for beta (bprior) which does not have the right number of elements!"
return
}
comp lftime = %if(%defined(ftime),ftime,end+1)
*** Define model and produce OLS forecasts
system(model=VARmodel)
variables varlist
lags 1 to lags
det detlist
end(system)
estimate(noprint) start end
forecast(noprint,model=VARmodel,nostatic,from=lftime,steps=fsteps,results=fc_ols)
*** Construct data matrices
dim y2(end-start+1,nvar) x2(end-start+1,nvar*lags+ndet)
ewise y2(i,j) = ([series] varlist(j))(start-1+i)
comp y = %vec(y2)
do i=start,end
do j=1,nvar
do l=1,lags
comp x2(i-start+1,(j-1)*lags+l) = ([series] varlist(j))(i-l)
end do l
end do j
do j=1,ndet
comp x2(i-start+1,nvar*lags+j) = ([series] dettable(1,j))(i-dettable(2,j))
end do j
end do i
comp x = %kroneker(%identity(nvar),x2)
*** Set up arrays to hold the simulated parameters
gset beta_draws 1 ndraws = %zeros(lags*nvar*nvar+nvar*ndet,1)
gset gamma_draws 1 ndraws = %zeros(lags*nvar*nvar+nvar*ndet,1)
gset theta_draws 1 ndraws = %zeros(lags*nvar*nvar+nvar*ndet,1)
gset sigma_draws 1 ndraws = %zeros(nvar,nvar)
*** Define priors
if %defined(bprior)
comp b0 = bprior; * prior mean of beta
else {
if prior==1.or.prior==3
comp b0 = %zeros(lags*nvar*nvar+nvar*ndet,1)
else if prior==2 {
comp b0 = %zeros(lags*nvar*nvar+nvar*ndet,1)
ewise b0(i) = %if(%clock(i,lags*nvar+ndet)==1,1.,.0)
}
}
comp D0 = %identity(lags*nvar*nvar+nvar*ndet); * prior precision of beta
if prior==1
ewise D0(i,j) = D0(i,j)*%if(%clock(i,lags*nvar+ndet)<=lags*nvar,rlprior,rldetprior)
else if prior==2 {
* Getting error variance from univariate models
dim uni_var(nvar)
do k=1,nvar
linreg(noprint) varlist(k) start end
# constant varlist(k){1 to lags}
comp uni_var(k) = %SEESQ
end do k
* Filling prior matrix
do i=1,nvar
do j=1,lags*nvar+ndet
if %clock(j,lags*nvar+ndet)>lags*nvar
comp D0((i-1)*(lags*nvar+ndet)+j,(i-1)*(lags*nvar+ndet)+j) = D0((i-1)*(lags*nvar+ndet)+j,(i-1)*(lags*nvar+ndet)+j)*100*uni_var(i)
else if j>(i-1)*lags.and.j<=i*lags
comp D0((i-1)*(lags*nvar+ndet)+j,(i-1)*(lags*nvar+ndet)+j) = D0((i-1)*(lags*nvar+ndet)+j,(i-1)*(lags*nvar+ndet)+j)*%clock(j,lags)^2
else
comp D0((i-1)*(lags*nvar+ndet)+j,(i-1)*(lags*nvar+ndet)+j) = D0((i-1)*(lags*nvar+ndet)+j,(i-1)*(lags*nvar+ndet)+j)*(mlprior/%clock(j,lags)^2)*(uni_var(i)/uni_var(fix((j-1)/lags)+1))
end do j
end do i
}
comp D0_inv = inv(D0)
comp p_j = gprior*%ones(lags*nvar*nvar+nvar*ndet,1); * prior probability for inclusion of parameter
comp eta = %identity(nvar)
comp eta_inv = inv(eta)
*** Initialize parameters
comp gamma = %ones(lags*nvar*nvar+nvar*ndet,1)
comp beta_OLS = inv(tr(x)*x)*tr(x)*y
comp beta_OLS2 = inv(tr(x2)*x2)*tr(x2)*y2
comp sse = tr(y2-x2*beta_OLS2)*(y2-x2*beta_OLS2)
comp sigma = sse/((end-start+1)-lags*nvar-1)
dim ux(lags*nvar*nvar+nvar*ndet)
dim fc_gibbs(nvar)
*** Running Gibbs sampler
infobox(action=define,progress,lower=-nburns,upper=ndraws) "Running Gibbs sampler"
do draw=-nburns,ndraws
infobox(action=modify,current=draw)
*
comp V = %kroneker(inv(sigma),%identity(end-start+1))
comp x_star = %ones((end-start+1)*nvar,1)*tr(gamma).*x
*** Update beta
comp beta_var = inv(D0_inv+tr(x_star)*V*x_star)
comp beta = beta_var*(D0*b0+tr(x_star)*V*y) + %decomp(beta_var)*(ux=%ran(1.0))
comp beta2 = %vectorect(beta,lags*nvar+ndet)
*** Update gamma
boot(noreplace) permutation 1 lags*nvar*nvar+nvar*ndet 1 lags*nvar*nvar+nvar*ndet
do i=1,lags*nvar*nvar+nvar*ndet
comp theta = beta.*gamma
comp theta_star = theta
comp theta_star_star = theta
comp theta_star(permutation(i)) = beta(permutation(i))
comp theta_star_star(permutation(i)) = .0
comp c = p_j(permutation(i))*exp(%scalar(-gprior*tr(y-x*theta_star)*V*(y-x*theta_star)))
comp d = (1-p_j(permutation(i)))*exp(%scalar(-gprior*tr(y-x*theta_star_star)*V*(y-x*theta_star_star)))
comp p_j_tilde = c/(c+d)
comp gamma(permutation(i)) = %ranbernoulli(p_j_tilde)
end do i
comp gamma2 = %vectorect(gamma,lags*nvar+ndet)
comp theta2 = beta2.*gamma2
*** Update sigma
comp sigma = %ranwisharti(%decomp(inv(eta_inv+tr(y2-x2*theta2)*(y2-x2*theta2))),end-start+1-(lags*nvar+ndet))
*** Save draws and make forecasts
if draw>0 {
*** Save draws
comp beta_draws(draw) = beta
comp gamma_draws(draw) = gamma
comp theta_draws(draw) = theta
comp sigma_draws(draw) = sigma
*** Make forecasts
comp %modelsetcoeffs(VARmodel,theta2)
simulate(model=VARmodel,cv=sigma,results=simfc,from=lftime,steps=fsteps)
do i=1,nvar
set fc_gibbs(i) lftime lftime+fsteps-1 = %if(draw==1,simfc(i)(t),fc_gibbs(i)(t)+simfc(i)(t))
end do i
}
end do draw
infobox(action=remove)
*** Post simulation processing
if report {
do i=1,nvar
set fc_gibbs(i) lftime lftime+fsteps-1 = fc_gibbs(i)/ndraws
end do i
@MCMCPostProc(ndraws=ndraws,mean=theta_mean,cd=theta_cd) theta_draws
@MCMCPostProc(ndraws=ndraws,mean=beta_mean,stderrs=beta_stderrs) beta_draws
@MCMCPostProc(ndraws=ndraws,mean=gamma_mean) gamma_draws
}
*** Report output
if report {
dim rowlabels(lags*nvar*nvar+nvar*ndet)
do k=1,nvar
do j=1,nvar
do l=1,lags
comp rowlabels((k-1)*(lags*nvar+ndet)+(j-1)*lags+l) = %if(j==1.and.l==1,%l(varlist(k))+"_"+%l(varlist(j))+"{"+%string(l)+"}"," "+%l(varlist(j))+"{"+%string(l)+"}")
end do l
end do j
do j=1,ndet
comp rowlabels((k-1)*(lags*nvar+ndet)+nvar*lags+j) = " "+%l(dettable(1,j))+"{"+%string(dettable(2,j))+"}"
end do j
end do k
disp; disp "BVAR has been estimated over the period" %datelabel(start) "-" %datelabel(end) "."
report(action=define,vlabels=rowlabels,hlabels=||"","OLS-Beta","Beta*","Stderrs","Gamma**","Theta***","Geweke CD****"||)
report(action=modify,fillby=columns,atrow=1,atcol=2) beta_OLS
report(action=modify,fillby=columns,atrow=1,atcol=3) beta_mean
report(action=modify,fillby=columns,atrow=1,atcol=4) beta_stderrs
report(action=modify,fillby=columns,atrow=1,atcol=5) gamma_mean
report(action=modify,fillby=columns,atrow=1,atcol=6) theta_mean
report(action=modify,fillby=columns,atrow=1,atcol=7) theta_cd
report(action=format,picture="*.##")
report(action=show,title="Outputs from Gibbs sampler")
disp "* Mean of posterior density of beta."
disp "** Average posterior probability of the restiction index."
disp "*** Mean of posterior density of beta*gamma."
disp "**** Geweke CD for theta."
}
*** Plot forecasts
if graph {
inquire(seasonal) periods
@sethfiandvfi(number=nvar,hfi=hfi,vfi=vfi)
spgraph(hfi=hfi,vfi=vfi)
do k=1,nvar
graph(header=%l(varlist(k)),key=upleft,klabels=||"Data","OLS","Restr. B-VAR"||) 3
# varlist(k) %imax(start,end-periods*3) end; # fc_ols(k) lftime lftime+fsteps-1; # fc_gibbs(k) lftime lftime+fsteps-1
end do k
spgraph(done)
}
*** Return variables if requested
if %defined(forecasts) {
dim forecasts(nvar)
do k=1,nvar
set forecasts(k) lftime lftime+fsteps-1 = fc_gibbs(k)(t)
end do k
}
end procedure BVARSELECTION