The RATS Software Forum

Code: Select all

*
* Replication File for Uhlig(2005), "What are the effects of monetary policy on
* output? Results from an agnostic identification procedure," Journal of Monetary
* Economics, 52, pp 381-419. Penalty function approach.
*
open data uhligdata.xls
calendar(m) 1965
data(format=xls,org=columns) 1965:1 2003:12 gdpc1 gdpdef cprindex totresns bognonbr fedfunds
*
set gdpc1    = log(gdpc1)*100.0
set gdpdef   = log(gdpdef)*100.0
set cprindex = log(cprindex)*100.0
set totresns = log(totresns)*100.0
set bognonbr = log(bognonbr)*100.0
*
system(model=varmodel)
variables gdpc1 gdpdef cprindex fedfunds bognonbr totresns
lags 1 to 12
end(system)
estimate(noprint)
*
dec vect[strings] vl(6)
compute vl=||"real GDP","GDP price defl","Comm. Price Ind.",$
             "Fed Funds Rate","Nonborr. Reserv.","Total Reserves"||
*
* Get the scale factors to be used for impulse responses
*
dec vect scales(6)
dec real func
compute fill=0
dofor i = gdpc1 gdpdef cprindex fedfunds bognonbr totresns
   diff i / diff1
   stats(noprint) diff1
   compute fill=fill+1,scales(fill)=sqrt(%variance)
end dofor i
*
* ndraws is the desired number of draws from the posterior of the VAR
* nvar is the number of variables
* nstep is the number of IRF steps to compute
* KMAX is the "K" value for the number of steps constrained
*
*
compute ndraws=1000
compute nvar=6
compute nstep=60
compute KMAX=5
declare vect[rect] goodresp(ndraws)
declare vector ik a(nvar)
*
dec vect g(nvar-1)
compute g=%fill(nvar-1,1,1.0)
nonlin g
*
* This is the standard setup for MC integration of an OLS VAR
*
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)
*
infobox(action=define,progress,lower=1,upper=ndraws) "Monte Carlo Integration"
*
* Plan for double the number of draws to allow for the ones which we
* reject.
*
compute accept=0
do draws=1,ndraws*2
   *
   * Make a draw from the posterior for the VAR and compute its impulse
   * responses.
   *
   compute sigmad  =%ranwisharti(svt,wishdof)
   compute swish   =%decomp(sigmad)
   compute ranc    =%ran(1.0)
   compute betau   =sxx*ranc*tr(swish)
   compute betadraw=betaols+betau
   compute %modelsetcoeffs(varmodel,betadraw)
   impulse(noprint,model=varmodel,decomp=swish,results=impulses,steps=nstep)
   *
   * Minimize the penalty function, starting from the last set of minimizers
   *
   find(noprint) min func
      compute a=%stereo(g)
      compute func=0.0
      do k=1,KMAX+1
         compute ik=(%xt(impulses,k)*a)./scales
         compute func=func+%if(ik(4)>0,-ik(4),-100*ik(4))+$
                           %if(ik(3)<0,ik(3),100*ik(3))+$
                           %if(ik(2)<0,ik(2),100*ik(2))+$
                           %if(ik(5)<0,ik(5),100*ik(5))
      end do k
   end find
   compute testfunc=func,testbeta=%beta
   *
   * Try the minimization again, starting from a standard set of values
   *
   compute g=%fill(nvar-1,1,1.0)
   find(noprint) min func
      compute a=%stereo(g)
      compute func=0.0
      do k=1,KMAX+1
         compute ik=(%xt(impulses,k)*a)./scales
         compute func=func+%if(ik(4)>0,-ik(4),-100*ik(4))+$
                           %if(ik(3)<0,ik(3),100*ik(3))+$
                           %if(ik(2)<0,ik(2),100*ik(2))+$
                           %if(ik(5)<0,ik(5),100*ik(5))
      end do k
   end find
   *
   * If the two estimates don't match, reject the draw. If they do, copy out the
   * impulse responses.
   *
   if abs(testfunc-func)<.01 {
      compute accept=accept+1
      dim goodresp(accept)(nstep,nvar)
      ewise goodresp(accept)(i,j)=ik=%xt(impulses,i)*a,ik(j)
   }
   *
   * If we've hit out desired number of accepted draws, break
   *
   if accept>=ndraws
      break
   infobox(current=accept)
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="Figure 14. Impulse Responses with Penalty Function Approach")
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(ticks,number=0,picture="##.##",header="Impulse Responses for "+vl(i)) 3
   # resp
   # upper / 2
   # lower / 2
end do i
*
spgraph(done)

Code: Select all

DEC SERIES[RECT] irfsquared

Code: Select all

*
*  Replication File for Uhlig (2005), "What are the effects of monetary policy on output?
*  Results from an agnostic identification procedure." Journal of Monetary Economics, 52, pp
*  381-419. Pure sign restriction approach
*
open data levelsa1.xls
calendar 1970 1 4
data(format=xls,org=columns) 1970:01 2008:03 lrgspen1 lrntax lrgdp pgdpch rint lrpcon lrpinv
*

system(model=varmodel)
variables lrgspen1 lrntax lrgdp pgdpch rint lrpcon
lags 1 to 4
end(system)
estimate(noprint)
*
dec vect[strings] vl(6)
compute vl=||'Govt spending','net tax','real gdp','GDP deflator','real interest rate','Private onsumption'||
*
* n1 is the number of draws from the posterior of the VAR
* n2 is the number of draws from the unit sphere for each draw for the VAR
* nvar is the number of variables
* nstep is the number of IRF steps to compute
* KMAX is the "K" value for the number of steps constrained
*
compute n1=200
compute n2=200
compute nvar=6
compute nstep=60
compute KMAX=5
*
* This is the standard setup for MC integration of an OLS VAR
*
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 vector ik a(nvar)
*
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 swish   =%decomp(sigmad)
   compute ranc    =%ran(1.0)
   compute betau   =sxx*ranc*tr(swish)
   compute betadraw=betaols+betau
   compute %modelsetcoeffs(varmodel,betadraw)
   impulse(noprint,model=varmodel,decomp=swish,results=impulses,steps=nstep)
   *
   * Do the subdraws over the unit sphere. These give the weights on the
   * orthogonal components.
   *
   do subdraws=1,n2
      compute rperp=$%perp(||1.0,-1.0,0.0,0.0,0.0,0.0||*swish)
      compute a=rperp*%ranmat(5,1),a=a/sqrt(%normsqr(a))
      *
      * Check that the responses have the correct signs for steps 1 to KMAX+1
      * (+1 because in the paper, the steps used 0-based subscripts, rather than
      * the 1 based subscripts used by RATS).
      *
      do k=1,KMAX+1
         compute ik=%xt(impulses,k)*a
         if ik(1)<0
         branch 105
      end do k
      *
      * This is an accepted draw. Copy the information out. If we have our 1000
      * good ones, drop out of the overall loop.
      *
      compute accept=accept+1
      dim goodresp(accept)(nstep,nvar)
      ewise goodresp(accept)(i,j)=ik=%xt(impulses,i)*a,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='Figure 6. Impulse Responses with Pure-Sign Approach',subhea='monetary Shock')
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(ticks,number=0,picture='##.##',header='Impulse Responses for '+vl(i)) 3
   # resp
   # upper / 2
   # lower / 2
end do i
*
spgraph(done)

Code: Select all

*
*  Replication File for Uhlig (2005), "What are the effects of monetary policy on output?
*  Results from an agnostic identification procedure." Journal of Monetary Economics, 52, pp
*  381-419. Pure sign restriction approach
*
open data levelsa1.xls
calendar 1970 1 4
data(format=xls,org=columns) 1970:01 2008:03 lrgspen1 lrntax lrgdp pgdpch rint lrpcon lrpinv
*

system(model=varmodel)
variables lrgspen1 lrntax lrgdp pgdpch rint lrpcon
lags 1 to 4
end(system)
estimate(noprint)
*
dec vect[strings] vl(6)
compute vl=||'Govt spending','net tax','real gdp','GDP deflator','real interest rate','Private onsumption'||
*
* n1 is the number of draws from the posterior of the VAR
* n2 is the number of draws from the unit sphere for each draw for the VAR
* nvar is the number of variables
* nstep is the number of IRF steps to compute
* KMAX is the "K" value for the number of steps constrained
*
compute n1=200
compute n2=200
compute nvar=6
compute nstep=60
compute KMAX=5
*
* This is the standard setup for MC integration of an OLS VAR
*
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 vector ik a(nvar)
*
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 swish   =%decomp(sigmad)
   compute ranc    =%ran(1.0)
   compute betau   =sxx*ranc*tr(swish)
   compute betadraw=betaols+betau
   compute %modelsetcoeffs(varmodel,betadraw)
   impulse(noprint,model=varmodel,decomp=swish,results=impulses,steps=nstep)
   *
   * Do the subdraws over the unit sphere. These give the weights on the
   * orthogonal components.
   *
   do subdraws=1,n2
      compute rperp=%perp(||1.0,-1.0,0.0,0.0,0.0,0.0||*swish)
      compute a=rperp*%ranmat(5,1),a=a/sqrt(%normsqr(a))
      *
      * Check that the responses have the correct signs for steps 1 to KMAX+1
      * (+1 because in the paper, the steps used 0-based subscripts, rather than
      * the 1 based subscripts used by RATS).
      *
      do k=1,KMAX+1
         compute ik=%xt(impulses,k)*a
         if ik(1)<0
         branch 105
      end do k
      *
      * This is an accepted draw. Copy the information out. If we have our 1000
      * good ones, drop out of the overall loop.
      *
      compute accept=accept+1
      dim goodresp(accept)(nstep,nvar)
      ewise goodresp(accept)(i,j)=ik=%xt(impulses,i)*a,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='Figure 6. Impulse Responses with Pure-Sign Approach',subhea='monetary Shock')
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(ticks,number=0,picture='##.##',header='Impulse Responses for '+vl(i)) 3
   # resp
   # upper / 2
   # lower / 2
end do i
*
spgraph(done)

Code: Select all

    OPEN DATA uhligdata.xls
    CALENDAR 1965 1 12
    compute missc=1.0e+32
    data(format=xlS,org=columns) 1965:1 2003:12 gdpc1 gdpdef cprindex totresns bognonbr fedfunds

    set gdpc1    = log(gdpc1)*100.0
    set gdpdef   = log(gdpdef)*100.0
    set cprindex = log(cprindex)*100.0
    set totresns = log(totresns)*100.0
    set bognonbr = log(bognonbr)*100.0
    *
    system(model=varmodel)
    variables gdpc1 gdpdef cprindex fedfunds bognonbr totresns
    lags 1 to 12
    end(system)
    estimate(noprint,resid=resids)

    dec vect[strings] vl(6)
    compute vl=||'real GDP','GDP price defl','Comm. Price Ind.',$
         'Fed Funds Rate','Nonborr. Reserv.','Total Reserves'||

    compute n1=4000
    compute n2=4000
    compute nvar=6
    compute nstep=60
    compute KMAX=5

    * 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) goodfevd(1000)
    declare vect[rect] goodrespa(1000) goodfevda(1000)

    declare vector ik a(nvar) ones(nvar)
    declare series[rect] irfsquared
    compute ones=%const(1.0)
    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)
       gset irfsquared 1 1     = %xt(impulses,t).^2
       gset irfsquared 2 nstep = irfsquared{1}+%xt(impulses,t).^2
       *
       * 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 - Demand Shock
          *****************************************************
          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 both restrictions
          *
          compute accept=accept+1
          dim goodrespa(accept)(nstep,nvar) goodfevda(accept)(nstep,nvar)
          dim goodresp(accept)(nstep,nvar) goodfevd(accept)(nstep,nvar)
          ewise goodresp(accept)(i,j)=(ik=%xt(impulses,i)*i1),ik(j)
          ewise goodfevd(accept)(i,j)=(ik=(irfsquared(i)*(i1.^2))./(irfsquared(i)*ones)),ik(j)
          ewise goodrespa(accept)(i,j)=(ik=%xt(impulses,i)*i2),ik(j)
          ewise goodfevda(accept)(i,j)=(ik=(irfsquared(i)*(i2.^2))./(irfsquared(i)*ones)),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 Pure-Sign Approach',subhea='Subhea')
    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)
          display resp(k)
       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)
    smpl
    *
    clear upper lower resp
    *
    spgraph(vfields=3,hfields=2,hlabel='SECOND Impulse Responses with Pure-Sign Approach',subhea='Subhea')
    do i=1,nvar
       compute minlower=maxupper=0.0
       smpl 1 accept
       do k=1,nstep
          set work = goodrespa(t)(k,i)
          compute frac=%fractiles(work,||.16,.84||)
          compute lower(k)=frac(1)
          compute upper(k)=frac(2)
          compute resp(k)=%avg(work)
          display resp(k)
       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)
    smpl
    *
    clear upper lower resp
    *
    spgraph(vfields=3,hfields=2,hlabel='Fraction of Variance Explained with Pure-Sign Approach',subhea='Subhea')
    do i=1,nvar
       compute minlower=maxupper=0.0
       smpl 1 accept
       do k=1,nstep
          set work = goodfevd(t)(k,i)
          compute frac=%fractiles(work,||.16,.50,.84||)
          compute lower(k)=frac(1)
          compute upper(k)=frac(3)
          compute resp(k)=frac(2)
          display resp(k)
       end do k
    *
       smpl 1 nstep
       graph(pattern,ticks,number=0,min=0.0,picture="##.##",header="Fraction Explained for "+vl(i)) 3
       # resp
       # upper / 2
       # lower / 2
    end do i
    *
    spgraph(done)
    smpl
    *
    clear upper lower resp
    *
    spgraph(vfields=3,hfields=2,hlabel='SECOND Fraction of Variance Explained with Pure-Sign Approach',subhea='Subhea')
    do i=1,nvar
       compute minlower=maxupper=0.0
       smpl 1 accept
       do k=1,nstep
          set work = goodfevda(t)(k,i)
          compute frac=%fractiles(work,||.16,.50,.84||)
          compute lower(k)=frac(1)
          compute upper(k)=frac(3)
          compute resp(k)=frac(2)
          display resp(k)
       end do k
    *
       smpl 1 nstep
       graph(pattern,ticks,number=0,min=0.0,picture="##.##",header="Fraction Explained for "+vl(i)) 3
       # resp
       # upper / 2
       # lower / 2
    end do i
    *
    spgraph(done)
    smpl

Code: Select all

   impulse(noprint,model=varmodel,decomp=swish,results=impulses,steps=nstep)
   compute step0=%xt(impulses,1)
   compute step1=%xt(impulses,2)
   compute step2=%xt(impulses,3)
   compute step3=%xt(impulses,4)
   dec rect r(4,6)
   compute %psubmat(r,1,1,||1.0,-1.0,0.0,0.0,0.0,0.0||*step0)
   compute %psubmat(r,2,1,||1.0,-1.0,0.0,0.0,0.0,0.0||*step1)
   compute %psubmat(r,3,1,||1.0,-1.0,0.0,0.0,0.0,0.0||*step2)
   compute %psubmat(r,4,1,||1.0,-1.0,0.0,0.0,0.0,0.0||*step3)
   compute rfperp=%perp(r)
   compute alpha=rfperp*%ranmat(%cols(rfperp),1)
   compute alpha=alpha/sqrt(%normsqr(alpha))

Code: Select all

   impulse(noprint,model=varmodel,decomp=swish,results=impulses,steps=nstep)
   gset irfsquared 1 1     = %xt(impulses,t).^2
   gset irfsquared 2 nstep = irfsquared{1}+%xt(impulses,t).^2

Code: Select all

      ewise goodfevd(accept)(i,j)=ik=(irfsquared(i)*(a.^2))./(irfsquared(i)*ones),ik(j)

Code: Select all

   impulse(noprint,model=varmodel,decomp=%identity(nvar),results=impulses,steps=nstep)

Code: Select all

   gset irfsquared 1 1     = %xt(impulses,t).^2
   gset irfsquared 2 nstep = irfsquared{1}+%xt(impulses,t).^2

Code: Select all

       compute v1=%ransphere(6)
       compute i1=p*v1
       do k=1,KMAX+1
           compute ik=%xt(impulses,k)*v1
           if ik(1)>0.or.ik(2)>0.or.ik(3)<0.or.ik(6)>0
              branch 105
       end do k

Code: Select all

       @forcedfactor(force=column) sigmad i1 f
       compute v2=%ransphere(5)
       compute i2=%xsubmat(f,1,6,2,6)*v2
       compute v2=inv(p)*i2
       *****************************************************
       *  Second Set of Restrictions - Demand Shock
       *****************************************************
       do k=1,KMAX+1
          compute ik=%xt(impulses,k)*v2
          if ik(1)<0.or.ik(2)<0.or.ik(3)<0.or.ik(6)<0
             branch 105
       end do k

Code: Select all

       ewise goodfevd(accept)(i,j)=(ik=(irfsquared(i)*(v1.^2))./(irfsquared(i)*ones)),ik(j)
       ewise goodfevda(accept)(i,j)=(ik=(irfsquared(i)*(v2.^2))./(irfsquared(i)*ones)),ik(j)

Code: Select all

 OPEN DATA uhligdata.xls
 CALENDAR 1965 1 12
 compute missc=1.0e+32
 data(format=xlS,org=columns) 1965:1 2003:12 gdpc1 gdpdef cprindex totresns bognonbr fedfunds

 set gdpc1    = log(gdpc1)*100.0
 set gdpdef   = log(gdpdef)*100.0
 set cprindex = log(cprindex)*100.0
 set totresns = log(totresns)*100.0
 set bognonbr = log(bognonbr)*100.0
 *
 system(model=varmodel)
 variables gdpc1 gdpdef cprindex fedfunds bognonbr totresns
 lags 1 to 12
 end(system)
 estimate(noprint,resid=resids)

 dec vect[strings] vl(6)
 compute vl=||'real GDP','GDP price defl','Comm. Price Ind.',$
      'Fed Funds Rate','Nonborr. Reserv.','Total Reserves'||

 compute n1=1000
 compute n2=1000
 compute nvar=6
 compute nstep=60
 compute KMAX=5

 * 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) goodfevd(1000)
 declare vect[rect] goodrespa(1000) goodfevda(1000)

 declare vector ik a(nvar) ones(nvar)
 declare series[rect] irfsquared
 compute ones=%const(1.0)
 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=p,results=impulses,steps=nstep)
    gset irfsquared 1 1     = %xt(impulses,t).^2
    gset irfsquared 2 nstep = irfsquared{1}+%xt(impulses,t).^2
    *
    * 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=%ransphere(6)
       compute i1=p*v1
       do k=1,KMAX+1
           compute ik=%xt(impulses,k)*v1
           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=%ransphere(5)
       compute i2=%xsubmat(f,1,6,2,6)*v2
       compute v2=inv(p)*i2
       *****************************************************
       *  Second Set of Restrictions - Demand Shock
       *****************************************************
       do k=1,KMAX+1
          compute ik=%xt(impulses,k)*v2
          if ik(1)<0.or.ik(2)<0.or.ik(3)<0.or.ik(6)<0
             branch 105
       end do k
       *
       * Meets both restrictions
       *
       compute accept=accept+1
       dim goodrespa(accept)(nstep,nvar) goodfevda(accept)(nstep,nvar)
       dim goodresp(accept)(nstep,nvar) goodfevd(accept)(nstep,nvar)
       ewise goodresp(accept)(i,j)=(ik=%xt(impulses,i)*v1),ik(j)
       ewise goodfevd(accept)(i,j)=(ik=(irfsquared(i)*(v1.^2))./(irfsquared(i)*ones)),ik(j)
       ewise goodrespa(accept)(i,j)=(ik=%xt(impulses,i)*v2),ik(j)
       ewise goodfevda(accept)(i,j)=(ik=(irfsquared(i)*(v2.^2))./(irfsquared(i)*ones)),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 Pure-Sign Approach',subhea='Subhea')
 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)
 smpl
 *
 clear upper lower resp
 *
 spgraph(vfields=3,hfields=2,hlabel='SECOND Impulse Responses with Pure-Sign Approach',subhea='Subhea')
 do i=1,nvar
    compute minlower=maxupper=0.0
    smpl 1 accept
    do k=1,nstep
       set work = goodrespa(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)
 smpl
 *
 clear upper lower resp
 *
 spgraph(vfields=3,hfields=2,hlabel='Fraction of Variance Explained with Pure-Sign Approach',subhea='Subhea')
 do i=1,nvar
    compute minlower=maxupper=0.0
    smpl 1 accept
    do k=1,nstep
       set work = goodfevd(t)(k,i)
       compute frac=%fractiles(work,||.16,.50,.84||)
       compute lower(k)=frac(1)
       compute upper(k)=frac(3)
       compute resp(k)=frac(2)
    end do k
 *
    smpl 1 nstep
    graph(pattern,ticks,number=0,min=0.0,picture="##.##",header="Fraction Explained for "+vl(i)) 3
    # resp
    # upper / 2
    # lower / 2
 end do i
 *
 spgraph(done)
 smpl
 *
 clear upper lower resp
 *
 spgraph(vfields=3,hfields=2,hlabel='SECOND Fraction of Variance Explained with Pure-Sign Approach',subhea='Subhea')
 do i=1,nvar
    compute minlower=maxupper=0.0
    smpl 1 accept
    do k=1,nstep
       set work = goodfevda(t)(k,i)
       compute frac=%fractiles(work,||.16,.50,.84||)
       compute lower(k)=frac(1)
       compute upper(k)=frac(3)
       compute resp(k)=frac(2)
    end do k
 *
    smpl 1 nstep
    graph(pattern,ticks,number=0,min=0.0,picture="##.##",header="Fraction Explained for "+vl(i)) 3
    # resp
    # upper / 2
    # lower / 2
 end do i
 *
 spgraph(done)
 smpl