BVAR Minnesota Prior IRF

[Jules89]

Dear Tom,

I tried to adjust the gibbsvar.rpf, with an independent Minnesota-Wishart prior, to do impulse response analysis. I deleted everything, which had to do with forecasting since I am not interested in that. To draw and store the impulse responses within the Gibbs loop I took parts of the codes in Bayesian Econometrics workbook example 6.2. Is this in general the right way to do this task? After the Gibbs sampling I tried to graph the IRFs with confidence bands by using @mcgraphirf. It graphs the IRFs, but no error bands. Why is this the case?

The resulting code is given by:

Code: Select all

*** Data

open data haversample.rat
cal(q) 1959
data(format=rats) 1959:1 2006:4 ftb3 gdph ih cbhm

set loggdp = log(gdph)
set loginv = log(ih)
set logc   = log(cbhm)


*** VAR Setup and Storage Matrices

compute lags=4			;*Number of lags
compute nstep=40		;*Number of response steps
compute nburn=500   	;*Number of burn-in draws
compute ndraws=2500	    ;*Number of keeper draws



compute tight=.1
compute other=.5

system(model=varmodel)
   variables loggdp loginv logc ftb3
   lags 1 to lags
   det constant
end(system)

estimate(ols,sigma)  ;*for initial draw of VAR coefficients
impulse(noprint,model=varmodel,steps=nstep,results=impulses)  ;*initial set of of impulse responses

compute nvar=%nvar   ;* Number of equations estimated
dec vect olssee(nvar)
ewise olssee(i)=sqrt(%sigma(i,i))   ;* standard deviations for rescaling

declare vect[rect] %%responses(ndraws);
declare rect[series] impulses(nvar,nvar)


*** Minnesota Prior 


* Create the prior mean and (square root of) precision. The olssee
* values are used to scale the precision. The prior mean is 1 for the
* first own lag in each equation, and 0 otherwise. The constant is given
* a mean of zero with a zero precision (that is, a flat prior).

compute ncoef=lags*nvar+1  ;* Number of coefficients per equation
dec rect minnmean(ncoef,nvar) minnprec(ncoef,nvar)

do i=1,nvar
   do j=1,nvar
      do l=1,lags
         compute minnmean((j-1)*lags+l,i)=(i==j.and.l==1)
         compute minnprec((j-1)*lags+l,i)=olssee(j)/olssee(i)*%if(i==j,1.0/tight,1.0/(other*tight))
      end do l
   end do j
   compute minnmean(ncoef,i)=0.0,minnprec(ncoef,i)=0.0
end do i

ewise minnprec(i,j)=minnprec(i,j)^2
compute hprior=%diag(%vec(minnprec))
compute bprior=minnmean

compute sigmad=%sigma
cmom(model=varmodel)



infobox(action=define,progress,lower=-nburn,upper=ndraws) "Gibbs Sampler"
do draw=-nburn,ndraws
   infobox(current=draw)
   *
   * Draw b given sigma
   *
   compute bdraw =%ranmvkroncmom(%cmom,inv(sigmad),hprior,bprior);* %cmom is the cross moment matrix and is always the same
   compute rssmat=%sigmacmom(%cmom,bdraw)
   *
   * Draw sigma given b
   *
   compute sigmad=%ranwisharti(%decomp(inv(rssmat)),%nobs)
   if draw<=0
      next

   compute %modelsetcoeffs(varmodel,bdraw)

   * Store the impulse responses
   *
   dim %%responses(draw)(nvar*nvar,nstep)
   local vector ix
   ewise %%responses(draw)(i,j)=ix=%vec(%xt(impulses,j)),ix(i)
 
end do draw
infobox(action=remove)

** Impulse response

@mcgraphirf(model=varmodel,center=median,percentiles=||.10,.90||)

Thank you in advance

Jules

[TomDoan]

You don't have an IMPULSE instruction inside the loop, so it's just repeating the base IRF's each time. Add

impulse(noprint,model=varmodel,steps=nstep,flatten=%%responses(draw)) ;*initial set of of impulse responses

The following can be commented out (it's for an older way of saving the responses):

*
* * Store the impulse responses
* *
* dim %%responses(draw)(nvar*nvar,nstep)
* local vector ix
* ewise %%responses(draw)(i,j)=ix=%vec(%xt(impulses,j)),ix(i)

[Jules89]

Dear Tom,

thank you for your quick reply. The new code is below for anybody who is interested:

Code: Select all

*** Data

open data haversample.rat
cal(q) 1959
data(format=rats) 1959:1 2006:4 ftb3 gdph ih cbhm

set loggdp = log(gdph)
set loginv = log(ih)
set logc   = log(cbhm)


*** VAR Setup and Storage Matrices

compute lags=4			;*Number of lags
compute nstep=40		;*Number of response steps
compute nburn=500   	;*Number of burn-in draws
compute ndraws=2500	    ;*Number of keeper draws



compute tight=.1
compute other=.5

system(model=varmodel)
   variables loggdp loginv logc ftb3
   lags 1 to lags
   det constant
end(system)

estimate(ols,sigma)  ;*for initial draw of VAR coefficients
impulse(noprint,model=varmodel,steps=nstep,results=impulses)  ;*initial set of of impulse responses


compute nvar=%nvar   ;* Number of equations estimated
dec vect olssee(nvar)   ;* vector of variances 
ewise olssee(i)=sqrt(%sigma(i,i))   ;* standard deviations for rescaling

declare vect[rect] %%responses(ndraws)
declare rect[series] impulses(nvar,nvar)



*** Minnesota Prior


* Create the prior mean and (square root of) precision. The olssee
* values are used to scale the precision. The prior mean is 1 for the
* first own lag in each equation, and 0 otherwise. The constant is given
* a mean of zero with a zero precision (that is, a flat prior).

compute ncoef=lags*nvar+1  ;* Number of coefficients per equation
dec rect minnmean(ncoef,nvar) minnprec(ncoef,nvar)

do i=1,nvar
   do j=1,nvar
      do l=1,lags
         compute minnmean((j-1)*lags+l,i)=(i==j.and.l==1)
         compute minnprec((j-1)*lags+l,i)=olssee(j)/olssee(i)*%if(i==j,1.0/tight,1.0/(other*tight))
      end do l
   end do j
   compute minnmean(ncoef,i)=0.0,minnprec(ncoef,i)=0.0
end do i

ewise minnprec(i,j)=minnprec(i,j)^2
compute hprior=%diag(%vec(minnprec))
compute bprior=minnmean

compute sigmad=%sigma
cmom(model=varmodel)



infobox(action=define,progress,lower=-nburn,upper=ndraws) "Gibbs Sampler"
do draw=-nburn,ndraws
   infobox(current=draw)
   *
   * Draw b given sigma
   *
   compute bdraw =%ranmvkroncmom(%cmom,inv(sigmad),hprior,bprior)
   compute rssmat=%sigmacmom(%cmom,bdraw)
   *
   * Draw sigma given b
   *
   compute sigmad=%ranwisharti(%decomp(inv(rssmat)),%nobs)

   if draw<=0
      next

   compute %modelsetcoeffs(varmodel,bdraw)
   impulse(noprint,model=varmodel,results=impulses,steps=nstep,flatten=%%responses(draw),cv=sigmad)

end do draw
infobox(action=remove)

** Impulse response

@mcgraphirf(model=varmodel,center=median,percentiles=||.10,.90||)

Would you in general say that this is the right way to do estimates an independent Minnesota Wishart Prior VAR with Cholesky shocks? Or would you recommend some other way?
I have two further questions.
1) In case that this is the right way to estimate such a recursive SVAR, how would I do a FEVD with the code above?
2) Do I still need the "declare rect[series] impulses(nvar,nvar)" or is that done by the "flatten" option within the impulse command?

Thank you in advance

Best

Jules

[TomDoan]

Would you in general say that this is the right way to do estimates an independent Minnesota Wishart Prior VAR with Cholesky shocks? Or would you recommend some other way?

Yes.

I have two further questions.
1) In case that this is the right way to estimate such a recursive SVAR, how would I do a FEVD with the code above?

You can use @MCFEVDTABLE applied to the saved responses.

2) Do I still need the "declare rect[series] impulses(nvar,nvar)" or is that done by the "flatten" option within the impulse command?

No. That's for the old method.

[Jules89]

I forgot one final detail I need to know.

As far as I understand the Minnesota prior given by

Code: Select all

compute ncoef=lags*nvar+1  ;* Number of coefficients per equation
dec rect minnmean(ncoef,nvar) minnprec(ncoef,nvar)

do i=1,nvar
   do j=1,nvar
      do l=1,lags
         compute minnmean((j-1)*lags+l,i)=(i==j.and.l==1)
         compute minnprec((j-1)*lags+l,i)=olssee(j)/olssee(i)*%if(i==j,1.0/tight,1.0/(other*tight))
      end do l
   end do j
   compute minnmean(ncoef,i)=0.0,minnprec(ncoef,i)=0.0
end do i

then this specification yields non decaying variances. When I want to have a Minnesota prior with decaying variances for increasing lags, because for example the more recent lags are assumed to be more likely to have nonzero values would the right code then be the following one?

Code: Select all

compute ncoef=lags*nvar+1  ;* Number of coefficients per equation
dec rect minnmean(ncoef,nvar) minnprec(ncoef,nvar)

do i=1,nvar
   do j=1,nvar
      do l=1,lags
         compute minnmean((j-1)*lags+l,i)=(i==j.and.l==1)
         compute minnprec((j-1)*lags+l,i)=olssee(j)/olssee(i)*%if(i==j,l^(-1)/tight,l^(-1)/(other*tight))
      end do l
   end do j
   compute minnmean(ncoef,i)=0.0,minnprec(ncoef,i)=0.0
end do i

Thank you in advance

[TomDoan]

That's correct.

[Jules89]

Dear Tom,

would it be compuatationally feasible to identify the Gibbs VAR above with sign restrictions?
In general would I have to do the sign restriction draws on each draw of the Gibbs loop?

Best Jules

[TomDoan]

Yes, that's what you would do. I believe Uhlig used a very weak prior in his first paper on sign restrictions.

[Jules89]

Dear Tom,

I have rewitten the Uhlig replication code, such that I am able to use a Gibbs sampler with an independent Normal-Wishart Prior, where the prior for the VAR slope coefficients is of the Minnesota type. The code is the following:

Code: Select all

open data uhligdata.xls
calendar(m) 1965:1
data(format=xls,org=columns) 1965:1 2003:12 $
gdpc1 gdpdef cprindex totresns bognonbr fedfunds
*
set gdpc1 = log(gdpc1)*100.0
set gdpdef = log(gdpdef)*100.0
set cprindex = log(cprindex)*100.0
set totresns = log(totresns)*100.0
set bognonbr = log(bognonbr)*100.0
*


source uhligfuncs.src



*** VAR Setup and Storage Matrices

compute lags=12                     ;*Number of lags
compute nstep=60                    ;*Number of response steps
compute n1= 10000                   ;*Number of draws for Gibbs-Loop
compute n2= 50000                   ;*Number of draws from the unit sphere for each VAR draw
compute nburn = 9000                ;*Number of burn-in-draws
compute keep = n1 - nburn           ;*Number of kept draws
compute nvari = 6                   ;*Number of variables
compute kmin = 1                    ;*Number of Steps constrained (Lower Bouund)
compute kmax = 5                    ;*Number of Steps constrained (Upper Bouund)
compute ncoef = nvari*lags          ;*Number of coefficients per equation
compute nkeep = 10000               ;*Number of excepted draws that we actually want


compute tight=.1
compute other=.5


system(model=varmodel)
variables gdpc1 gdpdef cprindex fedfunds bognonbr totresns
lags 1 to 12
end(system)

estimate(noprint,ols,sigma)                                       ;*initial set of estimations
impulse(noprint,model=varmodel,steps=nstep,results=impulses)      ;*initial set of of impulse responses

dec vect[strings] vl(6)
compute vl=||"real GDP","GDP price defl","Comm. Price Ind.",$
"Fed Funds Rate","Nonborr. Reserv.","Total Reserves"||


compute nvar=%nvar                                                ;* Number of equations estimated
dec vect olssee(nvar)                                             ;* vector of variances
ewise olssee(i)=sqrt(%sigma(i,i))                                 ;* standard deviations for rescaling

declare series[rect] bgibbs                                       ;* Storage Series of Rectangulars to save the VAR Coefficient Draws

gset bgibbs 1 keep = %zeros(ncoef,nvari)                          ;* To dimension the storage matrix inside the series



***********************************
********* Minnesota Prior *********
***********************************

source BVARBuildPriorMN.src

@BVARBuildPriorMN(model=varmodel,tight=tight,other=other) HBPRIORS HPRIORS
@BVARFinishPrior hbpriors hpriors bprior hprior

compute sigmad=%sigma
cmom(noprint,model=varmodel)



*********************************
********* Gibbs Sampler *********
*********************************

declare vect[rect] goodresp(nkeep) goodfevd(nkeep)
declare vector ik a(nvar) ones(nvar)
declare series[rect] irfsquared
compute ones=%const(1.0)                                          ;* Fills matrix with a constant value

compute accept = 0

infobox(action=define,progress,lower=1,upper=n1) "Gibbs Sampler"
do draw= 1,n1
   infobox(current=draw)
   *
   * Draw b given sigma
   *
   compute bdraw =%ranmvkroncmom(%cmom,inv(sigmad),hprior,bprior)
   compute rssmat=%sigmacmom(%cmom,bdraw)
   *
   * Draw sigma given b
   *
   compute sigmad = %ranwisharti(%decomp(inv(rssmat)),%nobs)
   compute swish = %decomp(sigmad)
   if draw<=nburn
      next

   compute %modelsetcoeffs(varmodel,bdraw)
   impulse(noprint,model=varmodel,results=impulses,steps=nstep,decomp=swish)

   gset irfsquared 1 1 = %xt(impulses,t).^2
   gset irfsquared 2 nstep = irfsquared{1}+%xt(impulses,t).^2


      do subdraws=1,n2
         compute a=%ransphere(nvar)
         *
         * Check that the responses have the correct signs for steps
         * 1 to KMAX+1 (+1 because in the paper, the steps used 0-based
         * subscripts, rather than the 1 based subscripts used by RATS).
         *
         if UhligAccept(a,KMIN,KMAX+1,||+4,-3,-2,-5||)==0
         goto reject
         *
         * This is an accepted draw. Copy the information out. If we
         * have enough good ones, drop out of the overall loop.
         *
         compute accept=accept+1

         dim goodresp(accept)(nstep,nvari)
         dim goodfevd(accept)(nstep,nvari)

         ewise goodresp(accept)(i,j)=ik=%xt(impulses,i)*a,ik(j)                              ;*%xt(impulses,k) pulls the [nvar x nvar] matrix of responses to the Cholesky shocks at step k.
         ewise goodfevd(accept)(i,j)=ik=(irfsquared(i)*(a.^2))./(irfsquared(i)*ones),ik(j)

         if accept>=nkeep
         break
         :reject
      end do subdraws

   if accept>=nkeep
   break

   * Bookkeeping of VAR coefficients:

   compute bgibbs(draw-nburn) = bdraw                                                       ;* Saves the Draws of the VAR Coefficients

end do draw
infobox(action=remove)


***********************************
********* Post Processing *********
***********************************

clear upper lower resp
*
spgraph(vfields=3,hfields=2,footer="Figure 6. Impulse Responses with Pure-Sign Approach")
   do i=1,nvar
      smpl 1 accept
      do k=1,nstep
         set work = goodresp(t)(k,i)
         compute frac=%fractiles(work,||.16,.50,.84||)
         compute lower(k)=frac(1)
         compute upper(k)=frac(3)
         compute resp(k)=frac(2)
      end do k

      smpl 1 nstep
      graph(ticks,number=0,picture="##.##",header="Impulse Responses for "+vl(i)) 3
      # resp
      # upper / 2
      # lower / 2
   end do i
*
spgraph(done)

When the Gibbs sampler iterates I often get the following error message:
"## MAT14. Non-invertible Matrix. Using Generalized Inverse for SYMMETRIC.
The Error Occurred At Location 61, Line 6 of loop/block
3660296 Position 2878"

Moreover I get very wired looking impulse response functions.

Would you kindly have a look at my code?

Thank you in advance

Jules

[TomDoan]

With that program and the original Uhlig data it seems to work fine.

[Jules89]

Dear Tom,

in the above posted sign restriction approach with Gibbs sampling. Is it possible so get the structural residuals for each draw in the subdraws, that fullfills the sign restrictions?

Best

Jules

[TomDoan]

You can, though they usually aren't very interesting:

https://estima.com/forum/viewtopic.php?p=11773#p11773

[Jules89]

Thanks for your response,

the porcedure calculates the shocks for four different structural shocks within the 10 variable VAR. Is it possible to calculate it also when you are just interested in one of the shocks?
For example in the above Gibbs sampling approach with sign restrictions, I calculate only a monetary policy shock. Is it possible to identify just one shock like in Uhlig (2005)(and like its done in the Gibbs sample example posted above) and then use a modification of the procedure you posted?


Best Jules

[TomDoan]

Sure. Why would you think you couldn't?

[Jules89]

I am just not sure how to code that up:

I use "forecast(model=varmodel,from=hstart,to=hend,errors=resids,static)" to get the structural residuals from the cholesky decomposition

How do I have to adjust the following part, when I just have one shock like in the Gibbs sample approach above:

Code: Select all

dim shocks(accept)(nshocks)
      compute vweights=v1~v2~v3r~v3e
      do j=1,nshocks
         set shocks(accept)(j) = %dot(%xcol(vweights,j),inv(sigmap)*%xt(resids,t))
      end do j

Thanks

Jules