Chapter 3. Multivariate Volatility Models (in R/Python)


Copyright 2011 - 2019 Jon Danielsson. This code is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This code is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. The GNU General Public License is available at: https://www.gnu.org/licenses/.

The original 2011 R code will not fully work on a recent R because there have been some changes to libraries. At least two R packages support estimating GARCH style models, some are old and not maintained, rmgarch by Alexios Ghalanos is regularly maintained, and what is used below,


Listing 3.1/3.2: Download stock prices in R
Last updated August 2019

p = read.csv('stocks.csv')
y=apply(log(p),2,diff)     # calculate returns
y = y[,1:2]                # consider first two stocks
y[,1] = y[,1]-mean(y[,1])  # subtract mean
y[,2] = y[,2]-mean(y[,2])
TT = dim(y)[1]
		
Listing 3.1/3.2: Download stock prices in Python
Last updated June 2018

import numpy as np
p = np.loadtxt('stocks.csv',delimiter=',',skiprows=1)
p = p[:,[0,1]]                                        # consider first two stocks
y = np.diff(np.log(p), n=1, axis=0)*100               # calculate returns
y[:,0] = y[:,0]-np.mean(y[:,0])                       # subtract mean
y[:,1] = y[:,1]-np.mean(y[:,1])
T = len(y[:,0])
		

Listing 3.3/3.4: EWMA in R
Last updated August 2019

## create a matrix to hold covariance matrix for each t
EWMA = matrix(nrow=TT,ncol=3)
lambda = 0.94
S = cov(y)                                             # initial (t=1) covar matrix
EWMA[1,] = c(S)[c(1,4,2)]                              # extract var and covar
for (i in 2:dim(y)[1]){
      S = lambda*S+(1-lambda)*  y[i-1,] %*% t(y[i-1,])
      EWMA[i,] = c(S)[c(1,4,2)]
}
EWMArho = EWMA[,3]/sqrt(EWMA[,1]*EWMA[,2])             # calculate correlations
print(head(EWMArho))
print(tail(EWMArho))
		
Listing 3.3/3.4: EWMA in Python
Last updated June 2018

EWMA = np.full([T,3], np.nan)
lmbda = 0.94
S = np.cov(y, rowvar = False)
EWMA[0,] = S.flatten()[[0,3,1]]
for i in range(1,T):
    S = lmbda * S + (1-lmbda) * np.transpose(np.asmatrix(y[i-1]))* np.asmatrix(y[i-1])
    EWMA[i,] = [S[0,0], S[1,1], S[0,1]]
EWMArho = np.divide(EWMA[:,2], np.sqrt(np.multiply(EWMA[:,0],EWMA[:,1])))
print(EWMArho)
		

Listing 3.5/3.6: GOGARCH in R
Last updated August 2019

library(rmgarch)
spec = gogarchspec(mean.model = list(armaOrder = c(0, 0),
    include.mean =FALSE),
    variance.model = list(model = "sGARCH",
    garchOrder = c(1,1)) ,
    distribution.model =  "mvnorm"
)
fit = gogarchfit(spec = spec, data = y)
show(fit)
		
Listing 3.5/3.6: OGARCH in Python
Last updated June 2018

## Python does not have a proper OGARCH package at present
		

Listing 3.7/3.8: DCC in R
Last updated August 2019

library(rmgarch)
xspec = ugarchspec(mean.model = list(armaOrder = c(0, 0), include.mean = FALSE))
uspec = multispec(replicate(2, xspec))
spec = dccspec(uspec = uspec, dccOrder = c(1, 1), distribution = 'mvnorm')
res = dccfit(spec, data = y)
H=res@mfit$H
DCCrho=vector(length=dim(y)[1])
for(i in 1:dim(y)[1]){
    DCCrho[i] =  H[1,2,i]/sqrt(H[1,1,i]*H[2,2,i])
}
		
Listing 3.7/3.8: DCC in Python
Last updated June 2018

## Python does not have a proper DCC package at present
		

Listing 3.9/3.10: Sample statistics in R
Last updated August 2019

matplot(cbind(EWMArho,DCCrho),type='l',las=1,lty=1,col=2:3,ylab="")
mtext("Correlations",side=2,line=0.3,at=1,las=1,cex=0.8)
legend("bottomright",c("EWMA","DCC"),lty=1,col=2:3,bty="n",cex=0.7)
		
Listing 3.9/3.10: Correlation comparison in Python
Last updated June 2018

## Python does not have a proper OGARCH/DCC package at present