Minnesota Prior in Near VAR

[Jules89]

Dear Tom,

I would like to estimate the coefficients of a near VAR by Gibbs sampling and than do IRF analysis with Cholesky shocks. At first I thought about doing it just like in the MONTESUR.RPF. However, this uses a flat prior on the coefficients. I would like to use a Minnesota Prior for the slope coefficients. Therefore I think the Gibbs procedure for the seemingly unrelated regression is the right choice. My near-Var is contains two blocks. One is exogenous to the local economy and the other one is the block containing the local variables. An example would be a setup in which there are US variables (GDP and FFR) and German Variables (GDP, CPI, Interest Rate). The US variables shall have an impact on the German variables, but not the other way around.

Code: Select all

system(model=ExoBlock)
   variables GDPUS FFR
   lags 1 to 3
   det constant
end(system)

system(model=EndoBlock)
   variables GDPGER CPIGER INTEREST
   lags 1 to 3
   det GDPUS{1 to 3} FFR{1 to 3} constant
end(system)

compute nearvar = ExoBlock+EndoBlock.

Within the computed model "nearvar" the first two equations have 6 regressors, while the rest have 15 regressors. Therefore the standard code for the minnesota prior below does not work, because it is written for the case that the VAR equations all have the same number of regressors:

Code: Select all

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)=%if(i==j.and.l==1,1,0)                                         
         compute minnprec((j-1)*lags+l,i)=(olssee(j)/olssee(i))*%if(i==j,l/tight,l/(other*tight))           
      end do l
   end do j

   compute minnmean(ncoef,i)=0.0,minnprec(ncoef,i)=100*olssee(i)                                            

end do i

Could you give me an advice how to implement that prior into this "asymmetric" model setup? Maybe the @BVARBuildPriorMN peocedure might work?
After I have build the priors, I will simpliy use it similar to the Gibbs sampler for the seemingly unrelated regression (for example like baseball example in the Bayes Workbook)

Thank you in advance

best

Jules

[Jules89]

Dear Tom/ Community,

following my question regarding Minnesota Near VAR Gibbs, I have written some code, which I would like to share to the community. However, sofar I am not entirely sure, whether I coded up everything correctly. The data form Cushman & Zha (1997) are aused, to assure that everybody can use the code.

Code: Select all

*****************
*  Read in Data *
*****************

cal(m) 1974:1
all 1993:12

open data dc11cit.asc
data(format=free,org=columns) / exc m1boc tbill cpi ip $
   expo imphs ipus cpius ffrus wxcp

***************************************
************ General Setup ************
***************************************


compute lags=2                                 ;*Number of lags
compute nstep=50                               ;*Number of response steps
compute ndraws=5000                            ;*Number of keeper draws
compute nburn=4000                             ;*Number of burn-in draws
compute keep = ndraws - nburn                  ;*Number of kept draws
compute nvari = 11                              ;*Number of coefficients in endogenous equation
compute ncoef = nvari*lags+1                   ;*Number of coefficients per equation, +1 for intercept
compute nexo = 4                               ;*Number of variables in exogenous block
compute tight= 0.2
compute other= 0.5


*******************
*** Model Setup ***
*******************

system(model=varexo)
   variables ipus cpius ffrus wxcp
   lags 1 to lags
   det constant
end(system)

system(model=varend)
   variables exc m1boc tbill cpi ip expo imphs
   lags 1 to lags
   det ipus{1 to lags} cpius{1 to lags} ffrus{1 to lags} wxcp{1 to lags} constant
end(system)



compute fullmodel =  varexo + varend                                    ;* Creates the Near VAR, such that the exogenous Block is ordered last

@SURGibbsSetup fullmodel

sur(noprint,nosigma,model=fullmodel,cv=%identity(neqn),cmom)

compute bdraw=%vec(%modelgetcoeffs(fullmodel))                          ;*for initial draw of VAR coefficients, use vector of stacked slope coefficeints from OLS
compute nvar=%nvar

************************************************
*** Set up SDs from AR(1) Models for Scaling ***
************************************************


declare vect olsseen(nvar)                                   ;* vector of standard errors from AR(1), endogenous variables
compute i = 1
dofor x = exc m1boc tbill cpi ip expo imphs ipus cpius ffrus wxcp

   linreg(noprint) x
   # constant x{1}

   compute olsseen(i)=sqrt(%SEESQ)
   compite i=i+1

end do


declare vect olsseex(nexo)                                        ;* vector of standard errors from AR(1), exogenous variables
compute i = 1
dofor x = ipus cpius ffrus wxcp

   linreg(noprint) x
   # constant x{1}

   compute olsseex(i)=sqrt(%SEESQ)
   compite i=i+1

end dofor


*********************************
*** Set up some Storgae Space ***
*********************************


declare vect[rect] %%responses(keep)                              ;* Storage Vector of Rectanculars for Impulse Response Draws


***********************************************************************************
******** Minnesota Prior for VAR Coefficients Exogenous before Endogenous *********
***********************************************************************************

dec rect minnmean(ncoef,nvar) minnprec(ncoef,nvar)

compute minnmean = %zeros(ncoef,nvar)
compute minnprec = %zeros(ncoef,nvar)


** first do the part for the exogenous block:
do i=1,nexo
   do j=1,nexo
      do l=1,lags
         compute minnmean((j-1)*lags+l,i)=%if(i==j.and.l==1,1,0)                                                         ;* sets first on log to 1 and everything else to zero
         compute minnprec((j-1)*lags+l,i)=(olsseex(j)/olsseex(i))*%if(i==j,l/tight,l/(other*tight))                ;* computes the corresponding minnesota precisions
      end do l
   end do j

   compute minnmean(ncoef-(nvar-nexo)*lags,i)=0.0,minnprec(ncoef-(nvar-nexo)*lags,i)=100*olsseex(i)                                   ;* sets the mean and precision for the intercept

end do i


** second do the part for the endogenous block:
do i=nexo+1,nvar
   do j=1, nvar
      do l=1,lags
         compute minnmean((j-1)*lags+l,i)=%if(i-nexo==j.and.l==1,1,0)                                                         ;* sets first on log to 1 and everything else to zero
         compute minnprec((j-1)*lags+l,i)=(olsseen(j)/olsseen(i-nexo))*%if(i-nexo==j,l/tight,l/(other*tight))                ;* computes the corresponding minnesota precisions
      end do l
   end do j

   compute minnmean(ncoef,i)=0.0,minnprec(ncoef,i)=100*olsseen(i-nexo)                                   ;* sets the mean and precision for the intercept

end do i

ewise minnprec(i,j)=minnprec(i,j)^2


*** compress the precision and slope prior to get all rows out that are zero in the exogenous equations
***  to have matching dimensions in the sampler


declare vect help
compute help = %ones(%rows(%vec(minnprec)),1)

do i=1,%rows(%vec(minnprec))

   compute help(i)=%if(%vec(minnprec)(i)==0,0,1)

end do i

compute vecminnprec = %compress(%vec(minnprec),help)
compute vecminnmean = %compress(%vec(minnmean),help)


compute hprior=%diag(vecminnprec)
compute bprior=vecminnmean

********************************************
********* Flat Prior for Precision *********
********************************************

compute hxxprior=%zeros(neqn,neqn)
compute nuprior =0.0

compute wishdof=nuprior+%nobs


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



infobox(action=define,progress,lower=1,upper=ndraws) "Gibbs Sampler"
do draw= 1,ndraws
   infobox(current=draw)


   ******************************
   *** Draw Precision given b ***
   ******************************

   compute covmat=SURGibbsSigma(bdraw)+hxxprior

   compute hdraw=%ranwishartf(%decomp(inv(covmat)),wishdof)


   **************************
   *** Draw b given sigma ***
   **************************

   @SURGibbsDataInfo hdraw hdata hbdata


   compute hpost=hdata+hprior
   compute vpost=inv(hpost)
   compute bpost=vpost*(hbdata+hprior*bprior)
   compute bdraw=bpost+%ranmvnormal(%decomp(vpost))


   if draw<=nburn
      next

   *******************************************
   *** Set Modelcoefficients and draw IRFs ***
   *******************************************



   compute %modelsetcoeffs(fullmodel,bdraw)
   compute fact=%decomp(inv(hdraw))

   impulse(noprint,model=fullmodel,results=impulses,steps=nstep,flatten=%%responses(draw-nburn),factor=fact)


end do draw
infobox(action=remove)

display %modelgetcoeffs(fullmodel)


@mcgraphirf(model=fullmodel,center=median,percent=||0.16,0.84||,page=byshock)

The VAR is estimated like a SUR.
I still have some questions:

1) When I type in "display %modelgetcoeffs(fullmodel)" I get for the endogenous block variables zeros in the equations of the exogenous block variables

Code: Select all

      1.24072       0.01385       1.03698       0.07823       1.07464       0.04782      -2.33703      -0.00781       0.01577      -0.19422      -0.26280
     -0.26147   5.96474e-04      -2.41643      -0.00302      -0.06803      -0.05613       2.81767      -0.00106      -0.08700       0.14589       0.21388
      0.18931       1.09332      -9.90841      -0.18247       0.00143       0.80422       3.07745       0.01217      -0.00117      -0.06534      -0.13205
     -0.18502      -0.10546      10.01533       0.16083       0.00130       0.13264      -2.12586  -7.84323e-04       0.06404       0.10133       0.29192
 -6.96298e-05   1.49607e-04       1.15070      -0.00105   4.42312e-04      -0.00292       0.89651   1.27891e-04   4.43314e-04      -0.00162      -0.00274
 -9.91628e-04   1.69525e-04      -0.21637       0.00140   8.29788e-04       0.00133      -0.03677   9.52162e-05  -6.71114e-05  -4.81388e-04  -9.17110e-04
      0.03445       0.00966       2.89109       1.28742       0.04626       0.20593      15.32523       0.97133      -0.35942      -0.43150      -0.27949
     -0.02232      -0.00565      -1.52980      -0.32648      -0.00696      -0.23940     -10.83981      -0.01041       0.27840       0.30974       0.00366
      0.02778      -0.02847       0.01266      -0.06521       0.09301       0.03806       1.75577      -0.00147       0.78453       0.10840       0.13670
      0.00000       0.00000       0.00000       0.00000      -0.04352      -0.02862       0.21689       0.00462       0.11134      -0.15316      -0.01881
      0.00000       0.00000       0.00000       0.00000      -0.00779      -0.06572       1.35918  -6.51365e-04      -0.03249       1.00747       0.05918
      0.00000       0.00000       0.00000       0.00000      -0.00496       0.01350      -1.20777      -0.00422       0.03785      -0.03495       0.07508
      0.00000       0.00000       0.00000       0.00000       0.00192       0.01889      -0.82737      -0.00994      -0.03206      -0.08473       0.82165
      0.00000       0.00000       0.00000       0.00000       0.02726       0.03348       0.20086       0.00839       0.00531       0.03009      -0.05235
      0.00000       0.00000       0.00000       0.00000       0.07587       0.01371       2.87824       0.00978       0.04221       0.22686      -0.19240
      0.00000       0.00000       0.00000       0.00000      -0.11411      -0.01811      -6.69188      -0.00928       0.05279      -0.11289       0.10268
      0.00000       0.00000       0.00000       0.00000       0.06977       0.10793       0.13234       0.00544      -0.35410       0.03554      -0.31752
      0.00000       0.00000       0.00000       0.00000      -0.11342       0.00102      -4.99108       0.02147       0.35215       0.04720       0.44730
      0.00000       0.00000       0.00000       0.00000      -0.00141  -8.32498e-04       0.09103   3.35145e-04      -0.00174   5.13385e-05       0.00302
      0.00000       0.00000       0.00000       0.00000   6.17902e-04  -3.83159e-05       0.00884  -1.27061e-04   1.99758e-04   3.77815e-04  -3.99306e-05
      0.00000       0.00000       0.00000       0.00000      -0.00424       0.01610       2.17363      -0.00498       0.01776       0.03983   9.95922e-05
      0.00000       0.00000       0.00000       0.00000      -0.00451      -0.02452      -0.89379       0.00761       0.01028      -0.00632       0.04543
      0.00000       0.00000       0.00000       0.00000      -0.31445       0.41087      -0.00586       0.04696       0.28936       0.86634       0.35873

However I need to get rid of all those enties which are zero within the Gibbs sampler, because otherwise "hdata+hprior" have different dimensions. I did this by vectorizing the meanprecision and set up a help variable, which takes values of zero whenever the the vectorized meanprec is equal to zero. Then I used "%compress(%vec(minnprec),help)" and "%compress(%vec(minnmean),help)" to take all the zero elements out such that the dimensions for "hdata" and "hprior" match.
Is this the correct way of getting the Minnesota prior for the near VAR?

2) I checked it 100 times, but its good to have some cross-check. Are precision parameter assingned to the right slots?

Code: Select all

dec rect minnmean(ncoef,nvar) minnprec(ncoef,nvar)

compute minnmean = %zeros(ncoef,nvar)
compute minnprec = %zeros(ncoef,nvar)


** first do the part for the exogenous block:
do i=1,nexo
   do j=1,nexo
      do l=1,lags
         compute minnmean((j-1)*lags+l,i)=%if(i==j.and.l==1,1,0)                                                         ;* sets first on log to 1 and everything else to zero
         compute minnprec((j-1)*lags+l,i)=(olsseex(j)/olsseex(i))*%if(i==j,l/tight,l/(other*tight))                ;* computes the corresponding minnesota precisions
      end do l
   end do j

   compute minnmean(ncoef-(nvar-nexo)*lags,i)=0.0,minnprec(ncoef-(nvar-nexo)*lags,i)=100*olsseex(i)                                   ;* sets the mean and precision for the intercept

end do i


** second do the part for the endogenous block:
do i=nexo+1,nvar
   do j=1, nvar
      do l=1,lags
         compute minnmean((j-1)*lags+l,i)=%if(i-nexo==j.and.l==1,1,0)                                                         ;* sets first on log to 1 and everything else to zero
         compute minnprec((j-1)*lags+l,i)=(olsseen(j)/olsseen(i-nexo))*%if(i-nexo==j,l/tight,l/(other*tight))                ;* computes the corresponding minnesota precisions
      end do l
   end do j

   compute minnmean(ncoef,i)=0.0,minnprec(ncoef,i)=100*olsseen(i-nexo)                                   ;* sets the mean and precision for the intercept

end do i

ewise minnprec(i,j)=minnprec(i,j)^2

3) The code slows down extreamly when I increase the number of lags, do you see some inefficiency, which can be altered such that the code gets a bit faster?

Thank you in advance

Best


Jules

[TomDoan]

compute hpost=hdata+hprior
compute vpost=inv(hpost)
compute bpost=vpost*(hbdata+hprior*bprior)
compute bdraw=bpost+%ranmvnormal(%decomp(vpost))

can be made more efficient by using %RANMVPOSTHB, which is designed for this. It will probably cut the time by about 50%.

compute bdraw=%ranmvposthb(hdata+hprior,hbdata+hprior*bprior)

(and you can do the hprior * bprior outside the loop, though that isn't the big time-user).

Yes. It gets a lot slower as you increase the lags (or the number of variables). Inverting or factoring an N x N matrix is an O(N^3) operation. This requires dealing with a matrix the size of the full model, so if you double the size of the model (by doubling the number of lags), you multiply the time by about 8.

%MODELGETCOEFFS doesn't get used, so how it pads out the columns isn't relevant. %MODELSETCOEFFS does get used, and since it knows the structure of the target, it knows how to handle the stacked coefficient vector that the simulation process creates.

[Jules89]

Thank you very much for the answer, thats a huge time saver !

%MODELGETCOEFFS doesn't get used, so how it pads out the columns isn't relevant. %MODELSETCOEFFS does get used, and since it knows the structure of the target, it knows how to handle the stacked coefficient vector that the simulation process creates.

Sorry but I do not understand what you mean... I did not use %modelgetcoeffs in the code. I just included it in the Gibbs loop to check whether the drawn coefficients are actually zero. However, that was only to check if everything goes right.
But what I actually meant was that when I set up the minnesota prior like I do I have zero values for the precision (and mean coefficients) for the endogenous block variables in the exogenous block. Since the model, which I set up with "compute fullmodel = varexo + varend" is "asymmetric" in the sense that the the first four equations have less regressors than the rest, the vectorized minnesota mean and precision have larger dimensions than the "hdata" such that you cannot add them. Therefore I have to get rid of all those zeros, which are for the endogenous block variables in the exogenous block. Is the procedure I applied correct? I was unsure about just deleting those lines...

Again thank you very much!

Best

Jules