This script fits the classical models, reads the Chronos forecasts, and scores all seven models through one identical MASE/RMSSE code path. It reads the train/test splits produced by Scripts/horserace/export_aus_retail.R (which exports the real aus_retail data from tsibbledata).

# =============================================================================
# score_horserace.R — Forecasting horse race scoring (single code path)
# =============================================================================
# Fits classical models on the aus_retail train split, reads the Chronos
# foundation-model forecasts, and scores EVERY model through ONE MASE/RMSSE
# function so no model gets an unfair code path.
#
#   Classical (fable): SNAIVE, NAIVE, ETS (auto), ARIMA (auto), ETS+ARIMA combo
#   Foundation (read): amazon/chronos-bolt-small, amazon/chronos-bolt-base
#
# Scoring (Hyndman & Koehler 2006, seasonal m = 12):
#   scale (MASE denom) = mean_{t=m+1..n} |Y_t - Y_{t-m}|   on TRAIN
#   MASE  = mean_{h}|actual_h - fcast_h| / scale
#   RMSSE = sqrt(mean_{h}(actual_h - fcast_h)^2) / sqrt(mean_{t}(Y_t-Y_{t-m})^2)
#
# Run from project root:
#   Rscript Scripts/horserace/score_horserace.R
# =============================================================================

suppressPackageStartupMessages({
  library(fpp3)
  library(dplyr)
  library(tidyr)
  library(readr)
  library(jsonlite)
})

set.seed(42)

res_dir <- "Research/horserace_results"
M       <- 12L   # seasonal period (monthly)

# --- 1. Load data -----------------------------------------------------------
to_ym <- function(x) yearmonth(as.Date(paste0(x, "-01")))

train_raw <- read_csv(file.path(res_dir, "train.csv"), show_col_types = FALSE)
test_raw  <- read_csv(file.path(res_dir, "test.csv"),  show_col_types = FALSE)
chronos   <- read_csv(file.path(res_dir, "chronos_forecasts.csv"),
                      show_col_types = FALSE)

train_ts <- train_raw |>
  mutate(month = to_ym(month)) |>
  as_tsibble(index = month, key = series_id)

n_series <- n_distinct(train_ts$series_id)
cat(sprintf("Loaded %d series | train rows %d | test rows %d | chronos rows %d\n",
            n_series, nrow(train_raw), nrow(test_raw), nrow(chronos)))

# --- 2. Fit classical models (timed separately for the effort comparison) ---
fit_forecast <- function(formula_expr, label) {
  t <- system.time({
    mbl <- train_ts |> model(m = formula_expr)
    fc  <- mbl |> forecast(h = 24)
  })
  cat(sprintf("  %-8s fitted+forecast in %5.1f s\n", label, t[["elapsed"]]))
  list(
    forecast = fc |>
      as_tibble() |>
      transmute(series_id, month, model = label, forecast = .mean),
    seconds = unname(t[["elapsed"]])
  )
}

cat("Fitting classical models on train (through 2016-12), h = 24:\n")
snaive_r <- fit_forecast(SNAIVE(turnover ~ lag("year")), "snaive")
naive_r  <- fit_forecast(NAIVE(turnover),                 "naive")
ets_r    <- fit_forecast(ETS(turnover),                   "ets")
arima_r  <- fit_forecast(ARIMA(turnover),                 "arima")

# ETS + ARIMA average combination — point forecast is the mean of the two
# point forecasts, so we derive it directly (timing = ets + arima cost).
combo_fc <- ets_r$forecast |>
  rename(f_ets = forecast) |>
  inner_join(arima_r$forecast |> rename(f_arima = forecast),
             by = c("series_id", "month")) |>
  transmute(series_id, month, model = "ets_arima_combo",
            forecast = (f_ets + f_arima) / 2)
combo_seconds <- ets_r$seconds + arima_r$seconds
cat(sprintf("  %-8s derived from ets+arima (%5.1f s combined)\n",
            "combo", combo_seconds))

# --- 3. Build ONE long forecast table (fable + chronos) ---------------------
# Normalise everyone to (series_id, month=char "YYYY-MM", model, forecast).
fable_long <- bind_rows(
  snaive_r$forecast, naive_r$forecast, ets_r$forecast,
  arima_r$forecast, combo_fc
) |>
  mutate(month = format(as.Date(month), "%Y-%m"))

chronos_long <- chronos |>
  transmute(series_id, month = as.character(month), model, forecast)

all_fc <- bind_rows(fable_long, chronos_long)

stopifnot(!anyNA(all_fc$forecast))
cat(sprintf("Combined forecast table: %d rows, %d models, %d series\n",
            nrow(all_fc), n_distinct(all_fc$model), n_distinct(all_fc$series_id)))

# --- 4. Per-series scale factors from TRAIN (seasonal naive, m = 12) --------
scales_tbl <- train_raw |>
  arrange(series_id, month) |>
  group_by(series_id) |>
  summarise(
    mae_scale  = mean(abs(turnover - lag(turnover, M)),   na.rm = TRUE),
    rmse_scale = sqrt(mean((turnover - lag(turnover, M))^2, na.rm = TRUE)),
    .groups = "drop"
  )
stopifnot(all(scales_tbl$mae_scale > 0))

# --- 5. THE single scoring code path ----------------------------------------
# Join forecasts to test actuals + per-series scale, then every model is
# scored by exactly this pipeline. No branch is model-specific.
scored <- all_fc |>
  inner_join(test_raw, by = c("series_id", "month")) |>   # brings in turnover (actual)
  inner_join(scales_tbl, by = "series_id") |>
  rename(actual = turnover) |>
  arrange(series_id, model, month) |>
  group_by(series_id, model) |>
  mutate(h = row_number()) |>                              # horizon index 1..24
  ungroup()

# sanity: each (series, model) must have all 24 test points
per_pair <- scored |> count(series_id, model)
stopifnot(all(per_pair$n == 24))

mase_by_series <- scored |>
  group_by(series_id, model) |>
  summarise(
    mase  = mean(abs(actual - forecast)) / first(mae_scale),
    rmsse = sqrt(mean((actual - forecast)^2)) / first(rmse_scale),
    .groups = "drop"
  )

# horizon buckets: months 1-6, 7-12, 13-24
bucket_of <- function(h) dplyr::case_when(
  h <= 6  ~ "h01_06",
  h <= 12 ~ "h07_12",
  TRUE    ~ "h13_24"
)
mase_by_bucket <- scored |>
  mutate(bucket = bucket_of(h)) |>
  group_by(series_id, model, bucket) |>
  summarise(mase = mean(abs(actual - forecast)) / first(mae_scale),
            .groups = "drop")

# --- 6. Summaries -----------------------------------------------------------
summary_tbl <- mase_by_series |>
  group_by(model) |>
  summarise(
    mean_mase   = mean(mase),
    median_mase = median(mase),
    mean_rmsse  = mean(rmsse),
    .groups = "drop"
  ) |>
  arrange(mean_mase)

# win rate: per series, model(s) with the minimum MASE. Ties split evenly.
wins <- mase_by_series |>
  group_by(series_id) |>
  mutate(is_min = mase == min(mase),
         win_weight = is_min / sum(is_min)) |>
  ungroup()
n_ties <- wins |>
  group_by(series_id) |>
  summarise(tie = sum(is_min) > 1, .groups = "drop") |>
  summarise(n = sum(tie)) |> pull(n)
winrate_tbl <- wins |>
  group_by(model) |>
  summarise(win_rate = sum(win_weight) / n_series, .groups = "drop") |>
  arrange(desc(win_rate))

bucket_tbl <- mase_by_bucket |>
  group_by(model, bucket) |>
  summarise(mean_mase = mean(mase), .groups = "drop") |>
  pivot_wider(names_from = bucket, values_from = mean_mase) |>
  arrange(h01_06)

timing_tbl <- tibble(
  model   = c("snaive", "naive", "ets", "arima", "ets_arima_combo"),
  seconds = c(snaive_r$seconds, naive_r$seconds, ets_r$seconds,
              arima_r$seconds, combo_seconds)
)

# --- 7. Verification print --------------------------------------------------
cat("\n================ LEADERBOARD (mean MASE, ascending) ================\n")
print(as.data.frame(summary_tbl), digits = 4, row.names = FALSE)
cat("\n================ WIN RATE (share of 148 series) ====================\n")
print(as.data.frame(winrate_tbl |> mutate(win_rate = round(win_rate, 4))),
      row.names = FALSE)
cat(sprintf("Series with ties at the top: %d\n", n_ties))
cat("\n================ MASE BY HORIZON BUCKET ============================\n")
print(as.data.frame(bucket_tbl), digits = 4, row.names = FALSE)
cat("\n================ CLASSICAL FITTING WALL-CLOCK (s) =================\n")
print(as.data.frame(timing_tbl), digits = 3, row.names = FALSE)

# --- 8. Cross-check vs the team lead's rough Python numbers ------------------
rough <- tibble(
  model       = c("chronos_bolt_small", "chronos_bolt_base", "snaive"),
  rough_mean  = c(1.63, 1.47, 1.46),
  rough_median= c(1.45, 1.30, 1.20)
)
check <- summary_tbl |>
  inner_join(rough, by = "model") |>
  mutate(d_mean = round(mean_mase - rough_mean, 3),
         d_median = round(median_mase - rough_median, 3),
         flag = ifelse(abs(d_mean) > 0.02 | abs(d_median) > 0.02, "CHECK", "ok"))
cat("\n================ CROSS-CHECK vs rough Python sanity =================\n")
print(as.data.frame(check |>
  select(model, mean_mase, rough_mean, d_mean,
         median_mase, rough_median, d_median, flag)),
  digits = 4, row.names = FALSE)

# --- 9. Persist outputs -----------------------------------------------------
write_csv(mase_by_series |> arrange(series_id, model),
          file.path(res_dir, "mase_by_series.csv"))

# Full forecast table (fable + chronos) so the images script never refits.
write_csv(all_fc |> arrange(series_id, model, month),
          file.path(res_dir, "forecasts_all.csv"))

summary_out <- summary_tbl |>
  left_join(winrate_tbl, by = "model") |>
  left_join(bucket_tbl,  by = "model") |>
  left_join(timing_tbl,  by = "model")

json_out <- list(
  meta = list(
    n_series          = n_series,
    forecast_rows      = nrow(all_fc),
    rows_per_model     = nrow(all_fc) / n_distinct(all_fc$model),
    h                  = 24L,
    train_span         = paste(min(train_raw$month), max(train_raw$month), sep = " .. "),
    test_span          = paste(min(test_raw$month),  max(test_raw$month),  sep = " .. "),
    mase_definition    = "Hyndman & Koehler 2006, seasonal m=12; denom = in-sample MAE of 12-month seasonal naive on train",
    fable_point        = "mean (.mean, fable default)",
    chronos_point      = "median (q=0.5)",
    n_series_with_ties = n_ties
  ),
  models = summary_out |>
    transmute(model, mean_mase, median_mase, mean_rmsse, win_rate,
              mase_h01_06 = h01_06, mase_h07_12 = h07_12, mase_h13_24 = h13_24,
              fit_seconds = seconds) |>
    purrr::transpose(),
  cross_check = check |>
    select(model, mean_mase, rough_mean, d_mean,
           median_mase, rough_median, d_median, flag) |>
    purrr::transpose()
)
write_json(json_out, file.path(res_dir, "scores_summary.json"),
           pretty = TRUE, auto_unbox = TRUE, digits = 6)

cat(sprintf("\nWrote %s (%d rows)\n",
            file.path(res_dir, "mase_by_series.csv"), nrow(mase_by_series)))
cat(sprintf("Wrote %s\n", file.path(res_dir, "scores_summary.json")))
cat("Done.\n")

House rule: Python is welcome when the subject library is Python-native, and this post is the poster child. Chronos-Bolt ships as a Python package, so the zero-shot forecasts live here. Point forecast is the median (the 0.5 quantile).

#!/usr/bin/env python3
"""
Zero-shot Chronos-Bolt forecasts for the horse race.

Reads Research/horserace_results/train.csv (148 aus_retail series through
2016-12), forecasts h=24 months (2017-01..2018-12) with amazon/chronos-bolt-small
and amazon/chronos-bolt-base, entirely zero-shot (no fine-tuning, no covariates).
Point forecast = median (0.5 quantile). Writes chronos_forecasts.csv and
timing_chronos.json. Run from project root.
"""

import json
import time
from pathlib import Path

import pandas as pd
import torch
from chronos import BaseChronosPipeline

H = 24
OUT = Path("Research/horserace_results")
MODELS = {
    "chronos_bolt_small": "amazon/chronos-bolt-small",
    "chronos_bolt_base": "amazon/chronos-bolt-base",
}

device = "mps" if torch.backends.mps.is_available() else "cpu"
print(f"device: {device}")

train = pd.read_csv(OUT / "train.csv")
series_ids = sorted(train["series_id"].unique())
print(f"series: {len(series_ids)}, train rows: {len(train)}")

contexts = []
for sid in series_ids:
    vals = train.loc[train["series_id"] == sid].sort_values("month")["turnover"].to_numpy()
    contexts.append(torch.tensor(vals, dtype=torch.float32))

months = pd.period_range("2017-01", periods=H, freq="M").strftime("%Y-%m").tolist()

rows = []
timing = {}
for model_key, model_name in MODELS.items():
    pipeline = BaseChronosPipeline.from_pretrained(model_name, device_map=device,
                                                   torch_dtype=torch.float32)
    t0 = time.time()
    quantiles, _ = pipeline.predict_quantiles(inputs=contexts, prediction_length=H,
                                              quantile_levels=[0.5])
    elapsed = time.time() - t0
    timing[model_key] = {"model": model_name, "seconds_forecast": round(elapsed, 2),
                         "device": device, "n_series": len(series_ids), "h": H}
    print(f"{model_key}: {elapsed:.1f}s for {len(series_ids)} series")
    med = quantiles[:, :, 0]  # (n_series, H)
    for i, sid in enumerate(series_ids):
        for j, m in enumerate(months):
            rows.append({"series_id": sid, "month": m, "model": model_key,
                         "forecast": float(med[i, j])})

pd.DataFrame(rows).to_csv(OUT / "chronos_forecasts.csv", index=False)
(OUT / "timing_chronos.json").write_text(json.dumps(timing, indent=2))
print(f"wrote {OUT/'chronos_forecasts.csv'} ({len(rows)} rows) and timing_chronos.json")