{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Waiting for the *n*-th Trade"
      ],
      "id": "c6c6fb37-6c84-441c-9310-58e950ce45f5"
    },
    {
      "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": "d02d9ce2-a17c-4234-a83a-acf6541c3919"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Introduction\n",
        "\n",
        "If trades in a bond arrive as a Poisson process with rate $\\lambda$, the\n",
        "times between trades are independent exponential random variables, and\n",
        "the waiting time $W_n$ until the $n$-th trade, the sum of $n$ of them,\n",
        "has a $\\text{Gamma}(n, \\lambda)$ distribution. This notebook takes the\n",
        "gamma model to one day of corporate bond transactions from Enhanced\n",
        "TRACE. We keep $\\alpha$ and $\\lambda$ the same for every bond, but we\n",
        "estimate the shape $\\alpha$ instead of fixing it at $n$, so the estimate\n",
        "tells us how far the data are from the Poisson benchmark. We estimate\n",
        "the model for the first through the tenth trade, under three assumptions\n",
        "about when trades happen during the day. Under the Poisson model the\n",
        "estimated shapes line up with $\\alpha = n$, and the rates are the same\n",
        "for every $n$.\n",
        "\n",
        "The density of a $\\text{Gamma}(\\alpha, \\lambda)$ random variable is\n",
        "<span id=\"eq-gamma-density\">$$\n",
        "  f(w) = \\frac{\\lambda^\\alpha}{\\Gamma(\\alpha)} w^{\\alpha - 1} e^{-\\lambda w}, \\qquad w > 0,\n",
        " \\qquad(1)$$</span> with mean $\\operatorname{E}(W) = \\alpha / \\lambda$\n",
        "and variance $\\alpha / \\lambda^2$."
      ],
      "id": "adaddfd2-350d-4f1a-ab96-ca61ac6fc8da"
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [],
      "source": [
        "import os\n",
        "import platform\n",
        "from pathlib import Path\n",
        "\n",
        "import matplotlib as mpl\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import psycopg\n",
        "from psycopg.types.numeric import FloatLoader\n",
        "from scipy.optimize import minimize\n",
        "from scipy.special import gammainc, gammaincc, gammaln\n",
        "\n",
        "# Make SVG figures byte-identical across reruns, so the frozen outputs\n",
        "# change in git only when a figure does\n",
        "os.environ.setdefault(\"SOURCE_DATE_EPOCH\", \"0\")\n",
        "mpl.rcParams[\"svg.hashsalt\"] = \"waiting-for-the-nth-trade\"\n",
        "\n",
        "DATA_DIR = Path(\"../_data\")\n",
        "\n",
        "\n",
        "def wrds_username():\n",
        "    \"\"\"The WRDS username, read from the pgpass file.\"\"\"\n",
        "    if platform.system() == \"Windows\":\n",
        "        pgpass = Path(os.environ[\"APPDATA\"]) / \"postgresql\" / \"pgpass.conf\"\n",
        "    else:\n",
        "        pgpass = Path.home() / \".pgpass\"\n",
        "    for line in pgpass.read_text().splitlines():\n",
        "        if line.startswith(\"wrds-pgdata\"):\n",
        "            return line.split(\":\")[3]\n",
        "\n",
        "\n",
        "def load_wrds(name, sql):\n",
        "    \"\"\"Run a query on WRDS and cache the result in DATA_DIR/<name>.parquet.\"\"\"\n",
        "    file = DATA_DIR / f\"{name}.parquet\"\n",
        "    if not file.exists():\n",
        "        with psycopg.connect(host=\"wrds-pgdata.wharton.upenn.edu\", port=9737,\n",
        "                             dbname=\"wrds\", user=wrds_username(),\n",
        "                             sslmode=\"require\") as conn:\n",
        "            # Return SQL numbers as floats rather than Python Decimals\n",
        "            conn.adapters.register_loader(\"numeric\", FloatLoader)\n",
        "            cur = conn.execute(sql)\n",
        "            df = pd.DataFrame(cur.fetchall(), columns=[c.name for c in cur.description])\n",
        "        DATA_DIR.mkdir(exist_ok=True)\n",
        "        df.to_parquet(file)\n",
        "    return pd.read_parquet(file)"
      ],
      "id": "2a87b94d"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Data\n",
        "\n",
        "We use all trades executed on Wednesday, October 15, 2025. The sample\n",
        "consists of corporate bonds with at least one customer trade in the\n",
        "month before that day, which keeps out the thousands of bonds\n",
        "outstanding that almost never trade. Bonds in the sample that do not\n",
        "trade on October 15 stay in, because they tell us how long the wait can\n",
        "be."
      ],
      "id": "796f92a0-958c-4c04-bbfa-5d6a925c5aa8"
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {},
      "outputs": [],
      "source": [
        "day = load_wrds(\"trace_2025-10-15\", \"\"\"\n",
        "    SELECT cusip_id, trd_exctn_tm, trc_st, msg_seq_nb, orig_msg_seq_nb,\n",
        "           entrd_vol_qt, rptd_pr, rpt_side_cd, cntra_mp_id, sub_prdct, wis_fl\n",
        "    FROM trace_enhanced.trace_enhanced\n",
        "    WHERE trd_exctn_dt = '2025-10-15'\n",
        "\"\"\")\n",
        "\n",
        "universe = load_wrds(\"trace_universe_2025-10-15\", \"\"\"\n",
        "    SELECT DISTINCT cusip_id\n",
        "    FROM trace_enhanced.trace_enhanced\n",
        "    WHERE trd_exctn_dt BETWEEN '2025-09-15' AND '2025-10-14'\n",
        "      AND cntra_mp_id = 'C'\n",
        "      AND trc_st = 'T'\n",
        "      AND sub_prdct = 'CORP'\n",
        "\"\"\")"
      ],
      "id": "16a0b381"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "TRACE records errors as separate rows that point back to the original\n",
        "trade report. We follow the filter of Dick-Nielsen (2014) for data\n",
        "reported after February 2012. A cancellation (`X`) or correction (`C`)\n",
        "repeats the original report with the same message sequence number, and\n",
        "both the original and the cancellation or correction are deleted. A\n",
        "correction also files a replacement report (`R`) with the correct terms,\n",
        "which stays in the data. A reversal (`Y`) points to the original through\n",
        "`orig_msg_seq_nb`, and both are deleted. Besides the sequence number, a\n",
        "report is matched to the original on CUSIP, volume, price, execution\n",
        "time, side and counterparty.\n",
        "\n",
        "We then keep customer trades only (`cntra_mp_id = 'C'`). A trade between\n",
        "two dealers is reported by both, and a dealer who buys from a customer\n",
        "often sells to another dealer seconds later, so counting interdealer\n",
        "trades would produce many artificially short waits. Time is measured in\n",
        "hours from 8:00 a.m. Eastern, and the day ends at 5:00 p.m., nine hours\n",
        "later."
      ],
      "id": "7b853c4a-af88-4e8c-9db4-eb64a28be5be"
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {},
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Bonds in the sample:          31,603\n",
            "Customer trades on the day:   69,085\n",
            "Bonds with no trade:          65.7%\n",
            "Bonds with  1 or more:        34.3%\n",
            "Bonds with  2 or more:        25.4%\n",
            "Bonds with  5 or more:        12.8%\n",
            "Bonds with 10 or more:        5.1%"
          ]
        }
      ],
      "source": [
        "match = [\"cusip_id\", \"entrd_vol_qt\", \"rptd_pr\", \"trd_exctn_tm\", \"rpt_side_cd\",\n",
        "         \"cntra_mp_id\", \"msg_seq_nb\"]\n",
        "cancelled = day[day.trc_st.isin([\"X\", \"C\"])]\n",
        "reversed_ = day[day.trc_st == \"Y\"].drop(columns=\"msg_seq_nb\").rename(\n",
        "    columns={\"orig_msg_seq_nb\": \"msg_seq_nb\"})\n",
        "dead = pd.concat([cancelled[match], reversed_[match]]).drop_duplicates()\n",
        "\n",
        "trades = day[day.trc_st.isin([\"T\", \"R\"])].merge(dead, on=match, how=\"left\", indicator=True)\n",
        "trades = trades[trades._merge == \"left_only\"].drop(columns=\"_merge\")\n",
        "\n",
        "trades = trades[\n",
        "    (trades.cntra_mp_id == \"C\")\n",
        "    & (trades.sub_prdct == \"CORP\")\n",
        "    & (trades.wis_fl != \"Y\")\n",
        "    & trades.cusip_id.isin(universe.cusip_id)\n",
        "]\n",
        "\n",
        "hours = pd.to_timedelta(trades.trd_exctn_tm.astype(str)).dt.total_seconds() / 3600 - 8\n",
        "trades = trades.assign(t=hours)[(hours > 0) & (hours <= 9)]\n",
        "\n",
        "counts = trades.groupby(\"cusip_id\").size().reindex(universe.cusip_id, fill_value=0)\n",
        "print(f\"Bonds in the sample:          {len(counts):,}\")\n",
        "print(f\"Customer trades on the day:   {counts.sum():,}\")\n",
        "print(f\"Bonds with no trade:          {(counts == 0).mean():.1%}\")\n",
        "for k in [1, 2, 5, 10]:\n",
        "    print(f\"Bonds with {k:>2} or more:        {(counts >= k).mean():.1%}\")"
      ],
      "id": "7b702369"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For each bond that trades at least $n$ times, the waiting time $w_{n,i}$\n",
        "is the time of its $n$-th customer trade."
      ],
      "id": "5df029fe-d395-45a5-90c6-868bb6cb528d"
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {},
      "outputs": [],
      "source": [
        "by_bond = trades.sort_values(\"t\").groupby(\"cusip_id\").t\n",
        "n_bonds = len(counts)\n",
        "\n",
        "\n",
        "def waits(n):\n",
        "    \"\"\"Time of the n-th customer trade of each bond that trades n times.\"\"\"\n",
        "    return by_bond.nth(n - 1).to_numpy()"
      ],
      "id": "a0abf38f"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Maximum likelihood with censoring\n",
        "\n",
        "We observe each bond until some time $T$. A bond whose $n$-th trade\n",
        "happens at $w_{n,i} \\leq T$ contributes the density\n",
        "(<a href=\"#eq-gamma-density\" class=\"quarto-xref\">1</a>) to the\n",
        "likelihood. A bond that has not traded $n$ times by $T$ contributes the\n",
        "probability of that event, $\\operatorname{P}(W_n > T) = 1 - F(T)$.\n",
        "Writing $n_c$ for the number of censored bonds, the log-likelihood is\n",
        "<span id=\"eq-loglik\">$$\n",
        "  \\ell(\\alpha, \\lambda) = \\sum_{w_{n,i} \\leq T} \\Big( \\alpha \\ln \\lambda - \\ln \\Gamma(\\alpha) + (\\alpha - 1) \\ln w_{n,i} - \\lambda w_{n,i} \\Big) + n_c \\ln Q(\\alpha, \\lambda T),\n",
        " \\qquad(2)$$</span> where\n",
        "$Q(\\alpha, x) = \\Gamma(\\alpha, x) / \\Gamma(\\alpha)$ is the regularized\n",
        "upper incomplete gamma function, so that\n",
        "$1 - F(T) = Q(\\alpha, \\lambda T)$. Dropping the censored bonds would\n",
        "keep only the bonds that trade often and bias $\\lambda$ upward. Because\n",
        "every bond is observed over the same window, the censoring term is the\n",
        "only adjustment needed.\n",
        "\n",
        "Later we let the rate change during the day. The code is written for\n",
        "that case from the start, using the expected number of trades by time\n",
        "$t$, $\\Lambda(t) = \\int_0^t \\lambda(s) \\, ds$, which equals $\\lambda t$\n",
        "when the rate is constant. We take $\\lambda(t)$ to be constant between\n",
        "break points, so a single break at $\\tau$ gives\n",
        "$\\Lambda(t) = \\lambda_1 \\min(t, \\tau) + \\lambda_2 \\max(t - \\tau, 0)$. We\n",
        "maximize over $\\ln \\alpha$ and the log rates to keep all parameters\n",
        "positive, and compute standard errors from a numerical Hessian of the\n",
        "log-likelihood."
      ],
      "id": "879e775c-3811-4e79-9083-609102d1bafd"
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {},
      "outputs": [],
      "source": [
        "def clock(t, rates, breaks):\n",
        "    edges = np.concatenate([[0.0], breaks, [np.inf]])\n",
        "    widths = np.diff(edges)\n",
        "    return sum(r * np.clip(t - e, 0, d) for r, e, d in zip(rates, edges, widths))\n",
        "\n",
        "\n",
        "def intensity(t, rates, breaks):\n",
        "    return np.asarray(rates)[np.searchsorted(breaks, t)]\n",
        "\n",
        "\n",
        "def neg_loglik(theta, w, T, breaks):\n",
        "    alpha, rates = np.exp(theta[0]), np.exp(theta[1:])\n",
        "    obs = w[w <= T]\n",
        "    L = clock(obs, rates, breaks)\n",
        "    observed = (alpha - 1) * np.log(L) - L - gammaln(alpha) + np.log(intensity(obs, rates, breaks))\n",
        "    censored = (n_bonds - len(obs)) * np.log(gammaincc(alpha, clock(T, rates, breaks)))\n",
        "    return -(observed.sum() + censored)\n",
        "\n",
        "\n",
        "def hessian(f, x, h=1e-4):\n",
        "    k = len(x)\n",
        "    H = np.empty((k, k))\n",
        "    E = h * np.eye(k)\n",
        "    for i in range(k):\n",
        "        for j in range(k):\n",
        "            H[i, j] = (f(x + E[i] + E[j]) - f(x + E[i] - E[j])\n",
        "                       - f(x - E[i] + E[j]) + f(x - E[i] - E[j])) / (4 * h**2)\n",
        "    return H\n",
        "\n",
        "\n",
        "def fit(w, T, breaks=(), alpha0=1.0):\n",
        "    breaks = np.asarray(breaks, dtype=float)\n",
        "    f = lambda theta: neg_loglik(theta, w, T, breaks)\n",
        "    start = np.concatenate([[np.log(alpha0)], np.zeros(len(breaks) + 1)])\n",
        "    res = minimize(f, start, method=\"Nelder-Mead\",\n",
        "                   options={\"xatol\": 1e-8, \"fatol\": 1e-8, \"maxiter\": 20000})\n",
        "    se = np.sqrt(np.diag(np.linalg.inv(hessian(f, res.x))))\n",
        "    lo, hi = np.exp(res.x[0] - 1.96 * se[0]), np.exp(res.x[0] + 1.96 * se[0])\n",
        "    return {\"alpha\": np.exp(res.x[0]), \"alpha 95% CI\": f\"[{lo:.2f}, {hi:.2f}]\",\n",
        "            \"alpha lo\": lo, \"alpha hi\": hi,\n",
        "            \"rates\": np.exp(res.x[1:]), \"breaks\": breaks, \"loglik\": -res.fun}"
      ],
      "id": "94d17041"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## The opening hours\n",
        "\n",
        "Corporate bonds trade over the counter, so there is no opening bell, but\n",
        "almost no customer trades happen before 8:00 a.m. Eastern, an hour and a\n",
        "half before the stock market opens. We first fit the model to that\n",
        "window alone, so $T = 1.5$ and a bond is censored if it has not traded\n",
        "$n$ times by 9:30. Here and below we start each search at the Poisson\n",
        "value $\\alpha = n$."
      ],
      "id": "03ff895e-a35a-466a-965d-287d5947e730"
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {},
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Bonds reaching trade  1 by 9:30: 2,243\n",
            "Bonds reaching trade  2 by 9:30:   518\n",
            "Bonds reaching trade  3 by 9:30:   148\n",
            "Bonds reaching trade  4 by 9:30:    64\n",
            "Bonds reaching trade  5 by 9:30:    31\n",
            "Bonds reaching trade  6 by 9:30:    21\n",
            "Bonds reaching trade  7 by 9:30:    12\n",
            "Bonds reaching trade  8 by 9:30:     9\n",
            "Bonds reaching trade  9 by 9:30:     6\n",
            "Bonds reaching trade 10 by 9:30:     6"
          ]
        }
      ],
      "source": [
        "ns = range(1, 11)\n",
        "\n",
        "\n",
        "def fit_all(T, breaks=(), ns=ns):\n",
        "    return {n: fit(waits(n), T, breaks, alpha0=n) for n in ns}\n",
        "\n",
        "\n",
        "for n in ns:\n",
        "    print(f\"Bonds reaching trade {n:>2} by 9:30: {(waits(n) <= 1.5).sum():>5,}\")"
      ],
      "id": "5118d4db"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Few bonds trade often enough before 9:30 to estimate the shape for large\n",
        "$n$: only 31 reach their fifth trade and 6 their tenth. We therefore fit\n",
        "the opening window for $n \\leq 4$, where at least 64 bonds reach the\n",
        "$n$-th trade."
      ],
      "id": "b44b2eb9-be8a-478d-83f9-aec4f1e01e58"
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {},
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "n = 1: alpha = 0.99 [0.95, 1.03], lambda = 0.0475\n",
            "n = 2: alpha = 1.66 [1.51, 1.82], lambda = 0.0741\n",
            "n = 3: alpha = 2.47 [2.07, 2.94], lambda = 0.1282\n",
            "n = 4: alpha = 2.64 [2.02, 3.44], lambda = 0.1112"
          ]
        }
      ],
      "source": [
        "opening = fit_all(T=1.5, ns=range(1, 5))\n",
        "for n, m in opening.items():\n",
        "    print(f\"n = {n}: alpha = {m['alpha']:.2f} {m['alpha 95% CI']}, lambda = {m['rates'][0]:.4f}\")"
      ],
      "id": "43688275"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The wait for the first trade has $\\hat\\alpha_1 = 0.99$, the exponential\n",
        "distribution the Poisson model predicts. The estimates for the second\n",
        "and third trades, 1.66 and 2.47, fall short of $n$ but not by much, and\n",
        "the fourth, 2.64, falls further short. Before the stock market opens,\n",
        "the first few bond trades arrive roughly the way the Poisson model\n",
        "describes.\n",
        "\n",
        "## The whole day\n",
        "\n",
        "Extending the window to the close, $T = 9$, and keeping one rate for the\n",
        "whole day gives a very different shape."
      ],
      "id": "7cb449ed-822a-4645-b194-df2f2751acaf"
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {},
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "n =  1: alpha = 0.85 [0.83, 0.86], lambda = 0.0353\n",
            "n =  2: alpha = 1.30 [1.26, 1.33], lambda = 0.0538\n",
            "n =  3: alpha = 1.65 [1.60, 1.70], lambda = 0.0666\n",
            "n =  4: alpha = 1.89 [1.83, 1.96], lambda = 0.0721\n",
            "n =  5: alpha = 2.09 [2.02, 2.17], lambda = 0.0749\n",
            "n =  6: alpha = 2.25 [2.16, 2.35], lambda = 0.0759\n",
            "n =  7: alpha = 2.43 [2.32, 2.55], lambda = 0.0793\n",
            "n =  8: alpha = 2.53 [2.41, 2.66], lambda = 0.0776\n",
            "n =  9: alpha = 2.66 [2.52, 2.81], lambda = 0.0789\n",
            "n = 10: alpha = 2.75 [2.59, 2.92], lambda = 0.0772"
          ]
        }
      ],
      "source": [
        "constant = fit_all(T=9)\n",
        "for n in ns:\n",
        "    m = constant[n]\n",
        "    print(f\"n = {n:>2}: alpha = {m['alpha']:.2f} {m['alpha 95% CI']}, lambda = {m['rates'][0]:.4f}\")"
      ],
      "id": "e94839fa"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<a href=\"#fig-constant\" class=\"quarto-xref\">Figure 1</a> compares the\n",
        "fitted distribution functions with the empirical ones for a few values\n",
        "of $n$. Since every bond is observed for the same nine hours, the\n",
        "empirical $F(t)$ is simply the fraction of all bonds in the sample whose\n",
        "$n$-th trade came by time $t$."
      ],
      "id": "93072f5e-8479-42cb-972b-95129f62c2a7"
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {},
      "outputs": [],
      "source": [
        "grid = np.linspace(0, 9, 400)\n",
        "hour_labels = [f\"{h}:00\" for h in range(8, 13)] + [f\"{h}:00\" for h in range(1, 6)]\n",
        "shown = [1, 2, 5, 10]\n",
        "\n",
        "\n",
        "def empirical_cdf(w):\n",
        "    return np.searchsorted(np.sort(w), grid, side=\"right\") / n_bonds\n",
        "\n",
        "\n",
        "def fitted_cdf(model):\n",
        "    return gammainc(model[\"alpha\"], clock(grid, model[\"rates\"], model[\"breaks\"]))\n",
        "\n",
        "\n",
        "def plot_fits(models):\n",
        "    fig, ax = plt.subplots(figsize=(7, 4.5))\n",
        "    for n, color in zip(shown, [\"C0\", \"C1\", \"C2\", \"C3\"]):\n",
        "        ax.plot(grid, empirical_cdf(waits(n)), color=color, label=rf\"$n = {n}$\")\n",
        "        ax.plot(grid, fitted_cdf(models[n]), color=color, linestyle=\"--\")\n",
        "    ax.set_xticks(range(10), hour_labels)\n",
        "    ax.set_xlabel(\"Time of day (Eastern)\")\n",
        "    ax.set_ylabel(r\"$P(W_n \\leq t)$\")\n",
        "    ax.set_xlim(0, 9)\n",
        "    ax.set_ylim(0)\n",
        "    ax.legend(frameon=False)\n",
        "    plt.tight_layout()\n",
        "    plt.show()"
      ],
      "id": "1eb6c56e"
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {}
        }
      ],
      "source": [
        "plot_fits(constant)"
      ],
      "id": "cell-fig-constant"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Over the whole day the estimates fall well below $n$, from 0.85 for the\n",
        "first trade to 2.75 for the tenth. In\n",
        "<a href=\"#fig-constant\" class=\"quarto-xref\">Figure 1</a> the fitted\n",
        "curves for the first and second trades rise too early and too slowly.\n",
        "The shape parameter governs how a curve starts, since\n",
        "$F(t) \\approx (\\lambda t)^\\alpha / \\Gamma(\\alpha + 1)$ for small $t$. A\n",
        "single rate has to match both the quiet opening and the busy late\n",
        "morning, and lowering $\\alpha$ is the only way the model can make the\n",
        "curve rise faster later in the day, at the cost of rising too fast at\n",
        "the start. The fifth and tenth trades come mostly after 9:30, and there\n",
        "the fitted curves follow the empirical ones closely.\n",
        "\n",
        "## A rate that changes at 9:30\n",
        "\n",
        "The empirical curves bend when the stock market opens. The simplest way\n",
        "to capture this is to let the rate jump at $\\tau = 1.5$, so that the\n",
        "bond market has one rate before 9:30 and another after. With the clock\n",
        "$\\Lambda(t)$ in place of $\\lambda t$, the waiting time is still gamma\n",
        "distributed in the sense that\n",
        "$\\operatorname{P}(W_n \\leq t) = P(\\alpha, \\Lambda(t))$, where\n",
        "$P = 1 - Q$ is the regularized lower incomplete gamma function. The\n",
        "model adds one parameter."
      ],
      "id": "0f0db2a9-cbe6-4355-b1dc-85c4e12d8bf1"
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {},
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "n =  1: alpha = 0.58 [0.57, 0.60], rates = 0.0053, 0.0184, ratio = 3.5, gain = 335.6\n",
            "n =  2: alpha = 0.79 [0.76, 0.81], rates = 0.0033, 0.0230, ratio = 7.0, gain = 440.8\n",
            "n =  3: alpha = 1.03 [0.99, 1.06], rates = 0.0036, 0.0307, ratio = 8.5, gain = 291.2\n",
            "n =  4: alpha = 1.23 [1.18, 1.28], rates = 0.0047, 0.0362, ratio = 7.8, gain = 170.8\n",
            "n =  5: alpha = 1.41 [1.34, 1.47], rates = 0.0054, 0.0399, ratio = 7.4, gain = 108.6\n",
            "n =  6: alpha = 1.57 [1.49, 1.66], rates = 0.0072, 0.0431, ratio = 5.9, gain = 66.8\n",
            "n =  7: alpha = 1.74 [1.64, 1.85], rates = 0.0092, 0.0471, ratio = 5.1, gain = 41.7\n",
            "n =  8: alpha = 1.87 [1.74, 2.01], rates = 0.0116, 0.0481, ratio = 4.1, gain = 26.0\n",
            "n =  9: alpha = 2.01 [1.85, 2.19], rates = 0.0143, 0.0509, ratio = 3.6, gain = 16.5\n",
            "n = 10: alpha = 2.15 [1.95, 2.37], rates = 0.0180, 0.0525, ratio = 2.9, gain = 10.2"
          ]
        }
      ],
      "source": [
        "jump = fit_all(T=9, breaks=[1.5])\n",
        "for n in ns:\n",
        "    m = jump[n]\n",
        "    print(f\"n = {n:>2}: alpha = {m['alpha']:.2f} {m['alpha 95% CI']}, \"\n",
        "          f\"rates = {m['rates'][0]:.4f}, {m['rates'][1]:.4f}, \"\n",
        "          f\"ratio = {m['rates'][1] / m['rates'][0]:.1f}, \"\n",
        "          f\"gain = {m['loglik'] - constant[n]['loglik']:.1f}\")"
      ],
      "id": "1e5cf9bc"
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {}
        }
      ],
      "source": [
        "plot_fits(jump)"
      ],
      "id": "cell-fig-jump"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<a href=\"#fig-jump\" class=\"quarto-xref\">Figure 2</a> shows that the jump\n",
        "at 9:30 reproduces the bend in the curves for the first and second\n",
        "trades. The first trade is fitted worst: the model puts too many first\n",
        "trades before 9:30 and too few in the late morning. For one extra\n",
        "parameter, the log-likelihood rises by 336 for the first trade, 441 for\n",
        "the second and 291 for the third. The gain then shrinks steadily, to 10\n",
        "for the tenth trade, which a single rate already fits well in\n",
        "<a href=\"#fig-constant\" class=\"quarto-xref\">Figure 1</a>.\n",
        "\n",
        "The rate after 9:30 is 3.5 times the rate before it for the first trade,\n",
        "8.5 times for the third and 2.9 times for the tenth. Under a single\n",
        "Poisson process this ratio would be the same whichever trade we wait\n",
        "for. The levels of the rates cannot be compared across models or across\n",
        "$n$, because a different $\\alpha$ needs a different clock to produce the\n",
        "same curve; only their ratio is informative.\n",
        "\n",
        "## Comparing the models\n",
        "\n",
        "<a href=\"#fig-alpha\" class=\"quarto-xref\">Figure 3</a> collects the\n",
        "estimates of the shape from the three models, and\n",
        "<a href=\"#tbl-alpha\" class=\"quarto-xref\">Table 1</a> lists them with\n",
        "their confidence intervals. The fits for different $n$ use the same\n",
        "bonds, and a bond’s tenth trade comes after its second, so the\n",
        "confidence intervals hold one at a time, not jointly."
      ],
      "id": "d7742802-c594-4c6e-b7b4-489f61175c13"
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {}
        }
      ],
      "source": [
        "models = {\n",
        "    \"Opening hours, 8:00 to 9:30\": opening,\n",
        "    \"Whole day, constant rate\": constant,\n",
        "    \"Whole day, rate jumps at 9:30\": jump,\n",
        "}\n",
        "fig, ax = plt.subplots(figsize=(7, 4.5))\n",
        "ax.plot([0.5, 10.5], [0.5, 10.5], color=\"gray\", linestyle=\"--\", label=r\"Poisson, $\\alpha = n$\")\n",
        "for (label, fits), shift, marker in zip(models.items(), [-0.15, 0, 0.15], [\"s\", \"o\", \"^\"]):\n",
        "    x = np.array(list(fits)) + shift\n",
        "    a = np.array([m[\"alpha\"] for m in fits.values()])\n",
        "    lo = np.array([m[\"alpha lo\"] for m in fits.values()])\n",
        "    hi = np.array([m[\"alpha hi\"] for m in fits.values()])\n",
        "    ax.errorbar(x, a, yerr=[a - lo, hi - a], fmt=marker, capsize=2, label=label)\n",
        "ax.set_xticks(list(ns))\n",
        "ax.set_xlim(0.5, 10.5)\n",
        "ax.set_ylim(0)\n",
        "ax.set_xlabel(\"Trade number $n$\")\n",
        "ax.set_ylabel(r\"Shape $\\alpha$\")\n",
        "ax.legend(frameon=False)\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ],
      "id": "cell-fig-alpha"
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "metadata": {},
      "outputs": [],
      "source": [
        "pd.DataFrame(\n",
        "    {label: [f\"{fits[n]['alpha']:.2f} {fits[n]['alpha 95% CI']}\" if n in fits else \"\"\n",
        "             for n in ns]\n",
        "     for label, fits in models.items()},\n",
        "    index=pd.Index(ns, name=\"n\"),\n",
        ")"
      ],
      "id": "2e6e74d5-346c-4ee0-a144-4b3d1a173b67"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Before the stock market opens, the estimates for the first three trades\n",
        "lie close to the Poisson line. Over the whole day they grow far more\n",
        "slowly than $n$: by about 0.2 per trade with a constant rate and 0.17\n",
        "with the jump at 9:30. Removing the intraday pattern moves the estimates\n",
        "further from the Poisson line, not closer, since the jump model gives a\n",
        "lower $\\hat\\alpha_n$ than the constant rate at every $n$.\n",
        "\n",
        "## What the estimates say\n",
        "\n",
        "Once the intraday pattern is removed, the remaining departure from the\n",
        "Poisson model is a shape parameter that starts below 1 and grows far\n",
        "more slowly than $n$. That is what a mixture of slow and fast bonds\n",
        "produces: a mixture of gamma distributions with different rates has more\n",
        "mass at both short and long waits than any single member of the family.\n",
        "The counts point the same way. With a common rate, the number of trades\n",
        "a bond makes during the day would be Poisson, with variance equal to its\n",
        "mean."
      ],
      "id": "05dbf2f0-78d3-412d-ad22-b200a1e50750"
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "metadata": {},
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Mean trades per bond: 2.19\n",
            "Variance:             194.31"
          ]
        }
      ],
      "source": [
        "print(f\"Mean trades per bond: {counts.mean():.2f}\")\n",
        "print(f\"Variance:             {counts.var():.2f}\")"
      ],
      "id": "360b2779"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The variance is about ninety times the mean: a few bonds trade hundreds\n",
        "of times a day while most never trade. Some of the low $\\hat\\alpha_n$\n",
        "may also come from trades clustering in time within a bond, which the\n",
        "customer-trade filter weakens but does not remove. With one day of data\n",
        "and one rate for all bonds, the two explanations cannot be told apart.\n",
        "\n",
        "A model that covers every $n$ at once gives each bond its own rate,\n",
        "drawn from a gamma distribution across bonds. The number of trades by\n",
        "time $t$ is then negative binomial rather than Poisson, so its variance\n",
        "can exceed its mean, and the waiting time to the $n$-th trade has a\n",
        "closed-form distribution with a single extra parameter, the shape of the\n",
        "distribution of rates. We leave that model for another notebook. A\n",
        "further step lets the rate depend on bond characteristics,\n",
        "$\\lambda_i(t) = \\lambda(t) \\exp(\\mathbf{x}_i' \\pmb{\\beta})$, with amount\n",
        "outstanding, age, maturity and rating in $\\mathbf{x}_i$, and checks\n",
        "whether $\\hat\\alpha_n$ moves back toward $n$ once that heterogeneity is\n",
        "accounted for. Time variation within a bond, beyond the time of day,\n",
        "leads to the autoregressive conditional duration model of Engle and\n",
        "Russell (1998).\n",
        "\n",
        "Dick-Nielsen, Jens. 2014. *How to Clean Enhanced TRACE Data*. Working\n",
        "Paper. Copenhagen Business School.\n",
        "\n",
        "Engle, Robert F., and Jeffrey R. Russell. 1998. “Autoregressive\n",
        "Conditional Duration: A New Model for Irregularly Spaced Transaction\n",
        "Data.” *Econometrica* 66 (5): 1127–62."
      ],
      "id": "77c2bac5-a930-4fc7-9b66-434b8d09a014"
    }
  ],
  "nbformat": 4,
  "nbformat_minor": 5,
  "metadata": {
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3 (ipykernel)",
      "language": "python",
      "path": "/home/lnaranjo/.pyenv/versions/3.15.2/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.15.2"
    }
  }
}