sell.py: optimal-ask model via empirical first-passage, not random walk

A driftless-random-walk closed form was tried first and backtested wrong
against 3 years of real exchange history (overstated touch probability
2-4x, confounded by discrete daily-VWAP sampling and real mean-reversion).
Replaced with puga.history.empirical_touch(): a non-parametric estimator
measuring first-passage directly from a ticker's own history, conditioned
on current state (trailing-60d z-score) since mean reversion's effect
depends on where price is now, not just move size. Queue-depth (existing
resting quantity) folds in as a deterministic shift of the hitting time,
not a max, ranking by the resulting hit_factor rather than a point
estimate of median days (which created a false cliff at the 50%-touched
threshold). Optional fixed-horizon right-censoring avoids recency bias
in "never reached" samples.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>

puga/history.py

@@ -1,6 +1,6 @@
-"""Daily exchange price/volume series from FIO's cxpc candles, shared by tools/history.py and
-tools/undervalued.py."""
-import datetime
+"""Daily exchange price/volume series from FIO's cxpc candles, shared by tools/history.py,
+tools/undervalued.py, tools/sell.py and puga/parkcash.py."""
+import datetime, math, statistics
 from . import fio
 
 
@@ -11,3 +11,126 @@ def daily_prices(tk: str, cx: str) -> dict:
         if e.get("Interval") == "DAY_ONE" and e.get("Traded"):
             out[datetime.datetime.fromtimestamp(e["DateEpochMs"] / 1000, datetime.timezone.utc).date()] = (e["Volume"] / e["Traded"], e["Traded"])
     return out
+
+
+def trading_days(tk: str, cx: str, days: int = 30) -> list[tuple[datetime.date, float]]:
+    """(date, units traded) for actual trading days in the last `days` CALENDAR days - fio.cxpc
+    only emits a row when a trade happened, so a sparse ticker needs a real calendar cutoff here,
+    not "last N rows" (that would silently reach back far more than `days` actual days)."""
+    cutoff = datetime.datetime.now(datetime.timezone.utc).date() - datetime.timedelta(days=days)
+    return sorted((d, v) for d, (_, v) in daily_prices(tk, cx).items() if d >= cutoff)
+
+
+def trade_gaps_days(dates: list[datetime.date]) -> list[int]:
+    return [(b - a).days for a, b in zip(dates, dates[1:])]
+
+
+def full_series(tk: str, cx: str) -> list[tuple[datetime.date, float]]:
+    """Full traded-price history (all available days, not windowed) ascending by date - for the
+    empirical first-passage estimator below, which needs as much real history as exists rather
+    than a 30d window, especially for distant target prices."""
+    return sorted((d, p) for d, (p, _) in daily_prices(tk, cx).items())
+
+
+def _trailing_zscores(logp: list[float], trail: int = 60, min_window: int = 10) -> list[float | None]:
+    """z-score of each day's log-price against its own trailing `trail`-day mean/stdev (None where
+    there isn't `min_window` days of history yet). A flat window (stdev 0) scores 0, not None -
+    the price IS exactly at its own flat mean, that's a real, defined state, not missing data."""
+    z = [None] * len(logp)
+    for i in range(len(logp)):
+        window = logp[max(0, i - trail + 1):i + 1]
+        if len(window) < min_window:
+            continue
+        m = statistics.mean(window)
+        sd = statistics.stdev(window) if len(window) > 1 else 0.0
+        z[i] = (logp[i] - m) / sd if sd > 1e-12 else 0.0
+    return z
+
+
+def empirical_touch(series: list[tuple[datetime.date, float]], y: float, rate: float,
+                     condition_on_state: bool = True, trail: int = 60, z_tol: float = 0.5,
+                     horizon: int | None = None) -> dict:
+    """Non-parametric first-passage stats for a target log-move `y` = ln(target_price/start_price),
+    estimated directly from how often and how fast this ticker's REAL price history has made moves
+    of that size - no random-walk/GBM assumption, so real mean-reversion shows up directly in the
+    result instead of being assumed away (see docs/decisions.md for the evidence that real prices
+    here do mean-revert).
+
+    Mean reversion means the answer depends on WHERE price is now relative to its own recent
+    trend, not just on the move size - a 20% rise is more/less likely depending on whether today
+    is already above or below the recent mean. With `condition_on_state` (default), only historical
+    start days whose z-score against a trailing `trail`-day mean are within `z_tol` of TODAY's
+    z-score (the series' last point) are used - keeping the estimate non-parametric (no fitted
+    mean-reversion speed/strength) while still making it state-aware. Falls back to using every
+    start day if there isn't enough history to compute today's z-score, or if the z_tol band is
+    empty (too little data to condition on, not a real "always/never" signal).
+
+    For each eligible historical day i, scan forward for the first day j>i where the cumulative
+    log-return reaches y (>= y if y>0, <= y if y<0). Contributes e^(-rho*dt) to the discount-factor
+    average, where dt is the elapsed calendar days and rho = ln(1+rate); a day that NEVER reaches y
+    within the observed history contributes 0 (not excluded - excluding right-censored samples
+    would bias the estimate toward only the easy cases, since the hard/slow/never cases are exactly
+    the ones that get cut off by the data running out).
+
+    CAVEAT: `series` is daily VWAP (one point per day), a discrete sample of a continuous price
+    path. Discrete monitoring can only miss touches, never invent them - an intraday spike that
+    crosses the target and reverts before the next day's VWAP is computed shows up as "not
+    touched" here even though a resting ask would genuinely have filled on it. So this estimator
+    is systematically conservative: real touch probability/hit_factor is at least this high, likely
+    somewhat higher, and the resulting "sell now" lean (tools/sell.py --tiers) is a bit stronger
+    than fully warranted. Standard discrete-monitoring barrier corrections (e.g. Broadie-Glasserman-
+    Kou) put the gap at roughly 0.58*sigma*sqrt(step) per sampling interval for a short horizon -
+    real at short horizons, shrinking as the horizon grows. Not corrected for here.
+
+    Without `horizon`, "never reached" means "not reached before the series ends" - which gives
+    start days near the end of history much less time to ever register a hit than start days near
+    the beginning, purely from running out of data, not economics. With `horizon` (days), only
+    start days with at least that much real history remaining are used, and the search for day j
+    is capped at `horizon` days past i for all of them - every sample gets the same fixed window,
+    so touched_frac/hit_factor are comparable across samples instead of recency-biased. Off by
+    default (every caller currently omits it) since it discards the most recent `horizon` days of
+    start-day samples entirely - a real tradeoff, not a strict improvement on small histories.
+
+    Returns: hit_factor (empirical E[e^-rho*tau], the discounted-value multiplier), touched_frac
+    (share of start days that ever reached y at all within their remaining history), median_days
+    (None if touched_frac < 0.5 - can't report a median that was never reached more than half the
+    time), n (sample count - note: these are heavily OVERLAPPING windows, not independent draws,
+    so n overstates the effective sample size; same caveat as the backtest this reuses).
+    """
+    n = len(series)
+    if y == 0 or n < 2:
+        return {"hit_factor": 1.0, "touched_frac": 1.0, "median_days": 0.0, "n": n}
+    dates = [d for d, _ in series]
+    logp = [math.log(p) for _, p in series if p > 0]
+    rho = math.log(1 + rate)
+
+    eligible = range(n - 1)
+    if horizon is not None:
+        eligible = [i for i in eligible if (dates[-1] - dates[i]).days >= horizon]
+    if condition_on_state:
+        z = _trailing_zscores(logp, trail)
+        z_today = z[-1]
+        if z_today is not None:
+            matched = [i for i in eligible if z[i] is not None and abs(z[i] - z_today) <= z_tol]
+            if matched:
+                eligible = matched
+
+    discount_sum, hit_days, count = 0.0, [], 0
+    for i in eligible:
+        target = logp[i] + y
+        found = None
+        for j in range(i + 1, n):
+            if horizon is not None and (dates[j] - dates[i]).days > horizon:
+                break
+            if (y > 0 and logp[j] >= target) or (y < 0 and logp[j] <= target):
+                found = j
+                break
+        count += 1
+        if found is not None:
+            dt = (dates[found] - dates[i]).days
+            discount_sum += math.exp(-rho * dt)
+            hit_days.append(dt)
+    touched_frac = len(hit_days) / count if count else 0.0
+    median_days = statistics.median(hit_days) if touched_frac >= 0.5 else None
+    return {"hit_factor": discount_sum / count if count else 0.0, "touched_frac": touched_frac,
+            "median_days": median_days, "n": count}