The RATS Software Forum

Code: Select all

compute qr=%qrdecomp(%ranmat(n,n))

Code: Select all

*
* This program applies the sign restriction algorithm of "Structural Vector Autoregressions: Theory of Identification 
* and Algorithms for Inference" and "Markov-Switching Structural Vector Autoregressions: Theory and Application"
* by J. Rubio-Ramirez, D. Waggoner, and T. Zha to estimate impulse responses as in 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. 
*
* Todd Clark wrote the program by modifying Tom Doan's Uhlig code.
*
* This version omits the computation of the forecast error variance but does inclue an historical decomposition 
* of the data showing the cumulative effects
* of the monetary policy shocks. Note that this is not part of the paper.
*
grparm(portrait) header 15 subheader 12 keylabeling 16 footer 10 vlabel 12 hlabel 12
env nowshowgraphs gsave="Myresults*.rgf" gformat=rgf

open data uhligdata.xls
calendar(m) 1965
data(format=xls,org=columns) 1965:01 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"||
*
* 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=10000
compute n2=10000
compute nkeep=20000
compute nvar=6
compute nstep=60
compute KMAX=5

* set up the number of shocks considered, along with signs
* Note that the computations portion of the program is set up to allow
*  for multiple shocks, but the portion of the program that constructs the charts
*  assumes only one shock of interest
comp nshock = 1
dec rec signs(nvar,nshock)   ;* a matrix with the signs to be imposed, for # shocks = nshock
                              * use a 1 for >= 0, -1 for <= 0, and %NA for unrestricted
comp signs = ||%NA|-1|-1|1|-1|%NA||

*************************************************************************************
compute hstart = 1998:1
compute hend = 2002:12
compute nhor = hend-hstart+1
dec vect[series] upaths(nvar) u(nvar) onesteps(nvar) base(nvar) baseplus(nvar)
clear upaths
*************************************************************************************
*
* 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
*
* Most draws are going to get rejected. We allow for up to nkeep good ones. The
* variable accept will count the number of accepted draws. GOODRESP will be a
* RECT(nvar,nshock) at each accepted draw.
*
declare rect[rect] goodresp(nkeep,nstep) ;* goodfevd(nkeep)
declare vector ik a ones(nvar)
dec rec ikmat
declare series[rect] irfsquared
compute ones=%const(1.0)

*************************************************************************************
declare vector hk(nvar)
declare rect[rect] goodhist(nkeep,nshock)
*************************************************************************************
* set up a vector of selection matrices that correspond to e(i) e(i)' in Tom Doan's notes on hist. decomp.
dec vec[rec] selmat(nshock)
do i = 1,nshock
 comp selmat(i) = %zeros(nvar,nvar)
 comp selmat(i)(i,i) = 1.
end do i

*
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)
  comp [symm] invsigmad = inv(sigmad)
  compute swish =%decomp(sigmad)
  compute betau =%ranmvkron(swish,sxx)
  compute betadraw=betaols+betau
  compute %modelsetcoeffs(varmodel,betadraw)
  **********************************************************************************
  *
  * Use static forecasts over the period for historical decomposition to get the
  * residuals. Use dynamic forecasts over the same period to get the base forecast
  * series for the historical decomposition.
  *
  forecast(model=varmodel,from=hstart,to=hend,results=onestep,errors=u,static)
  forecast(model=varmodel,from=hstart,to=hend,results=base)
  **********************************************************************************
  *
  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
  *
  * Do the subdraws over the unit sphere. These give the weights on the
  * orthogonal components.
  *
  do subdraws = 1,n2
    * first compute the QR decomposition and normalize signs as in RWZ Matlab code
    comp qr=%qrdecomp(%ranmat(nvar,nvar))
    do i = 1,nvar
     if qr(2)(i,i)<0.
      comp %psubmat(qr(1),1,i,-1.*%xcol(qr(1),i))
    end do i
    comp pmat = qr(1)
    *
    * 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
     comp ikmat = %xt(impulses,k)*pmat  ;* matrix of responses of i (row) to shock to j (col)
     do i = 1,nshock  ;* loop across shocks to check signs
      comp ik = %xcol(ikmat,i).*%xcol(signs,i)  ;* if the sign is right, the product should be >= 0
      if %minvalue(ik)<0
       branch 105
     end do i
    end do k
    *
    * This is an accepted draw. Copy the information out. If we have enough
    * good ones, drop out of the overall loop.
    *
    compute accept=accept+1
    dim goodresp(accept,nstep)(nvar,nshock) ;*goodfevd(accept)(nstep,nvar)
    do k = 1,nstep
      comp ikmat = %xt(impulses,k)*pmat  ;* matrix of responses of i (row) to shock to j (col), at horizon k
      comp goodresp(accept,k) = %xsubmat(ikmat,1,nvar,1,nshock)  ;* pull out just responses to shocks of interest (subset of cols)
    end do i
    *ewise goodresp(accept)(i,j)=ik=%xt(impulses,i)*a,ik(j)
    *ewise goodfevd(accept)(i,j)=ik=(irfsquared(i)*(a.^2))./(irfsquared(i)*ones),ik(j)
    *******************************************************************************
    *
    * Compute the historical paths of shocks to the system corresponding to the
    * <<a>> weights on the Choleski factor.
    *
    comp PAmat = swish*pmat
    do k = 1,nshock
     dim goodhist(accept,k)(nhor,nvar)
     compute remap = %mqform(selmat(k),tr(PAmat))*invsigmad ;* this uses the middle term of equation 4 of Tom Doan's notes
                                                             * on historical decompositions, with F = PAmat
     do t=hstart,hend
      compute %pt(upaths,t,remap*%xt(u,t))
     end do t
     *
     * Forecast over the decomposition period with those added shocks
     *
     forecast(model=varmodel,from=hstart,to=hend,results=baseplus,paths)
     # upaths
     *
     * Save the difference between the forecasts with and without and added shocks
     * (and shift them into the slots 1,...,nhor) of <<goodhist>>.
     *
     ewise goodhist(accept,k)(i,j)=baseplus(j)(i+hstart-1)-base(j)(i+hstart-1)
    end do k
    *******************************************************************************
    if accept>=nkeep
     break
  :105
  end do subdraws
  if accept>=nkeep
    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")

do i=1,nvar
 compute minlower=maxupper=0.0
 smpl 1 accept
 do k=1,nstep
  set work = goodresp(t,k)(i,1)   ;* only 1 shock of interest, so we only pull out the first column
  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)

/*
spgraph(vfields=3,hfields=2,hlabel="Figure 8. Fraction of Variance Explained with Pure-Sign Approach")

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(ticks,number=0,min=0.0,picture="##.##",header="Fraction Explained for "+vl(i)) 3
 # resp
 # upper / 2
 # lower / 2
end do i

spgraph(done)
*/

*************************************************************************************

declare series upperH lowerH respH
*
spgraph(vfields=3,hfields=2,hlabel="Figure X. Historical Decomposition with Pure-Sign Approach")

do i=1,nvar
 compute minlower=maxupper=0.0
 smpl hstart hend
 clear upperH lowerH respH
 smpl 1 accept
 do k=hstart,hend
  set work 1 accept = goodhist(t,1)(k-hstart+1,i) ;* only 1 shock of interest, so we only pull out the first column
  compute frac=%fractiles(work,||.16,.50,.84||)
  compute lowerH(k)=frac(1)
  compute upperH(k)=frac(3)
  compute respH(k)=frac(2)
 end do k

 smpl hstart hend
 graph(ticks,dates,picture="##.##",header="Contribution of Monetary Policy shock to "+vl(i)) 3
 # respH hstart hend
 # upperH / 2
 # lowerH / 2
end do i

spgraph(done)

Code: Select all

*******************************************************************************
         ** Forecast over the decomposition period with those added shocks
*******************************************************************************

         forecast(model=varmodel,from=hstart,to=hend,results=baseplus,paths)
         # upaths

         * Save the difference between the forecasts with and without and added shocks
         * (and shift them into the slots 1,...,nhor) of <<goodhist>>.

         ewise goodhist(accept,k)(i,j) = baseplus(j)(i+hstart-1)-base(j)(i+hstart-1)

         end do k

         *******************************************************************************
         if accept>=nkeep
         break
         :105
      end do subdraws

      if accept>=nkeep
        break
      infobox(current=draws)
}
end do draws

infobox(action=remove)