"""Forecast engine: statistical model race (statsforecast) + a global LightGBM quantile model. Weeks run Wednesday -> Tuesday (ds = the closing Tuesday). Evaluation is rolling one-week-ahead: each week is forecast from data up to the previous week. Windows (target weeks): calibration 7 Jul - 25 Aug 2026 (prediction bands) validation 1 Sep - 22 Sep 2026 (model selection per series) hold-out 29 Sep - 20 Oct 2026 (reported accuracy; never used for fitting or selection) """ import warnings, numpy as np, pandas as pd, lightgbm as lgb from statsforecast import StatsForecast from statsforecast.models import (Naive, WindowAverage, SimpleExponentialSmoothingOptimized, AutoETS, AutoTheta, CrostonOptimized, CrostonSBA, TSB, ADIDA, IMAPA) from pcommon import * warnings.filterwarnings('ignore') WEEKS = pd.date_range(FIRST_FULL_WEEK_END, DATA_END, freq='W-TUE') # 33 full weeks CAL_T = WEEKS[(WEEKS >= '2026-07-07') & (WEEKS <= '2026-08-25')] VAL_T = WEEKS[(WEEKS >= '2026-09-01') & (WEEKS <= TRAIN_END)] TEST_T = pd.DatetimeIndex(HOLDOUT_WEEKS) FWD = pd.date_range(DATA_END + pd.Timedelta(days=7), periods=12, freq='W-TUE') FAMILY = {'Naive': 'Baseline (repeat last week)', 'MA4': 'Statistical', 'SES': 'Statistical', 'AutoETS': 'Statistical', 'AutoTheta': 'Statistical', 'Croston': 'Intermittent', 'CrostonSBA': 'Intermittent', 'TSB': 'Intermittent', 'ADIDA': 'Intermittent', 'IMAPA': 'Intermittent', 'LightGBM': 'Machine learning'} def stat_models(light=False): m = [Naive(), WindowAverage(window_size=4, alias='MA4')] if not light: m += [SimpleExponentialSmoothingOptimized(alias='SES'), AutoETS(season_length=1, alias='AutoETS'), AutoTheta(season_length=1, alias='AutoTheta'), CrostonOptimized(alias='Croston'), CrostonSBA(alias='CrostonSBA'), TSB(alpha_d=0.1, alpha_p=0.1, alias='TSB'), ADIDA(alias='ADIDA'), IMAPA(alias='IMAPA')] return m # ------------------------------------------------------------------ calendar features per week def week_calendar(weeks): ev = {'ramadan': pd.date_range('2026-02-18', '2026-03-19'), 'eid': pd.DatetimeIndex(list(pd.date_range('2026-03-19', '2026-03-23')) + list(pd.date_range('2026-05-26', '2026-05-29'))), 'natday': pd.DatetimeIndex(['2026-09-23', '2026-02-22'])} rows = [] for w in weeks: days = pd.date_range(w - pd.Timedelta(days=6), w) rows.append({'ds': w, 'ramadan_days': days.isin(ev['ramadan']).sum(), 'eid_days': days.isin(ev['eid']).sum(), 'natday': int(days.isin(ev['natday']).any()), 'payday': int((days.day == 27).any())}) return pd.DataFrame(rows).set_index('ds') ALLW = WEEKS.append(FWD) CAL = week_calendar(ALLW) # ------------------------------------------------------------------ statistical race def stats_cv(df, light=False, n_jobs=2): """Rolling 1-step backtest over the last 16 weeks + fixed-origin 12-week backtests (5 origins) + 12-week forward forecast.""" sf = StatsForecast(models=stat_models(light), freq='W-TUE', n_jobs=n_jobs, fallback_model=Naive()) cv1 = sf.cross_validation(df=df, h=1, n_windows=16, step_size=1) cv12 = sf.cross_validation(df=df, h=12, n_windows=5, step_size=1) fwd = sf.forecast(df=df, h=12) return cv1, cv12, fwd # ------------------------------------------------------------------ LightGBM global model def _panel(df): w = df.pivot(index='unique_id', columns='ds', values='y').reindex(columns=WEEKS).fillna(0.0) return w.index.to_numpy(), w.to_numpy(float) def _features(Y, t, h, cats=None): """Features at origin index t (inclusive) for horizon h. Returns DataFrame + scale + mask of series already started.""" n = Y.shape[0] hist = Y[:, :t + 1] started = (hist > 0).any(axis=1) first = np.where(started, (hist > 0).argmax(axis=1), t) def mean_last(k): lo = max(0, t + 1 - k); return hist[:, lo:].mean(axis=1) ma4, ma8, ma13 = mean_last(4), mean_last(8), mean_last(13) lo8 = max(0, t - 7) std8 = hist[:, lo8:].std(axis=1) nz8 = (hist[:, lo8:] > 0).mean(axis=1) rev = (hist[:, ::-1] > 0) since = np.where(rev.any(axis=1), rev.argmax(axis=1), 99) histmean = np.where(started, hist.sum(axis=1) / np.maximum(t + 1 - first, 1), 0.0) scale = np.where(ma13 > 0, ma13, np.where(histmean > 0, histmean, 1.0)) lag = lambda k: (hist[:, t - k] if t - k >= 0 else np.zeros(n)) / scale tw = ALLW[t + h] if t + h < len(ALLW) else None ow = ALLW[t] X = pd.DataFrame({'l0': lag(0), 'l1': lag(1), 'l2': lag(2), 'l3': lag(3), 'ma4': ma4 / scale, 'ma8': ma8 / scale, 'ma13': ma13 / scale, 'std8': std8 / scale, 'nz8': nz8, 'since_nz': since, 'age': t - first, 'histmean': histmean / scale, 'trend': np.where(ma13 > 0, ma4 / np.maximum(ma13, 1e-9), 1.0), 'log_scale': np.log1p(scale), 'h': h, 't_ramadan': CAL.loc[tw, 'ramadan_days'], 't_eid': CAL.loc[tw, 'eid_days'], 't_natday': CAL.loc[tw, 'natday'], 't_payday': CAL.loc[tw, 'payday'], 'o_ramadan': CAL.loc[ow, 'ramadan_days'], 'o_eid': CAL.loc[ow, 'eid_days'], 'o_natday': CAL.loc[ow, 'natday']}) if cats is not None: X['cat'] = pd.Categorical(cats) return X, scale, started def lgb_fit(Y, cats, train_last_t, horizons, alphas=(0.1, 0.5, 0.9), min_t=3, rounds=350): """Train one quantile model per alpha on all (origin, h) pairs whose target week index <= train_last_t.""" Xs, ys, ws = [], [], [] for h in horizons: for t in range(min_t, train_last_t - h + 1): X, sc, st = _features(Y, t, h, cats) X = X[st]; y = Y[st, t + h] / sc[st] Xs.append(X); ys.append(y); ws.append(sc[st]) X = pd.concat(Xs, ignore_index=True); y = np.concatenate(ys); w = np.concatenate(ws) if 'cat' in X: X['cat'] = pd.Categorical(X['cat'].astype(str)) models = {} for a in alphas: m = lgb.LGBMRegressor(objective='quantile', alpha=a, n_estimators=rounds, learning_rate=0.05, num_leaves=31, min_child_samples=50, subsample=0.8, subsample_freq=1, colsample_bytree=0.8, reg_lambda=1.0, verbose=-1, n_jobs=2) m.fit(X, np.clip(y, 0, 50), sample_weight=w) models[a] = m return models, (X['cat'].cat.categories if 'cat' in X else None) def lgb_predict(models, catlev, Y, cats, t, h): X, sc, st = _features(Y, t, h, cats) if catlev is not None: X['cat'] = pd.Categorical(X['cat'].astype(str), categories=catlev) out = {a: np.clip(m.predict(X), 0, None) * sc for a, m in models.items()} q = np.sort(np.vstack([out[a] for a in sorted(out)]), axis=0) res = {a: q[i] for i, a in enumerate(sorted(out))} for a in res: res[a] = np.where(st, res[a], 0.0) return res, sc def widx(d): return int(np.where(ALLW == pd.Timestamp(d))[0][0]) def lgb_backtests(df, cats_by_id=None): """LightGBM predictions matching the statistical backtests: rolling 1-step for calibration / validation / hold-out, fixed-origin 12-step from 30 Jun - 28 Jul origins, and the 12-week forward forecast.""" ids, Y = _panel(df) cats = None if cats_by_id is None else np.array([cats_by_id.get(i, 'NA') for i in ids]) rows = [] # rolling one-step: three training cut-offs (each fold trained only on targets up to its own start) for fold, targets in (('cal', CAL_T), ('val', VAL_T), ('test', TEST_T)): last_train = widx(targets[0]) - 1 mods, lev = lgb_fit(Y, cats, last_train, horizons=[1]) for tw in targets: t = widx(tw) - 1 p, sc = lgb_predict(mods, lev, Y, cats, t, 1) rows.append(pd.DataFrame({'unique_id': ids, 'ds': tw, 'cutoff': ALLW[t], 'LightGBM': p[0.5], 'LightGBM_q10': p[0.1], 'LightGBM_q90': p[0.9], 'scale': sc})) cv1 = pd.concat(rows, ignore_index=True) # fixed-origin 12-step (for horizon-specific bands of the forward forecast) o0 = widx('2026-06-30') mods, lev = lgb_fit(Y, cats, o0, horizons=list(range(1, 13))) rows = [] for t in range(o0, widx('2026-07-28') + 1): for h in range(1, 13): p, sc = lgb_predict(mods, lev, Y, cats, t, h) rows.append(pd.DataFrame({'unique_id': ids, 'ds': ALLW[t + h], 'cutoff': ALLW[t], 'h': h, 'LightGBM': p[0.5], 'scale': sc})) cv12 = pd.concat(rows, ignore_index=True) # forward tl = widx(DATA_END) mods, lev = lgb_fit(Y, cats, tl, horizons=list(range(1, 13))) rows = [] for h in range(1, 13): p, sc = lgb_predict(mods, lev, Y, cats, tl, h) rows.append(pd.DataFrame({'unique_id': ids, 'ds': ALLW[tl + h], 'h': h, 'LightGBM': p[0.5], 'LightGBM_q10': p[0.1], 'LightGBM_q90': p[0.9], 'scale': sc})) fwd = pd.concat(rows, ignore_index=True) imp = pd.Series(mods[0.5].feature_importances_, index=mods[0.5].feature_name_).sort_values(ascending=False) return cv1, cv12, fwd, imp def scale_at(df, cutoffs): """Mean of the last 13 weeks up to each cutoff (the band scale), per series.""" ids, Y = _panel(df) out = [] for c in cutoffs: t = widx(c) lo = max(0, t - 12) m13 = Y[:, lo:t + 1].mean(axis=1) hm = np.array([Y[i, :t + 1][Y[i, :t + 1] > 0].mean() if (Y[i, :t + 1] > 0).any() else 1.0 for i in range(len(ids))]) out.append(pd.DataFrame({'unique_id': ids, 'cutoff': c, 'scale': np.where(m13 > 0, m13, hm)})) return pd.concat(out, ignore_index=True)