*
* MONTESUR.PRG
* Example of MC integration for a near-VAR
*
* This uses importance sampling for the covariance matrix. This is likely to be
* quite efficient for a near-VAR. For a more general SUR model with explanatory
* variables which differ more substantially from equation to equation, a Gibbs
* sampling procedure might work better.
*
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
*
open data haversample.rat
calendar(q) 1959
data(format=rats) 1959:1 2006:4 fm1 gdph
*
log gdph
log fm1
*
equation gdpeq gdph
# gdph{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.
*
compute logfx  =.5*%logdetxx(%xx)-%logdetxx(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,cv=sigmad)
      *
      * Compute the true log marginal posterior for sigma, and the weight on the draw
      *
      compute logfd  = .5*%logdetxx(%xx)-.5*%nobs*%dot(%sigma,inv(sigmad))-.5*(%nobs+nvar+1)*%logdetxx(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,factor=%decomp(sigmad),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)
   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.
*

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,.50||)
         compute lower(j)(k)=frac(1)
         compute upper(j)(k)=frac(2)
         compute resp(j)(k)=frac(3)
      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)