The Square Root Impact Law — Empirical Evidence and Derivation

Companion notebook to the Square Root Impact Law article (Market Impact, Module 2). All computations use math and numpy only — no external market data or financial libraries required.

§1 — Setup: imports, named constants, formula definition
§2 — Impact table: bps and dollars for seven order sizes; concavity check
§3 — Linear vs. square root model comparison; ratio analysis
§4 — Dimensional analysis: participation rate sweep; doubling invariant
§5 — Validation summary: ✓/✗ checks for all quantitative claims

Press Shift+Enter or click ▶ Run on each cell. Run cells in order — §2–§5 depend on functions defined in §1.

§1 — Setup: Imports and Named Constants

We define all model parameters as named constants. The square root impact formula is:

impact_bps = Y * sigma_bps * sqrt(Q / ADV)

where Y ≈ 1.0 is the universal constant, sigma_bps is the daily volatility in basis points, Q is the order size in shares, and ADV is the average daily volume in shares.

All magic numbers are prohibited. Every parameter has a named constant and a comment explaining its provenance and units.

import math
import numpy as np

# ── Model parameters ─────────────────────────────────────────────────────────
# Universal impact constant; empirical range [0.5, 1.5]; canonical baseline = 1.0
# Source: Almgren et al. (2005), Torre & Ferrari (1997)
Y = 1.0

# Daily return volatility (decimal, not percent). 0.02 = 2% daily vol.
# Convention: close-to-close log returns, annualised vol / sqrt(252).
SIGMA_DAILY = 0.02

# Average daily volume in shares. Representative large-cap US equity.
ADV = 5_000_000  # shares

# Reference mid-price in dollars at the start of execution.
S0 = 100.0  # USD per share

# Derived: daily vol in basis points (1 bp = 0.01% = 1e-4)
SIGMA_BPS = SIGMA_DAILY * 1e4  # 200 bps for 2% daily vol

# ── Core formula ─────────────────────────────────────────────────────────────
def sqrt_impact_bps(order_shares: float, adv: float, sigma_bps: float, y: float) -> float:
    """
    Square root impact law: impact in basis points.

    Parameters
    ----------
    order_shares : float
        Order size in shares. Must be > 0.
    adv : float
        Average daily volume in shares. Must be > 0.
    sigma_bps : float
        Daily return volatility in basis points (e.g. 200 for 2%).
    y : float
        Universal impact constant (dimensionless). Canonical value: 1.0.

    Returns
    -------
    float
        Expected price impact in basis points (10^4 * delta_P / S).

    Notes
    -----
    Source: Almgren et al. (2005) "Direct estimation of equity market impact".
    Assumes: liquid equity, institutional order size, daily timescale.
    """
    if order_shares <= 0 or adv <= 0:
        raise ValueError("order_shares and adv must be strictly positive")
    participation_rate = order_shares / adv
    return y * sigma_bps * math.sqrt(participation_rate)


def sqrt_impact_dollars(order_shares: float, adv: float, sigma_bps: float,
                        y: float, price: float) -> float:
    """
    Square root impact law: total dollar cost of impact.

    Returns
    -------
    float
        Total expected impact cost in dollars (delta_P * order_shares).
    """
    impact_bps = sqrt_impact_bps(order_shares, adv, sigma_bps, y)
    delta_p_per_share = (impact_bps / 1e4) * price
    return delta_p_per_share * order_shares


print("=== Square Root Impact Law — Setup ===")
print(f"Universal constant Y         : {Y:.2f}")
print(f"Daily volatility sigma       : {SIGMA_DAILY:.2%}  ({SIGMA_BPS:.0f} bps)")
print(f"Average daily volume ADV     : {ADV:,} shares")
print(f"Reference price S0           : ${S0:.2f}")
print()

# Quick sanity check on the worked example from the article:
# Q = 100,000 shares, ADV = 5,000,000, sigma = 2%, Y = 1.0
Q_EXAMPLE = 100_000
impact_example = sqrt_impact_bps(Q_EXAMPLE, ADV, SIGMA_BPS, Y)
cost_example = sqrt_impact_dollars(Q_EXAMPLE, ADV, SIGMA_BPS, Y, S0)
print(f"Article example: Q={Q_EXAMPLE:,}")
print(f"  Participation rate Q/ADV   : {Q_EXAMPLE/ADV:.2%}")
print(f"  sqrt(Q/ADV)                : {math.sqrt(Q_EXAMPLE/ADV):.4f}")
print(f"  Impact                     : {impact_example:.2f} bps")
print(f"  Total impact cost          : ${cost_example:,.2f}")
print()
print("Functions defined: sqrt_impact_bps(), sqrt_impact_dollars()")
print("Ready for §2.")

§2 — Impact Table and Concavity Check

We compute impact for seven representative order sizes spanning 10,000 to 2,000,000 shares. For each order, we report: the participation rate, impact in bps, impact per share (dollars), and total cost. We then verify concavity: per-share impact must decrease monotonically as order size grows.

Concavity is the central property of the square root law. It is the economic justification for algorithmic execution: slicing a large order into smaller pieces reduces the per-share impact on each piece, lowering total cost versus a single block trade.

# Order sizes to analyse (shares)
Q_VALUES = [10_000, 50_000, 100_000, 250_000, 500_000, 1_000_000, 2_000_000]

# ── Print impact table ────────────────────────────────────────────────────────
col_q       = "Order (shares)"
col_part    = "Q/ADV"
col_bps     = "Impact (bps)"
col_per_sh  = "Impact/share ($)"
col_total   = "Total cost ($)"

sep = "-" * 80
print("=== Square Root Impact Law — Impact Table ===")
print(f"Parameters: ADV={ADV:,}  sigma={SIGMA_DAILY:.2%}  Y={Y:.1f}  S0=${S0:.2f}")
print(sep)
print(f"{'Order (shares)':>16}  {'Q/ADV':>8}  {'Impact (bps)':>13}  "
      f"{'Impact/share ($)':>17}  {'Total cost ($)':>15}")
print(sep)

rows = []
for q in Q_VALUES:
    part = q / ADV
    imp_bps = sqrt_impact_bps(q, ADV, SIGMA_BPS, Y)
    imp_per_share = (imp_bps / 1e4) * S0
    total_cost = imp_per_share * q
    rows.append((q, part, imp_bps, imp_per_share, total_cost))
    print(f"{q:>16,}  {part:>8.2%}  {imp_bps:>13.2f}  "
          f"{imp_per_share:>17.4f}  {total_cost:>15,.2f}")

print(sep)
print()

# ── Concavity check ───────────────────────────────────────────────────────────
print("=== Concavity Check: per-share impact must decrease with order size ===")
per_share_impacts = [r[3] for r in rows]
concave = all(per_share_impacts[i] > per_share_impacts[i+1]
              for i in range(len(per_share_impacts) - 1))

print(f"{'Order':>16}   {'Impact/share ($)':>17}   {'vs. previous':>14}")
for i, (q, _, _, imp_per_share, _) in enumerate(rows):
    if i == 0:
        change_str = "  (baseline)"
    else:
        ratio = imp_per_share / rows[i-1][3]
        change_str = f"  x{ratio:.4f}"
    print(f"{q:>16,}   {imp_per_share:>17.4f}   {change_str}")

print()
status = "PASS" if concave else "FAIL"
print(f"Concavity (per-share impact strictly decreasing): {status}")
print()

# ── Marginal impact of the last unit ─────────────────────────────────────────
print("=== Marginal Cost of Adding 10,000 Shares to Each Order ===")
MARGINAL_INCREMENT = 10_000  # shares
print(f"Marginal increment: {MARGINAL_INCREMENT:,} shares")
print(f"{'Base order':>16}   {'Marginal cost ($)':>18}   {'Avg cost ($)':>13}")
for q, _, _, imp_per_share, _ in rows:
    total_base = sqrt_impact_dollars(q, ADV, SIGMA_BPS, Y, S0)
    total_ext  = sqrt_impact_dollars(q + MARGINAL_INCREMENT, ADV, SIGMA_BPS, Y, S0)
    marginal_cost = total_ext - total_base
    print(f"{q:>16,}   {marginal_cost:>18.2f}   {imp_per_share * q:>13,.2f}")

§3 — Linear vs. Square Root Model Comparison

The linear impact model assumes impact_bps = gamma_linear * Q where gamma_linear is a constant price impact per share (in basis points per share). We calibrate the linear model to match the square root model at Q = 100,000 shares, then compare predictions across all order sizes.

The linear model is inconsistent with the empirical record for institutional orders. It systematically overestimates impact at large sizes and underestimates at small sizes relative to the square root benchmark.

# Calibration point: match linear model to sqrt model at Q_CALIB
Q_CALIB = 100_000  # shares — the article's worked example

# Sqrt impact at calibration point
impact_calib_bps = sqrt_impact_bps(Q_CALIB, ADV, SIGMA_BPS, Y)

# Linear model: impact_bps = gamma_linear * Q
# Calibrate: gamma_linear = impact_calib_bps / Q_CALIB
GAMMA_LINEAR = impact_calib_bps / Q_CALIB  # bps per share

print("=== Linear vs. Square Root Impact Model ===")
print(f"Calibration point: Q = {Q_CALIB:,} shares")
print(f"Sqrt impact at calibration: {impact_calib_bps:.4f} bps")
print(f"Linear gamma (bps/share)  : {GAMMA_LINEAR:.6f}")
print(f"  => Linear model: impact = {GAMMA_LINEAR:.6f} * Q")
print()

# Compare across Q_VALUES
sep = "-" * 90
col_q       = "Order (shares)"
col_sqrt    = "Sqrt (bps)"
col_lin     = "Linear (bps)"
col_ratio   = "Ratio Lin/Sqrt"
col_ovrest  = "Overestimate ($)"

print(sep)
print(f"{'Order (shares)':>16}  {'Sqrt (bps)':>11}  {'Linear (bps)':>13}  "
      f"{'Ratio Lin/Sqrt':>15}  {'Overestimate ($)':>17}")
print(sep)

ratios = []
for q in Q_VALUES:
    imp_sqrt = sqrt_impact_bps(q, ADV, SIGMA_BPS, Y)
    imp_lin  = GAMMA_LINEAR * q
    ratio    = imp_lin / imp_sqrt
    # Dollar overestimate: (linear - sqrt) * Q * S0/10000
    overest  = (imp_lin - imp_sqrt) / 1e4 * S0 * q
    ratios.append(ratio)
    marker = " <-- calibration" if q == Q_CALIB else ""
    print(f"{q:>16,}  {imp_sqrt:>11.2f}  {imp_lin:>13.2f}  "
          f"{ratio:>15.4f}  {overest:>17,.2f}{marker}")

print(sep)
print()

# Key insight: ratio = linear/sqrt = gamma_linear * Q / (Y * sigma_bps * sqrt(Q/ADV))
#            = (gamma_linear * sqrt(ADV) / (Y * sigma_bps)) * sqrt(Q)
# So the ratio grows as sqrt(Q): linear model diverges from sqrt at large sizes.
print("=== Ratio Analysis: Linear/Sqrt grows as sqrt(Q) ===")
print("Predicted ratio formula: ratio = (gamma * sqrt(ADV) / (Y * sigma_bps)) * sqrt(Q/ADV)^(1/2)")
print()
print(f"{'Order (shares)':>16}   {'Actual ratio':>12}   {'sqrt(Q/ADV)':>12}   {'Consistent?':>12}")
for i, q in enumerate(Q_VALUES):
    actual_ratio = ratios[i]
    sqrt_part_rate = math.sqrt(q / ADV)
    # Predicted: ratio = sqrt_impact_bps(Q_CALIB,...) / SIGMA_BPS * sqrt(Q/Q_CALIB)
    predicted_ratio = math.sqrt(q / Q_CALIB)
    consistent = abs(actual_ratio - predicted_ratio) < 1e-10
    print(f"{q:>16,}   {actual_ratio:>12.4f}   {sqrt_part_rate:>12.4f}   "
          f"{'YES' if consistent else 'NO':>12}")

print()
print("Conclusion: for Q > calibration point, the linear model OVERESTIMATES impact.")
print("At Q = 2,000,000 (40% ADV), linear overstates by ratio =", f"{ratios[-1]:.2f}x")

§4 — Dimensional Analysis and Doubling Invariant

The square root law must satisfy two invariants that follow directly from its functional form:

Doubling invariant: doubling Q multiplies impact by √2 ≈ 1.4142, not 2.
Scaling sweep: as Q/ADV increases from 0.1% to 20%, impact grows as √(Q/ADV).

We verify both numerically and demonstrate dimensional consistency of the formula.

# ── Doubling invariant ────────────────────────────────────────────────────────
# Prediction: impact(2Q) / impact(Q) = sqrt(2) for all Q.
SQRT2 = math.sqrt(2.0)
TOL = 1e-10  # numerical tolerance for floating-point equality

print("=== Doubling Invariant: impact(2Q) / impact(Q) should equal sqrt(2) ===")
print(f"Expected ratio: sqrt(2) = {SQRT2:.10f}")
print()

header_q    = "Q (shares)"
header_imp  = "impact(Q) bps"
header_imp2 = "impact(2Q) bps"
header_rat  = "Ratio"
header_chk  = "Match sqrt(2)?"

sep = "-" * 78
print(f"{'Q (shares)':>14}  {'impact(Q) bps':>14}  {'impact(2Q) bps':>15}  "
      f"{'Ratio':>10}  {'Match sqrt(2)?':>15}")
print(sep)

doubling_checks = []
for q in Q_VALUES:
    imp_q  = sqrt_impact_bps(q,     ADV, SIGMA_BPS, Y)
    imp_2q = sqrt_impact_bps(2 * q, ADV, SIGMA_BPS, Y)
    ratio  = imp_2q / imp_q
    match  = abs(ratio - SQRT2) < TOL
    doubling_checks.append(match)
    print(f"{q:>14,}  {imp_q:>14.6f}  {imp_2q:>15.6f}  "
          f"{ratio:>10.8f}  {'YES' if match else 'NO':>15}")

print(sep)
all_match = all(doubling_checks)
print(f"All ratios equal sqrt(2) within tolerance {TOL}: {'PASS' if all_match else 'FAIL'}")
print()

# ── Participation rate sweep ──────────────────────────────────────────────────
# Sweep Q/ADV from 0.1% to 20%; verify impact grows as sqrt(participation rate).
PARTICIPATION_RATES = [0.001, 0.01, 0.02, 0.05, 0.10, 0.15, 0.20]

print("=== Participation Rate Sweep ===")
print(f"{'Q/ADV':>8}  {'Q (shares)':>12}  {'Impact (bps)':>13}  "
      f"{'sqrt(Q/ADV)':>12}  {'Normalised':>11}")
print("-" * 68)

# Normalise by impact at first point to verify sqrt scaling
baseline_part = PARTICIPATION_RATES[0]
baseline_q    = int(baseline_part * ADV)
baseline_imp  = sqrt_impact_bps(baseline_q, ADV, SIGMA_BPS, Y)

for part in PARTICIPATION_RATES:
    q_shares = int(part * ADV)
    imp      = sqrt_impact_bps(q_shares, ADV, SIGMA_BPS, Y)
    sqrt_pr  = math.sqrt(part)
    # Normalised ratio: if impact ~ sqrt(Q/ADV), then ratio = sqrt(part / baseline_part)
    normalised = imp / baseline_imp
    expected   = math.sqrt(part / baseline_part)
    consistent = abs(normalised - expected) < 1e-8
    chk = "OK" if consistent else "FAIL"
    print(f"{part:>8.1%}  {q_shares:>12,}  {imp:>13.4f}  "
          f"{sqrt_pr:>12.6f}  {normalised:>11.4f}  ({chk})")

print()

# ── Dimensional consistency ───────────────────────────────────────────────────
print("=== Dimensional Consistency Check ===")
print("Formula: Delta_P (dollars) = Y * sigma * S * sqrt(Q / ADV)")
print()
print("Units check:")
print(f"  Y            : dimensionless = {Y:.1f}")
print(f"  sigma        : dimensionless (decimal return) = {SIGMA_DAILY:.4f}")
print(f"  S0           : dollars per share = {S0:.2f}")
print(f"  sqrt(Q/ADV)  : dimensionless = sqrt(shares / shares)")
print(f"  Product units: dimensionless * dimensionless * (dollars/share) * dimensionless")
print(f"               = dollars per share => Delta_P has units: $/share   [CORRECT]")
print()
print(f"  Total cost = Delta_P * Q: ($/share) * shares = dollars          [CORRECT]")
print()

# Numeric verification: direct dollar formula vs bps formula
Q_TEST = 250_000
direct_dollars = Y * SIGMA_DAILY * S0 * math.sqrt(Q_TEST / ADV) * Q_TEST
via_bps_dollars = sqrt_impact_dollars(Q_TEST, ADV, SIGMA_BPS, Y, S0)
print(f"Numeric consistency at Q={Q_TEST:,}:")
print(f"  Direct dollar formula : ${direct_dollars:,.4f}")
print(f"  Via bps formula       : ${via_bps_dollars:,.4f}")
print(f"  Difference            : ${abs(direct_dollars - via_bps_dollars):.2e}  (should be ~0)")

§5 — Validation Summary

Final cell: automated ✓/✗ checks against all quantitative claims made in the article and preceding cells. All checks must pass for the notebook to be complete.

# ── Validation suite ─────────────────────────────────────────────────────────
# Each check is a (description, boolean) pair.
# All must be True for the module to pass.

results = []

# ── Check 1: Article worked example ──────────────────────────────────────────
# Q=100k, ADV=5M, sigma=2%, Y=1.0 => impact = 28.28 bps (article states 28.3 bps)
Q_ART  = 100_000
imp_art = sqrt_impact_bps(Q_ART, ADV, SIGMA_BPS, Y)
TARGET_ARTICLE_BPS = 200.0 * math.sqrt(100_000 / 5_000_000)  # = 28.2843...
check1 = abs(imp_art - TARGET_ARTICLE_BPS) < 1e-8
results.append(("Formula matches article worked example (28.28 bps)", check1))

# ── Check 2: Concavity — per-share impact strictly decreasing ─────────────────
Q_TEST_LIST = [10_000, 50_000, 100_000, 250_000, 500_000, 1_000_000, 2_000_000]
per_share_list = [(sqrt_impact_bps(q, ADV, SIGMA_BPS, Y) / 1e4) * S0 for q in Q_TEST_LIST]
check2 = all(per_share_list[i] > per_share_list[i+1] for i in range(len(per_share_list)-1))
results.append(("Per-share impact strictly decreasing (concavity verified)", check2))

# ── Check 3: Doubling invariant — impact(2Q)/impact(Q) = sqrt(2) ─────────────
doubling_ratios = []
for q in Q_TEST_LIST:
    r = sqrt_impact_bps(2*q, ADV, SIGMA_BPS, Y) / sqrt_impact_bps(q, ADV, SIGMA_BPS, Y)
    doubling_ratios.append(abs(r - math.sqrt(2)) < 1e-10)
check3 = all(doubling_ratios)
results.append(("Doubling Q multiplies impact by sqrt(2) for all tested sizes", check3))

# ── Check 4: Dimensional consistency ─────────────────────────────────────────
# Direct dollar formula and bps-route formula must agree to machine precision.
Q_DIM = 300_000
direct = Y * SIGMA_DAILY * S0 * math.sqrt(Q_DIM / ADV) * Q_DIM
via_bps = sqrt_impact_dollars(Q_DIM, ADV, SIGMA_BPS, Y, S0)
check4 = abs(direct - via_bps) < 1e-6
results.append(("Dimensional consistency: direct $ formula = bps route formula", check4))

# ── Check 5: Linear model overestimates large orders ─────────────────────────
# At Q = 2*ADV (200% ADV), the linear model (calibrated at Q_CALIB=100k) must exceed sqrt.
Q_LARGE = 2 * ADV
gamma_lin_calib = sqrt_impact_bps(Q_CALIB, ADV, SIGMA_BPS, Y) / Q_CALIB
imp_sqrt_large  = sqrt_impact_bps(Q_LARGE, ADV, SIGMA_BPS, Y)
imp_lin_large   = gamma_lin_calib * Q_LARGE
check5 = imp_lin_large > imp_sqrt_large
results.append(("Linear model overestimates impact for large orders vs. sqrt", check5))

# ── Check 6: Impact at 100% ADV (Q = V) ──────────────────────────────────────
# Q = ADV => Q/ADV = 1 => sqrt(Q/ADV) = 1 => impact = Y * sigma_bps
Q_FULL_ADV = ADV
imp_full = sqrt_impact_bps(Q_FULL_ADV, ADV, SIGMA_BPS, Y)
check6 = abs(imp_full - Y * SIGMA_BPS) < 1e-8
results.append(("At Q=ADV: impact = Y * sigma_bps (boundary check)", check6))

# ── Check 7: Participation rate scaling ───────────────────────────────────────
# impact(4Q) / impact(Q) should equal sqrt(4) = 2.0
Q_BASE = 50_000
ratio_4x = sqrt_impact_bps(4*Q_BASE, ADV, SIGMA_BPS, Y) / sqrt_impact_bps(Q_BASE, ADV, SIGMA_BPS, Y)
check7 = abs(ratio_4x - 2.0) < 1e-10
results.append(("4x order size => 2x impact (sqrt(4) = 2.0 invariant)", check7))

# ── Check 8: Non-negativity and monotonicity ──────────────────────────────────
impacts = [sqrt_impact_bps(q, ADV, SIGMA_BPS, Y) for q in Q_TEST_LIST]
check8 = all(imp > 0 for imp in impacts) and all(
    impacts[i] < impacts[i+1] for i in range(len(impacts)-1)
)
results.append(("Total impact positive and strictly increasing with order size", check8))

# ── Print results ──────────────────────────────────────────────────────────────
print("=" * 70)
print("  VALIDATION SUMMARY — Square Root Impact Law")
print("=" * 70)
all_pass = True
for desc, passed in results:
    mark = "\u2713" if passed else "\u2717"
    print(f"  [{mark}] {desc}")
    if not passed:
        all_pass = False

print("=" * 70)
if all_pass:
    print("  RESULT: ALL CHECKS PASSED")
else:
    failed = sum(1 for _, p in results if not p)
    print(f"  RESULT: {failed} CHECK(S) FAILED — review output above")
print("=" * 70)
print()
print(f"  Impact formula: {Y:.1f} x {SIGMA_BPS:.0f} bps x sqrt(Q / {ADV:,})")
print(f"  Article example (Q={Q_ART:,}): {imp_art:.4f} bps")
print(f"  At full ADV  (Q={ADV:,}):  {sqrt_impact_bps(ADV, ADV, SIGMA_BPS, Y):.4f} bps")