Executive Summary
Volatility, VaR, expected shortfall and drawdown for Index Level
The short answer
Over 7.66 years, Index Level realized an annualized volatility of 25.69% with the most recent 21-period rolling estimate at 20.11%. One-period Value-at-Risk is 2.68% at 95% confidence and 3.60% at 99%; the expected shortfall—average loss when VaR is breached—is 3.42% and 4.85% respectively. The deepest drawdown was −55.10% from 2021-05-18 to 2023-06-19 and had not recovered by the sample end.
The detail
Volatility clusters significantly (Ljung-Box test on squared returns to lag 10: 662.70, p < 0.001), meaning quiet and turbulent stretches group together. The normal distribution assumption underperforms: excess kurtosis is 0.82 (fatter-tailed than normal), and at 95% the historical VaR was breached 107 times against 100.0 expected (ratio 1.07). At 99%, 18 breaches occurred against 20.0 expected (ratio 0.90). The worst single-period loss was 7.62% on 2022-03-31. Historical Value-at-Risk cannot exceed the largest loss in its own sample, and past volatility does not forecast future volatility: the first half (2016–2019-11-01) realized 23.92% annualized while the second half realized 27.36%, a difference of 3.45 percentage points (factor of 1.14).
What this can't tell you
The square-root-of-time annualization is defensible but depends on low autocorrelation; the lag-1 autocorrelation of −0.008 supports it here. The normal model's breach counts appear contradictory but are not: volatility clustering widens the unconditional normal fit, causing over-prediction in calm periods and under-prediction in turbulent ones, with counts canceling across the full sample. Breach counts themselves have sampling error; with 20.0 expected, the standard deviation is approximately 4.5.
Analysis Overview
Volatility, Value-at-Risk and drawdown for Index Level across 2,000 time points.
The short answer
Index Level is a price series spanning 7.66 years (2,000 observations), with log returns used for volatility modeling and simple returns for capital-loss measures. The annualization rate of 252 periods per year was derived from the observed rate of 261.2 per year and snapped to convention.
The detail
The series type was detected algorithmically: 0.0% of values are below 1 in absolute size, average absolute magnitude is 1229.111, and zero observations are non-positive—all consistent with a price series, not returns. Log returns (1,999 of them) are used for volatility, EWMA, and GARCH because they add across periods, which is necessary for square-root-of-time annualization. Simple returns carry Value-at-Risk, expected shortfall, and drawdown because those measure percentage of capital lost. The observed rate of 261.2 periods per year was snapped to 252, the standard convention. Sample-wide annualized volatility is 25.69%, and the 21-period rolling estimate ranges from 10.41% to 64.09%.
What this can't tell you
The series is not resampled to a regular calendar: gaps in the data produce returns spanning multiple periods, biasing volatility estimates upward. Every risk figure here is historical, not forward-looking; the sample's own two halves differ by 3.45 percentage points of annualized volatility (23.92% vs 27.36%), showing the error cost of treating past volatility as a forecast.
Data Preparation
How the series was assembled, and what was dropped or averaged.
The short answer
All 2,000 rows were retained with no missing dates or values, forming 1,999 returns. The series was not resampled to a regular calendar, so returns spanning data gaps inherit extra volatility.
The detail
Initial rows: 2,000. Final rows: 2,000. Rows dropped for incomplete Session Date or Index Level: 0. This produced 2,000 distinct time points ordered by Session Date and 1,999 returns. The series is NOT resampled to a regular calendar: returns are taken between consecutive available observations, so a gap in the data produces a return that spans it and is correspondingly more volatile than a single period's. The observed observation rate of 261.2 per year is what the annualization was derived from.
What this can't tell you
Any return spanning a gap in the trading calendar or data collection will appear more volatile than a single-period return. This biases every volatility figure (rolling, EWMA, GARCH, annualized) upward relative to what a calendar-resampled series would show.
Rolling Volatility
Annualized volatility of Index Level over time, by three estimators.
The short answer
Rolling 21-period volatility ranged from 10.41% to 64.09% over the sample, a 53.68 percentage point spread. The latest estimate stands at 20.11%, having peaked at 64.09%. All three estimators (rolling, EWMA, GARCH) trace the same turbulent periods but smooth differently.
The detail
Rolling 21-period standard deviation, annualized by √252, ranges from 10.41% (minimum) to 64.09% (maximum)—a spread of 53.68 percentage points showing volatility is an average over regimes, not a single state. Latest rolling estimate: 20.11%. EWMA (decay 0.924, fitted by MLE) reacts immediately and forgets smoothly, latest 22.30%. GARCH(1,1) conditional volatility is pulled toward a long-run level of 25.44% by persistence 0.962, latest 23.44%. Square-root-of-time annualization is defensible: lag-1 return autocorrelation is −0.008 (two-standard-error band ±0.044), and 5-period variance ratio is 1.02 (z = 0.15, independence = 1).
What this can't tell you
The rolling window produces step changes when large returns enter and leave; neither EWMA nor GARCH forecasts next period's return direction, only the width of the distribution. GARCH persistence of 0.962 means volatility shocks decay slowly and the model approaches the boundary where unconditional variance ceases to exist, making the long-run estimate unstable.
Drawdown
How far the series fell below its own running peak, and for how long.
The short answer
The deepest drawdown was −55.10%, spanning 762.0 days from peak 2021-05-18 to trough 2023-06-19. It had not recovered by the sample end on 2023-09-01. Breaking even from a −55.10% loss requires a 122.70% gain.
The detail
Drawdown is measured peak-to-trough on the path. Maximum drawdown: −55.10%, reached 2023-06-19 after peaking 2021-05-18 (762.0 days from peak to trough). Recovery status: not recovered by 2023-09-01 (at least 74.0 days further). The break-even gain from −55.10% is 122.70%, demonstrating that drawdown is the risk a holder actually experiences. Drawdown is the one risk number on this report measured on the path rather than on a distribution.
What this can't tell you
Drawdown is historical—the worst this sample experienced, not the worst that can happen. The unrecovered drawdown may recover after the sample end or may not; the sample window alone does not determine future recovery.
Value-at-Risk, Expected Shortfall and Drawdown
One-period loss measures, historical and parametric, with the drawdown.
| Measure | Method | Estimate | Units | Horizon | Interpretation |
|---|---|---|---|---|---|
| Value-at-Risk 95% | Historical | 2.684 | percent loss | one period | The loss exceeded on 5.0% of periods in this sample. |
| Value-at-Risk 99% | Historical | 3.597 | percent loss | one period | The loss exceeded on 1.0% of periods in this sample. |
| Value-at-Risk 95% | Parametric (normal) | 2.655 | percent loss | one period | What a normal distribution with this mean and standard deviation predicts. |
| Value-at-Risk 99% | Parametric (normal) | 3.758 | percent loss | one period | What a normal distribution with this mean and standard deviation predicts. |
| Expected shortfall 95% | Historical | 3.42 | percent loss | one period | The AVERAGE loss on the 5.0% of periods that breached VaR — the size of the bad day, not its frequency. |
| Expected shortfall 99% | Historical | 4.846 | percent loss | one period | The AVERAGE loss on the 1.0% of periods that breached VaR. |
| Expected shortfall 95% | Parametric (normal) | 3.331 | percent loss | one period | Normal-theory expected shortfall for the same fitted distribution. |
| Expected shortfall 99% | Parametric (normal) | 4.306 | percent loss | one period | Normal-theory expected shortfall for the same fitted distribution. |
| Worst single-period loss | Observed | 7.624 | percent loss | one period | The largest single-period fall in this sample, on 2022-03-31. |
| Maximum drawdown | Observed | 55.1 | percent loss | peak to trough | Peak 2021-05-18 to trough 2023-06-19; not recovered by the end of the sample on 2023-09-01. |
The short answer
One-period Value-at-Risk is 2.68% (95%) and 3.60% (99%); expected shortfall—the average loss when breached—is 3.42% and 4.85%, respectively. The worst single period lost 7.62%, which is 2.12 times the 99% VaR, showing VaR is a threshold, not a bound.
The detail
Historical 95% VaR: 2.68% (5% of periods exceeded this loss); expected shortfall 3.42% (0.74 pp worse). Historical 99% VaR: 3.60% (1% of periods exceeded); expected shortfall 4.85% (1.25 pp worse). Parametric (normal) 95% VaR: 2.6554%; 99% VaR: 3.7578%. Worst single-period loss: 7.62% on 2022-03-31, which is 2.12 times the historical 99% VaR. The gap between VaR and expected shortfall is the part of tail risk a VaR number hides. Historical figures assume no distribution and count what happened; parametric assume normality. At 99%, the historical VaR sits 0.16 percentage points below the normal one.
What this can't tell you
VaR is a quantile, not a worst case; losses can and do exceed it. The expected shortfall is the average loss given breach, not the maximum. Historical VaR cannot represent losses larger than 7.62%, the worst in this 7.66-year sample.
What The Normal Assumption Costs
Excess kurtosis, a normality test, and a count of the normal VaR's breaches.
| Diagnostic | Observed | Expected | Units | Verdict |
|---|---|---|---|---|
| Excess kurtosis of returns | 0.8175 | 0 | dimensionless | fatter-tailed than normal by 0.82 |
| Skewness of returns | -0.0865 | 0 | dimensionless | roughly symmetric |
| Jarque-Bera normality statistic | 58.15 | 0 | chi-square, 2 df | normality rejected (p < 0.001) |
| Breaches of the normal 95% VaR | 107 | 99.95 | periods | 107 observed against 100.0 expected, a ratio of 1.07 |
| Breaches of the normal 99% VaR | 18 | 19.99 | periods | 18 observed against 20.0 expected, a ratio of 0.90 |
| Historical minus normal VaR 95% | 0.0287 | 0 | percentage points | the normal model understates the historical loss |
| Historical minus normal VaR 99% | -0.1612 | 0 | percentage points | the normal model does not understate the historical loss |
The short answer
Returns are fatter-tailed than normal (excess kurtosis 0.82), yet the normal model's breach counts are near prediction: 107 observed vs 100.0 expected at 95%, and 18 observed vs 20.0 expected at 99%. This apparent contradiction reflects volatility clustering, which widens the unconditional fit.
The detail
Excess kurtosis: 0.82 (normal = 0). Skewness: −0.0865 (roughly symmetric). Jarque-Bera test: 58.15 on 2 df, p < 0.001 (normality rejected). Normal 95% VaR breaches: 107 observed, 100.0 expected (ratio 1.07). Normal 99% VaR breaches: 18 observed, 20.0 expected (ratio 0.90). Historical minus normal VaR at 95%: +0.0287 pp (historical is larger). At 99%: −0.1612 pp (normal is larger). The fat tails are real, yet an unconditional normal fitted to the whole sample is widened by turbulent stretches, so it over-predicts losses during calm periods and under-predicts during turbulent ones; counts cancel across the sample while daily risk does not. With 20.0 expected breaches, the sampling standard deviation is about 4.5, so ratios near 1 should not be over-read.
What this can't tell you
Autocorrelation in squared returns can arise from clustering, isolated extremes, or a shift in variance level. The breach count is itself a small number subject to sampling variation. The normal model's near-accurate count masks regime-dependent misprediction.
Volatility Clustering
Autocorrelation of squared returns, with a Ljung-Box test.
The short answer
All 10 lags of squared-return autocorrelation exceed the two-standard-error band (0.044), and the Ljung-Box test gives 662.70 (p < 0.001). Volatility clusters: turbulent and calm stretches group together, and today's volatility is informative about tomorrow's.
The detail
Squared-return autocorrelation at lag 1: 0.146; lag 2: 0.2479; lag 5: 0.2101; lag 10: 0.2603. All 10 lags shown fall outside the ±0.044 two-standard-error band. Ljung-Box test on squared returns to lag 10: 662.70, p < 0.001. Because volatility clusters, a single unconditional volatility number understates risk during turbulent stretches and overstates it during calm ones. The EWMA and GARCH on the models card are designed to capture this time-varying behavior.
What this can't tell you
Autocorrelation in squared returns is consistent with volatility clustering, but can also arise from a few isolated extreme observations or a shift in variance level partway through the sample. The test does not distinguish these mechanisms.
Volatility Models
EWMA and GARCH(1,1), both fitted by direct maximum likelihood.
| Parameter | Estimate | Display |
|---|---|---|
| EWMA decay (lambda, fitted by MLE) | 0.924 | 0.924 |
| EWMA latest annualized volatility | 22.3 | 22.30% |
| RiskMetrics EWMA (lambda 0.94) latest | 22.24 | 22.24% |
| GARCH omega | 9.73e-06 | 0.00000973 |
| GARCH alpha (news impact) | 0.0903 | 0.090 |
| GARCH beta (persistence of past variance) | 0.8719 | 0.872 |
| GARCH alpha + beta | 0.9621 | 0.962 |
| GARCH shock half-life (periods) | 17.95 | 18.0 |
| GARCH long-run annualized volatility | 25.44 | 25.44% |
| GARCH latest conditional volatility | 23.44 | 23.44% |
| GARCH one-period-ahead forecast | 22.44 | 22.44% |
| Likelihood ratio, GARCH against constant variance | 212.1 | 212.14 (p < 0.001) |
The short answer
EWMA decay fitted at 0.924 (interior to search range, not on boundary) puts latest volatility at 22.30%. GARCH(1,1) has persistence 0.962 (alpha 0.090, beta 0.872), a half-life of 18.0 periods, and improves on constant variance by likelihood-ratio 212.14 (p < 0.001). Both models beat a single volatility number.
The detail
EWMA: decay 0.924 (fitted by MLE, vs RiskMetrics convention 0.94). Latest annualized volatility 22.30% (convention would give 22.24%). GARCH(1,1): omega 0.00000973, alpha 0.090, beta 0.872, persistence 0.962. Shock half-life: 18.0 periods. Long-run annualized volatility: 25.44% (realized: 25.69%). Latest conditional volatility: 23.44%. One-period-ahead forecast: 22.44%. Likelihood-ratio test GARCH vs constant variance: 212.14 on 2 df, p < 0.001. Clustering is a real feature, not decoration. Persistence 0.962 means volatility shocks fade slowly; the model is close to the boundary where unconditional variance ceases to exist, making the long-run figure unstable well before reaching 1.
What this can't tell you
A one-period-ahead conditional volatility is a forecast of variance (distribution width), never of return direction. Persistence near 1 makes the long-run volatility estimate unstable and unreliable as a forecast horizon lengthens.
Method & Limits
Every formula, every setting, and the four limits on what these numbers mean.
| Item | Detail |
|---|---|
| Series type | The series type was detected, not assumed: 0.0% of the Index Level values are smaller than 1 in absolute size, their average absolute size is 1229.111, and 0 of them are zero or negative. A return series straddles zero and is small on both counts; a price or value series is not. On that evidence Index Level was read as a price or value series. |
| Returns | Log returns (the difference of logs) are used for volatility, the EWMA and the GARCH because they add across periods, which is what annualizing by a square root requires. Simple returns (the proportional change) are used for Value-at-Risk, expected shortfall and drawdown, because those numbers describe a percentage of capital lost. 1,999 return(s) were formed from 2,000 time points. |
| Annualization | The series carries 2,000 observations across 7.66 years, an observed rate of 261.2 per year, which was snapped to the nearest standard convention of 252. Volatility is annualized by multiplying the per-period standard deviation by the square root of 252. That step assumes returns are independent across periods; the independence check below is what decides whether it holds here. |
| Rolling volatility | Standard deviation of log returns over a moving window of 21 period(s), annualized the same way. Across this sample it ranged from 10.41% to 64.09% and ended at 20.11%. |
| Historical VaR | The empirical quantile of the simple returns: the 99.0% VaR is the loss that 1.0% of the periods in this sample exceeded. No distribution is assumed. It is bounded by the sample — see the row below. |
| Parametric VaR | Mean plus the normal quantile times the standard deviation of the simple returns (mean 0.005, standard deviation 1.618 per period). This is the standard textbook VaR and it is the one the tail diagnostics test. |
| Expected shortfall | The average loss GIVEN that VaR was breached. Historically it is the mean of the returns at or below the quantile; parametrically it is the normal-theory closed form. Expected shortfall is reported beside every VaR because a quantile says nothing about how far past it the loss goes. |
| Drawdown | Computed on the wealth index (the price series itself, or the compounded return series), as the largest fall from a running peak. Peak 2021-05-18, trough 2023-06-19, depth 55.10%, never recovered within this sample. |
| Clustering test | Ljung-Box test on the squared returns to lag 10: statistic 662.70, p < 0.001. 10 of the 10 lag autocorrelations sit outside the two-standard-error band of 0.044. |
| EWMA | Exponentially weighted variance with the decay found by maximizing the Gaussian likelihood over lambda in the range 0.70 to 0.995 — a one-dimensional search, not a fixed convention. Fitted lambda 0.924; the RiskMetrics convention of 0.94 is reported beside it for comparison. The fitted decay of 0.924 is interior to the search range, so it is a genuine maximum-likelihood estimate rather than a boundary. |
| GARCH(1,1) | Variance recursion h_t = omega + alpha e_{t-1} squared + beta h_{t-1}, fitted by direct maximum likelihood with base R's optim on an unconstrained reparameterisation (omega through a log, and the persistence and alpha's share of it through logistic transforms) so that positivity and stationarity hold by construction. The GARCH(1,1) improves on a constant variance by a likelihood-ratio statistic of 212.14 on 2 degrees of freedom (p < 0.001), so the clustering is a real feature of this series rather than a fitted decoration. |
| Packages | Base R plus stats only. No volatility or finance package is used: rugarch, fGarch, PerformanceAnalytics and quantmod are all absent from the analysis image, so the EWMA and the GARCH are implemented here directly. |
| What VaR is not | VaR is a quantile, not a worst case. It answers how bad a loss you clear on a given fraction of periods, and says nothing about the size of the losses beyond it — which is why expected shortfall is printed next to it everywhere. In this sample the worst single period lost 7.62%, which is 2.12 times the historical 99% VaR of 3.60%. |
| What history cannot see | Historical VaR cannot see a loss larger than the worst one in its sample. This sample covers 7.66 year(s), from 2016-01-04 to 2023-09-01, and the worst single period in it lost 7.62%. A window of that length has never observed an event rarer than roughly one in 1,999. |
| What the normal assumption costs | Excess kurtosis of the returns is 0.82, and the historical series breached the normal 99% VaR 18 time(s) against the 20.0 the normal model expects. The historical 99% VaR sits 0.16 percentage points below the parametric one. |
| Forecast limit | Split this sample in half at 2019-11-01: the first half realized 23.92% annualized and the second half 27.36%, a change of 3.45 percentage points (a factor of 1.14). An estimate made at the midpoint would have been wrong about the second half by that much, which is the honest size of the error in treating any of these numbers as a forecast. |
Full disclosure of every choice made above, and of the four limits that bound what these numbers can mean. First, parametric Value-at-Risk understates tail risk when returns are fat-tailed, and the size of that understatement here is measured rather than asserted: excess kurtosis 0.82, with 18 historical breaches of the normal 99% VaR against 20.0 expected. Second, historical Value-at-Risk cannot see a loss bigger than the worst one in its sample: this sample covers 7.66 years and its worst period lost 7.62%, so no event rarer than that is represented anywhere in these figures. Third, Value-at-Risk is a quantile and not a worst case — the loss beyond it is unbounded, which is why expected shortfall is reported beside every VaR on this report. Fourth, annualizing by the square root of 252 assumes returns are independent across periods. The square-root-of-time annualization is defensible here: the lag-1 autocorrelation of returns is -0.008 against a two-standard-error band of 0.044, and the 5-period variance ratio is 1.02 (z = 0.15) where independence implies 1. And underneath all four: past volatility is not future volatility. Split this sample in half at 2019-11-01: the first half realized 23.92% annualized and the second half 27.36%, a change of 3.45 percentage points (a factor of 1.14). An estimate made at the midpoint would have been wrong about the second half by that much, which is the honest size of the error in treating any of these numbers as a forecast.
Methodology
Statistical methodology and diagnostics for Volatility & Value-at-Risk
Statistical Method
Standard-library analysis: how risky is this series? Map a date column and a price, value or return series and get the full risk picture — log and simple returns with each used where it belongs, rolling volatility with a stated window and annualization, historical and parametric Value-at-Risk at 95% and 99% with expected shortfall beside each one, the maximum drawdown with its peak date, trough date and recovery time, and a volatility-clustering diagnostic backed by an EWMA and a GARCH(1,1) fitted by direct maximum likelihood. The limits are computed rather than boilerplated: the excess kurtosis and the actual count of historical breaches beyond the normal model's prediction, the worst loss the sample has ever seen and the span it covers, and the measured gap between the first and second halves of the same history.
- The series is ordered in time and each date carries one observation (several rows on a date are averaged)
- A price or value series is strictly positive, so that log returns exist; a return series is not required to be
- Observations are spaced regularly enough that one observed interval means the same thing throughout — gaps are not resampled and produce returns that span them
- Returns are independent across periods for the square-root-of-time annualization to hold; this is tested rather than assumed
- Returns are normally distributed for the parametric VaR and expected shortfall to be accurate; this is tested rather than assumed
- Parametric normal Value-at-Risk understates tail risk when returns are fat-tailed; the analysis reports the excess kurtosis and the actual count of historical breaches beyond the normal prediction rather than leaving this as a warning
- Historical Value-at-Risk cannot see a loss larger than the worst one in its sample — a two-year window has never observed a ten-year event, and the report states the span covered and the worst loss in it
- Value-at-Risk is a quantile, not a worst case: the loss beyond it is unbounded, which is why expected shortfall is reported alongside every VaR figure
- Annualizing volatility by the square root of the period count assumes returns are independent; when the autocorrelation and variance-ratio diagnostics say otherwise, the report states that the annualized figures are unreliable
Analysis Code
Complete R source code for this analysis
Volatility & Value-at-Risk — How Risky Is This Series?
Takes a dated price/value series (or returns directly — it detects which) and quantifies risk: rolling volatility with a stated window and annualization, historical and parametric Value-at-Risk with expected shortfall beside them, the maximum drawdown with its peak, trough and recovery, and a volatility-clustering diagnostic backed by an EWMA and a GARCH(1,1) fitted by direct maximum likelihood.
Why This Method?
A single standard deviation hides everything that matters about financial risk: that volatility clusters, that the loss distribution has fatter tails than a normal, and that the number people quote (VaR) is a quantile rather than a worst case. This module computes each of those gaps as a number instead of printing a warning about them.
What This Analysis Covers
- Log and simple returns, with each used where it belongs
- Rolling volatility, annualized by the square root of the inferred period
- Historical and parametric VaR at 95% and 99%, each with expected shortfall
- The measured cost of the normal assumption: excess kurtosis and the count
of historical breaches beyond the normal prediction
- Maximum drawdown with peak date, trough date and recovery time
- Volatility clustering: autocorrelation of squared returns, a Ljung-Box
test, an EWMA fitted by maximum likelihood, and a GARCH(1,1) fitted by direct optim maximum likelihood
Standard Library
Platform standard-library module (LAT-1441): runs on ANY dataset via the semantic mapping {date, value}. All narrative is derived from the user's own column names and computed values.
suppressPackageStartupMessages(library(DT))
suppressPackageStartupMessages(library(htmlwidgets))
suppressPackageStartupMessages(library(arrow))
suppressPackageStartupMessages(library(knitr))
suppressPackageStartupMessages(library(rmarkdown))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(tidyr))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(stringr))
suppressPackageStartupMessages(library(lubridate))
suppressPackageStartupMessages(library(broom))
suppressPackageStartupMessages(library(Matrix))
suppressPackageStartupMessages(library(cluster))
suppressPackageStartupMessages(library(data.table))Core Analysis Pipeline
Step 1: Resolve the mapped columns (humanized everywhere in prose)
initial_rows <- nrow(df)
date_h <- humanize_semantic("date", col_map)
value_h <- humanize_semantic("value", col_map)
for (req in c("date", "value")) {
if (!(req %in% names(df))) {
stop(sprintf("Volatility analysis needs '%s' (%s) mapped.",
humanize_semantic(req, col_map),
c(date = "the date column that orders the series",
value = "the numeric price, value or return series")[[req]]))
}
}Step 2: Coerce the value to numeric (95% rule)
v_raw <- df$value
if (!is.numeric(v_raw)) {
ch <- as.character(v_raw)
non_blank <- !is.na(ch) & trimws(ch) != ""
conv <- suppressWarnings(as.numeric(ch))
if (sum(non_blank) == 0 ||
sum(!is.na(conv[non_blank])) < 0.95 * sum(non_blank)) {
stop(sprintf("The value column '%s' does not look numeric — fewer than 95%% of its values parse as numbers. Map the numeric price, value or return series whose risk you want measured.",
value_h))
}
v_raw <- conv
}Step 3: Parse the date column (ISO, ymd, mdy, dmy, or integer years)
parse_dates_vec <- function(ch) {
ok <- function(dd) sum(!is.na(dd)) >= 0.95 * sum(!is.na(ch) & trimws(ch) != "")
d <- suppressWarnings(as.Date(ch, format = "%Y-%m-%d"))
if (!ok(d)) d <- suppressWarnings(as.Date(lubridate::ymd(ch, quiet = TRUE)))
if (!ok(d)) d <- suppressWarnings(as.Date(lubridate::mdy(ch, quiet = TRUE)))
if (!ok(d)) d <- suppressWarnings(as.Date(lubridate::dmy(ch, quiet = TRUE)))
if (!ok(d)) {
nv <- suppressWarnings(as.numeric(ch))
nv_ok <- stats::na.omit(nv)
if (length(nv_ok) >= 0.95 * sum(!is.na(ch) & trimws(ch) != "") &&
length(nv_ok) > 0 && all(nv_ok == round(nv_ok)) &&
all(nv_ok >= 1900 & nv_ok <= 2100)) {
d <- as.Date(ifelse(is.na(nv), NA, sprintf("%04d-01-01", nv)),
format = "%Y-%m-%d")
}
}
if (ok(d)) d else NULL
}
d_ch <- trimws(as.character(df$date))
dts <- parse_dates_vec(d_ch)
if (is.null(dts)) {
stop(sprintf("The date column '%s' could not be read as dates — fewer than 95%% of its values parse as calendar dates (or integer years). Map a date column so the series can be ordered in time.",
date_h))
}Step 4: Keep complete rows; average several rows on the same date
keep <- !is.na(v_raw) & !is.na(dts)
n_dropped <- sum(!keep)
v_k <- v_raw[keep]
d_k <- dts[keep]
final_rows <- length(v_k)
rows_removed <- initial_rows - final_rows
if (final_rows < 2) {
stop(sprintf("Only %d usable row(s) remained after dropping rows missing %s or %s — there is no series to analyse.",
final_rows, date_h, value_h))
}
agg <- stats::aggregate(list(v = v_k),
by = list(date_iso = format(d_k, "%Y-%m-%d")),
FUN = mean)
agg <- agg[order(agg$date_iso), , drop = FALSE]
n_points <- nrow(agg)
agg_note <- if (n_points < final_rows) {
sprintf("Several rows share the same %s, so the %s values on each date were averaged into one point(%s rows became %s time points).",
date_h, value_h, format(final_rows, big.mark = ","),
format(n_points, big.mark = ","))
} else ""
MIN_POINTS <- 30L
if (n_points < MIN_POINTS) {
stop(sprintf("The series in '%s' has only %d distinct time point(s) — volatility and Value-at-Risk need at least %d so that a return distribution and a rolling window exist at all.",
date_h, n_points, MIN_POINTS))
}
dates_iso <- agg$date_iso
vals <- agg$vStep 5: Is this a level series or is it already returns?
The test is arithmetic, not a guess, and its inputs are reported: returns straddle zero and are small; prices and values do not.
frac_small <- mean(abs(vals) < 1)
n_nonpos <- sum(vals <= 0)
mean_abs <- mean(abs(vals))
st_req <- tolower(as.character(params$series_type %||% "auto"))
if (!(st_req %in% c("auto", "price", "return", "returns", "level", "value"))) {
stop(sprintf("series_type must be one of auto, price or return; received '%s'.",
as.character(params$series_type)))
}
auto_type <- if (n_nonpos > 0 && frac_small >= 0.95 && mean_abs < 0.5) "return" else "price"
series_type <- switch(st_req,
auto = auto_type,
price = "price", level = "price", value = "price",
return = "return", returns = "return")
detect_note <- if (st_req == "auto") {
sprintf("The series type was detected, not assumed: %s%% of the %s values are smaller than 1 in absolute size, their average absolute size is %s, and %s of them are zero or negative. A return series straddles zero and is small on both counts; a price or value series is not. On that evidence %s was read as a %s series.",
r1f(100 * frac_small), value_h, r3(mean_abs),
format(n_nonpos, big.mark = ","), value_h,
if (series_type == "return") "return" else "price or value")
} else {
sprintf("The series type was set explicitly by the request rather than detected: %s was read as a %s series.",
value_h, if (series_type == "return") "return" else "price or value")
}Returns. Log returns are used for volatility, EWMA and GARCH because they add across periods, which is what the square-root-of-time rule needs. Simple returns are used for VaR, expected shortfall and drawdown because a loss of a given percentage of capital is what those numbers are meant to describe.
if (series_type == "price") {
if (n_nonpos > 0) {
stop(sprintf("The value column '%s' contains %s value(s) that are zero or negative, so it cannot be treated as a price or value level — a log return needs a positive series. If these numbers are already returns, set series_type to 'return'.",
value_h, format(n_nonpos, big.mark = ",")))
}
ret_dates <- dates_iso[-1]
log_ret <- diff(log(vals))
simple_ret <- exp(log_ret) - 1
level_index <- vals
level_dates <- dates_iso
} else {
simple_ret <- vals
bad <- (1 + simple_ret) <= 0
if (any(bad)) {
simple_ret <- simple_ret[!bad]
ret_dates <- dates_iso[!bad]
} else {
ret_dates <- dates_iso
}
log_ret <- log(1 + simple_ret)
level_index <- cumprod(1 + simple_ret)
level_dates <- ret_dates
}
n_ret <- length(log_ret)
if (n_ret < MIN_POINTS - 1L) {
stop(sprintf("Only %d usable return(s) could be formed from '%s' — at least %d are needed to estimate volatility.",
n_ret, value_h, MIN_POINTS - 1L))
}
sd_log <- stats::sd(log_ret)
if (!is.finite(sd_log) || sd_log <= 0) {
stop(sprintf("The value column '%s' is constant across all %s time points — a series that never moves has no volatility and no loss distribution to measure.",
value_h, format(n_points, big.mark = ",")))
}Step 6: Periods per year — observed, then snapped to a convention
span_days <- as.numeric(as.Date(dates_iso[n_points]) - as.Date(dates_iso[1]))
span_years <- span_days / 365.25
ppy_raw <- if (span_years > 0) n_points / span_years else NA_real_
ppy_req <- params$periods_per_year %||% NULL
conventions <- c(1, 4, 12, 26, 52, 252, 365)
if (!is.null(ppy_req)) {
ppy <- suppressWarnings(as.numeric(ppy_req))
if (is.na(ppy) || ppy <= 0) {
stop(sprintf("periods_per_year must be a positive number; received '%s'.",
as.character(ppy_req)))
}
ppy_note <- sprintf("The annualization rate was set explicitly to %s periods per year.", r2(ppy))
} else if (is.na(ppy_raw) || ppy_raw <= 0) {
ppy <- 252
ppy_note <- "The dates did not yield a usable observation rate, so the daily-trading convention of 252 periods per year was used."
} else {
near <- conventions[which.min(abs(log(ppy_raw / conventions)))]
if (abs(log(ppy_raw / near)) <= log(1.20)) {
ppy <- near
ppy_note <- sprintf("The series carries %s observations across %s years, an observed rate of %s per year, which was snapped to the nearest standard convention of %s.",
format(n_points, big.mark = ","), r2(span_years),
r1f(ppy_raw), format(ppy, big.mark = ","))
} else {
ppy <- ppy_raw
ppy_note <- sprintf("The series carries %s observations across %s years, an observed rate of %s per year. No standard convention was within 20%% of that, so the observed rate itself was used to annualize.",
format(n_points, big.mark = ","), r2(span_years), r1f(ppy_raw))
}
}
ann_factor <- sqrt(ppy)Step 7: Rolling volatility
w_req <- params$window %||% NULL
window <- if (!is.null(w_req)) {
wv <- suppressWarnings(as.integer(w_req))
if (is.na(wv) || wv < 3L) {
stop(sprintf("window must be a whole number of at least 3; received '%s'.",
as.character(w_req)))
}
wv
} else {
max(5L, min(21L, as.integer(floor(n_ret / 4))))
}
if (window > n_ret) window <- max(3L, as.integer(floor(n_ret / 2)))
roll_vol <- rep(NA_real_, n_ret)
csum <- c(0, cumsum(log_ret))
csum2 <- c(0, cumsum(log_ret * log_ret))
for (t in window:n_ret) {
i <- t - window + 1L
s1 <- csum[t + 1] - csum[i]
s2 <- csum2[t + 1] - csum2[i]
vv <- (s2 - s1 * s1 / window) / (window - 1)
roll_vol[t] <- if (is.finite(vv) && vv >= 0) sqrt(vv) * ann_factor else NA_real_
}
ann_vol <- 100 * sd_log * ann_factor
roll_fin <- which(is.finite(roll_vol))
roll_last <- if (length(roll_fin) > 0) 100 * roll_vol[max(roll_fin)] else NA_real_
roll_min <- if (length(roll_fin) > 0) 100 * min(roll_vol[roll_fin]) else NA_real_
roll_max <- if (length(roll_fin) > 0) 100 * max(roll_vol[roll_fin]) else NA_real_Step 8: Value-at-Risk and expected shortfall, historical and normal
Both are stated as POSITIVE percentage losses of the position over one period, computed on simple returns.
levels_conf <- c(0.95, 0.99)
mu_s <- mean(simple_ret)
sd_s <- stats::sd(simple_ret)
var_hist <- setNames(numeric(length(levels_conf)), as.character(levels_conf))
es_hist <- var_hist
var_norm <- var_hist
es_norm <- var_hist
for (i in seq_along(levels_conf)) {
a <- 1 - levels_conf[i]
q <- as.numeric(stats::quantile(simple_ret, probs = a, names = FALSE))
tail_vals <- simple_ret[simple_ret <= q]
var_hist[i] <- -100 * q
es_hist[i] <- if (length(tail_vals) > 0) -100 * mean(tail_vals) else NA_real_
z <- stats::qnorm(a)
var_norm[i] <- -100 * (mu_s + z * sd_s)
es_norm[i] <- -100 * (mu_s - sd_s * stats::dnorm(z) / a)
}Step 9: What the normal assumption actually costs here
Not a warning — a measurement. Excess kurtosis, a normality test, and a straight count of how often the historical series breached the loss the normal model said it would breach only rarely.
m_c <- log_ret - mean(log_ret)
m2 <- mean(m_c^2); m3 <- mean(m_c^3); m4 <- mean(m_c^4)
skew <- if (m2 > 0) m3 / m2^1.5 else NA_real_
excess_kurt <- if (m2 > 0) m4 / (m2 * m2) - 3 else NA_real_
jb_stat <- if (!is.na(skew) && !is.na(excess_kurt))
n_ret / 6 * (skew^2 + (excess_kurt^2) / 4) else NA_real_
jb_p <- if (!is.na(jb_stat)) stats::pchisq(jb_stat, df = 2, lower.tail = FALSE) else NA_real_
breach_obs <- setNames(integer(length(levels_conf)), as.character(levels_conf))
breach_exp <- setNames(numeric(length(levels_conf)), as.character(levels_conf))
for (i in seq_along(levels_conf)) {
breach_obs[i] <- sum(simple_ret < -var_norm[i] / 100)
breach_exp[i] <- (1 - levels_conf[i]) * n_ret
}
breach_ratio <- ifelse(breach_exp > 0, breach_obs / breach_exp, NA_real_)
fat_flag <- (!is.na(excess_kurt) && excess_kurt > 0.5) ||
any(!is.na(breach_ratio) & breach_ratio > 1.5)Step 10: What the historical method cannot see
worst_i <- {
fin <- which(is.finite(simple_ret))
if (length(fin) > 0) fin[which.min(simple_ret[fin])] else NA_integer_
}
worst_loss <- if (!is.na(worst_i)) -100 * simple_ret[worst_i] else NA_real_
worst_date <- if (!is.na(worst_i)) ret_dates[worst_i] else NA_character_
worst_over_var <- if (!is.na(worst_loss) && var_hist[["0.99"]] > 0)
worst_loss / var_hist[["0.99"]] else NA_real_Step 11: Drawdown on the wealth index
run_max <- cummax(level_index)
dd_path <- 100 * (level_index / run_max - 1)
dd_i <- {
fin <- which(is.finite(dd_path))
if (length(fin) > 0) fin[which.min(dd_path[fin])] else NA_integer_
}
dd_depth <- if (!is.na(dd_i)) dd_path[dd_i] else NA_real_
dd_trough_date <- if (!is.na(dd_i)) level_dates[dd_i] else NA_character_
dd_peak_i <- if (!is.na(dd_i)) {
cand <- which(level_index[seq_len(dd_i)] == run_max[dd_i])
if (length(cand) > 0) cand[1] else NA_integer_
} else NA_integer_
dd_peak_date <- if (!is.na(dd_peak_i)) level_dates[dd_peak_i] else NA_character_
dd_rec_i <- if (!is.na(dd_i) && dd_i < length(level_index)) {
cand <- which(level_index[(dd_i + 1):length(level_index)] >= run_max[dd_i])
if (length(cand) > 0) dd_i + cand[1] else NA_integer_
} else NA_integer_
dd_recovery_date <- if (!is.na(dd_rec_i)) level_dates[dd_rec_i] else NA_character_
dd_recovery_days <- if (!is.na(dd_rec_i))
as.numeric(as.Date(dd_recovery_date) - as.Date(dd_trough_date)) else NA_real_
dd_peak_to_trough_days <- if (!is.na(dd_peak_i) && !is.na(dd_i))
as.numeric(as.Date(dd_trough_date) - as.Date(dd_peak_date)) else NA_real_
dd_recovered <- !is.na(dd_rec_i)Step 12: Volatility clustering
Squared returns are the standard proxy for realised variance; if variance clusters, they are autocorrelated even when the returns themselves are not.
n_lag <- min(10L, max(2L, as.integer(floor(n_ret / 5))))
l2 <- (log_ret - mean(log_ret))^2
acf2 <- vapply(seq_len(n_lag), function(k) {
a <- l2[-seq_len(k)]; b <- l2[seq_len(n_ret - k)]
if (stats::sd(a) <= 0 || stats::sd(b) <= 0) return(NA_real_)
suppressWarnings(stats::cor(a, b))
}, numeric(1))
acf_band <- 1.96 / sqrt(n_ret)
lb <- tryCatch(stats::Box.test(l2, lag = n_lag, type = "Ljung-Box"),
error = function(e) NULL)
lb_stat <- if (!is.null(lb)) as.numeric(lb$statistic) else NA_real_
lb_p <- if (!is.null(lb)) as.numeric(lb$p.value) else NA_real_
cluster_flag <- !is.na(lb_p) && lb_p < 0.05
n_acf_out <- sum(!is.na(acf2) & abs(acf2) > acf_band)Step 13: Is the square-root-of-time rule safe here?
r1_ret <- if (n_ret >= 10) suppressWarnings(stats::cor(log_ret[-1], log_ret[-n_ret])) else NA_real_
r1_band <- 1.96 / sqrt(n_ret)
vr_q <- 5L
vr_m <- as.integer(floor(n_ret / vr_q))
vr_ratio <- NA_real_; vr_z <- NA_real_
if (vr_m >= 10) {
qr <- vapply(seq_len(vr_m), function(i)
sum(log_ret[((i - 1) * vr_q + 1):(i * vr_q)]), numeric(1))
denom <- vr_q * stats::var(log_ret)
if (is.finite(denom) && denom > 0) {
vr_ratio <- stats::var(qr) / denom
vr_se <- sqrt(2 * (vr_q - 1) / vr_m)
vr_z <- (vr_ratio - 1) / vr_se
}
}
indep_flag <- (!is.na(r1_ret) && abs(r1_ret) > r1_band) ||
(!is.na(vr_z) && abs(vr_z) > 2)Step 14: EWMA and GARCH(1,1), both fitted by maximum likelihood
Neither uses a volatility package. The EWMA decay is found by a one-dimensional likelihood search; the GARCH parameters come from optim on an unconstrained reparameterisation that keeps the model stationary.
e_c <- log_ret - mean(log_ret)
e2 <- e_c * e_c
ew <- tryCatch(stats::optimize(ewma_nll, interval = c(0.70, 0.995), e2 = e2),
error = function(e) NULL)
ewma_lambda <- if (!is.null(ew) && is.finite(ew$minimum)) ew$minimum else NA_real_
ewma_h <- if (!is.na(ewma_lambda)) ewma_filter(e2, ewma_lambda) else rep(NA_real_, n_ret)
ewma_last <- if (!is.na(ewma_lambda)) 100 * sqrt(ewma_h[n_ret]) * ann_factor else NA_real_
ewma_rm_h <- ewma_filter(e2, 0.94)
ewma_rm_last <- 100 * sqrt(ewma_rm_h[n_ret]) * ann_factorA decay that lands on an endpoint of the search range is not an estimate, it is the search giving up in that direction, and it is reported as such.
ewma_bound <- !is.na(ewma_lambda) &&
(ewma_lambda >= 0.995 - 1e-4 || ewma_lambda <= 0.70 + 1e-4)
ewma_note <- if (is.na(ewma_lambda)) {
"The EWMA decay could not be fitted on this series."
} else if (ewma_lambda >= 0.995 - 1e-4) {
sprintf("The fitted decay sits ON the upper end of the search range(%s), which is the likelihood saying it wants as much smoothing as it can have: no recent-past weighting beats treating the variance as constant. Read it as 'no usable EWMA signal here', not as a decay of %s.",
r3(ewma_lambda), r3(ewma_lambda))
} else if (ewma_lambda <= 0.70 + 1e-4) {
sprintf("The fitted decay sits ON the lower end of the search range(%s), meaning the likelihood wants to weight almost only the single most recent observation. That is usually a sign of a near-deterministic or piecewise-constant series rather than a genuine volatility process.",
r3(ewma_lambda))
} else {
sprintf("The fitted decay of %s is interior to the search range, so it is a genuine maximum-likelihood estimate rather than a boundary.",
r3(ewma_lambda))
}
GARCH_MIN <- 100L
gf <- if (n_ret >= GARCH_MIN) garch_fit(e2) else NULL
garch_ok <- !is.null(gf)
garch_omega <- if (garch_ok) gf$omega else NA_real_
garch_alpha <- if (garch_ok) gf$alpha else NA_real_
garch_beta <- if (garch_ok) gf$beta else NA_real_
garch_persist <- if (garch_ok) gf$persistence else NA_real_
garch_h <- if (garch_ok) gf$h else rep(NA_real_, n_ret)
garch_last <- if (garch_ok) 100 * sqrt(garch_h[n_ret]) * ann_factor else NA_real_
garch_fc <- if (garch_ok) {
hn <- garch_omega + garch_alpha * e2[n_ret] + garch_beta * garch_h[n_ret]
100 * sqrt(hn) * ann_factor
} else NA_real_
garch_lr_vol <- if (garch_ok && garch_persist < 1)
100 * sqrt(garch_omega / (1 - garch_persist)) * ann_factor else NA_real_
garch_halflife <- if (garch_ok && garch_persist > 0 && garch_persist < 1)
log(0.5) / log(garch_persist) else NA_real_
const_nll <- 0.5 * sum(log(mean(e2)) + e2 / mean(e2))Constant variance is the GARCH with alpha = beta = 0, so the two models are nested and the likelihood ratio cannot truly be negative. A small negative value only means the optimizer stopped just short of the nested optimum, so it is clamped to zero rather than reported as a negative statistic.
garch_lr_raw <- if (garch_ok) 2 * (const_nll - gf$nll) else NA_real_
garch_lr <- if (is.na(garch_lr_raw)) NA_real_ else max(0, garch_lr_raw)
garch_lr_p <- if (!is.na(garch_lr))
stats::pchisq(garch_lr, df = 2, lower.tail = FALSE) else NA_real_
garch_supported <- !is.na(garch_lr_p) && garch_lr_p < 0.05
garch_note <- if (!garch_ok) {
if (n_ret < GARCH_MIN)
sprintf("A GARCH(1,1) was not fitted: %s returns are available and at least %d are needed before its three parameters can be estimated with any stability.",
format(n_ret, big.mark = ","), GARCH_MIN)
else
"A GARCH(1,1) was attempted but no starting point produced a usable maximum-likelihood fit on this series, so no conditional-variance model is reported."
} else if (garch_supported) {
sprintf("The GARCH(1,1) improves on a constant variance by a likelihood-ratio statistic of %s on 2 degrees of freedom(%s), so the clustering is a real feature of this series rather than a fitted decoration.",
r2(garch_lr), fmt_pp(garch_lr_p))
} else {
sprintf("The GARCH(1,1) does NOT improve significantly on a constant variance here — the likelihood-ratio statistic is %s on 2 degrees of freedom(%s) — so its parameters should be read as a fit to noise rather than as evidence of clustering.",
r2(garch_lr), fmt_pp(garch_lr_p))
}Step 15: Past volatility is not future volatility — measured, not said
half <- as.integer(floor(n_ret / 2))
vol_h1 <- 100 * stats::sd(log_ret[seq_len(half)]) * ann_factor
vol_h2 <- 100 * stats::sd(log_ret[(half + 1):n_ret]) * ann_factor
half_gap <- vol_h2 - vol_h1
half_ratio <- if (is.finite(vol_h1) && vol_h1 > 0) vol_h2 / vol_h1 else NA_real_
h1_end <- ret_dates[half]
stability_note <- sprintf(
"Split this sample in half at %s: the first half realized %s annualized and the second half %s, a change of %s percentage points(a factor of %s). An estimate made at the midpoint would have been wrong about the second half by that much, which is the honest size of the error in treating any of these numbers as a forecast.",
h1_end, pct(vol_h1), pct(vol_h2), r2(abs(half_gap)),
if (is.na(half_ratio)) "n/a" else r2(half_ratio))Step 16: Chart datasets
thin <- function(n, cap) {
if (n <= cap) return(seq_len(n))
sort(unique(c(1L, n, as.integer(round(seq(1, n, length.out = cap))))))
}
vi <- window:n_ret
keep_v <- vi[thin(length(vi), 400L)]
vol_parts <- list(
data.frame(period_date = ret_dates[keep_v],
annualized_volatility = round(100 * roll_vol[keep_v], 3),
vol_series_label = sprintf("Rolling %d-period", window),
stringsAsFactors = FALSE)
)
if (!is.na(ewma_lambda)) {
vol_parts[[length(vol_parts) + 1]] <- data.frame(
period_date = ret_dates[keep_v],
annualized_volatility = round(100 * sqrt(ewma_h[keep_v]) * ann_factor, 3),
vol_series_label = sprintf("EWMA(lambda %s)", r3(ewma_lambda)),
stringsAsFactors = FALSE)
}
if (garch_ok) {
vol_parts[[length(vol_parts) + 1]] <- data.frame(
period_date = ret_dates[keep_v],
annualized_volatility = round(100 * sqrt(garch_h[keep_v]) * ann_factor, 3),
vol_series_label = "GARCH(1,1) conditional",
stringsAsFactors = FALSE)
}
volatility_series_df <- do.call(rbind, vol_parts)
volatility_series_df <-
volatility_series_df[is.finite(volatility_series_df$annualized_volatility), , drop = FALSE]
rownames(volatility_series_df) <- NULL
keep_d <- thin(length(dd_path), 700L)
drawdown_series_df <- data.frame(
period_date = level_dates[keep_d],
drawdown_pct = round(dd_path[keep_d], 3),
stringsAsFactors = FALSE)
rownames(drawdown_series_df) <- NULLStep 17: Result tables
lvl_lab <- paste0(format(100 * levels_conf, trim = TRUE), "%")
risk_measures_df <- data.frame(
measure = c(paste("Value-at-Risk", lvl_lab[1]), paste("Value-at-Risk", lvl_lab[2]),
paste("Value-at-Risk", lvl_lab[1]), paste("Value-at-Risk", lvl_lab[2]),
paste("Expected shortfall", lvl_lab[1]), paste("Expected shortfall", lvl_lab[2]),
paste("Expected shortfall", lvl_lab[1]), paste("Expected shortfall", lvl_lab[2]),
"Worst single-period loss", "Maximum drawdown"),
method = c(rep("Historical", 2), rep("Parametric(normal)", 2),
rep("Historical", 2), rep("Parametric(normal)", 2),
"Observed", "Observed"),
estimate = round(c(var_hist[1], var_hist[2], var_norm[1], var_norm[2],
es_hist[1], es_hist[2], es_norm[1], es_norm[2],
worst_loss, abs(dd_depth)), 4),
units = "percent loss",
horizon = c(rep("one period", 8), "one period", "peak to trough"),
stringsAsFactors = FALSE)
risk_measures_df$interpretation <- c(
sprintf("The loss exceeded on %s%% of periods in this sample.", r1f(100 * (1 - levels_conf[1]))),
sprintf("The loss exceeded on %s%% of periods in this sample.", r1f(100 * (1 - levels_conf[2]))),
"What a normal distribution with this mean and standard deviation predicts.",
"What a normal distribution with this mean and standard deviation predicts.",
sprintf("The AVERAGE loss on the %s%% of periods that breached VaR — the size of the bad day, not its frequency.", r1f(100 * (1 - levels_conf[1]))),
sprintf("The AVERAGE loss on the %s%% of periods that breached VaR.", r1f(100 * (1 - levels_conf[2]))),
"Normal-theory expected shortfall for the same fitted distribution.",
"Normal-theory expected shortfall for the same fitted distribution.",
sprintf("The largest single-period fall in this sample, on %s.",
if (is.na(worst_date)) "n/a" else worst_date),
sprintf("Peak %s to trough %s; %s.",
if (is.na(dd_peak_date)) "n/a" else dd_peak_date,
if (is.na(dd_trough_date)) "n/a" else dd_trough_date,
if (dd_recovered) sprintf("recovered by %s", dd_recovery_date)
else sprintf("not recovered by the end of the sample on %s", level_dates[length(level_dates)])))
rownames(risk_measures_df) <- NULL
gap95 <- var_hist[["0.95"]] - var_norm[["0.95"]]
gap99 <- var_hist[["0.99"]] - var_norm[["0.99"]]
tail_diagnostics_df <- data.frame(
diagnostic = c("Excess kurtosis of returns",
"Skewness of returns",
"Jarque-Bera normality statistic",
sprintf("Breaches of the normal %s VaR", lvl_lab[1]),
sprintf("Breaches of the normal %s VaR", lvl_lab[2]),
sprintf("Historical minus normal VaR %s", lvl_lab[1]),
sprintf("Historical minus normal VaR %s", lvl_lab[2])),
observed = round(c(excess_kurt, skew, jb_stat,
breach_obs[1], breach_obs[2], gap95, gap99), 4),
expected = round(c(0, 0, 0, breach_exp[1], breach_exp[2], 0, 0), 3),
units = c("dimensionless", "dimensionless", "chi-square, 2 df",
"periods", "periods", "percentage points", "percentage points"),
verdict = c(
if (is.na(excess_kurt)) "not computed"
else if (excess_kurt > 0.5) sprintf("fatter-tailed than normal by %s", r2(excess_kurt))
else if (excess_kurt < -0.5) sprintf("thinner-tailed than normal by %s", r2(abs(excess_kurt)))
else "close to normal",
if (is.na(skew)) "not computed"
else if (skew < -0.3) "left-skewed — large falls outweigh large rises"
else if (skew > 0.3) "right-skewed — large rises outweigh large falls"
else "roughly symmetric",
if (is.na(jb_p)) "not computed"
else if (jb_p < 0.05) sprintf("normality rejected(%s)", fmt_pp(jb_p))
else sprintf("normality not rejected(%s)", fmt_pp(jb_p)),
sprintf("%s observed against %s expected, a ratio of %s",
format(breach_obs[1], big.mark = ","), r1f(breach_exp[1]), r2(breach_ratio[1])),
sprintf("%s observed against %s expected, a ratio of %s",
format(breach_obs[2], big.mark = ","), r1f(breach_exp[2]), r2(breach_ratio[2])),
if (gap95 > 0) "the normal model understates the historical loss"
else "the normal model does not understate the historical loss",
if (gap99 > 0) "the normal model understates the historical loss"
else "the normal model does not understate the historical loss"),
stringsAsFactors = FALSE)
rownames(tail_diagnostics_df) <- NULL
squared_return_acf_df <- data.frame(
lag = as.character(seq_len(n_lag)),
autocorrelation = round(acf2, 4),
band_95 = round(acf_band, 4),
outside_band = ifelse(is.na(acf2), "not computed",
ifelse(abs(acf2) > acf_band, "outside", "inside")),
stringsAsFactors = FALSE)
rownames(squared_return_acf_df) <- NULL
model_fits_df <- data.frame(
parameter = c("EWMA decay(lambda, fitted by MLE)",
"EWMA latest annualized volatility",
"RiskMetrics EWMA(lambda 0.94) latest",
"GARCH omega",
"GARCH alpha(news impact)",
"GARCH beta(persistence of past variance)",
"GARCH alpha + beta",
"GARCH shock half-life(periods)",
"GARCH long-run annualized volatility",
"GARCH latest conditional volatility",
"GARCH one-period-ahead forecast",
"Likelihood ratio, GARCH against constant variance"),
estimate = round(c(ewma_lambda, ewma_last, ewma_rm_last,
garch_omega, garch_alpha, garch_beta, garch_persist,
garch_halflife, garch_lr_vol, garch_last, garch_fc,
garch_lr), 8),
display = c(r3(ewma_lambda), pct(ewma_last), pct(ewma_rm_last),
r8(garch_omega), r3(garch_alpha), r3(garch_beta),
r3(garch_persist), r1f(garch_halflife), pct(garch_lr_vol),
pct(garch_last), pct(garch_fc),
sprintf("%s(%s)", r2(garch_lr), fmt_pp(garch_lr_p))),
stringsAsFactors = FALSE)
rownames(model_fits_df) <- NULL
methods_df <- data.frame(
item = c("Series type", "Returns", "Annualization", "Rolling volatility",
"Historical VaR", "Parametric VaR", "Expected shortfall",
"Drawdown", "Clustering test", "EWMA", "GARCH(1,1)",
"Packages", "What VaR is not", "What history cannot see",
"What the normal assumption costs", "Forecast limit"),
detail = c(
detect_note,
sprintf("Log returns(the difference of logs) are used for volatility, the EWMA and the GARCH because they add across periods, which is what annualizing by a square root requires. Simple returns(the proportional change) are used for Value-at-Risk, expected shortfall and drawdown, because those numbers describe a percentage of capital lost. %s return(s) were formed from %s time points.",
format(n_ret, big.mark = ","), format(n_points, big.mark = ",")),
sprintf("%s Volatility is annualized by multiplying the per-period standard deviation by the square root of %s. That step assumes returns are independent across periods; the independence check below is what decides whether it holds here.",
ppy_note, format(ppy, big.mark = ",")),
sprintf("Standard deviation of log returns over a moving window of %d period(s), annualized the same way. Across this sample it ranged from %s to %s and ended at %s.",
window, pct(roll_min), pct(roll_max), pct(roll_last)),
sprintf("The empirical quantile of the simple returns: the %s%% VaR is the loss that %s%% of the periods in this sample exceeded. No distribution is assumed. It is bounded by the sample — see the row below.",
r1f(100 * levels_conf[2]), r1f(100 * (1 - levels_conf[2]))),
sprintf("Mean plus the normal quantile times the standard deviation of the simple returns(mean %s, standard deviation %s per period). This is the standard textbook VaR and it is the one the tail diagnostics test.",
r3(100 * mu_s), r3(100 * sd_s)),
sprintf("The average loss GIVEN that VaR was breached. Historically it is the mean of the returns at or below the quantile; parametrically it is the normal-theory closed form. Expected shortfall is reported beside every VaR because a quantile says nothing about how far past it the loss goes."),
sprintf("Computed on the wealth index(the price series itself, or the compounded return series), as the largest fall from a running peak. Peak %s, trough %s, depth %s, %s.",
if (is.na(dd_peak_date)) "n/a" else dd_peak_date,
if (is.na(dd_trough_date)) "n/a" else dd_trough_date,
pct(abs(dd_depth)),
if (dd_recovered) sprintf("recovered %s after %s period-days", dd_recovery_date, r1f(dd_recovery_days))
else "never recovered within this sample"),
sprintf("Ljung-Box test on the squared returns to lag %d: statistic %s, %s. %s of the %d lag autocorrelations sit outside the two-standard-error band of %s.",
n_lag, r2(lb_stat), fmt_pp(lb_p), n_acf_out, n_lag, r3(acf_band)),
sprintf("Exponentially weighted variance with the decay found by maximizing the Gaussian likelihood over lambda in the range 0.70 to 0.995 — a one-dimensional search, not a fixed convention. Fitted lambda %s; the RiskMetrics convention of 0.94 is reported beside it for comparison. %s",
r3(ewma_lambda), ewma_note),
sprintf("Variance recursion h_t = omega + alpha e_{t-1} squared + beta h_{t-1}, fitted by direct maximum likelihood with base R's optim on an unconstrained reparameterisation (omega through a log, and the persistence and alpha's share of it through logistic transforms) so that positivity and stationarity hold by construction. %s",
garch_note),
"Base R plus stats only. No volatility or finance package is used: rugarch, fGarch, PerformanceAnalytics and quantmod are all absent from the analysis image, so the EWMA and the GARCH are implemented here directly.",
sprintf("VaR is a quantile, not a worst case. It answers how bad a loss you clear on a given fraction of periods, and says nothing about the size of the losses beyond it — which is why expected shortfall is printed next to it everywhere. In this sample the worst single period lost %s, which is %s times the historical %s VaR of %s.",
pct(worst_loss),
if (is.na(worst_over_var)) "an undetermined number of" else r2(worst_over_var),
lvl_lab[2], pct(var_hist[["0.99"]])),
sprintf("Historical VaR cannot see a loss larger than the worst one in its sample. This sample covers %s year(s), from %s to %s, and the worst single period in it lost %s. A window of that length has never observed an event rarer than roughly one in %s.",
r2(span_years), dates_iso[1], dates_iso[n_points], pct(worst_loss),
format(n_ret, big.mark = ",")),
sprintf("Excess kurtosis of the returns is %s, and the historical series breached the normal %s VaR %s time(s) against the %s the normal model expects. The historical %s VaR sits %s percentage points %s the parametric one.",
r2(excess_kurt), lvl_lab[2], format(breach_obs[2], big.mark = ","),
r1f(breach_exp[2]), lvl_lab[2], r2(abs(gap99)),
if (gap99 > 0) "above" else "below"),
stability_note),
stringsAsFactors = FALSE)
rownames(methods_df) <- NULLStep 18: Headline metrics and the computed answer
metrics <- list(
`Time Points` = n_points,
`Returns` = as.integer(n_ret),
`Periods Per Year` = as.numeric(ppy),
`Annualized Volatility` = pct(ann_vol),
`Latest Rolling Volatility` = pct(roll_last),
`Historical VaR 95%` = pct(var_hist[["0.95"]]),
`Expected Shortfall 95%` = pct(es_hist[["0.95"]]),
`Historical VaR 99%` = pct(var_hist[["0.99"]]),
`Maximum Drawdown` = pct(dd_depth),
`Excess Kurtosis` = if (is.na(excess_kurt)) "n/a" else r2(excess_kurt),
`GARCH Persistence` = if (is.na(garch_persist)) "not fitted" else r3(garch_persist)
)
indep_clause <- if (indep_flag) {
sprintf("CAUTION: annualizing by the square root of %s assumes returns are independent across periods, and this series is not — the lag-1 autocorrelation of returns is %s against a two-standard-error band of %s, and the %d-period variance ratio is %s(z = %s) where independence implies 1. Every annualized figure here is therefore unreliable and should be read at the raw per-period scale instead.",
format(ppy, big.mark = ","), r3(r1_ret), r3(r1_band), vr_q,
r2(vr_ratio), r2(vr_z))
} else {
sprintf("The square-root-of-time annualization is defensible here: the lag-1 autocorrelation of returns is %s against a two-standard-error band of %s, and the %d-period variance ratio is %s(z = %s) where independence implies 1.",
r3(r1_ret), r3(r1_band), vr_q, r2(vr_ratio), r2(vr_z))
}
cluster_clause <- if (cluster_flag) {
sprintf("Volatility clusters: the Ljung-Box test on squared returns to lag %d gives %s(%s), so quiet and turbulent stretches group together and today's volatility is informative about tomorrow's.",
n_lag, r2(lb_stat), fmt_pp(lb_p))
} else {
sprintf("No volatility clustering was detected: the Ljung-Box test on squared returns to lag %d gives %s(%s), so this sample gives no evidence that turbulent periods group together.",
n_lag, r2(lb_stat), fmt_pp(lb_p))
}The cost of the normal assumption is stated as the measured numbers at BOTH levels, so the sentence can never assert a direction its own figures contradict.
breach_bits <- paste(vapply(seq_along(levels_conf), function(i)
sprintf("at %s, %s breach(es) against %s expected(a ratio of %s)",
lvl_lab[i], format(breach_obs[i], big.mark = ","),
r1f(breach_exp[i]), r2(breach_ratio[i])),
character(1)), collapse = "; ")
kurt_word <- if (is.na(excess_kurt)) "could not be computed"
else if (excess_kurt > 0.5) sprintf("%s, fatter-tailed than a normal", r2(excess_kurt))
else if (excess_kurt < -0.5) sprintf("%s, thinner-tailed than a normal", r2(excess_kurt))
else sprintf("%s, close to a normal", r2(excess_kurt))Fat tails and a low breach count can appear together, and when they do the reason is worth saying rather than papering over.
n_lvl <- length(levels_conf)
tail_nuance <- if (!is.na(excess_kurt) && excess_kurt > 0.5 &&
!is.na(breach_ratio[n_lvl]) && breach_ratio[n_lvl] <= 1) {
" Those two findings look contradictory and are not: the returns are fat-tailed, yet the normal VaR was breached no more often than predicted. That is what volatility clustering does to an unconditional fit — one normal fitted to the whole sample is widened by the turbulent stretches, so it over-predicts losses during the calm ones and under-predicts them during the turbulent ones, and the counts cancel out across the sample while the risk on any given day does not."
} else ""
fat_verdict <- if (fat_flag) {
"On that evidence the normal model under-predicts losses on at least one of these levels, so the historical figures are the ones to trust here."
} else {
"On that evidence the normal model is not badly wrong on this particular sample — which is a statement about this history, not a licence to trust the normal assumption on the next one."
}
fat_clause <- sprintf("What the normal assumption costs is measured here rather than assumed: excess kurtosis is %s, and the historical series breached the parametric normal VaR %s. %s%s",
kurt_word, breach_bits, fat_verdict, tail_nuance)
json_output <- list(
answer = paste0(
"Risk of ", value_h, " across ", format(n_points, big.mark = ","),
" time points in ", date_h, " (", dates_iso[1], " to ",
dates_iso[n_points], ", ", r2(span_years), " years): annualized volatility ",
pct(ann_vol), " (rolling ", window, "-period, latest ", pct(roll_last),
"); one-period historical Value-at-Risk ", pct(var_hist[["0.95"]]),
" at 95% and ", pct(var_hist[["0.99"]]), " at 99%, with expected shortfall ",
pct(es_hist[["0.95"]]), " and ", pct(es_hist[["0.99"]]),
" — the average loss GIVEN the VaR is breached, which the VaR alone hides. ",
"The parametric normal VaR at 99% is ", pct(var_norm[["0.99"]]), ". ",
"Maximum drawdown ", pct(dd_depth), " from ",
if (is.na(dd_peak_date)) "n/a" else dd_peak_date, " to ",
if (is.na(dd_trough_date)) "n/a" else dd_trough_date,
if (dd_recovered) paste0(", recovered ", dd_recovery_date) else ", not yet recovered",
". ", cluster_clause, " ", fat_clause, " ", indep_clause,
" Historical VaR cannot see a loss larger than the worst in its sample, which here lost ",
pct(worst_loss), " on ", if (is.na(worst_date)) "n/a" else worst_date,
"; and past volatility is not future volatility. ", stability_note
),
cards = lapply(
c("tldr", "overview", "preprocessing", "rolling_volatility", "drawdown",
"risk_measures", "tail_gap", "clustering", "models", "methods"),
function(cid) list(id = cid, metrics = metrics)
)
)
list(
initial_rows = initial_rows, final_rows = final_rows,
rows_removed = rows_removed, n_dropped = n_dropped,
date_h = date_h, value_h = value_h,
n_points = n_points, n_ret = n_ret, agg_note = agg_note,
dates_iso = dates_iso, ret_dates = ret_dates, level_dates = level_dates,
series_type = series_type, detect_note = detect_note,
frac_small = frac_small, mean_abs = mean_abs, n_nonpos = n_nonpos,
ppy = ppy, ppy_raw = ppy_raw, ppy_note = ppy_note,
ann_factor = ann_factor, window = window,
span_days = span_days, span_years = span_years,
ann_vol = ann_vol, roll_last = roll_last,
roll_min = roll_min, roll_max = roll_max,
mu_s = mu_s, sd_s = sd_s, levels_conf = levels_conf, lvl_lab = lvl_lab,
var_hist = var_hist, var_norm = var_norm,
es_hist = es_hist, es_norm = es_norm,
gap95 = gap95, gap99 = gap99,
skew = skew, excess_kurt = excess_kurt, jb_stat = jb_stat, jb_p = jb_p,
breach_obs = breach_obs, breach_exp = breach_exp,
breach_ratio = breach_ratio, fat_flag = fat_flag,
worst_loss = worst_loss, worst_date = worst_date,
worst_over_var = worst_over_var,
dd_depth = dd_depth, dd_peak_date = dd_peak_date,
dd_trough_date = dd_trough_date, dd_recovery_date = dd_recovery_date,
dd_recovery_days = dd_recovery_days, dd_recovered = dd_recovered,
dd_peak_to_trough_days = dd_peak_to_trough_days,
n_lag = n_lag, acf2 = acf2, acf_band = acf_band, n_acf_out = n_acf_out,
lb_stat = lb_stat, lb_p = lb_p, cluster_flag = cluster_flag,
r1_ret = r1_ret, r1_band = r1_band, vr_q = vr_q,
vr_ratio = vr_ratio, vr_z = vr_z, indep_flag = indep_flag,
ewma_lambda = ewma_lambda, ewma_last = ewma_last, ewma_rm_last = ewma_rm_last,
ewma_bound = ewma_bound, ewma_note = ewma_note,
kurt_word = kurt_word, tail_nuance = tail_nuance, breach_bits = breach_bits,
garch_ok = garch_ok, garch_omega = garch_omega, garch_alpha = garch_alpha,
garch_beta = garch_beta, garch_persist = garch_persist,
garch_last = garch_last, garch_fc = garch_fc, garch_lr_vol = garch_lr_vol,
garch_halflife = garch_halflife, garch_lr = garch_lr,
garch_lr_p = garch_lr_p, garch_supported = garch_supported,
garch_note = garch_note,
vol_h1 = vol_h1, vol_h2 = vol_h2, half_gap = half_gap,
half_ratio = half_ratio, stability_note = stability_note,
indep_clause = indep_clause, cluster_clause = cluster_clause,
fat_clause = fat_clause,
volatility_series_df = volatility_series_df,
drawdown_series_df = drawdown_series_df,
risk_measures_df = risk_measures_df,
tail_diagnostics_df = tail_diagnostics_df,
squared_return_acf_df = squared_return_acf_df,
model_fits_df = model_fits_df, methods_df = methods_df,
metrics = metrics, json_output = json_output
)
}