Seminar 5
In this seminar, we will use Y.RData file that we created in the Seminar 2. We focus on coding multivariate volatility models using rmgarch package. Please review the theory of bivariate EWMA and DCC models that we discussed in the lecture.
The plan for this week:
1. Introduce multivariate volatility models
2. Build a bivariate EWMA model
3. Run DCC models with different specifications
4. Compare models
The following sections and subsections are materials for students to learn by themselves: DCC apARCH and tapARCH and DCC with several assets.
library(rmgarch)
library(lubridate)
Multivariate volatility models¶
In a univariate volatility model, we consider a stock with returns $Y_t$ that can be written as:
$$Y_t = \sigma_{t} Z_{t}$$Where $\sigma_t$ is the conditional volatility and $Z_t$ are random shocks.
For many financial applications, we have to consider a vector of assets, instead of a single asset. So, if we have $K>1$ assets, it is necessary to indeicate which asset and paremeters are being referred to. The notation we will use is:
$$Y_{t,i} = \sigma_{t,i} Z_{t,i}$$Where the first subscript indicates the date, and the second subscript the asset.
Now our conditional covariance is a symmetric square conditional covariance matrix, where the dimension is the number of assets. For example, for three assets it would be:
$$ \begin{equation*} \Sigma_{t} = \begin{pmatrix} \sigma_{t,11} & & \\ \sigma_{t,12} & \sigma_{t,22} & \\ \sigma_{t,13} & \sigma_{t,23} & \sigma_{t,33} \end{pmatrix} \end{equation*} $$If we have a portfolio wight a vector $w$ of portfolio weights, then the portfolio variance would be:
$$\sigma^2_{\textrm{portfolio}} = w'\Sigma w$$As in the case of univariate models, we need to ensure that the variance is not negative. In this case, this means that the covariance matrix $\Sigma$ is positive semi-definite:
$$ |\Sigma| \geq 0$$Working with several assets presents the curse of dimensionality. This means that the number of variance and covariance terms will explode as the number of assets increase, which makes it challenging to estimate the covariance matrix and ensure its positive semi-definiteness.
For example, if we try to estimate a GARCH model for two assets, we would have 21 parameters to estimate, which is almost impossible in practice.
The concept of stationarity is more important in multivariate volatility models, since violating it could lead to numerical problems. In this seminar we will focus in implementing two different multivariate models:
- Exponentially-Weighted Moving Average
- Dynamic Conditional Correlation Models
To do so, we will work with the returns from JP Morgan and Citigroup.
# Load the data
load("Y.RData")
# Extract the returns for JPM and C
y <- Y[c("JPM", "C")]
EWMA¶
In the Exponentially-Weighted Moving Average model, the estimated conditional volatility at time $t$ is a convex combination of the estimation at time $t-1$ and the squared returns at time $t-1$. In general, if we have a vector of returns for each time $t$ that includes all assets $K$:
$$\textbf{y}_{t} = [y_{t,1}, y_{t,2}, \ldots, y_{t,K}]$$Then the multivariate EWMA is:
$$\hat{\Sigma}_{t} = \lambda \hat{\Sigma}_{t-1} + (1-\lambda) \textbf{y}'_{t-1} \textbf{y}_{t-1}$$The properties of this mdoel are:
- The same pre-specified weight, $\lambda$ is used for all assets
- The variance of any particular asset only depends on its own lags
Implementation¶
Since we will perform linear algebra operations to compute the EWMA, it is more convenient to work with a matrix instead of a dataframe, so let's first turn y from a dataframe into a matrix:
# Check the class of y
class(y)
# Transform into matrix
y <- as.matrix(y)
# Check class again
class(y)
Let's take a look at the returns:
# Plotting the returns in a 2x1 grid
par(mfrow = c(2,1))
plot(y[,1], type = "l", main = "Returns for JPM", col = 1)
plot(y[,2], type = "l", main = "Returns for C", col = 2)
# Reset the grid
par(mfrow = c(1,1))
We will create a EWMA matrix that will hold our estimations for the elements of $\Sigma$. First let's establish the number of entries the matrix should have:
# Determine number of entries
n <- dim(y)[1]
Now we create the matrix. We could initialize the matrix with a number, like 0, but it is recommended to initialize it with NA to avoid any potential mistakes.
To determine the number of columns in the matrix, we need to figure out how many parameters we are estimating in every time period. For any given time $t$, we will be estimating the conditional variance of each asset, and the conditional covariances. For $K$ assets, this is:
$$K + K(K-1)/2$$In our case, it is 3. Column 1 will hold the conditional variance of JPM, column 2 will hold the conditional covariance between JPM and C, and column 3 will hold the conditional variance for Citigroup:
# Initializing the EWMA matrix
EWMA <- matrix(NA, nrow = n, ncol = 3)
Why do we initialize a matrix/vector?¶
It is convenient to create the vector or matrix we want to populate before doing so for two reasons:
- Good programming practice. It is better to have a fixed-length object than to re-size an object in every iteration of the loop.
- Be able to spot mistakes. In case there is an error in your loop and try to allocate an element to a non-existing row, you will get an error message. If you are dynamically modifying the size of the object, you can add more rows/columns by accident without realizing it
We could fill the new matrix with any values, and it is common to fill them with zeros. However, it is useful to fill the matrix with NA. This way, if something in the loop doesn't work, our matrix will have some empty spots instead of holding zeros, and whenever we try to make any operation we will get an error message.
We determine the value for $\lambda$:
# Determine lambda
lambda <- 0.94
It is necessary is to determine how to estimate the conditional covariances of the first period. In other words, how do we estimate $\hat{\Sigma}_1$ without $\hat{\Sigma}_0$ or $\textbf{y}_0$. For this, we will use the unconditional covariance of the sample and "burn" the first few observations. The effect of a given conditional covariance from a past period quickly dies out as time passes, because $\lambda^n \rightarrow 0$ as $n \rightarrow \infty$, so the effect of initializing the EWMA matrix with the unconditional covariance will not be a problem after a few time periods. A rule of thumb is to burn the first 30 observations.
# Get the sample covariance
S <- cov(y)
S
We want to get the three unique values of the matrix. This can be done in three different ways shown below. We recommend using the last, since it is easily reproducible for any number of assets:
# Getting the unique values in three ways:
# 1. Creating a vector with the distinct elements
c(S[1,1], S[1,2], S[2,2])
# 2. Vectorizing the matrix and getting distinct elements
c(S)[c(1,2,4)]
# 3. Using the fact that S is symmetric and using upper.tri()/lower.tri()
# upper.tri() returns a matrix of logicals
# with entries TRUE in the upper triangle
S[upper.tri(S, diag = TRUE)]
# Fill the initial row of EWMA with the sample covariances
EWMA[1,] <- S[upper.tri(S, diag = TRUE)]
# View the matrix
head(EWMA)
The first column is the sample variance of JPM, the second column is the sample covariance, and the third column is the sample variance of C.
As an example, let's manually compute how we would estimate the variances and covariance for $t=2$:
# Manually computing EWMA elements for t = 2
# Apply the formula for EWMA
S_2 <- lambda * S + (1-lambda) * y[1,] %*% t(y[1,])
# Get the variances and covariances
S_2[upper.tri(S_2, diag = TRUE)]
Now that we have seen how to get the elements for a given time, we can write a for loop to populate the entire EWMA matrix:
# Populating the EWMA matrix
# Create a loop for rows 2 to n
for (i in 2:n) {
# Update S with the new weighted moving average
S <- lambda * S + (1-lambda) * y[i-1,] %*% t(y[i-1,])
# Fill the following EWMA row with the covariances
EWMA[i,] <- S[upper.tri(S, diag = TRUE)]
}
Let's see the head() of our EWMA matrix:
head(EWMA)
Important note on Matrix Operations¶
Since we are working with matrices, it is very important to pay attention to the dimensions of the elements, to make sure we don't do any mistake when multiplying vectors. Here is a short review on the order of vectors for multiplication:
# Matrix operations - The order is important
# This is a 2x1 vector
y[1,]
# This is a scalar (a 1x2 row vector multiples a 2x1 column vector)
t(y[1,]) %*% y[1,]
# This is a matrix (a 2x1 column vector multiples a 1x2 row vector)
y[1,] %*% t(y[1,])
Let's create a wrong_EWMA matrix to show what would happen if we were to mix up the transposes:
# Creating a wrong EWMA
# Initialize the matrix the same way
wrong_EWMA <- matrix(NA, nrow = n, ncol = 3)
S <- cov(y)
wrong_EWMA[1,] <- S[upper.tri(S, diag = TRUE)]
# Do the loop but interchange the transpose
for (i in 2:n) {
# Update S with the new weighted moving average
S <- lambda * S + (1-lambda) * t(y[i-1,]) %*% y[i-1,]
# Fill the following EWMA row with the covariances
wrong_EWMA[i,] <- S[upper.tri(S, diag = TRUE)]
}
Plotting the conditional variances and covariances¶
# Plotting estimated variances and covariances
matplot(EWMA, type = "l", main = "EWMA", lty = 1)
legend("topleft", legend = c("JPM", "Covar", "C"), col = 1:3, lty = 1)
We can calculate the correlation coefficient of the two stocks, which is the covariance over the square root of variances:
# Correlation coefficient
EWMArho <- EWMA[,2]/sqrt(EWMA[,1]*EWMA[,3])
# Plot
plot(EWMArho, type = "l", main = "Correlation coefficient of JPM and C", col = "red")
Let's replicate the plots including the date as the x-axis. We can get it from our Y data frame:
# Plots for conditional volatility
plot(Y$date, sqrt(EWMA[,1]), type = "l", main = "Conditional volatility JPM")
plot(Y$date, sqrt(EWMA[,3]), type = "l", main = "Conditional volatility C")
# Correlation Coefficient plot
plot(Y$date, EWMArho, type = "l", main = "Correlation coefficient of JPM and C", col = "red")
Dynamic Conditional Correlation models¶
To discuss and implement Dynamic Conditional Correlation (DCC) models, we first have to discuss Constant Conditional Correlation (CCC) models. In this type of model, we separate out correlation modeling from volatility model. For the former, we use a correlation matrix, and for the latter we use GARCH or some standard method.
We have discussed the theory of DCC models in lecture. We will implement the DCC model using the rmgarch package, which is the multivariate version of rugarch.
Implementation¶
Similarly to a univariate model from rugarch, we need to specify the model we are going to use. In this case, we use the dccspec() function. We will replicate a univariate specification for each asset and implement both in the DCC model:
# DCC model
# Specify the default univariate GARCH model with no mean
xspec <- ugarchspec(
mean.model = list(armaOrder = c(0,0), include.mean = FALSE)
)
# Replicate it into a multispec() element
uspec <- multispec(replicate(2, xspec))
# Define the specification for the DCC model
spec <- dccspec(
# GARCH specification
uspec = uspec,
# DCC specification
dccOrder = c(1,1),
# Distribution, here multivariate normal
distribution = "mvnorm"
)
# Fit the specification to the data
res <- dccfit(spec, data = y)
We can call the variable to see the fitted model:
# Call the variable
res
The output includes the GARCH parameters $\alpha$ and $\beta$ for each one of the assets, and the joint DCC parameters $\zeta$ and $\xi$ respectively.
coef(res)
Using plot( ) on a DCC object¶
Calling the function plot() on a DCC object gives us a menu of options:
Make a plot selection (or 0 to exit):
1. Conditional Mean (vs Realized Returns)
2. Conditional Sigma (vs Realized Absolute Returns)
3. Conditional Covariance
4. Conditional Correlation
5. EW Portfolio Plot with conditional density VaR limits
Some of these plots, like Plot 2, 3 and 4, include information similar to what we plotted in the case of EWMA (The plot is different for different time span):
Plot 2¶
Plot 3¶
Plot 4¶
Exploring a DCC object¶
We can use the @ operator on the DCC object to explore what output is provided:
# Explore DCC object
names(res@mfit)
Let's extract the H matrix, which includes the covariances:
# Extracting the H matrix
H <- res@mfit$H
# Dimensions of the H matrix
dim(H)
The H matrix has three dimensions. You can think of it as a set that contains 8126 2x2 matrices. Each 2x2 matrix is a variance-covariance matrix, and we have one for each time period. We can access it the same way as we would access a two-dimensional matrix, but adding one more element within the square brackets.
If we want to see the first period's matrix:
# First period's covariances
H[,,1]
Let's see the first four covariance matrix estimations:
# First four covariance matrix estimations
print(H[,,1:4])
If we want to compute the conditional correlations of the DCC model, we can extract the variances and covariances from the H matrix. Recall the covariance matrix is a symmetric matrix and the variances are found in the diagonal of the matrix while the covariance will be found in the off-diagonal element:
# Computing the conditional correlations
# Initializing the vector
rhoDCC <- vector(length = n)
# Populate with the correlations
rhoDCC <- H[1,2,] / sqrt(H[1,1,]*H[2,2,])
DCC apARCH and tapARCH (Not covered in class)¶
Let's fit more DCC models using different univariate specification for the estimations of the assets' returns.
First, consider a DCC apARCH:
# DCC apARCH model
# Univariate specification
xspec <- ugarchspec(variance.model = list(model = "apARCH"),
mean.model = list(armaOrder = c(0,0), include.mean = FALSE)
)
# Duplicate the specification using multispec
uspec <- multispec(replicate(2, xspec))
# Create the DCC specification
spec <- dccspec(uspec = uspec,
dccOrder = c(1,1),
distribution = "mvnorm")
# Fit it to the data
res_aparch <- dccfit(spec, data = y)
# Call the object
res_aparch
The output of the DCC apARCH includes the extra apARCH parameters for each individual stock. We can enrich this model further by using a multivariate T distribution as the joint distribution between the returns instead of a multivariate normal:
# DCC tapARCH model
# Univariate specification
xspec <- ugarchspec(variance.model = list(model = "apARCH"),
mean.model = list(armaOrder = c(0,0), include.mean = FALSE)
)
# Duplicate the specification using multispec
uspec <- multispec(replicate(2, xspec))
# Create the DCC specification, replace the multivariate normal by a multivariate T
spec <- dccspec(uspec = uspec,
dccOrder = c(1,1),
distribution = "mvt")
# Fit it to the data
res_taparch <- dccfit(spec, data = y)
# Call the object
res_taparch
Each stock has five individual parameters, three that are part of the standard GARCH estimation (omega, alpha, beta), and two that are particular to the apARCH model (gamma and delta). Additionally, the joint parameters include mshape, which is the estimation for the degrees of freedom of the multivariate T distribution.
Let's extract the estimated variances and covariances to create the conditional correlations vector:
# Creating the conditional correlations vector from the tapARCH model
# Extracting covariances
H <- res_taparch@mfit$H
# Initializing vector
rhoDCCRich <- vector(length = n)
# Fill the covariance matri
rhoDCCRich <- H[1,2,] / sqrt(H[1,1,] * H[2,2,])
DCC with several assets (Not covered in class)¶
Now instead of using two assets, let's build a DCC model using all the assets in our Y data frame:
# DCC model for all assets
# Extracting all the columns from Y except the date
y_all <- subset(Y, select = -c(date, int_date))
# See y
head(y_all)
# Transform to matrix
y_all <- as.matrix(y_all)
# Check dimensions
dim(y_all)
# Build DCC model
# Univariate spec
xspec <- ugarchspec(mean.model = list(armaOrder = c(0,0), include.mean = FALSE))
# Replicate for each asset
uspec <- multispec(replicate(dim(y_all)[2], xspec))
# Build DCC spec
spec <- dccspec(uspec = uspec,
dccOrder = c(1,1),
distribution = "mvnorm")
# Fit the model
dcc_all <- dccfit(spec, data = y_all)
# Call the object
dcc_all
Each stock has three parameters, and there are two joint ones. In total we have 20 parameters. Let's now extract the H matrix and see how it looks:
# Extracting H matrix
H <- dcc_all@mfit$H
# Check the dimensions
dim(H)
# See the first one - 6x6 matrix
H[,,1]
To build the conditional correlations, we have to be careful on specifying which two assets we are interested in. We can see which asset belongs to which column by calling:
# Check which asset is in which column
colnames(y_all)
If we wanted the conditional correlation between 'MSFT' and 'JPM', we would have to write this:
# Conditional correlation between two stocks
r_msft_jpm <- H[1,4,]/sqrt(H[1,1,]*H[4,4,])
This can be annoying since we have to keep track of the numbers for each stock. Alternatively, we can write a short function that solves this problem:
# Writing a function that prevents us from keeping track of numbers
cond_corr <- function(stock1, stock2) {
# Finds the index of each ticker in colnames
index1 <- which(colnames(y_all) == stock1)
index2 <- which(colnames(y_all) == stock2)
# Return the vector operation
return(H[index1, index2, ]/sqrt(H[index1, index1,]*H[index2, index2,]))
}
# Call it on MSFT and JPM
r_msft_jpm_2 <- cond_corr("MSFT", "JPM")
# Compare the two:
head(r_msft_jpm)
head(r_msft_jpm_2)
Now we can easily find the conditional correlation of stocks without worrying of the index. Let's get the conditional correlation for JP Morgan and Citigroup, and plot it over time:
# JPM and C
r_jpm_c <- cond_corr("JPM", "C")
# Plot it
plot(Y$date, r_jpm_c, main = "Conditional correlation for JPM and C",
type = "l", col = "red", las = 1)
Add a horizontal line at the mean:
# Horizontal line at mean
plot(Y$date, r_jpm_c, main = "Conditional correlation for JPM and C",
type = "l", col = "red", las = 1)
abline(h = mean(r_jpm_c), lwd = 2)
Model comparison¶
In this section we will compare both visually and numerically the output of the models we have ran. We will focus in comparing the EWMA and DCC models.
Plotting the Conditional Correlations¶
Let's use matplot() to plot the conditional correlations obtained in EWMA, DCC, and the DCC using apARCH:
# Comparing DCC and EWMA visually
matplot(cbind(EWMArho, rhoDCC, rhoDCCRich), type = "l", lty = 1,
main = "Conditional Correlations for EWMA, DCC, DCC and tapARCH")
legend("bottomright", legend = c("EWMA", "DCC", "DCC tapARCH"), lty = 1, col = 1:3)
We can see the estimations are quite similar. Computing the correlation between the conditional correlations of the three methods shows us a near-perfect correlation:
# Correlation
cor(cbind(EWMArho, rhoDCC, rhoDCCRich))
We see a strong linear relationship between EWMArho and rhoDCC:
# Relationship between EWMArho and rhoDCC
plot(rhoDCC, EWMArho, main = "Relationship between EWMArho and rhoDCC")
Let's use colMeans() to obtain the mean of each series:
# Mean of each
colMeans(cbind(EWMArho, rhoDCC, rhoDCCRich))
# Standard deviation of each
sd(EWMArho)
sd(rhoDCC)
sd(rhoDCCRich)
The mean of the conditional correlation for the three models is very similar. We see that the EWMA model presents the largest variability.
Now we plot only the EWMA and standard DCC models:
# EWMA and standard DCC
matplot(cbind(EWMArho, rhoDCC), type = "l", lty = 1,
main = "Conditional Correlations for EWMA and DCC")
legend("bottomright", legend = c("EWMA", "DCC"), lty = 1, col = 1:2)
Recap¶
In this seminar we have covered:
- Multivariate volatility modelling in R
- Implementing an EWMA model
- Extracting and plotting conditional variances, covariances, and correlations for EWMA
- Fitting different specifications of DCC models
- Extracting and plotting conditional variances, covariances, and correlations for DCC
- Comparing and plotting models
Some new functions used:
cov()upper.tri()multispec()dccspec()dccfit()
For more discussion on the material covered in this seminar, refer to Chapter 3: Multivariate volatility modeling on Financial Risk Forecasting by Jon Danielsson.
Acknowledgements: Thanks to Alvaro Aguirre and Yuyang Lin for creating these notebooks
© Jon Danielsson, 2022