Case study: fitting a Student t to daily returns

The data, and the best normal

The series is the one of case studies 5.3.5 and 5.4.6: 1,000 invented daily log returns. They were drawn with seed 50 from a Student t with 4 degrees of freedom and scaled to a volatility near 1%. They are illustrative, not market data.

Work in percent: $x_t = 100\,r_t$. Then the scale is near 1, like the degrees of freedom, which suits an optimiser and a single Hessian bump (lesson 4.3.1). The desk notes below show what changes in decimals.

Lesson 5.9.1 gives the normal's maximum likelihood estimates in closed form: the mean and the standard deviation dividing by $n$. Here they are $\hat{\mu} = -0.0061\%$ and $\hat{\sigma} = 1.0073\%$. That is case study 5.3.5's volatility: it was the maximum likelihood estimate all along.

At the peak the normal's log-likelihood is $-\frac{n}{2}\big(\ln(2\pi\hat{\sigma}^2) + 1\big) = -1426.25$. If the days really were normal, lesson 5.9.2 would give standard errors of $\frac{\hat{\sigma}}{\sqrt{n}} = 0.0319\%$ for the mean and $\frac{\hat{\sigma}}{\sqrt{2n}} = 0.0225\%$ for the volatility.

The Student t model

The model takes a standard t with $\nu$ degrees of freedom (lesson 5.4.4). It stretches it by a scale $s$ and shifts it to a location $m$: $X = m + sT$. Its density is the t's, changed as lesson 5.6.1 changes a density under $aX + b$:

Here $c_\nu$ is the t's normalising constant. SciPy computes $\ln f$ as t.logpdf(x, nu, m, s), so the log-likelihood is one sum.

Case study 5.4.6 fixed $\nu = 4$ and chose $s = \hat{\sigma}/\sqrt{2}$, so that the variance matched the sample's. Here all three parameters are free, and the likelihood picks them.

Fitting by maximum likelihood

There is no formula for the peak, so minimise the negative log-likelihood numerically, as in lesson 5.9.1. The scale and the degrees of freedom must stay positive, so search over $\ln s$ and $\ln\nu$. Start from the normal fit and $\nu = 5$.

scipy.stats.t.fit finds the same peak to the digits shown. The fitted t's volatility is $0.7386 \times \sqrt{4.353/2.353} = 1.0046\%$, close to the sample's 1.0073%.

The hand-chosen t of case study 5.4.6 sits only 1.03 below the peak. The data barely tell $\nu = 4$ from $\nu = 4.35$; the standard error below says how barely.

Why the location moved

The fitted location is $+0.0221\%$, while the sample mean is $-0.0061\%$. In the normal's log-likelihood a day $z$ scales from the centre costs $\frac{z^2}{2}$; in the t's it costs $\frac{\nu + 1}{2}\ln(1 + \frac{z^2}{\nu})$, which grows only like a log.

A day 5 scales out costs 12.5 under the normal but 5.11 under the fitted t. So the extreme days pull the t's location much less. Here that lands it nearer the recipe's $0.04\%$, though one sample proves nothing.

Standard errors from the Hessian

Lesson 5.9.2's recipe: the Hessian of the negative log-likelihood at the peak is the observed information $J$. The standard errors are the square roots of the diagonal of $J^{-1}$. Take the Hessian in the natural parameters $(m, s, \nu)$, by lesson 4.3.1's central differences.

The interval for $\nu$ contains 4, the value the recipe used and case study 5.4.6 chose. It stays well clear of 2, below which the fitted t would have no variance.

The location's standard error, 0.0274%, is below the normal's $\frac{\hat{\sigma}}{\sqrt{n}} = 0.0319\%$. When the data really are t, the t's location is the more precise estimate, because it does not lean on a few wild days.

The scale and the degrees of freedom trade off

The estimates of $s$ and $\nu$ have correlation $0.657$. A larger $\nu$ means thinner tails and a smaller factor $\frac{\nu}{\nu - 2}$, so a larger $s$ is needed for the same spread. As in lesson 5.9.2's line, the two move together.

So 1 over the square root of the diagonal entry for $\nu$ gives 0.45, a quarter below the correct 0.597. Always invert the whole matrix.

A profile likelihood for the degrees of freedom

The interval for $\nu$ rests on the quadratic copy of the log-likelihood. Lesson 5.9.2's likelihood-ratio interval needs no copy.

For each fixed $\nu$, maximise over $m$ and $s$. The result is the profile log-likelihood $\ell_p(\nu)$: how well the best t with that $\nu$ explains the data.

The 95% likelihood-ratio interval keeps every $\nu$ within 1.92 of the peak: 3.38 to 5.81. It is lopsided, reaching 1.45 above $\hat{\nu}$ and 0.97 below it. A change of $\nu$ matters less to the shape when $\nu$ is large, so the profile falls more slowly to the right.

The Hessian interval, 3.18 to 5.52, is symmetric by construction and misses that lean. For a tail parameter the likelihood-ratio interval is the one to report.

Is the t really better? Comparing log-likelihoods

The t gains 69.58 in log-likelihood over the normal, with one more parameter. The normal is the t with $\nu \to \infty$ (lesson 5.4.4), so the question is whether letting $\nu$ be finite helps by more than chance.

The likelihood-ratio test compares $2 \times 69.58 = 139.15$ with the 0.95-quantile of $\chi^2_1$, 3.841, the same number as the interval's 1.92 doubled. In the language of lesson 5.8.3, the p-value is about $4 \times 10^{-32}$: the normal is rejected.

One caution: $\nu = \infty$ is the edge of the range, where the $\chi^2_1$ rule is only rough. With a statistic 36 times the threshold, the verdict does not depend on it.

The fitted tails against the sample's

A fitted model is only as good as its tails, where risk limits bite. Case study 5.3.5 read value at risk (VaR) per $1,000,000 from the sample and from a normal. The fitted t gives a third answer, from all 1,000 days at once.

At 99% the fitted t gives about $25,900, close to the sample's $26,212, where the normal says $23,221. At 99.9% the normal is far too small. The t is below the sample there, but the sample's figure rests on the worst one or two days.

Counts of large days tell the same story. Case study 5.4.6 set them against the normal and its hand-chosen t; the fitted t joins them.

The two t's agree with each other and with the sample to within the noise of a few days. The normal expects 0.063 days beyond 4 volatilities, where the sample has 6.

import numpy as np


def _returns():
    rng = np.random.default_rng(50)
    return 0.0004 + 0.01 * rng.standard_t(4, size=1000) / np.sqrt(2)


def test_shared_returns():
    """The shared daily returns fit a t with nu about 4.353, scale about 0.007386 and location about 0.000221"""
    loc, scale, nu, loglik = fit_t(_returns())
    expect_close(nu, 4.352924, what="nu", rel=1e-3)
    expect_close(scale, 0.0073859132, what="scale", rel=1e-4)
    expect_close(loc, 0.00022087671, what="loc", rel=0, abs_tol=2e-7)


def test_loglik_decimals():
    """Their log-likelihood at the peak, in decimals, is about 3248.4946"""
    expect_close(fit_t(_returns())[3], 3248.4946329553, what="loglik", rel=0, abs_tol=1e-3)


def test_percent():
    """In percent the fit has location 0.0221, scale 0.7386 and log-likelihood -1356.6756"""
    loc, scale, nu, loglik = fit_t(100 * _returns())
    expect_close(loc, 0.0220876745, what="loc in percent", rel=0, abs_tol=2e-5)
    expect_close(scale, 0.7385913, what="scale in percent", rel=1e-4)
    expect_close(nu, 4.352924, what="nu in percent", rel=1e-3)
    expect_close(loglik, -1356.6755530328, what="loglik in percent", rel=0, abs_tol=1e-3)


def test_units():
    """Changing units scales the location and scale by 100 and leaves nu alone"""
    r = _returns()
    a, b = fit_t(r), fit_t(100 * r)
    expect_close(b[1], 100 * a[1], what="scale in percent / 100 * scale in decimals", rel=1e-4)
    expect_close(b[2], a[2], what="nu in percent against nu in decimals", rel=1e-3)
    expect_close(a[3] - b[3], 1000 * np.log(100), what="loglik difference, n ln 100", rel=0, abs_tol=2e-3)


def test_one_year():
    """The first 250 days alone give nu about 3.37"""
    expect_close(fit_t(_returns()[:250])[2], 3.3666, what="nu from 250 days", rel=1e-2)


def test_normal_data():
    """On 1,000 seeded normal returns nu comes out above 30: no fat tails"""
    z = np.random.default_rng(594).normal(0, 0.01, 1000)
    nu = fit_t(z)[2]
    assert nu > 30, f"expected nu above 30 for normal data, got {nu}"


def test_bad_returns():
    """Fewer than 20 returns, or a nan or inf among them, raises ValueError"""
    r = _returns()
    bad = [r[:19], np.append(r[:100], np.nan), np.append(r[:100], np.inf)]
    for x in bad:
        try:
            fit_t(x)
        except ValueError:
            continue
        raise AssertionError(f"expected ValueError for an input of {len(x)} values ending in {x[-1]}")


def test_var_99():
    """The 99% one-day VaR per $1,000,000 from the fitted t is about $25,912.55"""
    expect_close(t_var(_returns()), 25912.548, what="t_var(r)", rel=1e-4)


def test_var_levels():
    """At 99.9% it is about $47,465 and at 95% about $15,051"""
    r = _returns()
    expect_close(t_var(r, 0.999), 47465.04, what="t_var(r, 0.999)", rel=1e-4)
    expect_close(t_var(r, 0.95), 15050.61, what="t_var(r, 0.95)", rel=1e-4)


def test_var_value():
    """A $2,000,000 position has twice the VaR"""
    r = _returns()
    expect_close(t_var(r, 0.99, 2_000_000), 2 * 25912.548, what="t_var(r, 0.99, 2_000_000)", rel=1e-4)


def test_bad_level():
    """A level of 0, 1 or 1.5 raises ValueError"""
    r = _returns()
    for level in [0, 1, 1.5]:
        try:
            t_var(r, level)
        except ValueError:
            continue
        raise AssertionError(f"expected ValueError for level={level}")