Forecasting Producer Price Inflation with Chronos-2

Univariate, multivariate, probabilistic, and covariate-informed forecasting

Author

Mike Aguilar

1 Overview

Quantitative question: Can a pretrained time-series foundation model forecast producer price inflation more accurately than a conventional autoregressive integrated moving-average (ARIMA) model? Does accuracy improve when the model jointly forecasts related series or incorporates historical upstream input-price inflation?

Economic question: Can we forecast final-demand producer price inflation, and do earlier production stages contain useful leading information?

  • Final versus intermediate demand: The Bureau of Labor Statistics (BLS) publishes many producer price indexes. The headline final-demand index comprises goods, services, and construction. Intermediate-demand indexes track prices for goods and services used as production inputs rather than purchased for final demand.
  • Production-flow stages: The BLS assigns industries to four stages according to predominant production flows. Stage 4 industries primarily sell to final demand; stage 3 industries predominantly supply stage 4, and stage 2 industries predominantly supply stage 3. These are broad relationships, not a rigid sequence followed by every product. The stage indexes used here describe prices of inputs consumed by industries in each stage, not prices of their output.
  • Forecasting intuition: Upstream price changes may precede changes in final-demand inflation. Our analysis tests whether their historical values improve forecasts, not whether they identify causal pass-through.

We will focus on monthly inflation as our primary object of interest: \[ \pi_{j,t}=100\bigl(\log P_{j,t}-\log P_{j,t-1}\bigr), \]

which is approximately a percentage change, expressed in percentage points. The forecast horizon is six consecutive months, and the forecast origin advances by one month at a time. The analysis begins with an origin corresponding to information available through December 2017, allowing us to evaluate January 2018 and subsequent months.

Specification Information used Forecasts
Auto-ARIMA benchmark Own inflation history Three component forecasts and a direct total-PPI forecast
Chronos-2, univariate Own inflation history Three separate component forecasts
Chronos-2, direct total Own history of published total final-demand PPI One direct aggregate forecast
Chronos-2, multivariate Joint history of three final-demand series Three targets jointly
Chronos-2, multivariate + covariates Three target histories + eight historical intermediate-demand series Three targets jointly

Each model supplies point predictions and 10th/90th percentiles (an 80% prediction interval). ARIMA reports a model-based point forecast; Chronos-2’s predictions output represents its median (50th percentile). Thus, RMSE compares the point forecasts supplied by the methods, whereas pinball loss also evaluates forecast quantiles. Chronos-2 is used zero-shot: no pretrained parameters are updated or fine-tuned on this dataset. Background on the model: Chronos-2 model card and Chronos GitHub repository.

ImportantInterpretation of the experiment

This is a pseudo-real-time, retrospective backtest. At each origin we hide subsequent rows, but we use the latest revised historical FRED data rather than contemporaneous data vintages.

Moreover, Chronos-2’s pretraining may have included historical PPI observations. Therefore, the main 2018-onward comparison demonstrates the forecasting workflow but cannot establish clean real-time superiority over ARIMA. A smaller post-release-origin check appears near the end.

2 Reproducible environment and settings

The analysis is contained in this Quarto document. On first use, it downloads 12 PPI series through the official FRED API using a personal key stored in a local text file, then saves a fixed snapshot as data/ppi_fred_snapshot.csv. Subsequent renders reuse that snapshot. Forecasts are cached separately under cache/, so editing the lesson’s prose and figures does not normally require refitting the models. To refresh either dataset or forecasts, use the settings and cache procedures described below.

Install the packages once in the Python environment used by Quarto (Python 3.11 or 3.12 recommended). Install the Quarto CLI separately. The following command is a shell command, not a Python notebook cell:

python -m pip install "chronos-forecasting>=2.0,<3" "statsforecast>=2.0,<3" \
    "pandas>=2.2" "numpy>=1.26" "matplotlib>=3.8" \
    "requests>=2.31" "openpyxl>=3.1" "ipython>=8" "ipykernel>=6" "torch>=2.2"

Render from the project folder using quarto render Chronos2_PPI_Quarto_Lesson.qmd --to html (adjust the filename if necessary). A GPU is optional: the script uses CUDA when available, otherwise the CPU. Downloading the pretrained model weights initially requires internet access. To reproduce archived results, retain the saved FRED snapshot and package-version record. The commands below use shell-style line continuation; on Windows, they can also be entered on one line.

from pathlib import Path
import datetime as dt
import hashlib
import io
import json
import re
import time
import warnings

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import requests
import torch
from IPython.display import display, Markdown
from statsforecast.models import AutoARIMA
from chronos import Chronos2Pipeline
from importlib.metadata import version

# EDIT THIS ONE SETTING to switch between a quick test and the full backtest.
RUN_MODE = "full"              # "quick" or "full"
REFRESH_FRED = False           # True deliberately replaces the frozen snapshot
REFRESH_WEIGHTS = False        # True re-downloads archived BLS weights if possible
FORECAST_HORIZON = 6
FIRST_TARGET_DATE = pd.Timestamp("2018-01-01")
POST_RELEASE_ORIGIN = pd.Timestamp("2025-10-01")
QUANTILES = [0.1, 0.5, 0.9]
CODE_VERSION = "ppi-chronos2-lesson-v2-fred-api"

assert RUN_MODE in {"quick", "full"}
DATA_DIR = Path("data")
CACHE_DIR = Path("cache")
DATA_DIR.mkdir(exist_ok=True)
CACHE_DIR.mkdir(exist_ok=True)
SNAPSHOT_PATH = DATA_DIR / "ppi_fred_snapshot.csv"
METADATA_PATH = DATA_DIR / "ppi_fred_snapshot_metadata.json"

VERSIONS = {p: version(p) for p in
            ["chronos-forecasting", "statsforecast", "pandas", "numpy", "torch"]}
DEVICE = "cuda" if torch.cuda.is_available() else "cpu"
(CACHE_DIR / "package_versions.json").write_text(json.dumps(VERSIONS, indent=2))
Package Version
0 chronos-forecasting 2.3.2
1 statsforecast 2.1.1
2 pandas 2.3.3
3 numpy 2.5.4
4 torch 2.14.1

Execution mode: full · Device: cpu · Forecast horizon: 6 months.

3 Data

The Bureau of Labor Statistics groups final-demand producer prices into goods, services, and construction. Their official aggregate is separately available as total final-demand PPI. In the BLS production-flow system, upstream intermediate-demand indexes organize inputs to stage 2, 3, and 4 producers. We test whether their observed historical inflation rates improve our forecasts of final demand.

The input series are the seasonally adjusted versions of the BLS/FRED intermediate-demand production-flow groupings. The eight covariates do not mean we know prices for the six future months.

Role Economic series FRED code
Target Final-demand goods PPIDGS
Target Final-demand services PPIDSS
Target Final-demand construction PPIDCS
Aggregate, for later evaluation Total final-demand PPI PPIFIS
Historical covariate Inputs to stage 2 goods producers WPSID521
Historical covariate Inputs to stage 2 services producers WPSID522
Historical covariate Inputs to stage 3 goods producers WPSID531
Historical covariate Inputs to stage 3 services producers WPSID532
Historical covariate Inputs to stage 3 construction producers WPSID533
Historical covariate Inputs to stage 4 goods producers WPSID541
Historical covariate Inputs to stage 4 services producers WPSID542
Historical covariate Inputs to stage 4 construction producers WPSID543

3.1 Fetch once, then reuse a fixed snapshot

If you are trying to recreate this exercise for yourself, save your FRED API key, with no quotes or other text, in a local file named fred_api_key.txt beside this .qmd file. Do not share or commit this file. By default, the next cell reuses the existing frozen snapshot. Only the first execution, or an intentional refresh, contacts the official FRED API. It downloads all 12 series, aligns them to monthly dates, and checks for missing values rather than silently filling them.

from requests.adapters import HTTPAdapter
from urllib3.util.retry import Retry

TARGETS = {"Goods": "PPIDGS", "Services": "PPIDSS", "Construction": "PPIDCS"}
AGGREGATE = "TotalFD"
AGGREGATE_ID = "PPIFIS"
COVARIATES = {
    "Stage2_Goods": "WPSID521", "Stage2_Services": "WPSID522",
    "Stage3_Goods": "WPSID531", "Stage3_Services": "WPSID532",
    "Stage3_Construction": "WPSID533", "Stage4_Goods": "WPSID541",
    "Stage4_Services": "WPSID542", "Stage4_Construction": "WPSID543",
}

# Keys are the friendly variable names used by every downstream analysis.
FRED_CODES = {**TARGETS, AGGREGATE: AGGREGATE_ID, **COVARIATES}
assert len(FRED_CODES) == 12 and len(set(FRED_CODES.values())) == 12

FRED_API_URL = "https://api.stlouisfed.org/fred/series/observations"
FRED_START = "2009-11-01"
API_KEY_FILE = Path("fred_api_key.txt")

# A persistent connection and retries make multiple consecutive requests sturdier.
fred_session = requests.Session()
fred_retries = Retry(
    total=4, connect=4, read=4, status=4, backoff_factor=1,
    status_forcelist=[429, 500, 502, 503, 504],
    allowed_methods=["GET"],
)
fred_session.mount("https://", HTTPAdapter(max_retries=fred_retries))


def fetch_fred_series(code, api_key):
    """Return one monthly FRED PPI index as a numeric, date-indexed Series."""
    params = {
        "series_id": code,
        "api_key": api_key,
        "file_type": "json",
        "observation_start": FRED_START,
    }
    try:
        response = fred_session.get(FRED_API_URL, params=params, timeout=(10, 45))
        response.raise_for_status()
        observations = response.json().get("observations", [])
    except requests.RequestException:
        # Avoid exposing the API key embedded in a request URL in an exception.
        raise RuntimeError(f"Could not retrieve FRED series {code}. Check API access.") from None
    except ValueError:
        raise RuntimeError(f"FRED did not return valid JSON for series {code}.") from None

    if not observations:
        raise ValueError(f"No FRED observations returned for {code}.")

    frame = pd.DataFrame(observations)
    if not {"date", "value"}.issubset(frame.columns):
        raise ValueError(f"Missing date/value fields in FRED data for {code}.")
    dates = pd.to_datetime(frame["date"]).dt.to_period("M").dt.to_timestamp()
    values = pd.to_numeric(frame["value"], errors="coerce")
    series = pd.Series(values.to_numpy(), index=pd.DatetimeIndex(dates), name=code)
    if series.index.has_duplicates:
        raise ValueError(f"Duplicate monthly observations returned for {code}.")
    return series.sort_index()


def build_snapshot():
    if not API_KEY_FILE.is_file():
        raise FileNotFoundError(
            "Create fred_api_key.txt beside the Quarto file, containing only your FRED API key."
        )
    api_key = API_KEY_FILE.read_text(encoding="utf-8").strip()
    if not api_key:
        raise ValueError("fred_api_key.txt is empty.")

    # Download every target, the aggregate, and all eight covariates.
    downloaded = {}
    for name, code in FRED_CODES.items():
        print(f"Retrieving {name}: {code}")
        downloaded[name] = fetch_fred_series(code, api_key)

    # Align dates, using exactly the variable names expected later in the QMD.
    wide = pd.concat(downloaded, axis=1).sort_index().loc[FRED_START:]
    wide.index.name = "date"

    incomplete = wide[wide.isna().any(axis=1)]
    if not incomplete.empty:
        print(f"Dropping {len(incomplete)} month(s) lacking at least one series: "
              + ", ".join(incomplete.index.strftime("%Y-%m")[:12]))
    wide = wide.dropna(how="any")
    if wide.empty:
        raise ValueError("There are no complete months across the 12 series.")
    expected = pd.date_range(wide.index.min(), wide.index.max(), freq="MS")
    if not wide.index.equals(expected):
        raise ValueError(
            "An internal monthly gap remains in the FRED data. "
            "Inspect each series instead of interpolating missing observations."
        )
    if not (wide > 0).all().all():
        raise ValueError("PPI index levels must be positive for log changes.")

    # Write only after all 12 series have been obtained and validated.
    wide.to_csv(SNAPSHOT_PATH, float_format="%.12g")
    METADATA_PATH.write_text(json.dumps({
        "downloaded_utc": dt.datetime.now(dt.timezone.utc).isoformat(),
        "source": "Official FRED API (BLS PPI)",
        "codes": FRED_CODES,
        "first_observation": str(wide.index.min().date()),
        "latest_observation": str(wide.index.max().date()),
        "number_of_series": len(wide.columns),
    }, indent=2), encoding="utf-8")
    return wide


if REFRESH_FRED or not SNAPSHOT_PATH.exists():
    build_snapshot()

# Keep these names and transformations consistent with all later sections.
levels = pd.read_csv(SNAPSHOT_PATH, parse_dates=["date"], index_col="date")
missing_columns = set(FRED_CODES) - set(levels.columns)
if missing_columns:
    raise ValueError(
        f"Snapshot is missing {sorted(missing_columns)}. "
        "Set REFRESH_FRED=True to rebuild it from all 12 API series."
    )
levels = levels.loc[:, list(FRED_CODES)].sort_index()
assert levels.index.equals(pd.date_range(levels.index.min(), levels.index.max(), freq="MS"))
assert levels.notna().all().all() and (levels > 0).all().all()
infl = 100 * np.log(levels).diff().dropna()
SNAPSHOT_HASH = hashlib.sha256(SNAPSHOT_PATH.read_bytes()).hexdigest()[:16]
assert len(levels.columns) == 12 and len(infl) > 0
Series FRED ID First inflation month Latest inflation month
0 Goods PPIDGS 2009-12 2026-09
1 Services PPIDSS 2009-12 2026-09
2 Construction PPIDCS 2009-12 2026-09
3 TotalFD PPIFIS 2009-12 2026-09
4 Stage2_Goods WPSID521 2009-12 2026-09
5 Stage2_Services WPSID522 2009-12 2026-09
6 Stage3_Goods WPSID531 2009-12 2026-09
7 Stage3_Services WPSID532 2009-12 2026-09
8 Stage3_Construction WPSID533 2009-12 2026-09
9 Stage4_Goods WPSID541 2009-12 2026-09
10 Stage4_Services WPSID542 2009-12 2026-09
11 Stage4_Construction WPSID543 2009-12 2026-09

Frozen snapshot downloaded: 2026-10-24T14:19:06.990751+00:00. SHA-256 prefix: 4e093ba09b707a19. Inflation observations: 202.

3.2 EDA

3.2.1 Sample statistics

We first examine the magnitude and volatility of monthly inflation in final-demand components and selected upstream input-price indexes. These descriptive statistics characterize the forecasting problem; they are not themselves evidence of predictive relationships.

def plot_history(data):
    fig, ax = plt.subplots(2, 1, figsize=(11, 7), sharex=True)
    for name in TARGETS:
        ax[0].plot(data.index, data[name], label=name, alpha=0.9)
    ax[0].axhline(0, color="gray", lw=0.7)
    ax[0].set(title="Final-demand monthly log inflation", ylabel="Percentage points")
    ax[0].legend(ncol=3)
    for name in ["Stage2_Goods", "Stage3_Goods", "Stage4_Goods"]:
        ax[1].plot(data.index, data[name], label=name, alpha=0.85)
    ax[1].axhline(0, color="gray", lw=0.7)
    ax[1].set(title="Illustrative upstream input-price inflation", ylabel="Percentage points")
    ax[1].legend(ncol=3)
    fig.tight_layout()
    plt.show()

summary = infl[[*TARGETS, *COVARIATES]].agg(["mean", "std", "min", "max"]).T
plot_history(infl)

mean std min max
Goods 0.24 0.78 -2.77 2.61
Services 0.22 0.31 -0.55 1.28
Construction 0.31 0.75 -1.20 5.77
Stage2_Goods 0.19 2.85 -9.10 9.67
Stage2_Services 0.23 0.42 -2.05 1.42
Stage3_Goods 0.27 1.17 -5.30 4.53
Stage3_Services 0.21 0.87 -5.38 3.16
Stage3_Construction 0.22 2.70 -13.35 10.48
Stage4_Goods 0.25 0.47 -0.96 2.17
Stage4_Services 0.25 0.37 -1.83 1.25
Stage4_Construction 0.31 0.79 -2.83 3.33

Interpretation. In this snapshot, Goods has the largest monthly standard deviation (0.78 percentage points), while Services has the smallest (0.31). Common inflation episodes are visible, but visual comovement alone does not show that upstream prices improve future forecasts.

3.2.2 PPI Persistence

Before forecasting producer price inflation, we should examine whether current inflation contains information about future inflation.

The autocorrelation function (ACF) measures the relationship between inflation and its own lagged values:

\[ \rho_j(k)=\operatorname{Corr}(\pi_{j,t},\pi_{j,t-k}) \]

where \(j\) indexes our three final-demand categories and \(k\) is the number of months separating observations.

  • Positive autocorrelation: Inflation tends to persist.
  • Negative autocorrelation: Inflation tends to reverse.
  • Little autocorrelation: Past inflation may offer limited linear predictive information.

We examine the first 12 monthly lags. This helps motivate ARIMA forecasting, which explicitly models temporal dependence.

Chronos-2 also uses historical observations, but does not require us to specify individual lag coefficients. Instead, it learns temporal relationships implicitly from the historical context.

MAX_ACF_LAG = 12

def plot_ppi_acf(data, targets, max_lag=12):
    fig, axes = plt.subplots(
        1, len(targets),
        figsize=(13, 3.8),
        sharey=True
    )

    results = []

    for ax, name in zip(axes, targets):
        series = data[name].dropna()
        lags = np.arange(1, max_lag + 1)

        acf_values = [
            series.autocorr(lag=int(k)) for k in lags
        ]

        # Approximate white-noise reference bands
        bound = 1.96 / np.sqrt(len(series))

        ax.vlines(lags, 0, acf_values, linewidth=2)
        ax.scatter(lags, acf_values, s=20)
        ax.axhline(0, color="black", linewidth=0.8)
        ax.axhline(bound, color="gray", linestyle="--")
        ax.axhline(-bound, color="gray", linestyle="--")

        ax.set_title(name)
        ax.set_xlabel("Lag (months)")
        ax.set_xticks(np.arange(1, max_lag + 1, 2))
        ax.set_ylim(-0.6, 1.0)

        for k, value in zip(lags, acf_values):
            results.append({
                "Target": name,
                "Lag": k,
                "ACF": value
            })

    axes[0].set_ylabel("Autocorrelation")
    fig.suptitle("Persistence in Final-Demand PPI Inflation")
    fig.tight_layout()
    plt.show()

    return pd.DataFrame(results)

acf_results = plot_ppi_acf(
    infl, list(TARGETS), MAX_ACF_LAG
)
Target Construction Goods Services
Lag
1 -0.004 0.405 0.028
2 -0.003 0.147 0.165
3 0.677 0.072 0.195
4 -0.001 0.129 0.188
5 0.018 0.200 0.139
6 0.604 0.156 0.173
7 -0.041 0.240 0.200
8 -0.030 0.138 0.142
9 0.450 0.036 0.214
10 -0.077 0.132 0.085
11 -0.069 0.011 0.048
12 0.216 -0.034 -0.011

Interpretation

  • Goods and services: Goods inflation shows pronounced first-order persistence (lag-one autocorrelation near 0.40). Services inflation exhibits weaker, more diffuse positive autocorrelation across several months.
  • Construction: Autocorrelation is pronounced at lags 3, 6, and 9, even though the series is seasonally adjusted. This suggests a three-month pricing or reporting pattern; the ACF alone cannot establish that it reflects residual calendar seasonality.
  • Modeling implications: ACF helps characterize persistence and can inform traditional lag-based specifications, but it does not select lag orders for Chronos-2. Chronos-2 instead uses the observed sequence in its historical context. The dashed bands are approximate white-noise reference bounds, not a formal forecastability test. Our ARIMA benchmark is deliberately nonseasonal, so construction’s three-month pattern is a reason to consider seasonal ARIMA in a separate robustness exercise.

3.2.3 Do Upstream Prices Lead?

The BLS organizes intermediate-demand PPI series according to production flows. This creates a natural hypothesis: price changes faced by upstream producers may precede price changes in final-demand goods, services, and construction.

We investigate this possibility using cross-correlations between intermediate-demand inflation and subsequent final-demand inflation.

Define:

\[ \rho_{XY}(k)= \operatorname{Corr}(X_{t-k},Y_t), \qquad k=0,\ldots,6 \]

where \(X\) represents intermediate-demand inflation and \(Y\) represents final-demand inflation.

Under this convention:

  • Lag 0: Contemporaneous correlation.
  • Lag 1: Upstream inflation one month earlier.
  • Lag 3: Upstream inflation three months earlier.
  • Lag 6: Upstream inflation six months earlier.

Positive correlations at longer lags would be consistent with upstream price changes preceding downstream inflation.

However, cross-correlation is descriptive. Common shocks, serial dependence, and overlapping economic influences may produce correlations even without a direct transmission mechanism.

MAX_CCF_LAG = 6

def compute_upstream_ccf(data, targets, covariates, max_lag=6):
    results = []

    for target in targets:
        for predictor in covariates:
            for lag in range(max_lag + 1):

                # Predictor observed lag months before target
                correlation = data[target].corr(
                    data[predictor].shift(lag)
                )

                results.append({
                    "Target": target,
                    "Predictor": predictor,
                    "Lag": lag,
                    "Correlation": correlation
                })

    return pd.DataFrame(results)


def plot_upstream_ccf(ccf_results, targets, covariates, max_lag=6):

    fig, axes = plt.subplots(
        1, len(targets),
        figsize=(16, 5),
        sharey=True
    )

    max_abs = max(
        0.1, ccf_results["Correlation"].abs().max()
    )

    for ax, target in zip(axes, targets):

        subset = ccf_results[
            ccf_results["Target"] == target
        ]

        matrix = subset.pivot(
            index="Predictor",
            columns="Lag",
            values="Correlation"
        ).reindex(
            index=list(covariates),
            columns=range(max_lag + 1)
        )

        im = ax.imshow(
            matrix.values,
            cmap="RdBu_r",
            vmin=-max_abs,
            vmax=max_abs,
            aspect="auto"
        )

        ax.set_title(f"Final Demand: {target}")
        ax.set_xlabel("Upstream lead (months)")
        ax.set_xticks(range(max_lag + 1))

        ax.set_yticks(range(len(covariates)))
        ax.set_yticklabels(list(covariates))

    fig.colorbar(
        im,
        ax=axes,
        label="Cross-correlation",
        shrink=0.75
    )

    fig.suptitle(
        "Do Intermediate-Demand Prices Lead Final-Demand Inflation?",
        y=1.03
    )

    plt.show()


ccf_results = compute_upstream_ccf(
    infl,
    list(TARGETS),
    list(COVARIATES),
    MAX_CCF_LAG
)

plot_upstream_ccf(
    ccf_results,
    TARGETS,
    COVARIATES,
    MAX_CCF_LAG
)

Interpretation

  • Goods: Strong contemporaneous correlations with nearly all intermediate-demand categories. Positive correlations persist for several months, particularly for Stage 4 goods, suggesting that upstream inflation may contain useful leading information.

  • Services: Contemporaneous correlations are generally weaker than for goods, but remain positive at several leads. Stage 4 goods inflation exhibits the most persistent positive relationship, suggesting possible predictive information beyond contemporaneous price movements.

  • Construction: Contemporaneous correlations are relatively weak. However, several intermediate-demand measures exhibit stronger positive correlations at longer leads (4–6 months), especially Stage 4 goods and construction inputs.

  • Implication for forecasting: These patterns motivate including historical intermediate-demand inflation as covariates in Chronos-2. However, cross-correlations do not establish incremental predictive value; our out-of-sample comparison will determine whether these variables improve forecasts beyond the information contained in past final-demand inflation.

We also summarize the strongest positive upstream correlation for each target–predictor combination. This is useful for identifying candidate lead-lag relationships, although selecting the largest observed correlation can exaggerate apparent predictability.

# Exclude lag zero to focus on leading relationships
positive_lags = ccf_results[
    ccf_results["Lag"] > 0
].copy()

# Identify the strongest positive correlation
idx = positive_lags.groupby(
    ["Target", "Predictor"]
)["Correlation"].idxmax()

lead_summary = (
    positive_lags.loc[
        idx,
        ["Target", "Predictor", "Lag", "Correlation"]
    ]
    .sort_values(
        ["Target", "Correlation"],
        ascending=[True, False]
    )
    .reset_index(drop=True)
)
Target Predictor Lead (months) Peak correlation
0 Construction Stage4_Goods 6 0.491
1 Construction Stage4_Construction 4 0.368
2 Construction Stage3_Goods 6 0.341
3 Construction Stage3_Services 4 0.316
4 Construction Stage2_Goods 5 0.310
5 Construction Stage2_Services 4 0.306
6 Construction Stage4_Services 6 0.224
7 Construction Stage3_Construction 5 0.176
8 Goods Stage4_Goods 1 0.445
9 Goods Stage2_Goods 1 0.419
10 Goods Stage4_Construction 1 0.417
11 Goods Stage3_Goods 1 0.403
12 Goods Stage3_Services 1 0.317
13 Goods Stage3_Construction 1 0.287
14 Goods Stage4_Services 2 0.247
15 Goods Stage2_Services 5 0.168
16 Services Stage4_Goods 4 0.410
17 Services Stage4_Construction 1 0.403
18 Services Stage3_Services 1 0.284
19 Services Stage4_Services 2 0.284
20 Services Stage3_Construction 1 0.268
21 Services Stage3_Goods 1 0.268
22 Services Stage2_Goods 1 0.220
23 Services Stage2_Services 1 0.217

Interpretation

  • Positive correlations at lags greater than zero are consistent with upstream price changes preceding downstream inflation.
  • Different peak lags across targets may suggest heterogeneous transmission patterns.
  • Weak or negative correlations suggest that an upstream predictor has a less straightforward positive leading relationship with that target.

These patterns motivate including all eight intermediate-demand histories as Chronos-2 covariates. The correlations are descriptive and use the full revised sample; they do not justify selecting a subset of predictors after observing test-period outcomes. Moreover, strong cross-correlations may reflect common shocks or information already present in lagged final-demand inflation. The subsequent rolling-origin comparison with and without upstream covariates tests their incremental forecast value.

4 Evaluation design

At forecast origin \(t\), the data include only monthly changes through month \(t\) (interpreted as available immediately after that month’s PPI release in calendar month \(t+1\)). Every method predicts \(\pi_{t+1},\ldots,\pi_{t+6}\). Because the PPI release occurs after month-end, the first target month may already be partly underway when the forecast is made; this timing convention uses no additional within-month information. The first forecast origin is December 2017, so January 2018 is the first target month. Forecast origins advance every month. Training histories expand rather than maintain a fixed length.

Let \(\widehat\pi_{t+h\mid t}\) be a forecast made at origin \(t\). We evaluate each horizon separately:

\[ \operatorname{RMSE}(h)=\sqrt{\frac{1}{N_h}\sum_t (\pi_{t+h}-\widehat\pi_{t+h\mid t})^2}. \]

We also compare empirical coverage and average width of nominal 80% intervals, along with pinball loss \(L_q(y,\widehat Q_q)=\bigl(q-\mathbf{1}\{y<\widehat Q_q\}\bigr)(y-\widehat Q_q)\) at \(q=0.1,0.5,0.9\). Prediction intervals and quantiles are marginal for each future month, not intervals for the six-month cumulative sum.

first_origin = FIRST_TARGET_DATE - pd.offsets.MonthBegin(1)
origins = infl.index[(infl.index >= first_origin) & (infl.index < infl.index.max())]
if RUN_MODE == "quick":
    origins = origins[-12:]
assert len(origins) >= 1, "Insufficient out-of-sample dates"

MODEL_ORDER = ["ARIMA", "Chronos-Uni", "Chronos-Multi", "Chronos-Multi+ID"]
MODEL_LABELS = {
    "ARIMA": "Automatic ARIMA", "Chronos-Uni": "Chronos-2 univariate",
    "Chronos-Multi": "Chronos-2 multivariate",
    "Chronos-Multi+ID": "Chronos-2 multivariate + upstream inputs",
    "Chronos-DirectTotal": "Chronos-2 direct total final demand",
}


def cache_path(model):
    config = f"{CODE_VERSION}|{SNAPSHOT_HASH}|{RUN_MODE}|{len(origins)}|" \
             f"{origins[0]}|{origins[-1]}|{FORECAST_HORIZON}|{model}|" \
             f"{VERSIONS.get('chronos-forecasting')}|{VERSIONS.get('statsforecast')}"
    digest = hashlib.sha256(config.encode()).hexdigest()[:16]
    return CACHE_DIR / f"forecast_{model}_{digest}.pkl"


def pack_predictions(model, origin, target, mean, lo, hi):
    forecast_dates = pd.date_range(origin + pd.offsets.MonthBegin(1),
                                   periods=FORECAST_HORIZON, freq="MS")
    n = len(forecast_dates)
    assert all(len(np.atleast_1d(v)) == n for v in [mean, lo, hi])
    out = pd.DataFrame({
        "model": model, "origin": origin, "date": forecast_dates,
        "h": np.arange(1, n + 1), "target": target,
        "pred": np.asarray(mean, dtype=float), "q10": np.asarray(lo, dtype=float),
        "q90": np.asarray(hi, dtype=float),
    })
    if not np.isfinite(out[["pred", "q10", "q90"]].values).all():
        raise ValueError(f"Non-finite forecast: {model}, {origin}, {target}")
    if (out.q10 > out.q90).any():
        raise ValueError(f"Crossing 10th and 90th percentiles: {model}, {origin}")
    return out


def run_or_restore(model, forecaster):
    file = cache_path(model)
    if file.exists():
        df = pd.read_pickle(file)
    else:
        parts = []
        for idx, origin in enumerate(origins):
            parts.append(forecaster(origin))
            if (idx + 1) % 12 == 0 or idx == len(origins) - 1:
                print(f"{model}: {idx + 1} / {len(origins)} forecast origins")
        df = pd.concat(parts, ignore_index=True)
        df.to_pickle(file)
    assert set(df["origin"]) == set(origins)
    return df


def attach_truth(df):
    observed = (infl[list(TARGETS) + [AGGREGATE]]
                .stack().rename("actual").rename_axis(["date", "target"]).reset_index())
    return df.merge(observed, on=["date", "target"], how="left", validate="many_to_one")


def evaluate(df):
    a = attach_truth(df).dropna(subset=["actual"]).copy()
    a["squared_error"] = (a.actual - a.pred) ** 2
    a["covered"] = (a.actual >= a.q10) & (a.actual <= a.q90)
    a["width"] = a.q90 - a.q10
    for q, col in [(0.1, "q10"), (0.5, "pred"), (0.9, "q90")]:
        err = a.actual - a[col]
        a[f"pinball_{q}"] = np.maximum(q * err, (q - 1) * err)
    def aggregate_scores(g):
        return pd.Series({
            "N": len(g), "RMSE": np.sqrt(g.squared_error.mean()),
            "Coverage80": g.covered.mean(), "Width80": g.width.mean(),
            "Pinball": g[["pinball_0.1", "pinball_0.5", "pinball_0.9"]].values.mean(),
        })
    scores = a.groupby(["model", "target", "h"], sort=False).apply(
        aggregate_scores, include_groups=False).reset_index()
    return scores, a


def show_one_model_snapshot(preds, model):
    # A historical origin with all six outcomes observed. Used only for *illustration*.
    origin = infl.index[-(FORECAST_HORIZON + 1)]
    one = preds[(preds.model == model) & (preds.origin == origin)]
    if one.empty:
        display(Markdown("No completed six-month path in the selected origins."))
        return
    fig, axes = plt.subplots(1, 3, figsize=(13, 3.6), sharey=False)
    for ax, target in zip(axes, TARGETS):
        seg = one[one.target == target].sort_values("date")
        hist = infl.loc[:origin, target].tail(12)
        ax.plot(hist.index, hist.values, label="History", color="gray", lw=1.5)
        ax.plot(seg.date, seg.pred, marker="o", label="Forecast")
        ax.fill_between(seg.date, seg.q10, seg.q90, alpha=0.18, label="80% interval")
        actual = infl[target].reindex(seg.date)
        ax.plot(seg.date, actual, ls="--", lw=1.6, label="Realized")
        ax.axvline(origin, ls=":", color="gray", lw=0.9)
        ax.set_title(target)
        ax.tick_params(axis="x", rotation=45)
    axes[0].set_ylabel("Monthly log inflation (pp)")
    axes[2].legend(loc="best", fontsize=8)
    fig.suptitle(f"{MODEL_LABELS[model]}: origin {origin:%b %Y}")
    fig.tight_layout()
    plt.show()


def describe_model(scores, model):
    x = scores[(scores.model == model) & scores.target.isin(TARGETS)]
    if x.empty:
        return
    by_target = x.groupby("target")["RMSE"].mean()
    coverage = x.Coverage80.mean()
    display(Markdown(
        f"**Interpretation.** Averaging horizon-specific RMSEs, the smallest error "
        f"for **{MODEL_LABELS[model]}** is for **{by_target.idxmin()}** "
        f"({by_target.min():.3f} percentage points). Mean empirical 80% interval "
        f"coverage across targets and horizons is **{coverage:.1%}**. "
        "Coverage near 80% is desirable, but should be considered together with interval width. "
        "Overlapping six-month forecasts are not independent observations."
    ))

Backtest origins: Dec 2017 through Aug 2026 (105 origins). The last origin has fewer than six realized future months; evaluation excludes unobserved outcomes for each horizon.

5 Automatic ARIMA

An ARIMA model provides a simple benchmark without any pretrained neural network. We let the algorithm select nonseasonal \(p\), \(d\), and \(q\) separately for each target at each forecast origin, using AICc and training data only. Because the input is inflation rather than the PPI level, many selected specifications may use \(d=0\), but we do not impose this.

For the aggregate extension, we also fit one direct ARIMA model of total final-demand inflation. That forecast is not an average of component forecasts. The benchmark is intentionally nonseasonal; the ACF for construction suggests that seasonal alternatives might be worth exploring beyond this demonstration.

def arima_at(origin):
    outputs = []
    for name in [*TARGETS, AGGREGATE]:
        y = infl.loc[:origin, name].to_numpy(dtype=np.float64)
        model = AutoARIMA(season_length=1, seasonal=False, ic="aicc",
                          max_p=3, max_q=3, max_d=1, stepwise=True)
        fit = model.forecast(y=y, h=FORECAST_HORIZON, level=[80])
        outputs.append(pack_predictions(
            "ARIMA", origin, name, fit["mean"], fit["lo-80"], fit["hi-80"]))
    return pd.concat(outputs, ignore_index=True)

arima_preds = run_or_restore("ARIMA", arima_at)
arima_scores, _ = evaluate(arima_preds)

ARIMA out-of-sample RMSE by target and forecast horizon (percentage points, pooled across backtest origins; each origin is fit using data available up to that date only).

target Construction Goods Services
h
1 0.770 0.896 0.365
2 0.767 0.977 0.351
3 0.766 0.995 0.348
4 0.919 1.018 0.364
5 0.916 1.014 0.382
6 0.925 1.023 0.381

Interpretation. Averaging horizon-specific RMSEs, the smallest error for Automatic ARIMA is for Services (0.365 percentage points). Mean empirical 80% interval coverage across targets and horizons is 71.4%. Coverage near 80% is desirable, but should be considered together with interval width. Overlapping six-month forecasts are not independent observations.

Benchmark interpretation. Services inflation is substantially less volatile than goods or construction, so its lower absolute RMSE does not automatically mean it is easier to forecast relative to a scale-adjusted benchmark. We compare models within each target and horizon before averaging across target categories.

6 Chronos-2: zero-shot univariate forecasting

Unlike ARIMA, Chronos-2 arrives with pretrained parameters. Its univariate configuration predicts each target using only that target’s own history. This is the clearest baseline for asking whether a foundation model adds anything relative to a traditional time-series model. We retain the same origin dates, target transformation, forecast horizons, and uncertainty summaries.

Implementation detail: Three rows of item_id values with a single target column represent three independent forecasting tasks when cross_learning=False. This is deliberately different from the multivariate model in the next section.

# One model download and initialization, reused by all Chronos variants.
pipeline = Chronos2Pipeline.from_pretrained("amazon/chronos-2", device_map=DEVICE)


def chronos_univariate_at(origin):
    train = infl.loc[:origin, list(TARGETS)]
    df = (train.rename_axis("timestamp").reset_index()
          .melt(id_vars="timestamp", var_name="item_id", value_name="target")
          .sort_values(["item_id", "timestamp"]).reset_index(drop=True))
    out = pipeline.predict_df(
        df, id_column="item_id", timestamp_column="timestamp", target="target",
        prediction_length=FORECAST_HORIZON, quantile_levels=QUANTILES,
        cross_learning=False, freq="MS")
    pieces = []
    for name in TARGETS:
        row = out[out.item_id == name].sort_values("timestamp")
        assert len(row) == FORECAST_HORIZON
        pieces.append(pack_predictions("Chronos-Uni", origin, name,
                                      row["predictions"].to_numpy(),
                                      row["0.1"].to_numpy(), row["0.9"].to_numpy()))
    return pd.concat(pieces, ignore_index=True)

uni_preds = run_or_restore("Chronos-Uni", chronos_univariate_at)
uni_scores, _ = evaluate(uni_preds)

Univariate Chronos-2 out-of-sample RMSE by target and forecast horizon (percentage points, pooled across backtest origins).

target Construction Goods Services
h
1 0.736 0.914 0.344
2 0.775 0.987 0.351
3 0.794 0.982 0.349
4 0.901 0.980 0.349
5 0.929 0.989 0.361
6 0.965 1.010 0.367

Interpretation. Averaging horizon-specific RMSEs, the smallest error for Chronos-2 univariate is for Services (0.353 percentage points). Mean empirical 80% interval coverage across targets and horizons is 69.4%. Coverage near 80% is desirable, but should be considered together with interval width. Overlapping six-month forecasts are not independent observations.

Against the ARIMA benchmark: Univariate Chronos-2 has lower RMSE in 11 of 18 target–horizon cells. This direct comparison holds both the target and historical information set fixed; zero-shot pretraining does not guarantee lower forecast error in every case.

7 Chronos-2: multivariate forecasting

Goods, services, and construction inflation respond to some of the same macroeconomic forces. A multivariate forecast asks whether the model can exploit these relationships rather than treating the target variables as unrelated series.

At each timestamp we supply one row with three target columns. We give the model one shared item ID and explicitly set target=["Goods", "Services", "Construction"], which invokes its multivariate pathway. This is not merely running three univariate models in a batch.

def chronos_multivariate_at(origin):
    df = infl.loc[:origin, list(TARGETS)].rename_axis("timestamp").reset_index()
    df.insert(0, "item_id", "FD_system")
    out = pipeline.predict_df(
        df, id_column="item_id", timestamp_column="timestamp",
        target=list(TARGETS), prediction_length=FORECAST_HORIZON,
        quantile_levels=QUANTILES, freq="MS")
    pieces = []
    for name in TARGETS:
        row = out[out.target_name == name].sort_values("timestamp")
        assert len(row) == FORECAST_HORIZON
        pieces.append(pack_predictions("Chronos-Multi", origin, name,
                                      row["predictions"].to_numpy(),
                                      row["0.1"].to_numpy(), row["0.9"].to_numpy()))
    return pd.concat(pieces, ignore_index=True)

multi_preds = run_or_restore("Chronos-Multi", chronos_multivariate_at)
multi_scores, _ = evaluate(multi_preds)

Multivariate Chronos-2 out-of-sample RMSE by target and forecast horizon (percentage points, pooled across backtest origins).

target Construction Goods Services
h
1 0.757 0.901 0.341
2 0.762 0.960 0.347
3 0.744 0.955 0.353
4 0.901 0.953 0.359
5 0.905 0.971 0.364
6 0.915 0.992 0.370

Interpretation. Multivariate Chronos-2 has lower RMSE than independent univariate Chronos-2 for 12 of 18 target–horizon combinations. Gains should be evaluated separately for goods, services, and construction: joint modeling may help one component while worsening another. Correlation among series alone does not guarantee a forecasting improvement.

8 Chronos-2: historical intermediate-demand covariates

The intermediate-demand production-flow indexes give the experiment an economic mechanism. Earlier-stage price changes may propagate through production, and they may be informative about future final-demand prices. But today’s upstream indexes do not reveal their values six months from now.

Chronos-2 distinguishes past-only covariates from known-future covariates. Here we intentionally provide all eight intermediate-demand inflation histories in the context DataFrame and no future_df. Thus, only already observed upstream prices can influence the six-month forecast. This is a no-leakage covariate specification, though the overall retrospective exercise retains the separate pretraining and revision concerns mentioned earlier.

def chronos_covariates_at(origin):
    cols = list(TARGETS) + list(COVARIATES)
    df = infl.loc[:origin, cols].rename_axis("timestamp").reset_index()
    df.insert(0, "item_id", "FD_system")
    # The eight non-target columns are inferred as past-only covariates.
    # No future_df is passed, because future intermediate-demand prices are unknown.
    out = pipeline.predict_df(
        df, id_column="item_id", timestamp_column="timestamp",
        target=list(TARGETS), prediction_length=FORECAST_HORIZON,
        quantile_levels=QUANTILES, freq="MS")
    pieces = []
    for name in TARGETS:
        row = out[out.target_name == name].sort_values("timestamp")
        assert len(row) == FORECAST_HORIZON
        pieces.append(pack_predictions("Chronos-Multi+ID", origin, name,
                                      row["predictions"].to_numpy(),
                                      row["0.1"].to_numpy(), row["0.9"].to_numpy()))
    return pd.concat(pieces, ignore_index=True)

cov_preds = run_or_restore("Chronos-Multi+ID", chronos_covariates_at)
cov_scores, _ = evaluate(cov_preds)

Chronos-2 with historical covariates: out-of-sample RMSE by target and forecast horizon (percentage points, pooled across backtest origins).

target Construction Goods Services
h
1 0.761 0.889 0.334
2 0.764 0.940 0.340
3 0.740 0.912 0.346
4 0.903 0.903 0.349
5 0.902 0.918 0.357
6 0.928 0.954 0.364

Interpretation. Historical upstream inflation improves RMSE relative to multivariate Chronos-2 without covariates in 14 of 18 target–horizon combinations. This supports testing the information value of upstream prices, but is not evidence that we have identified a causal pass-through mechanism. The model has no access to future values of these covariates.

9 Consolidated out-of-sample comparison

We now hold the target, origin dates, and horizon fixed while comparing the four approaches. For point forecasting, smaller RMSE is better. For probabilistic forecasting, interval coverage should be reasonably close to 80%, but extremely wide intervals can achieve high coverage without being useful. Lower pinball loss is preferred.

all_preds = pd.concat([arima_preds, uni_preds, multi_preds, cov_preds],
                      ignore_index=True)
all_scores, all_scored_rows = evaluate(all_preds)
target_scores = all_scores[all_scores.target.isin(TARGETS)].copy()

# Equal-weighted across the THREE target series for each model and horizon.
horizon_scores = (target_scores.groupby(["model", "h"], as_index=False)
                  [["RMSE", "Coverage80", "Width80", "Pinball"]].mean())


def plot_rmse_heatmaps(scores):
    fig, axes = plt.subplots(1, 3, figsize=(12, 3.9), sharey=True)
    # Use a common scale across all three panels, so colors are comparable.
    color_min = float(scores["RMSE"].min())
    color_max = float(scores["RMSE"].max())
    for ax, name in zip(axes, TARGETS):
        a = scores[scores.target == name].pivot(index="model", columns="h", values="RMSE")
        a = a.reindex(MODEL_ORDER)
        mat = a.to_numpy()
        im = ax.imshow(mat, aspect="auto", cmap="viridis_r",
                       vmin=color_min, vmax=color_max)
        ax.set_title(name)
        ax.set_xticks(range(FORECAST_HORIZON), range(1, FORECAST_HORIZON + 1))
        ax.set_yticks(range(len(MODEL_ORDER)), [MODEL_LABELS[m] for m in MODEL_ORDER])
        ax.set_xlabel("Forecast horizon (months)")
        for i in range(mat.shape[0]):
            for j in range(mat.shape[1]):
                ax.text(j, i, f"{mat[i,j]:.2f}", ha="center", va="center", fontsize=8,
                        color="black" if mat[i,j] < np.nanmedian(mat) else "white")
    fig.colorbar(im, ax=axes, label="RMSE (percentage points)", fraction=0.02)
    fig.subplots_adjust(left=0.26, right=0.88, bottom=0.18, top=0.86, wspace=0.12)
    plt.show()


def plot_metrics_by_horizon(summaries):
    fig, axes = plt.subplots(1, 2, figsize=(11, 3.6))
    for model in MODEL_ORDER:
        s = summaries[summaries.model == model].sort_values("h")
        axes[0].plot(s.h, s.RMSE, marker="o", label=MODEL_LABELS[model])
        axes[1].plot(s.h, s.Coverage80, marker="o", label=MODEL_LABELS[model])
    axes[0].set(xlabel="Months ahead", ylabel="RMSE (pp)", title="Point accuracy")
    axes[1].axhline(.8, ls="--", color="black", lw=1)
    axes[1].set(xlabel="Months ahead", ylabel="Empirical interval coverage",
                title="Coverage of nominal 80% intervals", ylim=(0, 1))
    axes[1].legend(fontsize=8, bbox_to_anchor=(1.04, 1), loc="upper left")
    fig.tight_layout()
    plt.show()

plot_rmse_heatmaps(target_scores)
plot_metrics_by_horizon(horizon_scores)

Out-of-sample RMSE by model and forecast horizon (percentage points, averaged equally across the three targets and pooled across backtest origins).

model ARIMA Chronos-Uni Chronos-Multi Chronos-Multi+ID
h
1 0.677 0.665 0.666 0.661
2 0.698 0.704 0.690 0.681
3 0.703 0.708 0.684 0.666
4 0.767 0.743 0.738 0.718
5 0.771 0.760 0.747 0.725
6 0.776 0.781 0.759 0.749
Coverage80 Width80 Pinball
model
ARIMA 0.714 1.226 0.175
Chronos-Uni 0.694 1.144 0.160
Chronos-Multi 0.702 1.139 0.158
Chronos-Multi+ID 0.699 1.068 0.155

Interpretation. Averaging the three target-specific RMSEs and six forecast horizons, Chronos-2 multivariate + upstream inputs has the lowest mean RMSE in this retrospective comparison. This is a descriptive ranking, not a significance test; longer horizons have fewer realized observations, and successive forecast errors overlap.

Uncertainty interpretation. All four methods under-cover the nominal 80% interval on average. The lowest mean pinball loss belongs to Chronos-2 multivariate + upstream inputs (0.155); consider this together with coverage and interval width. A narrower interval is not necessarily better if it systematically misses outcomes.

Incremental forecast value. Chronos-2 with upstream covariates improves on the ARIMA benchmark in 17 of 18 matched component–horizon comparisons. This pattern is more informative than comparing raw RMSE across series of very different volatilities, but it does not by itself establish statistical significance.

10 Cumulative inflation and reconstructed price levels

A forecast for month \(t+h\) is not the same as a forecast for inflation accumulated from \(t\) to \(t+h\). Given monthly log inflation forecasts, cumulative log inflation and the implied PPI level are

\[ \widehat\Pi_{t,h}=\sum_{k=1}^{h}\widehat\pi_{t+k\mid t},\qquad \widehat P_{t+h\mid t}=P_t\exp\!\left(\frac{\widehat\Pi_{t,h}}{100}\right). \]

To express the forecast as an ordinary cumulative percentage change, use \(100[\exp(\widehat\Pi_{t,h}/100)-1]\). We compare these point forecasts across horizons. We do not sum marginal prediction-interval endpoints: a valid interval for cumulative inflation requires the joint predictive distribution of future monthly changes (or another defensible approximation).

valid = attach_truth(all_preds)
valid = valid[valid.target.isin(TARGETS)].copy()
valid = valid.sort_values(["model", "target", "origin", "h"])
key = ["model", "target", "origin"]
valid["cum_pred_log"] = valid.groupby(key)["pred"].cumsum()
valid["cum_actual_log"] = valid.groupby(key)["actual"].cumsum()
valid["cum_valid"] = valid.groupby(key)["actual"].transform(
    lambda s: s.notna().cumprod().astype(bool))
valid.loc[~valid.cum_valid, "cum_actual_log"] = np.nan
valid["cum_pred_pct"] = 100 * np.expm1(valid.cum_pred_log / 100)
valid["cum_actual_pct"] = 100 * np.expm1(valid.cum_actual_log / 100)

# Last known index level at the origin converts cumulative predictions into levels.
level_base = (levels[list(TARGETS)].stack().rename("origin_level")
              .rename_axis(["origin", "target"]).reset_index())
valid = valid.merge(level_base, on=["origin", "target"], validate="many_to_one")
valid["implied_index"] = valid.origin_level * np.exp(valid.cum_pred_log / 100)

c = valid.dropna(subset=["cum_actual_pct"]).copy()
c["err2"] = (c.cum_pred_pct - c.cum_actual_pct) ** 2
cumulative_rmse = (c.groupby(["model", "h"])["err2"].mean()
                   .pipe(np.sqrt).unstack("model").reindex(columns=MODEL_ORDER))
model ARIMA Chronos-Uni Chronos-Multi Chronos-Multi+ID
h
1 0.721 0.713 0.714 0.710
2 1.168 1.177 1.144 1.126
3 1.542 1.566 1.491 1.438
4 1.999 1.987 1.881 1.798
5 2.448 2.405 2.273 2.137
6 2.905 2.868 2.692 2.528

At 1 month, the lowest pooled cumulative percentage-change RMSE belongs to Chronos-2 multivariate + upstream inputs (0.710 percentage points).

At 3 months, the lowest pooled cumulative percentage-change RMSE belongs to Chronos-2 multivariate + upstream inputs (1.438 percentage points).

At 6 months, the lowest pooled cumulative percentage-change RMSE belongs to Chronos-2 multivariate + upstream inputs (2.528 percentage points).

Interpretation. Cumulative forecast errors generally increase with the horizon because errors in successive monthly predictions accumulate. The cumulative metric summarizes changes in the price level, not the accuracy of any single future month’s inflation. Marginal monthly prediction intervals cannot simply be added to obtain a cumulative interval.

At the six-month horizon, the best approach reduces cumulative RMSE by approximately 13.0% relative to ARIMA when errors are pooled across the three component series.

11 Forecasting total final-demand inflation: direct and component-based

The three final-demand groups exhaust total final demand. But their index levels cannot simply be added or equally averaged to recover the official aggregate. As a teaching approximation, we combine the component inflation forecasts with contemporaneously available BLS relative-importance weights:

\[ \widehat\pi^{FD}_{t+h\mid t}\ \approx\ w_{G,t}\widehat\pi^G_{t+h\mid t} + w_{S,t}\widehat\pi^S_{t+h\mid t} + w_{C,t}\widehat\pi^C_{t+h\mid t}, \qquad \sum_jw_{j,t}=1. \]

This calculation is approximate: BLS uses its own index aggregation methodology and evolving effective weights; the simple weighted sum of seasonally adjusted component log changes is not an identity. We use one set of weights known at the forecast origin for the whole future six-month path, so no realized future shares enter the forecast.

11.1 Historical BLS weights, with a transparent fallback

BLS posts downloadable annual FDGROUPREL historical relative-importance spreadsheets. The following optional routine looks for the three exact broad categories in those Excel files; it accepts a year’s weights only if a common numeric column across all three categories sums to approximately 100%. Since the layout may vary across archived files, the code refuses to guess when the parser cannot identify the relevant column.

To avoid assuming that year-end weights were already published on December 31, we apply a conservative reporting lag: an annual archive from year \(y\) is first eligible for forecast origins April of year \(y+1\). If a suitable archived file is unavailable, we use the fixed weights from December 2016, publicly reported in February 2017: goods 33.837%, services 64.086%, construction 2.077%. This fallback was known before our December 2017 first origin. The rendered document shows exactly how often each source was used.

WEIGHT_ARCHIVE = DATA_DIR / "bls_annual_weights.csv"
FALLBACK_WEIGHTS = {"Goods": 0.33837, "Services": 0.64086, "Construction": 0.02077}
LABEL_MAP = {"final demand goods": "Goods", "final demand services": "Services",
             "final demand construction": "Construction"}


def norm_label(value):
    x = re.sub(r"\s+", " ", str(value).lower().replace("-", " ")).strip()
    return x


def parse_weight_book(content, archive_year):
    sheets = pd.read_excel(io.BytesIO(content), sheet_name=None, header=None)
    for sheet in sheets.values():
        matches = {}
        for row_i in range(len(sheet)):
            cells = [norm_label(v) for v in sheet.iloc[row_i].to_list()]
            for exact_label, name in LABEL_MAP.items():
                if exact_label in cells:
                    matches[name] = row_i
        if len(matches) != 3:
            continue
        candidates = []
        for column in sheet.columns:
            vals = pd.to_numeric(pd.Series(
                [sheet.iat[matches[name], column] for name in TARGETS]),
                errors="coerce").to_numpy(dtype=float)
            if np.isfinite(vals).all() and (vals > 0).all() and abs(vals.sum() - 100) < 1.5:
                candidates.append((int(column), vals / vals.sum()))
        if candidates:
            # Prefer the rightmost valid percentage column for annual month-end figures.
            _, vals = max(candidates, key=lambda p: p[0])
            return {"year": archive_year,
                    **{name: float(v) for name, v in zip(TARGETS, vals)}}
    return None


def build_weight_archive():
    rows = []
    for year in range(2013, dt.date.today().year + 1):
        url = ("https://www.bls.gov/ppi/tables/"
               f"final-demand-relative-importance-by-component-series-{year}.xlsx")
        try:
            r = requests.get(url, timeout=20, headers={"User-Agent": "Mozilla/5.0"})
            r.raise_for_status()
            record = parse_weight_book(r.content, year)
            if record is not None:
                rows.append(record)
        except (requests.RequestException, ValueError, OSError):
            pass  # The fallback is explicit; one broken BLS archive is not fatal.
    if rows:
        pd.DataFrame(rows).drop_duplicates("year").to_csv(WEIGHT_ARCHIVE, index=False)
    elif not WEIGHT_ARCHIVE.exists():
        pd.DataFrame(columns=["year", *TARGETS]).to_csv(WEIGHT_ARCHIVE, index=False)

if REFRESH_WEIGHTS or not WEIGHT_ARCHIVE.exists():
    build_weight_archive()
annual_weights = pd.read_csv(WEIGHT_ARCHIVE)


def origin_weights(origin):
    # From April of a given year, the PREVIOUS year archive is eligible;
    # January-March use the year before that. The rule is deliberately conservative.
    eligible = origin.year - (1 if origin.month >= 4 else 2)
    usable = annual_weights[annual_weights.year <= eligible] if not annual_weights.empty else pd.DataFrame()
    if len(usable):
        row = usable.sort_values("year").iloc[-1]
        weights = {k: float(row[k]) for k in TARGETS}
        return weights, f"BLS archive {int(row.year)}"
    return FALLBACK_WEIGHTS.copy(), "Fixed Dec 2016"

weight_sources = pd.DataFrame([{"origin": o, "source": origin_weights(o)[1]}
                               for o in origins])
Forecast origins
source
Fixed Dec 2016 105
Component Weight
0 Goods 0.3384
1 Services 0.6409
2 Construction 0.0208

Weight interpretation. The December 2016 fixed-weight fallback is used for 105 of 105 forecast origins; the final origin uses Fixed Dec 2016. When fixed weights are prevalent, differences between reconstructed and published total inflation reflect both the quality of component forecasts and aggregation approximation. Do not treat the reconstructed series as the exact BLS aggregate.

11.2 Chronos-2: direct forecast of total final-demand PPI

So far, Chronos-2 forecasts the three components, which we combine using BLS weights. A fairer comparison with ARIMA fitted directly to the published total also requires forecasting that total directly with Chronos-2. Otherwise, any accuracy difference could reflect the direct versus component-based forecasting strategy, not just the choice of model.

We now give Chronos-2 the historical monthly log inflation series for total final-demand PPI (PPIFIS) as a single univariate target. We retain the same zero-shot model, expanding history, monthly forecast origins, six horizons, and 10th/50th/90th quantiles used elsewhere. The future aggregate and component values are never included as inputs.

This additional specification has its own forecast cache, so introducing it does not force the original three Chronos-2 configurations to rerun.

def chronos_direct_total_at(origin):
    # One historical series, not a multivariate set or a weighted component sum.
    df = (infl.loc[:origin, [AGGREGATE]]
          .rename_axis("timestamp").reset_index()
          .rename(columns={AGGREGATE: "target"}))
    df.insert(0, "item_id", AGGREGATE)

    out = pipeline.predict_df(
        df,
        id_column="item_id",
        timestamp_column="timestamp",
        target="target",
        prediction_length=FORECAST_HORIZON,
        quantile_levels=QUANTILES,
        cross_learning=False,
        freq="MS",
    )
    row = out[out.item_id == AGGREGATE].sort_values("timestamp")
    assert len(row) == FORECAST_HORIZON
    return pack_predictions(
        "Chronos-DirectTotal", origin, AGGREGATE,
        row["predictions"].to_numpy(),
        row["0.1"].to_numpy(), row["0.9"].to_numpy(),
    )

# A separate cache key retains all previously computed component forecasts.
direct_chronos_preds = run_or_restore("Chronos-DirectTotal", chronos_direct_total_at)

direct_arima_preds = arima_preds[arima_preds.target == AGGREGATE].copy()
direct_total_preds = pd.concat(
    [direct_arima_preds, direct_chronos_preds], ignore_index=True
)
direct_total_scores, direct_total_scored = evaluate(direct_total_preds)

Both direct models produce marginal 80% prediction intervals for total inflation; no BLS aggregation approximation is involved. This supports a probabilistic comparison that is not available from simply averaging the components’ marginal quantiles.

Direct aggregate forecasts: point accuracy and uncertainty by horizon.

model h N RMSE Coverage80 Width80 Pinball
0 ARIMA 1 105.0 0.396 0.657 0.754 0.101
1 Chronos-DirectTotal 1 105.0 0.392 0.695 0.774 0.099
2 ARIMA 2 104.0 0.421 0.663 0.780 0.107
3 Chronos-DirectTotal 2 104.0 0.406 0.692 0.809 0.104
4 ARIMA 3 103.0 0.441 0.650 0.798 0.113
5 Chronos-DirectTotal 3 103.0 0.419 0.689 0.824 0.106
6 ARIMA 4 102.0 0.448 0.637 0.814 0.115
7 Chronos-DirectTotal 4 102.0 0.425 0.667 0.839 0.109
8 ARIMA 5 101.0 0.468 0.634 0.829 0.121
9 Chronos-DirectTotal 5 101.0 0.439 0.673 0.842 0.113
10 ARIMA 6 100.0 0.480 0.640 0.842 0.125
11 Chronos-DirectTotal 6 100.0 0.447 0.660 0.854 0.115

Interpretation. Direct Chronos-2 has lower RMSE than direct ARIMA at 6 of 6 horizons. Across horizons, average 80% interval coverage is 67.9% for Chronos-2 and 64.7% for ARIMA. Thus, point accuracy and interval calibration must be assessed separately; these forecasts do not depend on approximating component weights.

11.3 Compare component-based and direct aggregate forecasts

We compare two forecasting strategies: (i) forecast goods, services, and construction inflation, then approximate total inflation with BLS weights; or (ii) forecast published total final-demand PPI directly, using ARIMA or Chronos-2. All approaches are evaluated against the same realized total (PPIFIS) at the same origins and horizons.

The reconstructed component forecasts have point predictions only because weighted averages of marginal quantiles are not valid aggregate quantiles. The two direct models, however, have 80% intervals that we can evaluate separately.

component_preds = all_preds[all_preds.target.isin(TARGETS)].copy()
weighted = []
for (model, origin, h), grp in component_preds.groupby(["model", "origin", "h"]):
    w, source = origin_weights(origin)
    assert set(grp.target) == set(TARGETS)
    pred = sum(float(grp.loc[grp.target == k, "pred"].iloc[0]) * w[k]
               for k in TARGETS)
    d = grp.date.iloc[0]
    weighted.append({"model": model + " components", "origin": origin,
                     "date": d, "h": h, "pred": pred, "weight_source": source})
agg_components = pd.DataFrame(weighted)
# Compare direct forecasts with the same point-only fields used for components.
agg_direct = direct_total_preds[["model", "origin", "date", "h", "pred"]].copy()
agg_direct["model"] = agg_direct["model"].map({
    "ARIMA": "ARIMA direct total",
    "Chronos-DirectTotal": "Chronos-2 direct total",
})
assert agg_direct["model"].notna().all()
agg_direct["weight_source"] = "Not applicable"
agg_preds = pd.concat([agg_components, agg_direct], ignore_index=True)
agg_truth = infl[AGGREGATE].rename("actual").rename_axis("date")
agg_preds = agg_preds.join(agg_truth, on="date").dropna(subset=["actual"])
agg_preds["sqerr"] = (agg_preds.pred - agg_preds.actual) ** 2

# Every approach faces exactly the same observable origin/horizon outcomes.
AGG_MODEL_ORDER = [
    "ARIMA direct total", "Chronos-2 direct total",
    "ARIMA components", "Chronos-Uni components",
    "Chronos-Multi components", "Chronos-Multi+ID components",
]
agg_rmse = (agg_preds.groupby(["model", "h"])["sqerr"].mean()
            .pipe(np.sqrt).unstack("model").reindex(columns=AGG_MODEL_ORDER))

# Separate the direct methods' interval metrics from component reconstruction.
direct_interval_scores = direct_total_scores[
    ["model", "h", "N", "Coverage80", "Width80", "Pinball"]
].copy()
direct_interval_scores["model"] = direct_interval_scores["model"].map({
    "ARIMA": "ARIMA direct total",
    "Chronos-DirectTotal": "Chronos-2 direct total",
})

Monthly aggregate-inflation RMSE by forecast horizon. Direct models are shown first, followed by component-based forecasts using the same BLS weight approximation. Lower is better.

model ARIMA direct total Chronos-2 direct total ARIMA components Chronos-Uni components Chronos-Multi components Chronos-Multi+ID components
h
1 0.396 0.392 0.399 0.389 0.396 0.384
2 0.421 0.406 0.429 0.421 0.421 0.406
3 0.441 0.419 0.447 0.437 0.435 0.412
4 0.448 0.425 0.463 0.440 0.443 0.414
5 0.468 0.439 0.467 0.452 0.458 0.428
6 0.480 0.447 0.472 0.464 0.468 0.447

Probabilistic comparison, direct forecasts only. Empirical coverage can be compared with the nominal 80% level; narrower intervals and lower pinball loss are generally preferable, conditional on suitable calibration.

model h N Coverage80 Width80 Pinball
0 ARIMA direct total 1 105.0 0.657 0.754 0.101
1 Chronos-2 direct total 1 105.0 0.695 0.774 0.099
2 ARIMA direct total 2 104.0 0.663 0.780 0.107
3 Chronos-2 direct total 2 104.0 0.692 0.809 0.104
4 ARIMA direct total 3 103.0 0.650 0.798 0.113
5 Chronos-2 direct total 3 103.0 0.689 0.824 0.106
6 ARIMA direct total 4 102.0 0.637 0.814 0.115
7 Chronos-2 direct total 4 102.0 0.667 0.839 0.109
8 ARIMA direct total 5 101.0 0.634 0.829 0.121
9 Chronos-2 direct total 5 101.0 0.673 0.842 0.113
10 ARIMA direct total 6 100.0 0.640 0.842 0.125
11 Chronos-2 direct total 6 100.0 0.660 0.854 0.115

Interpretation. The lowest aggregate RMSE is achieved most often by Chronos-Multi+ID components (6 of 6 horizons). Compare ARIMA direct total with Chronos-2 direct total to isolate model-class differences; compare each direct approach with its component-based alternative to study the forecasting strategy. Caution: component aggregation uses fixed December 2016 weights for 105 of 105 origins. These are approximations, not exact BLS-consistent aggregate forecasts. Small RMSE differences should not be over-interpreted.

12 Post-release-origin robustness check

Chronos-2 was released in October 2025. Forecasts generated using October 2025 observations, released in November 2025, therefore have target months after the model’s release, when that future inflation could not have been in the original pretrained model. For this check, keep only forecast origins of October 2025 or later. The same monthly vintage-revision caveat remains. Because the dataset ends at a recent month, especially large horizons have very few realized observations and any rankings are unstable.

post_preds = all_preds[all_preds.origin >= POST_RELEASE_ORIGIN].copy()
if len(post_preds):
    post_scores, post_observations = evaluate(post_preds)
    post_summary = (post_scores[post_scores.target.isin(TARGETS)]
                    .groupby(["model", "h"], as_index=False)
                    [["RMSE", "Coverage80", "Width80", "Pinball", "N"]].mean())
else:
    post_summary = pd.DataFrame(columns=["model", "h", "RMSE", "N"])
model ARIMA Chronos-Uni Chronos-Multi Chronos-Multi+ID
h
1 0.786 0.762 0.772 0.783
2 0.865 0.860 0.844 0.852
3 0.897 0.860 0.859 0.836
4 0.849 0.796 0.784 0.756
5 0.883 0.831 0.817 0.809
6 0.907 0.865 0.872 0.875
model ARIMA Chronos-Uni Chronos-Multi Chronos-Multi+ID
h
1 11 11 11 11
2 10 10 10 10
3 9 9 9 9
4 8 8 8 8
5 7 7 7 7
6 6 6 6 6

Interpretation. The sample contains only 11 completed one-month-ahead outcomes and 6 completed six-month-ahead outcomes per target. Apparent rankings may reverse between the full retrospective sample and this much shorter, more temporally credible check. Post-release timing does not eliminate historical data revisions or provide a formal test of predictive superiority.

Post-release: direct total-PPI forecast comparison.

model h N RMSE Coverage80 Pinball
0 ARIMA 1 11.0 0.441 0.727 0.115
1 Chronos-DirectTotal 1 11.0 0.469 0.636 0.128
2 ARIMA 2 10.0 0.504 0.600 0.134
3 Chronos-DirectTotal 2 10.0 0.537 0.600 0.151
4 ARIMA 3 9.0 0.512 0.444 0.141
5 Chronos-DirectTotal 3 9.0 0.545 0.556 0.151
6 ARIMA 4 8.0 0.478 0.500 0.129
7 Chronos-DirectTotal 4 8.0 0.509 0.500 0.143
8 ARIMA 5 7.0 0.489 0.429 0.130
9 Chronos-DirectTotal 5 7.0 0.528 0.571 0.147
10 ARIMA 6 6.0 0.524 0.667 0.142
11 Chronos-DirectTotal 6 6.0 0.568 0.667 0.160

Direct-total post-release interpretation. ARIMA has lower RMSE than direct Chronos-2 at 6 of 6 horizons in this restricted period. This contrast with the full retrospective comparison should be treated cautiously: forecast origins overlap, the sample is very small, and later horizons have progressively fewer observed outcomes.

13 Current six-month forecast

To end on a practical forecasting task, re-run the original four specifications plus direct Chronos-2 total PPI using the latest month in the frozen snapshot as the origin. These predictions are genuinely unobserved within that snapshot and illustrate what the models forecast next. They are not included in the retrospective performance tables. The current origin changes only when the instructor explicitly refreshes the FRED snapshot.

latest = infl.index.max()
# The functions are the same ones used in the backtest.
current = pd.concat([
    arima_at(latest), chronos_univariate_at(latest),
    chronos_multivariate_at(latest), chronos_covariates_at(latest),
    chronos_direct_total_at(latest),
], ignore_index=True)


def plot_current(target):
    fig, ax = plt.subplots(figsize=(10, 3.8))
    hist = infl[target].tail(14)
    ax.plot(hist.index, hist.values, color="gray", lw=2, label="Observed")
    for model in MODEL_ORDER:
        f = current[(current.model == model) & (current.target == target)]
        ax.plot(f.date, f.pred, marker="o", label=MODEL_LABELS[model])
    f = current[(current.model == "Chronos-Multi+ID") & (current.target == target)]
    ax.fill_between(f.date, f.q10, f.q90, alpha=0.16,
                    label="Chronos + ID: 80% interval")
    ax.axvline(latest, ls="--", color="gray", lw=1)
    ax.set(title=f"{target}: six-month forecasts after {latest:%B %Y}",
           ylabel="Monthly log inflation (pp)")
    ax.legend(fontsize=8, ncol=2)
    fig.tight_layout()
    plt.show()


def plot_current_direct_total():
    fig, ax = plt.subplots(figsize=(10, 3.8))
    hist = infl[AGGREGATE].tail(14)
    ax.plot(hist.index, hist.values, color="gray", lw=2, label="Observed total PPI")
    for model, label in [("ARIMA", "ARIMA direct total"),
                         ("Chronos-DirectTotal", "Chronos-2 direct total")]:
        f = current[(current.model == model) & (current.target == AGGREGATE)]
        assert len(f) == FORECAST_HORIZON
        ax.plot(f.date, f.pred, marker="o", label=label)
        ax.fill_between(f.date, f.q10, f.q90, alpha=.13,
                        label=f"{label}: 80% interval")
    ax.axvline(latest, ls="--", color="gray", lw=1)
    ax.set(title=f"Published total PPI: direct forecasts after {latest:%B %Y}",
           ylabel="Monthly log inflation (pp)")
    ax.legend(fontsize=8, ncol=2)
    fig.tight_layout()
    plt.show()

Interpretation. Forecast paths above begin in October 2026. Differences among Chronos-2 variants show the effect of changing the information set, while comparison with ARIMA shows a more fundamental model-class difference. Current predictions should not be described as accurate or inaccurate until new PPI observations become available.

14 Conclusions and limitations

The progression isolates three questions: whether zero-shot forecasting competes with a conventional ARIMA benchmark, whether joint modeling of related final-demand series improves accuracy, and whether historical upstream prices improve forecasts beyond joint target information. Direct versus component-based forecasting is a separate modeling decision, so the aggregate section now compares ARIMA and Chronos-2 applied directly to total PPI, alongside the component-based alternatives. Evaluation of prediction intervals adds the distinct question of whether the model communicates uncertainty usefully.

The results should be read as evidence about this particular PPI forecasting exercise, not as a general ranking of foundation models versus econometric methods. In the full retrospective comparison, joint Chronos-2 forecasts with upstream information are often competitive for the components and their weighted aggregate; direct forecasting of published total PPI provides a more controlled head-to-head comparison with ARIMA. The shorter post-release period can produce a different ranking and has limited statistical power.

Interpretation boundaries:

  • Pretraining contamination: The 2018-onward retrospective window may overlap Chronos-2 pretraining data. The post-release-origin check is cleaner temporally, but very small.
  • Data vintages: FRED’s latest revised PPI data are not identical to what an analyst would have seen at every historical forecast origin.
  • Causal interpretation: Past intermediate-demand prices can be predictive without establishing a causal pass-through coefficient.
  • Probabilistic aggregation: Component marginal quantiles cannot be combined into a valid aggregate or six-month cumulative prediction interval without joint dependence information.
  • Aggregation weights: A contemporaneously weighted log-inflation sum approximates, but does not exactly replicate, BLS index construction. Archive parsing failures trigger a documented pre-2018 fixed-weight fallback.
  • Point summaries: ARIMA point forecasts approximate conditional means, whereas Chronos-2’s predictions output is the median (0.5 quantile). RMSE compares the methods’ supplied point forecasts, not necessarily the same conditional functional; quantile loss provides a complementary probabilistic comparison.
  • Benchmark scope: The automatic ARIMA benchmark is nonseasonal and deliberately parsimonious. The marked three-month autocorrelation in construction makes seasonal specifications a plausible robustness extension.
  • Inference: Monthly rolling origins produce overlapping forecast errors. Differences in RMSE and pinball loss are descriptive unless evaluated with an appropriate formal comparison.

Summary of the rendered experiment. Chronos-2 using upstream inputs beats component ARIMA in 17 of 18 matched comparisons. When both models forecast published total final-demand PPI directly, Chronos-2 has lower full-period RMSE in 6 of 6 horizons. In the shorter post-release sample, direct ARIMA instead has lower RMSE in 6 of 6 horizons. Together, these findings illustrate why both information sets and the evaluation period matter, and why the simplest model remains an essential benchmark.

14.1 Sources and implementation references