import pandas as pd
import numpy as np
from scipy import stats
from scipy.optimize import minimize
import plotnine as p9
from mizani.formatters import percent_format, comma_format
import warnings
warnings.filterwarnings("ignore")41 Tail Risk and Extreme Events
Finance is fundamentally a discipline of extremes. The central objects of concern (e.g., asset pricing, portfolio allocation, capital regulation, systemic stability) are all disproportionately shaped by rare but severe events. A portfolio manager who correctly estimates the mean and variance of returns but ignores the possibility of a \(-20\%\) daily move is, in the language of Taleb (2010), “picking up pennies in front of a steamroller.” The 1997 Asian Financial Crisis, the 2008 Global Financial Crisis, and the 2020 COVID-19 crash each demonstrated that the most consequential realizations in financial markets are precisely those that conventional Gaussian models assign near-zero probability.
This chapter develops the statistical and econometric toolkit for measuring, modeling, and diagnosing tail risk in Vietnamese equity markets. We focus on three interconnected layers. First, at the individual stock level, we estimate crash risk measures that capture the asymmetry and thickness of the left tail of return distributions. Second, at the portfolio or regulatory level, we implement Value at Risk (VaR) and Expected Shortfall (ES) under multiple estimation approaches and evaluate their accuracy via backtesting. Third, at the system level, we apply Extreme Value Theory (EVT) and tail dependence measures to characterize joint crash behavior across sectors and assess financial contagion.
Vietnamese markets present a particularly rich laboratory for studying tail phenomena. Daily price limits mechanically censor the tails of observed return distributions, creating latent tail risk that is not directly observable but must be inferred. State ownership concentrations, thin liquidity, and herding behavior among retail investors amplify crash dynamics. The absence of a developed derivatives market means that tail hedging instruments familiar in more developed markets (e.g., deep out-of-the-money puts, variance swaps) are largely unavailable, making accurate tail risk measurement all the more important for risk management.
41.1 Crash Risk in Equity Markets
41.1.1 The Left Tail Problem
The return distribution of individual stocks and market indices is not symmetric. Decades of evidence, beginning with Mandelbrot et al. (1963) and formalized in the asset pricing context by Harvey and Siddique (2000), establish that equity returns exhibit negative skewness (i.e., large negative returns occur more frequently than a Gaussian model predicts) and excess kurtosis (i.e., the tails are thicker than normal in both directions, but the left tail is of primary economic concern).
Let \(r_{i,t}\) denote the daily log return of stock \(i\) on day \(t\). The standardized third central moment (skewness) is:
\[ \text{Skew}(r_i) = \frac{E\left[(r_{i,t} - \mu_i)^3\right]}{\sigma_i^3} \tag{41.1}\]
Negative skewness implies that the distribution has a longer or fatter left tail relative to the right. In economic terms, negative skewness means that large losses are more likely than large gains of equal magnitude. The standardized fourth moment (excess kurtosis) is:
\[ \text{Kurt}(r_i) = \frac{E\left[(r_{i,t} - \mu_i)^4\right]}{\sigma_i^4} - 3 \tag{41.2}\]
Positive excess kurtosis indicates heavier tails than the normal distribution. Both moments have direct asset pricing implications: Harvey and Siddique (2000) and Kraus and Litzenberger (1976) argue that investors demand compensation for holding negatively skewed assets, implying that crash-prone stocks should earn higher expected returns, ceteris paribus.
41.1.2 Crash Risk Measures
The literature has developed several firm-specific crash risk measures. We implement the three most widely used, following Chen, Hong, and Stein (2001) and Hutton, Marcus, and Tehranian (2009).
Measure 1: Negative Coefficient of Skewness (NCSKEW)
For each firm \(i\) in fiscal year \(y\), compute firm-specific weekly returns \(W_{i,\tau}\) as residuals from a market model augmented with lead and lag market returns to account for nonsynchronous trading:
\[ r_{i,t} = \alpha_i + \beta_{1i} r_{m,t-1} + \beta_{2i} r_{m,t} + \beta_{3i} r_{m,t+1} + \varepsilon_{i,t} \tag{41.3}\]
where \(r_{m,t}\) is the value-weighted market return on day \(t\). The firm-specific weekly return is \(W_{i,\tau} = \ln(1 + \sum_{t \in \tau} \hat{\varepsilon}_{i,t})\), where \(\tau\) indexes weeks. NCSKEW is then:
\[ \text{NCSKEW}_{i,y} = -\frac{n(n-1)^{3/2} \sum W_{i,\tau}^3}{(n-1)(n-2)\left(\sum W_{i,\tau}^2\right)^{3/2}} \tag{41.4}\]
where \(n\) is the number of firm-specific weekly returns in the year. The negative sign ensures that higher NCSKEW corresponds to greater crash risk.
Measure 2: Down-to-Up Volatility (DUVOL)
Partition the firm-specific weekly returns into “down weeks” (\(W_{i,\tau} < \bar{W}_i\)) and “up weeks” (\(W_{i,\tau} \geq \bar{W}_i\)), where \(\bar{W}_i\) is the annual mean. Then:
\[ \text{DUVOL}_{i,y} = \ln\left(\frac{(n_u - 1)\sum_{\text{down}} W_{i,\tau}^2}{(n_d - 1)\sum_{\text{up}} W_{i,\tau}^2}\right) \tag{41.5}\]
where \(n_u\) and \(n_d\) are the number of up and down weeks, respectively. Higher DUVOL indicates greater crash risk. Chen, Hong, and Stein (2001) argue that DUVOL is less sensitive to outliers than NCSKEW because it does not involve the cube of returns.
Measure 3: Crash Count (COUNT)
Define a crash event as a firm-specific weekly return that falls below \(k\) standard deviations of its annual distribution:
\[ \text{CRASH}_{i,\tau} = \mathbb{1}\left(W_{i,\tau} < \bar{W}_i - k \cdot \hat{\sigma}_{W_i}\right) \tag{41.6}\]
The standard threshold in the literature is \(k = 3.09\), corresponding to the 0.1% tail of a standard normal distribution. The crash count for year \(y\) is \(\text{COUNT}_{i,y} = \sum_{\tau \in y} \text{CRASH}_{i,\tau}\).
# DataCore.vn API
from datacore import DataCore
dc = DataCore()
# Load daily stock returns and market returns
daily_returns = dc.get_daily_returns(
start_date="2010-01-01",
end_date="2024-12-31"
)
# Market return: value-weighted index
market_returns = dc.get_market_returns(
start_date="2010-01-01",
end_date="2024-12-31",
frequency="daily"
)
# Merge
df = daily_returns.merge(
market_returns[["date", "mkt_ret"]],
on="date",
how="left"
)
print(f"Universe: {df['ticker'].nunique()} firms, "
f"{df['date'].nunique()} trading days")from sklearn.linear_model import LinearRegression
def compute_firm_specific_returns(group):
"""
Estimate firm-specific daily returns from augmented market model
with lead and lag market returns (Chen, Hong, Stein 2001).
"""
g = group.sort_values("date").copy()
g["mkt_lag1"] = g["mkt_ret"].shift(1)
g["mkt_lead1"] = g["mkt_ret"].shift(-1)
g = g.dropna(subset=["ret", "mkt_ret", "mkt_lag1", "mkt_lead1"])
if len(g) < 30:
return pd.DataFrame()
X = g[["mkt_lag1", "mkt_ret", "mkt_lead1"]].values
y = g["ret"].values
model = LinearRegression().fit(X, y)
g["resid"] = y - model.predict(X)
return g[["ticker", "date", "resid"]]
# Estimate by firm-year
df["year"] = df["date"].dt.year
firm_resids = (
df.groupby(["ticker", "year"], group_keys=False)
.apply(compute_firm_specific_returns)
.reset_index(drop=True)
)# Aggregate daily residuals to weekly firm-specific returns
firm_resids["date"] = pd.to_datetime(firm_resids["date"])
firm_resids["week"] = firm_resids["date"].dt.isocalendar().week.astype(int)
firm_resids["year"] = firm_resids["date"].dt.year
weekly_returns = (
firm_resids.groupby(["ticker", "year", "week"])
.agg(W=("resid", lambda x: np.log(1 + x.sum())))
.reset_index()
)def compute_crash_measures(group):
"""Compute NCSKEW, DUVOL, and CRASH count for one firm-year."""
W = group["W"].values
n = len(W)
if n < 26: # Require at least 26 weeks
return pd.Series({
"ncskew": np.nan, "duvol": np.nan,
"crash_count": np.nan, "n_weeks": n
})
# NCSKEW
W_demean = W - W.mean()
num = n * (n - 1)**1.5 * np.sum(W_demean**3)
den = (n - 1) * (n - 2) * (np.sum(W_demean**2))**1.5
ncskew = -(num / den) if den != 0 else np.nan
# DUVOL
mean_W = W.mean()
down = W[W < mean_W]
up = W[W >= mean_W]
n_d, n_u = len(down), len(up)
if n_d > 1 and n_u > 1:
duvol = np.log(
((n_u - 1) * np.sum(down**2)) /
((n_d - 1) * np.sum(up**2))
)
else:
duvol = np.nan
# Crash count (k = 3.09)
threshold = W.mean() - 3.09 * W.std()
crash_count = int(np.sum(W < threshold))
return pd.Series({
"ncskew": ncskew,
"duvol": duvol,
"crash_count": crash_count,
"n_weeks": n
})
crash_risk = (
weekly_returns.groupby(["ticker", "year"])
.apply(compute_crash_measures)
.reset_index()
)
crash_risk = crash_risk.dropna(subset=["ncskew", "duvol"])
print(f"Crash risk panel: {len(crash_risk)} firm-years")41.1.3 Cross-Sectional Distribution of Crash Risk
Table 41.1 presents summary statistics for the three crash risk measures across the Vietnamese equity market.
summary = crash_risk[["ncskew", "duvol", "crash_count"]].describe(
percentiles=[0.05, 0.25, 0.5, 0.75, 0.95]
).T.round(4)
summary.columns = ["N", "Mean", "Std", "Min", "5%", "25%",
"Median", "75%", "95%", "Max"]
summaryplot_data = crash_risk[crash_risk["year"].between(2012, 2024)].copy()
(
p9.ggplot(plot_data, p9.aes(x="ncskew"))
+ p9.geom_histogram(bins=50, fill="#2E5090", alpha=0.7)
+ p9.facet_wrap("~year", scales="free_y", ncol=4)
+ p9.geom_vline(xintercept=0, linetype="dashed", color="red", size=0.5)
+ p9.labs(
x="NCSKEW",
y="Count",
title="Distribution of Firm-Level Crash Risk (NCSKEW)"
)
+ p9.theme_minimal()
+ p9.theme(figure_size=(12, 8))
)The rightward shift of the NCSKEW distribution during crisis episodes, particularly in 2020 and 2022, indicates a market-wide increase in crash risk. The red dashed line at zero separates firms with negative skewness (right of zero, i.e., high NCSKEW since the measure is negated) from those with positive skewness.
41.1.4 Interpretation in Asset Pricing
Crash risk has first-order implications for asset pricing. Harvey and Siddique (2000) demonstrate that coskewness (i.e., the contribution of an individual asset to the skewness of the aggregate market portfolio) is priced in the cross-section of expected returns. Specifically, assets that tend to crash when the market crashes (negative coskewness) should command a risk premium.
The coskewness of stock \(i\) with the market is:
\[ \text{CoSkew}_i = \frac{E\left[(r_i - \mu_i)(r_m - \mu_m)^2\right]}{\sigma_i \cdot \sigma_m^2} \tag{41.7}\]
Stocks with more negative coskewness exhibit disproportionately poor performance during market downturns. In Vietnam, this is especially relevant because: (1) retail investor herding amplifies co-movement during selloffs; (2) price limits delay crash completion, spreading crash risk across multiple days; and (3) state-owned enterprises, which constitute a large fraction of market capitalization, may exhibit coordinated crash behavior driven by common policy shocks.
def compute_coskewness(group):
"""Estimate coskewness for a firm-year."""
r_i = group["ret"].values
r_m = group["mkt_ret"].values
if len(r_i) < 60:
return np.nan
r_i_dm = r_i - r_i.mean()
r_m_dm = r_m - r_m.mean()
coskew = (
np.mean(r_i_dm * r_m_dm**2) /
(np.std(r_i) * np.std(r_m)**2)
)
return coskew
df_coskew = (
df.groupby(["ticker", "year"])
.apply(lambda g: pd.Series({"coskewness": compute_coskewness(g)}))
.reset_index()
.dropna()
)median_coskew = (
df_coskew.groupby("year")
.agg(
median_coskew=("coskewness", "median"),
q25=("coskewness", lambda x: x.quantile(0.25)),
q75=("coskewness", lambda x: x.quantile(0.75))
)
.reset_index()
)
(
p9.ggplot(median_coskew, p9.aes(x="year"))
+ p9.geom_ribbon(
p9.aes(ymin="q25", ymax="q75"),
fill="#2E5090", alpha=0.2
)
+ p9.geom_line(p9.aes(y="median_coskew"), color="#2E5090", size=1)
+ p9.geom_hline(yintercept=0, linetype="dashed", color="gray")
+ p9.labs(
x="Year",
y="Coskewness",
title="Cross-Sectional Coskewness: Median and Interquartile Range"
)
+ p9.theme_minimal()
+ p9.theme(figure_size=(10, 5))
)41.1.5 Downside Tail Concentration
Beyond skewness, the concentration of returns in the extreme left tail provides additional diagnostic information. We define the downside tail concentration ratio as the fraction of total return variance attributable to observations below the \(p\)-th percentile:
\[ \text{TCR}_p = \frac{\sum_{r_{i,t} < q_p} (r_{i,t} - \bar{r}_i)^2}{\sum_t (r_{i,t} - \bar{r}_i)^2} \tag{41.8}\]
where \(q_p\) is the \(p\)-th percentile of the firm’s return distribution. Under a symmetric distribution, \(\text{TCR}_{5\%}\) would be close to \(5\%\) (with a slight upward bias due to the squaring). Values substantially exceeding \(5\%\) indicate that the left tail contributes disproportionately to total risk.
def tail_concentration_ratio(returns, p=0.05):
"""Compute tail concentration ratio at percentile p."""
q = np.quantile(returns, p)
r_dm = returns - returns.mean()
total_var = np.sum(r_dm**2)
tail_var = np.sum(r_dm[returns < q]**2)
return tail_var / total_var if total_var > 0 else np.nan
tcr = (
df.groupby(["ticker", "year"])
.agg(
tcr_5=("ret", lambda x: tail_concentration_ratio(x.values, 0.05)),
tcr_1=("ret", lambda x: tail_concentration_ratio(x.values, 0.01)),
n_obs=("ret", "count")
)
.reset_index()
.query("n_obs >= 120")
)
print("Median TCR(5%):", round(tcr["tcr_5"].median(), 4))
print("Median TCR(1%):", round(tcr["tcr_1"].median(), 4))41.2 Value at Risk and Expected Shortfall
41.2.1 Definitions
Value at Risk (VaR) and Expected Shortfall (ES) are the two principal quantile-based risk measures used in financial regulation and internal risk management. For a portfolio return \(r\) with the cumulative distribution function \(F\), the VaR at confidence level \(\alpha\) is:
\[ \text{VaR}_\alpha = -F^{-1}(\alpha) = -\inf\{x : F(x) \geq \alpha\} \tag{41.9}\]
This is simply the negative of the \(\alpha\)-quantile of the return distribution. For \(\alpha = 0.01\), \(\text{VaR}_{1\%}\) answers: “What is the loss that is exceeded with only 1% probability?”
Expected Shortfall (also called Conditional VaR or CVaR) is the expected loss conditional on the loss exceeding VaR:
\[ \text{ES}_\alpha = -\frac{1}{\alpha}\int_0^\alpha F^{-1}(u) \, du = -E\left[r \mid r \leq -\text{VaR}_\alpha\right] \tag{41.10}\]
For continuous distributions, ES simplifies to the conditional expectation in Equation 41.10. ES is always at least as large as VaR (\(\text{ES}_\alpha \geq \text{VaR}_\alpha\)) and is strictly larger whenever the distribution has any mass beyond the VaR quantile.
41.2.2 Why Expected Shortfall Dominates VaR
Table 41.2 summarizes the fundamental differences.
| Property | VaR | Expected Shortfall |
|---|---|---|
| Definition | Quantile of the loss distribution | Conditional mean loss beyond VaR |
| Tail information | None beyond the quantile | Captures severity of tail losses |
| Subadditivity | Not subadditive in general | Subadditive (coherent risk measure) |
| Regulatory status | Basel II primary measure | Basel III/IV primary measure |
| Backtesting | Binary: breach or no breach | Continuous: requires severity assessment |
| Sensitivity to tail shape | None | Directly sensitive |
The critical deficiency of VaR is the absence of subadditivity. A risk measure \(\rho\) is subadditive if \(\rho(X + Y) \leq \rho(X) + \rho(Y)\) for all portfolios \(X, Y\). VaR violates this condition for non-elliptical distributions, meaning that diversification can appear to increase risk under VaR, which is a perverse result. Artzner et al. (1999) established the axiomatic framework for coherent risk measures, and ES satisfies all four axioms: monotonicity, translation invariance, positive homogeneity, and subadditivity. VaR satisfies only three.
41.2.3 Estimation Methods
We implement four approaches to VaR and ES estimation, each with distinct assumptions and data requirements.
Method 1: Historical Simulation
The simplest nonparametric approach uses the empirical distribution directly:
\[ \widehat{\text{VaR}}_\alpha^{\text{HS}} = -\hat{F}_n^{-1}(\alpha) = -r_{(\lceil n\alpha \rceil)} \tag{41.11}\]
where \(r_{(k)}\) is the \(k\)-th order statistic of the return sample. ES is the average of all returns below the VaR:
\[ \widehat{\text{ES}}_\alpha^{\text{HS}} = -\frac{1}{\lfloor n\alpha \rfloor}\sum_{k=1}^{\lfloor n\alpha \rfloor} r_{(k)} \tag{41.12}\]
Method 2: Parametric (Gaussian)
Assume \(r_t \sim N(\mu, \sigma^2)\). Then:
\[ \text{VaR}_\alpha^{\text{Gauss}} = -(\hat{\mu} + \hat{\sigma} \cdot \Phi^{-1}(\alpha)) \tag{41.13}\]
\[ \text{ES}_\alpha^{\text{Gauss}} = -\hat{\mu} + \hat{\sigma} \cdot \frac{\phi(\Phi^{-1}(\alpha))}{\alpha} \tag{41.14}\]
where \(\Phi\) and \(\phi\) are the standard normal CDF and PDF, respectively. The Gaussian assumption is known to underestimate tail risk for equity returns, but provides a useful baseline.
Method 3: Parametric (Student-\(t\))
The Student-\(t\) distribution with \(\nu\) degrees of freedom captures excess kurtosis. The VaR and ES are:
\[ \text{VaR}_\alpha^{t} = -\left(\hat{\mu} + \hat{\sigma}\sqrt{\frac{\nu - 2}{\nu}} \cdot t_\nu^{-1}(\alpha)\right) \tag{41.15}\]
\[ \text{ES}_\alpha^{t} = -\hat{\mu} + \hat{\sigma}\sqrt{\frac{\nu-2}{\nu}} \cdot \frac{f_{t_\nu}(t_\nu^{-1}(\alpha))}{\alpha} \cdot \frac{\nu + (t_\nu^{-1}(\alpha))^2}{\nu - 1} \tag{41.16}\]
where \(t_\nu^{-1}\) and \(f_{t_\nu}\) are the quantile function and PDF of the Student-\(t\) distribution.
Method 4: GARCH-Filtered EVT (Semi-Parametric)
This approach, developed by McNeil and Frey (2000), first filters returns through a GARCH(1,1) model to obtain standardized residuals \(z_t = \hat{\varepsilon}_t / \hat{\sigma}_t\), then fits a Generalized Pareto Distribution (GPD) to the lower tail of the standardized residuals. We detail the EVT component in Section 41.3.
def estimate_var_es(returns, alpha=0.01):
"""
Estimate VaR and ES at level alpha using four methods.
Parameters
----------
returns : array-like
Return series.
alpha : float
Tail probability (e.g., 0.01 for 1%).
Returns
-------
dict : VaR and ES estimates for each method.
"""
r = np.array(returns)
r = r[~np.isnan(r)]
n = len(r)
mu, sigma = r.mean(), r.std(ddof=1)
results = {}
# Method 1: Historical simulation
sorted_r = np.sort(r)
idx = int(np.ceil(n * alpha))
results["var_hs"] = -sorted_r[idx - 1]
results["es_hs"] = -sorted_r[:idx].mean()
# Method 2: Gaussian
z_alpha = stats.norm.ppf(alpha)
results["var_gauss"] = -(mu + sigma * z_alpha)
results["es_gauss"] = -mu + sigma * stats.norm.pdf(z_alpha) / alpha
# Method 3: Student-t (MLE for degrees of freedom)
try:
nu, loc, scale = stats.t.fit(r)
t_alpha = stats.t.ppf(alpha, df=nu)
results["var_t"] = -(loc + scale * t_alpha)
results["es_t"] = (
-loc + scale *
(stats.t.pdf(t_alpha, df=nu) / alpha) *
((nu + t_alpha**2) / (nu - 1))
)
results["nu_hat"] = nu
except Exception:
results["var_t"] = np.nan
results["es_t"] = np.nan
results["nu_hat"] = np.nan
return results# Compute VaR and ES for the aggregate market index
market_daily = dc.get_market_returns(
start_date="2010-01-01",
end_date="2024-12-31",
frequency="daily"
)
# Rolling 1-year VaR/ES
window = 252
results_list = []
for end_idx in range(window, len(market_daily)):
window_returns = market_daily["mkt_ret"].iloc[end_idx - window:end_idx].values
date = market_daily["date"].iloc[end_idx]
est = estimate_var_es(window_returns, alpha=0.01)
est["date"] = date
results_list.append(est)
var_es_ts = pd.DataFrame(results_list)plot_var = var_es_ts.melt(
id_vars="date",
value_vars=["var_hs", "var_gauss", "var_t"],
var_name="method",
value_name="var"
)
method_labels = {
"var_hs": "Historical Simulation",
"var_gauss": "Gaussian",
"var_t": "Student-t"
}
plot_var["method"] = plot_var["method"].map(method_labels)
(
p9.ggplot(plot_var, p9.aes(x="date", y="var", color="method"))
+ p9.geom_line(alpha=0.8, size=0.6)
+ p9.labs(
x="", y="1% VaR (Loss)",
title="Rolling 252-Day Value at Risk Estimates",
color="Method"
)
+ p9.scale_color_manual(
values=["#2E5090", "#C0392B", "#27AE60"]
)
+ p9.theme_minimal()
+ p9.theme(figure_size=(12, 5), legend_position="top")
)plot_es = var_es_ts.melt(
id_vars="date",
value_vars=["es_hs", "es_gauss", "es_t"],
var_name="method",
value_name="es"
)
method_labels_es = {
"es_hs": "Historical Simulation",
"es_gauss": "Gaussian",
"es_t": "Student-t"
}
plot_es["method"] = plot_es["method"].map(method_labels_es)
(
p9.ggplot(plot_es, p9.aes(x="date", y="es", color="method"))
+ p9.geom_line(alpha=0.8, size=0.6)
+ p9.labs(
x="", y="1% ES (Loss)",
title="Rolling 252-Day Expected Shortfall Estimates",
color="Method"
)
+ p9.scale_color_manual(
values=["#2E5090", "#C0392B", "#27AE60"]
)
+ p9.theme_minimal()
+ p9.theme(figure_size=(12, 5), legend_position="top")
)41.2.4 VaR Backtesting
A VaR model is only useful if its predictions are accurate. Backtesting evaluates whether the realized frequency of VaR breaches is consistent with the nominal confidence level. If VaR is estimated at the \(\alpha = 1\%\) level, we expect approximately 1% of days to exhibit losses exceeding VaR.
The Kupiec unconditional coverage test (Kupiec et al. 1995) evaluates whether the total number of breaches \(x\) in \(n\) observations is consistent with \(\alpha\):
\[ \text{LR}_{\text{UC}} = -2\ln\left(\frac{\alpha^x (1-\alpha)^{n-x}}{\hat{p}^x (1-\hat{p})^{n-x}}\right) \sim \chi^2(1) \tag{41.17}\]
where \(\hat{p} = x/n\) is the empirical breach frequency. The Christoffersen conditional coverage test (Christoffersen 1998) additionally evaluates whether breaches are independent (no clustering):
\[ \text{LR}_{\text{CC}} = \text{LR}_{\text{UC}} + \text{LR}_{\text{IND}} \tag{41.18}\]
where \(\text{LR}_{\text{IND}}\) tests the null of independence using a first-order Markov chain model for breach indicators.
def kupiec_test(returns, var_estimates, alpha=0.01):
"""
Kupiec unconditional coverage test for VaR backtesting.
Returns
-------
dict : breach count, expected, breach rate, LR statistic, p-value.
"""
breaches = returns < -var_estimates
x = breaches.sum()
n = len(returns)
p_hat = x / n
if p_hat == 0 or p_hat == 1:
return {"breaches": x, "expected": n * alpha,
"breach_rate": p_hat, "lr_stat": np.nan, "p_value": np.nan}
lr = -2 * (
x * np.log(alpha) + (n - x) * np.log(1 - alpha)
- x * np.log(p_hat) - (n - x) * np.log(1 - p_hat)
)
p_value = 1 - stats.chi2.cdf(lr, df=1)
return {
"breaches": int(x),
"expected": round(n * alpha, 1),
"breach_rate": round(p_hat, 4),
"lr_stat": round(lr, 4),
"p_value": round(p_value, 4)
}
# Backtest each method
actual_returns = market_daily["mkt_ret"].iloc[window:].values
backtest_results = {}
for method, col in [("Historical", "var_hs"),
("Gaussian", "var_gauss"),
("Student-t", "var_t")]:
backtest_results[method] = kupiec_test(
actual_returns, var_es_ts[col].values, alpha=0.01
)
bt_df = pd.DataFrame(backtest_results).T
bt_df.index.name = "Method"
bt_df41.2.5 Estimation Under Limited Data
In Vietnamese equity markets, many stocks have short trading histories, thin liquidity, and censored returns due to price limits. These features create challenges for tail risk estimation that are less severe in developed markets.
Price limit censoring. When a stock hits its daily price limit, the observed return is truncated. The true latent return may extend well beyond the limit, but is unobservable. This means that historical simulation mechanically understates VaR and ES for stocks that frequently hit limits. A Tobit-type correction can partially address this:
\[ r_{i,t}^{\text{latent}} \sim N(\mu_i, \sigma_i^2), \quad r_{i,t}^{\text{obs}} = \max(\underline{L}, \min(\bar{L}, r_{i,t}^{\text{latent}})) \tag{41.19}\]
The parameters \((\mu_i, \sigma_i^2)\) can be estimated via maximum likelihood on the censored sample, and VaR/ES can then be computed from the estimated latent distribution.
Short histories. For recently listed stocks, the available sample may be too short for reliable nonparametric estimation. In these cases, shrinkage toward a cross-sectional prior (e.g., the sector median VaR) or Bayesian approaches with informative priors can improve stability.
41.3 Extreme Value Theory
41.3.1 Motivation: What EVT Captures That GARCH Does Not
GARCH models capture volatility clustering (i.e., the tendency of large returns to be followed by large returns), but they make specific distributional assumptions about standardized residuals (typically Gaussian or Student-\(t\)). The tail behavior of the conditional distribution is thus determined by the parametric innovation assumption, not learned from the data.
Extreme Value Theory (EVT) takes a fundamentally different approach. It provides a statistical framework for estimating the distribution of extreme observations without specifying the entire return distribution. The key theoretical results (i.e., the Fisher-Tippett-Gnedenko theorem and the Pickands-Balkema-de Haan theorem) establish that the distribution of properly normalized maxima (or exceedances over high thresholds) converges to specific parametric families regardless of the parent distribution. This universality is what makes EVT powerful: the tail can be estimated even when the center of the distribution is misspecified.
41.3.2 Tail Index Estimation
The tail index \(\xi\) (also denoted \(\alpha^{-1}\) in some formulations) governs the rate of tail decay. For heavy-tailed distributions, the survival function satisfies:
\[ P(X > x) \sim L(x) \cdot x^{-1/\xi}, \quad x \to \infty \tag{41.20}\]
where \(L(x)\) is a slowly varying function. The tail index \(\xi > 0\) implies a Pareto-type tail; larger \(\xi\) means heavier tails. Financial return distributions typically have \(\xi \in (0.2, 0.5)\), corresponding to tail indices \(\alpha = 1/\xi \in (2, 5)\), which implies finite variance but potentially infinite higher moments.
The Hill estimator (Hill 1975) provides a simple nonparametric estimate of \(\xi\) from the upper order statistics:
\[ \hat{\xi}_{\text{Hill}}(k) = \frac{1}{k}\sum_{j=1}^{k} \ln X_{(n-j+1)} - \ln X_{(n-k)} \tag{41.21}\]
where \(X_{(1)} \leq \cdots \leq X_{(n)}\) are the order statistics and \(k\) is the number of upper order statistics used. The choice of \(k\) involves a bias-variance tradeoff: small \(k\) yields high variance; large \(k\) introduces bias from non-tail observations.
def hill_estimator(data, k_values=None):
"""
Compute Hill estimator of tail index for a range of k values.
Parameters
----------
data : array-like
Positive values (use absolute values of losses).
k_values : array-like, optional
Number of upper order statistics to use.
Returns
-------
DataFrame with columns [k, xi_hat, se].
"""
x = np.sort(np.abs(data))[::-1] # Descending order
n = len(x)
if k_values is None:
k_values = range(10, min(n // 2, 500), 5)
results = []
for k in k_values:
if k >= n or k < 2:
continue
log_excess = np.log(x[:k]) - np.log(x[k])
xi_hat = log_excess.mean()
se = xi_hat / np.sqrt(k)
results.append({"k": k, "xi_hat": xi_hat, "se": se})
return pd.DataFrame(results)# Use negative market returns (losses)
losses = -market_daily["mkt_ret"].dropna().values
losses_positive = losses[losses > 0]
hill_results = hill_estimator(losses_positive)
(
p9.ggplot(hill_results, p9.aes(x="k", y="xi_hat"))
+ p9.geom_ribbon(
p9.aes(ymin="xi_hat - 1.96 * se", ymax="xi_hat + 1.96 * se"),
fill="#2E5090", alpha=0.2
)
+ p9.geom_line(color="#2E5090", size=0.8)
+ p9.geom_hline(yintercept=0.33, linetype="dashed", color="red")
+ p9.annotate(
"text", x=hill_results["k"].max() * 0.7, y=0.35,
label="ξ = 1/3 (cubic tail)", color="red", size=8
)
+ p9.labs(
x="Number of Order Statistics (k)",
y="Hill Estimate ξ̂",
title="Hill Plot for Market Return Left Tail"
)
+ p9.theme_minimal()
+ p9.theme(figure_size=(10, 5))
)The Hill plot is a standard diagnostic in EVT. A stable region of the Hill estimate across a range of \(k\) suggests a reliable tail index estimate. The red dashed line at \(\xi = 1/3\) corresponds to a tail index of \(\alpha = 3\), which implies that the third moment (skewness) is infinite, which is a finding consistent with much of the empirical finance literature (Cont 2001; Gabaix et al. 2003).
41.3.3 Block Maxima and the GEV Distribution
The Fisher-Tippett-Gnedenko theorem establishes that if the maximum of \(n\) i.i.d. random variables, after proper normalization, converges in distribution, the limit must be a Generalized Extreme Value (GEV) distribution:
\[ G_\xi(x) = \exp\left\{-\left(1 + \xi \cdot \frac{x - \mu}{\sigma}\right)^{-1/\xi}\right\} \tag{41.22}\]
for \(1 + \xi(x - \mu)/\sigma > 0\), where \(\mu\) is the location, \(\sigma > 0\) the scale, and \(\xi\) the shape parameter. The three sub-families are in Table 41.3.
| Shape \(\xi\) | Distribution | Tail Type | Finance Example |
|---|---|---|---|
| \(\xi > 0\) | Fréchet | Heavy (polynomial decay) | Equity returns, credit losses |
| \(\xi = 0\) | Gumbel | Light (exponential decay) | Normal/lognormal models |
| \(\xi < 0\) | Weibull | Bounded upper tail | Bounded loss distributions |
For the block maxima approach, we divide the return series into non-overlapping blocks (typically months or quarters) and extract the minimum return (maximum loss) from each block. The GEV is then fitted to these block minima.
from scipy.stats import genextreme
# Monthly block minima (maximum losses)
market_daily["month"] = market_daily["date"].dt.to_period("M")
block_minima = (
market_daily.groupby("month")
.agg(min_ret=("mkt_ret", "min"))
.reset_index()
)
# Fit GEV to negated block minima (so we model maxima of losses)
block_losses = -block_minima["min_ret"].values
xi_hat, loc_hat, scale_hat = genextreme.fit(block_losses)
print(f"GEV fit to monthly block maxima of losses:")
print(f" Shape (ξ): {xi_hat:.4f}")
print(f" Location (μ): {loc_hat:.4f}")
print(f" Scale (σ): {scale_hat:.4f}")x_grid = np.linspace(
block_losses.min() * 0.8,
block_losses.max() * 1.2,
200
)
gev_pdf = genextreme.pdf(x_grid, xi_hat, loc=loc_hat, scale=scale_hat)
plot_data = pd.DataFrame({
"x": np.concatenate([block_losses, x_grid]),
"source": (["Empirical"] * len(block_losses) +
["GEV Fit"] * len(x_grid)),
"density": np.concatenate([
np.full(len(block_losses), np.nan),
gev_pdf
]),
"value": np.concatenate([block_losses, np.full(len(x_grid), np.nan)])
})
empirical = plot_data[plot_data["source"] == "Empirical"]
fitted = plot_data[plot_data["source"] == "GEV Fit"]
(
p9.ggplot()
+ p9.geom_histogram(
data=empirical, mapping=p9.aes(x="value", y="..density.."),
bins=30, fill="#2E5090", alpha=0.5
)
+ p9.geom_line(
data=fitted, mapping=p9.aes(x="x", y="density"),
color="#C0392B", size=1
)
+ p9.labs(
x="Monthly Maximum Loss",
y="Density",
title="GEV Distribution Fit to Monthly Block Maxima"
)
+ p9.theme_minimal()
+ p9.theme(figure_size=(10, 5))
)41.3.4 Peaks Over Threshold and the GPD
The Peaks Over Threshold (POT) approach is generally preferred over block maxima because it uses tail data more efficiently. The Pickands-Balkema-de Haan theorem establishes that for a sufficiently high threshold \(u\), the distribution of exceedances \(Y = X - u \mid X > u\) converges to a Generalized Pareto Distribution (GPD):
\[ H_{\xi,\beta}(y) = \begin{cases} 1 - \left(1 + \frac{\xi y}{\beta}\right)^{-1/\xi} & \text{if } \xi \neq 0 \\ 1 - \exp\left(-\frac{y}{\beta}\right) & \text{if } \xi = 0 \end{cases} \tag{41.23}\]
for \(y > 0\) and \(1 + \xi y / \beta > 0\), where \(\beta > 0\) is the scale parameter and \(\xi\) is the shape parameter (identical to the GEV shape). The POT-based VaR and ES at level \(\alpha\) are:
\[ \text{VaR}_\alpha^{\text{GPD}} = u + \frac{\hat{\beta}}{\hat{\xi}}\left[\left(\frac{n}{N_u} \cdot \alpha\right)^{-\hat{\xi}} - 1\right] \tag{41.24}\]
\[ \text{ES}_\alpha^{\text{GPD}} = \frac{\text{VaR}_\alpha^{\text{GPD}}}{1 - \hat{\xi}} + \frac{\hat{\beta} - \hat{\xi} u}{1 - \hat{\xi}} \tag{41.25}\]
where \(N_u\) is the number of observations exceeding \(u\) and \(n\) is the total sample size. These formulae are valid for \(\hat{\xi} < 1\).
from scipy.stats import genpareto
def fit_gpd_pot(returns, threshold_quantile=0.95):
"""
Fit GPD to exceedances over threshold using POT approach.
Parameters
----------
returns : array-like
Return series (losses should be positive).
threshold_quantile : float
Quantile for threshold selection.
Returns
-------
dict : Threshold, GPD parameters, VaR, ES estimates.
"""
losses = -np.array(returns)
losses = losses[~np.isnan(losses)]
n = len(losses)
u = np.quantile(losses, threshold_quantile)
exceedances = losses[losses > u] - u
N_u = len(exceedances)
# Fit GPD
xi_hat, _, beta_hat = genpareto.fit(exceedances, floc=0)
# VaR and ES at 1%
alpha = 0.01
var_gpd = u + (beta_hat / xi_hat) * (
(n / N_u * alpha)**(-xi_hat) - 1
)
es_gpd = var_gpd / (1 - xi_hat) + (beta_hat - xi_hat * u) / (1 - xi_hat)
return {
"threshold": u,
"n_exceedances": N_u,
"xi_hat": xi_hat,
"beta_hat": beta_hat,
"var_1pct": var_gpd,
"es_1pct": es_gpd
}
# Fit GPD to market losses
gpd_results = fit_gpd_pot(market_daily["mkt_ret"].dropna().values, 0.95)
print("GPD-POT Results (1% level):")
for k, v in gpd_results.items():
if isinstance(v, float):
print(f" {k}: {v:.6f}")
else:
print(f" {k}: {v}")41.3.5 Threshold Selection
The choice of threshold \(u\) is the key practical decision in POT analysis. Too low a threshold violates the asymptotic approximation; too high wastes data. The mean residual life plot provides a graphical diagnostic: if the GPD approximation holds above \(u_0\), the mean excess function \(e(u) = E[X - u \mid X > u]\) is linear in \(u\) for \(u > u_0\):
\[ e(u) = \frac{\beta + \xi u}{1 - \xi} \tag{41.26}\]
losses = -market_daily["mkt_ret"].dropna().values
losses_sorted = np.sort(losses)
thresholds = np.quantile(losses, np.linspace(0.80, 0.99, 50))
mean_excess = []
for u in thresholds:
exceedances = losses[losses > u] - u
if len(exceedances) > 5:
mean_excess.append({
"threshold": u,
"mean_excess": exceedances.mean(),
"se": exceedances.std() / np.sqrt(len(exceedances)),
"n_exceed": len(exceedances)
})
mrl_df = pd.DataFrame(mean_excess)
(
p9.ggplot(mrl_df, p9.aes(x="threshold", y="mean_excess"))
+ p9.geom_ribbon(
p9.aes(
ymin="mean_excess - 1.96 * se",
ymax="mean_excess + 1.96 * se"
),
fill="#2E5090", alpha=0.2
)
+ p9.geom_line(color="#2E5090", size=0.8)
+ p9.geom_point(color="#2E5090", size=1.5)
+ p9.labs(
x="Threshold (u)",
y="Mean Excess e(u)",
title="Mean Residual Life Plot"
)
+ p9.theme_minimal()
+ p9.theme(figure_size=(10, 5))
)A threshold in the region where the mean residual life plot is approximately linear is appropriate. In practice, we also compare GPD parameter stability across a range of thresholds and select the lowest threshold at which estimates stabilize.
41.3.6 What EVT Captures That GARCH Does Not
Table 41.4 summarizes the complementary roles of these two modeling frameworks.
| Feature | GARCH | EVT |
|---|---|---|
| What is modeled | Conditional variance dynamics | Unconditional tail distribution |
| Distributional assumption | Full distribution (Gaussian or \(t\)) | Only the tail (GPD or GEV) |
| Time dependence | Explicit (volatility clustering) | None (i.i.d. or filtered) |
| Tail shape | Determined by innovation distribution | Estimated from data |
| Extreme quantiles | Extrapolation from parametric assumption | Principled extrapolation via EVT theorems |
| Best use | Short-horizon risk forecasting | Extreme quantile estimation, stress testing |
The GARCH-filtered EVT approach of McNeil and Frey (2000) combines both: GARCH captures the time-varying volatility, and EVT models the residual tail. This hybrid approach is considered best practice for financial risk management (McNeil, Frey, and Embrechts 2015).
41.4 Tail Dependence and Contagion
41.4.1 Why Correlation Breaks Down in Crises
Perhaps the most important practical fact about multivariate tail risk is that correlations estimated from the body of the distribution systematically understate dependence in the tails. Longin and Solnik (2001) demonstrate this using extreme value theory: for bivariate normal returns, the correlation of extreme realizations converges to zero as the threshold moves further into the tail. Yet empirically, the correlation of crash returns remains large and often increases (a phenomenon incompatible with multivariate normality).
Formally, the coefficient of lower tail dependence between two return series \(X\) and \(Y\) is:
\[ \lambda_L = \lim_{q \to 0^+} P\left(Y \leq F_Y^{-1}(q) \mid X \leq F_X^{-1}(q)\right) \tag{41.27}\]
If \(\lambda_L > 0\), the two series exhibit asymptotic tail dependence (i.e., extreme losses tend to occur jointly). For the bivariate normal distribution, \(\lambda_L = 0\) for all \(|\rho| < 1\) (asymptotic tail independence). This means that any model built on multivariate normality will structurally underestimate the probability of joint crashes. The Student-\(t\) copula, by contrast, has \(\lambda_L > 0\) for \(\rho > -1\), which is why it has become the workhorse model for joint tail risk in finance (Demarta and McNeil 2005).
41.4.2 Co-Crash Measures
We implement several empirical measures of joint crash behavior.
Exceedance Correlation. Following Longin and Solnik (2001), compute the correlation of returns conditional on both returns falling below a threshold \(\theta\):
\[ \rho^-(\theta) = \text{Corr}(r_i, r_j \mid r_i < \theta, r_j < \theta) \tag{41.28}\]
If returns are bivariate normal, \(\rho^-(\theta) \to 0\) as \(\theta \to -\infty\). Empirically, \(\rho^-(\theta)\) typically increases as \(\theta\) decreases, indicating tail dependence beyond what the normal model implies.
Co-Exceedance Count. For each day \(t\), count the number of stocks or sectors whose returns fall below the \(q\)-th percentile of their respective marginal distributions:
\[ C_t(q) = \sum_{i=1}^{N} \mathbb{1}\left(r_{i,t} < F_i^{-1}(q)\right) \tag{41.29}\]
Days with high \(C_t(q)\) represent system-wide stress events.
def empirical_tail_dependence(x, y, quantiles=None):
"""
Estimate lower tail dependence coefficient at various quantiles.
Uses the empirical copula approach: transform to uniform margins,
then estimate P(V <= q | U <= q) for decreasing q.
"""
if quantiles is None:
quantiles = [0.20, 0.15, 0.10, 0.05, 0.03, 0.01]
# Transform to uniform margins via empirical CDF
from scipy.stats import rankdata
n = len(x)
u = rankdata(x) / (n + 1)
v = rankdata(y) / (n + 1)
results = []
for q in quantiles:
mask_u = u <= q
mask_both = mask_u & (v <= q)
n_u = mask_u.sum()
n_both = mask_both.sum()
lambda_hat = n_both / n_u if n_u > 0 else np.nan
results.append({
"quantile": q,
"lambda_L": lambda_hat,
"n_below_x": int(n_u),
"n_co_crash": int(n_both)
})
return pd.DataFrame(results)# Load sector-level daily returns
sector_returns = dc.get_sector_returns(
start_date="2010-01-01",
end_date="2024-12-31",
frequency="daily"
)
# Compute pairwise tail dependence for major sectors
sectors = sector_returns["sector"].unique()[:8] # Top 8 sectors
import itertools
pairwise_td = []
for s1, s2 in itertools.combinations(sectors, 2):
r1 = sector_returns.loc[
sector_returns["sector"] == s1, ["date", "ret"]
].set_index("date")["ret"]
r2 = sector_returns.loc[
sector_returns["sector"] == s2, ["date", "ret"]
].set_index("date")["ret"]
# Align dates
aligned = pd.concat([r1, r2], axis=1, join="inner").dropna()
if len(aligned) < 500:
continue
td = empirical_tail_dependence(
aligned.iloc[:, 0].values,
aligned.iloc[:, 1].values,
quantiles=[0.05]
)
pairwise_td.append({
"sector_1": s1,
"sector_2": s2,
"lambda_L_5pct": td["lambda_L"].iloc[0],
"pearson_corr": aligned.iloc[:, 0].corr(aligned.iloc[:, 1])
})
td_df = pd.DataFrame(pairwise_td)td_display = td_df.copy()
td_display["lambda_L_5pct"] = td_display["lambda_L_5pct"].round(4)
td_display["pearson_corr"] = td_display["pearson_corr"].round(4)
td_display["ratio"] = (
td_display["lambda_L_5pct"] / td_display["pearson_corr"]
).round(3)
td_display.sort_values("lambda_L_5pct", ascending=False).head(15)The ratio column in Table 41.5 is particularly informative. When \(\lambda_L / \rho \gg 1\), the sectors exhibit disproportionately high co-crash frequency relative to their overall correlation. This is the hallmark of tail dependence that Gaussian models miss.
41.4.3 System-Wide Stress Transmission
To assess system-wide contagion dynamics, we construct daily co-exceedance counts and examine their time-series properties.
# Compute daily co-exceedance count at the 5% level
sector_wide = sector_returns.pivot_table(
index="date", columns="sector", values="ret"
)
# Rolling 252-day quantile thresholds
quantile_thresholds = sector_wide.rolling(252, min_periods=60).quantile(0.05)
# Indicator: return below rolling 5th percentile
exceedance_indicators = (sector_wide < quantile_thresholds).astype(int)
co_exceedance = exceedance_indicators.sum(axis=1)
co_exceedance.name = "co_exceedance"
co_exceed_df = co_exceedance.reset_index()
co_exceed_df.columns = ["date", "co_exceedance"]
co_exceed_df["n_sectors"] = exceedance_indicators.notna().sum(axis=1).values(
p9.ggplot(co_exceed_df.dropna(), p9.aes(x="date", y="co_exceedance"))
+ p9.geom_line(color="#2E5090", alpha=0.5, size=0.3)
+ p9.geom_smooth(method="lowess", color="#C0392B", size=1, se=False)
+ p9.labs(
x="",
y="Number of Sectors in Left Tail",
title="System-Wide Stress: Daily Co-Exceedance Count"
)
+ p9.theme_minimal()
+ p9.theme(figure_size=(12, 5))
)Peaks in the co-exceedance series correspond to system-wide stress episodes. The smoothed trend (red line) reveals the evolving baseline level of systemic fragility.
41.4.4 Financial Contagion Mechanisms
The literature distinguishes several channels through which stress propagates across assets, sectors, and markets (Forbes and Rigobon 2002; Kaminsky, Reinhart, and Vegh 2002):
Direct linkages. Balance sheet interconnections, interbank lending, and cross-holdings create direct transmission channels. When institution \(A\) holds assets issued by institution \(B\), losses at \(B\) directly impair \(A\)’s balance sheet.
Information contagion. Adverse news about one firm triggers reassessment of similarly exposed firms. King and Wadhwani (1990) and Kodres and Pritsker (2002) formalize this via Bayesian updating: investors cannot distinguish idiosyncratic from systematic shocks, so a crash in one asset leads them to update beliefs about common risk factors.
Liquidity contagion. Fire sales by distressed institutions depress prices of commonly held assets, creating losses for other holders. Brunnermeier (2009) model this as a liquidity spiral: declining prices trigger margin calls, which force further sales, further depressing prices. In Vietnamese markets, the absence of circuit breakers beyond price limits and the concentration of holdings among a few large institutional investors amplify this channel.
Behavioral contagion. Herding behavior (i.e., investors mimicking others’ trades) generates correlated selling pressure even in the absence of fundamental linkages. Bikhchandani, Hirshleifer, and Welch (1992) and Banerjee (1992) provide theoretical foundations. Choe, Kho, and Stulz (2005) document herding in emerging markets specifically.
We can test for contagion (as distinct from interdependence) using the Forbes and Rigobon (2002) framework. The null hypothesis is that the co-movement observed during a crisis period is fully explained by the increase in volatility (which mechanically inflates correlations). The adjusted correlation is:
\[ \rho^* = \frac{\rho_{\text{crisis}}}{\sqrt{1 + \delta[1 - \rho_{\text{crisis}}^2]}} \tag{41.30}\]
where \(\delta = \sigma_{\text{crisis}}^2 / \sigma_{\text{tranquil}}^2 - 1\) is the relative increase in variance. If \(\rho^* > \rho_{\text{tranquil}}\), the increase in co-movement exceeds what volatility adjustment alone predicts, which is evidence of contagion.
def forbes_rigobon_test(returns_df, crisis_start, crisis_end):
"""
Forbes-Rigobon adjusted correlation test for contagion.
Parameters
----------
returns_df : DataFrame
Columns are sector returns, index is date.
crisis_start, crisis_end : str
Crisis period boundaries.
Returns
-------
DataFrame : Pairwise contagion test results.
"""
tranquil = returns_df[returns_df.index < crisis_start].dropna()
crisis = returns_df[
(returns_df.index >= crisis_start) &
(returns_df.index <= crisis_end)
].dropna()
sectors = returns_df.columns.tolist()
results = []
for s1, s2 in itertools.combinations(sectors, 2):
if s1 not in tranquil.columns or s2 not in crisis.columns:
continue
rho_t = tranquil[s1].corr(tranquil[s2])
rho_c = crisis[s1].corr(crisis[s2])
var_t = tranquil[s1].var()
var_c = crisis[s1].var()
if var_t == 0:
continue
delta = var_c / var_t - 1
rho_adj = rho_c / np.sqrt(1 + delta * (1 - rho_c**2))
results.append({
"sector_1": s1,
"sector_2": s2,
"rho_tranquil": round(rho_t, 4),
"rho_crisis": round(rho_c, 4),
"rho_adjusted": round(rho_adj, 4),
"contagion": rho_adj > rho_t
})
return pd.DataFrame(results)# Test for contagion during COVID-19 crash
contagion_covid = forbes_rigobon_test(
sector_wide,
crisis_start="2020-02-01",
crisis_end="2020-04-30"
)
n_pairs = len(contagion_covid)
n_contagion = contagion_covid["contagion"].sum()
print(f"COVID-19 contagion test: {n_contagion}/{n_pairs} sector pairs "
f"show evidence of contagion ({100*n_contagion/n_pairs:.1f}%)")contagion_covid.sort_values(
"rho_adjusted", ascending=False
).head(15)41.5 Applications to Financial Stability
41.5.1 Market-Wide Stress Indicators
We construct several market-wide stress indicators that aggregate the individual-level and sectoral tail risk measures developed above into a comprehensive system-level diagnostic.
Aggregate Crash Risk Index. The cross-sectional average of firm-level NCSKEW provides a market-wide measure of crash risk that is forward-looking (it captures the buildup of negative skewness before a crash materializes):
\[ \text{ACRI}_y = \frac{1}{N_y}\sum_{i=1}^{N_y} \text{NCSKEW}_{i,y} \tag{41.31}\]
Market Tail Risk (MTR). Following Kelly and Jiang (2014), we estimate the time-varying tail risk of the market portfolio using option-implied or realized tail measures. In the absence of a liquid options market in Vietnam, we use the realized version based on the Hill tail index estimated from rolling windows:
\[ \text{MTR}_t = \hat{\xi}_t^{\text{Hill}} \cdot \hat{\sigma}_t \tag{41.32}\]
where \(\hat{\xi}_t^{\text{Hill}}\) is the Hill estimator computed from a trailing 252-day window and \(\hat{\sigma}_t\) is the conditional volatility (e.g., from GARCH). This interaction captures both the thickness of the tail and the current scale of volatility.
# Aggregate Crash Risk Index
acri = (
crash_risk.groupby("year")
.agg(
acri=("ncskew", "mean"),
median_ncskew=("ncskew", "median"),
pct_high_crash=("ncskew", lambda x: (x > x.quantile(0.75)).mean()),
n_firms=("ncskew", "count")
)
.reset_index()
)
# Market Tail Risk: rolling Hill estimator x rolling volatility
mkt_returns = market_daily["mkt_ret"].dropna().values
dates = market_daily["date"].dropna().values
mtr_list = []
for end_idx in range(504, len(mkt_returns)):
window_ret = mkt_returns[end_idx - 252:end_idx]
losses_pos = -window_ret[window_ret < 0]
if len(losses_pos) < 30:
continue
# Hill estimator at k = 50
sorted_losses = np.sort(losses_pos)[::-1]
k = min(50, len(sorted_losses) - 1)
if k < 10:
continue
log_excess = np.log(sorted_losses[:k]) - np.log(sorted_losses[k])
xi_hat = log_excess.mean()
sigma_hat = window_ret.std()
mtr = xi_hat * sigma_hat
mtr_list.append({
"date": dates[end_idx],
"xi_hat": xi_hat,
"sigma_hat": sigma_hat,
"mtr": mtr
})
mtr_df = pd.DataFrame(mtr_list)
mtr_df["date"] = pd.to_datetime(mtr_df["date"])(
p9.ggplot(mtr_df, p9.aes(x="date", y="mtr"))
+ p9.geom_line(color="#2E5090", alpha=0.5, size=0.4)
+ p9.geom_smooth(method="lowess", color="#C0392B", size=1, se=False)
+ p9.labs(
x="",
y="MTR (ξ × σ)",
title="Market Tail Risk: Rolling Hill Index × Conditional Volatility"
)
+ p9.theme_minimal()
+ p9.theme(figure_size=(12, 5))
)41.5.2 Sectoral Fragility
Not all sectors are equally vulnerable to tail events. We decompose system-wide tail risk into sectoral contributions to identify the most fragile parts of the market.
# Sector-level tail statistics
sector_tail_stats = []
for sector in sectors:
sector_data = sector_returns[sector_returns["sector"] == sector]["ret"].dropna()
if len(sector_data) < 500:
continue
s_data = sector_data.values
# Skewness and kurtosis
skew = stats.skew(s_data)
kurt = stats.kurtosis(s_data)
# VaR and ES at 1%
var_1 = -np.quantile(s_data, 0.01)
es_1 = -s_data[s_data <= np.quantile(s_data, 0.01)].mean()
# Hill tail index (k=50)
losses_pos = -s_data[s_data < 0]
sorted_l = np.sort(losses_pos)[::-1]
k = min(50, len(sorted_l) - 1)
if k >= 10:
log_exc = np.log(sorted_l[:k]) - np.log(sorted_l[k])
xi = log_exc.mean()
else:
xi = np.nan
sector_tail_stats.append({
"sector": sector,
"skewness": round(skew, 3),
"excess_kurtosis": round(kurt, 3),
"var_1pct": round(var_1, 4),
"es_1pct": round(es_1, 4),
"tail_index_xi": round(xi, 3) if not np.isnan(xi) else np.nan,
"n_obs": len(s_data)
})
fragility_df = pd.DataFrame(sector_tail_stats)fragility_df.sort_values("es_1pct", ascending=False)(
p9.ggplot(
fragility_df.dropna(),
p9.aes(x="tail_index_xi", y="es_1pct", label="sector")
)
+ p9.geom_point(color="#2E5090", size=3)
+ p9.geom_text(nudge_y=0.002, size=7)
+ p9.labs(
x="Tail Index (ξ, Hill)",
y="1% Expected Shortfall",
title="Sector Fragility Map"
)
+ p9.theme_minimal()
+ p9.theme(figure_size=(10, 7))
)Sectors in the upper-right quadrant of Figure 41.10 exhibit both high expected shortfall (large average losses in the tail) and high tail index (heavier-than-typical tails). These represent the most fragile components of the market from a systemic stability perspective.
41.5.3 Crisis Diagnostics
We construct a comprehensive crisis diagnostic framework that combines the indicators developed above to identify, characterize, and compare stress episodes in Vietnamese financial markets.
# Identify extreme co-exceedance days (system-wide stress events)
threshold_high_stress = co_exceed_df["co_exceedance"].quantile(0.99)
stress_events = co_exceed_df[
co_exceed_df["co_exceedance"] >= threshold_high_stress
].copy()
# Add market return on stress days
stress_events = stress_events.merge(
market_daily[["date", "mkt_ret"]], on="date", how="left"
)
print(f"System-wide stress days (top 1%): {len(stress_events)}")
print(f"Average market return on stress days: "
f"{stress_events['mkt_ret'].mean():.4f}")
print(f"Average co-exceedance on stress days: "
f"{stress_events['co_exceedance'].mean():.1f}")# Cluster stress days into episodes (within 10 trading days)
stress_events = stress_events.sort_values("date").reset_index(drop=True)
stress_events["gap"] = stress_events["date"].diff().dt.days
stress_events["episode"] = (stress_events["gap"] > 15).cumsum()
episode_summary = (
stress_events.groupby("episode")
.agg(
start_date=("date", "min"),
end_date=("date", "max"),
n_stress_days=("date", "count"),
avg_market_return=("mkt_ret", "mean"),
min_market_return=("mkt_ret", "min"),
avg_coexceedance=("co_exceedance", "mean")
)
.reset_index(drop=True)
.sort_values("min_market_return")
)
episode_summary["avg_market_return"] = episode_summary["avg_market_return"].round(4)
episode_summary["min_market_return"] = episode_summary["min_market_return"].round(4)
episode_summary["avg_coexceedance"] = episode_summary["avg_coexceedance"].round(1)
episode_summary.head(10)41.6 Summary
This chapter developed a comprehensive toolkit for measuring and diagnosing tail risk in Vietnamese equity markets. The key objects estimated (crash risk measures, VaR and Expected Shortfall, tail indices, and tail dependence coefficients) capture distinct but complementary aspects of extreme event behavior.
At the individual stock level, the NCSKEW, DUVOL, and crash count measures quantify asymmetry in firm-specific return distributions, providing forward-looking indicators of crash vulnerability. At the portfolio level, EVT-based approaches to VaR and ES estimation avoid the parametric misspecification that plagues Gaussian and even Student-\(t\) models in the deep tails. At the system level, tail dependence and co-exceedance measures reveal the extent to which crashes propagate across sectors, and the Forbes-Rigobon contagion test distinguishes genuine contagion from the mechanical increase in co-movement driven by higher volatility.
Several features of Vietnamese markets deserve emphasis. Price limits censor the observed return distribution, causing systematic underestimation of true tail risk. The concentration of trading among retail investors amplifies herding-driven crash dynamics. State ownership of large firms creates correlated exposure to policy shocks that is not captured by standard diversification. And the absence of deep derivatives markets limits the hedging instruments available for tail risk management, making accurate measurement all the more critical.