Skip to content

The EWMA variance, live

RiskMetrics' exponentially weighted variance one sample on, v′ = λ v + (1 − λ) r²: 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 (lusine.finance.ewma_variance). Below, RiskMetrics' twenty daily returns of USD/DEM are folded on your GPU as you read, a day a step at a decay of 0.94, from the first day's return squared; and four series of returns drawn uniform with a spread of 1% are folded from nothing at decays of 0.5, 0.9, 0.94 and 0.97 for three hundred samples, each sample's square kept beside them with its weight, (1 − λ) when it is taken and λ times that each sample after.

This page runs the program on your GPU through WebGPU, so it needs a browser with WebGPU (a recent Chrome, Edge, Safari or Firefox). Without it the program does not run and nothing is drawn below.

What it shows

  • The strip above the controls is green while every check holds, and turns red for good once one fails.
  • RiskMetrics' Table 5.3: the chart draws each day's variance, in percent squared, the kernel's and the table's, one over the other. Below it, the day and the two.
  • The weighted sum: the chart draws the series at a decay of 0.94, its variance by the kernel and the weighted sum of its history, (1 − λ) Σ λⁱ r², in percent squared, one over the other. Below it, the sample.
  • The checks: RiskMetrics' Table 5.3, each day's variance the table's to 0.001 percent squared, the returns being printed to 0.001% (RiskMetrics Technical Document, 4th ed., 1996, 5.2.1); the forecast, 0.473%, the last day's variance's root RiskMetrics' forecast for 25 April within 0.0005; the weighted sum, every series' variance the weighted sum of its samples' squares, the estimator RiskMetrics derives the recursion from (Eq. [5.5]), within a part in 10⁵ over the run; a decay of one, a quiet market, a decay of one keeping a variance of 0.0004 bit for bit over ten samples of 5%, and ten samples of nothing letting it fall by 0.94 each, within 16 binary32 epsilons (2⁻²³) of 0.0004 × 0.94ᵏ in float64. Each is held at every step.

Every number here is computed by the program on the GPU; the page only reads the values back and draws them. The seed, the series' spread, the samples a run takes and RiskMetrics' decay are the program's settings, which the host may write; this page offers no control to change them yet.

The program

RiskMetrics' days, the window, the series, their history and the checks, the panels the page draws, and the one procedure, main, run each step: a day folded into the window, a sample into each series and its history, the weights moved on, and the checks.

// The EWMA variance's live page: RiskMetrics' twenty days of USD/DEM folded a day a step by the library's kernel
// against the variances its table prints; four series of drawn returns folded at four decays beside the exponentially
// weighted sum they stand for, kept sample by sample; and a decay of one and a quiet market, with the page's checks
// computed here on the GPU and shown by the panels the program declares.

namespace ewma_variance;

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

// The reader's settings: the seed the series are drawn from, their spread a sample, the samples a run takes before it
// starts again, and RiskMetrics' daily decay.
host struct { integer seed = 391; float spread = 0.01; integer length = 300; float decay = 0.94; } settings;

// RiskMetrics' Table 5.3: USD/DEM's daily returns from 28 March to 24 April 1996, and each day's variance in percent
// squared as the table prints it.
struct day { integer number; float sample; float printed; };
host table<day> days(capacity = 20) = {
  { number = 0, sample = 0.00633, printed = 0.401 }, { number = 1, sample = 0.00115, printed = 0.378 },
  { number = 2, sample = -0.00459, printed = 0.368 }, { number = 3, sample = 0.00093, printed = 0.346 },
  { number = 4, sample = 0.00176, printed = 0.327 }, { number = 5, sample = -0.00087, printed = 0.308 },
  { number = 6, sample = -0.00142, printed = 0.291 }, { number = 7, sample = 0.00324, printed = 0.28 },
  { number = 8, sample = -0.00943, printed = 0.316 }, { number = 9, sample = -0.00528, printed = 0.314 },
  { number = 10, sample = -0.00107, printed = 0.296 }, { number = 11, sample = -0.00159, printed = 0.28 },
  { number = 12, sample = -0.00445, printed = 0.275 }, { number = 13, sample = 0.00053, printed = 0.258 },
  { number = 14, sample = 0.00152, printed = 0.244 }, { number = 15, sample = -0.00318, printed = 0.236 },
  { number = 16, sample = 0.00424, printed = 0.232 }, { number = 17, sample = -0.00708, printed = 0.248 },
  { number = 18, sample = -0.00105, printed = 0.234 }, { number = 19, sample = -0.00257, printed = 0.224 },
};

// For the limits: a variance of 0.0004, a sample of 5%, and a decay of one.
host struct { float start = 0.0004; float sample = 0.05; float one = 1.0; } limit;

// The table's window: the day, from 0 to 19 and round again; the variance, from the first day's return squared and
// each day after it by the kernel at RiskMetrics' decay; and in percent squared, beside the table's.
struct { integer day = -1; float variance; float percent; float printed; } window;

// Four series, a sample each step uniform with the spread, each folded at its own decay from nothing; and the
// exponentially weighted sum of the samples' squares so far, (1 - lambda) sum lambda^i r^2, from the history.
struct series { integer index; float decay; float sample; float variance; float weighted; bool agrees; };
table<series> serieses(capacity = 4) = {
  { index = 0, decay = 0.5 }, { index = 1, decay = 0.9 }, { index = 2, decay = 0.94 }, { index = 3, decay = 0.97 },
};

// The history: each sample's square kept with its series and its weight, (1 - lambda) when it is taken and lambda times
// that each sample after.
struct taken { integer series; float decay; float square; float weight; };
table<taken> history(capacity = 1200);

// The run: the sample within it, and, for the chart, the third series' variance and weighted sum in percent squared.
struct { integer step = -1; float variance; float weighted; } run;

// The checks: each day's variance RiskMetrics' printed one to 0.001 percent squared, and the forecast for 25 April,
// the last variance's root, its 0.473%; each series' variance the weighted sum of its history within a part in 10^5; a
// decay of one keeping the variance bit for bit over ten samples; and a quiet market's variance falling by the decay
// each sample, within 16 binary32 epsilons of 4e-4 times 0.94^k in float64.
struct { bool table = true; bool forecast = true; bool weighted = true; bool limits = true; } checks;

// A strip of lamps the page draws, lit green while every check holds and red once one fails.
struct lamp { vector<float, 4> colour; };
lattice<lamp, 2> lamps(spacing = 1 m, boundary = wall);

vector<float, 4> lamp_colour(bool holds)
{
  return holds ? vector<float, 4>(0.15, 0.6, 0.3, 1.0) : vector<float, 4>(0.75, 0.15, 0.15, 1.0);
}

float eps() { return 1.0 / 8388608.0; }

float on_day(table<integer> numbers, table<float> values, integer day)
{
  let k = day;
  return sum(values) where numbers == k;
}

float drawn(integer seed, integer step, integer index, float spread)
{
  let u = lusine.error.uniform(drawn = lusine.error.draw(seed = bits(seed), step = bits(step), key = fnv1a("series"),
                                                         index = bits(index)));
  return sqrt(3.0) * spread * (2.0 * u - 1.0);
}

taken recorded(integer index, float decay, float sample)
{
  return taken { series = index, decay = decay, square = sample * sample, weight = 1.0 - decay };
}

bool restarting(integer step) { return step == 0; }

float weighting(integer index, table<integer> series, table<float> square, table<float> weight)
{
  let k = index;
  return sum(weight * square) where series == k;
}

integer disagreeing(table<bool> agrees) { return count(agrees) where !agrees; }

bool kept(float start, float sample, float one)
{
  let v1 = lusine.finance.ewma_variance(variance = start, sample, decay = one);
  let v2 = lusine.finance.ewma_variance(variance = v1, sample, decay = one);
  let v3 = lusine.finance.ewma_variance(variance = v2, sample, decay = one);
  let v4 = lusine.finance.ewma_variance(variance = v3, sample, decay = one);
  let v5 = lusine.finance.ewma_variance(variance = v4, sample, decay = one);
  let v6 = lusine.finance.ewma_variance(variance = v5, sample, decay = one);
  let v7 = lusine.finance.ewma_variance(variance = v6, sample, decay = one);
  let v8 = lusine.finance.ewma_variance(variance = v7, sample, decay = one);
  let v9 = lusine.finance.ewma_variance(variance = v8, sample, decay = one);
  let v10 = lusine.finance.ewma_variance(variance = v9, sample, decay = one);
  return v1 == start && v2 == start && v3 == start && v4 == start && v5 == start && v6 == start && v7 == start
      && v8 == start && v9 == start && v10 == start;
}

bool near(float got, float want) { return abs(got / want - 1.0) <= 16.0 * eps(); }

bool quiet(float start, float decay)
{
  let v1 = lusine.finance.ewma_variance(variance = start, sample = 0.0 * start, decay);
  let v2 = lusine.finance.ewma_variance(variance = v1, sample = 0.0 * start, decay);
  let v3 = lusine.finance.ewma_variance(variance = v2, sample = 0.0 * start, decay);
  let v4 = lusine.finance.ewma_variance(variance = v3, sample = 0.0 * start, decay);
  let v5 = lusine.finance.ewma_variance(variance = v4, sample = 0.0 * start, decay);
  let v6 = lusine.finance.ewma_variance(variance = v5, sample = 0.0 * start, decay);
  let v7 = lusine.finance.ewma_variance(variance = v6, sample = 0.0 * start, decay);
  let v8 = lusine.finance.ewma_variance(variance = v7, sample = 0.0 * start, decay);
  let v9 = lusine.finance.ewma_variance(variance = v8, sample = 0.0 * start, decay);
  let v10 = lusine.finance.ewma_variance(variance = v9, sample = 0.0 * start, decay);
  return near(got = v1, want = 0.000376) && near(got = v2, want = 0.00035344)
      && near(got = v3) && near(got = v4)
      && near(got = v5) && near(got = v6)
      && near(got = v7) && near(got = v8)
      && near(got = v9) && near(got = v10);
}

output table<lusine.widget> panels = {
  lusine.series { { panel = "RiskMetrics' Table 5.3", chart = "variance, %²", quantity = window.percent, steps = 20,
                    label = "the kernel's" } },
  lusine.series { { panel = "RiskMetrics' Table 5.3", chart = "variance, %²", quantity = window.printed, steps = 20,
                    label = "the table's" } },
  lusine.value { { panel = "RiskMetrics' Table 5.3", quantity = window.day, label = "day" } },
  lusine.value { { panel = "RiskMetrics' Table 5.3", quantity = window.percent, label = "the kernel's" } },
  lusine.value { { panel = "RiskMetrics' Table 5.3", quantity = window.printed, label = "the table's" } },
  lusine.series { { panel = "the weighted sum", chart = "decay 0.94, %²", quantity = run.variance, steps = 300,
                    label = "the kernel's" } },
  lusine.series { { panel = "the weighted sum", chart = "decay 0.94, %²", quantity = run.weighted, steps = 300,
                    label = "the history's" } },
  lusine.value { { panel = "the weighted sum", quantity = run.step, label = "sample" } },
  lusine.value { { panel = "the checks", quantity = checks.table, label = "RiskMetrics' Table 5.3" } },
  lusine.value { { panel = "the checks", quantity = checks.forecast, label = "the forecast, 0.473%" } },
  lusine.value { { panel = "the checks", quantity = checks.weighted, label = "the weighted sum" } },
  lusine.value { { panel = "the checks", quantity = checks.limits, label = "a decay of one, a quiet market" } },
};

void main()
{
  window.day = window.day >= 19 ? 0 : window.day + 1;
  window.variance = window.day == 0 ? on_day(numbers = days.number, values = days.sample, day = 0)
                                      * on_day(numbers = days.number, values = days.sample, day = 0)
                                    : lusine.finance.ewma_variance(window.variance,
                                                                   sample = on_day(numbers = days.number,
                                                                                   values = days.sample,
                                                                                   day = window.day),
                                                                   settings.decay);
  window.percent = 10000.0 * window.variance;
  window.printed = on_day(numbers = days.number, values = days.printed, day = window.day);
  checks.table = checks.table && abs(window.percent - window.printed) < 0.001;
  checks.forecast = checks.forecast && (window.day < 19 || abs(sqrt(window.percent) - 0.473) < 0.0005);

  run.step = run.step >= settings.length - 1 ? 0 : run.step + 1;
  delete from history where restarting(step = run.step);
  serieses.sample = drawn(settings.seed, step = run.step, serieses.index, settings.spread);
  serieses.variance = lusine.finance.ewma_variance(variance = run.step == 0 ? 0.0 : serieses.variance,
                                                   serieses.sample, serieses.decay);
  history.weight = history.weight * history.decay;
  insert into history select recorded(index, decay, sample) from serieses;
  serieses.weighted = weighting(serieses.index, series = history.series, square = history.square,
                                weight = history.weight);
  serieses.agrees = abs(serieses.variance / serieses.weighted - 1.0) <= 0.00001;
  run.variance = 10000.0 * (sum(serieses.variance) where serieses.index == 2);
  run.weighted = 10000.0 * (sum(serieses.weighted) where serieses.index == 2);
  checks.weighted = checks.weighted && disagreeing(agrees = serieses.agrees) == 0;
  checks.limits = checks.limits && kept(limit.start, limit.sample, limit.one)
                  && quiet(limit.start, decay = settings.decay);
  lamps.colour = lamp_colour(holds = checks.table && checks.forecast && checks.weighted && checks.limits);
}

The kernel's other proof (the same table, and four series of three hundred normal samples against the weighted sum in float64) is on its proof page.