*
* MONTESUR.PRG
* Example of MC integration for a near-VAR
*
* Revision Schedule:
*   06/06 Rewritten by Tom Doan to use importance sampling on the covariance
*      matrix
*
/*
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)  Create the equations needed for the near-VAR
*/
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
*

/*
**  The actual model is a near VAR. Set up the corresponding equations
**  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 m1eq gdpeq
**********************************************
*
* Do an iterated SUR to get the maximum likelihood estimates
*
sur(iters=10,model=surmodel)
compute sigmahat=%sigma
*
* Compute the parameters for drawing from the inverse Wishart distribution
* based upon the ML estimate of sigma. The degrees of freedom correction
* is the average number of parameters per equation
*
compute svt    =%decomp(inv(%nobs*sigmahat))
compute wishdof=%nobs-float(%nreg)/nvar
*
* This is to provide a rough approximation to the log constant of
* integration in the actual marginal posterior. For a true VAR, this
* will be log |X'X|**-(nvar/2)
*
compute logfx  =.5*log(%det(%xx))-log(%det(sigmahat))*.5*%nreg/nvar
*
* Run through the model to set up the "poke" array. %modelsetcoeffs needs
* the coefficients to be in a NREG x NEQN rectangular (as for a VAR), so
* we need to be able to create a matrix of that form out of a vector of
* coefficients from the SUR model
*
compute maxsize=0
do i=1,%modelsize(surmodel)
   compute j=%eqnsize(%modeleqn(surmodel,i))
   compute maxsize=%imax(maxsize,j)
end do j
dec rect fullcoeff(maxsize,%modelsize(surmodel))
dec rect[int] slots(2,%nreg)
compute fill=1
do i=1,%modelsize(surmodel)
   do j=1,%eqnsize(%modeleqn(surmodel,i))
      compute slots(2,fill)=i,slots(1,fill)=j,fill=fill+1
   end do j
end do i
*
dec vector ransur(%nreg) surdraw(%nreg)

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

infobox(action=define,progress,lower=1,upper=ndraws) 'Monte Carlo Integration'
*
compute sumwt=0.0,sumwt2=0.0
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. Evaluate the log density of sigma
*        and the importance weight
*
      compute sigmad = %ranwisharti(svt,wishdof)
      compute logwd  = logfx-.5*%nobs*%dot(sigmahat,inv(sigmad))-.5*(wishdof+nvar+1)*log(%det(sigmad))
      sur(noprint,nosigma,model=surmodel,isigma=sigmad)
      *
      * Compute the true log marginal posterior for sigma, and the weight on the draw
      *
      compute logfd  = .5*log(%det(%xx))-.5*%nobs*%dot(%sigma,inv(sigmad))-.5*(%nobs+nvar+1)*log(%det(sigmad))
      compute wt=exp(logfd-logwd)
*
*        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
   }
   else {
*
*        On even draws, just take the symmetric draw around the SUR estimator
*
      compute surdraw =%beta-betau
   }
   compute sumwt=sumwt+wt,sumwt2=sumwt2+wt**2
   compute %matpoke(fullcoeff,slots,surdraw)

   compute %modelsetcoeffs(surmodel,fullcoeff)
   impulse(noprint,model=surmodel,decomp=%decomp(sigmad),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)
   compute weights(draws)=wt
   infobox(current=draws)
end do draws
infobox(action=remove)
*
*  The efficacy of importance sampling depends upon function being estimated,
*  but the following is a simple estimate of the number of effective draws.
*
disp 'Effective sample size' sumwt**2/sumwt2
*
*  Normalize the weights to sum to 1.0
*
compute weights=weights/sumwt
*
dec vect[strings] xlabel(nvar) ylabel(nvar)
dec vect[integer] depvars
compute depvars=%modeldepvars(surmodel)
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
      set lower(j) 1 nstep = 0.0
      set upper(j) 1 nstep = 0.0
      set resp(j)  1 nstep = 0.0
      do k=1,nstep
         set work   1 ndraws = responses(t)((i-1)*nvar+j,k)
         compute frac=%wfractiles(work,weights,||.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)