valuation.py

"""Valuation and relative value primitives for the Elko land analysis.

Pure functions over DataFrames and arrays, no I/O, so every rule is testable without
executing the notebook. The notebook imports this module; the module never imports the
notebook.

Bootstrap functions take an explicit generator. The notebook's `RNG` seeds every published
confidence interval, so consuming it from new code would silently move existing results.
Pass a dedicated generator instead.
"""
import re

import numpy as np
import pandas as pd

# Copied from notebook cell 7 rather than shared: prep() drives the published analysis and
# must not change. tests/test_valuation.py asserts the two patterns stay identical.
TRS_PAT = re.compile(r"SEC\s+(\d+)\s+TWP\s+(\d+)N\s+R(?:GE|NG)\s+(\d+)E", re.I)


def trs_parse(site_address):
    """Township-range-section from an assessor site address. Non-matching rows give NA."""
    m = pd.Series(site_address).fillna("").astype(str).str.extract(TRS_PAT)
    out = pd.DataFrame(index=m.index)
    for i, name in enumerate(("sec", "twp", "rge")):
        out[name] = pd.to_numeric(m[i], errors="coerce").astype("Int64")
    return out


def _clean(x):
    return pd.Series(x, dtype="float64").replace([np.inf, -np.inf], np.nan).dropna()


def boot_median_ci(x, rng, n_boot=10_000, alpha=0.05, chunk=2_000):
    """Percentile bootstrap CI for a median.

    Resampling is chunked so memory stays bounded when this runs over the 2,400-parcel
    universe rather than the 170-trade frame.
    """
    v = _clean(x).to_numpy()
    n = len(v)
    if n == 0:
        return (np.nan, np.nan)
    if n == 1:
        return (float(v[0]), float(v[0]))
    meds = np.empty(n_boot, dtype="float64")
    done = 0
    while done < n_boot:
        k = min(chunk, n_boot - done)
        meds[done:done + k] = np.median(v[rng.integers(0, n, size=(k, n))], axis=1)
        done += k
    lo, hi = np.percentile(meds, [100 * alpha / 2, 100 * (1 - alpha / 2)])
    return (float(lo), float(hi))


def median_ci_table(pairs, rng, n_boot=10_000):
    """Median with bootstrap CI for each (label, values) pair, in one table."""
    rows = []
    for label, values in pairs:
        s = _clean(values)
        lo, hi = boot_median_ci(s, rng, n_boot=n_boot)
        rows.append({"segment": label, "n": len(s),
                     "median": s.median() if len(s) else np.nan, "lo": lo, "hi": hi})
    return pd.DataFrame(rows)


def boot_ratio_ci(num, den, rng, n_boot=10_000, alpha=0.05, chunk=2_000):
    """CI for median(num) / median(den), resampling both samples independently."""
    a, b = _clean(num).to_numpy(), _clean(den).to_numpy()
    if len(a) == 0 or len(b) == 0:
        return (np.nan, np.nan)
    out = np.empty(n_boot, dtype="float64")
    done = 0
    while done < n_boot:
        k = min(chunk, n_boot - done)
        ma = np.median(a[rng.integers(0, len(a), size=(k, len(a)))], axis=1)
        mb = np.median(b[rng.integers(0, len(b), size=(k, len(b)))], axis=1)
        out[done:done + k] = ma / mb
        done += k
    lo, hi = np.percentile(out, [100 * alpha / 2, 100 * (1 - alpha / 2)])
    return (float(lo), float(hi))


def half_label(dates):
    """Calendar half labels, '2024H1' style. Unparseable dates give NA."""
    d = pd.to_datetime(pd.Series(dates), errors="coerce")
    half = pd.Series(np.where(d.dt.month <= 6, "H1", "H2"), index=d.index, dtype="string")
    return (d.dt.year.astype("Int64").astype("string") + half).mask(d.isna())


def median_index(df, period_col, value_col, rng, min_n=8, n_boot=10_000):
    """Median of `value_col` per period with bootstrap CIs.

    Periods thinner than `min_n` keep their count but carry no median: a level computed
    from a handful of trades is noise, and plotting it would imply a precision the tape
    does not support.
    """
    rows = []
    for period, g in df.groupby(period_col, observed=True, sort=True):
        v = _clean(g[value_col])
        thin = len(v) < min_n
        lo, hi = (np.nan, np.nan) if thin else boot_median_ci(v, rng, n_boot=n_boot)
        rows.append({period_col: period, "n": len(v),
                     "median": np.nan if thin else v.median(),
                     "lo": lo, "hi": hi, "suppressed": thin})
    return pd.DataFrame(rows)


def time_adjust(prices, periods, index_df, base_period, period_col="period", value_col="median"):
    """Restate prices in `base_period` money using a median index.

    Trades in suppressed periods return NaN rather than being silently carried at their
    nominal value; the caller reports how many were dropped.
    """
    idx = index_df.dropna(subset=[value_col]).set_index(period_col)[value_col]
    if base_period not in idx.index:
        raise KeyError(f"base period {base_period!r} is absent or suppressed in the index")
    prices = pd.Series(prices, dtype="float64")
    periods = pd.Series(periods)
    periods.index = prices.index
    return prices * (idx[base_period] / periods.map(idx).astype("float64"))


def prb(assessed, price, alpha=0.05):
    """IAAO coefficient of price-related bias, with a confidence interval.

    Standard on Ratio Studies: regress the proportional deviation of each assessment
    ratio from the median on the base-2 log of a value proxy that averages the
    median-adjusted sale price and the assessed value. The slope is the change in
    assessment ratio per doubling of value. Negative is regressive. IAAO calls
    -0.05 to +0.05 acceptable.

    Note the direction: this takes the assessment ratio (assessed / price), the inverse
    of the notebook's price-to-assessed `ratio`.
    """
    import statsmodels.api as sm

    d = pd.DataFrame({"assessed": pd.to_numeric(pd.Series(assessed), errors="coerce"),
                      "price": pd.to_numeric(pd.Series(price), errors="coerce")}).dropna()
    d = d[(d.price > 0) & (d.assessed > 0)]
    if len(d) < 3:
        return (np.nan, np.nan, np.nan)
    ar = d.assessed / d.price
    med = ar.median()
    y = (ar - med) / med
    x = np.log(0.5 * (d.price * med + d.assessed)) / np.log(2)
    fit = sm.OLS(y.to_numpy(), sm.add_constant(x.to_numpy())).fit()
    lo, hi = fit.conf_int(alpha=alpha)[1]
    return (float(fit.params[1]), float(lo), float(hi))


def fair_value(d, results, price_col="price", log_price_col="log_price"):
    """Attach fitted value and residual from a fitted formula model.

    The published model is a median (quantile) regression, and quantiles are equivariant
    under monotone transforms, so exp(fitted log price) is the conditional median price
    directly. No smearing retransformation is needed, and applying one would bias the
    fitted values upward.
    """
    out = d.copy()
    out["fitted_log"] = results.predict(out)
    out["resid_log"] = out[log_price_col] - out["fitted_log"]
    out["fv_price"] = np.exp(out["fitted_log"])
    out["fv_gap_pct"] = out[price_col] / out["fv_price"] - 1
    return out


def loo_mae(d, formula, q=0.5, log_price_col="log_price"):
    """Leave-one-out mean absolute error, in log points.

    In-sample fit on 162 observations with this many terms flatters itself. Refitting
    without each observation and predicting it is the cheap honest check, and it is what
    justifies leaving the model specification frozen.
    """
    import statsmodels.formula.api as smf

    errors = []
    for i in d.index:
        fit = smf.quantreg(formula, data=d.drop(index=i)).fit(q=q)
        pred = float(fit.predict(d.loc[[i]]).iloc[0])
        errors.append(abs(float(d.loc[i, log_price_col]) - pred))
    return float(np.mean(errors))


def annual_average(obs):
    """Calendar-year mean of a date-indexed series, indexed by integer year."""
    s = obs.set_index("date")["value"] if isinstance(obs, pd.DataFrame) else pd.Series(obs)
    s = s.dropna()
    out = s.groupby(s.index.year).mean()
    out.index = out.index.astype(int)
    return out


def deflate(values, years, cpi_annual, base_year):
    """Restate nominal values in `base_year` money using an annual price index."""
    if base_year not in cpi_annual.index:
        raise KeyError(f"base year {base_year} is absent from the price index")
    values = pd.Series(values, dtype="float64")
    years = pd.Series(years)
    years.index = values.index
    return values * (float(cpi_annual[base_year]) / years.map(cpi_annual).astype("float64"))


def rebase(s, base_key, base=100.0):
    """Rescale a series so it equals `base` at `base_key`."""
    s = pd.Series(s, dtype="float64")
    anchor = float(s.loc[base_key])
    if not np.isfinite(anchor) or anchor == 0:
        raise ValueError(f"cannot rebase on {base_key!r}: anchor value is {anchor}")
    return s * (base / anchor)


def ratio_percentile(s):
    """Where the last observation sits in the distribution of the whole series."""
    v = _clean(s)
    if len(v) < 2:
        return np.nan
    return float((v <= v.iloc[-1]).mean())


def carry_rate(tax_df, universe_df, id_col="Property ID", nav_col="Net Assessed Value (h)"):
    """Median annual tax bill as a share of net assessed value, and the n behind it.

    Uses the most recent bill year present, joined to the parcel universe so the rate
    reflects this band of land rather than every property type in the county. Returns
    (rate_on_nav, n), where n counts the parcels behind the median: those that joined a
    bill and carry a positive net assessed value, since a non-positive one is dropped
    rather than divided by. Fewer than three of those is not a meaningful median, so the
    rate is nan and n carries however many there were, one of 0, 1 or 2; three parcels
    joining with no usable value among them therefore returns (nan, 0), not (nan, 3).
    n is 0 too when no priced bill row survives, the one nan case reached before any join.

    The export is bill-level, so one account can carry several bill rows in a year. They
    are summed into the account's annual bill before the join: left as separate rows they
    would weight that parcel's ratio more than once in the median and make n count bill
    rows rather than the parcels it is documented to count.
    """
    t = tax_df.drop_duplicates()
    t = t[t["Account Type"] == "Real Estate"].copy()
    t["bill"] = pd.to_numeric(t["Bill Total"], errors="coerce")
    t["bill_year"] = pd.to_numeric(t["Bill Year"], errors="coerce")
    t = t[(t.bill > 0) & t.bill_year.notna()]
    if t.empty:
        return (np.nan, 0)
    t = t[t.bill_year == t.bill_year.max()]
    bills = t.groupby("Account ID", as_index=False)["bill"].sum()

    u = universe_df.drop_duplicates(id_col).copy()
    u["_nav"] = pd.to_numeric(u[nav_col], errors="coerce")
    j = u.merge(bills, left_on=id_col, right_on="Account ID")
    j = j[j._nav > 0]
    if len(j) < 3:
        return (np.nan, len(j))
    return (float((j.bill / j._nav).median()), len(j))


def comps(frame, twp, rge, asof, rungs=((1, 24), (1, 36), (None, 36)), min_n=5,
          ppa_col="ppa_adj", date_col="sale_date", twp_col="twp", rge_col="rge"):
    """Comparable sales for a subject parcel, widening the search until min_n is met.

    `rungs` are (block_radius, months) pairs tried in order; a radius of None drops the
    location filter and keeps only the time window. The first rung reaching `min_n` wins;
    if none does, the widest rung tried is returned with whatever it found.

    Sales after `asof` are never included. A comparable that had not happened yet is not
    a comparable, and letting them in would flatter any backward-looking demonstration.

    Time adjustment is the caller's job: pass a frame whose `ppa_col` is already stated in
    a common period, so this function only decides which sales qualify.
    """
    asof = pd.Timestamp(asof)
    d = frame.dropna(subset=[ppa_col, date_col, twp_col, rge_col])
    d = d[d[date_col] <= asof]

    chosen, used = d.iloc[0:0], rungs[-1]
    for radius, months in rungs:
        window = d[d[date_col] >= asof - pd.DateOffset(months=months)]
        if radius is not None:
            window = window[(window[twp_col].sub(twp).abs() <= radius)
                            & (window[rge_col].sub(rge).abs() <= radius)]
        chosen, used = window, (radius, months)
        if len(window) >= min_n:
            break

    v = _clean(chosen[ppa_col])
    return {"comps": chosen, "rung": used, "n": len(chosen),
            "med_ppa": float(v.median()) if len(v) else np.nan,
            "iqr": ((float(v.quantile(0.25)), float(v.quantile(0.75)))
                    if len(v) else (np.nan, np.nan))}


def drawdown(level):
    """Peak-to-trough decline of a level series, with the recovery date.

    Returns the drawdown series alongside the worst episode: the peak it fell from, the
    trough, and the first period that regained the peak. `years_to_recover` is NaN while
    the series is still below its prior high, which is the honest answer rather than zero.
    """
    s = _clean(level)
    if len(s) < 2:
        return {"drawdown": s, "max_drawdown": np.nan, "peak": None, "trough": None,
                "recovery": None, "years_to_recover": np.nan}
    dd = s / s.cummax() - 1
    trough = dd.idxmin()
    peak = s.loc[:trough].idxmax()
    after = s.loc[trough:]
    regained = after[after >= s.loc[peak]]
    recovery = regained.index[0] if len(regained) else None
    return {"drawdown": dd, "max_drawdown": float(dd.min()), "peak": peak,
            "trough": trough, "recovery": recovery,
            "years_to_recover": (recovery - peak) if recovery is not None else np.nan}