CAL 1959 1 4
COMPUTE LAGS=4,NVAR=6,NSTEP=16,NDRAWS=10000
COMPUTE NCOEF=NVAR*LAGS+1
ALLOCATE 1989:3
*
DATA(FORMAT=DRI) / FYff GD GNPQ GNP GINQ LHUR FM1
SET FYFF = FYFF*.01
SET LHUR = LHUR*.01
LOG GNPQ
LOG GD
LOG GINQ
LOG FM1
*
*  In this section, we set up a system for a structural model.
*
SYSTEM(MODEL=VARMODEL)
VARIABLES FYFF FM1 GNPQ GD LHUR GINQ
LAGS 1 TO LAGS
DET CONSTANT
KFSET XXX
END(SYSTEM)
*
estimate(outsigma=vmat)
*
*  The nonlinear parameters are put into a vector. This is not the
*  most convenient way to code this for estimation (as it's harder
*  to adjust the model), but it considerably cleans up the process
*  of making draws.
*
*  nfree is the number of free coefficients
*
compute nfree=13
dec vect ax(nfree)
nonlin ax
*
dec frml[rect] afrml
frml afrml = ||1.0   ,ax(1),0.0   ,0.0   ,0.0 ,0.0|$
               ax(2) ,1.0  ,ax(3) ,ax(4) ,0.0 ,0.0|$
               ax(5) ,0.0  ,1.0   ,0.0   ,0.0 ,ax(6)|$
               ax(7) ,0.0  ,ax(8) ,1.0   ,0.0 ,ax(9)|$
               ax(10),0.0  ,ax(11),ax(12),1.0 ,ax(13)|$
               0.0   ,0.0  ,0.0   ,0.0   ,0.0 ,1.0||
*
compute ax=%const(0.0)
cvmodel(iters=50,method=simplex) vmat afrml
cvmodel(iters=300,method=bfgs) vmat afrml

*
declare rect sxx swish betaols betadraw
declare symm sigmad
*
*  Save a factor of the inv(X'X) matrix from the VAR, and the OLS coefficients
*
compute sxx    =%decomp(xxx)
compute betaols=%modelgetcoeffs(varmodel)
compute ncoef  =%rows(sxx)
*

*
*  Compute the maximum of the log of the marginal posterior density for the A coefficients
*  with a prior of the form |D|**(-delta). Delta should be at least (nvar+1)/2.0 to ensure
*  an integrable posterior distribution.
*
*  Since we need the inverse Hessian at the maximum for our importance density, we use the
*  bfgs and stderrs options on FIND.
*
dec vect deltas(6)
compute deltas=||4.0,1.0,3.0,2.0,0.0,3.0||
ewise deltas(i)=%nobs-ncoef+(deltas(i)+3)
compute delta=3.5
dec vect ddiag axbase
dec symm dhat
dec real pdensity
dec rect saxx
*
compute ax=%const(0.0)
find(iters=200,method=bfgs,stderrs) maximum pdensity
   compute a=afrml(1)
   compute dhat=a*vmat*tr(a)
   compute ddiag=%xdiag(dhat)
   compute pdensity=.5*(%nobs-ncoef)*log(%det(dhat))-.5*%dot(deltas,%log(ddiag))
end find
*
*  axbase is the maximizing vector of coefficients. saxx is a factor of the (estimated)
*  inverse Hessian at the final estimates. In this code, this gets scaled up slightly to fatten
*  up the tails a bit.
*
compute axbase=ax
compute saxx  =1.1*%decomp(%xx)
*
*  scladjust is used to prevent overflows when computing the weight function
*
compute scladjust=%funcval
*
*  nu is the degrees of freedom for the multivariate Student used in drawing A's
*
compute nu=50.0
*
*  Because this computes fractiles, it's necessary to save the full impulse responses
*  for each draw. We also need a vector to hold the weights.
*
declare vect[rect] responses(ndraws)
declare vect weights(ndraws)
declare rect[series] impulses(nvar,nvar)
***declare vect au(%nreg) axbase
* Change to NFREE because 5.03 does not set %NREG on the FIND command:
declare vect au(nfree) axbase
declare rect saxx
declare vect dbase d(nvar)
declare rect ranc(ncoef,nvar)

*
*  sumwt and sumwt2 will hold the sum of the weights and the sum of squared weights
*
compute sumwt=0.0
compute sumwt2=0.0
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 {
*
*         Do a draw for the ax coefficients from a multivariate t density. Compute (log kernels of)
*         the true marginal posterior density and the importance function.
*
      compute grandom =nu/%rangamma(nu)
      compute au      =%ran(sqrt(grandom))
      compute ax      =axbase+saxx*au
      compute a       =afrml(1)
      compute dhat    =a*vmat*tr(a)
      compute ddiag   =%xdiag(dhat)
      compute pdensity=.5*(%nobs-ncoef)*log(%det(dhat))-.5*%dot(deltas,%log(ddiag))
      compute idensity=-((nu+nfree)/2.0)*log(nu+%dot(au,au))
*
*         Compute the weight value by exp'ing the difference between the two densities, with
*         scale adjustment terms to prevent overflow.
*
      compute weight  =exp(pdensity-scladjust-idensity-((nu+nfree)/2.0)*log(nu))
*
*         Conditioned on A, make a draw for the D matrix
*
      ewise d(i)      =(%nobs/2.0)*ddiag(i)/%rangamma(.5*deltas(i))
*
*         Combine D and A to generate the draw for a factor of sigma.
*
      compute swish   =inv(a)*%diag(%sqrt(d))
*
*         Compute the + draw for coefficients
*
      compute ranc    =%ran(1.0)
      compute betau   =sxx*ranc*tr(swish)
      compute betadraw=betaols+betau
   }
   else
*
*         Compute the - draw for the coefficients
*
      compute betadraw=betaols-betau

   compute %modelsetcoeffs(varmodel,betadraw)

   impulse(noprint,model=varmodel,decomp=swish,results=impulses) nvar nstep
*
*     Store the impulse responses and the draw weight
*
   compute sumwt   =sumwt+weight
   compute sumwt2  =sumwt2+weight**2
   dim responses(draws)(nvar*nvar,nstep)
   compute weights(draws)=weight
   ewise responses(draws)(i,j)=impulses((i-1)/nvar+1,%clock(i,nvar))(j)
   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*(1.0/sumwt)

*
*  We need to compute fractiles of functions where the draws need to be
*  weighted. This requires sorting each component separately alongside
*  its weights, then locating the desired fractile in the series of
*  cumulative weights. This is most efficiently done by a binary search,
*  which is here done using the procedure "wfractiles"

procedure wfractiles work workwt ndraws desired fractiles
*
*  wfractiles returns a collection of fractiles of a weighted set of
*  data.
*
*  Syntax:
*
*  @wfractiles work workwt ndraws desired fractiles
*
*    Parameters:
*       work      series of sample values. Should be defined from entries 1 to ndraws
*       workwt    corresponding weights for the entries of work. This is altered within
*                 the procedure, so be sure to send a copy if you need it for other purposes.
*       ndraws    number of draws, ending entry of work and workwt
*       desired   VECTOR of percentiles that you want
*       fractiles output VECTOR of computed fractiles
*
*  Revision Schedule
*    11/02 Written by Estima
*
type series work workwt
type integer ndraws
type vector desired *fractiles
*
local integer lb ub trial
local integer n i
local real fractile
*
*  Order the work, workwt pair on the values of work. Compute a normalized
*  accumulation of the weights.
*
order work 1 ndraws workwt
acc(scale) workwt 1 ndraws
compute n=%rows(desired)
dim fractiles(n)
do i=1,n
   compute fractile=desired(i)
   *
   * Bracket the desired fractile
   *
   compute lb=0,ub=ndraws+1
   compute trial=fix(ndraws*fractile)
   loop
      if workwt(trial)<fractile
         compute lb=trial,trial=trial+ndraws/10
      else
         compute ub=trial,trial=trial-ndraws/10

      if trial>=ndraws
         compute ub=ndraws
      else
      if trial<=0
         compute lb=1
      if lb>0.and.ub<=ndraws
         break
   end loop
   if workwt(lb)>=fractile {
      compute fractiles(i)=work(lb)
      next
   }
   if workwt(ub)<=fractile {
      compute fractiles(i)=work(ub)
      next
   }
   *
   *  Now pinch it by bisection, until (lb,ub) includes
   * at most two points.
   *
   loop
      compute trial=(lb+ub)/2
      if workwt(trial)<fractile
         compute lb=trial
      else
         compute ub=trial
      if lb>=ub-1
         break
   end loop
   compute fractiles(i)=work(ub)
end do i
end

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
*
*  The xlabels are the ones which label the component being shocked.
*  For most structural VAR's, these aren't related 1-1 to the dependent
*  variables, so an adjustment to the line below will reset the labels
*
compute xlabel=||'MS','MD','Y','P','U','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 montevar.
*

dec vect[series] upper(nvar) lower(nvar) resp(nvar)
smpl 1 nstep
do i=1,nvar
   compute minlower=maxupper=0.0
   do j=1,nvar
      clear lower(j) upper(j) resp(j)
      smpl 1 ndraws
      do k=1,nstep
         set work   1 ndraws = responses(t)((i-1)*nvar+j,k)
         set workwt 1 ndraws = weights(t)
         compute resp(j)(k)=%dot(work,workwt)/%sum(workwt)
         @wfractiles work workwt ndraws ||.16,.84|| fracs
         compute lower(j)(k)=fracs(1),upper(j)(k)=fracs(2)
      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)