open data "estima.xls"
calendar(q) 1967
compute missc=1.0e+32
DATA(FORMAT=XLS,org=columns) 1967:1 2013:4 bro fed def sp spot rgdp

set ly = log(rgdp)
set trend = t
linreg ly 1967:1 2012:4 y
# constant trend

set lconsprice = log(def)
set pi  = (lconsprice-lconsprice{1})
diff pi / dpi

set lsp = log(sp/def)
set spgr  = (lsp-lsp{1})*100

set lbro = log(bro)
set dlbro = (lbro-lbro{1})*100


set lM4 = log(spot)
set M4gr  = (lM4-lM4{1})*100

*
system(model=varmodel)
variables y dpi M4gr dlbro spgr fed
lags 1 to 5
deterministic constant
end(system)
estimate(noprint,resid=resids)

dec vect[strings] vl(6)
compute vl=||"Output","Inf","Com","Bro","SP","Fed"||

compute n1=1000
compute n2=1000
compute nkeep=2000
compute nvar=6
compute nstep=21
compute KMAX=2

* This is the standard setup for MC integration of an OLS VAR
*
dec symm s(6,6)
dec vect v1(6)      ;* For the unit vector on the 1st draw
dec vect v2(5)      ;* For the unit vector on the 2nd draw
dec vect v(6)     ;* Working impulse vector

compute sxx    =%decomp(%xx)
compute svt    =%decomp(inv(%nobs*%sigma))
compute betaols=%modelgetcoeffs(varmodel)
compute ncoef  =%rows(sxx)
compute wishdof=%nobs-ncoef
dec rect ranc(ncoef,nvar)
*
* Most draws are going to get rejected. We allow for up to 1000
* good ones. The variable accept will count the number of accepted
* draws. GOODRESP will be a RECT(nsteps,nvar) at each accepted
* draw.
*
declare vect[rect] goodresp(1000)
declare vect[rect] goodrespa(1000)

declare vector ik a(nvar)
source forcedfactor.src
*
compute accept=0
infobox(action=define,progress,lower=1,upper=n1) 'Monte Carlo Integration'
do draws=1,n1
   *
   * Make a draw from the posterior for the VAR and compute its impulse
   * responses.
   *
   compute sigmad  =%ranwisharti(svt,wishdof)
   compute p       =%decomp(sigmad)
   compute ranc    =%ran(1.0)
   compute betau   =sxx*ranc*tr(p)
   compute betadraw=betaols+betau
   compute %modelsetcoeffs(varmodel,betadraw)
   *
   * This is changed to unit shocks rather than orthogonalized shocks.
   *
   impulse(noprint,model=varmodel,decomp=%identity(6),results=impulses,steps=nstep)
   *
   * Do the subdraws over the unit sphere. These give the weights on the
   * orthogonal components.
   *
   do subdraws=1,n2
   ************************************************
   * First Set of Restrictions - Unique Impulse Vector
   ************************************************
      compute v1=%ran(1.0),v1=v1/sqrt(%normsqr(v1))
      compute v=i1=p*v1
      do k=1,KMAX+1
          compute ik=%xt(impulses,k)*v
          if ik(1)>0.or.ik(2)>0.or.ik(3)>0.or.ik(6)<0
             branch 105
      end do k
      *
      * Meets the first restriction
      * Draw from the orthogonal complement of i1 (last five columns of the
      * factor "f").
      *
      @forcedfactor(force=column) sigmad i1 f
      compute v2=%ran(1.0),v2=v2/sqrt(%normsqr(v2))
      compute v=i2=%xsubmat(f,1,6,2,6)*v2
      *****************************************************
      *  Second Set of Restrictions
      *****************************************************
      do k=1,KMAX+1
         compute ik=%xt(impulses,k)*v
         if ik(4)<0
            branch 105
      end do k
      *
      * Meets both restrictions

*************Get fedrate shock***********************************************
      compute accept=accept+1
      *dim goodrespa(accept)(nstep,nvar)
      dim goodresp(accept)(nstep,nvar)
      ewise goodresp(accept)(i,j)=(ik=%xt(impulses,i)*i1),ik(j)
      *ewise goodrespa(accept)(i,j)=(ik=%xt(impulses,i)*i2),ik(j)
      if accept>=1000
         break
   :105
   end do subdraws
   if accept>=1000
      break
   infobox(current=draws)
end do draws
infobox(action=remove)
*
* Post-processing. Graph the mean of the responses along with the 16% and 84%-iles
*
clear upper lower resp
*
spgraph(vfields=3,hfields=2,hlabel='Impulse Responses with fed shock',subhea='Italy')
do i=1,nvar
   compute minlower=maxupper=0.0
   smpl 1 accept
   do k=1,nstep
      set work = goodresp(t)(k,i)
      compute frac=%fractiles(work,||.16,.84||)
      compute lower(k)=frac(1)
      compute upper(k)=frac(2)
      compute resp(k)=%avg(work)
   end do k
*
   smpl 1 nstep
   graph(pattern,ticks,number=0,picture='##.##',header='Impulse Responses for '+vl(i)) 3
   # resp
   # upper / 2
   # lower / 2
end do i
*
spgraph(done)