4.2 MGARCH modelling in R

TABLE 4.1: Commands for DCC – multivariate GARCH modelling within package `rmgarch`
Command	Description
`DCCtest()`	testing for non-constant correlations
`dccspec()`	creating a DCC specification object prior to fitting
`dccfit()`	estimating (fitting) specified DCC model
`show()`	displays the results of estimated DCC model (mGARCHfit object)
`fitted()`	produces fitted conditional means
`sigma()`	produces fitted conditional standard deviations
`residuals()`	extracts fitted conditional mean residuals
`dccforecast()`	forecasts the future variances and covariances for h-steps ahead
`infocriteria()`	exctrcts information criteria, such as AIC, BIC,…
`coef()`	displays coefficient vector
`likelihood()`	extracts the joint maximum likelihhod
`plot()`	displys $4$ possible plots by selecting numbers from $1$ do $4$
`nisurface()`	plots news impact surface
`rcor()`	extracts fitted dynamic conditional correlations
`rcov()`	extracts fitted dynamic conditional covariances
`rshape()`	obtains degrees of freedom from multivariate t–distribution
`rskew()`	obtains skewnees of multivariate t–distribution

Example 14. Estimate Engle’s DCC(

$1,1$ ) model to study the time–varying beta coefficient between Tesla’s returns and the S&P500 index. Install and load rmgarch package in R for the first time. Download daily prices for both assets, covering the period from “2021-01-01” to “2024-12-31”, and compute objects tesla.returns and sp500.returns (ensuring that any missing values are removed). Specify univariate GARCH(

$1,1$ ) models for both Tesla and S&P500 returns, which is necessary before specifying the multivariate DCC–GARCH model (assume univariate as well as multivariate normal distributions). Once the DCC(

$1,1$ )–GARCH model is fitted, extract and plot the dynamic conditional correlations. Test the joint null hypothesis

$H_0:~a_1=b_1=0$ using DCCtest() command. Finally, calculate the time–varying beta coefficients by multiplying the dynamic conditional correlation by the ratio of the Tesla and the S&P500 variances. After plotting time–varying beta’s, identify the dates when the beta coefficient was less than or equal to

$1$ , and summarize the beta coefficients using a five–number summary.

Solution

Copy the code lines below to the clipboard, paste them into an R Script file, and run them.
In this example estimate parameters

$\mu$ ,

$\omega$ ,

$\alpha_1$ and

$\gamma_1$ are interpreted in the same way as in univariate GARCH’s, whereas

$a_1$ and

$b_1$ are not interpretable but can serve to test the null hypothesis about constant correlations, which is not rejected in this example (

$p-$ value

$=0.40094$ ). This indicates that CCC model is more appropriate despite to significance of parameter

$b_1$ and non-significance of parameter

$a_1$ . However, plot exhibits time–varying correlations which are all positive in the range from

$0.21487$ to

$0.70336$ . Beta coefficients also varied from

$0.86669$ to

$3.38423$ , while in the first week of October

$2022$ recorded values below

$1$ .

# Installing and loading new package "rmgarch"
install.packages("rmgarch")
library(rmgarch) # also requires "rugarch" for univariate GARCH

# Downloading the data for Tesla and S&P500 index
getSymbols(c("TSLA","^GSPC"), # tickers for both equities
           src="yahoo", # specifying the source of the data
           from = as.Date("2021-01-01"), # starting date in the format "yyyy-mm-dd"
           to = as.Date("2024-12-31")) # ending date in the format "yyyy-mm-dd"

# Compute Tesla and S&P500 daily log returns and remove missig values (if any)
tesla.returns <- na.omit(dailyReturn(Cl(TSLA), type = "log"))
sp500.returns <- na.omit(dailyReturn(Cl(GSPC), type = "log"))

# Merge returns into a single data frame with renamed columns
# Keep only the dates that exist in both returns
returns_data <- merge(tesla.returns, sp500.returns, join = "inner") 
colnames(returns_data) <- c("Tesla", "SP500")

# Univariate GARCH(1,1) should be specified in the first step prior to DCC specification
univariate <- ugarchspec(
  variance.model = list(model = "sGARCH", garchOrder = c(1, 1)), # standard GARCH(1,1)
  mean.model = list(armaOrder = c(0, 0), include.mean = TRUE),    # constant mean
  distribution.model = "norm")                                    # normal distribution

# Multivariate DCC-GARCH specification in the second step
multivariate <- dccspec(
  uspec = multispec(replicate(2, univariate)), # same univariate GARCH for both assets
  dccOrder = c(1, 1),              # DCC(1,1) order
  distribution = "mvnorm")          # multivariate normal innovations

# Fitting the DCC(1,1)-GARCH model
dcc_fit <- dccfit(multivariate, data = as.matrix(returns_data))
print(dcc_fit) # display estimation results

# Performing DCC test to check for constant or dynamic correlation
DCCtest(as.matrix(returns_data))

# Extracting dynamic correlations
dcc_rho <- rcor(dcc_fit)            # 1004 correlation matrices 2x2
rho_t <- dcc_rho[1, 2, ]            # extracting off-diagonal elements 
rho_t <- xts(rho_t, order.by = index(returns_data))
colnames(rho_t) <- c("rho")
rho_t[1:10]                        # preview first 10 correlations

# Plotting dynamic correlations
library(ggplot2)
ggplot(rho_t, aes(x = Index, y = rho)) +
  geom_line(color = "orchid") + 
  labs(title = "Dynamic correlations between Tesla and S&P500", x = "Days", y = "Correlation") +
  theme_minimal()

# Calculating time-varying beta: beta_t = rho_t * (sigma_TSLA / sigma_sp500)
sigma <- sigma(dcc_fit)
index(sigma) <- as.Date(index(sigma))
sigma_tesla <- sigma[, "Tesla"]    # Tesla volatility
sigma_sp500 <- sigma[, "SP500"]  # S&P500 volatility
beta_t <- rho_t * (sigma_tesla / sigma_sp500)

# Plotting time-varying betas
colnames(beta_t) <- c("beta")
ggplot(beta_t, aes(x = Index, y = beta)) +
  geom_line(color = "orchid") + 
  labs(title = "Time-varying beta of Tesla", x = "Days", y = "Beta coefficient") +
  theme_minimal()

# Finding the dates when beta was less or equal to 1
beta_t[beta_t <= 1]

# Summarizing beta coefficients
quantile(beta_t) #  five-number summary

$~~~$

Example 15. Estimate both a standard (symmetric) and an asymmetric DCC(

$1,1$ ) model to analyze the dynamic correlation between BitCoin and the S&P500 index. Specify a univariate GJR–GARCH(

$1,1$ ) models with t–distribution of innovations (with fixed degrees of freedom

$5$ ) and use it to specify DCC(

$1,1$ ) models as objects multivariate1 and multivariate2. Estimate two models as dcc_fit1 and dcc_fit2 objects. Extract and plot the time–varying correlations between BitCoin and the S&P500 index based on dcc_fit1 and identify periods where the correlation was zero or negative. Finally, summarize and compare both models using custom–defined tidy.DCCfit and glance.DCCfit functions, and present the results using the modelsummary() command. Plot news impact surface from object dcc_fit2.

Solution

Copy the code lines below to the clipboard, paste them into an R Script file, and run them.

# Downloading the data for BitCoin and S&P500 index starting from January 2020
getSymbols(c("^GSPC","BTC-USD"), # tickers for equity and cryptocurrency
           src="yahoo", # specifying the source of the data
           from = as.Date("2020-01-01"), # also covers the pandemic period
           to = as.Date("2024-12-31"))

# Compute S&P500 and BitCoin daily log returns and removing missing values (if any)
sp500_returns <- na.omit(dailyReturn(Cl(GSPC), type = "log"))
bitcoin_returns <- na.omit(dailyReturn(Cl(`BTC-USD`), type = "log")) # use backticks `` when there is dash (-) in objects name


# Merge returns into a single data frame with renamed columns
# Keep only the dates that exist in both returns
returns_data <- merge(sp500_returns, bitcoin_returns, join = "inner") 
colnames(returns_data) <- c("SP500","BitCoin")

# Estimating one correlation for entire period
cor(returns_data)

# Univariate GARCH(1,1) should be specified in the first step prior do DCC specification
univariate <- ugarchspec(
  variance.model = list(model = "sGARCH", garchOrder = c(1, 1)), # standard GARCH(1,1)
  mean.model = list(armaOrder = c(0, 0), include.mean = TRUE),    # constant mean
  distribution.model = "std")                                    # normal distribution of innovations

# Multivariate DCC-GARCH specification in the second step
multivariate1 <- dccspec(
  uspec = multispec(replicate(2, univariate)), # same univariate GARCH for both assets
  dccOrder = c(1, 1),              # DCC(1,1) order
  model="DCC", # standard DCC model 
  distribution = "mvt")          # multivariate normal innovations

# Fitting the DCC(1,1)-GARCH model
dcc_fit1 <- dccfit(multivariate1, data = as.matrix(returns_data))
print(dcc_fit1) # Display estimation results

# Extracting dynamic correlations
dcc_rho <- rcor(dcc_fit1)            # 1257 correlation matrices 2x2
rho_t <- dcc_rho[1, 2, ]            # extracting off-diagonal elements 
rho_t <- xts(rho_t, order.by = index(returns_data))
colnames(rho_t) <- c("rho")
rho_t[1:10]                        # preview first 10 correlations

# Plotting dynamic correlations
ggplot(rho_t, aes(x = Index, y = rho)) +
  geom_line(color = "orchid") + 
  labs(title = "Dynamic correlations between S&P500 and BitCoin", x = "Days", y = "Correlation") +
  theme_minimal()

# Alternative to dynamic correlation plot
plot(dcc_fit1,which=4)

# Covariance can be also plotted
plot(dcc_fit1,which=3)

# Finding the dates when correlation was less or equal to 0
rho_t[rho_t <= 0]

# Two functions are defined that return two data frames: tidy and glance
# This is required as modelsummary() command does not support MGARCH model
tidy.DCCfit<- function(x, ...) {
  s <- x@mfit$matcoef
  ret <- data.frame(
    term      = row.names(s),
    estimate  = s[, " Estimate"],
    std.error = s[, " Std. Error"],
    statistic = s[, " t value"],
    p.value =  s[,"Pr(>|t|)"])
  ret
}

glance.DCCfit <- function(x, ...) {
  s <- x@mfit
  ret <- data.frame(
    Likelihood = round(s$llh, digits=2),
    Akaike=round(infocriteria(x)[1], digits=4),
    Bayes=round(infocriteria(x)[2], digits=4))
  ret
}

# Fitting an asymmetric DCC(1,1) model by incorporating univariate GJR-GARCH models
univariate2 <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)), # asymmetric GARCH(1,1)
  mean.model = list(armaOrder = c(0, 0), include.mean = TRUE),    # constant mean
  distribution.model = "std",                                    # t- distribution of innovations
  fixed.pars=list(shape=5))

# Same univariate GARCH for both assets
multivariate2 <- dccspec(uspec = multispec(replicate(2, univariate2)),
                         dccOrder = c(1, 1),              # DCC(1,1) order
                         model="DCC", # asymmetric DCC model 
                         distribution = "mvt")          # multivariate t-distribution

dcc_fit2 <- dccfit(multivariate2, data = as.matrix(returns_data))

modelsummary(list("DCC(1,1) model"=dcc_fit1,"Asymmetric DCC(1,1)"=dcc_fit2),
             stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01), # customizing significance levels
             fmt=4)

# Plotting news impact surface from object dcc_fit2
nisurface(dcc_fit2)

$~~~$