*
* Example of some of the calculations done in Pardo and Torro(2007),
* "Trading with asymmetric volatility spillovers", JBFA, vol. 34(9-10),
* 1548-1568.
*
* Applied to different data set (U.S Stocks). These are closing weekly
* adjusted values for the DJIA and the Russell 2000.
*
open data "usstocks.xls"
calendar(w) 2006:1:2
*
* The data were taken from Yahoo Finance and are in reversed
* chronological order. RATS v9.1 automatically detects that and
* rearranges the data. (The REVERSE option has no effect, but creates an
* immediate error with earlier versions of RATS).
*
data(format=xls,org=columns,reverse) / djia russell2000
*
* Normalize both to 100 * log relative to the first entry
*
set logdjia = 100.0*(log(djia/djia(1)))
set logrut  = 100.0*(log(russell2000/russell2000(1)))
*
graph(footer="Log Stock Indexes",key=upleft,klabels=||"DJIA","Russell 2000"||) 2
# logdjia
# logrut
*
* AIC and SBC (or HQ) pick very different lag lengths. The 1 favored by
* SBC and HQ seems too short. AIC has a local minimum at 4, so we choose
* that. (None of these criteria are designed for heteroscedastic data,
* so they can be thrown off by outliers).
*
@varlagselect(lags=10,crit=aic)
# logdjia logrut
@varlagselect(lags=10,crit=sbc)
# logdjia logrut
*
compute nlags=4
*
* There's no strong reason to believe that the two stock indexes are
* cointegrated since they have different risk-return structures. The
* Johansen test marginally rejects cointegration, but we'll go ahead and
* assign cointegration rank one for illustrative purposes. This uses the
* estimated cointegrating vector. In other applications (for instance
* with spot vs futures prices), the cointegrating vector will be known
* theoretically, and you should use that, not an estimated one.
*
* This analyzes cointegration using DET=CONSTANT, which allows for a
* trend in the series.
*
@johmle(lags=nlags,det=constant,cv=cv)
# logdjia logrut
*
* This sets up a VECM for the DET=CONSTANT model (including a CONSTANT
* as a DET term in the VAR). For this, you put in the original dependent
* variables, not the differences. The rearrangement of the regression is
* done internally when you ESTIMATE the model.
*
equation(coeffs=cv) ecteq *
# logdjia logrut
*
system(model=basevecm)
variables logdjia logrut
lags 1 to nlags
det constant
ect ecteq
end(system)
*
estimate(resids=vecmresids)
*
* Note that the number of lags used in the diagnostic tests has nothing
* to do with the choice of the number of lags in the VAR model. (The
* paper does some tests later on the OLS VECM residuals and uses 5 for
* those). Somewhere in the 3-5 lag range generally seems appropriate for
* this.
*
@mvarchtest(lags=5)
# vecmresids
*
* The model form with ECT isn't directly usable in GARCH, so this
* rewrites it in terms of differences and the observed value of the
* cointegrating vector.
*
set dx1 = logdjia-logdjia{1}
set dx2 = logrut-logrut{1}
*
set ect = %eqnprj(ecteq,t)
*
system(model=garchvecm)
var dx1 dx2
Lags 1 to nlags-1
det constant ect{1}
end(system)
*
* Estimate an asymmetric BEKK GARCH model. These saves the jointly
* standardized residuals (as STDU), the conditional covariance matrices
* (HH) and (non-standardized) residuals (RR) and the VECH representation
* (VECHCOMPS).
*
garch(model=garchvecm,mv=bekk,asymmetric,stdresids=stdu,hmatrices=hh,rvectors=rr,$
  pmethod=simplex,piters=20,method=bfgs,iters=500,vechmat=vechcomps)
*
* These are (our) recommended diagnostic tests on the jointly standardized
* residuals. Tests on the univariate standardized residuals come later.
*
@mvarchtest(lags=5)
# stdu
@mvqstat(lags=5)
# stdu
*
* Causality tests in the mean. The first tests does a joint test for
* whether the dx2 lags and the ect lag are zero in the first equation
* (first block of mean model coefficients) and the second does a similar
* test for the dx1 lags and ect in the second block of mean model
* coefficients. (You can use the Statistics-Regression Tests, Exclusion
* Restrictions wizard to set these up.)
*
* (The paper doesn't include the ECT in the test, but it should be
* included since if it's non-zero, the lagged "other" does help to
* predict).
*
test(zeros,title='x2-->x1 causality')
# 4 5 6 8
test(zeros,title='x1-->x2 causality')
# 9 10 11 16
*
* Wald tests of BEKK coefficients
*
test(zeros,title='bekk cross effects')
# 21 22 25 26 29 30
test(zeros,title='asymmetry')
# 28 29 30 31
*
* Tests on the univariate standardized residuals
*
set std1 = rr(t)(1)/sqrt(hh(t)(1,1))
set std2 = rr(t)(2)/sqrt(hh(t)(2,2))
*
* Skewness, (excess) kurtosis and normality (Jarque-Bera) are all part
* of the standard STATISTICS output.
*
stats std1
stats std2
*
* The Q(20) (Ljung-Box) and Q^2(20) (McLeod-Li) are included in the
* BDINDTESTS procedure.
*
@bdindtests(number=20) std1
@bdindtests(number=20) std2
*
* You can also do the LB and McLeod-Li tests separately using @RegCorrs
* (for LB) and @McLeodLi. The results are very slightly different
* because they use a slightly different calculation for the
* autocorrelations.
*
@regcorrs(qstats,number=20) std1
@mcleodli(number=20) std1
@archtest(lags=20) std1
*
@regcorrs(qstats,number=20) std2
@mcleodli(number=20) std2
@archtest(lags=20) std2
*
* These graph the conditional volatility (in standard deviations to
* allow more detail to be seen), the conditional correlation and the
* conditional beta. The latter two also include a horizontal line at the
* corresponding value from the original VECM estimates.
*
set sig11 = sqrt(hh(t)(1,1))
set sig22 = sqrt(hh(t)(2,2))
graph(footer="Conditional Volatility in Std Deviations",$
   key=below,klabels=||"DJIA","Russell 2000"||) 2
# sig11
# sig22
*
set rho12 = %cvtocorr(hh(t))(1,2)
graph(footer="Conditional correlation",vgrid=%cvtocorr(%sigma)(1,2))
# rho12
*
set beta12 = hh(t)(1,2)/hh(t)(1,1)
graph(footer="Conditional beta of Russell 2000 with respect to DJIA",vgrid=%sigma(1,2)/%sigma(1,1))
# beta12
*
* Compute impact surfaces (done as contour graphs)
*
compute ngrid=21
*
* This picks the grids in each direction based upon the observed
* absolute values of the residuals from the VECM (going from - to + of
* the 99%-ile of the absolute value).
*
set plus1 = abs(vecmresids(1))
stats(fract,noprint) plus1
compute grid1=%seqrange(-%fract99,+%fract99,ngrid)
set plus2 = abs(vecmresids(2))
stats(fract,noprint) plus2
compute grid2=%seqrange(-%fract99,+%fract99,ngrid)
*
dec rect var1g(ngrid,ngrid) var2g(ngrid,ngrid) covarg(ngrid,ngrid)
*
do i=1,ngrid
   do j=1,ngrid
      compute [symm] response=%vectosymm($
            vechcomps("A")*%vec(%outerxx(||grid1(i)|grid2(j)||))+$
            vechcomps("D")*%vec(%outerxx(||%minus(grid1(i))|%minus(grid2(j))||)),$
            2)
      compute var1g(i,j)=response(1,1)
      compute var2g(i,j)=response(2,2)
      compute covarg(i,j)=response(1,2)
   end do j
end do i
*
* Do contour plots
*
gcontour(x=grid1,y=grid2,f=var1g,hlabel="DJIA",vlabel="Russell 2000",$
  footer="Impact Surface for DJIA Variance")
gcontour(x=grid1,y=grid2,f=var2g,hlabel="DJIA",vlabel="Russell 2000",$
  footer="Impact Surface for Russell 2000 Variance")
gcontour(x=grid1,y=grid2,f=covarg,hlabel="DJIA",vlabel="Russell 2000",$
  footer="Impact Surface for Covariance")