Forecast a FRED Series and Publish a Challenger Only If It Beats Seasonal-Naive

id: fred-forecast-champion-challenger-gate
namespace: company.team

inputs:
  - id: series_id
    type: STRING
    displayName: FRED series ID
    description: Any series on fred.stlouisfed.org. The default is US electric and
      gas utilities output, monthly and strongly seasonal.
    defaults: IPG2211A2N
    validator: ^[A-Za-z0-9_.-]+$

  - id: season_length
    type: INT
    displayName: Season length
    description: Observations per seasonal cycle, for example 12 for monthly data or
      4 for quarterly data.
    defaults: 12
    min: 1

  - id: horizon
    type: INT
    displayName: Forecast horizon
    description: How many steps ahead to forecast. Each backtest fold is scored over
      the same horizon.
    defaults: 12
    min: 1

  - id: n_folds
    type: INT
    displayName: Backtest folds
    description: Number of rolling-origin folds. Each fold moves the cutoff forward
      by one horizon.
    defaults: 5
    min: 1
    max: 20

  - id: min_improvement
    type: FLOAT
    displayName: Minimum MASE improvement
    description: The champion is published only if its MASE is at least this share
      below the seasonal-naive baseline, for example 0.10 for 10%.
    defaults: 0.10
    min: 0
    max: 1

tasks:
  - id: download
    type: io.kestra.plugin.core.http.Download
    description: Download the full history of the series as CSV from FRED's public
      graph endpoint. No API key is needed.
    uri: "https://fred.stlouisfed.org/graph/fredgraph.csv?id={{ inputs.series_id }}"
    # FRED drops requests from Java's default User-Agent and from bare names like "kestra-blueprint",
    # and the task then hangs until it times out. A name with a project URL is accepted.
    headers:
      User-Agent: kestra-blueprint (https://kestra.io)
    options:
      # Let a 404 through so the next task can explain it, instead of failing here with FRED's HTML page.
      allowedResponseCodes: [ 200, 404 ]

  - id: unknown_series
    type: io.kestra.plugin.core.execution.Fail
    description: Stop with a readable message when FRED does not know the series ID.
    condition: "{{ outputs.download.code == 404 }}"
    errorMessage: "FRED has no series '{{ inputs.series_id }}'. Check the ID on
      https://fred.stlouisfed.org."

  - id: clean
    type: io.kestra.plugin.jdbc.duckdb.Queries
    description: >
      Read the CSV as text, turn empty and "." values into NULL, and write the
      remaining observations to clean.csv. Infer the frequency from the median
      spacing between dates, and for weekly or coarser series find the first
      date that is missing from the regular grid.
    inputFiles:
      series.csv: "{{ outputs.download.uri }}"
    outputFiles:
      - clean.csv
    sql: |
      CREATE TABLE raw AS
      SELECT CAST(ds AS DATE) AS ds,
             TRY_CAST(NULLIF(TRIM(value), '.') AS DOUBLE) AS y
      FROM read_csv('series.csv', header = true, all_varchar = true, names = ['ds', 'value']);

      CREATE TABLE obs AS
      SELECT ds, y, lag(ds) OVER (ORDER BY ds) AS prev_ds
      FROM raw
      WHERE y IS NOT NULL;

      -- Daily series skip weekends and holidays, so only weekly and coarser series must be regular.
      CREATE TABLE spacing AS
      SELECT freq,
             CASE freq WHEN 'week' THEN INTERVAL 7 DAY
                       WHEN 'month' THEN INTERVAL 1 MONTH
                       WHEN 'quarter' THEN INTERVAL 3 MONTH
                       WHEN 'year' THEN INTERVAL 1 YEAR END AS step
      FROM (
        SELECT CASE WHEN median_gap <= 1 THEN 'day'
                    WHEN median_gap = 7 THEN 'week'
                    WHEN median_gap BETWEEN 28 AND 31 THEN 'month'
                    WHEN median_gap BETWEEN 89 AND 92 THEN 'quarter'
                    WHEN median_gap BETWEEN 365 AND 366 THEN 'year'
                    ELSE 'irregular' END AS freq
        FROM (SELECT median(ds - prev_ds) AS median_gap FROM obs)
      );

      -- The first statement that returns rows, so later tasks read it as outputs.clean.outputs[0].row
      SELECT count(*) AS n_obs,
             (SELECT count(*) FROM raw WHERE y IS NULL) AS n_missing,
             any_value(freq) AS freq,
             CAST(min(prev_ds + step) FILTER (WHERE ds > prev_ds + step) AS DATE) AS first_missing_date
      FROM obs, spacing;

      COPY (SELECT ds, y FROM obs ORDER BY ds)
        TO '{{ outputFiles["clean.csv"] }}' (HEADER);
    fetchType: FETCH_ONE

  - id: irregular_dates
    type: io.kestra.plugin.core.execution.Fail
    description: >
      Stop if a weekly, monthly, quarterly or yearly series has a hole, because
      dropping it would shift every later observation against the season. Also
      stop if no regular frequency fits.
    condition: "{{ outputs.clean.outputs[0].row.freq == 'irregular' or
      outputs.clean.outputs[0].row.first_missing_date != null }}"
    errorMessage: >-
      {% if outputs.clean.outputs[0].row.freq == 'irregular' %}Series '{{
      inputs.series_id }}' has no regular daily, weekly, monthly, quarterly or
      yearly spacing, so seasonal models cannot be applied. {% else %}Series '{{
      inputs.series_id }}' is {{ outputs.clean.outputs[0].row.freq }}ly but has
      no value for {{ outputs.clean.outputs[0].row.first_missing_date }}. Fill
      or trim the gap before forecasting.{% endif %}

  - id: too_short
    type: io.kestra.plugin.core.execution.Fail
    description: >
      Stop before any modelling if the series cannot fill every backtest fold
      and still leave one full season plus one point of training data for the
      earliest fold.
    condition: "{{ outputs.clean.outputs[0].row.n_obs < inputs.n_folds *
      inputs.horizon + inputs.season_length + 1 }}"
    errorMessage: >-
      Series '{{ inputs.series_id }}' has {{ outputs.clean.outputs[0].row.n_obs
      }} usable observations, but {{ inputs.n_folds }} folds x {{ inputs.horizon
      }} steps plus {{ inputs.season_length }} + 1 training points need {{
      inputs.n_folds * inputs.horizon + inputs.season_length + 1 }}. Lower
      n_folds or horizon, or pick a longer series.

  - id: backtest
    type: io.kestra.plugin.core.flow.Loop
    description: >
      Back-test the seasonal-naive baseline and the two challengers at the same
      time, one iteration per model. Each iteration also forecasts the full
      history, so the chosen forecast is never refitted.
    values: [ "SeasonalNaive", "ETS", "Theta" ]
    concurrencyLimit: 0
    outputs:
      - id: scores
        type: FILE
        value: "{{ outputs.model.outputFiles['scores.csv'] }}"
      - id: forecast
        type: FILE
        value: "{{ outputs.model.outputFiles['forecast.csv'] }}"
    tasks:
      - id: model
        type: io.kestra.plugin.scripts.python.Script
        description: >
          Rolling-origin backtest of one model: train on everything before each
          cutoff, forecast the next horizon, and score it with MASE. Then
          forecast the horizon after the last observation.
        # Fail with a clear error instead of running forever if an install or a fit hangs.
        timeout: PT15M
        containerImage: python:3.13.16-slim-trixie
        inputFiles:
          clean.csv: "{{ outputs.clean.outputFiles['clean.csv'] }}"
          # Every package in the resolved tree is pinned, so a later release cannot change the result.
          requirements.txt: |
            formulaic==1.2.2
            interface_meta==2.0.1
            narwhals==2.27.1
            numpy==2.5.4
            packaging==26.3
            pandas==3.0.6
            patsy==1.0.3
            python-dateutil==2.9.0.post0
            scipy==1.18.1
            six==1.17.0
            statsmodels==0.15.0
            typing_extensions==4.16.0
            wrapt==2.5.1
        beforeCommands:
          # --timeout and --retries make a stalled download retry instead of hanging.
          - pip install --timeout 60 --retries 5 --progress-bar off
            --root-user-action=ignore -r requirements.txt
        env:
          MODEL: "{{ item.value }}"
          SEASON_LENGTH: "{{ inputs.season_length }}"
          HORIZON: "{{ inputs.horizon }}"
          N_FOLDS: "{{ inputs.n_folds }}"
          FREQ: "{{ outputs.clean.outputs[0].row.freq }}"
        outputFiles:
          - scores.csv
          - forecast.csv
        script: |
          import os
          import warnings
          import numpy as np
          import pandas as pd
          from statsmodels.tsa.exponential_smoothing.ets import ETSModel
          from statsmodels.tsa.forecasting.theta import ThetaModel

          name = os.environ["MODEL"]
          m, h, k = int(os.environ["SEASON_LENGTH"]), int(os.environ["HORIZON"]), int(os.environ["N_FOLDS"])
          df = pd.read_csv("clean.csv", parse_dates=["ds"])
          y = df["y"].to_numpy()

          # ETS grid: no trend, additive trend or damped additive trend, each with additive or
          # multiplicative seasonality. Multiplicative needs strictly positive data, so it is skipped otherwise.
          TRENDS = [("none", None, False), ("additive", "add", False), ("damped", "add", True)]
          SEASONALS = ["add", "mul"] if (y > 0).all() else ["add"]

          def seasonal_naive(train):
              return np.resize(train[-m:], h), "last season repeated"

          def ets(train):
              fits = []
              for label, trend, damped in TRENDS:
                  for seasonal in SEASONALS:
                      with warnings.catch_warnings():
                          warnings.simplefilter("ignore")   # convergence warnings from weaker configurations
                          fit = ETSModel(train, error="add", trend=trend, damped_trend=damped,
                                         seasonal=seasonal, seasonal_periods=m).fit(disp=False)
                      fits.append((fit.aic, f"trend={label}, seasonal={seasonal}", fit))
              _, config, best = min(fits, key=lambda f: f[0])   # lowest AIC wins
              return np.asarray(best.forecast(h)), config

          def theta(train):
              # Deseasonalises only if the lag-m autocorrelation is significant at the 10% level
              return np.asarray(ThetaModel(train, period=m).fit().forecast(h)), f"period={m}"

          forecast = {"SeasonalNaive": seasonal_naive, "ETS": ets, "Theta": theta}[name]

          scores = []
          for fold in range(1, k + 1):
              cutoff = len(y) - (k - fold + 1) * h   # the last fold ends at the last observation
              train, test = y[:cutoff], y[cutoff:cutoff + h]
              pred, config = forecast(train)
              # MASE: out-of-sample MAE divided by the in-sample MAE of a seasonal-naive forecast
              scale = np.mean(np.abs(train[m:] - train[:-m]))
              if scale == 0:
                  raise SystemExit(f"Fold {fold}: the training data repeats every season exactly, so MASE is undefined.")
              mae = np.mean(np.abs(test - pred))
              scores.append({"model": name, "fold": fold, "cutoff_date": df["ds"].iloc[cutoff - 1].date(),
                             "config": config, "mae": mae, "mase": mae / scale})
              print(f"{name} fold {fold}: MASE {mae / scale:.3f} [{config}]")
          pd.DataFrame(scores).to_csv("scores.csv", index=False)

          # Next dates after the last observation. Daily series without weekend values use business days.
          step = {"week": pd.DateOffset(weeks=1), "month": pd.DateOffset(months=1),
                  "quarter": pd.DateOffset(months=3), "year": pd.DateOffset(years=1),
                  "day": pd.offsets.BDay() if (df["ds"].dt.dayofweek < 5).all() else pd.DateOffset(days=1)}[os.environ["FREQ"]]
          dates = pd.date_range(df["ds"].iloc[-1], periods=h + 1, freq=step)[1:]
          mean, config = forecast(y)
          pd.DataFrame({"model": name, "step": range(1, h + 1), "ds": dates.date, "yhat": mean}).to_csv("forecast.csv", index=False)
          print(f"{name}: mean MASE {np.mean([s['mase'] for s in scores]):.3f} over {k} folds; final fit [{config}]")

  - id: decide
    type: io.kestra.plugin.jdbc.duckdb.Queries
    description: >
      Average MASE per model, pick the champion, apply the threshold, and keep
      only the forecast of the model actually used. The decision, all fold
      scores and the forecast are saved as a DuckDB file.
    # One scores file and one forecast file per Loop iteration, named after the model
    inputFiles: |
      {
      {%- for it in outputs.backtest.outputs %}
        "{{ it.item.value }}_scores.csv": "{{ it.outputs.scores }}",
        "{{ it.item.value }}_forecast.csv": "{{ it.outputs.forecast }}"{% if not loop.last %},{% endif %}
      {%- endfor %}
      }
    sql: |
      CREATE TABLE backtest_scores AS
      SELECT * FROM read_csv(['SeasonalNaive_scores.csv', 'ETS_scores.csv', 'Theta_scores.csv']);

      CREATE TABLE candidate_forecasts AS
      SELECT * FROM read_csv(['SeasonalNaive_forecast.csv', 'ETS_forecast.csv', 'Theta_forecast.csv']);

      CREATE TABLE decision AS
      WITH model_mase AS (
        SELECT model, avg(mase) AS mase FROM backtest_scores GROUP BY model
      ), scored AS (
        SELECT c.model AS champion, c.mase AS champion_mase, b.mase AS baseline_mase,
               1 - c.mase / b.mase AS improvement,
               -- Publish only a challenger that beats the baseline by more than the threshold
               NOT (c.model <> 'SeasonalNaive' AND c.mase < b.mase * (1 - {{ inputs.min_improvement }})) AS is_fallback
        FROM (SELECT * FROM model_mase ORDER BY mase, model <> 'SeasonalNaive' LIMIT 1) AS c,  -- an exact tie goes to the baseline
             (SELECT mase FROM model_mase WHERE model = 'SeasonalNaive') AS b
      )
      SELECT '{{ inputs.series_id }}' AS series_id,
             TIMESTAMPTZ '{{ execution.startDate }}' AS run_at,
             champion,
             champion_mase,
             baseline_mase,
             CAST({{ inputs.min_improvement }} AS DOUBLE) AS min_improvement,
             CASE WHEN is_fallback THEN 'SeasonalNaive' ELSE champion END AS model_used,
             is_fallback,
             CASE
               WHEN champion = 'SeasonalNaive'
                 THEN printf('No challenger beat the seasonal-naive baseline (MASE %.3f).', baseline_mase)
               WHEN is_fallback
                 THEN printf('%s improved MASE by only %.1f%% (%.3f vs baseline %.3f); more than %.1f%% is required.',
                             champion, improvement * 100, champion_mase, baseline_mase, {{ inputs.min_improvement }} * 100)
               ELSE printf('%s improved MASE by %.1f%% (%.3f vs baseline %.3f); more than %.1f%% is required.',
                           champion, improvement * 100, champion_mase, baseline_mase, {{ inputs.min_improvement }} * 100)
             END AS reason
      FROM scored;

      CREATE TABLE forecast AS
      SELECT d.series_id, d.run_at, d.model_used AS model, d.is_fallback, f.step, f.ds, f.yhat
      FROM candidate_forecasts AS f
      JOIN decision AS d ON f.model = d.model_used
      ORDER BY f.step;

      -- The only statement that returns rows: outputs.decide.outputs[0].row
      SELECT * FROM decision;
    fetchType: FETCH_ONE
    outputDbFile: true

  - id: gate
    type: io.kestra.plugin.core.flow.If
    description: Report the decision. The forecast table already holds the model
      that was used.
    condition: "{{ outputs.decide.outputs[0].row.is_fallback == false }}"
    then:
      - id: publish
        type: io.kestra.plugin.core.log.Log
        description: The challenger beat the baseline by more than the threshold.
        level: INFO
        message: >-
          Publishing the {{ outputs.decide.outputs[0].row.model_used }} forecast
          for {{ inputs.series_id }}. {{ outputs.decide.outputs[0].row.reason }}
    else:
      - id: fallback
        type: io.kestra.plugin.core.log.Log
        description: Logged at WARN so fallbacks are easy to filter for.
        level: WARN
        message: >-
          Falling back to the SeasonalNaive forecast for {{ inputs.series_id }}.
          {{ outputs.decide.outputs[0].row.reason }}

outputs:
  # Not required: if the run stops before `decide`, only the real error is reported.
  - id: model_used
    type: STRING
    required: false
    value: "{{ outputs.decide.outputs[0].row.model_used }}"
  - id: is_fallback
    type: BOOL
    required: false
    value: "{{ outputs.decide.outputs[0].row.is_fallback }}"
  - id: reason
    type: STRING
    required: false
    value: "{{ outputs.decide.outputs[0].row.reason }}"
  - id: results
    type: FILE
    required: false
    description: DuckDB file with the decision, backtest_scores, candidate_forecasts
      and forecast tables.
    value: "{{ outputs.decide.databaseUri }}"

triggers:
  - id: monthly
    type: io.kestra.plugin.core.trigger.Schedule
    description: Run on the 20th of each month, after the Federal Reserve's
      mid-month industrial production release.
    cron: "0 9 20 * *"
    timezone: America/New_York