{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "48c1f622",
   "metadata": {},
   "source": [
    "# Notebook 02 — Geometric Brownian Motion & Itô's Lemma\n",
    "\n",
    "In Notebook 01 we modelled the **log-price** as a Brownian motion. But the\n",
    "question is about the **price** $S_t$ itself, and a price has one\n",
    "non-negotiable property: **it must stay positive**. An arithmetic Brownian\n",
    "motion $X_t = \\nu t + \\sigma W_t$ goes negative with positive probability, so\n",
    "it cannot *be* a price. The fix is the most important model in quantitative\n",
    "finance: **geometric Brownian motion (GBM)**, and the tool that lets us do\n",
    "calculus on it is **Itô's lemma**.\n",
    "\n",
    "### What this notebook delivers\n",
    "\n",
    "1. Why prices demand a *multiplicative* (geometric) model.\n",
    "2. The GBM stochastic differential equation (SDE).\n",
    "3. **Itô's lemma** — the chain rule of stochastic calculus — derived with `sympy`.\n",
    "4. The exact solution for $S_t$, and the surprising \"$-\\tfrac12\\sigma^2$\"\n",
    "   correction that quants call *volatility drag*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "25e9993c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-08-02T16:21:48.315815Z",
     "iopub.status.busy": "2026-08-02T16:21:48.315286Z",
     "iopub.status.idle": "2026-08-02T16:21:50.210877Z",
     "shell.execute_reply": "2026-08-02T16:21:50.209480Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import sympy as sp\n",
    "\n",
    "# =====================================================================\n",
    "# Helper functions -- self-contained small library,\n",
    "# included verbatim at the top of every notebook so each one runs standalone.\n",
    "# (This cell is identical across the series.)\n",
    "# =====================================================================\n",
    "import types\n",
    "from scipy.stats import norm\n",
    "import matplotlib as mpl\n",
    "\n",
    "# --- config: trading-time & volatility conventions (annualised sigma) ---\n",
    "TRADING_DAYS_PER_YEAR = 252\n",
    "TRADING_HOURS_PER_DAY = 6.5\n",
    "TRADING_MINUTES_PER_HOUR = 60\n",
    "TRADING_MINUTES_PER_DAY = TRADING_HOURS_PER_DAY * TRADING_MINUTES_PER_HOUR     # 390\n",
    "TRADING_MINUTES_PER_YEAR = TRADING_DAYS_PER_YEAR * TRADING_MINUTES_PER_DAY      # 98280\n",
    "DROP_FRACTION = 0.10                       # a 10% drop\n",
    "WINDOW_MINUTES = 5                         # over a 5-minute window\n",
    "WINDOW_YEARS = WINDOW_MINUTES / TRADING_MINUTES_PER_YEAR    # ~5.09e-5\n",
    "LOG_BARRIER = float(np.log(1.0 - DROP_FRACTION))            # ln(0.9) ~ -0.10536\n",
    "TYPICAL_ANNUAL_VOL = 0.30                  # representative liquid large-cap volatility\n",
    "\n",
    "def minutes_to_years(minutes):\n",
    "    return minutes / TRADING_MINUTES_PER_YEAR\n",
    "\n",
    "def horizon_vol(sigma_annual, minutes):\n",
    "    \"\"\"sigma over `minutes` of trading time via Brownian scaling sigma*sqrt(T).\"\"\"\n",
    "    return sigma_annual * np.sqrt(minutes_to_years(minutes))\n",
    "\n",
    "def trading_minutes_to_seconds(minutes):\n",
    "    return minutes * 60.0\n",
    "\n",
    "# --- processes: samplers; each returns (t, X) with X of shape (n_paths, n_steps+1) ---\n",
    "def _as_rng(rng):\n",
    "    return np.random.default_rng() if rng is None else rng\n",
    "\n",
    "def time_grid(T, n_steps):\n",
    "    return np.linspace(0.0, T, n_steps + 1)\n",
    "\n",
    "def brownian_motion(T, n_steps, n_paths=1, drift=0.0, sigma=1.0, rng=None):\n",
    "    \"\"\"Arithmetic BM, dX = drift*dt + sigma*dW, via independent Gaussian increments.\"\"\"\n",
    "    rng = _as_rng(rng)\n",
    "    t = time_grid(T, n_steps)\n",
    "    dt = T / n_steps\n",
    "    dW = rng.normal(loc=drift * dt, scale=sigma * np.sqrt(dt), size=(n_paths, n_steps))\n",
    "    X = np.zeros((n_paths, n_steps + 1))\n",
    "    X[:, 1:] = np.cumsum(dW, axis=1)\n",
    "    return t, X\n",
    "\n",
    "def brownian_motion_cholesky(T, n_steps, n_paths=1, rng=None):\n",
    "    \"\"\"Standard BM via the min(s,t) covariance (Gram) matrix and its Cholesky factor L (x = L z).\"\"\"\n",
    "    rng = _as_rng(rng)\n",
    "    t = time_grid(T, n_steps)\n",
    "    t_inner = t[1:]\n",
    "    K = np.minimum(t_inner[:, None], t_inner[None, :])     # min(s,t) Gram matrix\n",
    "    L = np.linalg.cholesky(K)\n",
    "    z = rng.standard_normal(size=(n_paths, n_steps))\n",
    "    X = np.zeros((n_paths, n_steps + 1))\n",
    "    X[:, 1:] = z @ L.T\n",
    "    return t, X\n",
    "\n",
    "def geometric_brownian_motion(S0, T, n_steps, n_paths=1, mu=0.0, sigma=1.0, rng=None):\n",
    "    \"\"\"GBM, dS = mu*S*dt + sigma*S*dW, via the exact log-space solution (no Euler error).\"\"\"\n",
    "    rng = _as_rng(rng)\n",
    "    t, X = brownian_motion(T, n_steps, n_paths, drift=mu - 0.5 * sigma**2, sigma=sigma, rng=rng)\n",
    "    return t, S0 * np.exp(X)\n",
    "\n",
    "def brownian_bridge(T, n_steps, n_paths=1, start=0.0, end=0.0, sigma=1.0, rng=None):\n",
    "    \"\"\"BM conditioned on B(0)=start, B(T)=end.\"\"\"\n",
    "    rng = _as_rng(rng)\n",
    "    t, W = brownian_motion(T, n_steps, n_paths, drift=0.0, sigma=sigma, rng=rng)\n",
    "    linear = start + (end - start) * (t / T)\n",
    "    bridge = W - np.outer(W[:, -1], t / T) + linear[None, :]\n",
    "    return t, bridge\n",
    "\n",
    "def merton_jump_diffusion(S0, T, n_steps, n_paths=1, mu=0.0, sigma=0.2,\n",
    "                          jump_intensity=0.0, jump_mean=-0.10, jump_sigma=0.15, rng=None):\n",
    "    \"\"\"Merton jump-diffusion; jump_intensity=0 reduces exactly to GBM.\"\"\"\n",
    "    rng = _as_rng(rng)\n",
    "    dt = T / n_steps\n",
    "    k = np.exp(jump_mean + 0.5 * jump_sigma**2) - 1.0          # E[J-1]\n",
    "    nu = mu - 0.5 * sigma**2 - jump_intensity * k               # compensated drift\n",
    "    t, X = brownian_motion(T, n_steps, n_paths, drift=nu, sigma=sigma, rng=rng)\n",
    "    counts = rng.poisson(jump_intensity * dt, size=(n_paths, n_steps))\n",
    "    jump_log = np.zeros((n_paths, n_steps))\n",
    "    K = int(counts.max()) if counts.size else 0\n",
    "    if K > 0:\n",
    "        logJ = rng.normal(jump_mean, jump_sigma, size=(n_paths, n_steps, K))\n",
    "        mask = np.arange(K) < counts[:, :, None]\n",
    "        jump_log = (logJ * mask).sum(axis=2)\n",
    "    jump_path = np.zeros((n_paths, n_steps + 1))\n",
    "    jump_path[:, 1:] = np.cumsum(jump_log, axis=1)\n",
    "    return t, S0 * np.exp(X + jump_path)\n",
    "\n",
    "# --- barriers: first-passage CDFs & the Brownian-bridge miss formula (log-price BM, B<0) ---\n",
    "def first_passage_cdf_zero_drift(B, T, sigma):\n",
    "    \"\"\"P(tau_B<=T) for zero-drift BM, via the reflection principle: 2*Phi(B/(sigma*sqrt(T))).\"\"\"\n",
    "    return 2.0 * norm.cdf(B / (sigma * np.sqrt(T)))\n",
    "\n",
    "def first_passage_cdf(B, T, nu, sigma):\n",
    "    \"\"\"P(tau_B<=T) for BM with drift nu (Bachelier-Levy).\"\"\"\n",
    "    sT = sigma * np.sqrt(T)\n",
    "    return norm.cdf((B - nu * T) / sT) + np.exp(2.0 * nu * B / sigma**2) * norm.cdf((B + nu * T) / sT)\n",
    "\n",
    "def barrier_miss_prob(x0, xT, B, sigma, dt):\n",
    "    \"\"\"P(BM crossed B inside (0,dt) | both endpoints above B): exp(-2(x0-B)(xT-B)/(sigma^2 dt)).\"\"\"\n",
    "    if x0 <= B or xT <= B:\n",
    "        return 1.0\n",
    "    return float(np.exp(-2.0 * (x0 - B) * (xT - B) / (sigma**2 * dt)))\n",
    "\n",
    "def max_safe_dt(D, sigma, epsilon):\n",
    "    \"\"\"Largest sampling interval keeping P(miss)<=epsilon: dt = 2*D^2/(sigma^2 * ln(1/epsilon)).\"\"\"\n",
    "    if not (0.0 < epsilon < 1.0):\n",
    "        raise ValueError(\"epsilon must lie in (0, 1)\")\n",
    "    return 2.0 * D**2 / (sigma**2 * np.log(1.0 / epsilon))\n",
    "\n",
    "def estimate_first_passage_prob(B, T, nu, sigma, n_paths, n_steps, rng=None):\n",
    "    \"\"\"Naive MC of P(tau_B<=T) via discrete grid checks (O(sqrt(dt)) downward bias).\"\"\"\n",
    "    rng = _as_rng(rng)\n",
    "    _, X = brownian_motion(T, n_steps, n_paths, drift=nu, sigma=sigma, rng=rng)\n",
    "    return float(np.mean(np.any(X <= B, axis=1)))\n",
    "\n",
    "def estimate_first_passage_prob_bb(B, T, nu, sigma, n_paths, n_steps, rng=None):\n",
    "    \"\"\"Bias-free MC of P(tau_B<=T) using the Brownian-bridge crossing probability per step.\"\"\"\n",
    "    rng = _as_rng(rng)\n",
    "    _, X = brownian_motion(T, n_steps, n_paths, drift=nu, sigma=sigma, rng=rng)\n",
    "    dt = T / n_steps\n",
    "    x0, x1 = X[:, :-1], X[:, 1:]\n",
    "    both_above = (x0 > B) & (x1 > B)\n",
    "    p_cross = np.where(both_above, np.exp(-2.0 * (x0 - B) * (x1 - B) / (sigma**2 * dt)), 1.0)\n",
    "    u = rng.uniform(size=p_cross.shape)\n",
    "    return float(np.mean(np.any(u < p_cross, axis=1)))\n",
    "\n",
    "def estimate_bridge_breach_prob(x0, xT, B, sigma, dt, n_paths, n_inner, rng=None):\n",
    "    \"\"\"MC check of barrier_miss_prob by simulating bridges pinned at x0, xT.\"\"\"\n",
    "    rng = _as_rng(rng)\n",
    "    _, bridge = brownian_bridge(dt, n_inner, n_paths, start=x0, end=xT, sigma=sigma, rng=rng)\n",
    "    return float(np.mean(np.any(bridge <= B, axis=1)))\n",
    "\n",
    "# --- plotting helpers ---\n",
    "def style():\n",
    "    \"\"\"Apply the project's matplotlib style (serif maths, clean spines).\"\"\"\n",
    "    mpl.rcParams.update({\n",
    "        \"figure.figsize\": (8.0, 4.5), \"figure.dpi\": 110, \"savefig.dpi\": 140,\n",
    "        \"axes.grid\": True, \"grid.alpha\": 0.3,\n",
    "        \"axes.spines.top\": False, \"axes.spines.right\": False,\n",
    "        \"font.size\": 11, \"mathtext.fontset\": \"cm\",\n",
    "        \"axes.titlesize\": 12, \"axes.labelsize\": 11, \"legend.frameon\": False,\n",
    "    })\n",
    "\n",
    "def plot_paths(t, X, ax=None, alpha=0.5, lw=0.8, color=\"C0\", label=None, **kwargs):\n",
    "    \"\"\"Plot each row of X (shape (n_paths, len(t))) against t.\"\"\"\n",
    "    if ax is None:\n",
    "        _, ax = plt.subplots()\n",
    "    X = np.atleast_2d(X)\n",
    "    for i, row in enumerate(X):\n",
    "        ax.plot(t, row, color=color, alpha=alpha, lw=lw, label=label if i == 0 else None, **kwargs)\n",
    "    ax.set_xlabel(\"time\")\n",
    "    return ax\n",
    "\n",
    "def hline(ax, y, color=\"crimson\", label=None, ls=\"--\", lw=1.4):\n",
    "    \"\"\"Draw a horizontal barrier line.\"\"\"\n",
    "    ax.axhline(y, color=color, ls=ls, lw=lw, label=label)\n",
    "    return ax\n",
    "\n",
    "# --- module namespaces (so existing call sites like processes.X / config.Y keep working) ---\n",
    "config = types.SimpleNamespace(\n",
    "    TRADING_DAYS_PER_YEAR=TRADING_DAYS_PER_YEAR, TRADING_HOURS_PER_DAY=TRADING_HOURS_PER_DAY,\n",
    "    TRADING_MINUTES_PER_HOUR=TRADING_MINUTES_PER_HOUR, TRADING_MINUTES_PER_DAY=TRADING_MINUTES_PER_DAY,\n",
    "    TRADING_MINUTES_PER_YEAR=TRADING_MINUTES_PER_YEAR, DROP_FRACTION=DROP_FRACTION,\n",
    "    WINDOW_MINUTES=WINDOW_MINUTES, WINDOW_YEARS=WINDOW_YEARS, LOG_BARRIER=LOG_BARRIER,\n",
    "    TYPICAL_ANNUAL_VOL=TYPICAL_ANNUAL_VOL, minutes_to_years=minutes_to_years,\n",
    "    horizon_vol=horizon_vol, trading_minutes_to_seconds=trading_minutes_to_seconds,\n",
    ")\n",
    "processes = types.SimpleNamespace(\n",
    "    brownian_motion=brownian_motion, brownian_motion_cholesky=brownian_motion_cholesky,\n",
    "    geometric_brownian_motion=geometric_brownian_motion, brownian_bridge=brownian_bridge,\n",
    "    merton_jump_diffusion=merton_jump_diffusion,\n",
    ")\n",
    "barriers = types.SimpleNamespace(\n",
    "    first_passage_cdf_zero_drift=first_passage_cdf_zero_drift, first_passage_cdf=first_passage_cdf,\n",
    "    barrier_miss_prob=barrier_miss_prob, max_safe_dt=max_safe_dt,\n",
    "    estimate_first_passage_prob=estimate_first_passage_prob,\n",
    "    estimate_first_passage_prob_bb=estimate_first_passage_prob_bb,\n",
    "    estimate_bridge_breach_prob=estimate_bridge_breach_prob,\n",
    ")\n",
    "plotting = types.SimpleNamespace(style=style, plot_paths=plot_paths, hline=hline)\n",
    "\n",
    "from pathlib import Path\n",
    "\n",
    "def _img_subdir(sub):\n",
    "    \"\"\"Resolve <repo>/images/<sub>, creating it; cwd-independent.\"\"\"\n",
    "    root = Path.cwd()\n",
    "    for cand in (root, *root.parents):\n",
    "        if (cand / \"notebooks\").is_dir():\n",
    "            root = cand\n",
    "            break\n",
    "    d = root / \"images\" / sub\n",
    "    d.mkdir(parents=True, exist_ok=True)\n",
    "    return d\n",
    "\n",
    "_FIG = _img_subdir(\"02_gbm_ito\")\n",
    "\n",
    "plotting.style()\n",
    "RNG = np.random.default_rng(2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c308f848",
   "metadata": {},
   "source": [
    "## 1. Why the model must be multiplicative\n",
    "\n",
    "If we tried $S_t = S_0 + \\nu t + \\sigma W_t$ (arithmetic Brownian motion), then\n",
    "for any horizon $T$ there is a strictly positive probability that\n",
    "$S_T < 0$. A negative stock price is absurd, so this model is *disqualified*.\n",
    "\n",
    "The economic intuition is that returns — not absolute price changes — should be\n",
    "stationary. A $\\$1$ move in a $\\$10$ stock and in a $\\$1000$ stock are very\n",
    "different events. We therefore insist that the **instantaneous return** be\n",
    "\n",
    "$$\\frac{dS_t}{S_t} \\;=\\; \\text{(drift)}\\,dt + \\text{(noise)}\\,dW_t.$$\n",
    "\n",
    "Because $dS_t$ is proportional to $S_t$, if $S_t > 0$ it can never reach zero\n",
    "(or below) in finite time — positivity is preserved.\n",
    "\n",
    "Let's see the difference in simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "369a5a72",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-08-02T16:21:50.213694Z",
     "iopub.status.busy": "2026-08-02T16:21:50.213328Z",
     "iopub.status.idle": "2026-08-02T16:21:50.926702Z",
     "shell.execute_reply": "2026-08-02T16:21:50.925526Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 340304 characters]",
      "text/plain": [
       "<Figure size 990x462 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "min arithmetic-BM price: 17.62   -> crosses zero: False\n",
      "min geometric-BM price : 48.81   -> always positive\n"
     ]
    }
   ],
   "source": [
    "# Same noise, two models: arithmetic BM (can go negative) vs geometric BM.\n",
    "T, n_steps, n_paths = 1.0, 500, 30\n",
    "mu, sigma = 0.10, 0.30\n",
    "t, W = processes.brownian_motion(T, n_steps, n_paths, sigma=sigma, rng=RNG)\n",
    "S_arith = 100 * (1 + 0.10 * t + W)                 # arithmetic, can go negative\n",
    "S_geom  = processes.geometric_brownian_motion(100, T, n_steps, n_paths, mu, sigma, rng=RNG)[1]\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(9, 4.2))\n",
    "plotting.plot_paths(t, S_arith, ax=ax, alpha=0.5, color=\"C3\", label=\"arithmetic BM (can be $<0$)\")\n",
    "plotting.plot_paths(t, S_geom,  ax=ax, alpha=0.5, color=\"C0\", label=\"geometric BM (always $>0$)\")\n",
    "ax.axhline(0, color=\"k\", lw=0.8)\n",
    "ax.set_title(\"Positivity: arithmetic vs geometric Brownian motion\")\n",
    "ax.set_xlabel(\"time (years)\"); ax.set_ylabel(\"price\"); ax.legend(loc=\"upper left\")\n",
    "plt.savefig(_FIG / \"arithmetic_vs_gbm.png\", dpi=150, bbox_inches=\"tight\"); plt.show()\n",
    "\n",
    "print(\"min arithmetic-BM price:\", round(float(S_arith.min()), 2),\n",
    "      \"  -> crosses zero:\", bool((S_arith <= 0).any()))\n",
    "print(\"min geometric-BM price :\", round(float(S_geom.min()), 2),  \"  -> always positive\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de6fe76d",
   "metadata": {},
   "source": [
    "## 2. The Geometric Brownian Motion SDE\n",
    "\n",
    "Writing the noise amplitude on the *return* as a constant $\\sigma$ (the\n",
    "**volatility**) and the drift of the return as $\\mu$, we obtain the **GBM SDE**:\n",
    "\n",
    "$$\\boxed{\\;dS_t \\;=\\; \\mu\\, S_t\\, dt \\;+\\; \\sigma\\, S_t\\, dW_t\\;}$$\n",
    "\n",
    "The two parameters have units:\n",
    "\n",
    "* $\\mu$ — **drift**, the expected instantaneous *return* per unit time (annual);\n",
    "* $\\sigma$ — **volatility**, the instantaneous standard deviation of the return\n",
    "  per $\\sqrt{\\text{time}}$.\n",
    "\n",
    "> **Quant's note.** GBM is the model underlying the Black–Scholes option\n",
    "> formula. Its great virtues: positivity, closed-form solutions, and stationary\n",
    "> log-returns. Its great sin: it cannot reproduce **jumps** or **volatility\n",
    "> clustering**, which is exactly why real extreme moves lie outside it. We\n",
    "> return to this in Notebook 05.\n",
    "\n",
    "To do anything useful with this SDE — e.g. turn it into a closed-form price,\n",
    "or move to the log-price — we need a chain rule. That is Itô's lemma."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "feb6ff1b",
   "metadata": {},
   "source": [
    "## 3. Itô's lemma: the chain rule of stochastic calculus\n",
    "\n",
    "Suppose $S_t$ satisfies $dS_t = a(S_t,t)\\,dt + b(S_t,t)\\,dW_t$ and let\n",
    "$f(S_t, t)$ be a smooth function. The deterministic chain rule $df = f_S\\,dS +\n",
    "f_t\\,dt$ is **wrong** here. Why? Because Brownian increments are not $O(dt)$\n",
    "but $O(\\sqrt{dt})$: over $dt$ the typical move is $\\sigma\\sqrt{dt}$, so\n",
    "$(dS_t)^2 \\sim dt$, which is *not* negligible in a second-order Taylor\n",
    "expansion.\n",
    "\n",
    "### The derivation\n",
    "\n",
    "Expand $f$ to second order:\n",
    "\n",
    "$$df \\;=\\; f_t\\,dt + f_S\\,dS \\;+\\; \\tfrac12 f_{SS}\\,(dS)^2.$$\n",
    "\n",
    "Now substitute $dS = a\\,dt + b\\,dW$ and use the quadratic-variation facts\n",
    "\n",
    "$$dt \\cdot dt = 0, \\qquad dt \\cdot dW = 0, \\qquad (dW)^2 = dt.$$\n",
    "\n",
    "Keeping only terms of order $dt$:\n",
    "\n",
    "$$(dS)^2 \\;=\\; b^2\\,(dW)^2 \\;=\\; b^2\\,dt.$$\n",
    "\n",
    "So\n",
    "\n",
    "$$\\boxed{\\;df \\;=\\; \\Big(f_t + a\\,f_S + \\tfrac12 b^2 f_{SS}\\Big)\\,dt \\;+\\; b\\,f_S\\,dW_t\\;}$$\n",
    "\n",
    "The new ingredient versus ordinary calculus is the **$\\tfrac12 b^2 f_{SS}$\n",
    "term** — the *Itô correction*."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46064b09",
   "metadata": {},
   "source": [
    "### Itô's lemma, done symbolically with `sympy`\n",
    "\n",
    "Let's verify the algebra mechanically. We expand\n",
    "$f(S_t + dS) - f(S_t)$ to second order and substitute the multiplication\n",
    "table $(dt)^2 = dt\\,dW = 0$, $(dW)^2 = dt$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5be1f637",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-08-02T16:21:50.930201Z",
     "iopub.status.busy": "2026-08-02T16:21:50.929978Z",
     "iopub.status.idle": "2026-08-02T16:21:51.104896Z",
     "shell.execute_reply": "2026-08-02T16:21:51.103865Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Itô drift   df/dt coefficient: a(S, t)*Derivative(f(S, t), S) + b(S, t)**2*Derivative(f(S, t), (S, 2))/2\n",
      "Itô diffusion  dW coefficient: b(S, t)*Derivative(f(S, t), S)\n"
     ]
    }
   ],
   "source": [
    "# Symbolic Itô expansion for a general f(S, t) and dS = a dt + b dW.\n",
    "S, t, dW, dt = sp.symbols(\"S t dW dt\")\n",
    "a, b = sp.symbols(\"a b\", cls=sp.Function)   # a(S,t), b(S,t)\n",
    "f = sp.Function(\"f\")\n",
    "\n",
    "# First- and second-order Taylor terms of f around (S, t).\n",
    "f_S  = sp.diff(f(S, t), S)\n",
    "f_SS = sp.diff(f(S, t), S, 2)\n",
    "\n",
    "# Substitute dS = a(S,t) dt + b(S,t) dW, then apply the quadratic-variation rules.\n",
    "dS = a(S, t)*dt + b(S, t)*dW\n",
    "dS_squared = sp.expand(dS**2)\n",
    "rules = {dt**2: 0, dt*dW: 0, dW**2: dt}\n",
    "dS_squared = dS_squared.subs(rules)\n",
    "\n",
    "# Itô drift and diffusion.\n",
    "drift   = f_S * a(S, t) + sp.Rational(1, 2) * f_SS * b(S, t)**2\n",
    "diffuse = f_S * b(S, t)\n",
    "\n",
    "print(\"Itô drift   df/dt coefficient:\", sp.simplify(drift))\n",
    "print(\"Itô diffusion  dW coefficient:\", sp.simplify(diffuse))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1149ef5",
   "metadata": {},
   "source": [
    "## 4. Applying Itô to the log-price\n",
    "\n",
    "This is the pivotal calculation for the whole project. Take\n",
    "$f(S) = \\log S$ (no explicit $t$ dependence) applied to GBM, where\n",
    "$a(S,t) = \\mu S$ and $b(S,t) = \\sigma S$.\n",
    "\n",
    "$$f_S = \\frac{1}{S}, \\qquad f_{SS} = -\\frac{1}{S^2}.$$\n",
    "\n",
    "Plugging into Itô's formula:\n",
    "\n",
    "$$d(\\log S_t) \\;=\\; \\Big(\\mu S \\cdot \\tfrac{1}{S} + \\tfrac12 \\sigma^2 S^2 \\cdot \\big(-\\tfrac{1}{S^2}\\big)\\Big)\\,dt \\;+\\; \\sigma S\\cdot\\tfrac{1}{S}\\,dW_t,$$\n",
    "\n",
    "$$\\boxed{\\;dX_t \\;=\\; \\Big(\\mu - \\tfrac12 \\sigma^2\\Big)\\,dt \\;+\\; \\sigma\\, dW_t, \\qquad X_t := \\log S_t.\\;}$$\n",
    "\n",
    "The log-price is an **arithmetic Brownian motion** with drift\n",
    "$\\nu := \\mu - \\tfrac{1}{2}\\sigma^2$. The $-\\tfrac12\\sigma^2$ is the *Itô\n",
    "correction*; it is the price of applying a nonlinear function (log) to a\n",
    "process with non-zero quadratic variation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "be49c4cb",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-08-02T16:21:51.108753Z",
     "iopub.status.busy": "2026-08-02T16:21:51.108310Z",
     "iopub.status.idle": "2026-08-02T16:21:51.132966Z",
     "shell.execute_reply": "2026-08-02T16:21:51.131920Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "d(log S) drift   = mu - sigma**2/2  * dt\n",
      "d(log S) diffusion = sigma * dW\n",
      "\n",
      "=> nu = mu - sigma^2 / 2 : True\n"
     ]
    }
   ],
   "source": [
    "# Same calculation, with sympy, for GBM specifically.\n",
    "S, t, dW, dt, mu, sigma = sp.symbols(\"S t dW dt mu sigma\", positive=True)\n",
    "f = sp.log(S)\n",
    "f_S  = sp.diff(f, S)\n",
    "f_SS = sp.diff(f, S, 2)\n",
    "a_S, b_S = mu * S, sigma * S\n",
    "\n",
    "drift_X   = sp.simplify(a_S * f_S + sp.Rational(1, 2) * b_S**2 * f_SS)\n",
    "diffusion_X = sp.simplify(b_S * f_S)\n",
    "print(\"d(log S) drift   =\", drift_X, \" * dt\")\n",
    "print(\"d(log S) diffusion =\", diffusion_X, \"* dW\")\n",
    "print(\"\\n=> nu = mu - sigma^2 / 2 :\", sp.simplify(drift_X - (mu - sigma**2 / 2)) == 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "562dc912",
   "metadata": {},
   "source": [
    "### Volatility drag — why the median drifts slower than the mean\n",
    "\n",
    "Integrating the log-price SDE gives the **exact solution**\n",
    "\n",
    "$$\\boxed{\\;S_t \\;=\\; S_0\\,\\exp\\!\\Big(\\big(\\mu - \\tfrac12 \\sigma^2\\big)\\,t + \\sigma W_t\\Big).\\;}$$\n",
    "\n",
    "From this:\n",
    "\n",
    "* **Mean:** $\\;\\mathbb E[S_t] = S_0\\,e^{\\mu t}$ — drift $\\mu$.\n",
    "* **Median:** $\\;\\operatorname{Median}(S_t) = S_0\\,e^{(\\mu - \\sigma^2/2)\\,t}$ — drift $\\mu - \\sigma^2/2$.\n",
    "\n",
    "The median grows *slower* than the mean, and the gap is $\\tfrac12\\sigma^2 t$.\n",
    "This is **volatility drag**: high volatility raises the (positive-skew) mean\n",
    "but lowers the typical (median) outcome. A stock that halves then doubles has\n",
    "averaged a $41\\%$ gain yet left you exactly where you started."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "fa1a5c9f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-08-02T16:21:51.135671Z",
     "iopub.status.busy": "2026-08-02T16:21:51.135454Z",
     "iopub.status.idle": "2026-08-02T16:21:52.311349Z",
     "shell.execute_reply": "2026-08-02T16:21:52.310237Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "simulated  E[S_T]   =  110.626   formula  110.517\n",
      "simulated  median   =  105.780   formula  105.654\n",
      "vol drag over 1y    = 4.60%  (= (e^(sigma^2/2)-1) = 4.60%)\n"
     ]
    }
   ],
   "source": [
    "# Simulate GBM and check the exact-solution moments against the formulas.\n",
    "S0, T, mu, sigma = 100.0, 1.0, 0.10, 0.30\n",
    "n_paths, n_steps = 200_000, 200\n",
    "t, S = processes.geometric_brownian_motion(S0, T, n_steps, n_paths, mu, sigma, rng=RNG)\n",
    "ST = S[:, -1]\n",
    "\n",
    "mean_formula   = S0 * np.exp(mu * T)\n",
    "median_formula = S0 * np.exp((mu - 0.5*sigma**2) * T)\n",
    "\n",
    "print(f\"simulated  E[S_T]   = {ST.mean():8.3f}   formula {mean_formula:8.3f}\")\n",
    "print(f\"simulated  median   = {np.median(ST):8.3f}   formula {median_formula:8.3f}\")\n",
    "print(f\"vol drag over 1y    = {(mean_formula/median_formula - 1)*100:.2f}%  \"\n",
    "      f\"(= (e^(sigma^2/2)-1) = {(np.exp(sigma**2/2)-1)*100:.2f}%)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "1ae6ac95",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-08-02T16:21:52.313754Z",
     "iopub.status.busy": "2026-08-02T16:21:52.313497Z",
     "iopub.status.idle": "2026-08-02T16:21:53.304540Z",
     "shell.execute_reply": "2026-08-02T16:21:53.303598Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 66768 characters]",
      "text/plain": [
       "<Figure size 935x462 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8.5, 4.2))\n",
    "ax.hist(ST, bins=200, density=True, alpha=0.5, label=\"simulated $S_T$\")\n",
    "xs = np.linspace(ST.min(), ST.max(), 400)\n",
    "from scipy.stats import lognorm\n",
    "# lognorm with s=sigma, scale=S0*exp((mu-sigma^2/2)T) matches the exact law.\n",
    "pdf = lognorm.pdf(xs, s=sigma, scale=S0*np.exp((mu - 0.5*sigma**2)*T))\n",
    "ax.plot(xs, pdf, \"k-\", lw=2, label=\"exact log-normal density\")\n",
    "ax.axvline(mean_formula,   color=\"C3\", ls=\"--\", lw=1.5, label=f\"mean  $S_0 e^{{\\\\mu t}}$\")\n",
    "ax.axvline(median_formula, color=\"C2\", ls=\"--\", lw=1.5, label=f\"median $S_0 e^{{(\\\\mu-\\\\sigma^2/2)t}}$\")\n",
    "ax.set_title(\"GBM terminal prices: log-normal, mean > median (volatility drag)\")\n",
    "ax.set_xlabel(\"price $S_T$\"); ax.legend(fontsize=9); plt.savefig(_FIG / \"terminal_wealth_hist.png\", dpi=150, bbox_inches=\"tight\"); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3164e4c6",
   "metadata": {},
   "source": [
    "### Log-returns are Gaussian; prices are not\n",
    "\n",
    "A direct consequence of the exact solution is that the log-return over any\n",
    "horizon $T$ is exactly Gaussian:\n",
    "\n",
    "$$\\log S_T - \\log S_0 \\;\\sim\\; \\mathcal N\\!\\big((\\mu - \\tfrac12\\sigma^2)\\,T,\\;\\sigma^2 T\\big).$$\n",
    "\n",
    "The price $S_T$ itself is **log-normal** (and right-skewed, as the histogram\n",
    "above shows). This is why all of our barrier arithmetic in Notebooks 03–04 will\n",
    "be done in log-space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "ecc90de2",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-08-02T16:21:53.307252Z",
     "iopub.status.busy": "2026-08-02T16:21:53.306946Z",
     "iopub.status.idle": "2026-08-02T16:21:53.339665Z",
     "shell.execute_reply": "2026-08-02T16:21:53.338347Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "log-return mean=+0.0559 (formula +0.0550)\n",
      "log-return std =0.3002  (formula 0.3000)\n",
      "log-return skew=+0.001  (Gaussian => 0)\n",
      "price      skew=+0.950  (log-normal => >0)\n"
     ]
    }
   ],
   "source": [
    "# Verify: log(S_T/S_0) is Gaussian; S_T is not (right-skewed log-normal).\n",
    "logR = np.log(ST / S0)\n",
    "from scipy.stats import norm, skew\n",
    "print(f\"log-return mean={logR.mean():+.4f} (formula {(mu-0.5*sigma**2)*T:+.4f})\")\n",
    "print(f\"log-return std ={logR.std():.4f}  (formula {sigma*np.sqrt(T):.4f})\")\n",
    "print(f\"log-return skew={skew(logR):+.3f}  (Gaussian => 0)\")\n",
    "print(f\"price      skew={skew(ST):+.3f}  (log-normal => >0)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57d3d649",
   "metadata": {},
   "source": [
    "## 5. Summary\n",
    "\n",
    "* Prices must stay positive $\\Rightarrow$ the model is **multiplicative**: GBM.\n",
    "* **Itô's lemma** is the chain rule; its extra term\n",
    "  $\\tfrac12 b^2 f_{SS}$ is forced by $(dW)^2 = dt$.\n",
    "* Applied to $f = \\log$, GBM becomes arithmetic BM on the log-price with drift\n",
    "  $\\nu = \\mu - \\tfrac12\\sigma^2$ (the *volatility drag*).\n",
    "* Exact solution:\n",
    "  $S_t = S_0 \\exp\\big((\\mu - \\tfrac12\\sigma^2)t + \\sigma W_t\\big)$.\n",
    "\n",
    "For the barrier-monitoring problem we now have what we need: **the log-price is a\n",
    "Brownian motion with known drift $\\nu$ and volatility $\\sigma$**, so a $10\\%$\n",
    "price drop is the event that this Brownian motion falls below the level\n",
    "$\\ln(0.9)$. Computing the probability of that event — *continuously in time,\n",
    "not just at the endpoint* — is the **first-passage problem** of Notebook 03.\n",
    "\n",
    "---\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.14.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
