## Theory

The principle behind score-driven models is that the linear update $y_{t} - \theta_{t}$, that the Kalman filter relies upon, can be robustified by replacing it with the conditional score of a non-normal distribution. For this reason, any class of traditional state space model has a score-driven equivalent. This includes local level models, local linear trend models and dynamic regression models. Consider the dynamic regression problem in the score-driven framework:

$p\left(y_{t}\mid\theta_{t}\right)$

$\theta_{t} = \boldsymbol{x}_{t}^{'}\boldsymbol{\beta}_{t}$

$\boldsymbol{\beta}_{t} = \boldsymbol{\beta}_{t-1} + \boldsymbol{\eta}H_{t-1}^{-1}S_{t-1}$

Where $S_{t}$ is the score of the distribution. For a GAS regression with a Poisson distribution:

$p\left(y_{t}\mid\lambda_{t}\right)$

$\lambda_{t} = \exp\left(\theta_{t}\right) = \exp\left(\boldsymbol{x}_{t}^{'}\boldsymbol{\beta}_{t}\right)$

$\boldsymbol{\beta}_{t} = \boldsymbol{\beta}_{t-1} + \boldsymbol{\eta}\boldsymbol{x}_{t}\left(y_{t}-\theta_{t}\right)$

Here we have set $H_{t} = I_{t}$ so we are only using first-order information (for reasons of stability). Below we demonstrate this model type and others in PyFlux.

## PyFlux

### Poisson Data with a Local Level

We will use a custom dataset on the number of goals scored by soccer teams Nottingham Forest and Derby in their head-to-head matches from the beginning of their competitive history. We are interested to know whether these games have become more or less high scoring over time.

import numpy as np
import pandas as pd
import pyflux as pf
import matplotlib.pyplot as plt
%matplotlib inline

total_goals = pd.DataFrame(eastmidlandsderby['Forest'] + eastmidlandsderby['Derby'])
total_goals.columns = ['Total Goals']
plt.figure(figsize=(15,5))
plt.title("Total Goals in the East Midlands Derby")
plt.xlabel("Games Played");
plt.ylabel("Total Goals");
plt.plot(total_goals);


We can fit a GAS Local Level model using GASLLEV. We will choose a Poisson family.

model = pf.GASLLEV(data=total_goals,family=pf.GASPoisson())
x = model.fit()
x.summary()

Poisson GAS LLM
======================================================= =================================================
Dependent Variable: Total Goals                         Method: MLE
Start Date: 1                                           Log Likelihood: -198.7993
End Date: 96                                            AIC: 399.5985
Number of observations: 96                              BIC: 402.1629
=========================================================================================================
Latent Variable                          Estimate   Std Error  z        P>|z|    95% C.I.
======================================== ========== ========== ======== ======== ========================
SC(1)                                    0.1929     0.0468     4.1252   0.0      (0.1013 | 0.2846)
=========================================================================================================


We can plot the local level using plot_fit:

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


And predict forward using plot_predict:

model.plot_predict(figsize=(15,5),past_values=50)


### Local Linear Trend Model for US GDP

We will construct a local linear trend model for US GDP using a skew t-distribution. Here is the data:

growthdata = pd.read_csv('http://www.pyflux.com/notebooks/GDPC1.csv')
USgrowth = pd.DataFrame(np.log(growthdata['VALUE']))
USgrowth.index = pd.to_datetime(growthdata['DATE'])
USgrowth.columns = ['Logged US Real GDP']
plt.figure(figsize=(15,5))
plt.plot(USgrowth.index, USgrowth)
plt.ylabel('Real GDP')
plt.title('US Logged Real GDP');


We will fit with BBVI – this tends to be more robust than MLE for non-linear problems (MLE with L-BFGS-B relies on good starting values and is less effective at escaping local optima).

model2 = pf.GASLLT(data=USgrowth-np.mean(USgrowth),family=pf.GASSkewt())
x = model2.fit('BBVI',iterations=10000)

10% done : ELBO is -1716.00774587
20% done : ELBO is -446.28785162
30% done : ELBO is 99.3711159141
40% done : ELBO is 151.136245019
50% done : ELBO is 158.04058071
60% done : ELBO is 161.975861311
70% done : ELBO is 173.044112882
80% done : ELBO is 193.694813145
90% done : ELBO is 211.479364626
100% done : ELBO is 220.04243758

Final model ELBO is 215.74246198

x.summary()

Skewt GAS LLT
======================================================= =================================================
Dependent Variable: Logged US Real GDP                  Method: BBVI
Start Date: 1947-04-01 00:00:00                         Unnormalized Log Posterior: 230.1653
End Date: 2015-10-01 00:00:00                           AIC: -450.330655617
Number of observations: 275                             BIC: -432.246800129
=========================================================================================================
Latent Variable                          Median             Mean               95% Credibility Interval
======================================== ================== ================== ==========================
Level Scale                              0.0445             0.0445             (0.0368 | 0.0522)
Trend Scale                              0.0185             0.0185             (0.0102 | 0.0267)
Skewness                                 1.2089             1.2089             (1.1227 | 1.3002)
Skewt Scale                              0.1438             0.1438             (0.1362 | 0.1519)
v                                        5.886              5.8853             (5.3012 | 6.5345)
=========================================================================================================

And we can plot the latent variables with plot_z():
model2.plot_z([0,1,3])

model2.plot_z([2,4])


We can access the states through the saved model.fit() object. We can look at the local linear trend which provides an idea of the underlying growth potential of the economy:

plt.figure(figsize=(15,5))
plt.title("Local Trend for US GDP")
plt.ylabel("Trend")
plt.plot(USgrowth.index[21:],x.states[1][20:]);


We can predict forward using predict. Let’s predict the average growth rate for the next 4 quarters:

print("Average growth rate for this period is")
print(str(round(100*np.mean(np.exp(np.diff(model2.predict(h=4)['Logged US Real GDP'].values)) - 1),3)) + "%")

Average growth rate for this period is
0.502%


### Dynamic t regression

We start with a problem posed by another model. In portfolio construction, we are often interested in extracting the $\beta$ for a stock. If we assume returns are normally distributed then we can use dynamic linear regression from the Kalman filter algorithm. Let’s plot some example financial data, and plot the fit of a model on this data using this filter.

from pandas.io.data import DataReader
from datetime import datetime

a = DataReader('AMZN',  'yahoo', datetime(2012,1,1), datetime(2016,6,1))
a_returns = pd.DataFrame(np.diff(np.log(a['Adj Close'].values)))
a_returns.index = a.index.values[1:a.index.values.shape[0]]
a_returns.columns = ["Amazon Returns"]

spy = DataReader('SPY',  'yahoo', datetime(2012,1,1), datetime(2016,6,1))
spy_returns = pd.DataFrame(np.diff(np.log(spy['Adj Close'].values)))
spy_returns.index = spy.index.values[1:spy.index.values.shape[0]]
spy_returns.columns = ['S&P500 Returns']

one_mon = DataReader('DGS1MO', 'fred',datetime(2012,1,1), datetime(2016,6,1))
one_day = np.log(1+one_mon)/365

returns = pd.concat([one_day,a_returns,spy_returns],axis=1).dropna()
excess_m = returns["Amazon Returns"].values - returns['DGS1MO'].values
excess_spy = returns["S&P500 Returns"].values - returns['DGS1MO'].values
final_returns = pd.DataFrame(np.transpose([excess_m,excess_spy, returns['DGS1MO'].values]))
final_returns.columns=["Amazon","SP500","Risk-free rate"]
final_returns.index = returns.index

plt.figure(figsize=(15,5))
plt.title("Excess Returns")
x = plt.plot(final_returns);
plt.legend(iter(x), final_returns.columns);


model = pf.DynReg('Amazon ~ SP500',data=final_returns)
model.fit()
model.plot_fit(figsize=(15,15),series_type='Filtered')


Because we are assuming normality, the Kalman filter overreact to tail events – see for example the third graph around February 2015. Our beta estimate for Amazon drops from around $1$ to $0.5$ in response to this tail event. We want a dynamic estimate of beta that is more robust to tail events. We can use a GASReg model with a t-distribution:

model3 = pf.GASReg('Amazon ~ SP500',data=final_returns,family=pf.GASt())
x = model3.fit()
x.summary()

t GAS Regression
======================================================= =================================================
Dependent Variable: Amazon                              Method: MLE
Start Date: 2012-01-04 00:00:00                         Log Likelihood: 3158.435
End Date: 2016-06-01 00:00:00                           AIC: -6308.87
Number of observations: 1101                            BIC: -6288.8541
=========================================================================================================
Latent Variable                          Estimate   Std Error  z        P>|z|    95% C.I.
======================================== ========== ========== ======== ======== ========================
Scale 1                                  0.0
Scale SP500                              0.0474
t Scale                                  0.0095
v                                        2.8518
=========================================================================================================


And we can plot the fit:

model3.plot_fit(intervals=False,figsize=(15,15))


Looking at the final graph, we do not have the discrete jumps in our beta estimate corresponding with the February and August shocks (like with the Kalman filter). This is because of the ‘trimming’ property of the score term of the t-distribution, which reduces the influence of outliers.