The decay a sample¶

A day's decay spread over its samples, lambda = d^(1/n): the decay per sample whose product over n samples is the day's, so an EWMA stepped a sample at a time weights each day as the daily one would.

Both dimensionless: a decay in (0, 1) stated a day (RiskMetrics' 0.94) and a count of samples a day (above nothing). A decay stated over any other period spreads the same way: d over n of its parts. Source: RiskMetrics Technical Document, 4th ed. (J.P. Morgan/Reuters, 1996), 5.3, the decay factor a day; the relation lambda^n = d. Checked against a reference in numpy, which has its own tests.

Source

RiskMetrics Technical Document, 4th ed. (J.P. Morgan/Reuters, 1996), 5.3, the decay factor a day; the relation lambda^n = d.

Proof

reference, invariant, limit.

sample_decay.trait

namespace lusine.finance;

float sample_decay(float daily, float samples)
{
  return pow(daily, 1.0 / samples);
}

A day's decay spread over its samples, d^(1/n), so an EWMA stepped a sample at a time weights each day as the daily estimator would (RiskMetrics, 5.3).

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 lusine import trait_run
F32 = np.float32
U = float(np.finfo(F32).eps) / 2        # binary32's unit roundoff, 2^-24

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 lusine import trait_run F32 = np.float32 U = float(np.finfo(F32).eps) / 2 # binary32's unit roundoff, 2^-24

The kernel¶

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

In [2]:

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

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

// lusine.finance.sample_decay: a day's decay spread over its samples, lambda = d^(1/n): the decay per sample whose
// product over n samples is the day's, so an EWMA stepped a sample at a time weights each day as the daily one would.
// Both dimensionless: a decay in (0, 1) stated a day (RiskMetrics' 0.94) and a count of samples a day (above
// nothing). A decay stated over any other period spreads the same way: d over n of its parts.
// Source: RiskMetrics Technical Document, 4th ed. (J.P. Morgan/Reuters, 1996), 5.3, the decay factor a day; the
// relation lambda^n = d.
// Checked against a reference in numpy, which has its own tests.

namespace lusine.finance;

float sample_decay(float daily, float samples)
{
  return pow(daily, 1.0 / samples);
}

The program¶

A lattice of schedules, each a daily decay and the samples a day, both the host's; main() binds the kernel over them.

In [3]:

DECAY = """namespace decay;

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

struct schedule { host float daily; host float samples; float lambda; };
lattice<schedule, 1> schedules(spacing = 1 m, boundary = wall);

void main() { schedules.lambda = lusine.finance.sample_decay(schedules.daily, schedules.samples); }
"""
print(DECAY)


def per_sample(daily, samples):
    """Each schedule's decay a sample."""
    n = max(np.size(daily), np.size(samples))
    run = trait_run.Program(['finance/sample_decay.trait'], DECAY, {'decay.schedules': [n]})
    run.store('decay.schedules', 'daily', np.broadcast_to(np.asarray(daily, dtype=float), (n,)))
    run.store('decay.schedules', 'samples', np.broadcast_to(np.asarray(samples, dtype=float), (n,)))
    return run.steps('decay.main', 1, [('decay.schedules', 'lambda')])[0][0]

DECAY = """namespace decay; lusine.version version = { { number = 1 } }; struct schedule { host float daily; host float samples; float lambda; }; lattice<schedule, 1> schedules(spacing = 1 m, boundary = wall); void main() { schedules.lambda = lusine.finance.sample_decay(schedules.daily, schedules.samples); } """ print(DECAY) def per_sample(daily, samples): """Each schedule's decay a sample.""" n = max(np.size(daily), np.size(samples)) run = trait_run.Program(['finance/sample_decay.trait'], DECAY, {'decay.schedules': [n]}) run.store('decay.schedules', 'daily', np.broadcast_to(np.asarray(daily, dtype=float), (n,))) run.store('decay.schedules', 'samples', np.broadcast_to(np.asarray(samples, dtype=float), (n,))) return run.steps('decay.main', 1, [('decay.schedules', 'lambda')])[0][0]

namespace decay;

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

struct schedule { host float daily; host float samples; float lambda; };
lattice<schedule, 1> schedules(spacing = 1 m, boundary = wall);

void main() { schedules.lambda = lusine.finance.sample_decay(schedules.daily, schedules.samples); }

Reference: RiskMetrics' 0.94 over a day's samples¶

RiskMetrics' daily decay of 0.94 (Technical Document, 4th ed., 5.3) spread over 24 hourly samples, over 390 minutes, and over 6.5 hours of 10-second samples: each within the library's pow of the float64 reference's root, and each compounded over its samples the day's 0.94 again.

In [4]:

from lusine.kernels.finance.sample_decay import sample_decay as reference
n = np.array([24.0, 390.0, 2340.0])
got = per_sample(0.94, n)
want = reference(F32(0.94).astype(float), n)
assert np.all(np.abs(got / want - 1) <= 8 * U), got / want - 1
assert np.all(np.abs(got.astype(float) ** n - 0.94) < 1e-4), got.astype(float) ** n
print(dict(zip(n.astype(int), np.round(got.astype(float), 6))))

from lusine.kernels.finance.sample_decay import sample_decay as reference n = np.array([24.0, 390.0, 2340.0]) got = per_sample(0.94, n) want = reference(F32(0.94).astype(float), n) assert np.all(np.abs(got / want - 1) <= 8 * U), got / want - 1 assert np.all(np.abs(got.astype(float) ** n - 0.94) < 1e-4), got.astype(float) ** n print(dict(zip(n.astype(int), np.round(got.astype(float), 6))))

{np.int64(24): np.float64(0.997425), np.int64(390): np.float64(0.999841), np.int64(2340): np.float64(0.999974)}

Invariant: compounded over the day, the day's¶

Two hundred daily decays from 0.5 to 0.999 over 1 to 1,000 samples: each decay a sample, compounded over its samples in float64, is the day's within the pow's few ulps carried through the power; and a sample's decay is nearer one than its day's.

In [5]:

rng = np.random.default_rng(419)
d = rng.uniform(0.5, 0.999, 200)
m = rng.integers(1, 1000, 200).astype(float)
lam = per_sample(d, m).astype(float)
back = lam ** m
assert np.all(np.abs(back / F32(d).astype(float) - 1) <= 8 * U * m), np.max(np.abs(back / F32(d).astype(float) - 1) / m)
assert np.all(lam[m > 1] > F32(d[m > 1]))
print('worst relative error of the compounded decay, per sample compounded:', float(np.max(np.abs(back / F32(d).astype(float) - 1) / m)))
fig, ax = plt.subplots(figsize=(7, 3))
nn = np.array([1, 2, 4, 8, 24, 96, 390, 1440], dtype=float)
for dd in (0.9, 0.94, 0.97):
    ax.plot(nn, per_sample(dd, nn), 'o-', label=f'{dd} a day')
ax.set_xscale('log')
ax.set_xlabel('samples a day')
ax.set_ylabel('decay a sample')
ax.legend()
plt.show()

rng = np.random.default_rng(419) d = rng.uniform(0.5, 0.999, 200) m = rng.integers(1, 1000, 200).astype(float) lam = per_sample(d, m).astype(float) back = lam ** m assert np.all(np.abs(back / F32(d).astype(float) - 1) <= 8 * U * m), np.max(np.abs(back / F32(d).astype(float) - 1) / m) assert np.all(lam[m > 1] > F32(d[m > 1])) print('worst relative error of the compounded decay, per sample compounded:', float(np.max(np.abs(back / F32(d).astype(float) - 1) / m))) fig, ax = plt.subplots(figsize=(7, 3)) nn = np.array([1, 2, 4, 8, 24, 96, 390, 1440], dtype=float) for dd in (0.9, 0.94, 0.97): ax.plot(nn, per_sample(dd, nn), 'o-', label=f'{dd} a day') ax.set_xscale('log') ax.set_xlabel('samples a day') ax.set_ylabel('decay a sample') ax.legend() plt.show()

worst relative error of the compounded decay, per sample compounded: 3.0403930176922025e-08

Limit: one sample a day, and a decay of one¶

One sample a day keeps the daily decay, and a decay of one spreads to one: binary32's pow(x, 1) is x and pow(1, y) is 1 (lusine.stochastic.step_chance leans on the first).

In [6]:

out = per_sample([0.94, 1.0, 0.5], [1.0, 390.0, 1.0])
assert out[0] == F32(0.94) and out[1] == 1.0 and out[2] == 0.5, out
print(out)

out = per_sample([0.94, 1.0, 0.5], [1.0, 390.0, 1.0]) assert out[0] == F32(0.94) and out[1] == 1.0 and out[2] == 0.5, out print(out)

[0.94 1.   0.5 ]