*
* IMPORT.PRG
* Manual example 13.6
*
cal 1959 1 4
compute lags=4,nvar=6,nstep=16,ndraws=10000
compute ncoef=nvar*lags+1
allocate 2000:4
*
open data haversmp.rat
data(format=rats) / ffed dgnp gnph fnh lr fm1
set ffed = ffed*.01
set lr   = lr*.01
log gnph
log dgnp
log fnh
log fm1
*
*  In this section, we set up a system for a structural model.
*  First, estimate the VAR by OLS
*
system(model=varmodel)
variables ffed fm1 gnph dgnp lr fnh
lags 1 to lags
det constant
kfset xxx
end(system)
*
estimate(outsigma=vmat)
*
* Create the set of non-linear parameters, and the "A" formula
* for CVMODEL.
*
nonlin(parmset=simszha) a12 a21 a23 a24 a31 a36 a41 a43 a46 a51 a53 a54 a56
compute nfree=13
*
dec frml[rect] afrml
frml afrml = ||1.0   ,a12  ,0.0   ,0.0   ,0.0 ,0.0|$
               a21   ,1.0  ,a23   ,a24   ,0.0 ,0.0|$
               a31   ,0.0  ,1.0   ,0.0   ,0.0 ,a36|$
               a41   ,0.0  ,a43   ,1.0   ,0.0 ,a46|$
               a51   ,0.0  ,a53   ,a54   ,1.0 ,a56|$
               0.0   ,0.0  ,0.0   ,0.0   ,0.0 ,1.0||
*
cvmodel(parmset=simszha,pmethod=simplex,piters=10,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.
*
compute delta=3.5
cvmodel(parmset=simszha,dfc=ncoef,pdf=delta,dmatrix=marginalized,method=bfgs) vmat afrml
dec rect saxx
*
*  axbase is the maximizing vector of coefficients.
*
compute [vector] axbase=%parmspeek(simszha)
*
*  saxx is a factor of the (estimated) inverse Hessian at the final estimates. This gets scaled
*  slightly to fatten up the tails a bit.
*
compute saxx=1.2*%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=30.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)
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'
*
*  Antithetic acceleration is used for the lag coefficients only, so on odd number
*  draws, we generate a new sigma matrix, and a set of changes to the lag coefficients,
*  then on the even ones, we generate a new draw with the sign switched on the deltas
*  to the lag coefficients.
*
do draws = 1,ndraws
   if %clock(draws,2)==1 {
*
*        Do a draw for the coefficients from a multivariate t density and “poke” it back
*        into the parmset so the AFRML can get the new values. Compute (log kernels of)
*        the true marginal posterior density and the importance function.
*
      compute grandom =nu/%rangamma(nu)
      compute au      =%ran(sqrt(grandom))
      compute %parmspoke(simszha,axbase+saxx*au)
      compute a       =afrml(1)
      compute dhat    =a*vmat*tr(a)
      compute ddiag   =%xdiag(dhat)
      compute pdensity=.5*(%nobs-ncoef)*log(%det(a*tr(a)))-(.5*(%nobs-ncoef)+delta+1)*%sum(%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*(%nobs-ncoef)+delta+1)
*
*         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/sumwt

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.
*

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)
         compute fracs=%wfractiles(work,weights,||.16,.84||)
         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)