OPEN DATA "F:\japan1.xlsx"
CALENDAR(Q) 1960:1
DATA(FORMAT=XLSX,ORG=COLUMNS) 1960:01 2011:04 taxes g y
*
*
compute lags=4
compute nvar=3
compute nstep=16
compute ndraws=10000
*
set trend = T
*
* In this section, we set up a system for a structural model. First,
* estimate the VAR by OLS.
*
system(model=bp)
variables taxes g y
lags 1 to lags
end(system)
*
estimate(cvout=vmat)
compute xxx=%xx
*
* Save a factor of the inv(X'X) matrix from the VAR, and the OLS
* coefficients.
*
compute fxx    =%decomp(xxx)
compute betaols=%modelgetcoeffs(bp)
compute ncoef  =%rows(fxx)
*
* Create the set of non-linear parameters, and the "A"and "B" formula for
* CVMODEL.
* 
nonlin(parmset=simszha) a31 a32 b21
compute nfree=3
*
dec frml[rect] afrml bfrml
frml afrml = ||1.0,0.0,-1.218|0.0,1.0,0.0|a31,a32,1.0||
frml bfrml = ||1.0,0.0,0.0|b21,1.0,0.0|0.0,0.0,1.0||
*
*
declare rect fsigmad betaols betadraw
declare symm sigmad
*
* 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=2
cvmodel(parmset=simszha,dfc=ncoef,pdf=delta,dmatrix=marginalized,method=bfgs) vmat afrml bfrml
dec rect faxx
*
* axbase is the maximizing vector of coefficients.
*
compute [vector] axbase=%parmspeek(simszha)
*
* <<faxx>> is a factor of the (estimated) inverse Hessian at the final
* estimates. This gets scaled slightly to fatten up the tails a bit.
*
compute faxx=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
*
* We need to save the draws in %%RESPONSES as described in the manual
* (for use with MCGRAPHIRF). We also need a series to hold the weights.
*
set weights 1 ndraws = 0.0
declare vect au(%nreg)
declare vect dbase d(nvar)
declare rect[series] impulses(nvar,nvar)
*
* This is the standard name used by the MCGraphIRF procedure
*
declare vect[rect] %%responses(ndraws)
*
* 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 draw=1,ndraws
   if %clock(draw,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. Save the log kernel of the draw into <<idensity>>.
      *
      compute au      =%rant(nu)
      compute idensity=%ranlogkernel()
      compute %parmspoke(simszha,axbase+faxx*au)
      *
      * Evaluate the model at the draw for the coefficients. Save the
      * log likelihood.
      *
      compute a       =afrml(1)
      compute b       =bfrml(1)
      compute dhat    =a*vmat*tr(a)
      compute ddiag   =%xdiag(dhat)
      cvmodel(parmset=simszha,dfc=ncoef,pdf=delta,dmatrix=marginalized,method=evaluate) vmat afrml bfrml
      compute pdensity=%funcval
      *
      * Compute the weight value by exp'ing the difference between the
      * two densities, with a scale adjustment term to prevent overflow.
      *
      compute weight  =exp(pdensity-scladjust-idensity)
      *
      * 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 fsigmad =inv(a)*%diag(%sqrt(d))
		*
		* Compute the + draw for coefficients
		*
      compute betau   =%ranmvkron(fsigmad,fxx)
      compute betadraw=betaols+betau
   }
   else
		*
      * Compute the - draw for the coefficients
		*
      compute betadraw=betaols-betau

   compute %modelsetcoeffs(bp,betadraw)
   compute nshock=1
   impulse(noprint,model=bp,factor=fsigmad,results=impulses,steps=nstep) nvar
	*
	* Store the impulse responses and the draw weight
	*
   compute sumwt   =sumwt+weight
   compute sumwt2  =sumwt2+weight^2
   dim %%responses(draw)(nvar*3,nstep)
   compute weights(draw)=weight
   ewise %%responses(draw)(i,j)=xt=%xt(impulses,j),ix=%vec(%xcol(xt,3)~%xcol(xt,5)~%xcol(xt,6)),ix(i)
   infobox(current=draw)
end do draw
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
*
set weights 1 ndraws = weights/sumwt
*
@MCGraphIRF(model=bp,center=mean,page=all,weights=weights,$
  shocklabels=||"T","G","Y"||)