Understanding Returns and Assessing Risks with Value at Risk

Before choosing portfolio weights, we need a consistent way to describe investment performance. How much did an investment earn? How did its value evolve? How variable were its returns, and how severe were its losses? These questions require different measurements.

This chapter begins with returns and the growth of invested wealth, then develops measures of variability, risk-adjusted performance, drawdown, and tail losses. The aim is to understand what each measure calculates and how to interpret it before applying it to historical data.

In the opening examples, returns are stored as decimals: 0.04 means a return of 4%. Observations are equally spaced, distributions are reinvested when measuring compounded total returns, and the investor makes no additional deposits or withdrawals. Prices, returns, and wealth are distinct quantities; we will keep their units explicit.

The first examples are self-contained and can be run independently. Later examples that build on an earlier dataset should be run in chapter order. Keep full precision in calculations and round only the displayed results.

Analyzing Returns

Percentage Returns Explained

Suppose an asset has price $P_t>0$ at the beginning of a period and price $P_{t+1}$ at its end. Its simple price return, denoted by $r_{t+1}$, measures the price change relative to the capital initially invested:

$$r_{t+1}=\frac{P_{t+1}-P_t}{P_t}=\frac{P_{t+1}}{P_t}-1.$$

Rearranging gives $P_{t+1}=P_t(1+r_{t+1})$. The quantity $1+r_{t+1}$ is the gross return, or the multiplier that converts beginning-of-period value into end-of-period value.

For a price increase from $100 to $104, the dollar gain is $4 and the return is $4/100=0.04=4\mathrm{\%}$. A return is a dimensionless ratio; it is not a dollar amount.

If the asset also pays a cash distribution $D_{t+1}$ at the end of the period, its total return includes that payment:

$$r_{t+1}^{\mathrm{total}}=\frac{P_{t+1}-P_t+D_{t+1}}{P_t}.$$

For the same price change and a $2 dividend, total return is $(104-100+2)/100=6\mathrm{\%}$. The example assumes one share and an end-of-period payment. When using a price series already adjusted to represent reinvested distributions, adding those distributions again would double-count them.

initial_price = 100.0
final_price = 104.0
cash_distribution = 2.0

price_return = (final_price - initial_price) / initial_price
total_return = (final_price - initial_price + cash_distribution) / initial_price

print(f"Price return: {price_return:.2%}")
print(f"Total return: {total_return:.2%}")

Try it: change final_price to 98.0 while keeping the distribution at 2.0. The price return becomes −2%, while total return is zero: the distribution offsets the price loss.

Understanding Compound Returns

How do returns combine across time? Let $V_0$ be initial invested wealth, and let $r_1$ and $r_2$ be the returns in two consecutive periods. With all proceeds reinvested and no external cash flows:

$$V_1=V_0(1+r_1),V_2=V_1(1+r_2)=V_0(1+r_1)(1+r_2).$$

The second return applies to the wealth available after the first period. Dividing final wealth by initial wealth therefore gives the cumulative return:

$$R_{0,2}=\frac{V_2}{V_0}-1=(1+r_1)(1+r_2)-1=r_1+r_2+r_1{r}_2.$$

The term $r_1{r}_2$ explains why adding simple returns does not generally give cumulative performance. Over $n$ periods:

$$V_n=V_0\prod _{t=1}^n(1+r_t),R_{0,n}=\prod _{t=1}^n(1+r_t)-1.$$

Only when the return is the same value $r$ in every period does this simplify to $R_{0,n}=(1+r{)}^n-1$.

Worked Example: A Gain Followed by a Loss

Start with $100, earn 10% in the first month, and lose 10% in the second:

$$V_1=100(1.10)=110,V_2=110(0.90)=99.$$

The investment loses 1% overall. Equal percentage gains and losses do not cancel because they apply to different amounts of capital.

import pandas as pd

monthly_returns = pd.Series(
    [0.10, -0.10],
    index=pd.period_range("2024-01", periods=2, freq="M"),
)
initial_wealth = 100.0
wealth = initial_wealth * (1 + monthly_returns).cumprod()
cumulative_return = (1 + monthly_returns).prod() - 1

print(pd.DataFrame({
    "Monthly return (%)": monthly_returns * 100,
    "End-of-month wealth ($)": wealth,
}).to_string(float_format=lambda value: f"{value:.2f}"))
print(f"Cumulative return: {cumulative_return:.2%}")

cumprod() keeps each intermediate gross-return product, giving the wealth path. prod() returns the final product, from which subtracting one gives the cumulative return. The toolkit functions pok.compound_returns() and pok.compound() implement these two calculations.

Try it: replace the second return with -0.03. Final wealth becomes $106.70 and cumulative return becomes 6.70%.

Arithmetic and Geometric Average Returns

An average return can answer two different questions. The arithmetic mean describes the average one-period observation:

$$\overline{r}=\frac{1}{n}\sum _{t=1}^n{r}_t.$$

It is a sample estimate of the expected one-period return when the observations are representative of the distribution we want to estimate. It does not describe the compounded growth achieved over the full sample.

The geometric mean, $g$, is the constant per-period return that would produce the same final wealth as the observed sequence. Solve:

$$(1+g{)}^n=\prod _{t=1}^n(1+r_t)⟹g={[\prod _{t=1}^n(1+r_t)]}^{1/n}-1.$$

For the two-month example, the arithmetic mean is zero, while the geometric mean is $\sqrt{0.99}-1\approx -0.5013\mathrm{\%}$ per month. Compounding that constant monthly loss twice recovers the 1% cumulative loss.

import pandas as pd

monthly_returns = pd.Series([0.10, -0.10])
arithmetic_mean = monthly_returns.mean()
geometric_mean = (1 + monthly_returns).prod() ** (1 / monthly_returns.size) - 1

print(f"Arithmetic mean per month: {arithmetic_mean:.4%}")
print(f"Geometric mean per month: {geometric_mean:.4%}")

For positive gross returns, the geometric mean cannot exceed the arithmetic mean; they are equal when all period returns are identical. This distinction will matter when we separate historical growth from expected-return inputs in portfolio optimization.

Monthly & Annual Returns

To compare growth measured over different sample lengths, we can express it as an equivalent annual rate. This changes the reporting horizon; it does not predict next year's performance.

Generalizing Annualized Returns

Let $n$ be the number of observed return periods and $p$ the number of those periods in one year. The sample spans $n/p$ years. The annualized geometric return, $R_{\mathrm{ann}}$, must satisfy:

$$(1+R_{\mathrm{ann}}{)}^{n/p}=1+R_{0,n}.$$

Solving for the annual rate gives:

$$R_{\mathrm{ann}}=(1+R_{0,n}{)}^{p/n}-1={[\prod _{t=1}^n(1+r_t)]}^{p/n}-1=(1+g{)}^p-1.$$

Observation frequency	Periods per year, $p$
Monthly	12
Quarterly	4
Weekly	52
Daily trading observations	252, by convention

For example, a cumulative gain of 21% over 24 monthly observations corresponds to ${1.21}^{12/24}-1=10\mathrm{\%}$ per year. A constant quarterly return of 1% corresponds to ${1.01}^4-1\approx 4.06\mathrm{\%}$ per year.

When the sample spans multiple years, this equivalent annual growth rate is commonly called the compound annual growth rate (CAGR). Annualizing a shorter sample uses the same formula but extrapolates its observed growth to a full year. Irregularly spaced observations require the actual elapsed time rather than an assumed number of equally spaced periods.

From Prices to Annualized Returns in Python

Thirteen consecutive month-end prices contain twelve monthly return intervals. The first price has no previous observation, so the return calculation produces an initial NaN. That entry must not count as an additional month of investment performance.

import numpy as np
import pandas as pd
import PortfolioOptimizationKit as pok

monthly_prices = pd.Series(
    100.0 * 1.01 ** np.arange(13),
    index=pd.period_range("2023-12", periods=13, freq="M"),
)
monthly_returns = pok.compute_returns(monthly_prices)
observed_returns = monthly_returns.dropna()
periods_per_year = 12

growth_factor = (1 + observed_returns).prod()
annualized_return = growth_factor ** (periods_per_year / observed_returns.size) - 1

print(f"Price observations: {monthly_prices.size}")
print(f"Observed monthly returns: {observed_returns.size}")
print(f"Annualized return: {annualized_return:.2%}")
print(f"Toolkit annualized return: {pok.annualize_rets(monthly_returns, periods_per_year):.2%}")

Both annualized-return calculations give approximately 12.68%. The toolkit counts non-missing returns separately for each Series, including each column of a DataFrame. An empty or entirely missing return series has no measurable growth rate and returns NaN.

The initial NaN is a consequence of converting prices to returns. Missing observations inside a market history require separate investigation: dropping them does not reconstruct the missing investment periods. Establish the frequency and continuity of the data before interpreting an observed-period growth rate as full-history performance.

Try it: change np.arange(13) to np.arange(25) and the index length to periods=25. You now observe two years of the same 1% monthly growth. Cumulative return increases, while annualized return remains 12.68%.

Assessing Volatility and Risk

Two investments can have similar compounded returns while following very different paths. Volatility measures the dispersion of period returns around their mean. It is one measure of risk; it does not describe the direction of returns or the full severity of potential losses.

Estimating One-Period Volatility

For $n>1$ observed returns, let $\overline{r}$ be their arithmetic mean. The sample standard deviation is:

$${\hat{\sigma }}_{\mathrm{period}}=\sqrt{\frac{1}{n-1}\sum _{t=1}^n(r_t-\overline{r}{)}^2}.$$

Squaring deviations makes both positive and negative departures from the mean contribute to variability. Taking the square root restores the units of returns.

The denominator $n-1$ accounts for estimating the mean from the same sample. It gives an unbiased estimator of variance under independent, identically distributed sampling; its square root is not generally an unbiased estimator of standard deviation. Pandas uses this convention with std(ddof=1). Setting ddof=0 instead divides by $n$, which is useful when describing the supplied observations as the complete population of interest.

Why Volatility Scales with the Square Root of Time

For period returns with a common variance ${\sigma }^2$ and zero covariance across different periods, the variance of their sum is:

$$\mathrm{Var}\left(\sum _{t=1}^p{r}_t\right)=\sum _{t=1}^p\mathrm{Var}(r_t)=p{\sigma }^2.$$

Taking the square root motivates the conventional annualization rule:

$${\hat{\sigma }}_{\mathrm{ann}}=\sqrt{p}{\hat{\sigma }}_{\mathrm{period}}.$$

Monthly returns use $\sqrt{12}$, weekly returns use $\sqrt{52}$, and daily trading returns conventionally use $\sqrt{252}$. The multiplier depends on the frequency of each return, not the number of observations used to estimate volatility. For example, monthly volatility of 4%, estimated from 60 monthly returns, annualizes to $4\mathrm{\%}\sqrt{12}\approx 13.86\mathrm{\%}$.

This derivation concerns a sum of returns. Annual simple returns compound, so square-root scaling is a reporting convention rather than an exact identity for the volatility of compounded annual simple returns. Log returns add across time, but their variance still follows the simple scaling rule only when the covariance and variance assumptions hold. Serial dependence introduces covariance terms; changing volatility also affects the interpretation of a constant annualized estimate.

import pandas as pd
import PortfolioOptimizationKit as pok

monthly_returns = pd.Series([0.02, -0.02, 0.02, -0.02])
periods_per_year = 12
monthly_volatility = monthly_returns.std(ddof=1)
annualized_volatility = pok.annualize_vol(
    monthly_returns, periods_per_year=periods_per_year, ddof=1
)

print(f"Observed months: {monthly_returns.size}")
print(f"Monthly sample volatility: {monthly_volatility:.2%}")
print(f"Annualized sample volatility: {annualized_volatility:.2%}")

The arithmetic mean is zero, and the sum of squared deviations is $4(0.02{)}^2=0.0016$. Monthly sample volatility is therefore $\sqrt{0.0016/3}\approx 2.31\mathrm{\%}$, which annualizes to 8.00%. This short synthetic sample makes the calculation checkable; estimating market risk requires a justified data sample.

Try it: change ddof=1 to ddof=0 in both calculations. The results become 2.00% monthly and approximately 6.93% annualized. The observations have not changed; the variance-estimation convention has.

Zero Volatility Concept

Consider two hypothetical assets over twelve months. Asset A loses exactly 1% every month, while Asset B gains exactly 1% every month. Which has greater return volatility?

Both have zero volatility: every return equals that asset's own mean, so every deviation in the variance formula is zero. Their wealth outcomes are very different:

$$R_A={0.99}^{12}-1\approx -11.36\mathrm{\%},R_B={1.01}^{12}-1\approx 12.68\mathrm{\%}.$$

import pandas as pd
import PortfolioOptimizationKit as pok

monthly_returns = pd.DataFrame({
    "Asset A": [-0.01] * 12,
    "Asset B": [0.01] * 12,
})
summary = pd.DataFrame({
    "Cumulative return (%)": pok.compound(monthly_returns) * 100,
    "Annualized volatility (%)": pok.annualize_vol(monthly_returns, 12) * 100,
})

print(summary.to_string(float_format=lambda value: f"{value:.2f}"))

Volatility describes return variability, not whether the investment preserves capital. A constant sequence of losses can have zero measured volatility. Historical zero volatility also does not establish that future returns are certain. This is why we will complement volatility with drawdown and downside-risk measures later in the chapter.

Python Example: Analyzing Stock Data

We can now combine these measurements in a comparison of two hypothetical stocks. Both begin at $100 and pay no dividends. We specify six monthly returns, construct the corresponding prices by compounding, and then recover the returns from those prices. The data are deliberately simple: Stock B's return is twice Stock A's return in every month.

This gives us a controlled way to compare growth and variability. The observations below are synthetic; the calendar dates label the periods rather than identify real market events. Run this example before the return-on-risk and Sharpe examples, which use its stock_returns DataFrame.

import json
import pandas as pd
import PortfolioOptimizationKit as pok

scenario_returns = pd.DataFrame({
    "Stock A": [0.03, -0.01, 0.02, 0.00, 0.04, -0.02],
    "Stock B": [0.06, -0.02, 0.04, 0.00, 0.08, -0.04],
}, index=pd.period_range("2024-01", periods=6, freq="M"))

initial_prices = pd.DataFrame(
    100.0, index=[scenario_returns.index[0] - 1], columns=scenario_returns.columns
)
stock_prices = pd.concat([initial_prices, 100.0 * (1 + scenario_returns).cumprod()])
stock_returns = pok.compute_returns(stock_prices).dropna()
periods_per_year = 12

summary = pd.DataFrame({
    "Cumulative (%)": pok.compound(stock_returns) * 100,
    "Ann. geometric (%)": pok.annualize_rets(stock_returns, periods_per_year) * 100,
    "Mean monthly (%)": stock_returns.mean() * 100,
    "Ann. volatility (%)": pok.annualize_vol(stock_returns, periods_per_year) * 100,
})
print(summary.to_string(float_format=lambda value: f"{value:.2f}"))

# Include the initial price, but leave its unobserved return blank in the chart.
plot_data = {
    "dates": stock_prices.index.astype(str).tolist(),
    "stockPrices": {
        "title": "Synthetic Stock Prices",
        "series": {name: stock_prices[name].tolist() for name in stock_prices},
        "type": "line",
        "yAxisName": "Price ($)",
    },
    "returns": {
        "title": "Monthly Stock Returns",
        "series": {
            name: [None] + (stock_returns[name] * 100).tolist()
            for name in stock_returns
        },
        "type": "bar",
        "yAxisName": "Monthly return (%)",
    },
}
print("\n<ECHARTS_DATA>" + json.dumps(plot_data, allow_nan=False))

Stock A ends at approximately $106.01 and Stock B at $112.01. Their cumulative returns are 6.01% and 12.01%, while their annualized sample volatilities are approximately 8.20% and 16.40%. Stock B has twice the arithmetic mean return and twice the volatility, but its compounded return is not exactly twice Stock A's: compounding includes products of period returns.

The line chart shows price paths; the bars show returns for individual months. The initial price has no return observation, so its bar is blank rather than an artificial zero. The annualized geometric returns extrapolate six months of observed growth to one year; they are not forecasts.

Try it: change Stock B's final return from -0.04 to -0.08 and rerun the three examples. Observe how the larger loss affects compounded growth, volatility, and the two ratios below.

Evaluating Return on Risk

Does a higher return compensate for higher variability? As a first comparison, define return on risk (ROR) here as the arithmetic mean return divided by return volatility. Both inputs must refer to the same observation frequency:

$${\widehat{\mathrm{ROR}}}_{\mathrm{period}}=\frac{\overline{r}}{{\hat{\sigma }}_{\mathrm{period}}}.$$

Using $p$ periods per year, the conventional annualized arithmetic mean is $p\overline{r}$ and annualized volatility is $\sqrt{p}{\hat{\sigma }}_{\mathrm{period}}$. Therefore:

$${\widehat{\mathrm{ROR}}}_{\mathrm{ann}}=\frac{p\overline{r}}{\sqrt{p}{\hat{\sigma }}_{\mathrm{period}}}=\sqrt{p}\frac{\overline{r}}{{\hat{\sigma }}_{\mathrm{period}}}.$$

This arithmetic annualization is different from the geometric growth rate calculated earlier. It gives a frequency-scaled mean for the ratio, not the investment's realized annual compounded return. Dividing a six-month cumulative return by monthly volatility would mix horizons and change the meaning of the comparison.

In our example, Stock A's monthly mean is 1% and its sample volatility is approximately 2.3664%. Its annualized ROR is:

$${\widehat{\mathrm{ROR}}}_{A,\mathrm{ann}}=\frac{12(0.01)}{\sqrt{12}(0.023664)}\approx 1.464.$$

Stock B has twice the mean and twice the volatility, so its ratio is the same. These ratios are reported as numbers such as 1.464, not as percentages.

import PortfolioOptimizationKit as pok

# Run the stock-data example first to create stock_returns.
periods_per_year = 12
annualized_mean = stock_returns.mean() * periods_per_year
annualized_volatility = pok.annualize_vol(stock_returns, periods_per_year)
return_on_risk = annualized_mean / annualized_volatility.where(annualized_volatility > 0)

print("Annualized return-on-risk ratios:")
for name, value in return_on_risk.items():
    print(f"{name}: {value:.3f}")

The comparison describes this sample using total volatility as the risk measure. It does not account for what an investor could have earned without taking the stock's risk. To include that opportunity cost, we introduce a risk-free benchmark.

Sharpe Ratio: Assessing Risk-Adjusted Returns

The Sharpe ratio measures the arithmetic mean excess return per unit of excess-return volatility. Let $r_{f,t}$ be the risk-free return over the same period as the asset return $r_t$, and define:

$$x_t=r_t-r_{f,t},\overline{x}=\frac{1}{n}\sum _{t=1}^n{x}_t,{\widehat{\mathrm{SR}}}_{\mathrm{period}}=\frac{\overline{x}}{{\hat{\sigma }}_x}.$$

Here ${\hat{\sigma }}_x$ is the sample standard deviation of the excess returns. Subtracting a constant risk-free return changes the mean but not the volatility. With a time-varying benchmark, calculate the volatility of the aligned excess-return series itself.

Matching the Risk-Free Rate to the Observation Period

Our examples use a constant effective annual risk-free rate $R_f$. To obtain the equivalent monthly rate, solve $(1+r_f{)}^{12}=1+R_f$. More generally:

$$r_f=(1+R_f{)}^{1/p}-1.$$

At $R_f=3\mathrm{\%}$ and $p=12$, the monthly risk-free return is approximately 0.2466%. Subtract this monthly return from each monthly stock return, rather than subtracting 3% from a six-month cumulative gain.

Under the constant-variance and zero-serial-covariance assumptions discussed earlier, the conventional annualized Sharpe ratio is:

$${\widehat{\mathrm{SR}}}_{\mathrm{ann}}=\sqrt{p}\frac{\overline{x}}{{\hat{\sigma }}_x}=\frac{p\overline{x}}{\sqrt{p}{\hat{\sigma }}_x}.$$

The numerator uses an arithmetic mean of period excess returns. Compounding the excess returns, or substituting CAGR, produces a different statistic.

Calculating and Interpreting Sharpe Ratios

import numpy as np
import pandas as pd
import PortfolioOptimizationKit as pok

# Run the stock-data example first; the return-on-risk cell is optional here.
periods_per_year = 12
risk_free_rate = 0.03
period_risk_free = (1 + risk_free_rate) ** (1 / periods_per_year) - 1
excess_returns = stock_returns - period_risk_free

excess_volatility = excess_returns.std(ddof=1)
sharpe = np.sqrt(periods_per_year) * excess_returns.mean() / excess_volatility.where(
    excess_volatility > 0
)
toolkit_sharpe = pok.sharpe_ratio(stock_returns, risk_free_rate, periods_per_year)

print(f"Monthly risk-free return: {period_risk_free:.4%}")
print(pd.DataFrame({
    "Annualized Sharpe": sharpe,
    "Toolkit Sharpe": toolkit_sharpe,
}).to_string(float_format=lambda value: f"{value:.3f}"))

Both calculations give approximately 1.103 for Stock A and 1.283 for Stock B. Before subtracting the benchmark, the stocks had identical ROR values. Subtracting the same positive risk-free return reduces Stock A's ratio more because that deduction is larger relative to its volatility.

A positive Sharpe ratio means the sample's arithmetic mean return exceeded the benchmark; a negative ratio means it fell short. A higher value indicates more average excess return per unit of measured variability in that sample. Compare assets over a common date range and with the same benchmark and annualization convention. A short sample or substantial serial dependence can make the estimate a poor guide to future performance.

The toolkit accepts a constant annual risk-free rate for its return-series calculation. It returns NaN when fewer than two returns are observed or all observed returns are identical, because a finite sample Sharpe ratio cannot be estimated in those cases.

Try it: set risk_free_rate = 0.0. Both Sharpe ratios become approximately 1.464, matching the return-on-risk ratios. Increasing the risk-free rate lowers both Sharpe ratios without changing either stock's realized return path.

Illustrating Financial Concepts Using a Real-World Dataset

We now apply the same measurements to historical research portfolios from the Kenneth French Data Library: Portfolios Formed on Size. The repository includes a fixed snapshot of monthly returns from July 1926 to December 2018 in Portfolios_Formed_on_ME_monthly_EW.csv.

The underlying portfolios are formed using market equity and NYSE size breakpoints. We use the equally weighted returns of two size groups:

• Lo 10: the smallest size decile, referred to here as small caps.
• Hi 10: the largest size decile, referred to here as large caps.

Equal weighting applies to the stocks within each research portfolio. These columns already contain portfolio returns, not prices. The raw file expresses them in percent, so 3.29 must become 0.0329 before compounding. The toolkit loader performs that conversion and creates a monthly PeriodIndex.

Data Analysis

Begin by checking the dates, observation count, and missing values. A missing market observation is not a zero return. For this comparison, both portfolios must have observations for the same consecutive months.

For the chart, combine monthly returns into calendar-year compounded returns:

$$R_y=\prod _{m=1}^{12}(1+r_{y,m})-1,$$

where $r_{y,m}$ is the return in month $m$ of year $y$. The arithmetic mean of those twelve returns would instead describe an average month. Because the snapshot starts in July, 1926 is a partial year and is excluded from this annual chart. Its six observed months remain available for the full-history statistics.

import json
import pandas as pd
import PortfolioOptimizationKit as pok

small_large_caps = pok.get_ffme_returns()
expected_months = pd.period_range(
    small_large_caps.index[0], small_large_caps.index[-1], freq="M"
)
if not small_large_caps.index.equals(expected_months):
    raise ValueError("Expected unique, consecutive monthly observations.")
if small_large_caps.isna().any().any():
    raise ValueError("Resolve missing portfolio returns before this comparison.")

print(f"Coverage: {small_large_caps.index[0]} to {small_large_caps.index[-1]}")
print(f"Monthly observations: {len(small_large_caps)}")
print("First monthly returns (%):")
print((small_large_caps.head() * 100).to_string(float_format=lambda value: f"{value:.2f}"))

year_groups = small_large_caps.groupby(small_large_caps.index.year)
complete_years = year_groups.count().eq(12).all(axis=1)
yearly_returns = year_groups.agg(lambda returns: (1 + returns).prod() - 1)
yearly_returns = yearly_returns.loc[complete_years]

plot_data = {
    "dates": yearly_returns.index.astype(str).tolist(),
    "smallLargeCaps": {
        "title": "Small- and Large-Cap Calendar-Year Returns",
        "series": {
            "Small caps (Lo 10)": (yearly_returns["Lo 10"] * 100).tolist(),
            "Large caps (Hi 10)": (yearly_returns["Hi 10"] * 100).tolist(),
        },
        "type": "bar",
        "yAxisName": "Calendar-year return (%)",
    },
}
print("\n<ECHARTS_DATA>" + json.dumps(plot_data, allow_nan=False))

The snapshot contains 1,110 monthly observations and 92 complete calendar years, from 1927 through 2018. For example, the compounded 1927 returns are approximately 45.82% for small caps and 29.12% for large caps. The chart expresses these returns as percentages, while the DataFrame retains decimal values for subsequent calculations.

Try it: replace the loading expression with pok.get_ffme_returns().loc["2000":] and rerun the following statistics. This changes the estimation sample while keeping the return frequency monthly.

Calculating Volatility

Use the monthly observations to estimate sample volatility, then multiply by $\sqrt{12}$. The 1,110 observations determine the estimation window; twelve months determine the annualization factor.

# Run the historical-data example first.
monthly_volatility = small_large_caps.std(ddof=1)
annualized_volatility = pok.annualize_vol(small_large_caps, periods_per_year=12)

volatility_comparison = pd.DataFrame({
    "Monthly volatility (%)": monthly_volatility * 100,
    "Annualized volatility (%)": annualized_volatility * 100,
})
print(volatility_comparison.to_string(float_format=lambda value: f"{value:.2f}"))

Over the full snapshot, annualized volatility is approximately 36.82% for small caps and 18.67% for large caps. Small-cap returns were more variable in this sample. These full-history estimates average across very different market conditions; they do not assert that volatility remained constant throughout the period.

Returns Analysis

Now distinguish the arithmetic mean of monthly returns from the geometric growth achieved by compounding them. With $n$ observed months, the equivalent monthly growth rate is $g=(1+R_{0,n}{)}^{1/n}-1$, and the annualized geometric return is $(1+g{)}^{12}-1$.

total_return = pok.compound(small_large_caps)
return_per_month = (1 + total_return) ** (1 / small_large_caps.count()) - 1
annualized_return = pok.annualize_rets(small_large_caps, periods_per_year=12)

return_comparison = pd.DataFrame({
    "Mean monthly (%)": small_large_caps.mean() * 100,
    "Geometric monthly (%)": return_per_month * 100,
    "Annualized geometric (%)": annualized_return * 100,
})
print(return_comparison.to_string(float_format=lambda value: f"{value:.2f}"))

The annualized geometric returns are approximately 16.75% for small caps and 9.28% for large caps. These are growth rates for the supplied research portfolios over the full sample. An investor's realized experience would also depend on implementation costs, taxes, and the ability to trade the constituent stocks.

For risk-adjusted performance, apply the arithmetic-mean ROR and excess-return Sharpe conventions introduced above. The next example uses a constant 3% effective annual risk-free rate as an illustrative benchmark. Actual Treasury returns varied substantially over this history; a historical benchmark study would use a period-matched risk-free series.

# Assuming a risk-free rate
risk_free_rate = 0.03

# Annualized arithmetic mean divided by annualized volatility
return_on_risk = (small_large_caps.mean() * 12) / annualized_volatility.where(annualized_volatility > 0)
sharpe_ratio = pok.sharpe_ratio(small_large_caps, risk_free_rate, periods_per_year=12)

print("\nReturn on Risk:")
for col in return_on_risk.index:
    print(f"{col}: {return_on_risk[col]:.3f}")

print("\nSharpe Ratio:")
for col in sharpe_ratio.index:
    print(f"{col}: {sharpe_ratio[col]:.3f}")

With that assumed benchmark, the annualized Sharpe ratios are approximately 0.493 for small caps and 0.412 for large caps. The higher small-cap growth rate came with higher volatility, and the Sharpe comparison depends on the sample and benchmark. Neither measure tells us how far wealth fell below a previous peak. That requires drawdown analysis.

Drawdown

Volatility measures the dispersion of returns. Drawdown answers a different question: how far is the investment's current wealth below its highest value so far?

Let $V_0>0$ be initial wealth and $r_t$ the return in period $t$. With reinvestment and no external cash flows, define the wealth path, running peak, and drawdown as:

$$V_t=V_0\prod _{i=1}^t(1+r_i),H_t=max\{V_0,V_1,\dots ,V_t\},D_t=\frac{V_t-H_t}{H_t}=\frac{V_t}{H_t}-1.$$

The running peak $H_t$ is also called the high-water mark. It includes initial wealth: if the first return is negative, the portfolio is already in drawdown. We report $D_t$ as a non-positive number, with zero indicating that wealth is at its running peak.

For observations through time $T$, the maximum drawdown as a positive loss magnitude is:

$${\mathrm{MDD}}_{0,T}=-\underset{0\le t\le T}{min}{D}_t=\underset{0\le t\le T}{max}(1-\frac{V_t}{H_t}).$$

This is a relative peak-to-trough loss. The dollar difference $H_t-V_t$ is an absolute loss and changes when initial capital changes. Some reports use the signed minimum drawdown instead; here the charts show negative drawdowns, while the maximum-drawdown tables show positive loss magnitudes.

Worked Example: Loss, Recovery, and a New Decline

Start with $100 and observe monthly returns of −20%, +25%, −10%, and +5%. Wealth evolves as:

$$100⟶80⟶100⟶90⟶94.50.$$

The high-water mark remains $100. Drawdowns are therefore −20%, 0%, −10%, and −5.5%, giving a maximum drawdown of 20%.

import pandas as pd
import PortfolioOptimizationKit as pok

example_returns = pd.Series(
    [-0.20, 0.25, -0.10, 0.05],
    index=pd.period_range("2024-01", periods=4, freq="M"),
)
example_drawdown = pok.drawdown(example_returns, start=100.0)
display_table = example_drawdown.copy()
display_table["Drawdown"] *= 100
display_table = display_table.rename(columns={"Drawdown": "Drawdown (%)"})

print("Initial wealth: $100.00")
print(display_table.to_string(float_format=lambda value: f"{value:.2f}"))
print(f"Maximum drawdown: {-example_drawdown['Drawdown'].min():.2%}")

A 20% loss requires a 25% gain to recover: after losing a fraction $0\le d<1$ of peak wealth, the required gain is $1/(1-d)-1=d/(1-d)$. The trough is the low point; recovery is the first later observation at which wealth regains the preceding peak. A sample can end before recovery occurs.

Try it: reorder the returns to [-0.20, -0.10, 0.05, 0.25]. Final wealth remains $94.50, but the maximum drawdown increases to 28%. Compounded terminal wealth is unchanged by reordering the same returns, while the drawdown depends on their sequence.

Practical example

Apply the same calculation to the historical portfolios. Run the historical-data example first so that small_large_caps contains the chosen monthly sample. The toolkit returns end-of-period wealth, peaks, and signed drawdowns for each Series; its peak calculation includes starting capital.

The wealth charts use a logarithmic vertical scale because the histories span many decades. Equal proportional changes then occupy equal vertical distances. The drawdown chart uses a linear percentage scale.

import json
import pandas as pd
import PortfolioOptimizationKit as pok

start_value = 100.0
rets = small_large_caps.rename(columns={"Lo 10": "Small Caps", "Hi 10": "Large Caps"})
historical_drawdowns = {
    name: pok.drawdown(rets[name], start=start_value) for name in rets
}
drawdowns = pd.DataFrame({
    name: result["Drawdown"] for name, result in historical_drawdowns.items()
})
drawdown_summary = pd.DataFrame({
    "Max drawdown (%)": -drawdowns.min() * 100,
    "Trough month": drawdowns.idxmin().astype(str),
})
print(drawdown_summary.to_string(float_format=lambda value: f"{value:.2f}"))

# Show starting capital at the month-end immediately before the first return.
chart_dates = [str(rets.index[0] - 1)] + rets.index.astype(str).tolist()
events = [
    {"date": "1929-10", "name": "1929 crash"},
    {"date": "2000-03", "name": "Dot-com peak"},
    {"date": "2008-09", "name": "Lehman failure"},
]
plot_data = {
    "dates": chart_dates,
    "crisis": [event for event in events if event["date"] in chart_dates],
}
for name, result in historical_drawdowns.items():
    plot_data[f"{name} wealth"] = {
        "title": f"{name}: Wealth and Running Peak",
        "series": {
            "Wealth": [start_value] + result["Wealth"].tolist(),
            "Running peak": [start_value] + result["Peaks"].tolist(),
        },
        "type": "line",
        "yAxis": {"type": "log", "name": "Wealth ($, log scale)"},
    }
plot_data["drawdown"] = {
    "title": "Historical Drawdowns",
    "series": {name: [0.0] + (drawdowns[name] * 100).tolist() for name in drawdowns},
    "type": "line",
    "yAxisName": "Drawdown (%)",
}
print("\n<ECHARTS_DATA>" + json.dumps(plot_data, allow_nan=False))

For the full snapshot, maximum drawdowns are approximately 83.30% for small caps and 84.00% for large caps, both reaching their trough in May 1932. These are losses from the running peaks, not losses measured from the original $100 investment.

Try it: change start_value to 1000.0. Wealth and peak values increase tenfold, while drawdowns and trough dates remain the same. Percentage drawdown is independent of the investment's scale.

Insights from Historical Crises

To examine particular episodes, select windows from the drawdown series already calculated. We use 1929–1933 for the Great Depression, 2000–2002 for the dot-com decline, and 2007–2009 for the global financial crisis.

The running peaks still come from the entire chosen return history. A peak can therefore precede the displayed window. This measures the worst drawdown observed during the window, which differs from restarting an investment at the window's first date.

# Run the historical drawdown example first.
windows = {
    "Great Depression (1929-1933)": ("1929-01", "1933-12"),
    "Dot-com decline (2000-2002)": ("2000-01", "2002-12"),
    "Global financial crisis (2007-2009)": ("2007-01", "2009-12"),
}
crisis_records = []
for label, (first_month, last_month) in windows.items():
    window = drawdowns.loc[first_month:last_month]
    if window.empty:
        continue
    for name in window:
        crisis_records.append({
            "Window": label,
            "Portfolio": name,
            "Worst drawdown (%)": -window[name].min() * 100,
            "Trough month": str(window[name].idxmin()),
        })

crisis_summary = pd.DataFrame(crisis_records, columns=[
    "Window", "Portfolio", "Worst drawdown (%)", "Trough month"
])
print(crisis_summary.to_string(index=False, float_format=lambda value: f"{value:.2f}"))

For the full-history calculation, the dot-com window reaches its worst observations at approximately 36.34% for small caps in December 2000 and 49.52% for large caps in September 2002. During the global financial crisis window, the corresponding losses are approximately 63.12% and 52.81%, both in February 2009. These dates identify troughs, not recoveries.

The episode labels provide historical context; they do not establish that one event caused every loss in the window. Monthly observations can also miss deeper intramonth drawdowns. Together, growth, volatility, and drawdown describe different aspects of the investment experience. Next we examine the shape of the return distribution and the information contained in its tails.

Gaussian Density & Distribution

Return and volatility estimates summarize a sample, but a distribution model lets us ask probability questions. For example: if monthly returns followed a specified distribution, what would be the probability of losing more than 5% in a month?

The Gaussian, or normal, distribution provides a useful reference model. Write:

$$X\sim \mathcal{N}(\mu ,{\sigma }^2),\sigma >0,$$

where $\mu$ is the mean, $\sigma$ is the standard deviation, and ${\sigma }^2$ is the variance. For a monthly return model, both $\mu$ and $\sigma$ refer to monthly observations. A 1% monthly mean and 4% monthly volatility correspond to $\mu =0.01$, $\sigma =0.04$, and variance $0.0016$.

Density Is Not a Probability

The probability density function is:

$$f_X(x)=\frac{1}{\sigma \sqrt{2\pi }}\exp [-\frac{1}{2}{\left(\frac{x-\mu }{\sigma }\right)}^2].$$

The standardized distance $(x-\mu )/\sigma$ determines how far $x$ lies from the mean in standard-deviation units. The squared distance makes the density symmetric around $\mu$, while the factor $1/(\sigma \sqrt{2\pi })$ ensures that the total area under the curve is one.

A probability is an area, not the height of the density at one point:

$$\mathbb{P}(a\le X\le b)={\int }_a^b{f}_X(u)du.$$

For a continuous distribution, $\mathbb{P}(X=x)=0$ at any single point. Density values can exceed one; that does not imply a probability above 100%.

Cumulative Probability and Standardization

The cumulative distribution function (CDF) gives the probability of an outcome at or below a threshold:

$$F_X(x)=\mathbb{P}(X\le x)={\int }_{-\mathrm{\infty }}^x{f}_X(u)du.$$

To evaluate it, standardize $X$:

$$Z=\frac{X-\mu }{\sigma }\sim \mathcal{N}(0,1).$$

Subtracting $\mu$ centers the distribution at zero; dividing by $\sigma$ gives unit variance. Denote the standard-normal density by $\varphi$ and its CDF by $\mathrm{\Phi }$. Because $\sigma >0$, standardization preserves the direction of the inequality:

$$F_X(x)=\mathbb{P}(Z\le \frac{x-\mu }{\sigma })=\mathrm{\Phi }\left(\frac{x-\mu }{\sigma }\right).$$

Thus $F_X(x)=\mathrm{\Phi }(x)$ only when $X$ itself is standard normal. Probabilities for intervals follow by subtraction: $\mathbb{P}(a\le X\le b)=F_X(b)-F_X(a)$.

Symmetry Property

The standard-normal density satisfies $\varphi (-z)=\varphi (z)$ because it depends on $z^2$. Equal areas on opposite sides of zero give:

$$\mathrm{\Phi }(-z)=1-\mathrm{\Phi }(z).$$

For example, the probability of being more than one standard deviation below the mean equals the probability of being more than one standard deviation above it: each tail has approximately 15.87% probability. The interval between those thresholds contains approximately 68.27% of the distribution.

Distribution of Negative Random Variable

Reflecting a standard-normal variable around zero leaves its distribution unchanged:

$$\mathbb{P}(-Z\le z)=\mathbb{P}(Z\ge -z)=1-\mathrm{\Phi }(-z)=\mathrm{\Phi }(z).$$

The equality uses continuity, so including or excluding the boundary does not change the probability. More generally:

$$X\sim \mathcal{N}(\mu ,{\sigma }^2)⟹-X\sim \mathcal{N}(-\mu ,{\sigma }^2).$$

The mean changes sign while variance is unchanged. This matters when moving from a return $R$ to a loss $L=-R$: large positive losses correspond to the left tail of returns.

Worked Example: Monthly Loss Probabilities

Assume monthly returns follow $R\sim \mathcal{N}(0.01,{0.04}^2)$. A return of −5% lies 1.5 standard deviations below the mean:

$$z=\frac{-0.05-0.01}{0.04}=-1.5,\mathbb{P}(R\le -0.05)=\mathrm{\Phi }(-1.5)\approx 6.68\mathrm{\%}.$$

from scipy.stats import norm

mu = 0.01
sigma = 0.04
loss_threshold = -0.05

standardized_threshold = (loss_threshold - mu) / sigma
loss_probability = norm.cdf(0.0, loc=mu, scale=sigma)
tail_probability = norm.cdf(loss_threshold, loc=mu, scale=sigma)
within_one_sigma = (
    norm.cdf(mu + sigma, loc=mu, scale=sigma)
    - norm.cdf(mu - sigma, loc=mu, scale=sigma)
)
density_at_mean = norm.pdf(mu, loc=mu, scale=sigma)

print(f"Standardized loss threshold: {standardized_threshold:.2f}")
print(f"Probability of a negative return: {loss_probability:.2%}")
print(f"Probability of return <= {loss_threshold:.2%}: {tail_probability:.2%}")
print(f"Probability within one standard deviation: {within_one_sigma:.2%}")
print(f"Density at the mean (not a probability): {density_at_mean:.4f}")

The probability of a negative return is approximately 40.13%, while the probability of a return at or below −5% is 6.68%. The density at the mean is approximately 9.9736, illustrating why a density height must not be interpreted as a probability. In SciPy, loc is the mean and scale is the standard deviation, not the variance.

Try it: increase sigma to 0.08. The probability of a return at or below −5% rises to approximately 22.66%. The interval defined by one standard deviation around the mean widens, but its probability remains approximately 68.27%.

These are probabilities under the assumed model. Estimating a mean and volatility from market data does not by itself establish that a normal distribution is an adequate description of the returns.

Quantiles

A CDF takes a return threshold and gives a probability. A quantile reverses the question: for a specified probability, what threshold marks that part of the distribution?

For $0<p<1$, a general definition of the $p$-quantile is:

$$q_p(X)=inf\{x:F_X(x)\ge p\}.$$

It is the smallest threshold at which the CDF reaches or exceeds $p$. For a continuous, strictly increasing CDF such as the Gaussian CDF, this simplifies to:

$$q_p(X)=F_X^{-1}(p),\mathbb{P}(X\le q_p(X))=p.$$

The equality need not hold for a discrete distribution because its CDF can jump over the requested probability. That distinction also affects how we estimate quantiles from a finite sample.

Standard-Normal Quantiles

Write $z_p={\mathrm{\Phi }}^{-1}(p)$ for a standard-normal quantile. Symmetry gives:

$$\mathrm{\Phi }(-z_p)=1-\mathrm{\Phi }(z_p)=1-p⟹z_{1-p}=-z_p.$$

For example, $z_{0.95}\approx 1.6449$ and $z_{0.05}\approx -1.6449$. The median is $z_{0.5}=0$.

Practical Application: Finding a Specific Quantile

To find the threshold below which 90% of a standard-normal distribution lies, calculate $z_{0.9}={\mathrm{\Phi }}^{-1}(0.9)$. SciPy calls the inverse CDF the percent point function, available as norm.ppf().

from scipy.stats import norm

p = 0.90
z_quantile = norm.ppf(p)
probability_check = norm.cdf(z_quantile)
complementary_quantile = norm.ppf(1 - p)

print(f"Standard-normal {p:.0%} quantile: {z_quantile:.4f}")
print(f"CDF evaluated at that quantile: {probability_check:.4f}")
print(f"Complementary {1 - p:.0%} quantile: {complementary_quantile:.4f}")

The results are 1.2816, 0.9000, and −1.2816. The two quantiles are distances from the mean in standard-deviation units. They are not return percentages until we specify a return distribution's mean and volatility.

From Standard-Normal Quantiles to Return Thresholds

For $R=\mu +\sigma Z$ with $\sigma >0$, transform the standard-normal quantile back into return units:

$$q_p(R)=\mu +\sigma z_p=\mu +\sigma {\mathrm{\Phi }}^{-1}(p).$$

The reason is that $\mathbb{P}(R\le \mu +\sigma z_p)=\mathbb{P}(Z\le z_p)=p$. With a 1% monthly mean and 4% monthly volatility, the 5% return quantile is $0.01+0.04(-1.6449)\approx -5.58\mathrm{\%}$.

A central probability interval uses a quantile at each end. For desired coverage $c$, allocate $(1-c)/2$ to each tail:

$$\mathbb{P}(q_{(1-c)/2}(R)\le R\le q_{(1+c)/2}(R))=c.$$

For 95% central coverage, the endpoints use the 2.5% and 97.5% quantiles, corresponding to approximately $\mu \pm 1.96\sigma$. Using the 5% and 95% quantiles would instead leave 5% in each tail and give 90% central coverage.

from scipy.stats import norm

mu = 0.01
sigma = 0.04
left_tail_quantile = norm.ppf(0.05, loc=mu, scale=sigma)

central_probability = 0.95
tail_probability = (1 - central_probability) / 2
lower = norm.ppf(tail_probability, loc=mu, scale=sigma)
upper = norm.ppf(1 - tail_probability, loc=mu, scale=sigma)
coverage = norm.cdf(upper, loc=mu, scale=sigma) - norm.cdf(lower, loc=mu, scale=sigma)

print(f"5% monthly return quantile: {left_tail_quantile:.2%}")
print(f"Central {central_probability:.0%} interval: [{lower:.2%}, {upper:.2%}]")
print(f"Probability between the endpoints: {coverage:.2%}")

The central interval is approximately [−6.84%, 8.84%]. It is a model-based range for an individual monthly return, assuming the specified parameters. It is not a confidence interval for an estimated mean, and it does not incorporate uncertainty in parameter estimates.

Try it: set central_probability = 0.90. The lower endpoint becomes the 5% return quantile, approximately −5.58%, and the upper endpoint becomes approximately 7.58%.

Quantiles Estimated from Observed Returns

For historical returns, we can estimate a quantile directly from the ordered observations instead of assuming a Gaussian distribution. The interpolation convention matters, especially with small samples.

Consider ten monthly returns whose two worst observations are −7% and −4%. NumPy's linear method places the 10% quantile at the zero-based position $h=(n-1)p=9(0.10)=0.9$. It interpolates 90% of the way from the first observation to the second:

$${\hat{q}}_{0.10}^{\mathrm{linear}}=(1-0.9)(-0.07)+0.9(-0.04)=-0.043.$$

The inverse empirical CDF method instead selects the smallest observed value at which at least 10% of the sample has accumulated. With ten observations, it selects the worst return, −7%.

import numpy as np

observed_returns = np.array([-0.04, 0.05, 0.02, -0.07, 0.01, 0.005, -0.02, -0.01, -0.02, 0.05])
p = 0.10
linear_quantile = np.quantile(observed_returns, p, method="linear")
empirical_quantile = np.quantile(observed_returns, p, method="inverted_cdf")

print(f"Linear-interpolated {p:.0%} return quantile: {linear_quantile:.2%}")
print(f"Inverse empirical CDF {p:.0%} return quantile: {empirical_quantile:.2%}")

The estimates are −4.30% and −7.00%, respectively. They use the same data but different finite-sample conventions. np.quantile() takes probabilities between zero and one; np.percentile() takes percentiles between zero and 100. Their default interpolation method is linear.

When comparing historical risk estimates, specify the quantile convention as well as the probability and observation frequency. Before using either historical or Gaussian quantiles to estimate tail risk, we next examine how observed return distributions can differ from the Gaussian benchmark.

Exploring Skewness and Kurtosis in Financial Data

Two return distributions can have the same mean and volatility but very different patterns of large gains and losses. Skewness and kurtosis describe aspects of that shape. They help us assess departures from a Gaussian benchmark before choosing how to estimate tail risk.

Skewness: Assessing Asymmetry in Distributions

For a return $R$ with mean $\mu$ and positive standard deviation $\sigma$, skewness is its third standardized central moment, assuming that moment exists:

$$S(R)=\frac{\mathbb{E}[(R-\mu {)}^3]}{{\sigma }^3}=\mathbb{E}\left[{\left(\frac{R-\mu }{\sigma }\right)}^3\right].$$

Centering measures deviations from the mean, and dividing by $\sigma$ expresses them in standard-deviation units. Cubing preserves their signs and gives greater weight to large deviations:

• Negative skewness: large negative standardized deviations dominate the third moment. In returns, this can reflect occasional severe losses alongside more frequent modest gains.
• Positive skewness: large positive standardized deviations dominate the third moment.
• Zero skewness: positive and negative third-moment contributions balance. This alone does not establish symmetry or normality.

Skewness is dimensionless. Converting decimal returns to percentages leaves it unchanged, while multiplying all returns by −1 reverses its sign.

Kurtosis: Analyzing Tails and Outliers

Kurtosis uses the fourth standardized central moment, assuming it exists:

$$K(R)=\frac{\mathbb{E}[(R-\mu {)}^4]}{{\sigma }^4}=\mathbb{E}\left[{\left(\frac{R-\mu }{\sigma }\right)}^4\right].$$

The fourth power makes both large positive and large negative deviations contribute strongly. Kurtosis measures the importance of extreme standardized observations; it is not simply a measure of how sharply peaked a histogram looks. It also does not identify which tail contains the large observations.

A Gaussian distribution has Pearson kurtosis $K=3$. Excess kurtosis subtracts that benchmark:

$$K_{\mathrm{excess}}=K-3.$$

Positive excess kurtosis means a larger fourth standardized moment than the Gaussian benchmark. Both skewness and kurtosis can be sensitive to a few observations, particularly in small samples; neither specifies an entire return distribution.

Estimating the Moments from Returns

For $n$ observed returns, define the sample central moments:

$$m_j=\frac{1}{n}\sum _{t=1}^n(r_t-\overline{r}{)}^j,\hat{S}=\frac{m_3}{m_2^{3/2}},\hat{K}=\frac{m_4}{m_2^2}.$$

The toolkit uses these moment ratios. Their scale is $\sqrt{m_2}$, calculated with ddof=0, so the same denominator $n$ is used throughout. This differs from the ddof=1 volatility estimate used earlier. Pandas' skew() and kurt() apply finite-sample corrections, and kurt() reports excess kurtosis. To match the toolkit in SciPy, use skew(..., bias=True) and kurtosis(..., fisher=False, bias=True) on the same observed values.

Consider returns of −4%, +1%, +1%, +1%, and +1%. Their mean is zero and $\sqrt{m_2}=2\mathrm{\%}$, giving standardized deviations of $(-2,0.5,0.5,0.5,0.5)$. Therefore:

$$\hat{S}=\frac{(-2{)}^3+4(0.5{)}^3}{5}=-1.5,\hat{K}=\frac{(-2{)}^4+4(0.5{)}^4}{5}=3.25.$$

import pandas as pd
import PortfolioOptimizationKit as pok

toy_returns = pd.Series([-0.04, 0.01, 0.01, 0.01, 0.01])
centered = toy_returns - toy_returns.mean()
m2 = (centered ** 2).mean()
sample_skewness = (centered ** 3).mean() / m2 ** 1.5
sample_kurtosis = (centered ** 4).mean() / m2 ** 2

print(f"Skewness: {sample_skewness:.3f} | Toolkit: {pok.skewness(toy_returns):.3f}")
print(f"Pearson kurtosis: {sample_kurtosis:.3f} | Toolkit: {pok.kurtosis(toy_returns):.3f}")
print(f"Excess kurtosis: {pok.exkurtosis(toy_returns):.3f}")

The results are −1.500, 3.250, and 0.250. The toolkit excludes missing observations separately for each Series or DataFrame column. Standardized moments are undefined when there is no variation; the toolkit returns NaN for constant or insufficient data.

Try it: insert toy_returns = -toy_returns after creating the Series. Skewness becomes +1.500, while kurtosis stays 3.250. Reflection changes the direction of asymmetry, not the size of the fourth-moment contributions.

Practical Analysis

Compare the historical large-cap monthly returns with a reproducible standard-normal sample of the same length. Standardize the market returns to remove their mean and scale; this preserves their skewness and kurtosis. Use common histogram bins so the comparison reflects distribution shape rather than different units or bin boundaries.

import json
import numpy as np
import pandas as pd
import PortfolioOptimizationKit as pok

market_returns = pok.get_ffme_returns()["Hi 10"].dropna()
market_z = (market_returns - market_returns.mean()) / market_returns.std(ddof=0)
rng = np.random.default_rng(42)
normal_sample = pd.Series(rng.normal(size=len(market_returns)), index=market_returns.index)
shape_samples = pd.DataFrame({
    "Gaussian sample": normal_sample,
    "Large caps (standardized)": market_z,
})
shape_comparison = pd.DataFrame({
    "Skewness": pok.skewness(shape_samples),
    "Pearson kurtosis": pok.kurtosis(shape_samples),
    "Excess kurtosis": pok.exkurtosis(shape_samples),
})
print(shape_comparison.to_string(float_format=lambda value: f"{value:.3f}"))

bin_edges = np.linspace(shape_samples.min().min(), shape_samples.max().max(), 41)
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
histogram_series = {}
for name in shape_samples:
    density, _ = np.histogram(shape_samples[name], bins=bin_edges, density=True)
    histogram_series[name] = np.column_stack([bin_centers, density]).tolist()

plot_data = {
    "distributionComparison": {
        "title": "Return Shape: Gaussian Sample and Large Caps",
        "type": "bar",
        "xAxis": {"type": "value", "name": "Standardized return"},
        "yAxisName": "Density",
        "series": histogram_series,
    },
}
print("\n<ECHARTS_DATA>" + json.dumps(plot_data, allow_nan=False))

With this seed, the Gaussian sample has skewness approximately −0.041 and kurtosis 3.097. The full-history large-cap returns have skewness approximately 0.233 and kurtosis 10.695. The simulated moments need not equal their population values exactly, while the market sample's large fourth moment indicates a substantial contribution from extreme standardized returns.

Histogram heights are densities: the sum of each height times its bin width is one. Their heights do not need to sum to one. The horizontal axis measures standard-deviation units, not percentage returns.

Try it: change the random seed from 42 to 7. The simulated moments and histogram change, but the historical market moments remain the same. A finite sample from a normal distribution will not generally have skewness exactly zero or kurtosis exactly three.

Comparing Hedge Fund Indices

The bundled edhec-hedgefundindices.csv contains 263 monthly observations for 13 strategy indices, from January 1997 through November 2018. Its dates use day/month/year notation, and its percentage returns are converted to decimals by pok.get_hfi_returns().

These are historical strategy-index returns. We use them to compare observed distribution shapes, rather than treating an index's moments as a forecast for an individual fund.

import pandas as pd
import PortfolioOptimizationKit as pok

hfi = pok.get_hfi_returns()
hfi_skew_kurt = pd.DataFrame({
    "Observations": hfi.count(),
    "Skewness": pok.skewness(hfi),
    "Pearson kurtosis": pok.kurtosis(hfi),
    "Excess kurtosis": pok.exkurtosis(hfi),
})
print(f"Coverage: {hfi.index[0]:%Y-%m} to {hfi.index[-1]:%Y-%m}")
print(hfi_skew_kurt.to_string(float_format=lambda value: f"{value:.3f}"))

CTA Global has skewness approximately 0.174 and kurtosis 2.953. Convertible Arbitrage has skewness approximately −2.640 and kurtosis 23.281. The latter sample has pronounced negative third-moment contributions and large fourth-moment contributions. CTA Global's moments are closer to the Gaussian benchmark, but two moment estimates alone cannot establish that its distribution is Gaussian.

Using the Jarque–Bera Test for Normality

The Jarque–Bera (JB) test combines deviations of sample skewness and kurtosis from their Gaussian values. Using the moment estimators defined above:

$$\mathrm{JB}=\frac{n}{6}[{\hat{S}}^2+\frac{(\hat{K}-3{)}^2}{4}].$$

Under the null hypothesis of independent, identically distributed Gaussian observations, the statistic has an asymptotic chi-squared distribution with two degrees of freedom. Larger values indicate stronger departures in these two moments.

Choose a significance level $\alpha$, here 1%. If the p-value is below $\alpha$, reject the normality null; otherwise, do not reject it. A p-value is a tail probability for the test statistic under the null hypothesis, not the probability that the return distribution is normal. Failure to reject is not proof of normality.

from scipy.stats import jarque_bera

# Run the hedge-fund index example first.
level = 0.01
jb_rows = {}
for name in ["CTA Global", "Convertible Arbitrage"]:
    observed = hfi[name].dropna()
    skew = pok.skewness(observed)
    pearson_kurtosis = pok.kurtosis(observed)
    jb_from_moments = len(observed) / 6 * (skew ** 2 + (pearson_kurtosis - 3) ** 2 / 4)
    statistic, p_value = jarque_bera(observed)
    decision = "Unavailable" if pd.isna(p_value) else ("Reject" if p_value < level else "Do not reject")
    jb_rows[name] = {
        "JB (formula)": jb_from_moments,
        "JB (SciPy)": statistic,
        "p-value": p_value,
        "Decision": decision,
    }

jb_comparison = pd.DataFrame.from_dict(jb_rows, orient="index")
print(f"Significance level: {100 * level:.3g}%")
print(jb_comparison.to_string(float_format=lambda value: f"{value:.4g}"))

The two implementations agree: CTA Global has JB approximately 1.347 and a p-value of 0.510, so normality is not rejected at 1%. Convertible Arbitrage has JB approximately 4,812.703, leading to rejection. Its p-value is so small that SciPy reports zero at floating-point precision.

The chi-squared calibration is a large-sample approximation. With 263 monthly observations, finite-sample behavior deserves attention; serial dependence or return smoothing can also affect the interpretation of the p-values.

Normality Across Indices

Apply the same decision rule separately to each strategy. Despite its name, pok.is_normal() returns True when the test does not reject normality at the chosen level. It returns NaN for an unavailable test, including insufficient or constant observations. For a DataFrame, each column is tested separately rather than pooling different strategies into one distribution.

# Run the hedge-fund and Jarque–Bera examples first.
normality_test_results = pok.is_normal(hfi, level=level)
p_values = hfi.apply(lambda returns: jarque_bera(returns.dropna()).pvalue)
decisions = normality_test_results.map({True: "Do not reject", False: "Reject"}).fillna("Unavailable")
normality_table = pd.DataFrame({
    "Observations": hfi.count(),
    "p-value": p_values,
    "Decision": decisions,
})
print(normality_table.to_string(formatters={"p-value": lambda value: f"{value:.3g}"}))

At the 1% level, CTA Global is the only index in this snapshot for which the test does not reject normality. This is a statement about the data, test, and threshold, not a ranking of which investment is safest or a guarantee of normally distributed future returns. These are individual tests without a multiple-testing adjustment.

Try it: reduce level to 1e-8 in the Jarque–Bera example and rerun both test examples. Long/Short Equity is then also not rejected. A smaller significance level requires stronger evidence to reject the null; the underlying returns and their moments have not changed.

Skewness, kurtosis, and normality tests help diagnose distributional assumptions. Next we focus directly on negative returns and loss thresholds through downside-risk measures.

Understanding Downside Risk Measures

Volatility and distribution moments describe variability around the mean. Downside-risk measures focus more directly on disappointing outcomes: variability among losses, shortfalls relative to a target, or the severity of a distribution's loss tail. These are different questions, so their definitions and conventions matter.

Semivolatility: Focusing on Negative Fluctuations

In this toolkit, semivolatility means the standard deviation of the negative returns around their own mean. Let $n_{-}$ be the number of negative observations and ${\overline{r}}_{-}$ their arithmetic mean:

$${\hat{\sigma }}_{-}=\sqrt{\frac{1}{n_{-}}\sum _{r_t<0}(r_t-{\overline{r}}_{-}{)}^2}.$$

This is conditional dispersion among the observed losses, using ddof=0. It does not measure their distance from zero or account directly for how often losses occur. A sequence of identical negative returns has zero semivolatility even though every observation is a loss.

Downside Deviation Relative to a Target

Another convention, often used in the Sortino ratio, is target downside deviation. For a per-period target return $\tau$, use:

$${\widehat{\mathrm{DD}}}_{\tau }=\sqrt{\frac{1}{n}\sum _{t=1}^n{[min(r_t-\tau ,0)]}^2}.$$

Returns above the target contribute zero; returns below it contribute their squared shortfall. The denominator is the number of all observed returns, so both the frequency and size of shortfalls matter. The target must have the same frequency and units as the returns.

For monthly returns of −4%, −2%, +1%, and +3%, the negative subset has mean −3% and semivolatility 1%. Relative to a zero target, downside deviation is $\sqrt{({0.04}^2+{0.02}^2)/4}\approx 2.24\mathrm{\%}$.

import numpy as np
import pandas as pd
import PortfolioOptimizationKit as pok

monthly_returns = pd.Series([-0.04, -0.02, 0.01, 0.03])
target = 0.0
negative_volatility = pok.semivolatility(monthly_returns)
shortfalls = (monthly_returns - target).clip(upper=0)
downside_deviation = np.sqrt((shortfalls ** 2).mean())

print(f"Negative-return semivolatility: {negative_volatility:.2%}")
print(f"Monthly target: {target:.2%}")
print(f"Target downside deviation: {downside_deviation:.2%}")

With no negative observations, conditional semivolatility is unavailable and the toolkit returns NaN. For a nonempty sample, downside deviation is zero if all observed returns meet the target. Neither result guarantees the absence of future losses.

Try it: change target to 0.01. Downside deviation increases to approximately 2.92%, because the investor now requires a 1% monthly return. The negative-return semivolatility remains 1% because its definition has not changed.

Value at Risk (VaR): A Loss Quantile

Let $R$ be the return over a specified horizon and $L=-R$ the corresponding loss as a fraction of starting portfolio value. At confidence level $0<\alpha <1$, Value at Risk is the $\alpha$-quantile of the loss distribution:

$${\mathrm{VaR}}_{\alpha }(L)=q_{\alpha }(L)=inf\{\ell :F_L(\ell )\ge \alpha \}.$$

When the return CDF is continuous and strictly increasing, write $\beta =1-\alpha$ for the lower-tail probability. Then:

$${\mathrm{VaR}}_{\alpha }(L)=-q_{\beta }(R).$$

A 95% monthly VaR of 3% identifies a loss threshold with 5% probability above it under the assumed continuous model. It is neither a maximum possible loss nor an average loss. Outcomes beyond that threshold can be substantially worse.

The toolkit takes the tail probability as level: use level=0.05 for 95% confidence. Return observations determine the horizon; monthly inputs produce monthly VaR. The reported result is a loss fraction, so 0.03 corresponds to 3% of starting portfolio value. A negative result is possible if the estimated loss quantile still represents a gain; the functions retain that sign rather than taking an absolute value.

Worked Example: Historical VaR with Linear Interpolation

Use the same ten monthly returns from the quantile example. For 90% confidence, $\beta =0.10$. With NumPy's linear interpolation convention, the estimated return quantile is −4.30%, giving a historical VaR estimate of 4.30%.

import numpy as np
import pandas as pd
import PortfolioOptimizationKit as pok

tail_returns = pd.Series([-0.04, 0.05, 0.02, -0.07, 0.01, 0.005, -0.02, -0.01, -0.02, 0.05])
confidence = 0.90
tail_probability = 1 - confidence

var_from_returns = -np.quantile(tail_returns, tail_probability, method="linear")
var_from_losses = np.quantile(-tail_returns, confidence, method="linear")
historical_var = pok.var_historic(tail_returns, level=tail_probability)
portfolio_value = 100_000.0

print(f"VaR from return quantile: {var_from_returns:.2%}")
print(f"VaR from loss quantile: {var_from_losses:.2%}")
print(f"Toolkit historical VaR: {historical_var:.2%}")
print(f"VaR for a $100,000 portfolio: ${portfolio_value * historical_var:,.2f}")

The three estimates agree, giving a monetary VaR of $4,300. In this finite-sample calculation, the return and loss forms agree because they use linear interpolation. Inverse empirical CDF conventions can differ at probability jumps, so the interpolation rule must be stated. The estimate summarizes this historical sample; it does not establish a known 10% future exceedance probability.

Exploring Conditional Value at Risk (CVaR)

VaR locates a tail boundary. Expected shortfall (ES), called CVaR here, measures the average loss in the worst $\beta =1-\alpha$ fraction of outcomes. A definition that also handles distributions with probability masses is:

$${\mathrm{ES}}_{\alpha }(L)=\frac{1}{1-\alpha }{\int }_{\alpha }^1{q}_u(L)du.$$

For a continuous distribution, this is the conditional mean loss beyond VaR:

$${\mathrm{ES}}_{\alpha }(L)=\mathbb{E}[L\mid L\ge {\mathrm{VaR}}_{\alpha }(L)]=-\mathbb{E}[R\mid R\le q_{\beta }(R)].$$

For discrete observations, simply selecting returns strictly below a cutoff can drop tied values or produce an empty tail. We instead assign probability $1/n$ to each observed return and average exactly the worst $\beta$ fraction, including a fractional observation at the boundary when necessary.

Illustrative Example

Order the $n$ returns from worst to best as $r_{(1)}\le \cdots \le r_{(n)}$. Let $m=n\beta$, $k=\lfloor m\rfloor$, and $f=m-k$. The empirical expected shortfall is:

$${\widehat{\mathrm{ES}}}_{\alpha }=-\frac{{\displaystyle \sum _{i=1}^k{r}_{(i)}+f{r}_{(k+1)}}}{m}.$$

At 80% confidence in our ten-observation sample, $m=2$: average the two worst returns, −7% and −4%, giving 5.50% ES. At 75% confidence, $m=2.5$: include half of the next-worst return, −2%, giving:

$${\widehat{\mathrm{ES}}}_{0.75}=-\frac{-0.07-0.04+0.5(-0.02)}{2.5}=4.80\mathrm{\%}.$$

import numpy as np
import pandas as pd
import PortfolioOptimizationKit as pok

tail_returns = pd.Series([-0.04, 0.05, 0.02, -0.07, 0.01, 0.005, -0.02, -0.01, -0.02, 0.05])
confidence = 0.75
tail_probability = 1 - confidence
ordered_returns = tail_returns.sort_values().to_numpy()
tail_mass = len(ordered_returns) * tail_probability
weights = np.clip(tail_mass - np.arange(len(ordered_returns)), 0, 1)
empirical_es = -np.dot(weights, ordered_returns) / tail_mass
toolkit_es = pok.cvar_historic(tail_returns, level=tail_probability)

print(f"Tail mass in observations: {tail_mass:.2f}")
print(f"Empirical expected shortfall: {empirical_es:.2%}")
print(f"Toolkit expected shortfall: {toolkit_es:.2%}")

The weights are $1$, $1$, $0.5$, then zeros. This implements the inverse-empirical-CDF tail average, while var_historic() uses linear interpolation for its threshold estimate. Those finite-sample conventions are explicit rather than silently treating every observation below an interpolated cutoff as an equally complete tail observation.

Try it: change confidence to 0.80, then 0.90. ES becomes 5.50% and 7.00%, respectively. At 90% confidence, the worst 10% of ten observations contains only the worst observation. With identical losses, ES equals that loss rather than becoming undefined.

Estimating VaR and CVaR: A Comparative Overview

We now compare estimates using the same monthly hedge-fund returns. Historical methods use observed outcomes; Gaussian methods impose a distributional model; Cornish–Fisher adjusts a Gaussian quantile using estimated skewness and kurtosis. Expected shortfall answers a different question from VaR: tail-average severity rather than a threshold.

Historical Approach (Non-Parametric)

Start with CTA Global. The histogram displays monthly returns in percent. Its vertical axis is density per percentage point, so density times bin width sums to one.

import json
import numpy as np
import PortfolioOptimizationKit as pok

hfi = pok.get_hfi_returns()
cta_returns = hfi["CTA Global"].dropna()
density, bin_edges = np.histogram(cta_returns * 100, bins=30, density=True)
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2

print(f"Observations: {len(cta_returns)}")
print(f"Mean monthly return: {cta_returns.mean():.2%}")
print(f"Monthly volatility (ddof=0): {cta_returns.std(ddof=0):.2%}")
print(f"Skewness: {pok.skewness(cta_returns):.3f}")
print(f"Pearson kurtosis: {pok.kurtosis(cta_returns):.3f}")

plot_data = {
    "ctaDistribution": {
        "title": "CTA Global Monthly Returns",
        "type": "bar",
        "xAxis": {"type": "value", "name": "Monthly return (%)"},
        "yAxisName": "Density per percentage point",
        "series": {"Density": np.column_stack([bin_centers, density]).tolist()},
    },
}
print("\n<ECHARTS_DATA>" + json.dumps(plot_data, allow_nan=False))

Calculate monthly VaR and ES at 90%, 95%, and 99% confidence. The corresponding tail probabilities are 10%, 5%, and 1%.

import pandas as pd

# Run the CTA Global histogram example first.
tail_estimates = []
for confidence in [0.90, 0.95, 0.99]:
    tail_probability = 1 - confidence
    tail_estimates.append({
        "Confidence": f"{confidence:.0%}",
        "Historical VaR (%)": pok.var_historic(cta_returns, level=tail_probability) * 100,
        "Historical ES (%)": pok.cvar_historic(cta_returns, level=tail_probability) * 100,
        "Tail mass (observations)": len(cta_returns) * tail_probability,
    })
tail_comparison = pd.DataFrame(tail_estimates).set_index("Confidence")
print(tail_comparison.to_string(float_format=lambda value: f"{value:.2f}"))

The historical VaRs are approximately 2.41%, 3.17%, and 4.95%; the empirical ES estimates are 3.50%, 4.19%, and 5.50%. These describe increasingly severe parts of the observed monthly loss distribution.

At 99% confidence, the 263-observation sample supplies only 2.63 observations' worth of tail probability. Raising confidence does not create additional tail information. Historical estimates cannot include losses absent from the sample, and these monthly estimates should not be mechanically interpreted as daily or annual measures.

Parametric Approach (Gaussian)

Assume $R\sim \mathcal{N}(\mu ,{\sigma }^2)$ over the chosen horizon. From the quantile derivation earlier, the lower-tail return threshold is $q_{\beta }(R)=\mu +\sigma {\mathrm{\Phi }}^{-1}(\beta )$. Negating it gives:

$${\mathrm{VaR}}_{\alpha }^{\mathrm{Gaussian}}=-[\mu +\sigma {\mathrm{\Phi }}^{-1}(1-\alpha )].$$

For a 1% monthly mean, 4% monthly volatility, and 95% confidence, this is $-[0.01+0.04(-1.6449)]\approx 5.58\mathrm{\%}$. norm.ppf() supplies the standardized quantile; the mean and volatility are separate inputs.

When fitting observed returns, the toolkit uses their arithmetic mean and std(ddof=0), the Gaussian maximum-likelihood estimates under iid normal sampling. This scale convention differs from the ddof=1 sample-volatility report earlier. Constant observed returns are treated as a zero-variance point mass for this fitted quantile, which does not establish certainty about future returns.

from scipy.stats import norm
import PortfolioOptimizationKit as pok

confidence = 0.95
tail_probability = 1 - confidence
mu, sigma = 0.01, 0.04
model_var = -(mu + sigma * norm.ppf(tail_probability))
print(f"Illustrative model VaR: {model_var:.2%}")

hfi = pok.get_hfi_returns()
gaussian_var = pok.var_gaussian(hfi, level=tail_probability)
print(f"\n{confidence:.0%} monthly Gaussian VaR for the indices (%):")
print((gaussian_var * 100).to_string(float_format=lambda value: f"{value:.2f}"))

Try it: raise confidence to 0.99. The illustrative model's VaR increases to approximately 8.31%. Higher confidence places the threshold farther into the loss tail under the same Gaussian model.

Cornish-Fisher Modification (Semi-Parametric)

The Cornish–Fisher expansion adjusts a standard-normal quantile using estimated skewness $S$ and Pearson kurtosis $K$. Set $\beta =1-\alpha$ and $z_{\beta }={\mathrm{\Phi }}^{-1}(\beta )$. The truncated expansion used here is:

$${\tilde{z}}_{\beta }=z_{\beta }+\frac{z_{\beta }^2-1}{6}S+\frac{z_{\beta }^3-3z_{\beta }}{24}(K-3)-\frac{2z_{\beta }^3-5z_{\beta }}{36}{S}^2.$$

The corresponding modified VaR estimate is:

$${\widehat{\mathrm{VaR}}}_{\alpha }^{\mathrm{CF}}=-(\hat{\mu }+\hat{\sigma }{\tilde{z}}_{\beta }).$$

When $S=0$ and $K=3$, the correction terms vanish. Otherwise, they approximate a change in the lower-tail quantile. This is not a guarantee of a better estimate: large or unstable higher moments can make the approximation unreliable, and the resulting quantiles can even violate the expected ordering across confidence levels.

from scipy.stats import norm
import PortfolioOptimizationKit as pok

index_name = "CTA Global"
returns = pok.get_hfi_returns()[index_name].dropna()
confidence = 0.95
tail_probability = 1 - confidence
z = norm.ppf(tail_probability)
S, K = pok.skewness(returns), pok.kurtosis(returns)
adjusted_z = (
    z + (z ** 2 - 1) * S / 6
    + (z ** 3 - 3 * z) * (K - 3) / 24
    - (2 * z ** 3 - 5 * z) * S ** 2 / 36
)
modified_var = -(returns.mean() + returns.std(ddof=0) * adjusted_z)
toolkit_modified_var = pok.var_gaussian(returns, level=tail_probability, cf=True)

print(f"{index_name}: skewness {S:.3f}, Pearson kurtosis {K:.3f}")
print(f"Gaussian quantile: {z:.4f} | Adjusted quantile: {adjusted_z:.4f}")
print(f"Cornish–Fisher VaR: {modified_var:.2%} | Toolkit: {toolkit_modified_var:.2%}")

For CTA Global, the modified 95% monthly VaR is approximately 3.31%, versus 3.42% under the fitted Gaussian model. Try it: change index_name to "Convertible Arbitrage" and examine the much larger estimated skewness and kurtosis. Consider whether a truncated adjustment remains a convincing tail model for that sample.

Comparison of VaR Methods

Use the same 95% confidence level and monthly observations for all four columns below. The first three columns estimate a VaR threshold; the fourth reports historical expected shortfall. Every value is expressed as a percentage of starting portfolio value.

import json
import pandas as pd
import PortfolioOptimizationKit as pok

hfi = pok.get_hfi_returns()
confidence = 0.95
tail_probability = 1 - confidence
risk_comparison = pd.DataFrame({
    "Historical VaR": pok.var_historic(hfi, level=tail_probability),
    "Gaussian VaR": pok.var_gaussian(hfi, level=tail_probability),
    "Cornish-Fisher VaR": pok.var_gaussian(hfi, level=tail_probability, cf=True),
    "Historical ES": pok.cvar_historic(hfi, level=tail_probability),
})
comparevars = risk_comparison * 100
print(f"{confidence:.0%} monthly loss measures (%):")
print(comparevars.to_string(float_format=lambda value: f"{value:.2f}"))

plot_data = {
    "varComparison": {
        "type": "bar",
        "title": f"{confidence:.0%} Monthly Loss Measures",
        "yAxisName": "Loss (% of portfolio value)",
        "series": {name: comparevars[name].tolist() for name in comparevars},
        "xAxis": {
            "type": "category",
            "name": "Strategy",
            "data": comparevars.index.tolist(),
            "axisLabel": {"rotate": 45, "interval": 0},
        },
    },
}
print("\n<ECHARTS_DATA>" + json.dumps(plot_data, allow_nan=False))

Within a single distribution, expected shortfall averages the upper loss tail and is at least as large as its VaR threshold. Comparisons across different fitted models and finite-sample estimators require more care: there is no universal ordering between historical ES and Gaussian or Cornish–Fisher VaR. Inspect the assumptions and tail observations rather than treating the largest number as automatically the most accurate forecast.

Chapter Summary

• Compounded growth measures the evolution of wealth; the arithmetic mean answers a different question about period returns.
• Volatility and Sharpe ratios require consistent frequencies and clearly defined return and benchmark conventions.
• Drawdown measures losses relative to a running peak and depends on the sequence of returns.
• Skewness, kurtosis, and normality tests diagnose distributional assumptions without proving a model correct.
• VaR identifies a loss threshold; expected shortfall measures average severity in the worst tail. Specify the horizon, confidence level, and estimation method for both.

These measurements provide the inputs and diagnostics for the next chapter, Portfolio Optimization, where we combine assets and study the trade-off between expected return and portfolio risk.