Stochastic volatility: quasi-maximum likelihood

This notebook describes estimating the basic univariate stochastic volatility model by quasi-maximum likelihood methods, as in Ruiz (1994) or Harvey et al. (1994), using the Python package statsmodels.

The code defining the model is also available in a Python script sv.py.

%matplotlib inline
from __future__ import division

import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt

A canonical univariate stochastic volatility model (Kim et al., 1998) is:

We can transform the model for $y_t$ into a linear model in $\log y_t^2$, we have:

here, $E(\log \varepsilon_t^2) \approx -1.27$ and $Var(\log \varepsilon_t^2) = \pi^2 / 2$.

Quasi-likelihood - state space model

The quasi-likelihood approach (see e.g. Ruiz [1994] or Harvey et al. [1994]) replaces $\log \varepsilon_t^2$ with $E(\log \varepsilon_t^2) + \xi_t$ where $\xi_t \sim N(0, \pi^2 / 2)$. With this replacement, the model above becomes a linear Gaussian state space model and the Kalman filter can be used to compute the log-likelihoods using the prediction error decomposition.

Since $\log \varepsilon_t^2 - E(\log \varepsilon_t^2) \not \sim N(0, \pi^2 / 2)$, this is a quasi-maximum likelihood estimation routine that will only be valid asymptotically. However, Ruiz (1994) suggests that it may perform well nonetheless.

Example: daily exchange rates

As an example, we replicate the results of Harvey et al (1994), in which daily exchange rates between October 1, 1981 and June 28, 1985 are examined via the stochastic volatility model above.

Data

Harvey et al. (1994) consider the following exchange rates:

  • Pound / Dollar (USXUK)
  • Deutschmark / Dollar (USXGER)
  • Yen / Dollar (USXJPN)
  • Swiss-Franc / Dollar (USXUI)

Tom Doan of the RATS econometrics software package has made available the replication data at https://estima.com/forum/viewtopic.php?f=8&t=1595.

ex = pd.read_excel('data/xrates.xls')
dta = np.log(ex).diff().iloc[1:]
dta.plot(subplots=True, figsize=(15, 8));
*** No CODEPAGE record, no encoding_override: will use 'ascii'

png

Model

We can relatively easily construct the state space model in Statsmodels. This code is also given in the associated file sv.py.

from statsmodels.tsa.statespace.tools import (
    constrain_stationary_univariate, unconstrain_stationary_univariate)

# "Quasi-likelihood stochastic volatility"
class QLSV(sm.tsa.statespace.MLEModel):
    def __init__(self, endog):
        # Convert to log squares
        endog = np.log(endog**2)

        # Initialize the base model
        super(QLSV, self).__init__(endog, k_states=1, k_posdef=1,
                                   initialization='stationary')
        
        # Setup the observation covariance
        self['obs_intercept', 0, 0] = -1.27
        self['design', 0, 0] = 1
        self['obs_cov', 0, 0] = np.pi**2 / 2
        self['selection', 0, 0] = 1.
        
    @property
    def param_names(self):
        return ['phi', 'sigma2_eta', 'mu']
    
    @property
    def start_params(self):
        return np.r_[0.9, 1., 0.]
    
    def transform_params(self, params):
        return np.r_[constrain_stationary_univariate(params[:1]), params[1]**2, params[2:]]
    
    def untransform_params(self, params):
        return np.r_[unconstrain_stationary_univariate(params[:1]), params[1]**0.5, params[2:]]

    def update(self, params, **kwargs):
        super(QLSV, self).update(params, **kwargs)
        
        gamma = params[2] * (1 - params[0])
        self['state_intercept', 0, 0] = gamma
        self['transition', 0, 0] = params[0]
        self['state_cov', 0, 0] = params[1]

Notes

There are two important notes that will be needed to compare our results to those of Harvey et al. (1994):

  1. Recall that Harvey et al. (1994) estimate the intercept $\gamma$ rather than the mean $\mu$, so to compare our estimated mean to their estimated intercept, convert with $\gamma = \mu * (1 - \phi)$.
  2. Harvey et al. (1994) do not include the constant part of the loglikelihood $(-0.5 T \log(2 \pi))$ in their reported results, so to compare our likelihood value to theirs we must substract that constant.

Pound / Dollar

# Harvey et al. demean the data before estimating to avoid problems
# with taking logs of zeros.
endog = dta['USXUK'] - dta['USXUK'].mean()

# Fit the model
mod_usxuk = QLSV(endog)
res_usxuk = mod_usxuk.fit(cov_type='robust')

# Print the intercept that is comparable to Harvey et al (1994)
print('Gamma = %.4f' % (res_usxuk.params[2] * (1 - res_usxuk.params[0])))
# Print the loglikelihood that is comparable to Harvey et al (1994)
print('log L = %.2f' % (res_usxuk.llf - (-np.log(2 * np.pi) * mod_usxuk.nobs / 2)))

# Print the estimation summary
print(res_usxuk.summary())
Gamma = -0.0878
log L = -1212.82
                           Statespace Model Results                           
==============================================================================
Dep. Variable:                  USXUK   No. Observations:                  945
Model:                           QLSV   Log Likelihood               -2081.220
Date:                Sat, 14 Apr 2018   AIC                           4168.439
Time:                        21:57:57   BIC                           4182.993
Sample:                             0   HQIC                          4173.986
                                - 945                                         
Covariance Type:               robust                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
phi            0.9912      0.007    134.755      0.000       0.977       1.006
sigma2_eta     0.0070      0.005      1.356      0.175      -0.003       0.017
mu           -10.0038      0.273    -36.624      0.000     -10.539      -9.468
===================================================================================
Ljung-Box (Q):                       20.30   Jarque-Bera (JB):               970.27
Prob(Q):                              1.00   Prob(JB):                         0.00
Heteroskedasticity (H):               1.08   Skew:                            -1.33
Prob(H) (two-sided):                  0.48   Kurtosis:                         7.19
===================================================================================

Warnings:
[1] Quasi-maximum likelihood covariance matrix used for robustness to some misspecifications; calculated using the observed information matrix (complex-step) described in Harvey (1989).
fig, ax = plt.subplots(figsize=(13, 5))
ax.plot(np.abs(dta['USXUK']))
ax.plot(dta.index, np.exp(res_usxuk.smoothed_state[0] / 2), linewidth=3)
ax.set(title='Figure 1 of Harvey et al. (1994)');

png

Fit remaining models

We repeat the above exercise for the Deutschmark, Yen, and Swiss-Franc.

# Deutschmark / Dollar
mod_usxger = QLSV(dta['USXGER'] - dta['USXGER'].mean())
res_usxger = mod_usxger.fit(cov_type='robust')
print('Deutschmark / Dollar')
print('  phi        = %.4f' % res_usxger.params[0])
print('  sigma2_eta = %.4f' % res_usxger.params[1])
print('  gamma      = %.4f' % (res_usxger.params[2] * (1 - res_usxger.params[0])))
print('  log L      = %.2f' % (res_usxger.llf - (-np.log(2 * np.pi) * mod_usxger.nobs / 2)))
Deutschmark / Dollar
  phi        = 0.9649
  sigma2_eta = 0.0313
  gamma      = -0.3529
  log L      = -1232.26
# Yen / Dollar
mod_usxjpn = QLSV(dta['USXJPN'] - dta['USXJPN'].mean())
res_usxjpn = mod_usxjpn.fit()
print('Yen / Dollar')
print('  phi        = %.4f' % res_usxjpn.params[0])
print('  sigma2_eta = %.4f' % res_usxjpn.params[1])
print('  gamma      = %.4f' % (res_usxjpn.params[2] * (1 - res_usxjpn.params[0])))
print('  log L      = %.2f' % (res_usxjpn.llf - (-np.log(2 * np.pi) * mod_usxjpn.nobs / 2)))
Yen / Dollar
  phi        = 0.9947
  sigma2_eta = 0.0049
  gamma      = -0.0556
  log L      = -1272.64
# Swiss-Franc / Dollar
mod_usxsui = QLSV(dta['USXSUI'] - dta['USXSUI'].mean())
res_usxsui = mod_usxsui.fit()
print('Swiss-Franc / Dollar')
print('  phi        = %.4f' % res_usxsui.params[0])
print('  sigma2_eta = %.4f' % res_usxsui.params[1])
print('  gamma      = %.4f' % (res_usxsui.params[2] * (1 - res_usxsui.params[0])))
print('  log L      = %.2f' % (res_usxsui.llf - (-np.log(2 * np.pi) * mod_usxjpn.nobs / 2)))
Swiss-Franc / Dollar
  phi        = 0.9581
  sigma2_eta = 0.0459
  gamma      = -0.4180
  log L      = -1288.51