Week 9: Simulation-Based Risk Measurement

Version 5.0 - August 2026

Author
Affiliation

Jon Danielsson

London School of Economics

9 Week 9

An option is nonlinear. When the underlying price moves by one dollar, the option price does not always move by the same amount. That makes a simple weighted-return calculation unsuitable for a portfolio containing options.

Simulation provides a direct solution. Generate possible future stock prices, reprice the option in every scenario, and inspect the resulting profit-and-loss distribution. We first validate the simulation against the Black-Scholes price, then use it to calculate risk for a stock-and-option portfolio.

NoteBefore you start

Open the main course folder and check that Prices.csv and Returns.csv, created in Week 2, are present. You will use random simulation from Week 5 and VaR and ES conventions from Week 8. This week introduces no additional libraries.

TipBy the end of this session

You should be able to simulate a stock price, validate a simulated option price against Black-Scholes, reprice an option in every risk scenario, and explain why an option makes portfolio risk nonlinear.

9.1 The plan for this week

  1. Price one option with a supplied Black-Scholes function
  2. Simulate the same option price and compare the answers
  3. Simulate a one-day stock-and-option profit and loss
  4. Calculate VaR and ES from that simulated distribution

9.2 Loading data and libraries

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats

Prices = pd.read_csv("Prices.csv", index_col="date", parse_dates=True)
Returns = pd.read_csv("Returns.csv", index_col="date", parse_dates=True)
rng = np.random.default_rng(888)
NoteCheck your result

Both files should contain a JPM column and have dates as their index. Prices should be positive, and returns are daily decimals.

assert "JPM" in Prices.columns and "JPM" in Returns.columns
assert (Prices["JPM"].dropna() > 0).all()
print(Prices.index.min(), "to", Prices.index.max())
2000-01-03 00:00:00 to 2025-09-23 00:00:00
CautionUnits and signs used this week
  • Stock prices and portfolio values are in dollars.
  • Historical returns are daily decimals.
  • Black-Scholes volatility and time are annual: 0.20 means 20% per year and 0.25 means one quarter of a year.
  • VaR and ES are positive loss amounts.
  • A one-day risk scenario reduces option maturity by 1 / 252.

9.3 A function for the Black-Scholes call price

Run this function unchanged. Its job is to turn stock price, strike, interest rate, volatility and time to maturity into a European call price.

def black_scholes_call(stock_price, strike, rate, volatility, maturity):
    d1 = (
        np.log(stock_price / strike)
        + (rate + 0.5 * volatility**2) * maturity
    ) / (volatility * np.sqrt(maturity))
    d2 = d1 - volatility * np.sqrt(maturity)
    return (
        stock_price * stats.norm.cdf(d1)
        - strike * np.exp(-rate * maturity) * stats.norm.cdf(d2)
    )
Financial quantity Python name Unit
Current stock price \(S_0\) stock_price Dollars
Strike \(K\) strike Dollars
Risk-free rate \(r\) rate Annual decimal
Volatility \(\sigma\) volatility Annual decimal
Time to maturity \(T\) maturity Years
\(d_1,d_2\) d1, d2 No units

9.4 Price one option

We use the latest JPM price, an at-the-money strike, an annualised historical volatility estimate and three months to maturity.

This is an illustrative valuation rather than an attempt to match a traded option price. The formula above is for a stock that pays no dividend, and JPM pays one. The volatility is estimated from past returns, whereas a market price reflects the volatility other participants expect over the life of the option. Both simplifications are usual in a first treatment and both would have to go in a valuation meant for trading.

S0 = Prices["JPM"].dropna().iloc[-1]
K = S0
r = 0.03
sigma_daily = Returns["JPM"].dropna().tail(1000).std()
sigma_annual = sigma_daily * np.sqrt(252)
T = 0.25

analytic_price = black_scholes_call(S0, K, r, sigma_annual, T)
print("Stock price: $", round(S0, 2))
print("Annual volatility:", round(sigma_annual, 3))
print("Black-Scholes call price: $", round(analytic_price, 2))
Stock price: $ 312.74
Annual volatility: 0.248
Black-Scholes call price: $ 16.61
NoteCheck your result

The call price should be positive and below the stock price. Check that you did not put daily volatility into a formula expecting annual volatility.

assert 0 < analytic_price < S0

9.5 Simulate the option price

Under the Black-Scholes assumptions, a risk-neutral future stock price can be simulated as

\[ S_T=S_0\exp\left((r-\tfrac12\sigma^2)T+\sigma\sqrt{T}Z\right), \qquad Z\sim N(0,1). \]

simulations = 100_000
z_option = rng.standard_normal(simulations)
simulated_ST = S0 * np.exp(
    (r - 0.5 * sigma_annual**2) * T
    + sigma_annual * np.sqrt(T) * z_option
)
call_payoff = np.maximum(simulated_ST - K, 0)
discounted_payoff = np.exp(-r * T) * call_payoff
simulation_price = discounted_payoff.mean()
simulation_se = discounted_payoff.std(ddof=1) / np.sqrt(simulations)

print("Black-Scholes price: $", round(analytic_price, 3))
print("Simulation price:    $", round(simulation_price, 3))
print("Simulation SE:       $", round(simulation_se, 3))
Black-Scholes price: $ 16.611
Simulation price:    $ 16.744
Simulation SE:       $ 0.08

The simulation is an approximation, so the two prices will not be identical. The standard error describes ordinary Monte Carlo sampling uncertainty.

NoteCheck your result

The simulated price should be close to the Black-Scholes price. With this fixed seed and sample size, their difference should be within four simulation standard errors.

assert abs(simulation_price - analytic_price) < 4 * simulation_se

9.6 One-day option and portfolio risk

Pricing used the risk-free drift because it valued the option. Risk measurement now asks about possible one-day market moves. Over one day the drift is set so that the expected future price equals today’s price, using the estimated daily volatility. The \(-\tfrac12\sigma^2\) term in the exponent is what achieves that: it leaves the expected simple return at zero, so the expected log return is slightly negative.

p = 0.01
risk_scenarios = 100_000
z_risk = rng.standard_normal(risk_scenarios)
S1 = S0 * np.exp(-0.5 * sigma_daily**2 + sigma_daily * z_risk)
T1 = T - 1 / 252
call_price_tomorrow = black_scholes_call(S1, K, r, sigma_annual, T1)

Suppose the portfolio contains ten JPM shares and ten calls, each written on one share. A listed equity option contract normally covers 100 shares, with the premium quoted per share, so a real position of ten contracts would be a hundred times this one. Repricing the option in every scenario captures its nonlinearity.

number_of_shares = 10
number_of_calls = 10

stock_pnl = number_of_shares * (S1 - S0)
option_pnl = number_of_calls * (call_price_tomorrow - analytic_price)
portfolio_pnl = stock_pnl + option_pnl

pnl_threshold = np.quantile(portfolio_pnl, p, method="lower")
portfolio_var = -pnl_threshold
portfolio_es = -portfolio_pnl[portfolio_pnl <= pnl_threshold].mean()

print("Stock-and-option VaR: $", round(portfolio_var, 2))
print("Stock-and-option ES:  $", round(portfolio_es, 2))
Stock-and-option VaR: $ 168.88
Stock-and-option ES:  $ 191.85
NoteCheck your result

VaR and ES should be positive, and ES should be at least as large as VaR.

assert portfolio_var > 0
assert portfolio_es >= portfolio_var

To see that the option is not a fixed multiple of the stock, replace the repricing with a single sensitivity. Delta is the derivative of the call price with respect to the stock price, \(N(d_1)\), and a linear approximation applies it to every scenario:

d1 = (
    np.log(S0 / K) + (r + 0.5 * sigma_annual**2) * T
) / (sigma_annual * np.sqrt(T))
delta = stats.norm.cdf(d1)

linear_pnl = stock_pnl + number_of_calls * delta * (S1 - S0)
linear_threshold = np.quantile(linear_pnl, p, method="lower")
linear_var = -linear_threshold

print("Full repricing VaR:      $", round(portfolio_var, 2))
print("Delta approximation VaR: $", round(linear_var, 2))
Full repricing VaR:      $ 168.88
Delta approximation VaR: $ 174.12

The two differ because delta itself changes as the stock price moves. A long call gains delta as the price rises and loses it as the price falls, so on the downside the option loses less than a fixed delta predicts, and the linear approximation overstates the loss. That curvature is what “nonlinear” means here, and it is the reason the option is repriced in every scenario rather than represented by one number. How large the gap is depends on the strike and on the horizon. A deep in-the-money call behaves almost like the stock over one day, and the two numbers then nearly agree.

fig, ax = plt.subplots()
ax.hist(portfolio_pnl, bins=60, edgecolor="white")
ax.axvline(pnl_threshold, color="red", label="1% P&L threshold")
ax.set_title("Simulated one-day stock-and-option P&L")
ax.set_xlabel("Profit and loss ($)")
ax.set_ylabel("Count")
ax.legend(frameon=False)
plt.show()

The nonlinear step is the repricing: black_scholes_call is applied to every possible S1. Replacing that with one constant option sensitivity would be a linear approximation and could miss curvature in larger market moves.

9.7 An equivalent lognormal draw

numpy can draw the same future-price distribution directly. This is an alternative simulation route, not another model.

direct_rng = np.random.default_rng(888)
direct_ST = S0 * direct_rng.lognormal(
    mean=(r - 0.5 * sigma_annual**2) * T,
    sigma=sigma_annual * np.sqrt(T),
    size=simulations,
)
direct_price = np.exp(-r * T) * np.maximum(direct_ST - K, 0).mean()
print("Direct-lognormal price: $", round(direct_price, 3))
Direct-lognormal price: $ 16.744

9.8 Simulation size

A single estimate at each sample size would not show anything. A small sample can land closer to the analytic price than a large one by luck. What changes with sample size is the spread of the estimates, so each size below is run twenty times and the standard deviation of its twenty prices reported alongside their average.

convergence_rng = np.random.default_rng(321)
repetitions = 20
convergence = {}
for sample_size in [1_000, 10_000, 100_000]:
    estimates = []
    for repetition in range(repetitions):
        z = convergence_rng.standard_normal(sample_size)
        future_price = S0 * np.exp(
            (r - 0.5 * sigma_annual**2) * T
            + sigma_annual * np.sqrt(T) * z
        )
        payoff = np.maximum(future_price - K, 0)
        estimates.append(np.exp(-r * T) * payoff.mean())
    convergence[sample_size] = {
        "Mean price": np.mean(estimates),
        "Standard deviation": np.std(estimates, ddof=1),
    }

pd.DataFrame(convergence).T.rename_axis("Draws").round(4)
Mean price Standard deviation
Draws
1000 16.1834 0.8509
10000 16.6574 0.1776
100000 16.6169 0.0627

The average sits near the analytic price at every size. The standard deviation falls by roughly a factor of ten for each hundredfold increase in draws, which is the \(1/\sqrt{n}\) rate of Monte Carlo error.

9.9 Two correlated assets

This extension uses the historical covariance matrix to generate two correlated one-day returns. It is useful when a portfolio contains more than one source of market risk, but it is not needed for the core option exercise.

The draws are log returns, and the P&L below multiplies them by the position values as though they were simple returns. Over one day the two are close enough for the purpose. The exact calculation would use np.exp(simulated_returns) - 1.

two_asset_returns = Returns[["JPM", "MCD"]].dropna().tail(1000)
covariance = two_asset_returns.cov().to_numpy()
multi_rng = np.random.default_rng(456)
simulated_returns = multi_rng.multivariate_normal(
    mean=np.zeros(2),
    cov=covariance,
    size=100_000,
)

two_asset_values = np.array([1000, 1000])
two_asset_pnl = simulated_returns @ two_asset_values
two_asset_threshold = np.quantile(two_asset_pnl, p, method="lower")
two_asset_var = -two_asset_threshold
print("Two-asset VaR: $", round(two_asset_var, 2))
Two-asset VaR: $ 50.39

9.10 Recap

  • Black-Scholes provides an analytic benchmark for a European call.
  • A risk-neutral simulation should reproduce that benchmark within simulation uncertainty.
  • Risk scenarios and pricing scenarios answer different questions and need not use the same drift.
  • An option makes portfolio risk nonlinear because it must be repriced as the underlying price changes.
  • Simulation-based VaR is a lower P&L quantile reported as a positive loss.
  • Larger simulated samples reduce, but do not eliminate, Monte Carlo noise.

9.11 Exercises

9.11.1 Apply the method

  1. Change the strike to 90% and 110% of the current stock price. Compare the Black-Scholes and simulated prices.
  2. Change the portfolio to ten shares and five calls. Recalculate VaR and ES.
  3. Explain why simply adding the standalone stock VaR and option VaR need not equal portfolio VaR.

9.11.2 Compare results

  1. Repeat the pricing simulation with three seeds. Compare the price difference with its reported simulation standard error.
  2. Compare sample sizes of 1,000, 10,000 and 100,000 across several repeats.

9.11.3 Extend the analysis

  1. Simulate a stock-and-put portfolio by supplying a Black-Scholes put-pricing function and repricing the put in every scenario.
  2. Extend the two-asset simulation to reprice an option on one asset while holding shares in the other.