The EWMA variance¶

The exponentially weighted variance one sample on, v' = lambda v + (1 - lambda) r^2: the variance of the log return over one sample's interval, its last value weighted by the decay and the new sample's square by the rest (the mean return taken as zero).

All three are dimensionless: a variance per sample, a log return, a decay in (0, 1). The decay is the game's (RiskMetrics' 0.94 is a day's), and so is spreading it over a day's samples and annualising the variance. Source: RiskMetrics Technical Document, 4th ed. (J.P. Morgan/Reuters, 1996), 5.2.1, Eq. [5.3] and Table 5.3. Checked against a reference in numpy, which has its own tests.

Source

RiskMetrics Technical Document, 4th ed. (J.P. Morgan/Reuters, 1996), 5.2.1, Eq. [5.3] and Table 5.3.

Proof

reference, invariant, limit.

ewma_variance.trait

namespace lusine.finance;

float ewma_variance(float variance, float sample, float decay)
{
  return decay * variance + (1.0 - decay) * sample * sample;
}

RiskMetrics' exponentially weighted variance one sample on, lambda v + (1 - lambda) r^2.

In [1]:

import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
    root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import ndtr
from scipy.stats import norm
from lusine import trait_run
EPS = float(np.finfo(np.float32).eps)   # Trait's float is binary32
YEAR = 31557600.0                        # lusine's year row, the Julian year in seconds: the host speaks base units

import os, sys root = os.getcwd() while not os.path.isfile(os.path.join(root, 'docsite.json')): root = os.path.dirname(root) sys.path.insert(0, os.path.join(root, 'tools')) import numpy as np import matplotlib.pyplot as plt from scipy.special import ndtr from scipy.stats import norm from lusine import trait_run EPS = float(np.finfo(np.float32).eps) # Trait's float is binary32 YEAR = 31557600.0 # lusine's year row, the Julian year in seconds: the host speaks base units

The kernel¶

The law as the library states it, and as this page runs it.

In [2]:

print(trait_run.source('finance/ewma_variance.trait'))

print(trait_run.source('finance/ewma_variance.trait'))

// lusine.finance.ewma_variance: the exponentially weighted variance one sample on, v' = lambda v + (1 - lambda) r^2:
// the variance of the log return over one sample's interval, its last value weighted by the decay and the new sample's
// square by the rest (the mean return taken as zero). All three are dimensionless: a variance per sample, a log return,
// a decay in (0, 1). The decay is the game's (RiskMetrics' 0.94 is a day's), and so is spreading it over a day's
// samples and annualising the variance.
// Source: RiskMetrics Technical Document, 4th ed. (J.P. Morgan/Reuters, 1996), 5.2.1, Eq. [5.3] and Table 5.3.
// Checked against a reference in numpy, which has its own tests.

namespace lusine.finance;

float ewma_variance(float variance, float sample, float decay)
{
  return decay * variance + (1.0 - decay) * sample * sample;
}

The program¶

One series a lattice item: each step the host writes the new sample, and main() folds it into the item's variance by the kernel. The decay is the host's too; RiskMetrics' 0.94 is a day's, and spreading it over a day's samples is the game's.

In [3]:

EWMA = """namespace vol;

lusine.version version = { { number = 1 } };

struct history { float variance; host float sample; host float decay; };
lattice<history, 1> series(spacing = 1 m, boundary = wall);

void main() { series.variance = lusine.finance.ewma_variance(series.variance, series.sample, series.decay); }
"""
print(EWMA)


def path(start, samples, decay):
    """The variance after each sample: start and decay one a series, samples [step, series]."""
    samples = np.atleast_2d(np.asarray(samples, dtype=float).T).T
    run = trait_run.Program(['finance/ewma_variance.trait'], EWMA, {'vol.series': [samples.shape[1]]})
    run.store('vol.series', 'variance', np.broadcast_to(start, samples.shape[1:]))
    run.store('vol.series', 'decay', np.broadcast_to(decay, samples.shape[1:]))
    out = []
    for row in samples:
        run.store('vol.series', 'sample', row)
        out.append(run.steps('vol.main', 1, [('vol.series', 'variance')])[0][0])
    return np.array(out)

EWMA = """namespace vol; lusine.version version = { { number = 1 } }; struct history { float variance; host float sample; host float decay; }; lattice<history, 1> series(spacing = 1 m, boundary = wall); void main() { series.variance = lusine.finance.ewma_variance(series.variance, series.sample, series.decay); } """ print(EWMA) def path(start, samples, decay): """The variance after each sample: start and decay one a series, samples [step, series].""" samples = np.atleast_2d(np.asarray(samples, dtype=float).T).T run = trait_run.Program(['finance/ewma_variance.trait'], EWMA, {'vol.series': [samples.shape[1]]}) run.store('vol.series', 'variance', np.broadcast_to(start, samples.shape[1:])) run.store('vol.series', 'decay', np.broadcast_to(decay, samples.shape[1:])) out = [] for row in samples: run.store('vol.series', 'sample', row) out.append(run.steps('vol.main', 1, [('vol.series', 'variance')])[0][0]) return np.array(out)

namespace vol;

lusine.version version = { { number = 1 } };

struct history { float variance; host float sample; host float decay; };
lattice<history, 1> series(spacing = 1 m, boundary = wall);

void main() { series.variance = lusine.finance.ewma_variance(series.variance, series.sample, series.decay); }

Reference: RiskMetrics' Table 5.3¶

USD/DEM's daily returns, 28 March to 24 April 1996, at a decay of 0.94, starting from the first return squared (RiskMetrics Technical Document, 4th ed., 1996, 5.2.1, Table 5.3): each day's variance, in percent squared, is the table's to 0.001 (the returns are printed to 0.001%, so the table's own carry more digits than these), and the forecast for 25 April is its 0.473%.

In [4]:

r = np.array([0.633, 0.115, -0.459, 0.093, 0.176, -0.087, -0.142, 0.324, -0.943, -0.528, -0.107,
              -0.159, -0.445, 0.053, 0.152, -0.318, 0.424, -0.708, -0.105, -0.257]) / 100
printed = np.array([0.401, 0.378, 0.368, 0.346, 0.327, 0.308, 0.291, 0.280, 0.316, 0.314, 0.296,
                    0.280, 0.275, 0.258, 0.244, 0.236, 0.232, 0.248, 0.234, 0.224])
v = np.concatenate([[r[0] ** 2], path(r[0] ** 2, r[1:], 0.94)[:, 0]]) * 1e4
assert np.max(np.abs(v - printed)) < 0.001, v
assert abs(np.sqrt(v[-1]) - 0.473) < 0.0005
print(np.round(v, 4))
fig, ax = plt.subplots(figsize=(7, 3.5))
ax.plot(v, label='the kernel')
ax.plot(printed, 'o', label='Table 5.3')
ax.set_xlabel('day')
ax.set_ylabel('variance (%^2)')
ax.legend()
plt.show()

r = np.array([0.633, 0.115, -0.459, 0.093, 0.176, -0.087, -0.142, 0.324, -0.943, -0.528, -0.107, -0.159, -0.445, 0.053, 0.152, -0.318, 0.424, -0.708, -0.105, -0.257]) / 100 printed = np.array([0.401, 0.378, 0.368, 0.346, 0.327, 0.308, 0.291, 0.280, 0.316, 0.314, 0.296, 0.280, 0.275, 0.258, 0.244, 0.236, 0.232, 0.248, 0.234, 0.224]) v = np.concatenate([[r[0] ** 2], path(r[0] ** 2, r[1:], 0.94)[:, 0]]) * 1e4 assert np.max(np.abs(v - printed)) < 0.001, v assert abs(np.sqrt(v[-1]) - 0.473) < 0.0005 print(np.round(v, 4)) fig, ax = plt.subplots(figsize=(7, 3.5)) ax.plot(v, label='the kernel') ax.plot(printed, 'o', label='Table 5.3') ax.set_xlabel('day') ax.set_ylabel('variance (%^2)') ax.legend() plt.show()

[0.4007 0.3774 0.3674 0.3459 0.327  0.3078 0.2906 0.2794 0.316  0.3138
 0.2957 0.2794 0.2746 0.2582 0.2441 0.2356 0.2322 0.2484 0.2341 0.224 ]

Invariant: the exponentially weighted sum¶

Run over 300 samples from nothing, the recursion is (1 - lambda) sum lambda^i r_(t-i)^2, the estimator RiskMetrics derives it from (Eq. [5.5]), within binary32's rounding over the run; four series at once, four decays.

In [5]:

rng = np.random.default_rng(391)
samples = rng.normal(0.0, 0.01, (300, 4))
decays = np.array([0.5, 0.9, 0.94, 0.97])
got = path(0.0, samples, np.broadcast_to(decays, (4,)))[-1]
lam = np.float32(decays).astype(float)
w = lam[None, :] ** np.arange(299)[::-1, None]
want = (1 - lam) * np.sum(w * np.float32(samples[1:]).astype(float) ** 2, axis=0) + (1 - lam) * lam ** 299 * np.float32(samples[0]).astype(float) ** 2
rel = np.abs(got / want - 1)
assert np.max(rel) <= 1e-5, rel
print('relative', rel)

rng = np.random.default_rng(391) samples = rng.normal(0.0, 0.01, (300, 4)) decays = np.array([0.5, 0.9, 0.94, 0.97]) got = path(0.0, samples, np.broadcast_to(decays, (4,)))[-1] lam = np.float32(decays).astype(float) w = lam[None, :] ** np.arange(299)[::-1, None] want = (1 - lam) * np.sum(w * np.float32(samples[1:]).astype(float) ** 2, axis=0) + (1 - lam) * lam ** 299 * np.float32(samples[0]).astype(float) ** 2 rel = np.abs(got / want - 1) assert np.max(rel) <= 1e-5, rel print('relative', rel)

relative [5.62363393e-08 6.53103664e-08 1.39158647e-07 8.76619708e-08]

Limit: a decay of one, a quiet market¶

A decay of one keeps the variance bit for bit; a run of samples of nothing lets it fall by the decay each sample.

In [6]:

kept = path(4e-4, np.full(10, 0.05), 1.0)[:, 0]
assert np.all(kept == np.float32(4e-4)), kept
quiet = path(4e-4, np.zeros(10), 0.94)[:, 0]
assert np.all(np.abs(quiet / (4e-4 * 0.94 ** np.arange(1, 11)) - 1) <= 16 * EPS), quiet
print(quiet)

kept = path(4e-4, np.full(10, 0.05), 1.0)[:, 0] assert np.all(kept == np.float32(4e-4)), kept quiet = path(4e-4, np.zeros(10), 0.94)[:, 0] assert np.all(np.abs(quiet / (4e-4 * 0.94 ** np.arange(1, 11)) - 1) <= 16 * EPS), quiet print(quiet)

[0.000376   0.00035344 0.00033223 0.0003123  0.00029356 0.00027595
 0.00025939 0.00024383 0.0002292  0.00021545]