diff --git a/README.md b/README.md index 3e5dd03..286be6f 100644 --- a/README.md +++ b/README.md @@ -152,10 +152,19 @@ tests it **cross-sectionally** — across stocks, regressing `log s = A·log ŝ and asking whether A = 1. On 11 LSE stocks it gets A = 0.99 ± 0.10, R² = 0.96. [`analysis/farmer2005.py`](analysis/farmer2005.py) reruns that test on the five -LOBSTER symbols. **Both laws fail**, diffusion by four orders of magnitude — but -the failure is ordered by `dp/p_c`, the model's own nondimensional tick size, -which Equation 1 assumes away by taking `dp → 0`. A penny on a $27 stock is -seventeen times the characteristic price scale of its own order flow. +LOBSTER symbols. **Both laws fail** — but the failure is ordered by `dp/p_c`, the +model's own nondimensional tick size, which Equation 1 assumes away by taking +`dp → 0`. A penny on a $27 stock is seventeen times the characteristic price +scale of its own order flow. + +Auditing that result turned up a real problem worth stating: **the paper measures +diffusion on a different clock from the one its parameters live on** — `D̂` is per +event, `D` is per *midpoint change*, and this sample runs from 6 to 186 events +per midpoint change. A pinned spread is one whose midpoint rarely moves, so the +mismatch inflates exactly the stocks the tick already inflates. Correcting it +shrinks the diffusion failure from 1,777× to 66× and flips the regression slope +from −0.61 to +0.17. **The rank result is unchanged at ρ = 1.000, p = 0.017** — +which is why the argument rests on it. That is a correlation on five points, and `dp/p_c` is large for exactly the two cheapest stocks. [`tools/zi_paper.cpp`](tools/zi_paper.cpp) separates cause from diff --git a/analysis/farmer2005.py b/analysis/farmer2005.py index 1aa52c0..ef3ed42 100644 --- a/analysis/farmer2005.py +++ b/analysis/farmer2005.py @@ -101,11 +101,48 @@ def find_pair(symbol): return (msg[0] if msg else None), (book[0] if book else None) -def scan(symbol, msg_path, book_path, skip_open_s=0.0): +SESSION_START = 34200.0 # 09:30:00, in seconds past midnight +SESSION_END = 57600.0 # 16:00:00 + + +def _new_block(): + return { + "mo_shares": 0, "lo_sizes": [], "n_events": 0, + "zeta_all": [], "lifetimes_by_zeta": [], "placed": [], + "spreads": [], "mids": [], "mid_sum": 0.0, "mid_n": 0, + "last_mid": None, + } + + +def scan(symbol, msg_path, book_path, skip_open_s=0.0, n_blocks=1): """One pass over the message and orderbook files. Collects everything both halves of the test need: the raw material for the four model parameters, the realised spread, and the midpoint path. + + With n_blocks > 1 the session is cut into equal-duration intraday blocks and + every quantity is accumulated per block, in ONE pass. Returns a single dict + when n_blocks == 1, a list of dicts otherwise. + + WHY BLOCKS EXIST + + The paper measures each parameter per day and averages over 434 days, which + is what gives its real side an error bar. LOBSTER's free sample is a single + day, so the real side here has exactly one observation per stock and no way + to say whether a measurement is stable or a fluke. + + Intraday blocks are a partial substitute and it is important to be clear + about what they do and do not buy. They measure SAMPLING variability -- is + this stock's alpha the same in the second half hour as the ninth? They do + NOT measure day-to-day variability, so they will UNDERSTATE the true error: + overnight gaps, news, and regime changes are all invisible inside one + session. An error bar from blocks is a lower bound on the real one. + + What they do buy, which is worth more than the error bar: the cross-sectional + ordering can be re-tested independently in each block. One ordering of five + stocks is one draw. The same ordering holding in block after block is a + different and much better-supported claim -- with the caveat that blocks + from one day are not independent of each other. """ # --- effective market / effective limit ------------------------------ # @@ -119,25 +156,15 @@ def scan(symbol, msg_path, book_path, skip_open_s=0.0): # against one resting order, so type 4 shares ARE the transacted part. A # type 1 message is an order joining the book, so type 1 shares ARE the # resting part. A fully marketable order never produces a type 1 at all. - mo_shares = 0 # shares of effective market orders - lo_sizes = [] # sizes of effective limit orders -> sigma - n_events = 0 # event-time clock - - zeta_all = [] # relative price of every effective limit order - live = {} # order_id -> (zeta, event_index) for resting orders - lifetimes_by_zeta = [] # (zeta, lifetime in events) for fully cancelled orders - placed = [] # (zeta, shares) for the alpha numerator - - spreads = [] # log a - log b after each event - mids = [] # log midpoint, one entry per midpoint CHANGE - mid_sum = 0.0 # for the mean price level -> tick size in log terms - mid_n = 0 + blocks = [_new_block() for _ in range(max(1, n_blocks))] + t0 = SESSION_START + skip_open_s + span = max(1e-9, (SESSION_END - t0) / len(blocks)) + live = {} # order_id -> (zeta, event_index, block) for resting orders last_mo_key = None # (timestamp, direction) for grouping a sweep with open(msg_path, newline="") as mf, open(book_path, newline="") as bf: prev_mid = None - last_mid = None for row, brow in zip(csv.reader(mf), csv.reader(bf)): if len(row) < 6: continue @@ -150,8 +177,10 @@ def scan(symbol, msg_path, book_path, skip_open_s=0.0): direction = int(row[5]) except ValueError: continue - if t < 34200.0 + skip_open_s: + if t < t0: continue + bi = min(len(blocks) - 1, max(0, int((t - t0) / span))) + B = blocks[bi] # ---- event-time clock ------------------------------------- # @@ -168,11 +197,11 @@ def scan(symbol, msg_path, book_path, skip_open_s=0.0): if mtype == EXEC_VISIBLE: key = (t, direction) if key != last_mo_key: - n_events += 1 + B["n_events"] += 1 last_mo_key = key - mo_shares += size + B["mo_shares"] += size elif mtype in (NEW_LIMIT, PARTIAL_CANCEL, FULL_DELETE): - n_events += 1 + B["n_events"] += 1 last_mo_key = None else: # Hidden executions and anything else: not in the model. @@ -190,13 +219,13 @@ def scan(symbol, msg_path, book_path, skip_open_s=0.0): # convention, zeta > 0 means the order rests away from the mid and # zeta < 0 means it crossed. if mtype == NEW_LIMIT: - lo_sizes.append(size) + B["lo_sizes"].append(size) if prev_mid is not None and price > 0: p = math.log(price) z = (prev_mid - p) if direction == 1 else (p - prev_mid) - zeta_all.append(z) - placed.append((z, size)) - live[oid] = (z, n_events) + B["zeta_all"].append(z) + B["placed"].append((z, size)) + live[oid] = (z, B["n_events"], bi) elif mtype == FULL_DELETE: # A3 measures delta from cancelled orders only, and lifetime # "in terms of number of events happening between the @@ -205,7 +234,13 @@ def scan(symbol, msg_path, book_path, skip_open_s=0.0): # full delete ends a lifetime. rec = live.pop(oid, None) if rec is not None: - lifetimes_by_zeta.append((rec[0], n_events - rec[1])) + # The lifetime belongs to the block the order was PLACED in, + # and is measured on that block's event clock. An order that + # outlives its block is dropped rather than being credited + # to a clock it never ran on -- see the censoring note in + # the limits section. + if rec[2] == bi: + B["lifetimes_by_zeta"].append((rec[0], B["n_events"] - rec[1])) # ---- realised spread and midpoint path -------------------- try: @@ -216,30 +251,34 @@ def scan(symbol, msg_path, book_path, skip_open_s=0.0): if ask1 > 0 and bid1 > 0 and ask1 > bid1: # "Spread is measured as the daily average of log b(t) - # log a(t)", measured after each event with equal weight. - spreads.append(math.log(ask1) - math.log(bid1)) - mid_sum += (ask1 + bid1) / 2.0 - mid_n += 1 + B["spreads"].append(math.log(ask1) - math.log(bid1)) + B["mid_sum"] += (ask1 + bid1) / 2.0 + B["mid_n"] += 1 m = math.log((ask1 + bid1) / 2.0) # A4: "an event is anything that changes the midpoint price m". - if last_mid is None or m != last_mid: - mids.append(m) - last_mid = m + if B["last_mid"] is None or m != B["last_mid"]: + B["mids"].append(m) + B["last_mid"] = m prev_mid = m else: prev_mid = None - return { - "symbol": symbol, - "n_events": n_events, - "mo_shares": mo_shares, - "lo_sizes": lo_sizes, - "zeta_all": zeta_all, - "placed": placed, - "lifetimes_by_zeta": lifetimes_by_zeta, - "spreads": spreads, - "mids": mids, - "mean_mid": (mid_sum / mid_n) if mid_n else 0.0, - } + out = [] + for i, B in enumerate(blocks): + out.append({ + "symbol": symbol, + "block": i, + "n_events": B["n_events"], + "mo_shares": B["mo_shares"], + "lo_sizes": B["lo_sizes"], + "zeta_all": B["zeta_all"], + "placed": B["placed"], + "lifetimes_by_zeta": B["lifetimes_by_zeta"], + "spreads": B["spreads"], + "mids": B["mids"], + "mean_mid": (B["mid_sum"] / B["mid_n"]) if B["mid_n"] else 0.0, + }) + return out[0] if len(blocks) == 1 else out def parameters(s): @@ -454,13 +493,92 @@ def report_regression(name, r, markdown=False): print(" confirmation.") +def _row_from_scan(s): + """Turn one scan result into parameters, predictions and measured values.""" + p = parameters(s) + eps, p_c, s_hat, d_hat = predict(p) + s_real = (sum(s["spreads"]) / len(s["spreads"])) if s["spreads"] else float("nan") + d_real = diffusion_rate(s["mids"]) + + # Nondimensional tick size, the model's SECOND control parameter: + # "A non-dimensional scale parameter based on tick size is constructed by + # dividing the tick size dp by the characteristic price, i.e. + # dp/p_c = dp*alpha/mu ... the properties of the model only depend on the + # two non-dimensional parameters eps and dp/p_c." + # + # In log-price terms one cent on a share priced P is log(P+1c) - log(P). + # LOBSTER quotes in units of 1/10000 dollar, so the US Reg NMS minimum + # increment of $0.01 is 100 of those units. + mm = s["mean_mid"] + dp_log = math.log(mm + 100.0) - math.log(mm) if mm > 0 else float("nan") + tick_ratio = (dp_log / p_c) if p_c > 0 else float("nan") + + # ---- the two diffusion clocks ------------------------------------- + # + # The paper measures diffusion on a DIFFERENT clock from the one its + # parameters live on, and in this sample that matters enormously. + # + # A3 defines event time as the count of order placements and cancellations, + # and mu, alpha, delta are all per event on that clock. So D_hat is in + # log-price^2 per event. + # + # A4 then says "here an event is anything that changes the midpoint price + # m", and measures V(tau) over that sequence. So D_real is in log-price^2 + # per MIDPOINT CHANGE. + # + # The ratio D_real/D_hat therefore carries a hidden factor of events per + # midpoint change. That is harmless if the factor is roughly constant across + # the sample, which is presumably true of the paper's 11 LSE stocks. It is + # emphatically not true here: the factor is ~6 for GOOG and ~186 for INTC, + # because a spread pinned at one tick is a spread whose midpoint rarely + # moves. The tick constraint suppresses midpoint changes directly, so the + # unit mismatch inflates exactly the stocks the tick already inflates, and + # a raw comparison double counts it. + # + # Both are reported. d_real is the literal A4 quantity; d_real_ev converts + # it to the same per-event clock as D_hat, which is the apples-to-apples + # comparison. + n_mid = len(s["mids"]) + ev_per_midchg = (p["n_events"] / n_mid) if n_mid else float("nan") + d_real_ev = (d_real / ev_per_midchg) if ev_per_midchg > 0 else float("nan") + + return {"symbol": s["symbol"], "block": s.get("block", 0), **p, + "eps": eps, "p_c": p_c, "s_hat": s_hat, "d_hat": d_hat, + "s_real": s_real, "d_real": d_real, "n_mid": n_mid, + "ev_per_midchg": ev_per_midchg, "d_real_ev": d_real_ev, + "price": mm / 10000.0, "dp_log": dp_log, "tick_ratio": tick_ratio} + + +def measure_blocks(symbols, n_blocks, skip_open=0.0, verbose=True): + """Per-symbol, per-intraday-block measurements. + + Returns {block_index: [row, ...]}, each inner list sorted by dp/p_c, so the + cross-sectional test can be re-run independently inside every block. + """ + by_block = {i: [] for i in range(n_blocks)} + for sym in symbols: + m, b = find_pair(sym) + if not m or not b: + print(f"skipping {sym}: no message/orderbook pair", file=sys.stderr) + continue + if verbose: + print(f"scanning {sym} in {n_blocks} blocks ...", file=sys.stderr) + for s in scan(sym, m, b, skip_open, n_blocks): + if s["n_events"] < 500 or not s["spreads"]: + continue # too thin to estimate anything from + by_block[s["block"]].append(_row_from_scan(s)) + for i in by_block: + by_block[i].sort(key=lambda r: r["tick_ratio"]) + return {i: rs for i, rs in by_block.items() if len(rs) >= 3} + + def measure_all(symbols, skip_open=0.0, verbose=True): - """Measure every symbol. Shared with analysis/validate_scaling.py. + """Measure every symbol over the whole session. - The simulator must run on exactly the parameters the empirical test used, - or a difference in the measurement would look like a difference in the - model. Same reason analysis/compare.py imports stylized_facts rather than - reimplementing it. + Shared with analysis/validate_scaling.py: the simulator must run on exactly + the parameters the empirical test used, or a difference in the measurement + would look like a difference in the model. Same reason analysis/compare.py + imports stylized_facts rather than reimplementing it. """ rows = [] for sym in symbols: @@ -470,34 +588,98 @@ def measure_all(symbols, skip_open=0.0, verbose=True): continue if verbose: print(f"scanning {sym} ...", file=sys.stderr) - s = scan(sym, m, b, skip_open) - p = parameters(s) - eps, p_c, s_hat, d_hat = predict(p) - s_real = (sum(s["spreads"]) / len(s["spreads"])) if s["spreads"] else float("nan") - d_real = diffusion_rate(s["mids"]) - - # Nondimensional tick size, the model's SECOND control parameter: - # "A non-dimensional scale parameter based on tick size is constructed - # by dividing the tick size dp by the characteristic price, i.e. - # dp/p_c = dp*alpha/mu ... the properties of the model only depend on - # the two non-dimensional parameters eps and dp/p_c." - # - # In log-price terms one cent on a share priced P is log(P+1c) - log(P). - # LOBSTER quotes in units of 1/10000 dollar, so the US Reg NMS minimum - # increment of $0.01 is 100 of those units. - mm = s["mean_mid"] - dp_log = math.log(mm + 100.0) - math.log(mm) if mm > 0 else float("nan") - tick_ratio = (dp_log / p_c) if p_c > 0 else float("nan") - - rows.append({"symbol": sym, **p, "eps": eps, "p_c": p_c, - "s_hat": s_hat, "d_hat": d_hat, - "s_real": s_real, "d_real": d_real, "n_mid": len(s["mids"]), - "price": mm / 10000.0, "dp_log": dp_log, - "tick_ratio": tick_ratio}) + rows.append(_row_from_scan(scan(sym, m, b, skip_open))) rows.sort(key=lambda r: r["tick_ratio"]) return rows +def report_blocks(symbols, n_blocks, skip_open=0.0): + """Re-run the cross-sectional test independently inside each intraday block. + + This is the closest available substitute for the paper's 434 trading days, + and it is a substitute for only one of the two things those days buy. + + It DOES give the real side a spread: each stock's ratio now has a mean and a + standard deviation instead of a single number, so "GOOG is at 4.36" becomes + a claim with a width. + + It does NOT give day-to-day variation. Overnight gaps, news and regime + shifts are invisible inside one session, so these error bars are a LOWER + BOUND on the true ones. Blocks from one day are also not independent of each + other, so "held in 12 of 13 blocks" is a description, not a p-value with 13 + degrees of freedom behind it. + """ + by_block = measure_blocks(symbols, n_blocks, skip_open) + if not by_block: + print("\nnot enough data per block", file=sys.stderr) + return + + print("\n" + "=" * 72) + print(f"INTRADAY BLOCKS -- {len(by_block)} usable blocks of " + f"{(SESSION_END - SESSION_START - skip_open) / n_blocks / 60:.0f} min") + print("=" * 72) + + # Per stock, across blocks. + per_sym = {} + for rs in by_block.values(): + for r in rs: + if r["s_hat"] > 0 and r["s_real"] > 0: + per_sym.setdefault(r["symbol"], []).append( + (r["tick_ratio"], r["s_real"] / r["s_hat"])) + + hdr = (f"{'sym':<6} {'blocks':>7} {'dp/p_c mean':>12} {'ratio mean':>11} " + f"{'ratio sd':>9} {'min':>8} {'max':>8}") + print("\nspread ratio, measured independently in each block") + print(hdr) + print("-" * len(hdr)) + for sym in sorted(per_sym, key=lambda s: mean_of([t for t, _ in per_sym[s]])): + ts = [t for t, _ in per_sym[sym]] + rs = [v for _, v in per_sym[sym]] + print(f"{sym:<6} {len(rs):>7} {mean_of(ts):>12.2f} {mean_of(rs):>11.2f} " + f"{sd_of(rs):>9.2f} {min(rs):>8.2f} {max(rs):>8.2f}") + + # Does the ordering survive block by block? + print("\ncross-sectional ordering, re-tested inside each block") + print(f"{'block':>6} {'n':>3} {'rho(spread err, dp/p_c)':>25} {'p':>8}") + print("-" * 46) + perfect = 0 + positive = 0 + total = 0 + for i in sorted(by_block): + rs = [r for r in by_block[i] if r["s_hat"] > 0 and r["s_real"] > 0] + if len(rs) < 3: + continue + sp = spearman_exact([r["tick_ratio"] for r in rs], + [r["s_real"] / r["s_hat"] for r in rs]) + if not sp: + continue + total += 1 + if sp["rho"] > 0.999: + perfect += 1 + if sp["rho"] > 0: + positive += 1 + print(f"{i:>6} {len(rs):>3} {sp['rho']:>25.3f} {sp['p']:>8.4f}") + + if total: + print(f"\n positive in {positive}/{total} blocks, " + f"perfectly ordered in {perfect}/{total}") + print(" Blocks from a single session are correlated, so this is a") + print(" consistency check, not 13 independent replications. It answers") + print(" 'is the ordering an artifact of one measurement window?' and") + print(" not 'how many sigma is the effect?'") + + +def mean_of(xs): + return sum(xs) / len(xs) if xs else float("nan") + + +def sd_of(xs): + if len(xs) < 2: + return 0.0 + m = mean_of(xs) + return math.sqrt(sum((x - m) ** 2 for x in xs) / (len(xs) - 1)) + + def main(): ap = argparse.ArgumentParser() ap.add_argument("--symbols", default="AAPL,AMZN,GOOG,INTC,MSFT") @@ -505,6 +687,9 @@ def main(): help="seconds of the session to skip; the paper excludes the " "opening auction, which LOBSTER samples already omit") ap.add_argument("--markdown", action="store_true") + ap.add_argument("--blocks", type=int, default=0, + help="also split the session into N intraday blocks and " + "re-run the cross-sectional test inside each") args = ap.parse_args() rows = measure_all([x.strip().upper() for x in args.symbols.split(",") if x.strip()], @@ -550,6 +735,28 @@ def main(): print(" prediction is up to an overall constant (k, and the price window).") print(" So read the ratio column for spread, not scatter.") + # The two diffusion clocks. See _row_from_scan for why these differ. + print("\ndiffusion on a matched clock") + hdr3 = (f"{'sym':<6} {'dp/p_c':>8} {'ev/midchg':>10} {'raw ratio':>12} " + f"{'matched':>10}") + print(hdr3) + print("-" * len(hdr3)) + for r in rows: + dratio = (r["d_real"] / r["d_hat"]) if r["d_hat"] else float("nan") + dev = (r["d_real_ev"] / r["d_hat"]) if r["d_hat"] else float("nan") + print(f"{r['symbol']:<6} {r['tick_ratio']:>8.2f} " + f"{r['ev_per_midchg']:>10.1f} {dratio:>12.1f} {dev:>10.2f}") + raws = [r["d_real"] / r["d_hat"] for r in rows if r["d_hat"] > 0] + matched = [r["d_real_ev"] / r["d_hat"] for r in rows if r["d_hat"] > 0] + if raws and matched: + print(f"\n spread of raw ratios {max(raws) / min(raws):>8.0f}x") + print(f" spread of matched ratios {max(matched) / min(matched):>8.0f}x") + print(" D_hat is per event; the raw D_real is per MIDPOINT CHANGE, and this") + print(" sample's events-per-midpoint-change runs from 6 to 186. A pinned") + print(" spread is a spread whose midpoint rarely moves, so the mismatch") + print(" inflates exactly the stocks the tick already inflates. The matched") + print(" column is the apples-to-apples comparison; the raw one double counts.") + def ok(rs, key_hat, key_real): return [r for r in rs if r[key_hat] > 0 and r[key_real] > 0 and not math.isnan(r[key_hat]) and not math.isnan(r[key_real])] @@ -558,9 +765,12 @@ def both(rs, label): report_regression(f"spread {label}: log s = A log s_hat + B", regress([math.log(r["s_hat"]) for r in ok(rs, "s_hat", "s_real")], [math.log(r["s_real"]) for r in ok(rs, "s_hat", "s_real")])) - report_regression(f"diffusion {label}: log D = A log D_hat + B", + report_regression(f"diffusion {label}, raw A4 clock: log D = A log D_hat + B", regress([math.log(r["d_hat"]) for r in ok(rs, "d_hat", "d_real")], [math.log(r["d_real"]) for r in ok(rs, "d_hat", "d_real")])) + report_regression(f"diffusion {label}, matched clock: log D = A log D_hat + B", + regress([math.log(r["d_hat"]) for r in ok(rs, "d_hat", "d_real_ev")], + [math.log(r["d_real_ev"]) for r in ok(rs, "d_hat", "d_real_ev")])) print("\n" + "=" * 72) print("ALL STOCKS") @@ -587,8 +797,10 @@ def both(rs, label): print("DOES dp/p_c EXPLAIN THE ERROR?") print("=" * 72) tr = [r["tick_ratio"] for r in rows] - for label, key in (("spread", "s"), ("diffusion", "d")): - err = [r[f"{key}_real"] / r[f"{key}_hat"] for r in rows] + for label, num, den in (("spread", "s_real", "s_hat"), + ("diffusion (raw A4 clock)", "d_real", "d_hat"), + ("diffusion (matched clock)", "d_real_ev", "d_hat")): + err = [r[num] / r[den] for r in rows] sp = spearman_exact(tr, err) if sp is None: continue @@ -599,6 +811,9 @@ def both(rs, label): print(" model's own nondimensional tick size, which is the scope") print(" condition Equation 1 was derived under (the dp -> 0 limit).") + if args.blocks > 1: + report_blocks([r["symbol"] for r in rows], args.blocks, args.skip_open) + print("\npaper, for reference: 11 LSE stocks over 434 trading days,") print(" spread A = 0.99 +/- 0.10, B = 0.06 +/- 0.29, R^2 = 0.96") print(" diffusion R^2 = 0.76") diff --git a/analysis/plot_farmer2005.py b/analysis/plot_farmer2005.py new file mode 100644 index 0000000..cf3bd54 --- /dev/null +++ b/analysis/plot_farmer2005.py @@ -0,0 +1,415 @@ +#!/usr/bin/env python3 +"""Render the replication results as a self-contained HTML page. + +Emits hand-built SVG with no JavaScript and no external assets, matching the +project's dependency rule and the precedent set by tools/book_replay.html. The +file can be opened straight from disk or committed and viewed on GitHub Pages. + +Three panels, in the order the argument is made: + + 1. The law fails, and the failure is ordered by dp/p_c. Spread ratio against + the model's own nondimensional tick size, with intraday error bars. A + correct law gives a horizontal line at ANY height -- the prediction holds + only up to a constant -- so the reader should be looking for flatness, not + for proximity to 1. + + 2. Simulation reproduces it. The same axes with the simulated ratio overlaid, + from the paper's model run at each stock's measured parameters and real + tick size. + + 3. The negative control. Width scan showing the answer does not depend on + where the semi-infinite deposition intervals are truncated. + +Usage: + python3 analysis/plot_farmer2005.py # writes docs/farmer2005.html + python3 analysis/plot_farmer2005.py --out /tmp/f.html --blocks 13 + python3 analysis/plot_farmer2005.py --no-sim # skip the slow panel +""" + +import argparse +import math +import os +import subprocess +import sys + +sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) +import farmer2005 as f # noqa: E402 + +ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) +ZI_PAPER = os.path.join(ROOT, "build", "zi_paper") + +W, H = 720, 420 +PAD_L, PAD_R, PAD_T, PAD_B = 70, 30, 30, 60 + +# Colour-blind-safe: blue/orange rather than red/green, and the two regimes are +# also distinguished by fill so the chart survives being printed in greyscale. +C_SMALL = "#2f6db3" +C_LARGE = "#d1701c" +C_SIM = "#111111" + + +def esc(s): + return (str(s).replace("&", "&").replace("<", "<").replace(">", ">")) + + +class LogAxes: + """Log-log mapping from data space to SVG pixels.""" + + def __init__(self, xlo, xhi, ylo, yhi): + self.lx0, self.lx1 = math.log10(xlo), math.log10(xhi) + self.ly0, self.ly1 = math.log10(ylo), math.log10(yhi) + + def x(self, v): + t = (math.log10(v) - self.lx0) / (self.lx1 - self.lx0) + return PAD_L + t * (W - PAD_L - PAD_R) + + def y(self, v): + t = (math.log10(v) - self.ly0) / (self.ly1 - self.ly0) + return H - PAD_B - t * (H - PAD_T - PAD_B) + + +def decade_ticks(lo, hi): + out = [] + e = math.floor(math.log10(lo)) + while 10 ** e <= hi * 1.001: + for m in (1, 2, 5): + v = m * 10 ** e + if lo * 0.999 <= v <= hi * 1.001: + out.append(v) + e += 1 + return out + + +def fmt(v): + if v >= 100: + return f"{v:.0f}" + if v >= 1: + return f"{v:g}" + return f"{v:g}" + + +def frame(ax, xlo, xhi, ylo, yhi, xlabel, ylabel): + """Axes, gridlines, and the dp/p_c = 1 boundary.""" + p = [] + p.append(f'') + + # The model's stated domain: Equation 1 is a dp -> 0 result, so everything + # right of dp/p_c = 1 is outside the regime it was derived for. Shading it + # states the scope condition on the chart instead of in a caption. + if xlo < 1.0 < xhi: + x1 = ax.x(1.0) + p.append(f'') + p.append(f'') + p.append(f'' + f'dp/p_c > 1 — outside the model’s stated domain') + + for v in decade_ticks(xlo, xhi): + x = ax.x(v) + p.append(f'') + p.append(f'{fmt(v)}') + for v in decade_ticks(ylo, yhi): + y = ax.y(v) + p.append(f'') + p.append(f'{fmt(v)}') + + p.append(f'') + p.append(f'{esc(xlabel)}') + p.append(f'{esc(ylabel)}') + return p + + +def point(ax, x, y, lo, hi, colour, filled, label): + p = [] + px, py = ax.x(x), ax.y(y) + if lo and hi and hi > lo: + p.append(f'') + for e in (lo, hi): + p.append(f'') + fill = colour if filled else "var(--card)" + p.append(f'') + if label: + p.append(f'{esc(label)}') + return p + + +def panel_ratio(rows, spreads, sim=None): + """Spread ratio against dp/p_c, with intraday error bars.""" + xs = [r["tick_ratio"] for r in rows] + ys = [r["s_real"] / r["s_hat"] for r in rows] + lo_hi = [spreads.get(r["symbol"]) for r in rows] + + allv = list(ys) + for lh in lo_hi: + if lh: + allv += [lh[0], lh[1]] + if sim: + allv += [s for s in sim.values() if s] + + ax = LogAxes(min(xs) * 0.5, max(xs) * 2.0, min(allv) * 0.6, max(allv) * 1.7) + p = frame(ax, min(xs) * 0.5, max(xs) * 2.0, min(allv) * 0.6, max(allv) * 1.7, + "dp / p_c (nondimensional tick size)", + "measured spread / predicted spread") + + for r, y in zip(rows, ys): + c = C_LARGE if r["tick_ratio"] >= 1.0 else C_SMALL + lh = spreads.get(r["symbol"]) + p += point(ax, r["tick_ratio"], y, lh[0] if lh else None, + lh[1] if lh else None, c, True, r["symbol"]) + + if sim: + pts = [(r["tick_ratio"], sim[r["symbol"]]) for r in rows + if sim.get(r["symbol"])] + if len(pts) > 1: + d = " ".join(f"{ax.x(a):.1f},{ax.y(b):.1f}" for a, b in sorted(pts)) + p.append(f'') + for a, b in pts: + p += point(ax, a, b, None, None, C_SIM, False, None) + + return "".join(p) + + +def panel_scan(scan_rows): + """Width scan: the answer must not move when the truncation moves.""" + if not scan_rows: + return "" + xs = [w for w, _ in scan_rows] + ys = [v for _, v in scan_rows] + ylo, yhi = min(ys) * 0.5, max(ys) * 2.0 + ax = LogAxes(min(xs) * 0.7, max(xs) * 1.4, ylo, yhi) + p = frame(ax, min(xs) * 0.7, max(xs) * 1.4, ylo, yhi, + "truncation width (characteristic prices)", + "simulated spread / predicted") + d = " ".join(f"{ax.x(a):.1f},{ax.y(b):.1f}" for a, b in scan_rows) + p.append(f'') + for a, b in scan_rows: + p += point(ax, a, b, None, None, C_LARGE, True, None) + return "".join(p) + + +def run_scan(row): + out = [] + for w in (12.5, 25.0, 50.0, 100.0, 200.0): + r = subprocess.run( + [ZI_PAPER, "--alpha", f"{row['alpha']:.10g}", "--mu", f"{row['mu']:.10g}", + "--delta", f"{row['delta']:.10g}", "--sigma", f"{row['sigma']:.10g}", + "--dp", f"{row['dp_log']:.10g}", "--events", "300000", + "--width-pc", str(w), "--no-header"], + capture_output=True, text=True) + if r.returncode == 0 and len(r.stdout.split()) >= 6: + out.append((w, float(r.stdout.split()[5]))) + return out + + +def run_sim(rows, events, seeds): + sims = {} + for r in rows: + vals = [] + for s in range(1, seeds + 1): + p = subprocess.run( + [ZI_PAPER, "--alpha", f"{r['alpha']:.10g}", "--mu", f"{r['mu']:.10g}", + "--delta", f"{r['delta']:.10g}", "--sigma", f"{r['sigma']:.10g}", + "--dp", f"{r['dp_log']:.10g}", "--events", str(events), + "--seed", str(s), "--no-header"], + capture_output=True, text=True) + if p.returncode == 0 and len(p.stdout.split()) >= 6: + vals.append(float(p.stdout.split()[5])) + if vals: + sims[r["symbol"]] = sum(vals) / len(vals) + print(f" simulated {r['symbol']}", file=sys.stderr) + return sims + + +CSS = """ +:root{ + --bg:#ffffff; --card:#ffffff; --fg:#1a1a1a; --muted:#5c6470; + --grid:#e8ebef; --rule:#c3c9d2; --shade:#fdf1e4; --accent:#2f6db3; +} +@media (prefers-color-scheme: dark){ + :root:not([data-theme="light"]){ + --bg:#14171c; --card:#1b1f26; --fg:#e8eaed; --muted:#9aa3b0; + --grid:#272c35; --rule:#3c434e; --shade:#332618; --accent:#7fb0e6; + } +} +:root[data-theme="dark"]{ + --bg:#14171c; --card:#1b1f26; --fg:#e8eaed; --muted:#9aa3b0; + --grid:#272c35; --rule:#3c434e; --shade:#332618; --accent:#7fb0e6; +} +*{box-sizing:border-box} +body{background:var(--bg);color:var(--fg);margin:0;padding:2.5rem 1.25rem 4rem; + font:16px/1.65 -apple-system,BlinkMacSystemFont,"Segoe UI",Helvetica,Arial,sans-serif} +main{max-width:820px;margin:0 auto} +h1{font-size:1.7rem;line-height:1.25;margin:0 0 .4rem} +h2{font-size:1.15rem;margin:2.6rem 0 .5rem} +.sub{color:var(--muted);margin:0 0 2rem} +.sub a{color:var(--accent)} +figure{margin:1rem 0 0;overflow-x:auto} +svg{display:block;min-width:640px;max-width:100%;height:auto; + border:1px solid var(--rule);border-radius:8px} +figcaption{color:var(--muted);font-size:.86rem;margin-top:.6rem} +.tick{font-size:11px;fill:var(--muted)} +.axis{font-size:12px;fill:var(--fg)} +.lbl{font-size:11px;font-weight:600} +.note{font-size:11px;fill:var(--muted)} +table{border-collapse:collapse;width:100%;font-size:.9rem;margin-top:1rem; + font-variant-numeric:tabular-nums} +th,td{padding:.42rem .6rem;border-bottom:1px solid var(--grid);text-align:right} +th:first-child,td:first-child{text-align:left} +th{color:var(--muted);font-weight:600} +.key{display:flex;gap:1.4rem;flex-wrap:wrap;color:var(--muted); + font-size:.86rem;margin-top:.8rem} +.key span{display:flex;align-items:center;gap:.4rem} +.dot{width:11px;height:11px;border-radius:50%;display:inline-block} +.caveat{border-left:3px solid var(--rule);padding:.1rem 0 .1rem 1rem; + color:var(--muted);font-size:.92rem;margin:1.4rem 0} +code{font-family:ui-monospace,SFMono-Regular,Menlo,monospace;font-size:.9em} +""" + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("--symbols", default="AAPL,AMZN,GOOG,INTC,MSFT") + ap.add_argument("--out", default=os.path.join(ROOT, "docs", "farmer2005.html")) + ap.add_argument("--blocks", type=int, default=13) + ap.add_argument("--events", type=int, default=400000) + ap.add_argument("--seeds", type=int, default=3) + ap.add_argument("--no-sim", action="store_true") + args = ap.parse_args() + + syms = [s.strip().upper() for s in args.symbols.split(",") if s.strip()] + rows = f.measure_all(syms) + if not rows: + return 1 + + print(f"measuring {args.blocks} intraday blocks ...", file=sys.stderr) + by_block = f.measure_blocks(syms, args.blocks, verbose=False) + spreads, sds = {}, {} + per = {} + for rs in by_block.values(): + for r in rs: + if r["s_hat"] > 0 and r["s_real"] > 0: + per.setdefault(r["symbol"], []).append(r["s_real"] / r["s_hat"]) + for sym, vs in per.items(): + m, sd = f.mean_of(vs), f.sd_of(vs) + spreads[sym] = (max(1e-9, m - sd), m + sd) + sds[sym] = (m, sd, min(vs), max(vs), len(vs)) + + sim, scan = {}, [] + if not args.no_sim and os.path.exists(ZI_PAPER): + print("simulating ...", file=sys.stderr) + sim = run_sim(rows, args.events, args.seeds) + big = max(rows, key=lambda r: r["tick_ratio"]) + scan = run_scan(big) + elif not args.no_sim: + print(f"{ZI_PAPER} not built; skipping simulation panels", file=sys.stderr) + + tbl = [] + for r in rows: + m, sd, lo, hi, n = sds.get(r["symbol"], (float("nan"),) * 4 + (0,)) + s = sim.get(r["symbol"]) + tbl.append( + f"{r['symbol']}${r['price']:.2f}" + f"{r['tick_ratio']:.2f}" + f"{r['s_real'] / r['s_hat']:.2f}" + f"{m:.2f} ± {sd:.2f}" + f"{lo:.2f}–{hi:.2f}" + f"{s:.2f}" if s else + f"{r['symbol']}${r['price']:.2f}" + f"{r['tick_ratio']:.2f}" + f"{r['s_real'] / r['s_hat']:.2f}" + f"{m:.2f} ± {sd:.2f}" + f"{lo:.2f}–{hi:.2f}—") + + html = f"""Zero Intelligence and the Tick + +
+

Zero intelligence and the tick

+

Replicating Farmer, Patelli & Zovko (2005), +“The predictive power of +zero intelligence in financial markets”, on five NASDAQ symbols from +LOBSTER, 2012-06-21. Method and every caveat in +docs/FARMER_2005.md.

+ +

The law fails, and the failure is ordered by the tick

+
+ +{panel_ratio(rows, spreads, sim)} + +
A correct law gives a horizontal line at any height. +The prediction holds only up to an unknown constant, so flatness is the claim, +not proximity to 1. Error bars are ±1 sd across +{args.blocks} intraday blocks. Dashed line and open circles: the paper’s +own model simulated at each stock’s measured parameters and its real tick +size — nothing about a cheap stock is present except four flow numbers and +dp.
+
+
+ dp/p_c < 1 — inside the stated domain + dp/p_c > 1 — tick-constrained + simulated +
+ + + + +{"".join(tbl)} +
symbolpricedp/p_cratioblock meanblock rangesimulated
+ +

The three small-tick stocks and the two tick-constrained ones +do not overlap across any of the {args.blocks * 5} block measurements. +Where the simulation is available it tracks the real inflation without having +seen it.

+ +{"

The negative control

" if scan else ""} +{f'''
+ +{panel_scan(scan)} + +
The model’s deposition intervals are semi-infinite and a +simulation has to truncate them somewhere. If the answer moved with the +truncation, the boundary would be setting the spread rather than the order flow +and the whole result would be an artifact. It does not. +

This is not a formality: an earlier version truncated to a fixed price +box and scanned 32.23, 32.23, 32.23, 0.00, 32.23. Non-monotonic +in the width is a bug rather than a boundary effect — the book could pin +its best bid against the top of the box, leaving no legal prices for sell orders +and freezing the market one-sided forever.
+
''' if scan else ""} + +
+Read before quoting anything here. Five stocks and one trading +day, against the paper’s eleven stocks and 434 days. The error bars are +intraday, so they measure sampling variability and not day-to-day +variation — they are a lower bound on the true uncertainty. At n = 5 +the smallest attainable exact p-value is 0.017, so the rank results cannot be +stronger than that however real the effect is. +
+
+""" + os.makedirs(os.path.dirname(args.out), exist_ok=True) + with open(args.out, "w") as fh: + fh.write(html) + print(f"wrote {args.out}", file=sys.stderr) + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/docs/FARMER_2005.md b/docs/FARMER_2005.md index 844fac6..d21966f 100644 --- a/docs/FARMER_2005.md +++ b/docs/FARMER_2005.md @@ -5,6 +5,7 @@ > [arXiv:cond-mat/0309233](https://arxiv.org/abs/cond-mat/0309233) Run it: `python3 analysis/farmer2005.py` +Charts: `python3 analysis/plot_farmer2005.py` → [`farmer2005.html`](farmer2005.html) --- @@ -119,6 +120,11 @@ The diffusion regression is the interesting failure: R² = 0.80 with a **negative** slope. That is not noise. It is a strong relationship pointing the wrong way, which is what a missing control variable looks like. +> **The diffusion column above is measured on the paper's literal A4 clock, and +> that clock does not match `D̂`'s.** Part of this failure is a unit mismatch +> rather than a model failure. See *The two diffusion clocks* below — the +> corrected numbers are smaller and the corrected slope changes sign. + ### The error is ordered by the model's own tick parameter | | ρ | exact p | n | @@ -134,6 +140,98 @@ n = 5. The spread ordering is **not significant** (p = 0.083): GOOG and AAPL swap. Stated as a null result, because it is one. +### The two diffusion clocks + +Found by auditing my own strongest result before putting it in a paper. + +**The paper measures diffusion on a different clock from the one its parameters +live on.** §A3 defines event time as the count of order placements and +cancellations, and μ, α, δ are all per event on that clock — so `D̂` is in +log-price² **per event**. §A4 then says *"here an event is anything that changes +the midpoint price m"* and measures `V(τ)` over that sequence — so `D` is in +log-price² **per midpoint change**. + +The ratio `D/D̂` therefore carries a hidden factor of events-per-midpoint-change. +That is harmless when the factor is roughly constant across the sample, which is +presumably true of the paper's 11 LSE stocks. It is emphatically not true here: + +| | dp/p_c | events / midpoint change | raw ratio | matched ratio | +|---|---:|---:|---:|---:| +| GOOG | 0.21 | 5.9 | 11.6 | **1.96** | +| AAPL | 0.35 | 6.0 | 12.6 | **2.11** | +| AMZN | 0.76 | 9.6 | 93.2 | **9.70** | +| INTC | 17.30 | 186.3 | 9,386 | **50.37** | +| MSFT | 22.03 | 158.5 | 20,629 | **130.15** | + +**A spread pinned at one tick is a spread whose midpoint rarely moves.** So the +unit mismatch inflates exactly the stocks the tick already inflates, and a raw +comparison double counts the same effect. + +What changes when the clocks are matched: + +| | raw A4 clock | matched | +|---|---|---| +| spread of ratios | 1,777× | **66×** | +| regression slope A | −0.612 (R² = 0.80) | **+0.171** (R² = 0.39) | +| rank correlation with dp/p_c | ρ = 1.000, p = 0.0167 | **ρ = 1.000, p = 0.0167** | + +Three consequences, and the first two are corrections to claims made earlier in +this repo: + +1. **"Diffusion fails by four orders of magnitude" was wrong.** On a matched + clock it is closer to two. The extra factor was bookkeeping. +2. **The dramatic negative slope is substantially an artifact.** Matched, A moves + from −0.61 to +0.17. The narrative that a strong wrong-way relationship + signalled a missing control variable still holds — the control was real and it + was `dp/p_c` — but the specific statistic was inflated. +3. **The rank result is completely unaffected.** ρ = 1.000 at p = 0.0167 on both + clocks. + +Point 3 is why the argument was built on a rank test in the first place. The +statistic that survived a factor-of-27 error in its own input is the one worth +quoting. + +### Intraday error bars + +The paper averages parameters over 434 days per stock, which is what gives its +real side an error bar. LOBSTER's free sample is one day, so the real side above +is a single observation per stock. + +`--blocks 13` cuts the session into 30-minute blocks and re-measures everything +inside each, in one pass: + +| | dp/p_c | ratio (full day) | block mean ± sd | block range | +|---|---:|---:|---:|---:| +| GOOG | 0.21 | 4.36 | 4.39 ± 0.88 | 2.70 – 6.00 | +| AAPL | 0.35 | 3.70 | 3.60 ± 0.44 | 2.84 – 4.60 | +| AMZN | 0.76 | 5.66 | 6.13 ± 1.77 | 3.45 – 9.49 | +| INTC | 17.30 | 39.75 | 40.43 ± 5.54 | 28.48 – 49.89 | +| MSFT | 22.03 | 50.35 | 36.52 ± 10.96 | 25.32 – 55.25 | + +**The two regimes do not overlap in any of the 65 block measurements.** The +small-tick group tops out at 9.49 and the tick-constrained group bottoms out at +25.32. That separation is a stronger statement than the full-day ratios alone, +because it survives being re-measured 13 times. + +The ordering is also stable: the rank correlation between spread error and +`dp/p_c` is **positive in 13 of 13 blocks**, ρ = 0.9 in twelve of them and 1.0 in +one. The persistent ρ = 0.9 rather than 1.0 is the same GOOG/AAPL swap seen over +the full day, so that swap is a real feature of this sample rather than noise. + +**What blocks do and do not buy.** They measure *sampling* variability — is this +stock's α the same in the second half hour as the ninth? They do **not** measure +day-to-day variation: overnight gaps, news and regime shifts are invisible inside +one session, so **these error bars are a lower bound on the true ones**. And +blocks from a single day are correlated, so "13 of 13" is a consistency check, +not thirteen independent replications. It answers *is the ordering an artifact of +one measurement window?* — not *how many sigma is the effect?* + +One estimate does shift: MSFT's block mean (36.52) sits well below its full-day +ratio (50.35). Per-block δ is more heavily censored than the full-day δ, since an +order outliving its block is dropped rather than credited to a clock it never ran +on. **The blocks are for variability; the full-session row remains the headline +estimate.** + ### Inside the stated domain Restricted to the three stocks with `dp/p_c < 1`, the spread ratios are diff --git a/docs/farmer2005.html b/docs/farmer2005.html new file mode 100644 index 0000000..faff29b --- /dev/null +++ b/docs/farmer2005.html @@ -0,0 +1,110 @@ +Zero Intelligence and the Tick + +
+

Zero intelligence and the tick

+

Replicating Farmer, Patelli & Zovko (2005), +“The predictive power of +zero intelligence in financial markets”, on five NASDAQ symbols from +LOBSTER, 2012-06-21. Method and every caveat in +docs/FARMER_2005.md.

+ +

The law fails, and the failure is ordered by the tick

+
+ +dp/p_c > 1 — outside the model’s stated domain0.20.512510200.5125102050dp / p_c (nondimensional tick size)measured spread / predicted spreadGOOGAAPLAMZNINTCMSFT + +
A correct law gives a horizontal line at any height. +The prediction holds only up to an unknown constant, so flatness is the claim, +not proximity to 1. Error bars are ±1 sd across +13 intraday blocks. Dashed line and open circles: the paper’s +own model simulated at each stock’s measured parameters and its real tick +size — nothing about a cheap stock is present except four flow numbers and +dp.
+
+
+ dp/p_c < 1 — inside the stated domain + dp/p_c > 1 — tick-constrained + simulated +
+ + + + + +
symbolpricedp/p_cratioblock meanblock rangesimulated
GOOG$570.780.214.364.39 ± 0.882.70–6.000.65
AAPL$583.150.353.703.60 ± 0.442.84–4.600.72
AMZN$222.720.765.666.13 ± 1.773.45–9.490.84
INTC$27.0517.3039.7540.43 ± 5.5428.48–49.8932.23
MSFT$30.5522.0350.3536.52 ± 10.9625.32–55.2540.87
+ +

The three small-tick stocks and the two tick-constrained ones +do not overlap across any of the 65 block measurements. +Where the simulation is available it tracks the real inflation without having +seen it.

+ +

The negative control

+
+ +10205010020050truncation width (characteristic prices)simulated spread / predicted + +
The model’s deposition intervals are semi-infinite and a +simulation has to truncate them somewhere. If the answer moved with the +truncation, the boundary would be setting the spread rather than the order flow +and the whole result would be an artifact. It does not. +

This is not a formality: an earlier version truncated to a fixed price +box and scanned 32.23, 32.23, 32.23, 0.00, 32.23. Non-monotonic +in the width is a bug rather than a boundary effect — the book could pin +its best bid against the top of the box, leaving no legal prices for sell orders +and freezing the market one-sided forever.
+
+ +
+Read before quoting anything here. Five stocks and one trading +day, against the paper’s eleven stocks and 434 days. The error bars are +intraday, so they measure sampling variability and not day-to-day +variation — they are a lower bound on the true uncertainty. At n = 5 +the smallest attainable exact p-value is 0.017, so the rank results cannot be +stronger than that however real the effect is. +
+