*
*  Example MONTEVAR.PRG
*  Updated as described in the RATS Newsletter from November 2002
*****************************************************************
/*
Between here and next line of ****
a)  Set the four variables (lags, nvar, nstep, ndraws)
b)  Read the data and perform any required transformations
c)  Replace the VARIABLES list in the system definition with your variables
*/
compute lags=4			;*Number of lags
compute nvar=2			;*Number of variables
compute nstep=16		;*Number of response steps
compute ndraws=2500		;*Number of draws
*
cal 1959 1 4
allocate 1998:4
*
open data drisampl.rat
data(format=rats) / gdpq fm1
log gdpq
log fm1
/*
Set up the system
*/
system(model=varmodel)
variables gdpq fm1
lags 1 to lags
kfset xxx
det constant
end(system)
******************************************************************

estimate(outsigma=vmat)

declare rect sxx svtr swish betaols betadraw
declare symm sigmad

compute sxx    =%decomp(xxx)
compute svtr   =tr(%decomp(vmat))
compute betaols=%modelgetcoeffs(varmodel)
compute ncoef  =%rows(sxx)
compute wishdof=%nobs-ncoef

declare rect ranc(ncoef,nvar)

*
*  This code handles the antithetic draws differently than the original MONTEVAR program.
*  Instead of doubling up on ndraws and combining results of the original draw with its
*  antithete, it just does a new draw on the odd values, and the antithetic draw on even ones.
*  Each is stored separately.
*
*  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.
*
*
declare vect[rect] responses(ndraws)
declare rect[series] impulses(nvar,nvar)

infobox(action=define,progress,lower=1,upper=ndraws) 'Monte Carlo Integration'
do draws = 1,ndraws
   if %clock(draws,2)==1 {
      compute sigmad  =%mqform(%nobs*inv(%ranwishart(nvar,wishdof)),svtr)
      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(varmodel,betadraw)
   impulse(noprint,model=varmodel,decomp=swish,results=impulses) nvar 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 ndraws
   do j=1,nvar
      clear lower(j) upper(j) resp(j)
      do k=1,nstep
         set work   1 ndraws = 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)