Case Study: Financial Volatility & Risk¶
Not every series is forecast for its level. A stock's daily return is close to unpredictable in direction from its own past, yet its risk, the size of the moves, is very forecastable, because volatility clusters. This case study runs the full risk workflow on five years of daily returns: confirm the stylized facts, fit a GARCH model, forecast volatility, turn it into a Value-at-Risk number a desk can act on, backtest that number, and check what leads the stock. We use statsmodels and the arch library.
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import seaborn as sns # seaborn = high-level statistical plots (heatmaps, pairplots, count/bar plots)
from matplotlib.colors import ListedColormap
EM="#0284c7"; DEEP="#075985"; LIGHT="#bae6fd"; INK="#1a2138"; GRID="#e6e9f2"; RED="#ef4444"; AMBER="#d97706"; GREEN="#059669"; BLUE="#2563eb"; PUR="#9333ea"; GREY="#94a3b8"; SLATE="#475569"; ORG="#0284c7"; CYAN="#0891b2"
plt.rcParams.update({"figure.facecolor":"white","axes.facecolor":"white","figure.dpi":110,"font.size":11,
"axes.edgecolor":GRID,"axes.grid":True,"grid.color":GRID,"axes.axisbelow":True,"axes.spines.top":False,
"axes.spines.right":False,"axes.titlesize":12,"axes.titleweight":"bold","legend.frameon":False})
sns.set_style("whitegrid")
BASE="https://raw.githubusercontent.com/johnfisher-ai/Statistics-Data-Science-AI-Visual-Book/main/data/"
from scipy import stats
from statsmodels.stats.diagnostic import acorr_ljungbox, het_arch
from statsmodels.tsa.stattools import grangercausalitytests, acf
from statsmodels.graphics.tsaplots import plot_acf
import warnings; warnings.filterwarnings('ignore')
try:
from arch import arch_model
except ImportError:
import subprocess, sys; subprocess.run([sys.executable,'-m','pip','install','-q','arch']); from arch import arch_model
A desk holds a $1,000,000 position in one stock and asks a concrete question: how much could we lose in a single day? The answer is Value-at-Risk (VaR): the 1-day 99% VaR is the loss we expect to exceed on only 1 day in 100. Because losses are driven by volatility, and volatility is forecastable, VaR is a forecasting problem in disguise.
try: raw = pd.read_excel('../../data/case-study-financial-volatility-and-risk--stock_returns.xlsx', sheet_name='Data')
except FileNotFoundError: raw = pd.read_excel(BASE + 'case-study-financial-volatility-and-risk--stock_returns.xlsx', sheet_name='Data')
raw['date']=pd.to_datetime(raw['date']); df = raw.set_index('date')
r, mkt = df['stock_return'], df['market_return']
print('days:', len(r), '| daily mean %.3f%% | daily std %.3f%%' % (r.mean(), r.std()))
days: 1260 | daily mean 0.071% | daily std 1.206%
fig,ax=plt.subplots(figsize=(11,3.8)); ax.plot(r.index, r.values, color=EM, lw=0.7)
ax.axhline(0,color=GREY); ax.set(title='Daily stock returns (%): quiet spells punctuated by bursts', ylabel='return %')
worst=r.idxmin(); ax.plot([worst],[r[worst]],'v',color=RED,ms=9,label='worst day %.1f%%'%r[worst]); ax.legend(); plt.tight_layout(); plt.show()
print('big moves cluster together -> volatility has memory, even though the direction does not')
big moves cluster together -> volatility has memory, even though the direction does not
The returns wander around zero with no visible trend, but the size of the moves comes in waves: quiet weeks, then bursts of large swings (including a crisis stretch). That is volatility clustering, the central fact of financial risk, and the reason a single average volatility number is not enough.
print('ACF of returns lag-1 = %.2f (direction is ~unpredictable)' % acf(r,nlags=1)[1])
print('ACF of returns^2 lag-1 = %.2f (VARIANCE has memory)' % acf((r-r.mean())**2,nlags=1)[1])
print('kurtosis = %.1f (a normal is 3; >3 means FAT TAILS -> extremes are common)' % stats.kurtosis(r,fisher=False))
fig,ax=plt.subplots(1,2,figsize=(11,3.8))
ax[0].hist(r, bins=60, density=True, color=LIGHT, edgecolor='white'); x=np.linspace(r.min(),r.max(),200)
ax[0].plot(x, stats.norm.pdf(x, r.mean(), r.std()), color=RED, lw=2, label='normal')
ax[0].set(title='Returns are peaked with fat tails (vs normal)', xlabel='return %'); ax[0].legend()
plot_acf((r-r.mean())**2, lags=20, ax=ax[1]); ax[1].set(title='ACF of squared returns: volatility clusters')
plt.tight_layout(); plt.show()
ACF of returns lag-1 = -0.04 (direction is ~unpredictable) ACF of returns^2 lag-1 = 0.12 (VARIANCE has memory) kurtosis = 6.3 (a normal is 3; >3 means FAT TAILS -> extremes are common)
Three facts, all visible: returns are close to white noise in the mean (a lag-1 autocorrelation of only -0.04, so their own past says almost nothing about tomorrow's direction), their squared values are autocorrelated (volatility is predictable), and the distribution has fat tails, kurtosis about 6.3 versus 3 for a normal, so crashes happen far more often than a bell curve says. A good risk model must honor all three.
lm, p, _, _ = het_arch(r, nlags=10)
print('ARCH-LM test: statistic %.1f, p = %.4f' % (lm, p))
print('p far below 0.05 -> significant ARCH effects -> volatility clustering is real and forecastable')
ARCH-LM test: statistic 86.5, p = 0.0000 p far below 0.05 -> significant ARCH effects -> volatility clustering is real and forecastable
The ARCH-LM test formally checks whether past squared returns help predict today's variance. The tiny p-value confirms it: there are strong ARCH effects, so a model of the variance (GARCH) will add real value over assuming constant volatility.
g = arch_model(r, mean='Constant', vol='GARCH', p=1, q=1, dist='t').fit(disp='off')
print('omega %.3f | alpha %.3f | beta %.3f | persistence (alpha+beta) %.3f' % (g.params['omega'], g.params['alpha[1]'], g.params['beta[1]'], g.params['alpha[1]']+g.params['beta[1]']))
print('Student-t degrees of freedom nu = %.1f (low -> confirms fat tails)' % g.params['nu'])
omega 0.037 | alpha 0.065 | beta 0.910 | persistence (alpha+beta) 0.976 Student-t degrees of freedom nu = 5.1 (low -> confirms fat tails)
GARCH(1,1) says today's variance depends on yesterday's squared shock and yesterday's variance. We use a Student-t distribution for the innovations, because the fat tails demand it. The fit shows persistence about 0.98 (a shock to volatility fades slowly, so turbulence lingers) and a low nu about 5, confirming the heavy tails.
sr = g.std_resid
lb = acorr_ljungbox((sr-sr.mean())**2, lags=[10])['lb_pvalue'].iloc[0]
print('standardized residual^2 Ljung-Box p = %.2f -> %s' % (lb, 'no clustering left (GARCH worked)' if lb>0.05 else 'clustering remains'))
fig,ax=plt.subplots(figsize=(11,3.8)); ax.plot(r.index, r.abs(), color=LIGHT, lw=0.6, label='|return|')
ax.plot(g.conditional_volatility.index, g.conditional_volatility, color=RED, lw=1.6, label='GARCH conditional volatility')
ax.set(title='GARCH conditional volatility tracks the clustering', ylabel='%'); ax.legend(); plt.tight_layout(); plt.show()
standardized residual^2 Ljung-Box p = 0.62 -> no clustering left (GARCH worked)
A good volatility model leaves standardized residuals with no clustering left. Dividing each return by its GARCH volatility, the squared residuals now pass a Ljung-Box test (p about 0.62), so the model absorbed the ARCH structure. The red line is the estimated conditional volatility, rising in the storms and falling in the calm.
fc = g.forecast(horizon=10); vol10 = np.sqrt(fc.variance.iloc[-1].values)
for h in [1,5,10]: print('volatility forecast, day +%2d: %.2f%%' % (h, vol10[h-1]))
print('the forecast drifts toward the long-run level -> after a shock it decays, in calm it drifts up')
volatility forecast, day + 1: 0.99% volatility forecast, day + 5: 1.02% volatility forecast, day +10: 1.05% the forecast drifts toward the long-run level -> after a shock it decays, in calm it drifts up
Because volatility is persistent but mean-reverting, GARCH forecasts a whole path: from today's level it drifts toward the long-run average. After a turbulent day it stays elevated for a while then subsides; in a calm spell it creeps back up. That path, not a single number, is what feeds a forward-looking risk limit.
mu, nu = g.params['mu'], g.params['nu']
q01 = stats.t.ppf(0.01, nu) * np.sqrt((nu-2)/nu) # 1% quantile of the standardized t
cv_next = np.sqrt(g.forecast(horizon=1).variance.iloc[-1, 0]) # GARCH's volatility forecast for tomorrow
long_run = np.sqrt(g.params['omega'] / (1 - g.params['alpha[1]'] - g.params['beta[1]']))
VaR_pct = -(mu + cv_next * q01)
POSITION = 1_000_000
print("tomorrow's forecast volatility: %.2f%% (long-run level %.2f%%)" % (cv_next, long_run))
print('1-day 99%% VaR = %.2f%% of the position = $%.0f' % (VaR_pct, VaR_pct/100*POSITION))
print('interpretation: on 99 days in 100 the loss stays under this; 1 day in 100 it is worse')
tomorrow's forecast volatility: 0.99% (long-run level 1.24%) 1-day 99% VaR = 2.52% of the position = $25239 interpretation: on 99 days in 100 the loss stays under this; 1 day in 100 it is worse
Value-at-Risk converts the volatility forecast into money. We take GARCH's volatility forecast for tomorrow, scale it by the 1% quantile of the Student-t (which respects the fat tails), and multiply by the position. The 1-day 99% VaR is about 2.5%, or 25,200 dollars, on the million-dollar position, the loss the desk should be prepared to exceed roughly two or three times a year.
var_line = -(mu + g.conditional_volatility * q01) # daily loss threshold
viol = r.values < -var_line.values
print('VaR violations: %d of %d days = %.2f%% (expected 1.00%%)' % (viol.sum(), len(r), 100*viol.mean()))
const_var = -(mu + r.std() * q01) # one threshold from the full-sample volatility
cviol = r.values < -const_var
print('constant-volatility VaR (%.2f%% every day): %d violations = %.2f%%' % (const_var, cviol.sum(), 100*cviol.mean()))
print('violations by year, GARCH: ', pd.Series(viol, index=r.index).groupby(r.index.year).sum().to_dict())
print('violations by year, constant:', pd.Series(cviol, index=r.index).groupby(r.index.year).sum().to_dict())
fig,ax=plt.subplots(figsize=(11,3.8)); ax.plot(r.index, r.values, color=LIGHT, lw=0.6, label='return')
ax.plot(var_line.index, -var_line.values, color=EM, lw=1.2, label='99% VaR threshold')
ax.plot(r.index[viol], r.values[viol], 'o', color=RED, ms=4, label='violations'); ax.set(title='VaR backtest: violations land near the expected 1%', ylabel='return %'); ax.legend(); plt.tight_layout(); plt.show()
VaR violations: 15 of 1260 days = 1.19% (expected 1.00%)
constant-volatility VaR (3.08% every day): 18 violations = 1.43%
violations by year, GARCH: {2019: 2, 2020: 4, 2021: 3, 2022: 6, 2023: 0}
violations by year, constant: {2019: 3, 2020: 4, 2021: 5, 2022: 6, 2023: 0}
A VaR number is only credible if it is backtested. Counting the days the actual loss breached the VaR line gives about 1.2%, close to the expected 1% (a formal Kupiec test is in Take It Further). The violations are spread across the first four years, with none in the calm final year. A constant VaR set from the full-sample volatility does worse on the count (18 breaches, 1.43%) with a similar spread by year, so GARCH's edge here is a rate closer to 1% and a threshold that moves with the risk.
gm = grangercausalitytests(df[['stock_return','market_return']], maxlag=2)
gs = grangercausalitytests(df[['market_return','stock_return']], maxlag=2)
print('market -> stock Granger p(lag1) = %.4f -> significant' % gm[1][0]['ssr_ftest'][1])
print('stock -> market Granger p(lag1) = %.4f -> not significant' % gs[1][0]['ssr_ftest'][1])
import statsmodels.api as sm
lagged = pd.concat([r, r.shift(1).rename('stock_lag1'), mkt.shift(1).rename('market_lag1')], axis=1).dropna()
ols = sm.OLS(lagged['stock_return'], sm.add_constant(lagged[['stock_lag1', 'market_lag1']])).fit()
print("size of the effect: %.2f x yesterday's market return (t = %.1f), R-squared %.3f" % (ols.params['market_lag1'], ols.tvalues['market_lag1'], ols.rsquared))
market -> stock Granger p(lag1) = 0.0000 -> significant stock -> market Granger p(lag1) = 0.0606 -> not significant
size of the effect: 0.22 x yesterday's market return (t = 6.8), R-squared 0.037
Risk is not only about one stock in isolation. A Granger causality test shows the market Granger-causes the stock (yesterday's market return helps predict today's stock return) but not the reverse at the 5% level. The regression behind the test puts the effect at about 0.22 of yesterday's market move, clearly significant yet small, since the two lags explain under 4% of the daily variation. So the stock's direction is slightly predictable from the market even though its own past says almost nothing, and a sharp market fall is a mild warning for the next day. The test is about returns, so it says nothing directly about volatility. (As always, Granger is predictive precedence, not proof of cause.)
def daily_var(returns, position=1_000_000, level=0.99):
m = arch_model(returns, mean='Constant', vol='GARCH', p=1, q=1, dist='t').fit(disp='off')
nu = m.params['nu']; q = stats.t.ppf(1-level, nu)*np.sqrt((nu-2)/nu)
sigma = np.sqrt(m.forecast(horizon=1).variance.iloc[-1, 0]) # tomorrow's volatility forecast
v = -(m.params['mu'] + sigma*q)
return v, v/100*position
# score today's risk on the latest data (in production, rerun each afternoon after the close)
pct, dollars = daily_var(r)
print('DEPLOYED RISK REPORT (today):')
print(' 1-day 99%% VaR = %.2f%% = $%.0f on a $1,000,000 position' % (pct, dollars))
print(' drivers: volatility is clustering (persistence 0.98), and yesterday\'s market move carries a small signal')
DEPLOYED RISK REPORT (today): 1-day 99% VaR = 2.52% = $25239 on a $1,000,000 position drivers: volatility is clustering (persistence 0.98), and yesterday's market move carries a small signal
In production the risk function reruns each afternoon: refit GARCH on the latest returns, forecast tomorrow's volatility, and publish the VaR. The plain-English report writes itself, “we could lose more than about 25,200 dollars tomorrow with 1% probability. Risk is a little below its long-run level, and a big move would raise it fast.” That single, backtested, forward-looking number, with the market's lead as a mild warning, is what turns a wall of returns into a decision.
The risk workflow, in one view¶
- Define the risk question in dollars (Value-at-Risk).
- Collect the daily returns.
- Inspect: volatility clusters into calm and turbulent stretches.
- Stylized facts: mean unpredictable, variance has memory, fat tails.
- Test for ARCH effects to confirm volatility is forecastable.
- Fit GARCH with fat-tailed (Student-t) innovations.
- Diagnose: standardized residuals should have no clustering left.
- Forecast volatility; it mean-reverts toward the long run.
- Value-at-Risk: scale volatility by the fat-tailed quantile times the position.
- Backtest the VaR: violations should land near the target rate.
- Relationships: the market Granger-leads the stock (systemic risk).
- Deploy & communicate a daily risk number in plain English.
The one idea to keep: you cannot forecast a return, but you can forecast its risk, and GARCH turns clustered, fat-tailed volatility into a backtested Value-at-Risk the business can actually plan and set limits against.