*
* Computes and graphs error bands for impulse response functions for a VAR using
* Monte Carlo simulation. This assumes the model is a symmetric VAR estimated
* using the ESTIMATE instruction. This uses a Choleski factorization, though it
* can be modified fairly easily to do any just-identified factorization.
*
* Syntax:
*  @MonteVar( options )
*
* Parameters:
*  model=model to analyze
*  step=number of response steps[48]
*  draws=number of Monte Carlo draws[1000]
*  shocklabels=VECT[STRINGS] of labels for the shocks [dependent variables of model]
*  varlabels=VECT[STRING] of labels for the variables [dependent variables of model]
*  header=' title for graph '
*
* Revision Schedule:
*  04/07 Converted to procedure by Tom Doan, Estima
*
procedure MonteVar
*
option model   model
option symm    vmat       %sigma
option integer steps      48
option integer draws      1000
option vect[strings] shocklabels
option vect[strings] varlabels
option string  header
*
local integer i j k nvar
local integer ncoef wishdof draw
local rect sxx svt swish betaols betau betadraw
local symm sigmad
local rect ranc
local vect[rect] responses
local rect[series] impulses
local vect[strings] xlabel(nvar) ylabel(nvar)
local vect[integer] depvars
local vect[series] upper lower resp
local vect[integer] depvars
local label ll
local vector frac work
local real minupper maxupper minlower
*
if .not.%defined(model) {
   disp "Syntax: @MonteVar(model=model to analyze,other options)"
   return
}
compute nvar=%modelsize(model)
*
dim responses(draws) impulses(nvar,nvar)
dim upper(nvar) lower(nvar) resp(nvar)
*
compute sxx    =%decomp(%xx)
compute svt    =%decomp(inv(%nobs*vmat))
compute betaols=%modelgetcoeffs(model)
compute ncoef  =%rows(sxx)
compute wishdof=%nobs-ncoef
*
dim ranc(ncoef,nvar)
*
* This uses antithetic acceleration to improve convergence of the Monte Carlo
* draws. 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)).
*
* A draw for the coefficients is computed by postmultiplying the Kronecker product
* of the decompositions of the sigma draw and X’X inverse by a vector of standard
* Normal draws. For even values of "draws", we take the reflection of the
* original draw around the OLS vector.
*
dim responses(draws) impulses(nvar,nvar)
*
infobox(action=define,progress,lower=1,upper=draws) "Monte Carlo Integration"
do draw = 1,draws
   if %clock(draw,2)==1 {
      compute sigmad  =%ranwisharti(svt,wishdof)
      compute swish   =%decomp(sigmad)
      compute ranc    =%ran(1.0)
      compute betau   =sxx*ranc*tr(swish)
      compute betadraw=betaols+betau
   }
   else
      compute betadraw=betaols-betau

   compute %modelsetcoeffs(model,betadraw)
   *
   *  The decomp=swish is where the particular factorization of sigmad is input into
   *  the calculation. (While swish=%decomp(sigmad) is used in drawing the
   *  coefficients, that is because all we need at that step is any factorization of
   *  the covariance matrix, and the Choleski is the simplest one).
   *
   impulse(noprint,model=model,decomp=swish,results=impulses,steps=steps)
   *
   *  Store the impulse responses
   *
   dim responses(draw)(nvar*nvar,steps)
   ewise responses(draw)(i,j)=impulses((i-1)/nvar+1,%clock(i,nvar))(j)
   infobox(current=draw)
end do draws
infobox(action=remove)
*
* Restore the original coefficients
*
compute %modelsetcoeffs(model,betaols)
*
* Set up the labels for the matrix of graphs
*
dim xlabel(nvar) ylabel(nvar)

compute depvars=%modeldepvars(model)

if %defined(shocklabels)
   ewise xlabel(i)=shocklabels(i)
else
   ewise xlabel(i)=%l(depvars(i))

if %defined(varlabels)
   ewise ylabel(i)=varlabels(i)
else
   ewise ylabel(i)=%l(depvars(i))

grparm(bold) hlabel 18 matrixlabels 14
grparm  axislabel 24
spgraph(header=header,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.
*
dim work(draws)
do i=1,nvar
   compute minlower=maxupper=0.0
   do j=1,nvar
      clear lower(j) upper(j) resp(j)
      do k=1,steps
         ewise work(t)=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

   do j=1,nvar
      graph(ticks,min=minlower,max=maxupper,number=0) 3 j i
      # resp(j)  1 steps
      # upper(j) 1 steps 2
      # lower(j) 1 steps 2
   end do j
end do i
*
spgraph(done)
end







