*
* Panel Data Workbook, Example 13.3
* Panel VAR with shrinkage prior
*
compute nindiv=79           ;*Number of individuals
compute nlags=1            ;*Number of lags

dec vect[strings] country(nindiv)
input country
Sri Lanka
India
Bangladesh
Nepal
Pakistan
Maldives
China
Thailand
Vietnam
Malaysia
Indonesia
Albania
Algeria
Angola
Argentina
Armania
Azerbaijan
Belize
Benin
Bolivia
Botswana
Brazil
Bulgaria
Burkina Faso
Cambodia
Cameroon
Colombia
Costa Rica
Cuba
Dominican Rep.
Ecuador
Egypt
El Salvador
Ethiopia
Fiji
Georgia
Ghana
Guatemala
Guyana
Honduras
Iran
Jameica
Jordan
Kazakhstan
Kenya
Kyrgys Republic
Lesotho
Macedonia
Madagascar
Malawi
Mali
Mauritiana
Mauritius
Moldova
Mongolia
Morocco
Mozambique
Namibia
Nicaragua
Nigeria
Panama
Paraguay
Peru
Philippines
Romania
Russian Federation
Senegal
South Africa
Suriname
SWAZILAND
Syria
Tjikistan
Tanzania
Tunisia
Uganda
Ukraine
Yeman
Zambia
Zimbabwe



*
open data abc2.xls
cal(panelobs=10) 2006
data(org=columns,format=xls) 1//2006:01 791//2015:01 gap	rgdppercapita	rgdp

dec vect[strings] longvars
compute longvars=||"Gender Gap","GDP"||
 *
set logRGDPpercapita / = log(rgdppercapita(t))
set logRGDP / = log(rgdp(t))

* Set up the model
*
system(model=varmodel)
variables gap logRGDPpercapita
lags 1 to nlags
det constant
end(system)
*
* Shrinkage based upon Minnesota prior. Same standard errors but no mean.
*
compute tight =.1
compute other =.5
compute decay =.5
*
* Compute the standard errors from OLS autoregression for use in scaling
* the prior.
*
compute dvlist=%modeldepvars(varmodel)
compute nvar  =%size(dvlist)
dec vect see(nvar)
do j=1,nvar
   linreg(noprint) dvlist(j)
   # constant dvlist(j){1 to nlags}
   compute see(j)=sqrt(%seesq)
end do j
*
* Set up the precision form for a Minnesota prior
*
compute dvtable =%tablefromrl(dvlist)
compute sampleeq=%modeleqn(varmodel,1)
compute nsize   =%eqnsize(sampleeq)
compute etable  =%eqntable(sampleeq)
dec vect hdiag(nsize*nvar)
do j=1,nvar
   compute jvar=dvlist(j)
   do i=1,nsize
      compute ivar =etable(1,i)
      compute ilag =etable(2,i)
      compute ipos =%tablefind(dvtable,ivar,0)
      compute slot =(j-1)*nsize+i
      if ipos==0
         compute hdiag(slot)=0.0
      else
         compute hdiag(slot)=ilag^(2*decay)*%if(ivar==jvar,$
            1.0/tight,see(ipos)/(tight*other*see(j)))^2
   end do i
end do j
*
compute [symm] hdelta=%diag(hdiag)
*
* Gibbs sampling for random coefficients model.
*
compute nburn =1000
compute ndraws=5000
*
* Initialize with results from pooled least squares
*
dec vect[vect] betas(nindiv)
dec vect[symm] sigmas(nindiv)
*
estimate
*
compute bdraw=%betasys
ewise sigmas(i)=%sigma
*
dec vect[symm] cmom(nindiv)
dec vect inobs(nindiv)
do i=1,nindiv
   cmom(smpl=%indiv(t)==i,model=varmodel)
   compute cmom(i)=%cmom
   compute inobs(i)=%nobs
   compute betas(i)=bdraw
end do i
*
compute varsize=%size(bdraw)
nonlin(parmset=gibbsparms) bdraw betas
dec series[vect] bgibbs
gset bgibbs 1 ndraws = %zeros(%size(%parmspeek(gibbsparms)),1)
*
infobox(action=define,lower=-nburn,upper=ndraws,progress) "Gibbs Sampling"
do draw=-nburn,ndraws
   *
   * Draw beta(i)'s given beta and sigma(i)'s given beta(i)
   *
   do i=1,nindiv
      compute betas(i)=%ranmvkroncmom(cmom(i),inv(sigmas(i)),$
         hdelta,bdraw)
      compute rssmat=%sigmacmom(%cmom,%reshape(betas(i),nsize,nvar))
      compute sigmas(i)=%ranwisharti(%decomp(inv(rssmat)),inobs(i))
   end do i
   *
   * Draw beta (called bdraw) given beta(i)'s
   *
   compute hbeta=%zeros(varsize,1)
   compute hsum =%zeros(varsize,varsize)
   do i=1,nindiv
      compute hbeta=hbeta+hdelta*betas(i),hsum=hsum+hdelta
   end do i
   compute hsum=%ginv(hsum)
   compute bdraw=hsum*hbeta+%ranmvnormal(%decomp(hsum))
   infobox(current=draw)
   if draw<=0
      next
   compute bgibbs(draw)=%parmspeek(gibbsparms)
end do draw
infobox(action=remove)
*
@mcmcpostproc(means=bmeans,cv=bcv,ndraws=ndraws) bgibbs
*
dec vect[rect[series]] irf(nindiv)
do i=1,nindiv
   compute coeffsi=%xsubvec(bmeans,i*varsize+1,(i+1)*varsize)
   compute %modelsetcoeffs(varmodel,coeffsi)
   impulse(model=varmodel,results=irf(i),$
     cv=%identity(nvar),steps=24)
end do i
list iser = 1 to nindiv
spgraph(vfields=nvar,hfields=nvar,footer="Shrinkage IRF's",$
  xlabels=longvars,ylabels=longvars)
do i=1,nvar
   do j=1,nvar
      graph(number=0,nodates) nindiv
      cards irf(iser)(i,j)
   end do j
end do i
spgraph(done)