{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Calibrating the Heston Model" ], "id": "17a74a4a-cbcb-4758-be54-64f7d06af8c3" }, { "cell_type": "raw", "metadata": { "raw_mimetype": "tex" }, "source": [ "\\renewcommand{\\perp}{\\mathrel\\bot}\n", "\\newcommand{\\ev}{\\operatorname{E}}\n", "\\newcommand{\\var}{\\operatorname{V}}\n", "\\newcommand{\\cov}{\\operatorname{Cov}}\n", "\\newcommand{\\Normal}{\\mathcal{N}}\n", "\\newcommand{\\qedf}{\\qed{\\parfillskip=0pt \\par}}\n", "\\renewcommand{\\vec}[1]{\\mathbf{#1}}\n", "\\newcommand{\\gvec}[1]{\\pmb{#1}}\n", "\\newcommand{\\inner}[1]{\\langle{#1}\\rangle}\n", "\\newcommand{\\norm}[1]{\\lVert{#1}\\rVert}\n", "\\newcommand{\\prob}{\\operatorname{P}}\n", "\\newcommand{\\1}[1]{1\\kern-0.25em\\text{l}_{\\{#1\\}}}" ], "id": "6705799c-cba7-4fca-bdaa-2f5c117bd788" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction\n", "\n", "The [Heston Model](../derivatives-and-bond-pricing/heston-model.qmd)\n", "notebook derives the closed-form option pricing formula from the\n", "stochastic volatility dynamics and the characteristic function of the\n", "log stock price. This notebook focuses entirely on implementation and\n", "calibration: we translate that formula into Python, and then find the\n", "model parameters that best match observed call prices on Apple (AAPL) as\n", "of October 1, 2024.\n", "\n", "Calibration means solving an inverse problem. Given a set of observed\n", "option prices $C^{\\text{obs}}(K_i, T_j)$ across strikes\n", "$K_1, \\ldots, K_m$ and maturities $T_1, \\ldots, T_n$, we set $v_0$ from\n", "the short-dated ATM implied volatility and search for the parameter\n", "vector $(\\kappa, \\theta, \\sigma, \\rho)$ that minimizes the root\n", "mean-squared pricing error $$\n", " \\text{RMSE} = \\sqrt{\\frac{1}{mn} \\sum_{i=1}^{m} \\sum_{j=1}^{n} \\left(C^{\\text{obs}}(K_i, T_j) - C^{\\text{model}}(K_i, T_j)\\right)^2}.\n", "$$ Once calibrated, the model can be used to price exotic derivatives,\n", "compute sensitivities, or simulate risk scenarios with a volatility\n", "structure that is consistent with the market.\n", "\n", "This notebook presents the simplest possible calibration exercise: we\n", "fit all parameters to a single cross-section of option prices on one\n", "date. A more rigorous approach would split the data into an estimation\n", "sample and a hold-out sample, calibrate $(\\kappa, \\theta, \\sigma, \\rho)$\n", "on the training set, and then backtest the model by repricing the\n", "hold-out contracts with those fixed parameters. Using the short-dated\n", "ATM implied volatility as a proxy for $\\sqrt{v_0}$ is a natural choice\n", "in any such exercise, because $v_0$ is the only parameter that should be\n", "updated daily as market conditions change, while the structural\n", "parameters $(\\kappa, \\theta, \\sigma, \\rho)$ are more stable and can be\n", "estimated over a longer horizon.\n", "\n", "One-factor stochastic volatility models like Heston face inherent\n", "limitations in fitting the full surface simultaneously across strikes\n", "and maturities. Multifactor extensions that allow the variance process\n", "to have richer dynamics can substantially improve the fit, as shown in\n", "Cortazar et al. (2017) in the context of commodity markets.\n", "\n", "## Implementing the Pricing Formula\n", "\n", "The Heston call price is $$\n", " C = S e^{-q T} P_1 - K e^{-r T} P_2,\n", "$$ where $P_2$ is the risk-neutral probability of expiring in the money\n", "and $P_1$ is the corresponding probability under the stock-numeraire\n", "measure. Both are recovered by Fourier inversion of the joint\n", "characteristic function of $x(T) = \\ln S(T)$ and $v(T)$.\n", "\n", "We begin by loading the required libraries." ], "id": "1f91803c-816d-4ad7-97c3-d5a374ff0801" }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from numpy import log, exp, pi, real\n", "from scipy.integrate import quad\n", "from scipy.stats import norm\n", "from scipy.optimize import brentq, minimize\n", "import matplotlib.pyplot as plt\n", "import pandas as pd" ], "id": "d7823052" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The characteristic function\n", "\n", "The joint characteristic function of $x(T)$ and $v(T)$ is $$\n", "\\begin{gathered}\n", " f(\\phi, \\varphi) = \\exp\\Bigg(i (r - q) \\phi T + \\left(\\delta - \\frac{2 \\gamma}{\\sigma^{2}}\\right) \\kappa \\theta^{*} T + i \\phi x + \\delta v \\\\\n", " \\qquad - \\frac{2 \\kappa \\theta^{*}}{\\sigma^{2}} \\ln\\!\\left(e^{-\\gamma T} - \\frac{\\sigma^{2}}{2 \\gamma} (1 - e^{-\\gamma T}) (i \\varphi - \\delta)\\right) \\\\\n", " \\qquad + \\frac{v (i \\varphi - \\delta)}{e^{- \\gamma T} - \\frac{\\sigma^{2}}{2 \\gamma} (1 - e^{- \\gamma T}) (i \\varphi - \\delta)}\\Bigg),\n", "\\end{gathered}\n", "$$ where $$\n", "\\begin{aligned}\n", " \\gamma &= \\sqrt{\\kappa^{2} + (1 - \\rho^{2}) \\sigma^{2} \\phi^{2} + i (\\sigma - 2 \\kappa \\rho) \\sigma \\phi}, \\\\\n", " \\delta &= \\frac{\\kappa + \\gamma - i \\rho \\sigma \\phi}{\\sigma^{2}}.\n", "\\end{aligned}\n", "$$ The implementation puts every factor inside the exponential to avoid\n", "numerical instabilities from branch cuts in the complex square root." ], "id": "12fc2594-a34e-4b2a-90c5-225e8c6478f4" }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def CF_Heston(phi1, phi2, S0, V0, T, r, q, kappa, theta, sigma, rho):\n", " gamma = (kappa**2 + (1 - rho**2) * sigma**2 * phi1**2\n", " + 1j * (sigma - 2 * kappa * rho) * sigma * phi1)**0.5\n", " delta = (kappa + gamma - 1j * rho * sigma * phi1) / sigma**2\n", "\n", " y = exp(\n", " 1j * (r - q) * phi1 * T\n", " + (delta - 2 * gamma / sigma**2) * kappa * theta * T\n", " + 1j * phi1 * log(S0)\n", " + delta * V0\n", " + (-2 * kappa * theta / sigma**2)\n", " * np.log(exp(-gamma * T) - sigma**2 / (2 * gamma)\n", " * (1 - exp(-gamma * T)) * (1j * phi2 - delta))\n", " + V0 * (1j * phi2 - delta)\n", " / (exp(-gamma * T) - sigma**2 / (2 * gamma)\n", " * (1 - exp(-gamma * T)) * (1j * phi2 - delta))\n", " )\n", " return y" ], "id": "6f5b73d4" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### In-the-money probabilities\n", "\n", "The risk-neutral probability $P_2$ that the call expires in the money is\n", "recovered by Fourier inversion of $f(\\phi, 0)$." ], "id": "34b6ab77-b1ae-47be-8b87-143f1a180dd9" }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def P2(S0, V0, K, T, r, q, kappa, theta, sigma, rho):\n", " integrand = lambda phi1: real(\n", " exp(-1j * phi1 * log(K))\n", " * CF_Heston(phi1, 0, S0, V0, T, r, q, kappa, theta, sigma, rho)\n", " / (1j * phi1)\n", " )\n", " y = 0.5 + 1 / pi * quad(integrand, 0, np.inf, full_output=1)[0]\n", " return y" ], "id": "b2a1a241" }, { "cell_type": "markdown", "metadata": {}, "source": [ "The probability $P_1$ under the stock-numeraire measure is obtained by\n", "evaluating the characteristic function at $\\phi - i$ instead of $\\phi$,\n", "then dividing by the forward price $F = S e^{(r-q)T}$." ], "id": "6397ae18-b451-443b-a37f-4778bfcecc56" }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def P1(S0, V0, K, T, r, q, kappa, theta, sigma, rho):\n", " F = S0 * exp((r - q) * T)\n", " integrand = lambda phi1: real(\n", " exp(-1j * phi1 * log(K))\n", " * CF_Heston(phi1 - 1j, 0, S0, V0, T, r, q, kappa, theta, sigma, rho)\n", " / (1j * phi1)\n", " )\n", " y = 0.5 + 1 / pi / F * quad(integrand, 0, np.inf, full_output=1)[0]\n", " return y" ], "id": "230868f1" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Call price\n", "\n", "The European call price follows directly from the two probabilities." ], "id": "58d5c0b2-7f37-480f-90fa-1f74ce750fe8" }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def heston_call_price(S0, V0, K, T, r, q, kappa, theta, sigma, rho):\n", " call = (S0 * exp(-q * T) * P1(S0, V0, K, T, r, q, kappa, theta, sigma, rho)\n", " - K * exp(-r * T) * P2(S0, V0, K, T, r, q, kappa, theta, sigma, rho))\n", " return call" ], "id": "77177b58" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Black-Scholes implied volatility\n", "\n", "We also need the Black-Scholes call price and its numerical inversion,\n", "both for computing the initial variance $v_0$ before calibration and for\n", "plotting the implied volatility surface afterward." ], "id": "a814f953-b4dc-4ded-9c0c-3bda74a148e1" }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "def bs_call(S, K, T, r, q, sigma):\n", " d1 = (log(S / K) + (r - q + 0.5 * sigma**2) * T) / (sigma * T**0.5)\n", " d2 = d1 - sigma * T**0.5\n", " return S * exp(-q * T) * norm.cdf(d1) - K * exp(-r * T) * norm.cdf(d2)\n", "\n", "def bs_implied_vol(S, K, T, r, q, market_price):\n", " try:\n", " return brentq(lambda sig: bs_call(S, K, T, r, q, sig) - market_price,\n", " 1e-6, 5.0)\n", " except ValueError:\n", " return np.nan" ], "id": "aa3d56de" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Verification\n", "\n", "We verify the implementation with a standard set of parameter values\n", "before moving to real data." ], "id": "09926cf4-5382-4a69-bbbc-31e2bf2b549b" }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Analytical call price: 8.7526" ] } ], "source": [ "S0_test = 100\n", "K_test = 100\n", "T_test = 1.0\n", "r_test = 0.05\n", "q_test = 0.02\n", "kappa_t = 2.0\n", "theta_t = 0.04\n", "sigma_t = 0.5\n", "rho_t = -0.7\n", "V0_test = 0.04\n", "\n", "call_test = heston_call_price(S0_test, V0_test, K_test, T_test,\n", " r_test, q_test, kappa_t, theta_t, sigma_t, rho_t)\n", "print(f\"Analytical call price: {call_test:.4f}\")" ], "id": "d0b79d19" }, { "cell_type": "markdown", "metadata": {}, "source": [ "We cross-check the analytical price against a Monte Carlo estimate\n", "obtained by simulating the Heston dynamics under the risk-neutral\n", "measure. The stock price and variance follow $$\n", "\\begin{aligned}\n", " dS &= (r - q)\\, S\\, dt + \\sqrt{v}\\, S\\, dW_1, \\\\\n", " dv &= \\kappa(\\theta - v)\\, dt + \\sigma\\sqrt{v}\\, dW_2,\n", "\\end{aligned}\n", "$$ where $dW_1\\, dW_2 = \\rho\\, dt$. We discretize with $n$ equal time\n", "steps using the Euler–Maruyama scheme, reflect negative variance to zero\n", "at each step, and estimate the call price as the discounted average\n", "payoff $e^{-rT} \\mathbb{E}[\\max(S(T)-K,0)]$." ], "id": "b95167c0-18c8-438d-984b-faab09322028" }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Analytical call price: 8.7526\n", "Monte Carlo estimate: 8.7832 (±0.0456 at 95%)\n", "Difference: 0.0306" ] } ], "source": [ "def heston_mc_call(S0, V0, K, T, r, q, kappa, theta, sigma, rho,\n", " n_paths=200_000, n_steps=200, seed=42):\n", " rng = np.random.default_rng(seed)\n", " dt = T / n_steps\n", " sqrt_dt = dt**0.5\n", "\n", " S = np.full(n_paths, float(S0))\n", " v = np.full(n_paths, float(V0))\n", "\n", " for _ in range(n_steps):\n", " Z1 = rng.standard_normal(n_paths)\n", " Z2 = rng.standard_normal(n_paths)\n", " W1 = Z1\n", " W2 = rho * Z1 + (1 - rho**2)**0.5 * Z2\n", "\n", " v_pos = np.maximum(v, 0.0)\n", " sv = v_pos**0.5\n", "\n", " S = S * np.exp((r - q - 0.5 * v_pos) * dt + sv * sqrt_dt * W1)\n", " v = v + kappa * (theta - v_pos) * dt + sigma * sv * sqrt_dt * W2\n", "\n", " payoff = np.maximum(S - K, 0.0)\n", " price = np.exp(-r * T) * payoff.mean()\n", " se = np.exp(-r * T) * payoff.std() / n_paths**0.5\n", " return price, se\n", "\n", "mc_price, mc_se = heston_mc_call(S0_test, V0_test, K_test, T_test,\n", " r_test, q_test, kappa_t, theta_t, sigma_t, rho_t)\n", "print(f\"Analytical call price: {call_test:.4f}\")\n", "print(f\"Monte Carlo estimate: {mc_price:.4f} (±{2*mc_se:.4f} at 95%)\")\n", "print(f\"Difference: {abs(call_test - mc_price):.4f}\")" ], "id": "90f846ab" }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Monte Carlo estimate should be within two standard errors of the\n", "analytical price, confirming that the characteristic-function formula\n", "and the SDE discretization are consistent.\n", "\n", "## AAPL Option Data\n", "\n", "The option prices used here come from the OptionMetrics (OPTIONM)\n", "database. OPTIONM records the highest closing bid and lowest closing ask\n", "across all exchanges for every listed contract. We use the midpoint\n", "$\\frac{\\text{bid} + \\text{ask}}{2}$ as the market price for each\n", "contract. The midpoint is the standard choice in empirical options\n", "research because it is free from the direction of trading and minimizes\n", "the impact of the bid-ask bounce on estimated parameters.\n", "\n", "AAPL options are American-style. For call options with a small dividend\n", "yield the early-exercise premium is negligible, so we treat the observed\n", "prices as European call prices without adjustment. On October 1, 2024,\n", "AAPL closed at \\$226.21. We use a dividend yield of 0.5% per year and a\n", "continuously compounded risk-free rate of 5%." ], "id": "b17700ba-cc58-43b8-beed-c47e37de800e" }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "S0 = 226.21\n", "r = 0.05\n", "q = 0.005" ], "id": "0c217526" }, { "cell_type": "markdown", "metadata": {}, "source": [ "We select seven strikes spanning roughly $\\pm$ 15% around the current\n", "stock price and six standard expiries ranging from 45 days to 198 days,\n", "giving 42 contracts in total. Time-to-maturity is measured in years from\n", "October 1, 2024." ], "id": "9da9cd0a-643b-456c-b404-bac4e8918c49" }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "strikes = [200, 210, 220, 230, 240, 250, 260]\n", "\n", "# Expiration dates and exact days to expiry from 2024-10-01\n", "# Nov 15 2024: 45 days; Dec 20 2024: 80 days; Jan 17 2025: 108 days;\n", "# Feb 21 2025: 143 days; Mar 21 2025: 171 days; Apr 17 2025: 198 days\n", "maturities = [45/365, 80/365, 108/365, 143/365, 171/365, 198/365]\n", "maturity_labels = ['Nov-24', 'Dec-24', 'Jan-25', 'Feb-25', 'Mar-25', 'Apr-25']" ], "id": "501981df" }, { "cell_type": "markdown", "metadata": {}, "source": [ "The midpoint prices below are taken directly from OPTIONM for the\n", "contracts described above." ], "id": "d2c69aa3-7b76-474c-9b4e-1db6ae65ea0d" }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "# Midpoint prices: rows = strikes (200 to 260), columns = maturities (Nov-24 to Apr-25)\n", "market_prices = np.array([\n", " [28.950, 30.775, 32.125, 33.925, 35.325, 36.475], # K = 200\n", " [20.450, 22.625, 24.175, 26.300, 27.725, 29.025], # K = 210\n", " [13.050, 15.475, 17.075, 19.475, 21.000, 22.400], # K = 220\n", " [ 7.275, 9.625, 11.250, 13.725, 15.225, 16.625], # K = 230\n", " [ 3.400, 5.350, 6.800, 9.075, 10.525, 11.850], # K = 240\n", " [ 1.395, 2.680, 3.750, 5.700, 6.975, 8.175], # K = 250\n", " [ 0.540, 1.260, 1.970, 3.425, 4.450, 5.450], # K = 260\n", "])" ], "id": "33729889" }, { "cell_type": "markdown", "metadata": {}, "source": [ "We display the data as a table before calibrating." ], "id": "cbb06dc7-e064-47c7-b083-a98df0605146" }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "df_market = pd.DataFrame(\n", " market_prices,\n", " index=pd.Index(strikes, name='Strike'),\n", " columns=maturity_labels,\n", ")\n", "df_market" ], "id": "68ce6d42-3c1e-48c1-9f03-2fa0badbfa22" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Calibration\n", "\n", "We calibrate four Heston parameters $(\\kappa, \\theta, \\sigma, \\rho)$ by\n", "minimizing the RMSE between observed and model call prices, holding\n", "$v_0$ fixed at a value extracted directly from market data.\n", "\n", "### Fixing the initial variance\n", "\n", "Treating $v_0$ as a free parameter creates a degeneracy: a very large\n", "$\\kappa$ paired with a $v_0$ far from the short-term implied variance\n", "can produce nearly the same option prices as a moderate $\\kappa$ with a\n", "well-anchored $v_0$, because the two parameters compensate each other.\n", "Fixing $v_0$ breaks this degeneracy, reduces the problem from five to\n", "four free parameters, and anchors the short end of the volatility\n", "surface to a directly observable market quantity.\n", "\n", "We set $v_0$ equal to the squared Black-Scholes implied volatility of\n", "the shortest-maturity contract whose strike is nearest to the forward\n", "price $F = S e^{(r-q)T}$." ], "id": "921354a7-942c-45a7-87b0-d81d3e786214" }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "ATM strike (Nov-24): 230\n", "ATM implied vol: 0.2662 (26.62%)\n", "V0 fixed to: 0.070862" ] } ], "source": [ "T_short = maturities[0]\n", "F_short = S0 * exp((r - q) * T_short)\n", "atm_idx = int(np.argmin(np.abs(np.array(strikes) - F_short)))\n", "iv_atm = bs_implied_vol(S0, strikes[atm_idx], T_short, r, q, market_prices[atm_idx, 0])\n", "V0_fixed = iv_atm**2\n", "\n", "print(f\"ATM strike (Nov-24): {strikes[atm_idx]}\")\n", "print(f\"ATM implied vol: {iv_atm:.4f} ({iv_atm*100:.2f}%)\")\n", "print(f\"V0 fixed to: {V0_fixed:.6f}\")" ], "id": "b2c86ae1" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Running the optimizer\n", "\n", "The objective function loops over all 42 contracts using the fixed\n", "$v_0$." ], "id": "7cb438b3-3f4f-4439-9c13-37a77379146b" }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "def rmse_heston(params):\n", " kappa, theta, sigma, rho = params\n", " total = 0.0\n", " for i, K in enumerate(strikes):\n", " for j, T in enumerate(maturities):\n", " model = heston_call_price(S0, V0_fixed, K, T, r, q, kappa, theta, sigma, rho)\n", " total += (market_prices[i, j] - model) ** 2\n", " return np.sqrt(total / (len(strikes) * len(maturities)))" ], "id": "b969352e" }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "kappa = 3.4879\n", "theta = 0.0627 (sqrt: 0.2504)\n", "sigma = 0.8174\n", "rho = -0.4823\n", "V0 = 0.0709 (sqrt: 0.2662) [fixed]\n", "RMSE = 0.2070" ] } ], "source": [ "bounds = [\n", " (0.001, 20.0), # kappa\n", " (0.001, 1.0), # theta\n", " (0.001, 5.0), # sigma\n", " (-1.0, 1.0), # rho\n", "]\n", "\n", "x0 = [2.0, V0_fixed, 0.5, -0.7]\n", "\n", "result = minimize(rmse_heston, x0, method='L-BFGS-B', bounds=bounds)\n", "\n", "kappa_cal, theta_cal, sigma_cal, rho_cal = result.x\n", "\n", "print(f\"kappa = {kappa_cal:.4f}\")\n", "print(f\"theta = {theta_cal:.4f} (sqrt: {theta_cal**0.5:.4f})\")\n", "print(f\"sigma = {sigma_cal:.4f}\")\n", "print(f\"rho = {rho_cal:.4f}\")\n", "print(f\"V0 = {V0_fixed:.4f} (sqrt: {V0_fixed**0.5:.4f}) [fixed]\")\n", "print(f\"RMSE = {result.fun:.4f}\")" ], "id": "d517b138" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Interpreting the calibrated parameters\n", "\n", "The four calibrated parameters have natural economic interpretations.\n", "\n", "The speed of mean reversion $\\kappa$ controls how quickly the\n", "instantaneous variance returns to its long-run level $\\theta$. Larger\n", "values imply faster reversion; as $\\kappa \\to \\infty$ the variance\n", "process becomes almost deterministic and the model collapses toward a\n", "constant-volatility world.\n", "\n", "The long-run variance $\\theta$ determines the volatility level that\n", "options with long maturities are priced against. Its square root\n", "$\\sqrt{\\theta}$ is directly comparable to the implied volatility of\n", "long-dated at-the-money options.\n", "\n", "The volatility of volatility $\\sigma$ governs the curvature of the\n", "implied volatility surface. A larger $\\sigma$ produces a more pronounced\n", "smile because the variance process itself is more uncertain.\n", "\n", "The correlation $\\rho$ captures the leverage effect: negative values\n", "make price drops and variance increases co-move, which tilts the implied\n", "volatility surface into the downward-sloping skew that is characteristic\n", "of equity markets. The sign should be negative for AAPL.\n", "\n", "## Implied Volatility Analysis\n", "\n", "We compute the Black-Scholes implied volatility for every contract in\n", "the data set." ], "id": "c7e5badf-769b-4ab1-8864-1ea95c55e9ad" }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "iv = np.zeros_like(market_prices)\n", "for i, K in enumerate(strikes):\n", " for j, T in enumerate(maturities):\n", " iv[i, j] = bs_implied_vol(S0, K, T, r, q, market_prices[i, j])" ], "id": "e8cf384d" }, { "cell_type": "markdown", "metadata": {}, "source": [ "The figure below shows the implied volatility surface. The axes are\n", "strike and days to expiry; the height and color show the level of\n", "implied volatility." ], "id": "8cded878-0300-419d-9506-f323a6be8199" }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": {} } ], "source": [ "from mpl_toolkits.mplot3d import Axes3D\n", "\n", "K_grid = np.array(strikes)\n", "T_grid = np.array(maturities) * 365 # days to expiry for the axis label\n", "\n", "KK, TT = np.meshgrid(K_grid, T_grid)\n", "\n", "fig = plt.figure(figsize=(8, 5))\n", "ax = fig.add_subplot(111, projection='3d')\n", "surf = ax.plot_surface(KK, TT, iv.T * 100, cmap='inferno', alpha=0.9,\n", " linewidth=0.3, edgecolor='0.4')\n", "fig.colorbar(surf, ax=ax, shrink=0.5, aspect=10, pad=0.1, label='IV (%)')\n", "ax.set_xlabel('Strike')\n", "ax.set_ylabel('Days to Expiry')\n", "ax.set_zlabel('Implied Vol (%)')\n", "ax.set_title('AAPL Implied Volatility Surface — October 1, 2024')\n", "plt.tight_layout()\n", "plt.show()" ], "id": "cell-fig-iv-surface" }, { "cell_type": "markdown", "metadata": {}, "source": [ "The surface shows the typical equity volatility skew: implied\n", "volatilities are highest for low strikes and decrease monotonically as\n", "the strike rises. The skew is steeper for short maturities and flattens\n", "as maturity increases, which is exactly the pattern produced by the\n", "Heston model with a negative $\\rho$.\n", "\n", "By construction, $\\sqrt{v_0}$ equals the implied volatility of the\n", "nearest-ATM Nov-24 option, so the short-dated ATM level is matched\n", "exactly. Long-dated options are priced closer to $\\sqrt{\\theta}$,\n", "reflecting the mean reversion of $v$ toward its long-run level.\n", "\n", "## Pricing Errors\n", "\n", "We compute model prices and pricing errors for all 42 contracts." ], "id": "98acce99-7884-4506-9188-85ec5a72dbc0" }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "model_prices = np.zeros_like(market_prices)\n", "for i, K in enumerate(strikes):\n", " for j, T in enumerate(maturities):\n", " model_prices[i, j] = heston_call_price(\n", " S0, V0_fixed, K, T, r, q, kappa_cal, theta_cal, sigma_cal, rho_cal\n", " )\n", "\n", "abs_errors = np.abs(market_prices - model_prices)" ], "id": "9e74d478" }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "rows = []\n", "for i, K in enumerate(strikes):\n", " for j, (label, T) in enumerate(zip(maturity_labels, maturities)):\n", " rows.append({\n", " 'Strike': K,\n", " 'Expiry': label,\n", " 'Observed': f\"{market_prices[i, j]:.3f}\",\n", " 'Model': f\"{model_prices[i, j]:.3f}\",\n", " 'Abs. Error': f\"{abs_errors[i, j]:.3f}\",\n", " })\n", "\n", "df_errors = pd.DataFrame(rows).set_index(['Strike', 'Expiry'])\n", "df_errors" ], "id": "c584c4ef-f2ab-4b1a-badc-5b5baf22f50b" }, { "cell_type": "markdown", "metadata": {}, "source": [ "We summarize the mean absolute error for each strike-expiry combination,\n", "with marginal averages by row and column." ], "id": "1c56fac2-0c2b-4147-98ca-6611ff295c41" }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "df_mae = pd.DataFrame(\n", " abs_errors.round(4),\n", " index=pd.Index(strikes, name='Strike'),\n", " columns=maturity_labels,\n", ")\n", "df_mae['Avg'] = abs_errors.mean(axis=1).round(4)\n", "df_mae.loc['Avg'] = abs_errors.mean(axis=0).tolist() + [abs_errors.mean().round(4)]\n", "df_mae" ], "id": "4ab36b58-c598-4553-ae15-9967bb577bc8" }, { "cell_type": "markdown", "metadata": {}, "source": [ "The figure below plots the pricing errors by strike for each expiry." ], "id": "2f924915-403b-4ba8-8154-16b96b06ad79" }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": {} } ], "source": [ "pricing_errors = model_prices - market_prices\n", "\n", "fig, axes = plt.subplots(2, 3, figsize=(9, 5), sharey=True)\n", "for j, (label, ax) in enumerate(zip(maturity_labels, axes.flat)):\n", " ax.bar(strikes, pricing_errors[:, j], color=['#d62728' if e > 0 else '#1f77b4'\n", " for e in pricing_errors[:, j]])\n", " ax.axhline(0, color='0.3', lw=0.8, ls='--')\n", " ax.set_title(label, fontsize=10)\n", " ax.set_xticks(strikes)\n", " ax.set_xticklabels(strikes, fontsize=7)\n", " ax.spines[['top', 'right']].set_visible(False)\n", "\n", "fig.supxlabel('Strike', fontsize=9)\n", "fig.supylabel('Error ($)', fontsize=9)\n", "plt.tight_layout()\n", "plt.show()" ], "id": "cell-fig-pricing-errors" }, { "cell_type": "markdown", "metadata": {}, "source": [ "The largest errors occur in short-maturity contracts. This is a\n", "structural limitation of the Heston model: the model-implied skew is\n", "proportional to $\\sqrt{T}$ as maturity shrinks, so it flattens too\n", "quickly near expiry and cannot match the steep short-dated skew observed\n", "in the market. The RMSE treats all contracts equally, but near-the-money\n", "options have larger dollar prices, so their squared errors receive more\n", "implicit weight and the optimizer tilts the fit toward those contracts\n", "across all maturities.\n", "\n", "The overall fit demonstrates that the Heston model can reproduce the\n", "main features of the equity implied volatility surface—the level, the\n", "skew, and the term structure of volatility—with just four free\n", "parameters.\n", "\n", "Cortazar, Gonzalo, Matias Lopez, and Lorenzo Naranjo. 2017. “A\n", "Multifactor Stochastic Volatility Model of Commodity Prices.” *Energy\n", "Economics* 67: 182–201." ], "id": "a489c856-5a51-40a2-9d94-3d18a0e99da1" } ], "nbformat": 4, "nbformat_minor": 5, "metadata": { "kernelspec": { "name": "python3", "display_name": "Python 3 (ipykernel)", "language": "python", "path": "/home/lnaranjo/.pyenv/versions/3.14.7/share/jupyter/kernels/python3" }, "language_info": { "name": "python", "codemirror_mode": { "name": "ipython", "version": "3" }, "file_extension": ".py", "mimetype": "text/x-python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.14.7" } } }