# Realized Garch model with dynamical conditional correlation

The accompanying report is here

The package rGARCH_DCC is included here and it is priliminary and hopefully I can make it user friendly and will publish to my github.

## Realized GARCH-DCC Model

In this project we will study the realized GARCH-DCC model. This model is designed to do forecast of volatility and correlation among symbols. The model takes into high frequency data.

Specifically, we let $r^i_t$ to denote the log return $\log P^i_t - \log P^i_{t-1}$ for stock $i$, $\sigma^i_t$ to be the conditional variance of stock $i$, $z_t^i$ to be the shock for stock $i$ and $C_t$ be the conditional correlation matrix among $z_t^i$. The model can be described as

\begin{aligned} &r_t^i = \mu^i + \sigma_t^i z_t^i \\ &\log(\text{RV}_t^i) = \xi^i + \phi^i \log(\sigma_{t}^i)^2 + \tau^i_1 z_t^i + \tau^i_2 ((z_t^i)^2-1) + u_t^i\\ &(\sigma_{t+1}^i)^2 = \omega^i + \beta^i (\sigma_{t}^i)^2 + \gamma^i \text{RV}_t^i \\ &C_t = \text{diag}(Q_t)^{-1/2} Q_t \text{diag}(Q_t)^{-1/2} \label{eq:section2:normalization} \\ &Q_{t+1} = S(1-a-b) + a (z_t z_t^T) + bQ_t \label{eq:section2:covariance} \\ &z_t \sim \mathcal{N}(0, C_t) \\ &u_t \sim \mathcal{N}(0, \Sigma_u) . \end{aligned}

The first equation describes the return using the volatility scale and the constant expected mean $\mu^i$. The second equation describes the “measurement” mechanics: the realized variance depends on the hidden true volatility and the shocks as well as a measurement noise $u_t$. These $u_t$ have zero pairwise correlation but have different scales for different stock $i$ and thus $\Sigma_u$ is a diagonal matrix. The log transformation on the measurement equation is to make the noise $u_t$ similar to Normal distribution. The third equation is the updating equation for the volatility $\sigma_t^i$. Notice, the usual GARCH model uses return squared to update the volatility but here, we use the more accurate estimate, the realized variance. The $Q_{t}$ is the sample estimate of the conditional covariance matrix of $z_t$ and the update equation is mean revert to the sample estimation $S$. The normalization ensures that the covariance matrix for the shocks is normalized.

import numpy as np
import pandas as pd
import pandas_market_calendars as mcal
import datetime

from scipy.optimize import minimize
from sklearn.linear_model import LinearRegression
from arch import arch_model

import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator
import matplotlib.dates as mdates
from matplotlib.ticker import FormatStrFormatter
from statsmodels.graphics.gofplots import qqplot
import matplotlib

from rGARCH_DCC import *

font = { "family": "sans",
"weight": "normal",
'size'  : 12}

matplotlib.rc('font', **font)


We first load data from https://firstratedata.com/free-intraday-data. The data contains AAPL, AMZN, FB, MSFT, TSLA and SP500 one minute tick data in 2019.

FOLDER = "/Users/bunny_home/Documents/Cunwei/stockData/"
AAPL_FILE = FOLDER + "AAPL_FirstRateDatacom1.txt"
AMZN_FILE = FOLDER + "AMZN_FirstRateDatacom1.txt"
FB_FILE   = FOLDER + "FB_FirstRateDatacom1.txt"
MSFT_FILE = FOLDER + "MSFT_FirstRateDatacom1.txt"
TSLA_FILE = FOLDER + "TSLA_FirstRateDatacom1.txt"
SPX_FILE  = FOLDER + "SPX_Firstratedata1.txt"

FILES = {"AAPL" : AAPL_FILE,
"AMZN" : AMZN_FILE,
"FB"   : FB_FILE,
"MSFT" : MSFT_FILE,
"TSLA" : TSLA_FILE,
"SP500": SPX_FILE}


Since there are missing ticks in the data, we then construct time ticks frame in one minute frequency and left join with the data. This ensures the data all in one tick frequency and can be joined conveniently.

def get_trading_index(start, end):
nyse = mcal.get_calendar('NYSE')
days = nyse.valid_days(start_date = start, end_date = end)
date_index = []

for d in days:
#print (d)
early_close = [datetime.date(d.year, 7,3),
datetime.date(d.year, 11,29),
datetime.date(d.year, 12,24)]
if d.date() in early_close :
open_time = pd.Timestamp(year=d.year, month = d.month,
day = d.day, hour = 9, minute = 30)
close_time = pd.Timestamp(year=d.year, month = d.month,
day = d.day, hour = 13, minute = 0)
else:
schedule = nyse.schedule(start_date = d.date(), end_date = d.date(), tz = nyse.tz.zone)
open_time = schedule["market_open"].replace(tzinfo=None)
close_time = schedule["market_close"].replace(tzinfo=None)

date_index += pd.date_range(open_time, close_time, freq = "1T")

return pd.Series(date_index)

df.columns = ["time", "open", "high", "low", "close", "volume"][:len(df.columns)]
df["time"] = pd.to_datetime(df["time"], format= "%Y-%m-%d %H:%M:%S")
df = df.set_index("time")

dummy = dummy.set_index("time")

df = dummy.join(df)
df = df.drop("dummy", axis = 1)
return df


DATA = {}
print ("Generating time ticks ...")
print ("Ticks number: %d" % (len(TRADING_MINUTES)))
for F in FILES:
DATA[F] = data
print ("%s\t length: %d, done" %(F, data.shape))

Generating time ticks ...
Ticks number: 97601
AAPL	 length: 97601, done
AMZN	 length: 97601, done
FB	 length: 97601, done
MSFT	 length: 97601, done
TSLA	 length: 97601, done
SP500	 length: 97601, done

combine = []
keys = []
for s in DATA:
combine.append(DATA[s]["close"])
keys.append(s)
prices = pd.concat(combine, axis = 1, keys = keys)
prices

AAPL AMZN FB MSFT TSLA SP500
time
2019-01-02 09:30:00 154.7800 1466.9690 129.7950 99.010 305.9600 2470.40
2019-01-02 09:31:00 155.1597 1469.0000 130.1100 99.250 308.8050 2470.80
2019-01-02 09:32:00 154.8073 1470.6300 130.4329 99.220 306.9370 2471.26
2019-01-02 09:33:00 154.6700 1470.8000 130.4700 99.355 305.1100 2469.64
2019-01-02 09:34:00 154.7500 1473.9806 130.6524 99.410 303.8800 2470.11
... ... ... ... ... ... ...
2019-12-30 15:56:00 291.4700 1846.7900 204.4700 157.585 414.6150 3220.75
2019-12-30 15:57:00 291.4477 1845.5500 204.4200 157.520 414.4712 3219.83
2019-12-30 15:58:00 291.3800 1845.5100 204.3000 157.410 414.3600 3218.84
2019-12-30 15:59:00 291.6470 1847.1800 204.4100 157.680 414.6200 3222.06
2019-12-30 16:00:00 291.6100 1846.8900 204.3600 157.590 414.4700 NaN

97601 rows × 6 columns

### Clean Nans and Create Log Returns

Since some data are missing for specific ticks and we should clean them. From the nan counts, we see that most of the nans are in 2019-09-30 and since each day should contain 391 points from 9:30am to 4:00pm inclusive, this means the data are all missing in 2019-09-30. We just drop that day.

def get_Nan_days(data):
nans = data.isna()
nans = nans.groupby(nans.index.date).sum()
return nans.loc[(nans > 1).sum(axis = 1)>1]

nan_days = get_Nan_days(prices)
nan_days

AAPL AMZN FB MSFT TSLA SP500
2019-08-28 0 3 0 0 2 1
2019-09-30 0 391 391 391 391 1
2019-10-17 2 2 0 0 0 1
2019-10-21 4 10 4 4 4 1
delete_idx = (prices.index.date==datetime.date(2019,9,30))
prices = prices.loc[~delete_idx]
prices = prices.fillna(method="ffill")
get_Nan_days(prices)

AAPL AMZN FB MSFT TSLA SP500
prices.shape

(97210, 6)


we then just construct the log returns in one minute frequency.

# The log return:

returns = np.log(prices).groupby(prices.index.date).diff()
returns

AAPL AMZN FB MSFT TSLA SP500
time
2019-01-02 09:30:00 NaN NaN NaN NaN NaN NaN
2019-01-02 09:31:00 0.002450 0.001384 0.002424 0.002421 0.009256 0.000162
2019-01-02 09:32:00 -0.002274 0.001109 0.002479 -0.000302 -0.006067 0.000186
2019-01-02 09:33:00 -0.000887 0.000116 0.000284 0.001360 -0.005970 -0.000656
2019-01-02 09:34:00 0.000517 0.002160 0.001397 0.000553 -0.004039 0.000190
... ... ... ... ... ... ...
2019-12-30 15:56:00 0.000103 0.000846 0.000489 0.000667 0.000084 0.000096
2019-12-30 15:57:00 -0.000077 -0.000672 -0.000245 -0.000413 -0.000347 -0.000286
2019-12-30 15:58:00 -0.000232 -0.000022 -0.000587 -0.000699 -0.000268 -0.000308
2019-12-30 15:59:00 0.000916 0.000904 0.000538 0.001714 0.000627 0.001000
2019-12-30 16:00:00 -0.000127 -0.000157 -0.000245 -0.000571 -0.000362 0.000000

97210 rows × 6 columns

returns = returns.dropna()
returns.to_csv("logReturns.csv")


### Test Train split

test_index = returns.index.date >= datetime.date(2019,12,1)
train = returns[~test_index]
test = returns[test_index]


## Construct Realized Variance

In calculating the realized variance, the frequency of data points need to be chosen. Since our data is spaced in one minute frequence, we should choose intervals larger than that. We calculate realized Variance for different spacing and plot the realized variance vs the spacing. Then we would choose the spacing corresponding to the minimum value as the optimal spacing.

def RV_singleDay(returns, spacing):
res = 0
for start_tick in range(spacing):
res += RV_singleDay_test(returns, spacing, start_tick)
return res/spacing

def RV_singleDay_test(returns, spacing, start_tick):
sums = returns.rolling(spacing).sum()
if len(sums.shape) == 1:
return np.sum((sums.iloc[start_tick+spacing::spacing])**2, axis = 0)
return np.sum((sums.iloc[start_tick+spacing::spacing,:])**2, axis = 0)

def getRVs(returns, spacing, start_idx = None):

def RVsingle(x):
return RV_singleDay(x, spacing)
if start_idx is None:
return returns.groupby(returns.index.date).apply(RVsingle)

def RVsingle_test(x):
return RV_singleDay_test(x, spacing, start_idx)

return returns.groupby(returns.index.date).apply(RVsingle_test)

def multipleRVs(returns, spacing):
res = {}
for s in returns.columns:
res[s] = getRVs(returns[s], spacing[s])
res = pd.DataFrame.from_dict(res, orient = "columns")
return res

def test_spacing(data, max_spacing, step = 1,  start_idx = None):

res = {}
for s in range(1,max_spacing+1, step):
res[s] = getRVs(train, s, start_idx).mean(axis = 0)
res = pd.DataFrame.from_dict(res, orient = "index")
res["spacing"] = res.index

return res


RV_tuning = test_spacing(train, 100, start_idx = 0, step=2)
fig, axs = plt.subplots(2, 3, figsize=(16,10))
for i in range(2):
for j in range(3):
ax = axs[i][j]
ax.scatter(RV_tuning.iloc[:,-1], RV_tuning.iloc[:,i*3+j], s = 60)
#ax.set_yscale("log")
ax.yaxis.set_major_formatter(FormatStrFormatter('%1.1e'))
ax.title.set_text(RV_tuning.columns[i*3+j])
plt.savefig("RVspacing.png")
plt.show() From the above plot, we choose the smallest minute after which the realized variance stablize and this corresponds to the values listed below and we can trasnform the data to get daily realized variance.

SPACING = {
"AAPL" : 40,
"AMZN" : 30,
"FB"   : 10,
"MSFT" : 10,
"TSLA" : 40,
"SP500": 10
}

def collapse_daily(returns, spacing=SPACING):
RVs = multipleRVs(returns, spacing)
returns = returns.groupby(returns.index.date).sum()
index = pd.MultiIndex.from_product([returns.columns, ["RV", "return"] ] )
data = pd.DataFrame(np.zeros((RVs.shape, RVs.shape*2)),
index = RVs.index, columns = index)
for j in RVs.columns:
data.loc[:, (j,"return")] = returns[j]
data.loc[:, (j, "RV")] = RVs[j]
return data


train_daily = collapse_daily(train)
test_daily = collapse_daily(test)

index = (train_daily["TSLA"]["RV"] == 0)
train_daily = train_daily[~index]
train_daily

AAPL AMZN FB MSFT TSLA SP500
RV return RV return RV return RV return RV return RV return
2019-01-02 0.000276 0.018817 0.000566 0.047986 0.000478 0.044343 0.000252 0.021878 0.000517 0.012795 0.000143 0.014631
2019-01-03 0.000334 -0.014916 0.000346 -0.017998 0.000446 -0.021329 0.000419 -0.025744 0.000573 -0.024370 0.000222 -0.019205
2019-01-04 0.000240 0.025416 0.000369 0.031095 0.000224 0.023634 0.000381 0.025587 0.000389 0.032885 0.000171 0.022167
2019-01-07 0.000355 -0.002600 0.000272 0.019787 0.000296 0.001024 0.000194 0.005699 0.000679 0.040694 0.000068 0.006271
2019-01-08 0.000114 0.007528 0.000419 -0.002183 0.000299 0.012709 0.000191 -0.003128 0.001019 -0.021433 0.000087 0.001863
... ... ... ... ... ... ... ... ... ... ... ... ...
2019-11-22 0.000038 -0.002365 0.000027 0.004836 0.000064 -0.000905 0.000034 -0.003204 0.000223 -0.020516 0.000007 -0.000595
2019-11-25 0.000022 0.012085 0.000052 0.008161 0.000050 0.000200 0.000012 0.006900 0.000164 -0.014727 0.000004 0.003998
2019-11-26 0.000012 -0.006490 0.000035 0.009007 0.000034 -0.004639 0.000016 0.003198 0.000120 -0.016852 0.000006 0.002352
2019-11-27 0.000011 0.008658 0.000036 0.011019 0.000085 0.011351 0.000013 0.000722 0.000083 0.004661 0.000002 0.002683
2019-11-29 0.000007 0.001948 0.000030 -0.010237 0.000081 0.001410 0.000015 -0.004232 0.000027 -0.000545 0.000003 -0.001396

229 rows × 12 columns

We can visualize the realized variance and the squared log returns.

for s in DATA:
plt.figure()
plt.title(s)
plt.plot((train_daily[s]["return"]**2),label="r2")
plt.plot(train_daily[s]["RV"],label="RV")
plt.legend(loc = "best")
plt.show()      We then fit the individual realized Garch

GARCHs = {}
PARAMs = {}
scale = 1e3
for symbol in DATA:
print ("fit symbol "+ str(symbol) + " ...")
returns = train_daily[symbol]["return"]
RVs = train_daily[symbol]["RV"]
rl_garch = RealizedGARCH(returns, RVs, scale = scale)
rl_garch.fit(verbose = 0)
GARCHs[symbol] = rl_garch
params = pd.concat( (rl_garch.params, rl_garch.measure_params) )
PARAMs[symbol] = params

PARAMs = pd.DataFrame.from_dict(PARAMs, orient="columns")
PARAMs

fit symbol AAPL ...
best L1 = -650.558
fit symbol AMZN ...
best L1 = -668.673
fit symbol FB ...
best L1 = -711.490
fit symbol MSFT ...
best L1 = -636.488
fit symbol TSLA ...
best L1 = -816.116
fit symbol SP500 ...
best L1 = -498.618

AAPL AMZN FB MSFT TSLA SP500
omega 0.000035 0.000013 0.000012 0.000011 0.000241 0.000003
beta 0.306265 0.804552 0.932296 0.665771 0.000000 0.415299
gamma 0.498232 0.112582 0.000000 0.248973 0.762372 0.502016
mu 0.001039 0.000049 0.000246 0.000175 0.000450 0.000387
xi 0.221838 7.076395 153.648558 2.076153 -1.130535 -0.787398
phi 1.084924 1.846474 18.791195 1.256567 0.931592 0.943701
tau1 -0.101627 -0.092318 -0.105260 -0.072026 -0.026562 -0.174259
tau2 0.244832 0.176118 0.134127 0.156663 0.149789 0.153689
sigmaU 0.650246 0.584097 0.490212 0.444128 0.566098 0.554517
def realizedGarch_log_likelihood():
Likelihoods = {}
for symbol in GARCHs:
Likelihoods[symbol] = GARCHs[symbol].train_log_likelihood()
Likelihoods = pd.DataFrame.from_dict(Likelihoods, orient= "index")
Likelihoods.columns = ["log(L)"]
Likelihoods = Likelihoods.transpose()
return Likelihoods

def get_garchs():
original_garches = {}
for symbol in DATA:
am = arch_model(1000*train_daily[symbol]["return"], mean='Constant', vol='garch')
res = am.fit(disp="off")
original_garches[symbol] = res

return original_garches

def arch_log_likelihood(arches, scale = 1):
ans = {}
for symbol in arches:
arch_res = arches[symbol]
vol = (arch_res.conditional_volatility/scale)**2
resid = arch_res.std_resid
loglikelihood = -0.5*np.sum(np.log(2*np.pi*vol) + resid**2)
ans[symbol] = loglikelihood
ans = pd.DataFrame.from_dict(ans, orient="index")
ans.columns=["log(L)"]
ans = ans.transpose()
return ans



We fit our reailized Garch model and the classical Garch separately and compare them by the likelihoods. Our realized model achieves better value in AAPL, MSF, TSLA and SP500. The fit for AMZN and FB seems not working.

realizedGarch_log_likelihood()

AAPL AMZN FB MSFT TSLA SP500
log(L) 720.881208 702.766017 659.948803 734.951069 555.322595 872.821334
arch_log_likelihood(get_garchs(), scale = 1000)

AAPL AMZN FB MSFT TSLA SP500
log(L) 714.535985 708.639553 662.078304 734.394404 549.42429 867.143143

From the below plot its clear why the FB and AMZN do not show higher likelihood. The fitted conditional variance $\sigma^2$ doest follow the realized variance. Especially, the FB model has $\beta = 0$ and does not take into account of realized variance in the model. The AMZN has very small coefficient for $\beta$ and is thus very insensitive to the change of the market (the rise or fall of realized variance).

fig, axs = plt.subplots(2, 3, figsize=(16,10))
for i in range(2):
for j in range(3):
ax = axs[i][j]
symbol = RV_tuning.columns[i*3+j]
GARCHs[symbol].plot(ax= ax)
ax.xaxis.set_major_locator(MaxNLocator(5))
ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y-%m"))
ax.xaxis.set_minor_formatter(mdates.DateFormatter("%Y-%m"))
ax.yaxis.set_major_formatter(FormatStrFormatter('%1.1g'))
ax.title.set_text(symbol) def plot_z(Garchs):
fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(16,8))
fig.suptitle(r"QQ plot for $z_t$")
axs = axs.reshape(-1)
for j, stock in enumerate(Garchs):
ax1 = axs[j]
ax1.set_title(stock)
qqplot(Garchs[stock].std(), line = "s", ax = ax1)
ax1.set_xlabel("")
plt.savefig("zt.png")
def plot_u(Garchs):
fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(16,8))
fig.suptitle(r"QQ plot for $u_t$")
axs = axs.reshape(-1)
for j, stock in enumerate(Garchs):
ax1 = axs[j]
ax1.set_title(stock)
measure_RV = Garchs[stock].measure_RV()
diff = np.log(measure_RV) - np.log(Garchs[stock].RVs_original[1:])
qqplot(diff, line = "q", ax = ax1)
ax1.set_xlabel("")
plt.savefig("ut.png")


We plot fitted residual for $z_t$ and measurement linear regression residual $u_t$ in qq plot. The plot shows that they almost follow the normal distribution. As a result, the assumption is mildly satisfied.

plot_z(GARCHs)
plot_u(GARCHs)  ## DCC correlations

We fit the rGARCH-DCC model with the DCC employed to describe the correlation among the shocks $z_t$ among the symbols.

rGarch = RGARCH_DCC(train_daily, garch_scale = 1000)
rGarch.fit(verbose = 0)
rGarch.params()

fit AAPL ...
best L1 = -650.558
fit AMZN ...
best L1 = -668.673
fit FB ...
best L1 = -711.490
fit MSFT ...
best L1 = -636.488
fit TSLA ...
best L1 = -816.116
fit SP500 ...
best L1 = -498.618
iteration: 1 	 log(L): 73.840 	 a: 0.38595 	 b:0.61405
iteration: 2 	 log(L): 680.975 	 a: 0.3732 	 b:0.31926
iteration: 3 	 log(L): 672.804 	 a: 0.3509 	 b:0.58848
iteration: 4 	 log(L): 515.326 	 a: 0.31236 	 b:0.68764
iteration: 5 	 log(L): 543.880 	 a: 0.30289 	 b:0.69711
iteration: 6 	 log(L): 579.609 	 a: 0.2883 	 b:0.7117
iteration: 7 	 log(L): 572.286 	 a: 0.29162 	 b:0.70838
iteration: 8 	 log(L): 625.487 	 a: 0.25989 	 b:0.74011
iteration: 9 	 log(L): 533.779 	 a: 0.30642 	 b:0.69358
iteration: 10 	 log(L): 649.872 	 a: 0.19871 	 b:0.80129
iteration: 11 	 log(L): 687.172 	 a: 0.30417 	 b:0.62254
iteration: 12 	 log(L): 687.863 	 a: 0.30861 	 b:0.59438
iteration: 13 	 log(L): 687.865 	 a: 0.30947 	 b:0.59154
iteration: 14 	 log(L): 687.865 	 a: 0.30952 	 b:0.59176
iteration: 15 	 log(L): 687.865 	 a: 0.30948 	 b:0.59181
iteration: 16 	 log(L): 687.865 	 a: 0.30948 	 b:0.59181
rGarch parameters:
AAPL      AMZN          FB      MSFT      TSLA     SP500
omega   0.000035  0.000013    0.000012  0.000011  0.000241  0.000003
beta    0.306265  0.804552    0.932296  0.665771  0.000000  0.415299
gamma   0.498232  0.112582    0.000000  0.248973  0.762372  0.502016
mu      0.001039  0.000049    0.000246  0.000175  0.000450  0.000387
xi      0.221838  7.076395  153.648558  2.076153 -1.130535 -0.787398
phi     1.084924  1.846474   18.791195  1.256567  0.931592  0.943701
tau1   -0.101627 -0.092318   -0.105260 -0.072026 -0.026562 -0.174259
tau2    0.244832  0.176118    0.134127  0.156663  0.149789  0.153689
sigmaU  0.650246  0.584097    0.490212  0.444128  0.566098  0.554517

DCC parameters:
a    0.309478
b    0.591807
dtype: float64

def plot_predict(RGARCH, test, n_sample = 1000, low_q = 0.025, high_q = 0.975):
horizon = test.shape
sigma2_pred, RV_pred = RGARCH.predict_horizon(
horizon,
n_sample = n_sample,
low_q = low_q, high_q = high_q)
sigma2_mean = pd.DataFrame(sigma2_pred[1,...])
sigma2_high = pd.DataFrame(RV_pred[2,...])
sigma2_low  = pd.DataFrame(RV_pred[0,...])

sigma2_mean.index = test.index
sigma2_high.index = test.index
sigma2_low.index  = test.index

fig, axs = plt.subplots(2, 3, figsize=(16,10))
fig.suptitle("prediction on in-sample and out-of-sample data")

axs = axs.reshape(-1)

for j, stock in enumerate(RGARCH.stocks):
test_returns = test[stock]["RV"]
ax = axs[j]
RGARCH.GARCHs[stock].plot(ax=ax)
ax.plot(test_returns, color ='#1f77b4')
ax.plot(sigma2_mean.iloc[:,j], color = "darkorange")
ax.fill_between(sigma2_mean.index, sigma2_low.iloc[:,j],
sigma2_high.iloc[:,j], color = "darkgrey")
ax.axvline(x=test.index, color = "k", linestyle = '--')
ax.xaxis.set_major_locator(MaxNLocator(5))
ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y-%m"))
ax.xaxis.set_minor_formatter(mdates.DateFormatter("%Y-%m"))
ax.yaxis.set_major_formatter(FormatStrFormatter('%1.1g'))
ax.title.set_text(stock)
return fig


In this plot, we plot the fitted $\sigma^2$ and the predicted $\sigma^2$ in December. The gray area is the $95\%$ confidence interval for the realized variance. The plots shows that the predicted realized variance confidence interval contains the true hidden test data. This indicates a successful forecast.

fig = plot_predict(rGarch, test_daily)
fig.savefig("train_test.png") We design a function that could use rolling bootstrap to estimate the confidence interval of the estimated parameters.

def rolling_bootstrap(train, n_bootstrap):
length = train.shape - n_bootstrap + 1
GARCH_params = None
DCC_params = None
count = 0
for i in range(n_bootstrap):
garchDCC = RGARCH_DCC(train.iloc[i:i+length,:],
garch_scale = 1000)
garchDCC.fit(verbose = 0)
if GARCH_params is None:
shape = [n_bootstrap] + list(garchDCC.GARCH_PARAMs.shape)
GARCH_params = np.zeros(shape=shape)
if DCC_params is None:
DCC_params = np.zeros(shape = (n_bootstrap, 2))
if garchDCC.DCC_PARAMs["b"] <= 1.0e-5 :
continue
GARCH_params[i,:,:] = garchDCC.GARCH_PARAMs
DCC_params[count,:] = garchDCC.DCC_PARAMs
count +=1

GARCH_params_std = np.std(GARCH_params, axis = 0)
DCC_params_std = np.std(DCC_params[:count,:], axis = 0)

GARCH_params_std = pd.DataFrame(GARCH_params_std)
GARCH_params_std.index = garchDCC.GARCH_PARAMs.index
GARCH_params_std.columns = garchDCC.GARCH_PARAMs.columns

DCC_params_std = pd.Series(DCC_params_std)
DCC_params_std.index = garchDCC.DCC_PARAMs.index
print (count)
return GARCH_params_std, DCC_params_std

rGarch.GARCH_PARAMs

AAPL AMZN FB MSFT TSLA SP500
omega 0.000035 0.000013 0.000012 0.000011 0.000241 0.000003
beta 0.306265 0.804552 0.932296 0.665771 0.000000 0.415299
gamma 0.498232 0.112582 0.000000 0.248973 0.762372 0.502016
mu 0.001039 0.000049 0.000246 0.000175 0.000450 0.000387
xi 0.221838 7.076395 153.648558 2.076153 -1.130535 -0.787398
phi 1.084924 1.846474 18.791195 1.256567 0.931592 0.943701
tau1 -0.101627 -0.092318 -0.105260 -0.072026 -0.026562 -0.174259
tau2 0.244832 0.176118 0.134127 0.156663 0.149789 0.153689
sigmaU 0.650246 0.584097 0.490212 0.444128 0.566098 0.554517
rGarch.DCC_PARAMs

a    0.309478
b    0.591807
dtype: float64

def getRV_data(data, weights):
indexs = data.columns.get_level_values(0).unique()
returns = np.zeros( (data.shape, len(indexs)) )
for i,symbol in enumerate(indexs):
returns[:,i] = data[symbol]["return"]
sigma2 = returns.dot(weights)**2
sigma2 = pd.Series(sigma2)
sigma2.index = data.index
return sigma2

def plot_portfolio_RV(train, test, model, weights, low_q = 0.025,
high_q = 0.975, ax = None):
model_sigma2 = model.get_portfolio_sigma2(weights)
train_sigma2 = getRV_data(train_daily, weights)
test_sigma2  = getRV_data(test_daily, weights)
pred_sigma2  = rGarch.predict_horizon_portfolio(
test_sigma2.shape,
weights,
low_q = low_q, high_q = high_q)
pred_sigma2 = pd.DataFrame(pred_sigma2.T)
pred_sigma2.index = test_daily.index
if ax is None:
ax = plt.gca()
ax.plot(train_sigma2, color = "#1f77b4", label = r"$r^2$")
ax.plot(test_sigma2, color = "#1f77b4")
ax.plot(model_sigma2, color = "darkorange", label = r"model $\sigma^2$")
ax.plot(pred_sigma2.iloc[:,1], color = "darkorange")
ax.fill_between(pred_sigma2.index, pred_sigma2.iloc[:,0],
pred_sigma2.iloc[:,2], color = "darkgrey")
ax.axvline(x=test.index, color = "k", linestyle = '--')
ax.legend(loc="best")
return

def plot_six_pack(train, test, model, weights = None,
low_q = 0.025, high_q = 0.975):

symbols = train.columns.get_level_values(0).unique()
n_symbols = len(symbols)

if weights is None:
weights = np.random.uniform(size=(6, n_symbols))
weights = weights/np.sum(weights, axis=1)[:,None]

fig, axs = plt.subplots(2, 3, figsize=(16,10))
fig.suptitle(r"portfolio $\sigma^2$ estimation")

print (weights)
print (weights.sum(axis=1))
axs = axs.reshape(-1)
for i in range(6):
ax = axs[i]
plot_portfolio_RV(train, test, model, weights[i,:],
low_q, high_q, ax = ax)
ax.xaxis.set_major_locator(MaxNLocator(5))
ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y-%m"))
ax.xaxis.set_minor_formatter(mdates.DateFormatter("%Y-%m"))
ax.yaxis.set_major_formatter(FormatStrFormatter('%1.1g'))
ax.title.set_text(f"portfolio {i+1:d}")

return fig



We then construct six random portfolio consisting of the 6 symbols we have. We do this because we want to test whether the correlation matrix is estimated well. The portfolio volatility needs the correlation matrix and individual volatility. Below we plotted the fitted portfolio $\sigma^2$ and the prected $\sigma^2$ in December. The blue line is the portfolio realized return squared.

fig = plot_six_pack(train_daily, test_daily, rGarch)
plt.savefig("six_pack.png")

[[0.16405074 0.17298219 0.19924118 0.0351977  0.247057   0.18147119]
[0.12837752 0.19303799 0.01225164 0.04758079 0.17934214 0.43940992]
[0.10489163 0.25960325 0.0583599  0.26450911 0.0773851  0.23525101]
[0.09721203 0.31495958 0.05358826 0.04046787 0.3348108  0.15896146]
[0.14997454 0.26229132 0.18209948 0.11644926 0.19562806 0.09355733]
[0.0638383  0.15728596 0.12019817 0.25126211 0.2573531  0.15006237]]
[1. 1. 1. 1. 1. 1.] Categories:

Updated: