Qual VAR - Dueker, JBES 2005

[jonasdovern]

Does anybody have code to implement the Qual VAR proposed by Michael Dueker (2005), Dynamic Forecasts of Qualitative Variables: A Qual VAR Model of U.S. Recessions, Journal of Business and Economic Statistics 23(1), 96-104?

The model augments a standard VAR with a probit-equation for the probability of a recession. Due to this non-linear equation estimating the model involves evaluating an integral nummerically, which is implemented by MCMC techniques.

Any hint - also maybe about similar models that are already coded in RATS is very much appretiated!!!!

[TomDoan]

A revised program for this is now (Sept 2011) posted on the forum at:
http://www.estima.com/forum/viewtopic.php?f=8&t=1191


That's actually fairly simple. The VAR step is a standard MC draw. The draw for sigma is also standard, except you have to rescale the last column to a unit variance when you're done. That can be done with

Code: Select all

dec vector rescale(n)
ewise rescale(i)=%if(i==n,1.0/sqrt(sigmadraw(n,n)),1.0)
compute sigmadraw=%diag(rescale)*sigmadraw*%diag(rescale)

A draw for the y* from the probit can be done with the method in koop_bp216.rpf from

http://www.estima.com/Koop's%20Bayesian%20Econometrics.shtml

However, instead of estimating the probit in isolation (as is done there), you just include the y* variable in the VAR sampler.

[jonasdovern]

Thanks, Tom. I'll put some code together and probably come back to the Forum with specific questions.

[jonasdovern]

I started coding the model and got stuck with two questions. To compute the parameters of the distribution of the latent variable I have to compute the matrices C and D given below equation (10) in the paper by Dueker.

1. Is there a function to easily extract the coefficient matrices corresponding to lags 1,...,p from %modelgetcoeffs(qualVAR)?
2. I could not yet figure out what exacly the kappas are that are used to compute D. I think it's clear to me from equations (7) and (8) how to compute kappa_t. But how are kappa_t+1, ..., kappa_t+p defined? (Where are the future values of Y gone? kappa should only contain known information...) Does anybody have a hint on how to interpret the notation of the paper?

Code: Select all

* Initial estimation of VAR by OLS
system(model=qualVAR)
   var dlgdp dlcpi spread ffr nber
   det constant
   lag   1 to 5
end(system)
estimate(model=qualVAR) 1955:3 ftime-1

comp fxx    = %decomp(%xx)
comp fwish   = %decomp(inv(%nobs*%sigma))
comp wishdof = %nobs-%nreg
comp betaols = %modelgetcoeffs(qualVAR)
comp nvar   = %rows(%sigma)

comp ystarmu    = 0.
comp ystarsigma   = 0.1


* Estimation by MCMC method
comp draws    = 1000
comp burn   = 250

decl vec rescale(nvar)
decl rec Cdraw(nvar,nvar) Ddraw(nvar,nvar)

infobox(action=define,progress,lower=1,upper=draws) "Estimation by MCMC"
do draw=1,draws
   infobox(action=modify,current=draw)
   *
   * Draw VAR coefficients
   comp betadraw    = betaols+%ranmvkron(fsigma,fxx)
   *
   * Draw and rescale VAR error variance
   comp sigmadraw   = %ranwisharti(fwish,wishdof)
   ewise rescale(i) = %if(i==nvar,1.0/sqrt(sigmadraw(nvar,nvar)),1.0)
   comp sigmadraw   = %diag(rescale)*sigmadraw*%diag(rescale)
   *
   * Draw the latent "recession variable"
   
   * Get pstarmu and pstarsigma
   comp Cdraw = inv(sigmadraw)+...


   * Draw ystar for main sample
   set ystar 1955:3+5 ftime-1 = %if(nber(t),%rantruncate(ystarmu,ystarsigma,%na,0.0-ystarmu/ystarsigma),%rantruncate(ystarmu,ystarsigma,0.0-ystarmu/ystarsigma,%na))
   * Draw ystar for first p=5 observations



   * Save single draws
   if draw>burn {



   }



end do draw
infobox(action=remove)


[TomDoan]

Sorry. It's isn't quite so simple because of the dynamics in the probit equation. With a standard probit, you can draw the y*'s independently, but here you can't since lags of y* enter into the equation for y*. The approach that he's taking is to draw each y*(t) separately, taking all other y*'s (both future and past) as given, that is, each y*(t) is a separate partition for the Gibbs sampler. That tends to produce a rather poorly behaved Gibbs sampler since, in general, you don't like to sample separately parameters that are highly correlated. However, unless there's some sneaky forward filter-backward sampler algorithm, there might not be a good alternative.

[TomDoan]

Actually, what he's doing is re-deriving the Kalman smoother. It will not be the first paper, nor the last, that didn't recognize that. (I, myself, must plead guilty to this, though the Kalman smoother wasn't as well known when I wrote my part of Doan, Litterman and Sims).

This is the little brother of the VAR, but it shows the basic idea. You set up a state-space model, and smooth the data with a measurement equation that leaves out only one data point. (You have to run a full Kalman smoother for each data point, so this takes a while). That gives you the mean and variance of ystar(time) given all the other ystar's plus the rest of the model. You then sample from the appropriate truncated Normal for ystar(time).

For the VAR, everything except the lag coefficients on the ystar end up in the "Z" component of the state space model; the state equation only needs the dynamics involving the ystars.

Code: Select all

*
* Gibbs sampling estimate of a dynamic "probit" model, where the latent
* index follows an AR(1) equation:
*
* ystar(t)=alpha+rho*ystar(t-1)+eps(t), with eps(t)~N(0,1) i.i.d.
*
cal(q) 1960
@nbercycles(up=expansion) 1960:1 2009:4
*
* Start with ystar of +/-1
*
set ystar = %if(expansion,1.0,-1.0)
*
* Start the sampler at the OLS estimates, given the starting ystars
*
linreg ystar
# constant ystar{1}
compute rho=%beta(2),alpha=%beta(1)
*
* Set up state space model for ystar
*
frml af = rho
frml zf = alpha
*
* (diffuse) prior for coefficients
*
compute bprior=%zeros(2,1)
compute hprior=%zeros(2,2)
*
compute nburn=1000
compute ndraws=1000
*
dec series[vect] bgibbs(ndraws)
gset bgibbs 1 ndraws = %zeros(2,1)
*
infobox(action=define,progress,lower=-nburn,upper=ndraws) "Gibbs Sampler: Dynamic Probit"
do draw=-nburn,ndraws
   infobox(current=draw)
   *
   * Draw ystar(time) given all other ystar's. This is done by Gibbs
   * sampling on
   *
   * ystar(t)=alpha+rho*ystar(t-1)+eps(t)  (state equation)
   * ystar(t)=ystar(t)+0 if t<>time        (measurement equation)
   *
   * Kalman smoothing gives the conditional mean and variance for
   * ystar(time).
   *
   do time=1960:1,2009:4
      dlm(y=%if(t==time,%na,ystar),a=af,z=zf,c=1.0,$
         sw=1.0,presample=ergodic,type=smooth) 1960:1 2009:4 xstates vstates
      compute cmean=%scalar(xstates(time)),cstddev=sqrt(%scalar(vstates(time)))
      *
      * Given the conditional mean and the observed value for the 0-1
      * variable, draw a truncated Normal for ystar.
      *
      compute ystar(time)=%if(expansion(time),%rantruncate(cmean,cstddev,0.0,%na),$
                                              %rantruncate(cmean,cstddev,%na,0.0))
   end do time
   *
   * Draw coefficients from the index equation given ystars
   *
   cmom
   # constant ystar{1} ystar
   compute bdraw=%ranmvpostcmom(%cmom,1.0,hprior,bprior)
   compute alpha=bdraw(1),rho=bdraw(2)
   *
   * After burn-in, save draws
   *
   if draw>0
      compute bgibbs(draw)=bdraw
end do draw
infobox(action=remove)
*
@mcmcpostproc(ndraws=ndraws,mean=bmean,stderrs=bstderrs,cd=bcd) bgibbs

[jonasdovern]

Thanks for the example!

Here is what I coded for the VAR (the resulting recession-indicator looks similar to Figure 1 of the Dueker paper, allthough a bit less smooth since I used only 100 draws):

Code: Select all

* Initial estimation of VAR by OLS
system(model=qualVAR)
   var dlgdp dlcpi spread ffr ystar
   det constant
   lag   1 to 5
end(system)
estimate(model=qualVAR) 1955:3 ftime-1

comp fxx       = %decomp(%xx)
comp fwish      = %decomp(inv(%nobs*%sigma))
comp wishdof    = %nobs-%nreg
comp betaols    = %modelgetcoeffs(qualVAR)
comp nvar      = %rows(%sigma)
comp sigmadraw   = %ranwisharti(fwish,wishdof)
decl vec rescale(nvar)
ewise rescale(i)= %if(i==nvar,1.0/sqrt(sigmadraw(nvar,nvar)),1.0)
comp sigmadraw   = %diag(rescale)*sigmadraw*%diag(rescale)
comp betadraw    = betaols+%ranmvkron(sigmadraw,fxx)

* Estimation by MCMC method
comp draws    = 100
comp burn   = 50

decl ser[vec] betagibbs
decl ser[vec] ystargibbs
gset betagibbs  1 draws = %zeros(1,%nreg*5)
gset ystargibbs 1 draws = %zeros(1,ftime-55:3+1)

comp sw = %zeros(5,5)
comp sw(1,1) = 1.0
frml zf = betadraw(26,5)+betadraw(1,5)*dlgdp{1}+betadraw(2,5)*dlgdp{2}+betadraw(3,5)*dlgdp{3}+betadraw(4,5)*dlgdp{4}+betadraw(5,5)*dlgdp{5}+   $
         betadraw(6,5)*dlcpi{1}+betadraw(7,5)*dlcpi{2}+betadraw(8,5)*dlcpi{3}+betadraw(9,5)*dlcpi{4}+betadraw(10,5)*dlcpi{5}+   $
         betadraw(11,5)*spread{1}+betadraw(12,5)*spread{2}+betadraw(13,5)*spread{3}+betadraw(14,5)*spread{4}+betadraw(15,5)*spread{5}+   $
         betadraw(16,5)*ffr{1}+betadraw(17,5)*ffr{2}+betadraw(18,5)*ffr{3}+betadraw(19,5)*ffr{4}+betadraw(20,5)*ffr{5}

infobox(action=define,progress,lower=-burn,upper=draws) "Estimation by MCMC"
do draw=-burn,draws
   infobox(action=modify,current=draw)
   *
   * Re-estimate VAR given the new estimates for ystar
   estimate(noprint,model=qualVAR) 1955:3 ftime-1
   comp fxx       = %decomp(%xx)
   comp fwish   = %decomp(inv(%nobs*%sigma))
   *
   * Get VAR coefficients
   comp betadraw    = %modelgetcoeffs(qualVAR)+%ranmvkron(sigmadraw,fxx)
   *
   * Get and rescale VAR error variance
   comp sigmadraw   = %ranwisharti(fwish,wishdof)
   ewise rescale(i) = %if(i==nvar,1.0/sqrt(sigmadraw(nvar,nvar)),1.0)
   comp sigmadraw   = %diag(rescale)*sigmadraw*%diag(rescale)
   *
   * Draw the latent "recession variable"
   do time=1955:3,ftime-1
      dlm(y=%if(t==time,%na,ystar),   $
         a=||betadraw(21,5),betadraw(22,5),betadraw(23,5),betadraw(24,5),betadraw(25,5)||~~(%identity(4)~%zeros(4,1)),   $
         z=zf~~%zeros(4,1),   $
         c=1.0~~%zeros(4,1),   $
         sw=sw,   $
         presample=ergodic,type=smooth) 1955:3 ftime-1 xstates vstates
      *         
      comp cmean         = xstates(time)(1)
      comp cstddev       = sqrt(%scalar(vstates(time)(1,1)))
      compute ystar(time) = %if(nber(time),%rantruncate(cmean,cstddev,%na,0.0),%rantruncate(cmean,cstddev,0.0,%na))
   end do time
   disp draw @+3 betadraw(1,1) sigmadraw(1,1) ystar(2008:4)

   * Save single draws
   if draw>0 {
      comp betagibbs(draw)  = %vec(betadraw)
      ewise ystargibbs(draw)(i) = ystar(55:3+i-1)
   }
end do draw
infobox(action=remove)

* Evaluate the draws from the Gibbs sampler
@mcmcpostproc(ndraws=draws,mean=ystaravg,stderrs=ystarstderrs) ystargibbs
@mcmcpostproc(ndraws=draws,mean=betaavg,stderrs=betastderrs) betagibbs

set recindic 55:3 ftime-1 = ystaravg(t-55:3+1)
graph(shad=nber) 1; # recindic


I am not sure, whether I got the updating of the VAR parameters a) completely right and b) almost efficiently programmed! E.g. is there a more efficient way than reestimating the entire VAR with the estimate-command?

[TomDoan]

The procedure for drawing from the VAR is much the same as it is in the standard VAR Monte Carlos, except that the data change from one pass to the next. Hence, the FXX, FWISH and BETAOLS statistics change as well. I've attached an alternative coding to montevar.prg which does all the calculations for drawing a VAR inside the loop. This is what you have to do when you have data generated (at least in part) by another step in a Gibbs sampler.

gibbsvaroneoff.rpf

(2.12 KiB) Downloaded 106 times

[TomDoan]

[quote="jonasdovern"]

Code: Select all

frml zf = betadraw(26,5)+betadraw(1,5)*dlgdp{1}+betadraw(2,5)*dlgdp{2}+betadraw(3,5)*dlgdp{3}+betadraw(4,5)*dlgdp{4}+betadraw(5,5)*dlgdp{5}+   $
         betadraw(6,5)*dlcpi{1}+betadraw(7,5)*dlcpi{2}+betadraw(8,5)*dlcpi{3}+betadraw(9,5)*dlcpi{4}+betadraw(10,5)*dlcpi{5}+   $
         betadraw(11,5)*spread{1}+betadraw(12,5)*spread{2}+betadraw(13,5)*spread{3}+betadraw(14,5)*spread{4}+betadraw(15,5)*spread{5}+   $
         betadraw(16,5)*ffr{1}+betadraw(17,5)*ffr{2}+betadraw(18,5)*ffr{3}+betadraw(19,5)*ffr{4}+betadraw(20,5)*ffr{5}


That's neither pretty nor flexible. A quicker and easier way to compute that is to define an "equation" which has the list of variables:

Code: Select all

equation fixedvars *
# constant dlgdp{1 to 5} dlcpi{1 to 5} ...

After you create betadraw, pull out the coefficients that you need:

Code: Select all

dec vect probzco(21)
ewise probzco(i)=%if(i==1,betadraw(26,5),betadraw(i-1,5))

Then, you can do the ZF with

Code: Select all

FRML ZF = %dot(probzco,%eqnxvector(fixedvars,t))

[jonasdovern]

Looks much simpler; thanks. But would probzco be re-evaluated during each Gibbs iteration based on the new betadraw in this case?

[TomDoan]

Looks much simpler; thanks. But would probzco be re-evaluated during each Gibbs iteration based on the new betadraw in this case?

Yes. After drawing the VAR, you need to re-set the values into probzco.

[Marcus]

Hi,
there seems to be a slight glitch in the code you posted on Jul 20.

The line:
dec series[vect] bgibbs(ndraws)
produced the following error message:
## MAT11. You Can't DIMENSION a Series

best regards

Marcus

[jonasdovern]

One - hopfully final - question relates to the transformation of the latent process (ystar) into recession probabilities.
After running the KF in each iteration I have a series ystar and a series of means for the single elements of that process:

Code: Select all

   do time=1960:1,2009:4
      dlm(y=%if(t==time,%na,ystar),a=af,z=zf,c=1.0,$
         sw=1.0,presample=ergodic,type=smooth) 1960:1 2009:4 xstates vstates
      compute cmean=%scalar(xstates(time)),cstddev=sqrt(%scalar(vstates(time)))
      *
      * Given the conditional mean and the observed value for the 0-1
      * variable, draw a truncated Normal for ystar.
      *
      compute ystar(time)=%if(expansion(time),%rantruncate(cmean,cstddev,0.0,%na),$
                                              %rantruncate(cmean,cstddev,%na,0.0))
      compute ystarmean(time)=cmean
   end do time


This series ystarmean does not necessarily have a unit variance. Is it the right approach to obtain proper recession probabilities by scaling it using its variance before using %cdf to get the probabilities?

Code: Select all

statistics(noprint) ystarmean 1960:1 2009:4
set ystarmean 1960:1 2009:4 = ystarmean(t)/sqrt(%VARIANCE)
set recprob 1960:1 2009:4 = %cdf(ystarmean(t))


[TomDoan]

No. You want to use %CDF(ystar). Note, however, that the generated ystar series in-sample is a smoothed estimate, so it will always generate a probability >.5 (and usually quite close to 1.0) for data points which are expansions and <.5 for those that are recessions.

I've attached what I believe to be a correct program. I don't know if the data are on the JBES data archive (they password protected the bloody thing!). If they are, could someone post them? The 2000 vintage data probably look quite different than the ones that I put together from current sources.

As mentioned earlier in the thread, the ugly parts of Dueker's Gibbs sampler are actually a re-derivation of the Kalman smoother. The pre-sample part that he does with M-H can also be done by extending the Kalman smoother one period before the working start of the data and doing fairly standard conditional draws from a small multivariate Normal.

A revised program for this is now (Sept 2011) posted on the forum at:
http://www.estima.com/forum/viewtopic.php?f=8&t=1191

[kampscr]

Here are the files Dueker posted on the JBES archive:

USRATESQ.TXT

(6.12 KiB) Downloaded 139 times

ROMUSM.TXT

(29.78 KiB) Downloaded 140 times