Returns, Volatility and Risk Metrics
Turn a price series into returns, measure volatility without losing digits, and compute drawdown, Sharpe ratio, Value at Risk and expected shortfall, with every convention stated.
A price series is not yet a measurement
Every later part of this course produces a series of numbers over time: the value of a position, the price of a simulated asset, the profit and loss (PnL) of a strategy in a back-test. Sooner or later someone asks how risky that series is. "It went from 100 to 157 over four years" says nothing about the 33% fall in the middle, nor about how large a bad day can be.
Risk numbers answer those questions, and they are easy to compute badly. Three things go wrong again and again in real code. The data is in prices when the formula wants returns. The variance is computed with a formula that is exact on paper and loses every digit on a computer. And a single number such as "99% Value at Risk" is reported without the convention behind it, so two systems compute different numbers from the same data and both are "right".
This lesson builds the standard measurements one at a time, each in a short program, and states the convention of each one: returns and compounding, variance with Welford's algorithm, annualization, drawdown, the Sharpe and Sortino ratios, and historical Value at Risk and expected shortfall, including how to check a VaR model against what actually happened. All data is synthetic, generated by the course's deterministic Lcg, and nothing here says anything about any real market. It builds on the money arithmetic of Money in C++: from here on we work with ratios and statistics, where double is the right tool.
Simple returns, log returns and why log returns add
A return is the relative change of a value over one period. If a price goes from P0 to P1, there are two common definitions:
- the simple return
R = P1 / P0 - 1, which is the percentage change you would quote in words; - the log return
r = log(P1 / P0), the natural logarithm of the price ratio.
For small moves they are almost equal (log(1 + R) is close to R when R is small), but they combine differently over several periods. Simple returns compound: the growth over two days is (1 + R1) * (1 + R2) - 1, not R1 + R2. Log returns add: log(P2 / P0) = log(P2 / P1) + log(P1 / P0), because the logarithm of a product is the sum of the logarithms.
A concrete case shows why the difference matters. A price rises 10% and then falls 10%. The simple returns sum to zero, but 1.1 * 0.9 = 0.99: you are down 1%. The log returns, +0.0953 and -0.1054, sum to -0.0101, and exp(-0.0101) - 1 is exactly the 1% loss.
#include <cmath>
#include <print>
#include <vector>
int main() {
// Closing prices on five consecutive days (synthetic).
const std::vector<double> price{100.0, 110.0, 99.0, 104.0, 102.0};
double growth = 1.0; // product of (1 + simple return)
double log_sum = 0.0; // sum of log returns
double simple_sum = 0.0;
std::println("day price simple log");
for (std::size_t t = 1; t < price.size(); ++t) {
double simple = price[t] / price[t - 1] - 1.0;
double logret = std::log(price[t] / price[t - 1]);
growth *= 1.0 + simple;
log_sum += logret;
simple_sum += simple;
std::println("{:>3} {:>6.2f} {:>8.4f} {:>8.4f}", t, price[t], simple, logret);
}
std::println("sum of simple returns: {:.4f}", simple_sum);
std::println("compounded simple: {:.4f}", growth - 1.0);
std::println("sum of log returns: {:.4f}", log_sum);
std::println("exp(sum of log) - 1: {:.4f}", std::exp(log_sum) - 1.0);
std::println("log(last / first): {:.4f}",
std::log(price.back() / price.front()));
}Press Run to compile and run this example.
day price simple log
1 110.00 0.1000 0.0953
2 99.00 -0.1000 -0.1054
3 104.00 0.0505 0.0493
4 102.00 -0.0192 -0.0194
sum of simple returns: 0.0313
compounded simple: 0.0200
sum of log returns: 0.0198
exp(sum of log) - 1: 0.0200
log(last / first): 0.0198The price ends 2% above where it started. Adding the simple returns claims 3.13%, which is wrong. Compounding them gives the true 2%. The log returns add up to 0.0198, which is log(1.02), the log of the total growth, and converting back gives 2% again.
Which one to use depends on the question:
- To aggregate over time for one asset, log returns are convenient because they add. Most statistics in this lesson (mean, volatility, Sharpe ratio) are computed per period and then scaled, and that scaling is cleanest for quantities that add.
- To aggregate across assets at one moment, simple returns are right: a portfolio's simple return is the weighted average of its holdings' simple returns. Log returns do not have that property.
- For risk numbers such as VaR, which are about money lost, simple returns map directly to currency: a simple return of
-0.03on a position of 1,000,000 is a loss of 30,000.
Whichever you choose, write it down and use it consistently. A volatility computed from log returns and one computed from simple returns of the same series differ slightly, and comparing them is comparing two different things.
Variance in one pass without losing every digit
Volatility is the standard deviation of returns, the usual measure of how much a series moves around its average. The sample variance of n observations with mean m is s^2 = sum((x_i - m)^2) / (n - 1). Dividing by n - 1 instead of n (Bessel's correction) makes it an unbiased estimate of the variance of the distribution the data came from; dividing by n gives the population variance of the data itself. With thousands of daily returns the difference is negligible, with ten it is not, so say which one you use. This lesson uses n - 1 for sample statistics.
The definition needs two passes: one to find m, one to sum the squared deviations. Algebra turns it into one pass, s^2 = (sum(x^2) - sum(x)^2 / n) / (n - 1), and that version appears in many textbooks and spreadsheets. On paper the two are equal. In floating point the one-pass formula subtracts two huge, nearly equal numbers, and the subtraction cancels the leading digits that agree and leaves only rounding error. This is called catastrophic cancellation.
Welford's algorithm (B. P. Welford, 1962) keeps a running mean and a running sum of squared deviations M2, updating both with each new value:
n = n + 1; delta = x - mean; mean = mean + delta / n; M2 = M2 + delta * (x - mean_new)
The variance is M2 / (n - 1). Every quantity it adds is a deviation from the current mean, so large values never get squared directly. It is one pass, uses constant memory, and works on a stream whose length you do not know in advance, which is exactly the situation of a live risk monitor.
The data below is four values whose sample variance is exactly 30, shifted by a growing offset. Shifting data does not change its variance.
#include <print>
#include <span>
// Textbook one-pass formula: s^2 = (sum(x^2) - sum(x)^2 / n) / (n - 1).
double naive_variance(std::span<const double> xs) {
double sum = 0.0, sum_sq = 0.0;
for (double x : xs) {
double sq = x * x; // a separate statement, so no compiler fuses it into an FMA
sum += x;
sum_sq += sq;
}
double n = static_cast<double>(xs.size());
double mean_part = sum * sum / n;
return (sum_sq - mean_part) / (n - 1);
}
// Welford (1962): update the mean and the sum of squared deviations M2 together.
struct Welford {
long long n = 0;
double mean = 0.0, m2 = 0.0;
void add(double x) {
++n;
double delta = x - mean; // distance from the old mean
mean += delta / static_cast<double>(n);
m2 += delta * (x - mean); // old distance times new distance
}
double variance() const { return n > 1 ? m2 / static_cast<double>(n - 1) : 0.0; }
};
int main() {
// Four values whose sample variance is exactly 30, shifted by a large offset.
// Shifting data does not change its variance.
const double base[] = {4.0, 7.0, 13.0, 16.0};
std::println("{:>8} {:>14} {:>10}", "offset", "naive", "Welford");
for (double offset : {0.0, 1e4, 1e6, 1e8, 1e9}) {
double xs[4];
Welford w;
for (int i = 0; i < 4; ++i) {
xs[i] = base[i] + offset;
w.add(xs[i]);
}
std::println("{:>8.0e} {:>14.4f} {:>10.4f}", offset, naive_variance(xs),
w.variance());
}
}Press Run to compile and run this example.
offset naive Welford
0e+00 30.0000 30.0000
1e+04 30.0000 30.0000
1e+06 30.0000 30.0000
1e+08 29.3333 30.0000
1e+09 -170.6667 30.0000With an offset of 1e9 the textbook formula reports a negative variance, which is impossible, and its square root would be NaN. The values are around 1e9, so their squares are around 1e18 and their sum around 4e18. A double has 53 bits of significand, about 16 decimal digits, so near 4e18 the spacing between neighboring representable numbers is 512. The true difference between sum(x^2) and sum(x)^2 / n is 90, smaller than that spacing, so the result is noise. Welford's version computes deviations of size 10 and keeps full precision.
Daily returns have a mean near zero, so this failure seems remote for them. It is not remote for prices stored as large tick counts, account balances in cents, timestamps in nanoseconds or latencies, and those are exactly the series that also get averaged and "standard-deviated". Use Welford everywhere and the question never comes up. (Note the comment in naive_variance: some compilers fuse a * b + c into one fused multiply-add instruction, which rounds differently. Writing the square as its own statement keeps this demonstration identical on every compiler.)
Annualizing with the square root of time
Volatility is usually quoted per year, while data is often daily. The usual conversion is
sigma_annual = sigma_daily * sqrt(252) and mean_annual = mean_daily * 252.
The number 252 is a convention: it approximates the number of trading days in a year on many stock exchanges once weekends and holidays are removed. Other markets trade on more days (foreign exchange trades on weekdays, crypto markets every day), and systems use 252, 260 or 365 depending on what they measure. Nothing about it is a law. Use the number of observation periods per year in your data and state it.
The square root comes from an assumption, not from arithmetic. If daily returns are independent and identically distributed, the variance of a sum of 252 of them is 252 times the daily variance, so the standard deviation grows with sqrt(252), about 15.87. If returns are correlated from one day to the next (as they are for many smoothed or illiquid series), this rule is wrong: positive correlation makes the true annual volatility larger than the scaled number, negative correlation makes it smaller. William Sharpe's 1994 article on the ratio named after him states this assumption (zero serial correlation) explicitly, and Andrew Lo's 2002 paper shows how far off scaled ratios can be when it fails.
So sigma_annual = sigma_daily * sqrt(252) is fine for a quick comparison of series measured the same way, and it is a claim about the data that you could test, not a fact.
A synthetic price path and its summary statistics
To have a realistic-looking series to measure, we generate one. The standard model for a price that cannot go negative is geometric Brownian motion (GBM): over a step of length dt years,
S(t + dt) = S(t) * exp((mu - sigma^2 / 2) * dt + sigma * sqrt(dt) * Z),
where Z is a standard normal random number, mu is the drift (the expected growth rate per year) and sigma is the volatility per year. The - sigma^2 / 2 term makes the expected value of S(T) equal to S(0) * exp(mu * T); without it the median path would grow at mu and the average path faster. Log returns under GBM are normally distributed with mean (mu - sigma^2 / 2) * dt and standard deviation sigma * sqrt(dt). Lessons 8 to 10 rely on the same model.
The normals come from the course's Lcg through the Box-Muller transform: for independent uniforms u1, u2 in (0, 1), sqrt(-2 * log(u1)) * cos(2 * pi * u2) is a standard normal. (Box-Muller produces a second normal with sin that this lesson discards for brevity; Monte Carlo Simulation keeps both.) We never print values from std::normal_distribution, because the standard does not specify its algorithm, so different libraries give different numbers.
The program measures four things on the path:
- Excess return: the daily return minus the risk-free rate, the return available without risk (here 2% per year, divided by 252 per day). Performance ratios are about the reward for taking risk, so they subtract it.
- The Sharpe ratio: mean excess return divided by its standard deviation, annualized by
sqrt(252)(the mean scales by 252, the standard deviation bysqrt(252), so the ratio scales by252 / sqrt(252) = sqrt(252)). - The Sortino ratio: the same mean divided by the downside deviation instead, which counts only returns below the target (here the risk-free rate). Definitions vary; this one is
sqrt(sum(min(0, r_i)^2) / n)over allnexcess returnsr_i. The idea is that upside surprises are not what an investor fears. - The drawdown at day
tis the fraction lost since the highest value so far,1 - S(t) / max(S(0..t)). The maximum drawdown is the largest such fall, from a peak to a later trough.
#include <cmath>
#include <cstdint>
#include <numbers>
#include <print>
#include <vector>
struct Lcg {
std::uint64_t s;
std::uint32_t next() {
s = s * 6364136223846793005ULL + 1442695040888963407ULL;
return std::uint32_t(s >> 33);
}
double uniform() { return (next() + 0.5) / 2147483648.0; } // in (0, 1)
};
double std_normal(Lcg& g) { // Box-Muller; uses two uniforms, keeps one normal
double u1 = g.uniform(), u2 = g.uniform();
return std::sqrt(-2.0 * std::log(u1)) * std::cos(2.0 * std::numbers::pi * u2);
}
// Geometric Brownian motion sampled once per trading day.
std::vector<double> gbm_path(double s0, double mu, double sigma, int days, Lcg& g) {
const double dt = 1.0 / 252.0;
std::vector<double> path{s0};
for (int t = 0; t < days; ++t) {
double z = std_normal(g);
double step = (mu - 0.5 * sigma * sigma) * dt + sigma * std::sqrt(dt) * z;
path.push_back(path.back() * std::exp(step));
}
return path;
}
struct Welford {
long long n = 0;
double mean = 0.0, m2 = 0.0;
void add(double x) {
++n;
double delta = x - mean;
mean += delta / static_cast<double>(n);
m2 += delta * (x - mean);
}
double variance() const { return n > 1 ? m2 / static_cast<double>(n - 1) : 0.0; }
};
int main() {
Lcg g{2024};
const double mu = 0.06, sigma = 0.20; // annual drift and volatility of the model
auto path = gbm_path(100.0, mu, sigma, 4 * 252, g);
const double rf_daily = 0.02 / 252; // risk-free rate per day (simple)
Welford excess; // daily simple returns minus rf
double downside_sq = 0.0; // for Sortino: only shortfalls below rf count
double peak = path[0], max_dd = 0.0;
std::size_t peak_day = 0, dd_peak = 0, dd_trough = 0;
for (std::size_t t = 1; t < path.size(); ++t) {
double r = path[t] / path[t - 1] - 1.0 - rf_daily;
excess.add(r);
if (r < 0) downside_sq += r * r;
if (path[t] > peak) { peak = path[t]; peak_day = t; }
double dd = 1.0 - path[t] / peak; // fraction lost since the running peak
if (dd > max_dd) { max_dd = dd; dd_peak = peak_day; dd_trough = t; }
}
double n = static_cast<double>(excess.n);
double sd = std::sqrt(excess.variance());
double downside_dev = std::sqrt(downside_sq / n);
const double ann = std::sqrt(252.0);
std::println("days simulated: {}", excess.n);
std::println("final price: {:.2f}", path.back());
std::println("daily excess mean: {:.6f}", excess.mean);
std::println("daily stdev: {:.6f}", sd);
std::println("annualized volatility: {:.4f} (model: {:.4f})", sd * ann, sigma);
std::println("annualized Sharpe: {:.4f}", excess.mean / sd * ann);
std::println("annualized Sortino: {:.4f}", excess.mean / downside_dev * ann);
std::println("max drawdown: {:.2f}% (day {} to day {})", 100 * max_dd,
dd_peak, dd_trough);
}Press Run to compile and run this example.
days simulated: 1008
final price: 157.11
daily excess mean: 0.000450
daily stdev: 0.012760
annualized volatility: 0.2026 (model: 0.2000)
annualized Sharpe: 0.5601
annualized Sortino: 0.8084
max drawdown: 32.88% (day 466 to day 878)The estimated volatility, 0.2026, is close to the model's 0.20. That is a check on the generator and on the formula: with 1008 daily returns the standard error of a volatility estimate is about sigma / sqrt(2n), roughly 0.0045 here, so 0.2026 is well inside the expected range.
The Sharpe ratio is the opposite story. The model's true Sharpe ratio is about (0.06 - 0.02) / 0.20 = 0.2, and the estimate from four years is 0.56. Nothing is wrong with the code; four years of data simply cannot pin down a mean. The pitfall section below measures how bad this is. The maximum drawdown is 32.88%, from day 466 to day 878: a model with a positive drift still spent more than a year and a half falling a third from its peak.
Seeing the path and its drawdown
The same seed, drawn. The top panel is the price; the bottom panel is the drawdown, which is zero whenever the price is at a new high and dips below the axis while it is under water. SVG's y-axis points down, so the program maps higher prices to smaller y values itself.
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <format>
#include <numbers>
#include <print>
#include <string>
#include <vector>
struct Lcg {
std::uint64_t s;
std::uint32_t next() {
s = s * 6364136223846793005ULL + 1442695040888963407ULL;
return std::uint32_t(s >> 33);
}
double uniform() { return (next() + 0.5) / 2147483648.0; } // in (0, 1)
};
double std_normal(Lcg& g) {
double u1 = g.uniform(), u2 = g.uniform();
return std::sqrt(-2.0 * std::log(u1)) * std::cos(2.0 * std::numbers::pi * u2);
}
std::vector<double> gbm_path(double s0, double mu, double sigma, int days, Lcg& g) {
const double dt = 1.0 / 252.0;
std::vector<double> path{s0};
for (int t = 0; t < days; ++t) {
double z = std_normal(g);
double step = (mu - 0.5 * sigma * sigma) * dt + sigma * std::sqrt(dt) * z;
path.push_back(path.back() * std::exp(step));
}
return path;
}
int main() {
Lcg g{2024};
auto path = gbm_path(100.0, 0.06, 0.20, 4 * 252, g);
std::vector<double> dd; // drawdown at each day, as a fraction of the running peak
double peak = path[0];
for (double p : path) {
peak = std::max(peak, p);
dd.push_back(1.0 - p / peak);
}
const auto [lo, hi] = std::ranges::minmax(path);
const double max_dd = std::ranges::max(dd);
const double w = 400, x0 = 30, top_h = 130, dd_y0 = 160, dd_h = 70;
const double last = double(path.size() - 1);
auto x = [&](std::size_t t) { return x0 + (w - x0 - 5) * double(t) / last; };
// SVG's y axis points down: a higher price must get a smaller y.
auto y_price = [&](double p) { return 10 + top_h * (hi - p) / (hi - lo); };
auto y_dd = [&](double d) { return dd_y0 + dd_h * d / max_dd; };
std::string price_pts, dd_pts = std::format("{:.1f},{:.1f} ", x(0), dd_y0);
for (std::size_t t = 0; t < path.size(); ++t) {
price_pts += std::format("{:.1f},{:.1f} ", x(t), y_price(path[t]));
dd_pts += std::format("{:.1f},{:.1f} ", x(t), y_dd(dd[t]));
}
dd_pts += std::format("{:.1f},{:.1f}", x(path.size() - 1), dd_y0);
std::println(R"(<svg viewBox="0 0 400 240" font-size="9" fill="#1f2937">)");
std::println(R"(<polyline points="{}" fill="none" stroke="#2563eb"/>)", price_pts);
std::println(R"(<polygon points="{}" fill="#dc2626" fill-opacity="0.25" )"
R"(stroke="#dc2626" stroke-width="0.6"/>)", dd_pts);
std::println(R"(<line x1="{}" y1="{}" x2="395" y2="{}" stroke="#1f2937"/>)",
x0, dd_y0, dd_y0);
std::println(R"(<text x="0" y="{:.1f}">{:.0f}</text>)", y_price(hi) + 3, hi);
std::println(R"(<text x="0" y="{:.1f}">{:.0f}</text>)", y_price(lo), lo);
std::println(R"(<text x="0" y="{}">0%</text>)", dd_y0 + 3);
std::println(R"(<text x="0" y="{}">-{:.0f}%</text>)", dd_y0 + dd_h, 100 * max_dd);
std::println(R"(<text x="40" y="20">price</text>)");
std::println(R"(<text x="40" y="{}">drawdown</text>)", dd_y0 + dd_h + 8);
std::println("</svg>");
}Press Run to compile and run this example.
The deepest red trough sits near the right end of the chart, around day 878, and its start is the peak near the middle of the price panel. Most of the path is under water: drawdown is the normal state of a volatile series, and new highs are the exceptions. That is why a drawdown chart says more about how a series feels to hold than its volatility does.
Value at Risk and expected shortfall from history
Value at Risk (VaR) answers: "what loss will not be exceeded on a given fraction of days?" The 99% one-day VaR is a loss level L such that on 99% of days the loss is at most L. Statistically it is the 99th percentile (quantile) of the loss distribution, where loss is minus the return. Historical VaR takes that percentile directly from past returns, with no model of their distribution.
VaR says nothing about how large the losses beyond it are. Expected shortfall (ES, also called conditional VaR) does: it is the average loss on the days that are at least as bad as VaR. Acerbi and Tasche (2002) show that, unlike VaR, ES is a coherent risk measure; among other things, the ES of a combined portfolio is never larger than the sum of the parts' ES, which VaR does not guarantee.
"The 99th percentile of 200 numbers" has several reasonable definitions, so the convention must be stated. This lesson uses: sort the n losses from worst to best, let k = ceil(n * (1 - level)); VaR is the k-th worst loss and ES is the average of the k worst losses. With 1000 days at 99%, k = 10: VaR is the 10th worst day and ES the average of the 10 worst. Other systems interpolate between neighbors or use the (k+1)-th worst; the results differ slightly and all are defensible, as long as they are written down.
The program compares two synthetic return series with the same volatility. One is normal. The other is a mixture: on 95% of days the standard deviation is 1%, on the other 5% it is 4%, which gives the same overall standard deviation of about 1.32% but heavier tails. The last row computes the exact values for a normal distribution from the closed forms VaR = z * sd and ES = sd * pdf(z) / (1 - level), where z is the standard normal quantile (1.644854 at 95%, 2.326348 at 99%) and pdf the normal density.
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <functional>
#include <numbers>
#include <print>
#include <span>
#include <vector>
struct Lcg {
std::uint64_t s;
std::uint32_t next() {
s = s * 6364136223846793005ULL + 1442695040888963407ULL;
return std::uint32_t(s >> 33);
}
double uniform() { return (next() + 0.5) / 2147483648.0; } // in (0, 1)
};
double std_normal(Lcg& g) {
double u1 = g.uniform(), u2 = g.uniform();
return std::sqrt(-2.0 * std::log(u1)) * std::cos(2.0 * std::numbers::pi * u2);
}
struct Risk { double var, es; }; // both reported as positive losses
// Historical VaR and ES at confidence `level` (e.g. 0.99).
// Convention: k = ceil(n * (1 - level)); VaR is the k-th largest loss and
// ES is the average of the k largest losses (so ES >= VaR always).
Risk historical(std::span<const double> returns, double level) {
std::vector<double> loss(returns.size());
std::ranges::transform(returns, loss.begin(), std::negate<>{});
std::ranges::sort(loss, std::greater<>{}); // worst first
// 1 - 0.99 is 0.010000000000000009 in double: nudge down so that
// n = 200000 gives k = 2000, not 2001.
double tail = double(loss.size()) * (1.0 - level);
auto k = std::max<std::size_t>(1, static_cast<std::size_t>(std::ceil(tail - 1e-9)));
double tail_sum = 0.0;
for (std::size_t i = 0; i < k; ++i) tail_sum += loss[i];
return {loss[k - 1], tail_sum / double(k)};
}
void row(const char* name, double sd, Risk r95, Risk r99) {
std::println("{:<8} {:>5.2f}% {:>5.2f}% {:>5.2f}% {:>5.2f}% {:>5.2f}%", name,
100 * sd, 100 * r95.var, 100 * r95.es, 100 * r99.var, 100 * r99.es);
}
double rms(std::span<const double> r) { // the stdev when the mean is zero
double sum_sq = 0.0;
for (double x : r) sum_sq += x * x;
return std::sqrt(sum_sq / double(r.size()));
}
int main() {
Lcg g{11};
const std::size_t n = 200'000; // many days, so estimates sit near the true values
std::vector<double> calm, jumpy;
for (std::size_t i = 0; i < n; ++i) calm.push_back(0.013229 * std_normal(g));
for (std::size_t i = 0; i < n; ++i) {
// 95% of days are quiet (1% stdev); 5% of days are wild (4% stdev).
double sd = g.uniform() < 0.05 ? 0.04 : 0.01;
jumpy.push_back(sd * std_normal(g));
}
std::println("{:<8} {:>6} {:>6} {:>6} {:>6} {:>6}", "", "stdev", "VaR95", "ES95",
"VaR99", "ES99");
row("normal", rms(calm), historical(calm, 0.95), historical(calm, 0.99));
row("mixture", rms(jumpy), historical(jumpy, 0.95), historical(jumpy, 0.99));
// Exact values for a normal distribution: VaR = z * sd and
// ES = sd * pdf(z) / (1 - level), with z the standard normal quantile at `level`.
auto pdf = [](double z) {
return std::exp(-0.5 * z * z) / std::sqrt(2 * std::numbers::pi);
};
const double sd = 0.013229, z95 = 1.644854, z99 = 2.326348;
row("exact N", sd, {z95 * sd, sd * pdf(z95) / 0.05},
{z99 * sd, sd * pdf(z99) / 0.01});
}Press Run to compile and run this example.
stdev VaR95 ES95 VaR99 ES99
normal 1.32% 2.18% 2.74% 3.09% 3.55%
mixture 1.32% 1.80% 2.93% 3.44% 5.64%
exact N 1.32% 2.18% 2.73% 3.08% 3.53%The normal sample agrees with the exact row to within sampling noise, which checks the quantile code. Then the mixture: its 95% VaR is lower than the normal's (1.80% against 2.18%), so by that measure it looks safer. Its 99% VaR is a little higher (3.44% against 3.09%), and its 99% ES is much higher (5.64% against 3.55%). Volatility could not tell the two apart, the 95% VaR ranked them the wrong way round, and only the tail average shows how bad the bad days are. (For this mixture the exact 99% VaR and ES, computed separately by root-finding on its distribution function, are 3.445% and 5.607%.)
Notice the comment in historical. In double, 1 - 0.99 is 0.010000000000000009, so 200000 * (1 - 0.99) is slightly more than 2000 and ceil returns 2001. The first version of this program had that bug, and it moved the 99% numbers in the third decimal. Counting from a floating-point fraction needs a small nudge or an integer representation of the level.