*
* MONTEVAR.PRG
* Manual Example 13.2
*
compute lags=4			;*Number of lags
compute nstep=16		;*Number of response steps
compute nkeep=2500	;*Number of keeper draws
*
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)
*
system(model=varmodel)
variables loggdp loginv logc ftb3
lags 1 to lags
det constant
end(system)
******************************************************************
estimate
compute nvar   =%nvar
compute fxx    =%decomp(%xx)
compute fwish  =%decomp(inv(%nobs*%sigma))
compute wishdof=%nobs-%nreg
compute betaols=%modelgetcoeffs(varmodel)
*
* On the odd values for draws, a draw is made from the inverse Wishart
* distribution for the covariance matrix. This assumes use of the Jeffrey's prior
* |S|**-(n+1)/2 where n is the number of equations in the VAR. The posterior for
* S with that prior is inv(S)~Wishart with T-p d.f. (p = number of explanatory
* variables per equation) and covariance matrix inv(T(S-hat)).
*
declare vect[rect] responses(nkeep)
declare rect[series] impulses(nvar,nvar)

infobox(action=define,progress,lower=1,upper=nkeep) "Monte Carlo Integration"
do draws = 1,nkeep
   if %clock(draws,2)==1 {
      compute sigmad  =%ranwisharti(fwish,wishdof)
      compute fsigma  =%decomp(sigmad)
      compute betau   =%ranmvkron(fsigma,fxx)
      compute betadraw=betaols+betau
   }
   else
      compute betadraw=betaols-betau
   compute %modelsetcoeffs(varmodel,betadraw)
   impulse(noprint,model=varmodel,factor=fsigma,results=impulses,steps=nstep)
*
*     Store the impulse responses
*
   dim responses(draws)(nvar*nvar,nstep)
   ewise responses(draws)(i,j)=impulses((i-1)/nvar+1,%clock(i,nvar))(j)
   infobox(current=draws)
end do draws

dec vect[strings] xlabel(nvar) ylabel(nvar)
dec vect[integer] depvars
compute depvars=%modeldepvars(varmodel)
do i=1,nvar
   compute ll=%l(depvars(i))
   compute xlabel(i)=ll
   compute ylabel(i)=ll
end do i
infobox(action=remove)

grparm(bold) hlabel 18 matrixlabels 14
grparm  axislabel 24
spgraph(header="Impulse responses",xpos=both,xlab=xlabel, $
        ylab=ylabel,vlab="Responses of",vfields=nvar,hfields=nvar)

*
* Because we want a common scale for all responses of a single variable, we need
* to do all the calculations for a full row of graphs first. Note, by the way,
* that this graph transposes rows and columns from the arrangement used in the
* original montevar.
*

dec vect[series] upper(nvar) lower(nvar) resp(nvar)
do i=1,nvar
   compute minlower=maxupper=0.0
   smpl 1 nkeep
   do j=1,nvar
      clear lower(j) upper(j) resp(j)
      do k=1,nstep
         set work   1 nkeep = responses(t)((i-1)*nvar+j,k)
         compute frac=%fractiles(work,||.16,.84||)
         compute lower(j)(k)=frac(1)
         compute upper(j)(k)=frac(2)
         compute resp(j)(k)=%avg(work)
      end do k
      compute maxupper=%max(maxupper,%maxvalue(upper(j)))
      compute minlower=%min(minlower,%minvalue(lower(j)))
   end do j
*
   smpl 1 nstep
   do j=1,nvar
      graph(ticks,min=minlower,max=maxupper,number=0) 3 j i
      # resp(j)
      # upper(j) / 2
      # lower(j) / 2
   end do j
end do i
*
*
spgraph(done)