""" market_intelligence.py ---------------------- Worked example by Ahmed Beddi: a weekly agricultural market-intelligence system for one irrigated cereal (paddy and milled local rice) and one horticultural crop (onion) in an irrigated river valley in the Sahel. ALL DATA ARE SYNTHETIC. The generator below produces six years of weekly prices (312 weeks, 16 Nov 2020 to 2 Nov 2026) for four markets per crop, plus an illustrative import-parity price for rice. Prices are in a generic local currency unit (LCU) per kg. Harvest calendars, costs, margins and shocks are illustrative; they are not the prices of any real country, market, agency or programme. Pipeline 1. generate() synthetic weekly prices with seasonality, trend, shocks 2. seasonal_norms() robust STL on log prices over the 5-year reference period (weeks 0-259); norm frozen and projected for the monitoring year (weeks 260-311); band = norm x exp(+/- k sigma) 3. alerts() flag when the price is outside the band on the same side for n consecutive weeks (default k = 2, n = 2) 4. spreads_margins() valley-capital spread and farm-gate-to-wholesale margin 5. transmission() Engle-Granger cointegration and error-correction model, Valley market A on Capital wholesale market 6. storage() net margin per tonne of holding from a buy month to a sell month, after losses, financing and storage fees 7. backtest(), forecast() seasonal naive vs SARIMA vs STL+ETS, rolling-origin backtest (MAPE, interval coverage), 12-week prediction bands Run: python market_intelligence.py (writes CSVs and PNGs next to it) python market_intelligence.py --ts-out PATH (also writes the demo data file) The website demo (MarketPriceMonitor.tsx) reads the series, norms and forecasts exported by this script and recomputes the band, alert rule, spreads, margins and storage calculator with the same formulas and defaults (components/demos/lib/marketPriceMonitor.ts). """ from __future__ import annotations import argparse import json import math import warnings from dataclasses import dataclass, asdict from pathlib import Path import numpy as np import pandas as pd warnings.filterwarnings("ignore") SEED = 20261106 N_WEEKS = 312 # six years of weekly observations N_REF = 260 # five-year reference period (weeks 0-259) START = pd.Timestamp("2020-11-16") # Monday; the last week starts 2026-11-02 PERIOD = 52 SERIES = [ # key, crop, market, label ("rice_farmgate", "rice", "farmgate", "Paddy, farm gate (valley)"), ("rice_valley", "rice", "valley", "Milled local rice, Valley market A"), ("rice_regional", "rice", "regional", "Milled local rice, Regional wholesale market"), ("rice_capital", "rice", "capital", "Milled local rice, Capital wholesale market"), ("rice_parity", "rice", "parity", "Imported rice, port import parity (capital)"), ("onion_farmgate", "onion", "farmgate", "Onion, farm gate (valley)"), ("onion_valley", "onion", "valley", "Onion, Valley market A"), ("onion_regional", "onion", "regional", "Onion, Regional wholesale market"), ("onion_capital", "onion", "capital", "Onion, Capital wholesale market"), ] KEYS = [s[0] for s in SERIES] @dataclass(frozen=True) class GenParams: """Generator parameters. All illustrative unless stated otherwise in the notebook.""" inflation: float = 0.04 # general price drift per year (LCU) # Import parity for rice: IPP = (FOB + freight) * FX / 1000 * (1 + duty) * (1 + margin) + port + haul fob0: float = 480.0 # world FOB price, US$/t, week 0 freight: float = 55.0 # ocean freight and insurance, US$/t fx0: float = 100.0 # LCU per US$, week 0 fx_drift: float = 0.03 # depreciation per year duty: float = 0.10 # import duty and taxes, share of CIF importer_margin: float = 0.06 # importer and wholesale margin, share port_cost: float = 2.5 # port handling, LCU/kg (week-0 prices) haul_cost: float = 1.5 # port-to-capital haulage, LCU/kg (week-0 prices) # Local rice quality_discount: float = 0.85 # local milled rice at capital, ratio to import parity (long run) ipp_passthrough: float = 0.6 # share of import-parity movements passed to local rice (log) milling_ratio: float = 0.65 # kg milled rice per kg paddy milling_cost: float = 3.0 # LCU per kg paddy (week-0 prices) # Valley-capital arbitrage (error-correction generator) rice_transfer: float = 7.0 # transport + trader cost valley -> capital, LCU/kg (week 0) onion_transfer: float = 5.0 rice_alpha: float = 0.25 # speed of adjustment of Valley market A to capital onion_alpha: float = 0.35 rice_gamma: float = 0.30 # same-week pass-through of capital price changes onion_gamma: float = 0.40 # Onion onion_base: float = 26.0 # capital wholesale, LCU/kg, week 0 level (deseasonalised) onion_farm_share: float = 0.80 # farm gate as share of Valley market A price onion_farm_cost: float = 1.0 # field-to-market cost, LCU/kg (week 0) # Shock episodes injected by the generator (used to evaluate the alert rule). # start date, weeks, what it hits, size SHOCKS = [ dict(id="S1", start="2022-08-29", weeks=10, crop="onion", where="capital", kind="supply shortfall", size=+0.35), dict(id="S2", start="2023-02-06", weeks=6, crop="rice", where="parity", kind="exchange-rate depreciation (persistent)", size=+0.12), dict(id="S3", start="2023-08-07", weeks=14, crop="rice", where="parity", kind="world price spike", size=+0.30), dict(id="S4", start="2024-08-12", weeks=8, crop="both", where="valley", kind="road cut: transfer cost up", size=+8.0), dict(id="S5", start="2025-03-10", weeks=8, crop="onion", where="capital", kind="bumper harvest glut", size=-0.25), dict(id="S6", start="2025-12-08", weeks=6, crop="rice", where="valley", kind="harvest glut, late paddy payments", size=-0.10), dict(id="S7", start="2026-09-14", weeks=8, crop="onion", where="capital", kind="storage losses, early stock-out", size=+0.30), ] def dates() -> pd.DatetimeIndex: return pd.date_range(START, periods=N_WEEKS, freq="7D") def _seasonal_profile(harvests, smooth_days: float = 12.0) -> np.ndarray: """Daily log-price seasonal profile (365 values, mean zero). Each harvest k pulls the price down by A_k at its peak-arrival day H_k; the effect decays with time constant lam_k days as stocks are drawn down. A circular Gaussian smoother (sd smooth_days) removes the jump at harvest onset. """ d = np.arange(365) prof = np.zeros(365) for H, A, lam in harvests: since = (d - H) % 365 prof -= A * np.exp(-since / lam) k = np.arange(-60, 61) w = np.exp(-0.5 * (k / smooth_days) ** 2) w /= w.sum() sm = np.array([np.sum(w * prof[(i + k) % 365]) for i in range(365)]) return sm - sm.mean() # Illustrative harvest calendars (peak-arrival day of year, amplitude in log points, decay in days) RICE_HARVESTS = [(324, 0.12, 120.0), # wet-season harvest, peak arrivals ~20 November (171, 0.06, 90.0)] # off-season (hot dry season) harvest, ~20 June ONION_HARVESTS = [(74, 0.55, 200.0)] # cold dry-season harvest, peak arrivals ~15 March def _episode(idx: pd.DatetimeIndex, start: str, weeks: int, ramp: int = 3, persistent: bool = False) -> np.ndarray: """Shape of a shock: ramps up over `ramp` weeks, holds, ramps down (unless persistent).""" t0 = idx.get_indexer([pd.Timestamp(start)])[0] x = np.zeros(len(idx)) for i in range(len(idx)): j = i - t0 if j < 0: continue if persistent: x[i] = min(1.0, (j + 1) / weeks) elif j < weeks: up = min(1.0, (j + 1) / ramp) down = min(1.0, (weeks - j) / ramp) x[i] = min(up, down) return x def _ar1(rng, n, rho, sd) -> np.ndarray: e = rng.normal(0, sd, n) x = np.zeros(n) x[0] = e[0] / math.sqrt(1 - rho ** 2) for i in range(1, n): x[i] = rho * x[i - 1] + e[i] return x def generate(g: GenParams = GenParams(), seed: int = SEED) -> pd.DataFrame: """Synthetic weekly prices, LCU/kg, rounded to 2 decimals (the canonical data).""" rng = np.random.default_rng(seed) idx = dates() n = len(idx) t = np.arange(n) doy = ((idx + pd.Timedelta(days=3)).dayofyear.values - 1) % 365 # mid-week day of year infl = np.exp(np.log(1 + g.inflation) * t / 52.0) s_rice = _seasonal_profile(RICE_HARVESTS)[doy] s_onion = _seasonal_profile(ONION_HARVESTS)[doy] harvest_rice = np.clip(-s_rice / max(1e-9, -s_rice.min()), 0, 1) # 0..1 harvest intensity harvest_onion = np.clip(-s_onion / max(1e-9, -s_onion.min()), 0, 1) sh = {s["id"]: _episode(idx, s["start"], s["weeks"], persistent=(s["id"] == "S2")) for s in SHOCKS} # Exchange rate and world price log_fx = np.log(g.fx0) + np.log(1 + g.fx_drift) * t / 52.0 + np.cumsum(rng.normal(0, 0.004, n)) log_fx += np.log(1 + 0.12) * sh["S2"] fx = np.exp(log_fx) log_fob = np.log(g.fob0) + _ar1(rng, n, 0.97, 0.012) + np.log(1 + 0.30) * sh["S3"] fob = np.exp(log_fob) parity = ((fob + g.freight) * fx / 1000.0 * (1 + g.duty) * (1 + g.importer_margin) + (g.port_cost + g.haul_cost) * infl) # Local milled rice, capital: partial pass-through of import parity, quality discount, seasonality ipp_bar = parity[0] * infl log_cap_r = (np.log(g.quality_discount) + g.ipp_passthrough * np.log(parity) + (1 - g.ipp_passthrough) * np.log(ipp_bar) + 0.6 * s_rice + _ar1(rng, n, 0.8, 0.010)) cap_r = np.exp(log_cap_r) # Onion, capital: seasonal, persistent noise, supply shocks log_cap_o = (np.log(g.onion_base) + np.log(infl) + 0.85 * s_onion + _ar1(rng, n, 0.85, 0.028) + np.log(1 + 0.35) * sh["S1"] + np.log(1 - 0.25) * sh["S5"] + np.log(1 + 0.30) * sh["S7"]) cap_o = np.exp(log_cap_o) cap_r = cap_r * (1 + 0.04 * sh["S4"]) # road cut: capital a little dearer cap_o = cap_o * (1 + 0.06 * sh["S4"]) def valley(cap, transfer0, alpha, gamma, s_extra, harvest, road, noise_sd, local_shock): tau = transfer0 * infl * (1 + 0.25 * harvest) + road v = np.zeros(n) target0 = cap[0] * np.exp(s_extra[0]) - tau[0] v[0] = target0 eps = rng.normal(0, noise_sd, n) for i in range(1, n): target = cap[i - 1] * np.exp(s_extra[i]) * (1 + local_shock[i]) - tau[i] v[i] = v[i - 1] + alpha * (target - v[i - 1]) + gamma * (cap[i] - cap[i - 1]) + eps[i] * infl[i] return v road = 8.0 * sh["S4"] * infl val_r = valley(cap_r, g.rice_transfer, g.rice_alpha, g.rice_gamma, 0.4 * s_rice, harvest_rice, road, 0.45, -0.10 * sh["S6"]) val_o = valley(cap_o, g.onion_transfer, g.onion_alpha, g.onion_gamma, 0.15 * s_onion, harvest_onion, road, 0.45, np.zeros(n)) reg_r = 0.55 * cap_r + 0.45 * val_r + _ar1(rng, n, 0.5, 0.35) * infl reg_o = 0.55 * cap_o + 0.45 * val_o + _ar1(rng, n, 0.5, 0.40) * infl # Farm gate: paddy = milling ratio x valley milled price - milling cost, with a harvest discount fg_r = (g.milling_ratio * val_r * (1 - 0.05 * harvest_rice) - g.milling_cost * infl + _ar1(rng, n, 0.6, 0.25) * infl) fg_o = (g.onion_farm_share * val_o * (1 - 0.08 * harvest_onion) - g.onion_farm_cost * infl + _ar1(rng, n, 0.6, 0.30) * infl) df = pd.DataFrame({ "week": t, "date": idx.strftime("%Y-%m-%d"), "period": np.where(t < N_REF, "reference", "monitoring"), "rice_farmgate": fg_r, "rice_valley": val_r, "rice_regional": reg_r, "rice_capital": cap_r, "rice_parity": parity, "onion_farmgate": fg_o, "onion_valley": val_o, "onion_regional": reg_o, "onion_capital": cap_o, "fx_lcu_per_usd": fx, "world_fob_usd_t": fob, }) for k in KEYS: df[k] = df[k].round(2) df["fx_lcu_per_usd"] = df["fx_lcu_per_usd"].round(3) df["world_fob_usd_t"] = df["world_fob_usd_t"].round(1) return df # ---------------------------------------------------------------- seasonal norms and alerts @dataclass(frozen=True) class AlertRule: k: float = 2.0 # band half-width, in robust SDs of log(price / norm) over the reference period n: int = 2 # consecutive weeks outside the band (same side) to raise a flag window: int = 13 # weeks of deseasonalised prices in the level estimate gap: int = 4 # weeks between the end of the level window and the current week def seasonal_index(df: pd.DataFrame) -> dict: """Robust STL (period 52, seasonal window 7) on log prices of the 5-year reference period. Returns per series: the STL components, and a fixed 52-week seasonal index s_p (p = week mod 52), the mean STL seasonal at each cycle position, centred on zero and rounded to 5 decimals (the exported value used everywhere downstream). """ from statsmodels.tsa.seasonal import STL out = {} pos = np.arange(N_REF) % PERIOD for k in KEYS: y = np.log(df[k].values[:N_REF]) res = STL(y, period=PERIOD, seasonal=7, robust=True).fit() s = np.array([res.seasonal[pos == p].mean() for p in range(PERIOD)]) s = np.round(s - s.mean(), 5) out[k] = dict(index=s, trend=res.trend, seasonal=res.seasonal, remainder=res.resid, strength=float(max(0.0, 1 - np.var(res.resid) / np.var(res.resid + res.seasonal)))) return out def seasonal_norm(price: np.ndarray, s: np.ndarray, rule: AlertRule = AlertRule()) -> np.ndarray: """Causal seasonal norm: norm_t = exp(s_{t mod 52}) x L_t, where L_t is the median of the deseasonalised price p_j / exp(s_{j mod 52}) over weeks t-gap-window+1 .. t-gap. Undefined (NaN) for the first gap+window-1 weeks.""" n = len(price) sf = np.exp(s[np.arange(n) % PERIOD]) des = price / sf norm = np.full(n, np.nan) for t in range(rule.gap + rule.window - 1, n): lo, hi = t - rule.gap - rule.window + 1, t - rule.gap + 1 norm[t] = sf[t] * np.median(des[lo:hi]) return norm def band_sigma(price: np.ndarray, norm: np.ndarray) -> float: """Robust SD (1.4826 x MAD) of log(price / norm) over the reference period.""" r = np.log(price[:N_REF] / norm[:N_REF]) r = r[~np.isnan(r)] return float(1.4826 * np.median(np.abs(r - np.median(r)))) def band(norm: np.ndarray, sigma: float, k: float): return norm * np.exp(-k * sigma), norm * np.exp(k * sigma) def alerts(price: np.ndarray, norm: np.ndarray, sigma: float, rule: AlertRule = AlertRule()) -> np.ndarray: """+1 high flag, -1 low flag, 0 none. A week is flagged when it and the previous n-1 weeks are all outside the band on the same side (causal: uses only weeks <= t).""" lo, hi = band(norm, sigma, rule.k) side = np.where(price > hi, 1, np.where(price < lo, -1, 0)) # NaN norm compares False -> 0 flag = np.zeros(len(price), dtype=int) run, last = 0, 0 for i, s in enumerate(side): run = run + 1 if (s != 0 and s == last) else (1 if s != 0 else 0) last = s if s != 0 and run >= rule.n: flag[i] = s return flag def monthly_norm(price: np.ndarray, s: np.ndarray): """Seasonal-norm price by calendar month at the current level, for the storage calculator. level: mean log deseasonalised price over the last 52 weeks, moved forward 25.5 weeks to the latest week with the drift d (average weekly change in the log deseasonalised level between the first and the last 52 weeks). Month m uses the mean of s over cycle positions whose mid-week date falls in month m (dates of the reference period).""" n = len(price) des = np.log(price) - s[np.arange(n) % PERIOD] d = (des[-52:].mean() - des[:52].mean()) / (n - 52) level = des[-52:].mean() + d * 25.5 months = (dates() + pd.Timedelta(days=3)).month.values[:N_REF] pos = np.arange(N_REF) % PERIOD sm_ = np.array([s[pos[months == m]].mean() for m in range(1, 13)]) return np.round(np.exp(level + sm_), 2), round(float(d), 6) def seasonal_norms(df: pd.DataFrame, rule: AlertRule = AlertRule()) -> dict: si = seasonal_index(df) out = {} for k in KEYS: p = df[k].values norm = seasonal_norm(p, si[k]["index"], rule) mn, d = monthly_norm(p, si[k]["index"]) out[k] = dict(si[k], norm=norm, sigma=round(band_sigma(p, norm), 4), monthly_norm=mn, drift_week=d) return out # The series each injected shock hits first, and its direction (for the detection table) SHOCK_TARGETS = {"S1": ["onion_capital"], "S2": ["rice_parity"], "S3": ["rice_parity", "rice_capital"], "S4": ["rice_valley", "onion_valley"], "S5": ["onion_capital"], "S6": ["rice_valley"], "S7": ["onion_capital"]} def evaluate_alerts(df: pd.DataFrame, norms: dict, rule: AlertRule = AlertRule(), buffer: int = 8): """(1) Detection: for each shock and each series it hits first, the first week (counted from the shock start) with a flag in the shock's direction, within the shock window plus buffer. (2) Flags outside every shock window of that crop (plus buffer), by period.""" idx = dates() flags = {k: alerts(df[k].values, norms[k]["norm"], norms[k]["sigma"], rule) for k in KEYS} det = [] for s in SHOCKS: t0 = idx.get_indexer([pd.Timestamp(s["start"])])[0] t1 = min(N_WEEKS, t0 + s["weeks"] + buffer) sign = 1 if s["size"] > 0 else -1 if s["id"] == "S4": sign = -1 # road cut lowers valley prices for k in SHOCK_TARGETS[s["id"]]: hit = np.nonzero(flags[k][t0:t1] == sign)[0] det.append(dict(shock=s["id"], start=s["start"], kind=s["kind"], series=k, direction="high" if sign > 0 else "low", detected=bool(len(hit)), weeks_to_flag=(int(hit[0]) if len(hit) else None))) rows = [] for k in KEYS: crop = k.split("_")[0] in_window = np.zeros(N_WEEKS, dtype=bool) for s in SHOCKS: if s["crop"] in (crop, "both"): t0 = idx.get_indexer([pd.Timestamp(s["start"])])[0] in_window[t0:min(N_WEEKS, t0 + s["weeks"] + buffer)] = True valid = ~np.isnan(norms[k]["norm"]) for per, sl in (("reference", slice(0, N_REF)), ("monitoring", slice(N_REF, N_WEEKS))): f, w, v = flags[k][sl] != 0, in_window[sl], valid[sl] rows.append(dict(series=k, period=per, weeks=int(v.sum()), flagged=int(f.sum()), flagged_outside_windows=int((f & ~w).sum()), weeks_outside_windows=int((~w & v).sum()))) tab = pd.DataFrame(rows) tab["false_flag_rate"] = (tab["flagged_outside_windows"] / tab["weeks_outside_windows"]).round(3) return tab, pd.DataFrame(det) def rule_tradeoff(df, ks=(1.5, 2.0, 2.5), ns=(1, 2, 3)) -> pd.DataFrame: """Shocks detected (of 9 shock-series pairs) and false-flag rate outside shock windows, by k and n.""" norms = seasonal_norms(df) rows = [] for k in ks: for n in ns: r = AlertRule(k=k, n=n) nr = {kk: dict(v, sigma=v["sigma"]) for kk, v in norms.items()} tab, det = evaluate_alerts(df, nr, r) rows.append(dict(k=k, n=n, detected=int(det["detected"].sum()), pairs=len(det), median_weeks_to_flag=float(det["weeks_to_flag"].median()), false_flag_rate_reference=round(tab[tab.period == "reference"]["flagged_outside_windows"].sum() / tab[tab.period == "reference"]["weeks_outside_windows"].sum(), 3), false_flag_rate_monitoring=round(tab[tab.period == "monitoring"]["flagged_outside_windows"].sum() / tab[tab.period == "monitoring"]["weeks_outside_windows"].sum(), 3))) return pd.DataFrame(rows) # ---------------------------------------------------------------- spreads and margins MILLING_RATIO = GenParams().milling_ratio MILLING_COST_0 = GenParams().milling_cost def farmgate_equivalent(df: pd.DataFrame, crop: str) -> np.ndarray: """Farm-gate price in wholesale-product terms. Rice: (paddy + milling cost) / milling ratio, milling cost indexed at 4% a year; onion: the farm-gate price itself.""" if crop == "onion": return df["onion_farmgate"].values infl = np.exp(np.log(1.04) * df["week"].values / 52.0) return (df["rice_farmgate"].values + MILLING_COST_0 * infl) / MILLING_RATIO def spreads_margins(df: pd.DataFrame) -> pd.DataFrame: out = pd.DataFrame({"week": df["week"], "date": df["date"]}) for crop in ("rice", "onion"): cap = df[f"{crop}_capital"].values out[f"{crop}_spread"] = (cap - df[f"{crop}_valley"].values).round(2) out[f"{crop}_margin_share"] = ((cap - farmgate_equivalent(df, crop)) / cap).round(4) out["rice_local_to_parity"] = (df["rice_capital"] / df["rice_parity"]).round(4) return out def same_week_reference_mean(x: np.ndarray, t: int) -> float: """Mean of x at t-52, t-104, ... restricted to the reference period (the 5-year norm).""" vals = [x[j] for j in range(t - 52, -1, -52) if j < N_REF] return float(np.mean(vals)) # ---------------------------------------------------------------- price transmission def transmission(df: pd.DataFrame, crop: str, seasonal: bool = False) -> dict: """Engle-Granger two-step on Valley market A (v) and Capital wholesale (c), levels, reference period. Long run: v_t = b0 + b1 c_t [+ two Fourier pairs of period 52 if seasonal] + e_t ECM: dv_t = a0 + alpha e_{t-1} + gamma dc_t + delta dv_{t-1} + u_t (HAC standard errors, 4 lags) Half-life of a deviation = ln(0.5) / ln(1 + alpha) weeks. The Engle-Granger p-value (MacKinnon) is reported for the specification without Fourier terms. """ import statsmodels.api as sm from statsmodels.tsa.stattools import adfuller, coint v = df[f"{crop}_valley"].values[:N_REF] c = df[f"{crop}_capital"].values[:N_REF] t = np.arange(N_REF) X = c[:, None] if seasonal: F = np.column_stack([f(2 * np.pi * j * t / PERIOD) for j in (1, 2) for f in (np.sin, np.cos)]) X = np.column_stack([c, F]) lr = sm.OLS(v, sm.add_constant(X)).fit() e = lr.resid t_stat, p_coint, _ = coint(v, c, trend="c") dv, dc = np.diff(v), np.diff(c) Z = np.column_stack([e[1:-1], dc[1:], dv[:-1]]) ecm = sm.OLS(dv[1:], sm.add_constant(Z)).fit(cov_type="HAC", cov_kwds={"maxlags": 4}) alpha = float(ecm.params[1]) g = GenParams() return dict( crop=crop, spec="with seasonal terms" if seasonal else "plain", adf_p_valley=float(adfuller(v, regression="ct", autolag="AIC")[1]), adf_p_capital=float(adfuller(c, regression="ct", autolag="AIC")[1]), b0=float(lr.params[0]), b1=float(lr.params[1]), eg_t=float(t_stat), eg_p=float(p_coint), alpha=alpha, alpha_se=float(ecm.bse[1]), gamma=float(ecm.params[2]), gamma_se=float(ecm.bse[2]), half_life_weeks=float(math.log(0.5) / math.log(1 + alpha)) if -1 < alpha < 0 else float("nan"), true_alpha=-(g.rice_alpha if crop == "rice" else g.onion_alpha), true_gamma=g.rice_gamma if crop == "rice" else g.onion_gamma, ) # ---------------------------------------------------------------- storage decision @dataclass(frozen=True) class Storage: buy_month: int # 1..12 sell_month: int # 1..12, after buy_month (wraps into next year) loss: float # weight loss per month held (share) finance: float # financing cost per month, simple interest on purchase value fee: float # storage fee, LCU per tonne per month STORAGE_DEFAULTS = { "rice": Storage(buy_month=12, sell_month=8, loss=0.005, finance=0.015, fee=250.0), # paddy "onion": Storage(buy_month=4, sell_month=9, loss=0.05, finance=0.015, fee=400.0), } def months_held(buy: int, sell: int) -> int: m = (sell - buy) % 12 return 12 if m == 0 else m def storage_margin(monthly_norm, drift_week: float, s: Storage) -> dict: """Net margin per tonne stored (LCU/t), at seasonal-norm prices. P_buy = monthly_norm[buy]; P_sell = monthly_norm[sell] x exp(g m), g = drift_week x 52/12 sold = (1 - loss)^m tonnes per tonne stored net = 1000 P_sell sold - 1000 P_buy (1 + finance m) - fee m """ m = months_held(s.buy_month, s.sell_month) g = drift_week * 52.0 / 12.0 p_buy = float(monthly_norm[s.buy_month - 1]) p_sell = float(monthly_norm[s.sell_month - 1]) * math.exp(g * m) sold = (1 - s.loss) ** m revenue = 1000.0 * p_sell * sold cost_goods = 1000.0 * p_buy finance = cost_goods * s.finance * m fees = s.fee * m net = revenue - cost_goods - finance - fees return dict(months=m, p_buy=p_buy, p_sell=p_sell, sold=sold, revenue=revenue, cost_goods=cost_goods, finance=finance, fees=fees, lost_value=1000.0 * p_sell * (1 - sold), net=net, return_on_stock=net / cost_goods) def storage_curve(monthly_norm, drift_week: float, s: Storage) -> pd.DataFrame: """Net margin for every possible sell month after the buy month (1..11 months held).""" rows = [] for m in range(1, 12): sell = (s.buy_month - 1 + m) % 12 + 1 r = storage_margin(monthly_norm, drift_week, Storage(s.buy_month, sell, s.loss, s.finance, s.fee)) rows.append(dict(sell_month=sell, **r)) cur = pd.DataFrame(rows) return cur def break_even(cur: pd.DataFrame): pos = cur[cur["net"] >= 0] return None if pos.empty else int(pos.iloc[0]["sell_month"]) def break_even_finance(r: dict, s: Storage) -> float: """Monthly financing rate at which the plan breaks even: (revenue - goods - fees) / (goods x months).""" return (r["revenue"] - r["cost_goods"] - r["fees"]) / (r["cost_goods"] * r["months"]) def break_even_loss(monthly_norm, drift_week: float, s: Storage) -> float: """Monthly loss rate at which the default buy/sell plan just breaks even (bisection).""" lo, hi = 0.0, 0.5 f = lambda L: storage_margin(monthly_norm, drift_week, Storage(s.buy_month, s.sell_month, L, s.finance, s.fee))["net"] if f(lo) < 0: return float("nan") for _ in range(100): mid = 0.5 * (lo + hi) if f(mid) >= 0: lo = mid else: hi = mid return lo # ---------------------------------------------------------------- forecasting MODELS = ("seasonal_naive", "sarima", "stl_ets") HORIZON = 12 Z80, Z95 = 1.2815515655446004, 1.959963984540054 def _fc_seasonal_naive(y: np.ndarray, h: int): """log y_{T+j} = log y_{T+j-52}; interval from the SD of 52-week log differences.""" d = y[52:] - y[:-52] sd = d.std(ddof=1) T = len(y) mean = np.array([y[T - 52 + j] for j in range(h)]) sd_h = np.full(h, sd) return mean, sd_h def _fc_sarima(y: np.ndarray, h: int): """SARIMA(1,0,0)(0,1,0)[52] with drift on log prices.""" from statsmodels.tsa.statespace.sarimax import SARIMAX with warnings.catch_warnings(): warnings.simplefilter("ignore") res = SARIMAX(y, order=(1, 0, 0), seasonal_order=(0, 1, 0, 52), trend="c").fit(disp=False, maxiter=200) f = res.get_forecast(h) return np.asarray(f.predicted_mean), np.sqrt(np.asarray(f.var_pred_mean)) def _fc_stl_ets(y: np.ndarray, h: int): """STL (robust, period 52) + ETS(A,Ad,N) on the seasonally adjusted log series.""" from statsmodels.tsa.forecasting.stl import STLForecast from statsmodels.tsa.exponential_smoothing.ets import ETSModel y = pd.Series(y, index=pd.RangeIndex(len(y))) mod = STLForecast(y, ETSModel, period=PERIOD, robust=True, model_kwargs=dict(error="add", trend="add", damped_trend=True)) res = mod.fit(fit_kwargs=dict(disp=False)) f = res.get_prediction(start=len(y), end=len(y) + h - 1) sf = f.summary_frame(alpha=0.2) mean = np.asarray(sf["mean"]) sd = (np.asarray(sf["mean_ci_upper"]) - mean) / Z80 # 80% prediction interval -> SD return mean, sd FORECASTERS = {"seasonal_naive": _fc_seasonal_naive, "sarima": _fc_sarima, "stl_ets": _fc_stl_ets} def backtest(df: pd.DataFrame, origins=None, h: int = HORIZON) -> pd.DataFrame: """Rolling origin over the monitoring year: fit on weeks <= origin, forecast 1..h ahead. Errors in price levels (median forecast exp(mean log)); 80% interval coverage.""" if origins is None: origins = list(range(N_REF - 1, N_WEEKS - h, 4)) # 259, 263, ..., 299 rows = [] for k in KEYS: y = np.log(df[k].values) for o in origins: for m in MODELS: mean, sd = FORECASTERS[m](y[: o + 1], h) for j in range(h): actual = df[k].values[o + 1 + j] pred = math.exp(mean[j]) lo, hi = math.exp(mean[j] - Z80 * sd[j]), math.exp(mean[j] + Z80 * sd[j]) rows.append(dict(series=k, origin=o, model=m, h=j + 1, actual=actual, pred=pred, ape=abs(pred - actual) / actual, covered80=int(lo <= actual <= hi))) return pd.DataFrame(rows) def backtest_summary(bt: pd.DataFrame) -> pd.DataFrame: bt = bt.copy() bt["hgroup"] = np.where(bt["h"] <= 4, "h1_4", "h5_12") a = bt.pivot_table(index=["series", "model"], columns="hgroup", values="ape", aggfunc="mean") a["mape_all"] = bt.groupby(["series", "model"])["ape"].mean() a["cover80"] = bt.groupby(["series", "model"])["covered80"].mean() return (a * 100).round(1).reset_index() def best_models(bs: pd.DataFrame) -> dict: """Per crop, the model with the lowest backtest MAPE averaged over that crop's series (1-12 weeks).""" bs = bs.assign(crop=bs["series"].str.split("_").str[0]) return {c: g.groupby("model")["mape_all"].mean().idxmin() for c, g in bs.groupby("crop")} def forecast(df: pd.DataFrame, models: dict, h: int = HORIZON) -> dict: """12-week forecast bands for every series with its crop's model (dict crop -> model name).""" out = {} for k in KEYS: y = np.log(df[k].values) mean, sd = FORECASTERS[models[k.split("_")[0]]](y, h) out[k] = dict( median=np.round(np.exp(mean), 2), lo80=np.round(np.exp(mean - Z80 * sd), 2), hi80=np.round(np.exp(mean + Z80 * sd), 2), lo95=np.round(np.exp(mean - Z95 * sd), 2), hi95=np.round(np.exp(mean + Z95 * sd), 2), ) return out # ---------------------------------------------------------------- headline numbers def headline(df, norms, rule: AlertRule = AlertRule(), fc=None) -> dict: T = N_WEEKS - 1 sm_ = spreads_margins(df) out = {"last_date": df["date"].iloc[-1]} for k in KEYS: p = df[k].values flag = alerts(p, norms[k]["norm"], norms[k]["sigma"], rule) out[k] = dict(latest=float(p[T]), norm=float(norms[k]["norm"][T]), dev=float(p[T] / norms[k]["norm"][T] - 1), flagged_monitoring=int((flag[N_REF:] != 0).sum()), flag_now=int(flag[T])) if fc is not None: out[k]["fc4"] = (float(fc[k]["median"][3]), float(fc[k]["lo80"][3]), float(fc[k]["hi80"][3])) for crop in ("rice", "onion"): sp = sm_[f"{crop}_spread"].values mg = sm_[f"{crop}_margin_share"].values cap = df[f"{crop}_capital"].values out[f"{crop}_spread"] = dict(latest=float(sp[T]), share=float(sp[T] / cap[T]), norm=same_week_reference_mean(sp, T)) out[f"{crop}_margin"] = dict(latest=float(mg[T]), norm=same_week_reference_mean(mg, T)) s = STORAGE_DEFAULTS[crop] k = f"{crop}_farmgate" r = storage_margin(norms[k]["monthly_norm"], norms[k]["drift_week"], s) cur = storage_curve(norms[k]["monthly_norm"], norms[k]["drift_week"], s) best = cur.loc[cur["net"].idxmax()] out[f"{crop}_storage"] = dict(net=r["net"], ret=r["return_on_stock"], months=r["months"], p_buy=r["p_buy"], p_sell=r["p_sell"], break_even_month=break_even(cur), best_month=int(best["sell_month"]), best_net=float(best["net"]), paying_months=[int(x) for x in cur[cur["net"] >= 0]["sell_month"]], break_even_finance=break_even_finance(r, s), break_even_loss=break_even_loss(norms[k]["monthly_norm"], norms[k]["drift_week"], s)) return out # ---------------------------------------------------------------- figures def _plt(): import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt plt.rcParams.update({"font.size": 10, "axes.spines.top": False, "axes.spines.right": False}) return plt LABEL = {k: lab for k, _, _, lab in SERIES} def fig_bands(df, norms, out: Path, rule: AlertRule = AlertRule()): """Price, seasonal norm, band and flags for the two capital wholesale series and two farm-gate series.""" plt = _plt() x = pd.to_datetime(df["date"]) keys = ["rice_farmgate", "rice_capital", "onion_farmgate", "onion_capital"] fig, axs = plt.subplots(4, 1, figsize=(10, 11), sharex=True) for ax, k in zip(axs, keys): nm = norms[k] lo, hi = band(nm["norm"], nm["sigma"], rule.k) f = alerts(df[k].values, nm["norm"], nm["sigma"], rule) ax.fill_between(x, lo, hi, color="#1f5f63", alpha=0.15, lw=0, label=f"band (norm x exp(+/-{rule.k:g} sigma))") ax.plot(x, nm["norm"], color="#1f5f63", lw=1, ls="--", label="seasonal norm") ax.plot(x, df[k], color="#16181d", lw=1.2, label="price") for sgn, col in ((1, "#a3471c"), (-1, "#9a7420")): m = f == sgn ax.scatter(x[m], df[k].values[m], s=10, color=col, zorder=3, label="flag: high" if sgn > 0 else "flag: low") ax.axvline(x.iloc[N_REF], color="#4f545e", lw=0.8, ls=":") ax.set_title(LABEL[k], loc="left", fontsize=10) ax.set_ylabel("LCU/kg") axs[0].legend(ncol=5, fontsize=8, loc="upper left", frameon=False) axs[-1].text(x.iloc[N_REF], axs[-1].get_ylim()[0], " monitoring year", va="bottom", fontsize=8, color="#4f545e") fig.tight_layout() fig.savefig(out / "price_bands.png", dpi=130) plt.close(fig) def fig_stl(df, norms, out: Path): plt = _plt() x = pd.to_datetime(df["date"]).iloc[:N_REF] fig, axs = plt.subplots(3, 2, figsize=(10, 6.5), sharex=True) for j, k in enumerate(["rice_farmgate", "onion_farmgate"]): nm = norms[k] axs[0, j].plot(x, df[k].values[:N_REF], color="#16181d", lw=1) axs[0, j].plot(x, np.exp(nm["trend"]), color="#1f5f63", lw=1.5) axs[0, j].set_title(LABEL[k] + ": price and STL trend", loc="left", fontsize=9) axs[1, j].plot(x, 100 * nm["seasonal"], color="#9a7420", lw=1) axs[1, j].set_title("seasonal component (% of trend, log points x 100)", loc="left", fontsize=9) axs[2, j].plot(x, 100 * nm["remainder"], color="#a3471c", lw=0.8) axs[2, j].set_title("remainder (log points x 100)", loc="left", fontsize=9) fig.tight_layout() fig.savefig(out / "stl_decomposition.png", dpi=130) plt.close(fig) def fig_storage(norms, out: Path): plt = _plt() mon = ["Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"] fig, axs = plt.subplots(1, 2, figsize=(10, 3.8)) for ax, crop in zip(axs, ("rice", "onion")): s = STORAGE_DEFAULTS[crop] k = f"{crop}_farmgate" cur = storage_curve(norms[k]["monthly_norm"], norms[k]["drift_week"], s) cols = ["#1f5f63" if v >= 0 else "#a3471c" for v in cur["net"]] ax.bar([mon[m - 1] for m in cur["sell_month"]], cur["net"], color=cols) ax.axhline(0, color="#4f545e", lw=0.8) ax.set_title(f"{LABEL[k]}: bought in {mon[s.buy_month - 1]}, net margin by sell month (LCU/t)", loc="left", fontsize=9) ax.tick_params(axis="x", labelsize=8) fig.tight_layout() fig.savefig(out / "storage_margin.png", dpi=130) plt.close(fig) def fig_forecast(df, fc, out: Path, keys=("rice_capital", "onion_capital"), back: int = 52): plt = _plt() x = pd.to_datetime(df["date"]) xf = pd.date_range(x.iloc[-1] + pd.Timedelta(days=7), periods=HORIZON, freq="7D") fig, axs = plt.subplots(1, len(keys), figsize=(10, 3.6)) for ax, k in zip(axs, keys): ax.plot(x.iloc[-back:], df[k].values[-back:], color="#16181d", lw=1.2) ax.fill_between(xf, fc[k]["lo95"], fc[k]["hi95"], color="#9a7420", alpha=0.15, lw=0, label="95%") ax.fill_between(xf, fc[k]["lo80"], fc[k]["hi80"], color="#9a7420", alpha=0.35, lw=0, label="80%") ax.plot(xf, fc[k]["median"], color="#9a7420", lw=1.2, ls="--", label="median") ax.set_title(LABEL[k] + ": 12-week forecast band", loc="left", fontsize=9) ax.tick_params(axis="x", labelsize=8) axs[0].legend(fontsize=8, frameon=False, loc="upper left") fig.tight_layout() fig.savefig(out / "forecast_bands.png", dpi=130) plt.close(fig) # ---------------------------------------------------------------- export for the website demo def export_ts(path: Path, df, norms, fc, model: dict, bt_sum: pd.DataFrame, rule: AlertRule = AlertRule()): def arr(a): return "[" + ",".join(f"{x:.2f}".rstrip("0").rstrip(".") for x in a) + "]" lines = [ "// Synthetic data for the MarketPriceMonitor demo. GENERATED by", "// public/artifacts/agricultural-market-intelligence/market_intelligence.py --ts-out. Do not edit by hand.", "// Prices in LCU/kg (generic local currency unit). Not the prices of any real market.", "export type SeriesKey = " + " | ".join(f"'{k}'" for k in KEYS) + ";", "export type SeriesData = { price: number[]; sIndex: number[]; sigma: number; driftWeek: number; monthlyNorm: number[];", " fc: { median: number[]; lo80: number[]; hi80: number[]; lo95: number[]; hi95: number[] };", " mape: { seasonal_naive: number; sarima: number; stl_ets: number } };", f"export const START_DATE = '{START.strftime('%Y-%m-%d')}';", f"export const N_WEEKS = {N_WEEKS};", f"export const N_REF = {N_REF};", f"export const FORECAST_MODEL: Record<'rice' | 'onion', 'seasonal_naive' | 'sarima' | 'stl_ets'> = {json.dumps(model)};", f"export const ALERT_DEFAULTS = {{ k: {rule.k}, n: {rule.n}, window: {rule.window}, gap: {rule.gap} }};", "export const DATA: Record = {", ] for k in KEYS: nm = norms[k] f = fc[k] mp = {m: float(bt_sum[(bt_sum.series == k) & (bt_sum.model == m)]["mape_all"].iloc[0]) for m in MODELS} sidx = "[" + ",".join(f"{x:.5f}".rstrip("0").rstrip(".") for x in nm["index"]) + "]" lines.append(f" {k}: {{ price: {arr(df[k].values)}, sIndex: {sidx}, sigma: {nm['sigma']}, " f"driftWeek: {nm['drift_week']}, monthlyNorm: {arr(nm['monthly_norm'])},") lines.append(f" fc: {{ median: {arr(f['median'])}, lo80: {arr(f['lo80'])}, hi80: {arr(f['hi80'])}, " f"lo95: {arr(f['lo95'])}, hi95: {arr(f['hi95'])} }},") lines.append(f" mape: {json.dumps(mp)} }},") lines.append("};") path.write_text("\n".join(lines) + "\n") def parameter_table() -> pd.DataFrame: g = GenParams() rows = [(k, v) for k, v in asdict(g).items()] tab = pd.DataFrame(rows, columns=["parameter", "value"]) tab["status"] = "illustrative" return tab def main(argv=None): ap = argparse.ArgumentParser() ap.add_argument("--out", default=str(Path(__file__).resolve().parent)) ap.add_argument("--ts-out", default=None) a = ap.parse_args(argv) out = Path(a.out) df = generate() df.to_csv(out / "synthetic_prices_weekly.csv", index=False) norms = seasonal_norms(df) tab, det = evaluate_alerts(df, norms) tab.to_csv(out / "alert_false_flags.csv", index=False) det.to_csv(out / "alert_detection.csv", index=False) print(tab.groupby("period")[["flagged", "flagged_outside_windows", "weeks_outside_windows"]].sum()) print(det.to_string(index=False)) tr = pd.DataFrame([transmission(df, c, sz) for c in ("rice", "onion") for sz in (False, True)]) tr.round(4).to_csv(out / "price_transmission.csv", index=False) print(tr.round(3).to_string()) bt = backtest(df) bs = backtest_summary(bt) bs.to_csv(out / "forecast_backtest.csv", index=False) print(bs.groupby("model")[["h1_4", "h5_12", "mape_all", "cover80"]].mean().round(1)) model = best_models(bs) fc = forecast(df, model) hl = headline(df, norms, fc=fc) hl["forecast_model"] = model (out / "headline_numbers.json").write_text(json.dumps(hl, indent=1, default=float)) print(json.dumps(hl, indent=1, default=float)) fig_bands(df, norms, out) fig_stl(df, norms, out) fig_storage(norms, out) fig_forecast(df, fc, out) if a.ts_out: export_ts(Path(a.ts_out), df, norms, fc, model, bs) if __name__ == "__main__": main()