*
*  Example MONTESUR.PRG
*  Performs Monte Carlo integration for impulse responses for a near VAR
*  As written, this requires the 5.03 upgrade.
*
******************************************************************
/*
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
d)  Create the equations needed for the near-VAR

You'll also need to change the SUR instructions farther down in the program

*/
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 full VAR system
*/

system(model=varmodel)
variables gdpq fm1
lags 1 to lags
kfset xxx
det constant
end(system)

/*
**  The actual model is a near VAR. Set up the corresponding equation
**  and create a model from them.
*/

equation gdpeq gdpq
# gdpq{1 to lags} fm1{1 to lags} constant
equation m1eq fm1
# fm1{1 to lags} constant

group surmodel gdpeq m1eq

******************************************************************
*
*  Estimate the full VAR. This is really used only for the information
*  needed to get draws for the covariance matrix.
*
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
dim betadraw(ncoef,nvar)
compute betadraw=%const(0.0)
*
*  MatrixPeek and MatrixPoke are procedures which pull information out of
*  or put information into a rectangular array at the row,column combinations
*  provided by the 2 x p integer array coords. MatrixPeek isn't used in this
*  program; MatrixPoke is used to put the coefficients from the restricted
*  regression into the coefficient matrix for a full VAR.
*
procedure matrixpoke target source coords
type rect target
type vect source
type rect[int] coords
*
local integer i
do i=1,%rows(source)
   compute target(coords(1,i),coords(2,i))=source(i)
end do i
end
*
procedure matrixpeek source target coords
type rect source
type vect *target
type rect[int] coords
*
local integer i
dim target(%cols(coords))
do i=1,%cols(coords)
   compute target(i)=source(coords(1,i),coords(2,i))
end do i
end

procedure modelpreppoke base comp coords
*
*  ModelPrepPoke fills in a "poke" coordinate array for filling
*  in a coefficient matrix for the model "base" from a stacked
*  coefficient vector for the model "comp"
*
type model base comp
type rect[int] *coords
*
local equation eqnbase eqncomp
local rect[int] tblbase tblcomp
local vect coebase coecomp
local integer sizbase sizcomp i j k n m
*
compute n=%modelsize(base)
compute m=0
do i=1,n
   compute eqncomp=%modeleqn(comp,i)
   compute m=m+%eqnsize(eqncomp)
end do i
*
dim coords(2,m)
compute m=0
*
do i=1,n
   compute eqnbase=%modeleqn(base,1)
   compute tblbase=%eqntable(eqnbase)
   compute sizbase=%eqnsize(eqnbase)
   compute eqncomp=%modeleqn(comp,i)
   compute tblcomp=%eqntable(eqncomp)
   compute sizcomp=%eqnsize(eqncomp)
   do j=1,sizcomp
      do k=1,sizbase
         if tblcomp(1,j)==tblbase(1,k).and.tblcomp(2,j)==tblbase(2,k) {
            compute m=m+1
            compute coords(1,m)=k,coords(2,m)=i
            break
         }
      end do k
   end do j
end do i
end modelpreppoke
*
*  The call to ModelPrepPoke creates a "coords" map which will put the
*  coefficients from the near-VAR into the correct slots in the full VAR
*  model. Note, however, than ModelPrepPoke won't work correctly in
*  RATS 5.02 or earlier, because of a minor bug in the %modeleqn function,
*  which reordered the equation, instead of just making a straight copy.
*  This has been corrected in 5.03.
*

@ModelPrepPoke varmodel surmodel coords


sur 2
# gdpeq
# m1eq

dec vector ransur(%nreg) surdraw(%nreg)

declare vect[rect] responses(ndraws)
declare rect[series] impulses(nvar,nvar)

infobox(action=define,progress,lower=1,upper=ndraws) 'Monte Carlo Integration'
*
*  This code handles the antithetic draws differently than the standard 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.
*
do draws = 1,ndraws
   if %clock(draws,2)==1 {
*
*        On odd draws, draw a covariance matrix from the unconditional density. Do a SUR
*        estimation with that input covariance matrix.
*
      compute sigmad  =%mqform(%nobs*inv(%ranwishart(nvar,wishdof)),svtr)
      compute swish   =%decomp(sigmad)
      sur(noprint,nosigma,isigma=sigmad) 2
      # gdpeq
      # m1eq
*
*        Make a draw from the Normal distribution for the coefficients. Poke the restricted
*        coefficients into an array set up for the full VAR.
*
      compute ransur  =%ran(1.0)
      compute betau   =%decomp(%xx)*ransur
      compute surdraw =%beta+betau
      @matrixpoke betadraw surdraw coords
      compute betagap =betadraw-betaols
   }
   else {
*
*        On even draws, just take the symmetric draw around the SUR estimator
*
      compute surdraw =%beta-betau
      @matrixpoke betadraw surdraw coords
   }

   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
infobox(action=remove)
*
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

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)