*Initial regressions
*
* Tables IV
*
open data "H:\My Documents\Data1.xls"
calendar(m) 1971:09
data(format=xls,org=columns) 1971:09 2014:02 dbill1 Spread1 rr4 rb1 mb
*
*
*
system(model=onelag)
variables dbill1 Spread1 rr4 rb1
lags 1
det constant
end(system)
*
* Construct a lookup table of ranges. Each row has the start and end of one of
* the ranges of interest.
*
dec rect[integer] rangesMP(5,2)
compute rangesMP=||$
 1985:1,2014:2|$
 1985:1,2007:7|$
 1985:1,2008:9|$
 2007:7,2014:2|$
 2008:9,2014:2||
compute nrangesMP=%rows(rangesMP)
dec rect[integer] ranges(4,2)
compute ranges=||$
 1985:1,2014:2|$
 1985:1,2014:2|$
 1985:1,2014:2|$
 1985:1,2014:2||
compute nranges=%rows(ranges)
compute nstats=9
dec vect stat(nstats) basestat(nstats)
dec vect impact(4) baseimpact(4)
dec vect phi(4)
dec symm statcv(nstats,nstats)
dec symm statcv_mp(4,4)
dec rect statderive
dec rect statderive_mp
*
* Table IV
*
report(action=define)
report(atrow=1,atcol=1,fillby=cols) "VAR Lag Length" "Sample Period" "Return R^2" "Significance"
report(atrow=5,atcol=1,fillby=cols) "Shares of" "Var(xpi)" "" "2Cov(xpi,xr)" "" "2Cov(xpi,xx)" "" $
   "Var(xr)" "" "2Cov(xr,xx)" "" "Var(xx)" ""
report(atrow=18,atcol=1,fillby=cols) "R^2(xpi)" "" "R^2(xr)" "" "R^2(xx)" ""


*
compute n=120
*
* Compute the excess returns on bonds
*
*set eb = -(n-1)/12.0*dbill1-(n-1)/12.0*Spread1
do rangeMP=1,nrangesMP

do range=1,nranges
   compute ncoeffs=4*(4+1)
   dim statderive(nstats,ncoeffs)
   dim statderive_mp(4,ncoeffs)
   report(atrow=1,atcol=range+1,fillby=cols) "1" %datelabel(ranges(range,1))+"-"+%datelabel(ranges(range,2))
   do calc=0,ncoeffs
      if calc==0 {
         estimate(resids=u,coeffs=betafull,noprint) ranges(range,1) ranges(range,2)
         mcov(equation=%modeleqn(onelag,1),nodepvar) ranges(range,1) ranges(range,2) u
         compute xxfull=%mqform(%cmom,%kroneker(%identity(4),%xx))
      }
      else {
         compute fiddle=sqrt(xxfull(calc,calc))*.0001
         compute beta=%vec(betafull)+fiddle*%unitv(ncoeffs,calc)
         compute %modelsetcoeffs(onelag,beta)
      }

      forecast(static,model=onelag,from=ranges(range,1),to=ranges(range,2),errors=u)

      compute a=%modelcompanion(onelag)
      compute discrev=inv(%identity(4)-a)
      *
      *     The reported calculations were done with an incorrect formula.
      *     To match the published results, uncomment out the second definition
      *     of rrevise.
      *
      compute rrevise=%xrow(discrev*(a-a^n),3)
      compute previse=-1.0*rrevise+%xrow(discrev*((n-1)*%identity(4)+discrev*(a^n-a)),1)
      set target ranges(range,1) ranges(range,2) = -(n-1)*(u(1)+u(2))
      set rtilde ranges(range,1) ranges(range,2) = %dot(rrevise,%xt(u,t))
      set ptilde ranges(range,1) ranges(range,2) = %dot(previse,%xt(u,t))
      set xtilde ranges(range,1) ranges(range,2) = -target-rtilde-ptilde
      vcv(nocenter,noprint,matrix=breakdown) ranges(range,1) ranges(range,2)
      # ptilde rtilde xtilde target
      compute breakdown=breakdown/breakdown(4,4)
      compute stat(1)=    breakdown(1,1)
      compute stat(2)=2.0*breakdown(2,1)
      compute stat(3)=2.0*breakdown(3,1)
      compute stat(4)=    breakdown(2,2)
      compute stat(5)=2.0*breakdown(3,2)
      compute stat(6)=    breakdown(3,3)
      compute breakdown=%cvtocorr(breakdown)
      compute stat(7)=    breakdown(1,4)^2
      compute stat(8)=    breakdown(2,4)^2
      compute stat(9)=    breakdown(3,4)^2

****Here we compute the monetary policy impact for given news terms****
***Step 1: I regress the residuals from the forecast VAR on monetary policy proxy (mb in this case), save the coefficients into separate vectors***
LINREG(ROBUST,noprint) U(1) rangesMP(rangeMP,1) rangesMP(rangeMP,2)
# mb Constant
compute [vector] beta1=%beta

LINREG(ROBUST,noprint) U(2) rangesMP(rangeMP,1) rangesMP(rangeMP,2)
# mb Constant
compute [vector] beta2=%beta

LINREG(ROBUST,noprint) U(3) rangesMP(rangeMP,1) rangesMP(rangeMP,2)
# mb Constant
compute [vector] beta3=%beta

LINREG(ROBUST,noprint) U(4) rangesMP(rangeMP,1) rangesMP(rangeMP,2)
# mb Constant
compute [vector] beta4=%beta
compute phi=||beta1(1),beta2(1),beta3(1),beta4(1)||
***This vector will be used in the new derived formulas to compute the impact of the monetary policy shock on the news components of the original decomposition***
compute mp_target=-(n-1)*phi(1)-(n-1)*phi(2)
compute mp_rtilde=%dot(rrevise,phi)
compute mp_ptilde=-1.0*%dot(rrevise,phi)+%dot(%xrow(discrev*((n-1)*%identity(4)+discrev*(a^n-a)),1),phi)
compute mp_xtilde=-mp_target-mp_rtilde-mp_ptilde
compute impact(1)=mp_target
compute impact(2)=mp_rtilde
compute impact(3)=mp_ptilde
compute impact(4)=mp_xtilde

	 if calc==0 {
         do istat=1,9
            report(atrow=2*istat+4,atcol=range+1) stat(istat)

        end do istat
         compute basestat=stat
         compute baseimpact=impact
      }
      else {
         compute %psubmat(statderive,1,calc,(stat-basestat)/fiddle)
         compute %psubmat(statderive_mp,1,calc,(impact-baseimpact)/fiddle)
      }
   end do calc
   compute statcv=%mqform(xxfull,tr(statderive))
   compute statcv_mp=%mqform(xxfull,tr(statderive_mp))
   do istat=1,9
      report(atrow=2*istat+5,atcol=range+1,special=parens) sqrt(statcv(istat,istat))
      end do istat

end do range
report(action=format,atcol=2,align=decimal,picture="*.###")
*report(action=show)

report(action=define)
report(atrow=1,atcol=1,fillby=cols) "MB" " " "Parameter" " " "SE"
report(atrow=1,atcol=2,fillby=rows) "Target" "" "R" "" "P" "" "X" %datelabel(rangesMP(rangeMP,1))+"-"+%datelabel(rangesMP(rangeMP,2))
report(atrow=3,atcol=2) baseimpact(1)
report(atrow=3,atcol=4) baseimpact(2)
report(atrow=3,atcol=6) baseimpact(3)
report(atrow=3,atcol=8) baseimpact(4)
report(atrow=5,atcol=2,special=parens) sqrt(statcv_mp(1,1))
report(atrow=5,atcol=4,special=parens) sqrt(statcv_mp(2,2))
report(atrow=5,atcol=6,special=parens) sqrt(statcv_mp(3,3))
report(atrow=5,atcol=8,special=parens) sqrt(statcv_mp(4,4))
report(action=format,atcol=2,align=decimal,picture="*.###")
report(action=show)
end do rangeMP