import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from scipy import stats
import statsmodels.api as sm
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.stats.diagnostic import acorr_ljungboxWeek 3: Distributions and Statistical Analysis
Version 5.0 - August 2026
3 Week 3
Risk measurement rests on an assumption about the distribution of returns, and the usual first assumption — normality — is wrong in a way that matters. Returns have fatter tails than the normal distribution allows, so a model built on it understates the probability of exactly the large losses that risk management exists to anticipate.
This week uses the normal and Student-t distributions, their densities, distribution functions and quantiles, and QQ plots for comparing a sample with a distribution. It then applies the Jarque-Bera test for normality and the Ljung-Box test for autocorrelation, in returns and in squared returns.
Open the main course folder and check that the Prices.csv and Returns.csv files created and validated in Week 2 are present. This week adds scipy and statsmodels, both included with Anaconda.
You should be able to explain why normality is a poor description of financial tails, read a histogram and QQ plot, and interpret Jarque–Bera and Ljung–Box test results for one return series.
3.1 The plan for this week
- Compare normal and Student-t distributions
- Summarise one return series
- Read a histogram and QQ plot
- Interpret normality and autocorrelation tests
3.2 Loading data and packages
Prices = pd.read_csv("Prices.csv", index_col="date", parse_dates=True)
Returns = pd.read_csv("Returns.csv", index_col="date", parse_dates=True)Both frames should have the same increasing date index and the ten security columns created in Week 2. If either file is missing, return to the Week 2 course data checkpoint before continuing.
3.3 Distributions
We work extensively with statistical distributions, such as the normal, log-normal, Student-t, binomial, Bernoulli and chi-square. scipy.stats handles many more distributions than these, but these are the ones we mostly use here. Each distribution object comes with four functions:
- Density (pdf) —
.pdf(), or probability mass (pmf) —.pmf() - Cumulative distribution (cdf) —
.cdf() - Quantiles —
.ppf() - Random numbers —
.rvs()
The first differs between the two kinds of distribution. A continuous distribution, such as the normal or Student-t, has a density, .pdf(), and the probability of any single value is zero. A discrete distribution, such as the Bernoulli or binomial, instead has a probability mass function, .pmf(), giving the probability of an exact outcome. Week 10 uses stats.binom.logpmf() for that reason. The other three names are the same for every distribution:
| Method | Normal | Student-t | Meaning |
|---|---|---|---|
.pdf() |
stats.norm.pdf |
stats.t.pdf |
density |
.cdf() |
stats.norm.cdf |
stats.t.cdf |
cumulative distribution |
.ppf() |
stats.norm.ppf |
stats.t.ppf |
quantile function |
.rvs() |
stats.norm.rvs |
stats.t.rvs |
random numbers |
We deal with random numbers in a later seminar. Here is an example showing how to plot the density, cumulative distribution and quantile function over their domain. The quantile function is drawn from 0.001 to 0.999 rather than from 0 to 1, since a normal quantile is infinite at both ends.
x = np.linspace(-3, 3, 1000)
z = np.linspace(0.001, 0.999, 1000)
fig, axes = plt.subplots(2, 2)
axes[0, 0].plot(x, stats.norm.pdf(x))
axes[0, 0].set_title("Normal density")
axes[0, 1].plot(x, stats.norm.cdf(x))
axes[0, 1].set_title("Cumulative distribution")
axes[1, 0].plot(z, stats.norm.ppf(z))
axes[1, 0].set_title("Normal quantile")
axes[1, 1].set_visible(False)
fig.tight_layout()
plt.show()3.4 Comparing the normal distribution with the Student-t
The Student-t distribution has fatter tails than the normal, and the fewer degrees of freedom it has, the fatter those tails are.
A Student-t with \(\nu\) degrees of freedom has variance \(\nu/(\nu-2)\), so a plain stats.t is wider than a standard normal for two separate reasons, its larger variance and its heavier tails. To see the tails on their own, scale it to variance one, which is also what arch does when it fits a Student-t in Week 6.
nu = 3
t_scale = np.sqrt((nu - 2) / nu)
x = np.linspace(-3, 3, 1000)
normal = stats.norm.pdf(x)
student = stats.t.pdf(x, df=nu, scale=t_scale)
fig, ax = plt.subplots()
ax.plot(x, normal, label="Normal")
ax.plot(x, student, label=f"Student-t, {nu} df, variance 1")
ax.set_title("Comparing distributions")
ax.set_xlabel("x")
ax.set_ylabel("f(x)")
ax.legend(frameon=False)
plt.show()These distributions are important in financial analysis because:
- The normal distribution is the traditional assumption in many financial models and risk calculations.
- The Student-t distribution better captures the fat tails often seen in financial returns. Extreme losses and gains happen more frequently than the normal distribution predicts.
- The choice between these distributions affects risk calculations, portfolio optimisation and derivative pricing.
3.5 Applying distribution concepts to real financial data
Now that we understand the theoretical distributions commonly used in finance, we can apply these concepts to analyse real stock price data. We examine how actual returns compare to theoretical distributions and identify periods where normal distribution assumptions break down.
3.6 Visualising and commenting on prices
Prices.head()| GSPC | IXIC | AAPL | MSFT | JPM | C | XOM | MCD | GE | NVDA | |
|---|---|---|---|---|---|---|---|---|---|---|
| date | ||||||||||
| 2000-01-03 | 1455.22 | 4131.1499 | 0.8401 | 35.6344 | 23.1248 | 209.0070 | 17.9612 | 21.1582 | 129.9334 | 0.0894 |
| 2000-01-04 | 1399.42 | 3901.6899 | 0.7693 | 34.4306 | 22.6175 | 196.1906 | 17.6172 | 20.7244 | 124.7360 | 0.0870 |
| 2000-01-05 | 1402.11 | 3877.5400 | 0.7805 | 34.7937 | 22.4778 | 204.0776 | 18.5776 | 21.0581 | 124.5194 | 0.0842 |
| 2000-01-06 | 1403.45 | 3727.1299 | 0.7130 | 33.6282 | 22.7970 | 213.9364 | 19.5380 | 20.7577 | 126.1842 | 0.0787 |
| 2000-01-07 | 1441.47 | 3882.6201 | 0.7468 | 34.0676 | 23.2158 | 212.9506 | 19.4807 | 21.2917 | 131.0703 | 0.0800 |
3.6.1 GE case
Jack Welch, who was the CEO of GE for twenty years, retired on 7 September 2001. He was considered to be one of the most valuable CEOs of all time. Add a vertical line on our plot to reflect this, and find the highest price GE reached before he left. The course data begin in 2000, so this covers only the last two years of his tenure.
welch_retirement = "2001-09-07"
welch_prices = Prices["GE"][Prices.index <= welch_retirement]
max_price = welch_prices.max()
max_date = welch_prices.idxmax()
fig, ax = plt.subplots()
ax.plot(Prices.index, Prices["GE"])
ax.set_title("Price of GE")
ax.axvline(pd.Timestamp(welch_retirement), linewidth=2, color="red")
ax.axvline(max_date, linewidth=2, color="blue")
ax.annotate(
f"Highest price before\nhis retirement: {max_price:.2f}",
xy=(max_date, max_price),
xytext=(10, 0),
textcoords="offset points",
)
plt.show()3.7 Zoom into the COVID crisis
To compare how different stocks performed during the crisis, we normalise all prices to start at 1. This shows relative performance regardless of absolute price levels — which stocks declined most and which recovered fastest.
crisis_prices = Prices[(Prices.index >= "2020-01-01") & (Prices.index <= "2020-12-31")]
crisis_prices_normalised = crisis_prices / crisis_prices.iloc[0]
fig, ax = plt.subplots()
for ticker in crisis_prices_normalised.columns:
ax.plot(crisis_prices_normalised.index, crisis_prices_normalised[ticker], label=ticker)
ax.legend(loc="upper left", ncol=2, frameon=False)
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=3))
ax.xaxis.set_major_formatter(mdates.DateFormatter("%b %Y"))
ax.set_xlim(crisis_prices.index[0], crisis_prices.index[-1])
ax.tick_params(axis="x", labelrotation=45)
fig.tight_layout()
plt.show()The x-axis shows one label every three months. The labels use abbreviated month names and are rotated so that they do not overlap.
Now apply the same crisis period analysis to returns data to understand volatility patterns during the crisis:
crisis_returns = Returns[(Returns.index >= "2020-01-01") & (Returns.index <= "2020-12-31")]
fig, ax = plt.subplots()
for ticker in crisis_returns.columns:
ax.plot(crisis_returns.index, crisis_returns[ticker], label=ticker)
ax.set_title("Returns during the COVID crisis")
ax.set_xlabel("Date")
ax.set_ylabel("Returns")
ax.legend(loc="lower right", ncol=2, frameon=False)
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=3))
ax.xaxis.set_major_formatter(mdates.DateFormatter("%b %Y"))
ax.set_xlim(crisis_returns.index[0], crisis_returns.index[-1])
ax.tick_params(axis="x", labelrotation=45)
fig.tight_layout()
plt.show()This shows the volatility clustering during the crisis. Periods of high volatility tend to be followed by more high volatility.
3.8 Graphical analyses and statistical tests
We can do some of the basic statistical and graphical analysis shown at the start of the course. Pick JPMorgan. The missing values are dropped once, here, because several of the tests below do not ignore them.
y = Returns["JPM"].dropna()First, print some summary statistics. scipy.stats.kurtosis reports excess kurtosis by default. Passing fisher=False gives raw kurtosis instead, the convention used throughout this course, in which the normal distribution has a kurtosis of 3 rather than 0.
print("mean:", y.mean())
print("sd:", y.std())
print("skewness:", stats.skew(y))
print("kurtosis:", stats.kurtosis(y, fisher=False))mean: 0.00039266060916775365
sd: 0.02326061861157198
skewness: 0.21236522160404447
kurtosis: 17.57827645552578
The histogram shows the centre and tails directly:
fig, ax = plt.subplots()
ax.hist(y, bins=50, density=True, color="lightgrey")
ax.set_title("Distribution of JPM returns")
ax.set_xlabel("Daily log return")
ax.set_ylabel("Density")
plt.show()Then run the Jarque-Bera test for normality and the Ljung-Box test for autocorrelation, on both the returns and the squared returns.
jb_stat, jb_pvalue = stats.jarque_bera(y)
print("Jarque-Bera statistic:", jb_stat, "p-value:", jb_pvalue)Jarque-Bera statistic: 57342.13817289246 p-value: 0.0
acorr_ljungbox(y, lags=[10], return_df=True)| lb_stat | lb_pvalue | |
|---|---|---|
| 10 | 66.190495 | 2.397765e-10 |
acorr_ljungbox(y ** 2, lags=[10], return_df=True)| lb_stat | lb_pvalue | |
|---|---|---|
| 10 | 3714.372294 | 0.0 |
A small Jarque–Bera p-value is evidence against normality. Ljung–Box results for returns and squared returns answer different questions. Dependence in squared returns is evidence of volatility clustering even when dependence in returns is weak.
Left without a lags argument, acorr_ljungbox returns one row per lag up to a default of its own choosing. Passing lags=[10] asks for a single row instead — the cumulative test statistic up to lag 10, a common choice for daily return diagnostics.
Then, plot the autocorrelation function of returns and returns squared. What information does the latter plot provide?
fig, ax = plt.subplots()
plot_acf(y, ax=ax, title="Autocorrelation of returns")
plt.show()fig, ax = plt.subplots()
plot_acf(y ** 2, ax=ax, title="Autocorrelation of returns squared")
plt.show()Finally, the QQ plot is informative about the distribution of returns, especially in the tails. Points that bend away from the line at both ends show that the empirical tails differ from those of the comparison distribution.
sm.qqplot(y, line="q")
plt.gca().set_title("QQ plot against the normal")
plt.show()3.9 Recap
3.9.1 In this seminar, we have covered:
- Working with statistical distributions in
scipy.stats:- Density functions (
.pdf()) - Cumulative distribution functions (
.cdf()) - Quantile functions (
.ppf()) - Comparing normal and Student-t distributions
- Density functions (
- Visualising and analysing stock price data:
- Adding reference lines for important events
- Finding and marking maximum values
- Creating crisis period analysis
- Statistical testing and graphical analysis:
- Calculating summary statistics (mean, standard deviation, skewness, kurtosis)
- Testing for normality and autocorrelation
- Creating ACF plots for returns and squared returns
- Using QQ plots to assess distributional assumptions
3.9.2 Some new functions used:
stats.norm.pdf()/.cdf()/.ppf()/.rvs()— the density/cdf/quantile/random-number quartet for the normal distributionstats.t.pdf()— Student-t density functionnp.linspace()— generate evenly spaced sequences of numbersax.axvline()— add a vertical reference line to a plotax.annotate()— add text annotations to a plotstats.skew()/stats.kurtosis()— measure asymmetry and tail heavinessstats.jarque_bera()— test for normality based on skewness and kurtosisacorr_ljungbox()— test for autocorrelation in a time seriesplot_acf()— plot the autocorrelation functionsm.qqplot()— quantile-quantile plot for comparing distributions
3.10 Exercises
3.10.1 Apply the method
Change y = Returns["JPM"].dropna() to another security. Before running the code, predict whether its returns will look normally distributed. Use its summary statistics, histogram, QQ plot and test results to give a short evidence-based answer.
3.10.2 Compare results
Compare the JPM QQ plot against Student-t distributions with 4 and 3 degrees of freedom:
fig, axes = plt.subplots(1, 2, figsize=(9, 4))
sm.qqplot(y, dist=stats.t, distargs=(4,), line="q", ax=axes[0])
axes[0].set_title("QQ plot against t(4)")
sm.qqplot(y, dist=stats.t, distargs=(3,), line="q", ax=axes[1])
axes[1].set_title("QQ plot against t(3)")
fig.tight_layout()
plt.show()Which comparison follows the tails more closely, and what does that imply about the normal assumption?
3.10.3 Extend the analysis
- Make a table of mean, standard deviation, skewness, kurtosis, minimum and maximum returns for all ten securities.
- Create ACF plots of squared returns for all securities and compare the strength of volatility clustering.
- Repeat the crisis section for 2008, taking prices from January 2008 to June 2009, and compare the two episodes using normalised prices, volatility and maximum drawdown.
- Plot 60-day rolling skewness and kurtosis for one security and identify periods in which they change sharply.