Page 1 of 2
BVAR Minnesota Prior IRF
Posted: Wed Mar 29, 2017 6:00 am
by 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
Re: BVAR Minnesota Prior IRF
Posted: Wed Mar 29, 2017 7:28 am
by 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)
Re: BVAR Minnesota Prior IRF
Posted: Wed Mar 29, 2017 7:57 am
by 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
Re: BVAR Minnesota Prior IRF
Posted: Wed Mar 29, 2017 9:35 am
by TomDoan
Jules89 wrote:
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.
Jules89 wrote:
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.
Jules89 wrote:
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.
Re: BVAR Minnesota Prior IRF
Posted: Wed Mar 29, 2017 10:39 am
by 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
Re: BVAR Minnesota Prior IRF
Posted: Thu Mar 30, 2017 6:23 am
by TomDoan
That's correct.
Re: BVAR Minnesota Prior IRF
Posted: Fri Jul 14, 2017 9:33 am
by 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
Re: BVAR Minnesota Prior IRF
Posted: Fri Jul 14, 2017 9:51 am
by TomDoan
Yes, that's what you would do. I believe Uhlig used a very weak prior in his first paper on sign restrictions.
Re: BVAR Minnesota Prior IRF
Posted: Tue Jul 18, 2017 3:57 am
by 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
Re: BVAR Minnesota Prior IRF
Posted: Tue Jul 18, 2017 3:34 pm
by TomDoan
With that program and the original Uhlig data it seems to work fine.
Re: BVAR Minnesota Prior IRF
Posted: Fri Jul 28, 2017 5:29 am
by 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
Re: BVAR Minnesota Prior IRF
Posted: Fri Jul 28, 2017 8:32 am
by TomDoan
You can, though they usually aren't very interesting:
https://estima.com/forum/viewtopic.php?p=11773#p11773
Re: BVAR Minnesota Prior IRF
Posted: Tue Aug 01, 2017 4:31 am
by 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
Re: BVAR Minnesota Prior IRF
Posted: Tue Aug 01, 2017 7:39 am
by TomDoan
Sure. Why would you think you couldn't?
Re: BVAR Minnesota Prior IRF
Posted: Tue Aug 01, 2017 7:57 am
by 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