## Theory

GARCH-in-mean, or GARCH-M, models extend the GARCH framework by allowing the returns to depend upon a volatility component. The Beta-t-EGARCH M model includes this term $\phi$ as follows:

$y_{t} = \mu + \phi\exp\left(\lambda_{t\mid{t-1}}/2\right) + \exp\left(\lambda_{t\mid{t-1}}/2\right)\epsilon_{t}$

$\lambda_{t\mid{t-1}} = \alpha_{0} + \sum^{p}_{i=1}\alpha_{i}\lambda_{t-i} + \sum^{q}_{j=1}\beta_{j}u_{t-j}$

$\epsilon_{t} \sim t_{\nu}$

## PyFlux

First let us load some financial time series data from Yahoo Finance:

import numpy as np
import pyflux as pf
import pandas as pd
from pandas.io.data import DataReader
from datetime import datetime
import matplotlib.pyplot as plt
%matplotlib inline

aapl = DataReader('AAPL',  'yahoo', datetime(2006,1,1), datetime(2016,3,10))
returns = pd.DataFrame(np.diff(np.log(aapl['Adj Close'].values)))
returns.index = aapl.index.values[1:aapl.index.values.shape[0]]
returns.columns = ['AAPL Returns']

plt.figure(figsize=(15,5));
plt.plot(returns.index,returns);
plt.ylabel('Returns');
plt.title('AAPL Returns');


Let’s fit an EGARCH-M(1,1) model. Here we’ll use BBVI to fit the model. BBVI for this model type starts with a MAP estimate as a starting point, then follows stochastic gradients after this point (a BBVI estimate may differ from a modal MAP estimate):

model.adjust_prior([0,1,2,3,4,5], pf.Uniform())
model = pf.EGARCHM(returns,p=1,q=1)
x = model.fit('BBVI',iterations=400)
x.summary()

10% done : ELBO is 6496.60179525
20% done : ELBO is 6499.82624111
30% done : ELBO is 6509.97490232
40% done : ELBO is 6464.00566071
50% done : ELBO is 6435.41359852
60% done : ELBO is 6502.28665116
70% done : ELBO is 6493.50219742
80% done : ELBO is 6494.55567465
90% done : ELBO is 6487.46927781
100% done : ELBO is 6391.83716934

Final model ELBO is 6492.17584989
EGARCHM(1,1)
======================================================= =================================================
Dependent Variable: AAPL Returns                        Method: BBVI
Start Date: 2006-01-05 00:00:00                         Unnormalized Log Posterior: 6505.2305
End Date: 2016-03-10 00:00:00                           AIC: -12998.461048
Number of observations: 2562                            BIC: -12963.3697871
=========================================================================================================
Latent Variable                          Median             Mean               95% Credibility Interval
======================================== ================== ================== ========================
Vol Constant                             -0.1223            -0.122             (-0.1925 | -0.0506)
p(1)                                     0.9834             0.9834             (0.9821 | 0.9847)
q(1)                                     0.1005             0.1005             (0.0932 | 0.1084)
v                                        5.964              5.9632             (5.5074 | 6.4545)
Returns Constant                         0.0039             0.004              (-0.0524 | 0.0602)
GARCH-M                                  -0.0339            -0.0338            (-0.1153 | 0.0464)
=========================================================================================================


We can plot the latent variable $\phi$:

model.plot_z([5],figsize=(15,5))


There is quite a lot of uncertainty for this latent variable, but the bulk of the mass is in negative territory. Let’s plot the fit with plot_fit:

model.plot_fit(figsize=(15,5))


And obtain predictions of future conditional volatility with predict:

model.predict(h=10)

AAPL Returns
2016-03-07 0.014771
2016-03-10 0.014901
2016-03-11 0.015029
2016-03-12 0.015157
2016-03-13 0.015283
2016-03-14 0.015408
2016-03-17 0.015532
2016-03-18 0.015656
2016-03-19 0.015778
2016-03-20 0.015899

Or plot the results with plot_predict:

model.plot_predict(h=10,figsize=(15,5))