BVARSELECTION - Korobilis (2011)

[jonasdovern]

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