Executive Summary

Volatility, VaR, expected shortfall and drawdown for Close

Time Points

2416

Returns

2415

Periods Per Year

252

Annualized Volatility

51.63%

Latest Rolling Volatility

79.21%

Historical VaR 95%

4.61%

Expected Shortfall 95%

6.98%

Historical VaR 99%

7.91%

Maximum Drawdown

-53.51%

Excess Kurtosis

5.92

GARCH Persistence

0.999

Over 2,416 time points and 9.60 years, Close realized an annualized volatility of 51.63% — very high — with the 21-period rolling estimate ending at 79.21%. On a one-period horizon the historical Value-at-Risk is 4.61% at 95% and 7.91% at 99%; the expected shortfall beside them — the average loss GIVEN the VaR is breached — is 6.98% and 11.14%. Read the second number, not the first: VaR is a quantile and says nothing about how far past it a loss can go. The deepest peak-to-trough fall was -53.51%, from 2017-09-18 to 2019-06-03, recovered by 2019-12-18 after 198.0 days. Volatility clusters: the Ljung-Box test on squared returns to lag 10 gives 94.62 (p < 0.001), so quiet and turbulent stretches group together and today's volatility is informative about tomorrow's. What the normal assumption costs is measured here rather than assumed: excess kurtosis is 5.92, fatter-tailed than a normal, and the historical series breached the parametric normal VaR at 95%, 90 breach(es) against 120.8 expected (a ratio of 0.75); at 99%, 34 breach(es) against 24.2 expected (a ratio of 1.41). 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. The square-root-of-time annualization is defensible here: the lag-1 autocorrelation of returns is 0.011 against a two-standard-error band of 0.040, and the 5-period variance ratio is 1.03 (z = 0.25) where independence implies 1. Two limits bound everything above. Historical Value-at-Risk cannot see a loss larger than the worst one in its own sample: here that was 19.33% on 2012-01-13, and 9.60 years of history has never contained a rarer event than the rarest in it. And past volatility is not future volatility. Split this sample in half at 2015-04-16: the first half realized 56.37% annualized and the second half 46.44%, a change of 9.93 percentage points (a factor of 0.82). 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.

What this means— written from the computed numbers · yours to edit

Tesla stock realized annualized volatility of 51.63% over 9.60 years, with the 21-period rolling estimate ending at 79.21%—very high throughout. One-period historical Value-at-Risk is 4.61% at 95% confidence and 7.91% at 99%; expected shortfall (average loss when breached) is 6.98% and 11.14%. The deepest drawdown was -53.51%, from 2017-09-18 to trough on 2019-06-03, recovering by 2019-12-18 after 198 days. Volatility clusters strongly: the Ljung-Box test on squared returns to lag 10 yields 94.62 (p < 0.001), so quiet and turbulent stretches group together. The normal distribution assumption understates tail risk: excess kurtosis is 5.92, and the 99% parametric VaR was breached 34 times against 24.2 expected (ratio 1.41). The worst single-period loss was 19.33% on 2012-01-13. Split-sample analysis shows the first half (2010–2015) realized 56.37% volatility while the second half (2015–2020) realized 46.44%—a shift of 9.93 percentage points that would have made a forecast at the midpoint substantially wrong.

Overview

Analysis Overview

Volatility, Value-at-Risk and drawdown for Close across 2,416 time points.

N Points2416

N Returns2415

Periods Per Year252

Rolling Window21

Series Typeprice

What this means— written from the computed numbers · yours to edit

Tesla's Close price exhibits very high volatility of 51.63% annualized over 9.60 years (2,416 observations). The series was correctly identified as a price series—all 2,416 values are large and positive, averaging 186.404. The analysis uses log returns for volatility estimation (which add across periods for annualization) and simple returns for drawdown and risk measures (which describe percentage capital loss). The observation rate of 251.7 per year was snapped to the standard convention of 252 for annualization. Volatility estimates range from 16.82% to 122.39% across 21-period rolling windows, confirming that a single volatility figure masks substantial regime shifts.

Data Preparation

How the series was assembled, and what was dropped or averaged.

Initial Rows2416

Final Rows2416

Rows Removed0

Rows Dropped Incomplete0

Time Points2416

N Returns2415

What this means— written from the computed numbers · yours to edit

The short answer

All 2,416 observations loaded without missing data, yielding 2,415 returns for analysis. The series is not resampled to a regular calendar—returns spanning data gaps are included as-is, which inflates measured volatility compared to single-period returns. This upward bias affects every volatility figure downstream.

The detail

The dataset contains 2,416 distinct time points with no rows dropped for incomplete dates or values. The observation rate of 251.7 per year was used to annualize volatility estimates. Because the series preserves gaps in the raw data rather than resampling to a regular grid, any return bridging a gap captures more than one period's movement, mechanically raising its magnitude and thus raising the volatility estimate above what single-period returns would show.

What this can't tell you

The magnitude of the upward bias from gap-spanning returns cannot be quantified without knowing the frequency and duration of gaps in the trading calendar. Consider exporting the date gaps explicitly to measure their contribution to the volatility figure.

Visualization

Rolling Volatility

Annualized volatility of Close over time, by three estimators.

What this means— written from the computed numbers · yours to edit

Volatility moved dramatically across the sample, with the rolling 21-period estimate ranging from 16.82% to 122.39%—a spread of 105.58 percentage points. The series opened around 106.7% in July 2010, fell to a low of 25.96% in November 2010, then spiked to 84.35% in late 2010. It ended the sample at 79.21%. The EWMA, fitted with decay 0.984, smooths these moves rather than stepping abruptly. The GARCH(1,1) conditional volatility, pulled toward a long-run level of 75.12% at persistence 0.999, traces a middle path. The 105.58 percentage-point spread between rolling min and max is the clearest evidence that volatility is not stable—a single number describes an average across regimes, not any one state.

Visualization

Drawdown

How far the series fell below its own running peak, and for how long.

What this means— written from the computed numbers · yours to edit

The underwater curve shows how far below its running peak Tesla stood at every point. The worst drawdown was -53.51%, reached on 2019-06-03 after peaking on 2017-09-18—623 days from peak to trough. Recovery to the prior peak occurred on 2019-12-18, 198 days after the trough, for a full round-trip of 821 days. Drawdown is the one risk measure holders actually experience on the path; a 53.51% loss requires a 115.12% gain to break even. Early in the sample (July 2010) there was a -27.17% drawdown, and the series experienced multiple recoveries and new peaks throughout the window.

Data Table

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	4.606	percent loss	one period	The loss exceeded on 5.0% of periods in this sample.
Value-at-Risk 99%	Historical	7.909	percent loss	one period	The loss exceeded on 1.0% of periods in this sample.
Value-at-Risk 95%	Parametric (normal)	5.191	percent loss	one period	What a normal distribution with this mean and standard deviation predicts.
Value-at-Risk 99%	Parametric (normal)	7.424	percent loss	one period	What a normal distribution with this mean and standard deviation predicts.
Expected shortfall 95%	Historical	6.979	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	11.14	percent loss	one period	The AVERAGE loss on the 1.0% of periods that breached VaR.
Expected shortfall 95%	Parametric (normal)	6.56	percent loss	one period	Normal-theory expected shortfall for the same fitted distribution.
Expected shortfall 99%	Parametric (normal)	8.534	percent loss	one period	Normal-theory expected shortfall for the same fitted distribution.
Worst single-period loss	Observed	19.33	percent loss	one period	The largest single-period fall in this sample, on 2012-01-13.
Maximum drawdown	Observed	53.51	percent loss	peak to trough	Peak 2017-09-18 to trough 2019-06-03; recovered by 2019-12-18.

What this means— written from the computed numbers · yours to edit

One-period losses are quoted as percentages of capital. The historical 95% Value-at-Risk of 4.61% means 5% of periods lost more; the expected shortfall of 6.98% is the average loss when that threshold is breached—a gap of 2.37 percentage points that VaR alone hides. At 99%, the historical VaR is 7.91% and expected shortfall is 11.14%—a gap of 3.23 percentage points. The parametric normal VaR at 99% is 7.42%, sitting 0.48 percentage points below the historical figure. The worst single-period loss in the sample was 19.33% on 2012-01-13, which is 2.44 times the historical 99% VaR—a clear statement that VaR is a threshold, not a bound.

Data Table

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	5.916	0	dimensionless	fatter-tailed than normal by 5.92
Skewness of returns	0.1969	0	dimensionless	roughly symmetric
Jarque-Bera normality statistic	3538	0	chi-square, 2 df	normality rejected (p < 0.001)
Breaches of the normal 95% VaR	90	120.8	periods	90 observed against 120.8 expected, a ratio of 0.75
Breaches of the normal 99% VaR	34	24.15	periods	34 observed against 24.2 expected, a ratio of 1.41
Historical minus normal VaR 95%	-0.5849	0	percentage points	the normal model does not understate the historical loss
Historical minus normal VaR 99%	0.4849	0	percentage points	the normal model understates the historical loss

What this means— written from the computed numbers · yours to edit

The short answer

Tesla returns are substantially fatter-tailed than a normal distribution: excess kurtosis of 5.92 means the tails hold far more probability than the normal model predicts. The Jarque-Bera test rejects normality decisively (p < 0.001). At the 99% confidence level, extreme losses occur 1.41 times as often as the normal VaR model expects (34 observed breaches versus 24.2 expected).

The detail

Excess kurtosis is 5.9164 (normal = 0), skewness is 0.1969 (roughly symmetric), and the Jarque-Bera statistic is 3537.82 on 2 degrees of freedom (p < 0.001). The 99% VaR breaches show 34 observed against 24.15 expected, a ratio of 1.41. At 95%, the ratio is 0.75 (90 observed against 120.75 expected). The normal model understates the historical 99% loss by 0.4849 percentage points. With a sampling standard deviation of about 4.9 around the 24.2 expected breaches, the 99% ratio near 1.41 reflects genuine tail risk, not noise.

What this can't tell you

The breach counts are small numbers; ratios near 1 should not be over-read. The 95% under-prediction and 99% over-prediction together suggest the normal model's symmetric widening masks the actual tail shape. Historical quantiles are more reliable than the parametric model for this series.

Visualization

Volatility Clustering

Autocorrelation of squared returns, with a Ljung-Box test.

What this means— written from the computed numbers · yours to edit

Volatility clusters: turbulent and quiet periods arrive together. Four of 10 lags in the autocorrelation of squared returns fall outside the 0.040 two-standard-error band. The tallest bar is lag 1 at 0.1623, and the Ljung-Box test on squared returns to lag 10 gives 94.62 (p < 0.001). This means today's volatility is informative about tomorrow's. A single unconditional volatility number (51.63%) understates risk during turbulent stretches and overstates it during calm ones—precisely why the EWMA and GARCH conditional models matter. The test does not distinguish between true clustering, isolated extreme observations, or a shift in variance level partway through the sample.

Data Table

Volatility Models

EWMA and GARCH(1,1), both fitted by direct maximum likelihood.

Parameter	Estimate	Display
EWMA decay (lambda, fitted by MLE)	0.9842	0.984
EWMA latest annualized volatility	52.36	52.36%
RiskMetrics EWMA (lambda 0.94) latest	59.58	59.58%
GARCH omega	2.52e-06	0.00000252
GARCH alpha (news impact)	0.0193	0.019
GARCH beta (persistence of past variance)	0.9796	0.980
GARCH alpha + beta	0.9989	0.999
GARCH shock half-life (periods)	614.7	614.7
GARCH long-run annualized volatility	75.12	75.12%
GARCH latest conditional volatility	54.47	54.47%
GARCH one-period-ahead forecast	67.01	67.01%
Likelihood ratio, GARCH against constant variance	154.8	154.81 (p < 0.001)

What this means— written from the computed numbers · yours to edit

Two conditional-variance models were fitted by maximum likelihood. The EWMA decay was fitted at 0.984 (interior to the search range, so a genuine estimate rather than a boundary), compared to the RiskMetrics convention of 0.94. The fitted EWMA puts latest annualized volatility at 52.36%. The GARCH(1,1) found alpha 0.019 (news impact) and beta 0.980 (persistence of past variance), summing to 0.999. This persistence implies shocks decay by half in 614.7 periods and produce a long-run annualized volatility of 75.12% against the realized 51.63%. Latest conditional volatility is 54.47% and the one-period-ahead forecast is 67.01%. The GARCH improves on constant variance by a likelihood-ratio statistic of 154.81 on 2 degrees of freedom (p < 0.001), confirming clustering is real. Persistence near 1 makes the long-run figure unstable; a conditional volatility forecast describes how wide the next period will be, never which direction returns will move.

Data Table

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 Close values are smaller than 1 in absolute size, their average absolute size is 186.404, 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 Close 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. 2,415 return(s) were formed from 2,416 time points.
Annualization	The series carries 2,416 observations across 9.60 years, an observed rate of 251.7 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 16.82% to 122.39% and ended at 79.21%.
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.198, standard deviation 3.276 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 2017-09-18, trough 2019-06-03, depth 53.51%, recovered 2019-12-18 after 198.0 period-days.
Clustering test	Ljung-Box test on the squared returns to lag 10: statistic 94.62, p < 0.001. 4 of the 10 lag autocorrelations sit outside the two-standard-error band of 0.040.
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.984; the RiskMetrics convention of 0.94 is reported beside it for comparison. The fitted decay of 0.984 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 154.81 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 19.33%, which is 2.44 times the historical 99% VaR of 7.91%.
What history cannot see	Historical VaR cannot see a loss larger than the worst one in its sample. This sample covers 9.60 year(s), from 2010-06-29 to 2020-02-03, and the worst single period in it lost 19.33%. A window of that length has never observed an event rarer than roughly one in 2,415.
What the normal assumption costs	Excess kurtosis of the returns is 5.92, and the historical series breached the normal 99% VaR 34 time(s) against the 24.2 the normal model expects. The historical 99% VaR sits 0.48 percentage points above the parametric one.
Forecast limit	Split this sample in half at 2015-04-16: the first half realized 56.37% annualized and the second half 46.44%, a change of 9.93 percentage points (a factor of 0.82). 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.

What this means— written from the computed numbers · yours to edit

Log returns are used for volatility, EWMA, and GARCH because they add across periods (required for square-root-of-252 annualization); simple returns are used for Value-at-Risk, expected shortfall, and drawdown because they describe percentage capital loss. Historical VaR is the empirical quantile—no distribution assumed—and is bounded by the worst loss in the sample (19.33% on 2012-01-13 across 9.60 years). Parametric VaR assumes normality; the tail diagnostics (excess kurtosis 5.92, 34 breaches at 99% against 24.2 expected) show the cost. Expected shortfall is the average loss given VaR is breached. Square-root-of-252 annualization assumes independence: lag-1 autocorrelation of returns is 0.011 (within the 0.040 two-standard-error band) and the 5-period variance ratio is 1.03 (z = 0.25), both supporting independence. The hardest limit: split the sample at 2015-04-16, the first half realized 56.37% volatility and the second 46.44%—a difference of 9.93 percentage points (factor of 0.82)—showing past volatility is not future volatility.

Rate this report Was this the answer you needed?

The exact source that produced this report — yours to keep, read, and re-run.

Download PDF

How this was computed method · R source · citation

The code that did it

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 &#x27;%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 &#x27;%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 &#x27;%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 &#x27;%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$v

Step 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 &#x27;%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 &#x27;%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 &#x27;%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 &#x27;%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 &#x27;%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 &#x27;%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_factor

A 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 &#x27;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) <- NULL

Step 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&#x27;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) <- NULL

Step 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&#x27;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
  )
}

Your data has more stories to tell.Run any analysis on your own data — validated R modules, interactive reports, AI insights, and PDF export. 500 free credits on signup.

Try Free — No Signup Sign Up Free

Download PDF · Print · Share ·Refresh this analysis ·Edit this report

What's wrong with this card?

Executive Summary

Analysis Overview

Data Preparation

The short answer

The detail

What this can't tell you

Rolling Volatility

Drawdown

Value-at-Risk, Expected Shortfall and Drawdown

What The Normal Assumption Costs

The short answer

The detail

What this can't tell you

Volatility Clustering

Volatility Models

Method & Limits

Volatility & Value-at-Risk — How Risky Is This Series?

Why This Method?

What This Analysis Covers

Standard Library

Core Analysis Pipeline

Step 1: Resolve the mapped columns (humanized everywhere in prose)

Step 2: Coerce the value to numeric (95% rule)

Step 3: Parse the date column (ISO, ymd, mdy, dmy, or integer years)

Step 4: Keep complete rows; average several rows on the same date

Step 5: Is this a level series or is it already returns?

Step 6: Periods per year — observed, then snapped to a convention

Step 7: Rolling volatility

Step 8: Value-at-Risk and expected shortfall, historical and normal

Step 9: What the normal assumption actually costs here

Step 10: What the historical method cannot see

Step 11: Drawdown on the wealth index

Step 12: Volatility clustering

Step 13: Is the square-root-of-time rule safe here?

Step 14: EWMA and GARCH(1,1), both fitted by maximum likelihood

Step 15: Past volatility is not future volatility — measured, not said

Step 16: Chart datasets

Step 17: Result tables

Step 18: Headline metrics and the computed answer

Cite this analysis

Report an Issue