Other interesting models

References

Introduction

In this practical, we will be looking at some applications of GARCH models mapped onto the multivariate plane. Particularly, we will be considering applications to the parsimonious estimation of conditional correlation estimates that can be estimated for a matrix of n-series.

Comovement is of particular interest to investors and investment managers, as it provides insight into diversification potential. The higher the correlation between two assets, the lower the diversification. In this practical, we will outline some of the most widely used techniques to calculate this measure. We will be looking at the:

DCC-GARCH model
ADCC-GARCH model
GO-GARCH model

We can then use the time-varying correlation estimates so defined to provide insight into the underlying comovement structures of a portfolio of assets. We can also use our estimates to forecast the expected future comovement between two assets.

We will be using mainly two packages: MTS from the Tsay (2014) textbook (where the MV volatility models are discussed in chapter 7), and Ghalanos (2014) rmgarch, which is the multivariate extension of the rugarch package. MTS provides us with excellent goodness-of-fit estimates, while rmgarch provides unmatched ease of estimation and forecasting and plotting functionality.

See the brilliant summary papers by Bauwens, Laurent, and Rombouts (2006) and Silvennoinen and Teräsvirta (2009), as well as the canonical paper on DCC modelling by Engle (2002). Various other applications of these MV-GARCH techniques have been written (c.f. Horvath and Poldauf (2012), Christiansen (2007), Johansson (2008), etc.).

Import data

TOday we will be using the same data as in the previous practical to keep things simple. I will provide you with much more in-depth data for the projects in due course. For now, load and calculate and clean the returns similar to last practical:

## No+te: we will install the MTS package and rmgarch today
pacman::p_load("MTS", "robustbase")
pacman::p_load("tidyverse", "devtools", "rugarch", "rmgarch", 
    "forecast", "tbl2xts", "lubridate", "PerformanceAnalytics", 
    "ggthemes")
dailydata <- fmxdat::findata

Calculating Returns and Cleaning

rtn <- dailydata %>% gather(Tickers, TRI, -Date) %>% group_by(Tickers) %>% 
    mutate(dlogret = log(TRI) - log(lag(TRI))) %>% mutate(scaledret = (dlogret - 
    mean(dlogret, na.rm = T))) %>% filter(Date > dplyr::first(Date)) %>% 
    ungroup() %>% mutate(Tickers = gsub("JSE\\.", "", Tickers)) %>% 
    mutate(Tickers = gsub("\\.Close", "", Tickers))

So now we have our return data, rtn, with which we will conduct our MV-Volatility modelling.

MV Conditional Heteroskedasticity tests

Tsay (2014) summarizes how to fit MV Portmanteau tests in chapter 7.1, which motivates the use of heteroskedasticity models such as GARCH. See e.g. the MarchTest from MTS package (Tsay 2014, 401:407)

xts_rtn <- rtn %>% tbl_xts(., cols_to_xts = "dlogret", spread_by = "Tickers")
MarchTest(xts_rtn)

## Q(m) of squared series(LM test):  
## Test statistic:  4225.995  p-value:  0 
## Rank-based Test:  
## Test statistic:  3326.733  p-value:  0 
## Q_k(m) of squared series:  
## Test statistic:  4017.933  p-value:  0 
## Robust Test(5%) :  2101.564  p-value:  0

The MARCH test indicates that all the MV portmanteau tests reject the null of no conditional heteroskedasticity, motivating our use of MVGARCH models.

Now there are very many types of multi-variate volatility models that we can fit. Here follows a few that we can use to control for the remaining serial heteroskedasticity in our series- as well as uncover the conditional covariance structure to do further analyses on.

MV modelling: Conceptually

Generalization of univariate GARCH models to the multivariate domain is a conceptually simple exercise: Given the stochastic process, of financial returns with dimension and mean vector , given the information set :

where the residuals of the process are modelled as:

above is a positive definite matrix such that is the conditional covariance matrix of .

is a i.i.d. N(0,1) series.

Now several techniques have been proposed that map the matrix into the multivariate plain:

VECH, BEKK, CCC-GARCH, DCC-GARCH, GOGARCH, etc.
A large literature has developed that attempts to map into the multivariate space,
- But MV Volatility models typically suffer from curse of high dimensionality..
Studying correlation between a large amount of series, implies a potentially very large set of unique bivariate correlation combinations.
- This requires a parsimonious and effective means of isolating the conditional comovements.

EWMA Models

The EWMA offers a simplified means of estimating MVGARCH models.

From our earlier theoretical definitions of the residuals,

, we define the EWMA model as:

where denotes the decaying rate, or the persistence parameter (Tsay 2014, p414). In practice, the value of is typically fixed at around 0.96, and not estimated. This follows simply to avoid complicated and long estimation times to reach a value that is anyway quite similar. This makes the model extremely parsimonious and easy to estimate. The problem, though, is that this fit is assumed to hold across many time-series (something which has been shown in the literature to be an inaccurate assumption).

To fit the EWMA on our data, let’s use the MTS package.

Before we fit the data, however, let’s follow Tsay (2014) in fitting a VAR(1) mean model to control for remaining serial auto-persistence in the first moment.

After controlling for MV first order persistence, let’s use a EWMA first to control for MV second moment persistence (GARCH effects).

Note: for the EWMA estimation, a negative spec for lambda implies estimating lambda, while a positive value sets it. Thus, we could, e.g., to save estimation time simply set lambda to 0.96.

rtnv <- VAR(x = xts_rtn, p = 1)

## Constant term: 
## Estimates:  0.0001411384 0.0005675912 0.001094744 0.0004293827 0.0002673927 0.0004502421 0.0002400217 0.0004910471 
## Std.Error:  0.0003186641 0.0002862482 0.0003079427 0.0003310735 0.0003061103 0.0003322034 0.0003116739 0.0003038726 
## AR coefficient matrix 
## AR( 1 )-matrix 
##         [,1]     [,2]     [,3]    [,4]    [,5]     [,6]    [,7]     [,8]
## [1,] -0.0577 -0.02485 -0.03328  0.0447 -0.0238 -0.04446  0.0610  0.05517
## [2,]  0.0422 -0.10849  0.00937  0.0487 -0.0475 -0.01781  0.0571 -0.02540
## [3,]  0.0268 -0.01156 -0.06206  0.0468 -0.0316  0.00529  0.0252  0.02583
## [4,]  0.0817 -0.04293 -0.01325 -0.0600  0.0133  0.01197  0.0271 -0.01124
## [5,]  0.0649 -0.00593 -0.02121  0.0196 -0.1134 -0.02768  0.0649  0.00776
## [6,]  0.1234 -0.03326 -0.00413  0.1853  0.0140 -0.32409  0.0284  0.00145
## [7,]  0.0559 -0.02498 -0.01339  0.0858 -0.0125 -0.04049 -0.0801 -0.01386
## [8,]  0.0348 -0.01500 -0.01612  0.0937  0.0054 -0.07034  0.0461 -0.10339
## standard error 
##        [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]
## [1,] 0.0272 0.0234 0.0185 0.0332 0.0280 0.0306 0.0299 0.0236
## [2,] 0.0245 0.0210 0.0166 0.0298 0.0252 0.0275 0.0268 0.0212
## [3,] 0.0263 0.0226 0.0178 0.0320 0.0271 0.0296 0.0288 0.0228
## [4,] 0.0283 0.0243 0.0192 0.0345 0.0291 0.0318 0.0310 0.0245
## [5,] 0.0262 0.0225 0.0177 0.0319 0.0269 0.0294 0.0287 0.0226
## [6,] 0.0284 0.0244 0.0192 0.0346 0.0292 0.0319 0.0311 0.0246
## [7,] 0.0266 0.0229 0.0181 0.0324 0.0274 0.0300 0.0292 0.0231
## [8,] 0.0260 0.0223 0.0176 0.0316 0.0267 0.0292 0.0285 0.0225
##   
## Residuals cov-mtx: 
##              ABG          BVT          CPI          FSR          NED
## ABG 0.0003598300 1.591799e-04 1.012735e-04 0.0002651588 2.404987e-04
## BVT 0.0001591799 2.903466e-04 7.760638e-05 0.0001757061 1.563139e-04
## CPI 0.0001012735 7.760638e-05 3.360245e-04 0.0001127613 9.796866e-05
## FSR 0.0002651588 1.757061e-04 1.127613e-04 0.0003884005 2.497708e-04
## NED 0.0002404987 1.563139e-04 9.796866e-05 0.0002497708 3.320373e-04
## RMH 0.0002520055 1.793830e-04 1.065675e-04 0.0003316892 2.386856e-04
## SBK 0.0002555962 1.645707e-04 9.843702e-05 0.0002695124 2.421150e-04
## SLM 0.0001879859 1.559992e-04 9.259939e-05 0.0002098443 1.846055e-04
##              RMH          SBK          SLM
## ABG 0.0002520055 2.555962e-04 1.879859e-04
## BVT 0.0001793830 1.645707e-04 1.559992e-04
## CPI 0.0001065675 9.843702e-05 9.259939e-05
## FSR 0.0003316892 2.695124e-04 2.098443e-04
## NED 0.0002386856 2.421150e-04 1.846055e-04
## RMH 0.0003910563 2.638572e-04 2.111537e-04
## SBK 0.0002638572 3.442167e-04 1.925607e-04
## SLM 0.0002111537 1.925607e-04 3.272006e-04
##   
## det(SSE) =  1.180731e-30 
## AIC =  -68.87554 
## BIC =  -68.7647 
## HQ  =  -68.83602

et <- rtnv$residuals  ## Save the VAR(1) model's residuals.
# Now do a GARCH test on remaining series heteroskedasticity:
MarchTest(et)

## Q(m) of squared series(LM test):  
## Test statistic:  3843.756  p-value:  0 
## Rank-based Test:  
## Test statistic:  3316.768  p-value:  0 
## Q_k(m) of squared series:  
## Test statistic:  3948.179  p-value:  0 
## Robust Test(5%) :  2106.163  p-value:  0

# The MarchTest shows the presence of MV heteroskedasticity,
# and motivates the use of EWMA (or any other MV vol model).

# Let's fit the EWMA:
EWMA <- EWMAvol(et, lambda = 0.96)

# The available output:
names(EWMA)

## [1] "Sigma.t" "return"  "lambda"

# gives the conditional variance estimates:
Sigma.t <- EWMA$Sigma.t

# The mean returns series (which is the same as et above)
all.equal(EWMA$return, et)

## [1] TRUE

# And the estimate of lambda.

# Model checking could be done using the built-in MCHdiag
# from Tsay (2014):
modcheck <- MCHdiag(et, Sigma.t)

## Test results:  
## Q(m) of et: 
## Test and p-value:  184.9199 0 
## Rank-based test: 
## Test and p-value:  192.7904 0 
## Qk(m) of epsilon_t: 
## Test and p-value:  1006.7 0 
## Robust Qk(m):  
## Test and p-value:  897.4593 6.503142e-11

Note the MARCH test (section 7.1 in Tsay (2014)) p-values are all close to zero, indicating the standardized residuals (zt) of EWMA still have remaining conditional heteroskedasticity.

Let’s progress to a more generalized estimation model: BEKK

BEKK Models

The BEKK GARCH family of models (Engle and Kroner (1995)), as discussed in class, allow us to consider spill-over effects and study the impact of one variable on another.

We won’t spend too much time on this technique, as it typically suffers from very high dimensionality. It is mostly used for studying spill-overs in a constrained framework.

See Tsay (2014) p.417 for code to run this - health warning: it will take long to run on this data set…

DCC Models

DCC models offer a simple and more parsimonious means of doing MV-vol modelling. In particular, it relaxes the constraint of a fixed correlation structure (assumed by the CCC model), to allow for estimates of time-varying correlation (see Engle (2002) for the original paper).

The DCC model can be defined as:

Equation splits the varcovar matrix into identical diagonal matrices and an estimate of the time-varying correlation. Estimating requires it to be inverted at each estimated period, and thus a proxy equation is used (Engle (2002)):

Note the above equation is similar in form to a GARCH(1,1) process, with non-negative scalars and , and with:
- the unconditional (sample) variance estimate between series and ,
- the unconditional matrix of standardized residuals from each univariate pair estimate.}

We next use eq. to estimate as:

Which has bivariate elements:

The resulting DCC model is then formulated as:

Next, the ADCC allows for leverage effects in the above (think GJR GARCH..)

Now fitting these techniques to our earlier data implies a two-step approach:

Fitting univariate GARCH models to each of our VAR series’ residuals: . The volatility approximation series so estimated, , will then be used in step 2.
These vol series are then standardized and used in fitting a DCC model for .

# dccPre fits the univariate GARCH models to each series in our data frame of returns.
# Let's select a VAR order of zero for the mean equation, and simply use the mean of each series.
# The mean equation is thus in our case simply: Rtn_t = mean(R) + et

# Then, for every series, a standard univariate GARCH(1,1) is run - giving us:
# et and sigmat, which is then used to calculate the standardized resids, zt, which is used in DCC calcs after.
DCCPre <- dccPre(xts_rtn, include.mean = T, p = 0)

## Sample mean of the returns:  0.0001245249 0.0005285527 0.001052786 0.0003952611 0.000234691 0.0004003896 0.0002159652 0.000443293 
## Component:  1 
## Estimates:  1.1e-05 0.071618 0.898339 
## se.coef  :  2e-06 0.009621 0.014006 
## t-value  :  4.714197 7.443803 64.13791 
## Component:  2 
## Estimates:  5e-06 0.068833 0.914122 
## se.coef  :  1e-06 0.008174 0.010532 
## t-value  :  4.025368 8.421161 86.79259 
## Component:  3 
## Estimates:  6e-06 0.057782 0.924458 
## se.coef  :  2e-06 0.008299 0.011618 
## t-value  :  3.763205 6.962738 79.5718 
## Component:  4 
## Estimates:  7e-06 0.058119 0.923202 
## se.coef  :  2e-06 0.007903 0.0106 
## t-value  :  4.246679 7.353888 87.0911 
## Component:  5 
## Estimates:  8e-06 0.060769 0.915368 
## se.coef  :  2e-06 0.008881 0.012882 
## t-value  :  4.062661 6.8428 71.05764 
## Component:  6 
## Estimates:  8e-06 0.07214 0.906428 
## se.coef  :  2e-06 0.008714 0.010868 
## t-value  :  4.733472 8.278996 83.40404 
## Component:  7 
## Estimates:  6e-06 0.067161 0.917309 
## se.coef  :  1e-06 0.008357 0.010255 
## t-value  :  4.176252 8.036888 89.44735 
## Component:  8 
## Estimates:  6e-06 0.078714 0.903836 
## se.coef  :  1e-06 0.010034 0.012244 
## t-value  :  4.162209 7.844578 73.81934

# If you want to fit other univariate garch models for each series, use the fGarch package to do so.

names(DCCPre)

## [1] "marVol"  "sresi"   "est"     "se.coef"

# We now have the estimates of volatility for each series. 
# Follow my lead below in changing the output to a usable Xts series for each column in xts_rtn:
Vol <- DCCPre$marVol
colnames(Vol) <- colnames(xts_rtn)
Vol <- 
  data.frame( cbind( date = index(xts_rtn), Vol)) %>% # Add date column which dropped away...
  mutate(date = as.Date(date)) %>%  tbl_df()  # make date column a date column...
TidyVol <- Vol %>% gather(Stocks, Sigma, -date)
ggplot(TidyVol) + geom_line(aes(x = date, y = Sigma, colour = Stocks))

# From the figure, it is clear that RMH is the most volatile stock...

# ------- Back to DCC:

# After saving now the standardized residuals:
StdRes <- DCCPre$sresi
# We can now use these sresids to calculate the DCC model.

# BUT FIRST NOTE THIS:
# Now, here follows a CLASSIC example of bad names for package commands. 
# If you run the dccFit command, notice that it gives you an error for filter_ 

# Upon investigation (Cntrl + click on the word dccFit in Rstudio) - you will notice
# the dccFit function uses the command 'filter'. The problem is, dplyr also uses filter.
# The MTS authors should have wrapped the command as stats::filter and not filter, as it 
# produces ambiguity between stats::filter and dplyr::filter...
# Try it yourself - Run the following code:
# DCC <- dccFit(StdRes, type="Engle")

# SO... to solve this petty issue, let's detach the tidyr and dplyr packages, 
# then run dccFit and then reload tidyr and dplyr...  
# (Note this takes a few minutes - go get coffee in the meantime):

detach("package:tidyverse", unload=TRUE)
detach("package:tbl2xts", unload=TRUE)
DCC <- dccFit(StdRes, type="Engle")

## Estimates:  0.95 0.01556144 6.902376 
## st.errors:  0.01451932 0.002667741 0.2442795 
## t-values:   65.43008 5.833189 28.25607

pacman::p_load("tidyverse", "tbl2xts", "broom")

Finally, we have our DCC model estimated…

So what can we do with it? Note that it is again a S4 class (i.e., not a normal data.frame) - so we need to be guided by the authors what it is that we can gain from it…

From the equations above, we know that we can now calculate the bivariate correlation between all the pairs in our data set. As we have 7 series, this translates to 42 pairs…

Let’s plot all the bivariate time-varying correlations with RMH from our DCC model:

Rhot <- DCC$rho.t
# Right, so it gives us all the columns together in the form:
# X1,X1 ; X1,X2 ; X1,X3 ; ....

# So, let's be clever about defining more informative col names. 
# I will create a renaming function below:
ReturnSeries = xts_rtn
DCC.TV.Cor = Rhot

renamingdcc <- function(ReturnSeries, DCC.TV.Cor) {
  
ncolrtn <- ncol(ReturnSeries)
namesrtn <- colnames(ReturnSeries)
paste(namesrtn, collapse = "_")

nam <- c()
xx <- mapply(rep, times = ncolrtn:1, x = namesrtn)
# Now let's be creative in designing a nested for loop to save the names corresponding to the columns of interest.. 

# TIP: draw what you want to achieve on a paper first. Then apply code.

# See if you can do this on your own first.. Then check vs my solution:

nam <- c()
for (j in 1:(ncolrtn)) {
for (i in 1:(ncolrtn)) {
  nam[(i + (j-1)*(ncolrtn))] <- paste(xx[[j]][1], xx[[i]][1], sep="_")
}
}

colnames(DCC.TV.Cor) <- nam

# So to plot all the time-varying correlations wrt SBK:
 # First append the date column that has (again) been removed...
DCC.TV.Cor <- 
    data.frame( cbind( date = index(ReturnSeries), DCC.TV.Cor)) %>% # Add date column which dropped away...
    mutate(date = as.Date(date)) %>%  tbl_df() 

DCC.TV.Cor <- DCC.TV.Cor %>% gather(Pairs, Rho, -date)

DCC.TV.Cor

}

# Let's see if our function works! Excitement!
Rhot <- 
  renamingdcc(ReturnSeries = xts_rtn, DCC.TV.Cor = Rhot)

head(Rhot %>% arrange(date))

## # A tibble: 6 x 3
##   date       Pairs     Rho
##   <date>     <chr>   <dbl>
## 1 2006-01-03 ABG_ABG 1    
## 2 2006-01-03 ABG_BVT 0.473
## 3 2006-01-03 ABG_CPI 0.302
## 4 2006-01-03 ABG_FSR 0.678
## 5 2006-01-03 ABG_NED 0.666
## 6 2006-01-03 ABG_RMH 0.652

# Let's now create a plot for all the stocks relative to the other stocks...
g1 <- 
  ggplot(Rhot %>% filter(grepl("ABG_", Pairs ), !grepl("_ABG", Pairs)) ) + 
  geom_line(aes(x = date, y = Rho, colour = Pairs)) + 
  theme_hc() +
  ggtitle("Dynamic Conditional Correlations: ABG")

print(g1)

g2 <- 
  g1 %+% subset(Rhot %>% filter(grepl("CPI_", Pairs ), !grepl("_CPI", Pairs))) + ggtitle("Dynamic Conditional Correlations: CPI")
print(g2)

# So... try write a function that loops through all your unique ticker names and plots the above..

… And those are the time-varying conditional correlations between all the stock pairs in our data set. And after calculating all these series’ time-varying correlation estimates using DCC-specifications, you rub your fur clean first before carrying on with your analysis:

You can now use these measures of time-varying correlations to study the level changes during different periods - particularly assessing whether and by how much correlation changed.

What makes these estimates different from rolling correlations?

For the same token as univariate garch model estimates of provide us with a “noise-reduce” version of volatility (in that the signal, , has been separated from noise, , using the volatility model ), so a DCC estimate of time-varying correlation provides us with a cleaner, noise-reduced version of the bivariate comovement structure in considered series.

DCC: flexible Univariate specs

You can ofcourse be flexible ito your specifications (but be wary of complicating steps for the pure sake of it)…

# Using the rugarch package, let's specify our own univariate
# functions to be used in the dcc process:

# Step 1: Give the specficiations to be used first:

# A) Univariate GARCH specifications:
uspec <- ugarchspec(variance.model = list(model = "gjrGARCH", 
    garchOrder = c(1, 1)), mean.model = list(armaOrder = c(1, 
    0), include.mean = TRUE), distribution.model = "sstd")
# B) Repeat uspec n times. This specification should be
# self-explanatory...
multi_univ_garch_spec <- multispec(replicate(ncol(xts_rtn), uspec))

# Right, so now every series will have a GJR Garch univariate
# specification. (see ?ugarchspec for other options...)

# C) DCC Specs
spec.dcc = dccspec(multi_univ_garch_spec, dccOrder = c(1, 1), 
    distribution = "mvnorm", lag.criterion = c("AIC", "HQ", "SC", 
        "FPE")[1], model = c("DCC", "aDCC")[1])  # Change to aDCC e.g.

# D) Enable clustering for speed:
cl = makePSOCKcluster(10)

# ------------------------ Step 2: The specs are now saved.
# Let's now build our DCC models...  ------------------------

# First, fit the univariate series for each column:
multf = multifit(multi_univ_garch_spec, xts_rtn, cluster = cl)

# Now we can use multf to estimate the dcc model using our
# dcc.spec:
fit.dcc = dccfit(spec.dcc, data = xts_rtn, solver = "solnp", 
    cluster = cl, fit.control = list(eval.se = FALSE), fit = multf)

# And that is our DCC fitted model!

# We can now test the model's fit as follows: Let's use the
# covariance matrices to test the adequacy of MV model in
# fitting mean residual processes:
RcovList <- rcov(fit.dcc)  # This is now a list of the monthly covariances of our DCC model series.
covmat = matrix(RcovList, nrow(xts_rtn), ncol(xts_rtn) * ncol(xts_rtn), 
    byrow = TRUE)
mc1 = MCHdiag(xts_rtn, covmat)

## Test results:  
## Q(m) of et: 
## Test and p-value:  65.60914 3.098944e-10 
## Rank-based test: 
## Test and p-value:  169.6296 0 
## Qk(m) of epsilon_t: 
## Test and p-value:  700.1147 0.0495894 
## Robust Qk(m):  
## Test and p-value:  921.517 1.69087e-12

# ....Check Tsay (2014 on the interpretation of these
# Portmanteau tests)....

# Now to save the time-varying correlations as specified by
# the DCC model, it again requires some gymnastics from our
# side.  First consider what the list looks like:
dcc.time.var.cor <- rcor(fit.dcc)
print(dcc.time.var.cor[, , 1:3])

## , , 2006-01-03
## 
##           ABG       BVT       CPI       FSR       NED       RMH       SBK
## ABG 1.0000000 0.4762022 0.3099609 0.6787156 0.6675293 0.6555642 0.6956543
## BVT 0.4762022 1.0000000 0.2733245 0.5141451 0.4887010 0.5252806 0.5097593
## CPI 0.3099609 0.2733245 1.0000000 0.3363768 0.3171730 0.3290350 0.3119902
## FSR 0.6787156 0.5141451 0.3363768 1.0000000 0.6794219 0.8454864 0.7201598
## NED 0.6675293 0.4887010 0.3171730 0.6794219 1.0000000 0.6559834 0.6958803
## RMH 0.6555642 0.5252806 0.3290350 0.8454864 0.6559834 1.0000000 0.7012190
## SBK 0.6956543 0.5097593 0.3119902 0.7201598 0.6958803 0.7012190 1.0000000
## SLM 0.5324867 0.4947969 0.2930279 0.5699836 0.5401437 0.5872740 0.5556683
##           SLM
## ABG 0.5324867
## BVT 0.4947969
## CPI 0.2930279
## FSR 0.5699836
## NED 0.5401437
## RMH 0.5872740
## SBK 0.5556683
## SLM 1.0000000
## 
## , , 2006-01-04
## 
##           ABG       BVT       CPI       FSR       NED       RMH       SBK
## ABG 1.0000000 0.4742107 0.2920909 0.6849396 0.6517577 0.6571002 0.6987774
## BVT 0.4742107 1.0000000 0.2699428 0.5111949 0.4864373 0.5251849 0.5091313
## CPI 0.2920909 0.2699428 1.0000000 0.3164803 0.3210558 0.3198876 0.2997783
## FSR 0.6849396 0.5111949 0.3164803 1.0000000 0.6617337 0.8443044 0.7228568
## NED 0.6517577 0.4864373 0.3210558 0.6617337 1.0000000 0.6490365 0.6857606
## RMH 0.6571002 0.5251849 0.3198876 0.8443044 0.6490365 1.0000000 0.7027941
## SBK 0.6987774 0.5091313 0.2997783 0.7228568 0.6857606 0.7027941 1.0000000
## SLM 0.5219732 0.4935449 0.2950783 0.5577587 0.5410254 0.5828711 0.5492294
##           SLM
## ABG 0.5219732
## BVT 0.4935449
## CPI 0.2950783
## FSR 0.5577587
## NED 0.5410254
## RMH 0.5828711
## SBK 0.5492294
## SLM 1.0000000
## 
## , , 2006-01-05
## 
##           ABG       BVT       CPI       FSR       NED       RMH       SBK
## ABG 1.0000000 0.4732274 0.2920766 0.6832127 0.6483952 0.6537574 0.6985511
## BVT 0.4732274 1.0000000 0.2680315 0.5038797 0.4912005 0.5140473 0.5099090
## CPI 0.2920766 0.2680315 1.0000000 0.3166800 0.3174179 0.3195627 0.2992519
## FSR 0.6832127 0.5038797 0.3166800 1.0000000 0.6496762 0.8446447 0.7192487
## NED 0.6483952 0.4912005 0.3174179 0.6496762 1.0000000 0.6319254 0.6847797
## RMH 0.6537574 0.5140473 0.3195627 0.8446447 0.6319254 1.0000000 0.6966282
## SBK 0.6985511 0.5099090 0.2992519 0.7192487 0.6847797 0.6966282 1.0000000
## SLM 0.5219615 0.4913328 0.2949618 0.5572370 0.5365909 0.5810885 0.5487272
##           SLM
## ABG 0.5219615
## BVT 0.4913328
## CPI 0.2949618
## FSR 0.5572370
## NED 0.5365909
## RMH 0.5810885
## SBK 0.5487272
## SLM 1.0000000

# Okay, so this format presents a particular challenge in
# getting our time-varying series into a plottable format...
# Every date is a list, and the list a matrix... How to get
# it in a tidy and plottable format?!!

# If you're up for the challenge - see if you can solve this
# little conundrum of getting the lists as bivariate pair
# columns before using my solution.

# -------------SOLUTION:--------------------

# We will first use the base R function aperm - which
# transposes any array into a list by changing the dimensions
dcc.time.var.cor <- aperm(dcc.time.var.cor, c(3, 2, 1))
dim(dcc.time.var.cor) <- c(nrow(dcc.time.var.cor), ncol(dcc.time.var.cor)^2)

# And now we can rename our columns the same way as before.
# Luckily we wrote a function so we can use it again...

dcc.time.var.cor <- renamingdcc(ReturnSeries = xts_rtn, DCC.TV.Cor = dcc.time.var.cor)

# Note that the figure is very similar to our earlier DCC TV
# correlations.  So having a GJR univariate model didn't
# exactly change much...

g1 <- ggplot(dcc.time.var.cor %>% filter(grepl("SBK_", Pairs), 
    !grepl("_SBK", Pairs))) + geom_line(aes(x = date, y = Rho, 
    colour = Pairs)) + theme_hc() + ggtitle("Dynamic Conditional Correlations: SBK")

print(g1)

Go-GARCH: Orthogonalizing sources of volatility

Another different approach to estimating second order persistence models is based on the assumption that returns are generated by a set of unobserved orthogonal conditionally heteroskedastic factors (c.f. Van der Weide (2002)). See my presentation on this under the Research Tab, or click here.

These orthogonal components (think PCA on second moment…) are measured by identifying independent and uncorrelated factors that make up the var-covar matrix .

The statistical transformations are done as follows:

with linking the unobserved uncorrelated components with the observed residual process. : constant and invertible ; : normalized, unit variance.

Also: represents the unobserved independent factors assigned to each series (factor weights), such that:

with

and

as before and:

So that the conditional covariance matrix is given by

Intuitively…

What is happening is that we are using the residual series of a simple mean return series to decompose the var-covar matrix into a factor loading matrix, A, which is made up of:

a sample (unconditional) covariance estimate ;
an orthogonal matrix, , mapping the uncorrelated factors onto the covar matrix: estimated using Independent Component Analysis (ICA).
Thus it assumes the variance process is driven by independent orthogonal latent components.
The aim is then to model the independent factors as independent univariate GARCH processes.

I strongly suggest reading wake me up before you Go-GARCH, which offers a very simple intuitive overview of this technique, as well as a neat application section in section 6 and 7 - a very good source for project ideas…

An example of a Normally distributed GoGARCH(1,1) model is:

With each uncorrelated component driven by a GARCH(1,1) process:

with : conditional variances of the factors, , so that:

The conditional covariances of are then driven by:

Why should GO-GARCH be used?

The DCC model estimations allow for the direct modelling of the conditional correlation.
Although efficient, this method has some drawbacks, as e.g. it assumes a constant structure to the correlation dynamics.
The GoGARCH model uses a highly efficient and less parameter intensive estimation technique to decompose the var-covar matrix into orthogonal sources of volatility, which could then be calculated using an additive log-likelihood approach because of the independence of volatility factors.
Although computationally highly efficient, the assumption of linearly indendent volatility factors is strong.

Fitting the Go-GARCH

We will once again use the rmgarch package to fit the Go-GARCH model. Also see the vignette for technical details (specifically, see section 2.4 for Go-GARCH).

Tip: search for the rmgarch folder installed on your computer - go to rmgarch/doc, and opening The_rmgarch_models.Rnw you will find the source code for the vignette - and by definition all the math equations needed for these models!

So let’s fit these models on our data:

# Let's keep using our earlier first step individual garch definitions (GJR):
# Now we create a gogarch model specification:
spec.go <- gogarchspec(multi_univ_garch_spec, 
                       distribution.model = 'mvnorm', # or manig.
                       ica = 'fastica') # Note: we use the fastICA
cl <- makePSOCKcluster(10)
multf <- multifit(multi_univ_garch_spec, xts_rtn, cluster = cl)

fit.gogarch <- gogarchfit(spec.go, 
                      data = xts_rtn, 
                      solver = 'hybrid', 
                      cluster = cl, 
                      gfun = 'tanh', 
                      maxiter1 = 40000, 
                      epsilon = 1e-08, 
                      rseed = 100)

print(fit.gogarch)

## 
## *------------------------------*
## *        GO-GARCH Fit          *
## *------------------------------*
## 
## Mean Model       : CONSTANT
## GARCH Model      : sGARCH
## Distribution : mvnorm
## ICA Method       : fastica
## No. Factors      : 8
## No. Periods      : 3568
## Log-Likelihood   : 84138.58
## ------------------------------------
## 
## U (rotation matrix) : 
## 
##         [,1]      [,2]    [,3]    [,4]    [,5]    [,6]     [,7]    [,8]
## [1,] -0.3544 -0.421955 -0.0437 -0.4208  0.3902  0.2396  0.13250 -0.5387
## [2,]  0.0598  0.126882 -0.0105  0.1081 -0.1181  0.0548  0.97451 -0.0438
## [3,] -0.2857 -0.246671  0.0524 -0.3003  0.3033 -0.3885  0.17361  0.7011
## [4,] -0.7770 -0.018505 -0.3467  0.3938 -0.2931  0.1564 -0.03742  0.0942
## [5,] -0.1208  0.015008 -0.0763  0.0565 -0.1254 -0.8725  0.01218 -0.4462
## [6,] -0.0672  0.032895  0.0274 -0.6782 -0.7278  0.0490 -0.01395  0.0372
## [7,]  0.2759 -0.862183  0.0122  0.2744 -0.3188 -0.0148  0.02938  0.0482
## [8,]  0.3005 -0.000178 -0.9318 -0.1659  0.0938 -0.0193  0.00545  0.0684
## 
## A (mixing matrix) : 
## 
##          [,1]    [,2]      [,3]    [,4]     [,5]      [,6]      [,7]    [,8]
## [1,] 0.003064 0.00488 -0.001036 0.01592 -0.00521 -0.003006 -1.41e-03 0.00599
## [2,] 0.000285 0.00208  0.000465 0.00279 -0.00109  0.000503 -6.05e-05 0.01673
## [3,] 0.002238 0.00162  0.000687 0.00168 -0.00123 -0.001985 -1.75e-02 0.00405
## [4,] 0.013398 0.00725 -0.003129 0.00597 -0.00631 -0.002211 -1.05e-03 0.00819
## [5,] 0.001911 0.00436 -0.000532 0.00511 -0.01493 -0.003074 -1.51e-03 0.00712
## [6,] 0.011495 0.00780  0.008026 0.00631 -0.00615 -0.002067 -4.57e-04 0.00809
## [7,] 0.001942 0.01480 -0.000631 0.00609 -0.00567 -0.002658 -1.20e-03 0.00672
## [8,] 0.003568 0.00379  0.001047 0.00310 -0.00269 -0.014535 -3.51e-04 0.00862

# Extracting time-varying conditional correlations: You know the drill...
gog.time.var.cor <- rcor(fit.gogarch)
gog.time.var.cor <- aperm(gog.time.var.cor,c(3,2,1))
dim(gog.time.var.cor) <- c(nrow(gog.time.var.cor), ncol(gog.time.var.cor)^2)
# Finally:
gog.time.var.cor <-
renamingdcc(ReturnSeries = xts_rtn, DCC.TV.Cor = gog.time.var.cor)

Now, see the discussion on p.21 of Wake me up before you gogarch, for interpreting the model parameter outputs of the Go-GARCH, as well as the discussion of tv-correlation estimates between DCC and Go-GARCH.

g1 <- ggplot(dcc.time.var.cor %>% filter(grepl("SBK_", Pairs), 
    !grepl("_SBK", Pairs))) + geom_line(aes(x = date, y = Rho, 
    colour = Pairs)) + theme_hc() + ggtitle("DCC: SBK")

g2 <- ggplot(gog.time.var.cor %>% filter(grepl("SBK_", Pairs), 
    !grepl("_SBK", Pairs))) + geom_line(aes(x = date, y = Rho, 
    colour = Pairs)) + theme_hc() + ggtitle("Go-GARCH: SBK")

# DCC tv corr:
print(g1)

# Go-GARCH tv corr:
print(g2)

Time-varying third and fourth moments

In addition to providing us with an estimate of the time-varying second moments, the ACD model allows a closed form solution of the time-varying skewness and kurtosis measures. I would be very interested in someone studying the stability of the third and fourth moment in the domestic market. See here for a step-by-step example of setting up such a tv third and fourth moment model. Important: install the needed packages (racd package) here.

Other interesting models

The literature on second order moment modelling is fast expanding, and there are very many brilliant techniques that can be applied and used.

I am quite interested in developments in:

GAS models
Bayesian GARCH models

I recommend checking out the GAS model page for interesting papers and recommendations on studies. The recent gasmodel R package also provides easily accessible code for fitting these models.. so go play!

References

Bauwens, Luc, Sébastien Laurent, and Jeroen VK Rombouts. 2006. “Multivariate Garch Models: A Survey.” Journal of Applied Econometrics 21 (1): 79–109.

Christiansen, Charlotte. 2007. “Volatility-Spillover Effects in European Bond Markets.” European Financial Management 13 (5): 923–48.

Engle, Robert. 2002. “Dynamic Conditional Correlation: A Simple Class of Multivariate Generalized Autoregressive Conditional Heteroskedasticity Models.” Journal of Business & Economic Statistics 20 (3): 339–50.

Ghalanos, Alexios. 2014. rmgarch: Multivariate Garch Models.

Horvath, Roman, and Petr Poldauf. 2012. “International Stock Market Comovements: What Happened During the Financial Crisis?” Global Economy Journal 12 (1).

Johansson, Anders C. 2008. “Interdependencies Among Asian Bond Markets.” Journal of Asian Economics 19 (2): 101–16.

Silvennoinen, Annastiina, and Timo Teräsvirta. 2009. “Multivariate Garch Models.” In Handbook of Financial Time Series, 201–29. Springer.

Tsay, Ruey S. 2014. An Introduction to Analysis of Financial Data with R. John Wiley & Sons.

Van der Weide, Roy. 2002. “GO-Garch: A Multivariate Generalized Orthogonal Garch Model.” Journal of Applied Econometrics 17 (5): 549–64.

Financial Econometrics Practical

Practical 7: Multi-variate Volatility Modelling

NF Katzke