# Run cells from top to bottom: later cells reuse the data, functions, and fitted
# models created here. pandas/NumPy handle monthly tables and arrays, sklearn
# and XGBoost fit models, matplotlib draws charts, and joblib saves estimators.
from pathlib import Path
import json, platform, warnings
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import joblib, sklearn, xgboost
from sklearn.exceptions import ConvergenceWarning
from sklearn.base import BaseEstimator, RegressorMixin
from sklearn.linear_model import LinearRegression, Ridge, HuberRegressor
from sklearn.ensemble import RandomForestRegressor, VotingRegressor
from sklearn.neural_network import MLPRegressor
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler, PolynomialFeatures, SplineTransformer
from sklearn.compose import TransformedTargetRegressor
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
from sklearn.model_selection import GridSearchCV, TimeSeriesSplit
from xgboost import XGBRegressor
from IPython.display import display

# Only input data are external. No local Python module is imported.
# Support execution from either the product folder or the commodities workspace.
# Fail early if the saved EIA observation CSV cannot be found.
HERE=Path.cwd()
if not (HERE/'data/eia_observations.csv').exists():
    HERE=HERE/'oil'/'us_snd_crude'
assert (HERE/'data/eia_observations.csv').exists(), 'Keep the data folder beside this notebook.'
OUTPUT=HERE/'model_output'
# These settings affect chart/table presentation only, not stored model values.
plt.style.use('seaborn-v0_8-whitegrid')
pd.set_option('display.max_columns',20)
pd.set_option('display.float_format',lambda x:f'{x:,.2f}')

# Dates represent observation months using their first day. Both exclusion
# endpoints are inclusive; keep the raw calendar intact so lags retain their meaning.
COVID_START = pd.Timestamp('2020-03-01')
COVID_END = pd.Timestamp('2021-03-01')


# Return one Boolean per date: True means outside the excluded COVID interval.
def outside_covid(months):
    dates = pd.DatetimeIndex(months)
    return ~((dates >= COVID_START) & (dates <= COVID_END))


# Filter fitting rows after constructing lagged features on the intact calendar.
# A clean target is insufficient if one of its direct lag inputs is excluded.
def training_rows(frame):
    """Exclude COVID targets and any training row whose direct inputs touch COVID."""
    dates = pd.DatetimeIndex(frame.month)
    keep = outside_covid(dates)
    # lag1 stocks/flows, lag2 for the latest change, lag12 for annual stocks.
    for lag in [1, 2, 12]:
        keep &= outside_covid(dates - pd.DateOffset(months=lag))
    return frame.loc[keep]


# Create chronological expanding-window folds; never shuffle time-series rows.
# This filter concerns target dates; fitting applies the stricter lag-input filter.
def validation_splits(frame, folds=10, test_size=None):
    """Split eligible target dates, retaining real calendar positions and lag dates."""
    allowed = np.flatnonzero(outside_covid(frame.month))
    splitter = TimeSeriesSplit(n_splits=folds, test_size=test_size)
    # Map positions within eligible dates back to the original frame. A test block
    # of 12 eligible observations may span more than 12 calendar months.
    return [(allowed[tr], allowed[te]) for tr, te in splitter.split(allowed)]


# Search backward in annual steps for six eligible year-long development paths.
# Each origin leaves its next 12 forecast months inside development history.
def recursive_origins(development_end, count=6):
    """Use six full test years that do not overlap the excluded interval."""
    origins = []
    cutoff = development_end - pd.DateOffset(years=1)
    while len(origins) < count:
        months = pd.date_range(cutoff + pd.offsets.MonthBegin(), periods=12, freq='MS')
        if outside_covid(months).all() and outside_covid([cutoff]).all():
            origins.append(cutoff)
        cutoff -= pd.DateOffset(years=1)
    return sorted(origins)

# PADDs are the five Petroleum Administration for Defense Districts. Series
# templates substitute {p} with the district number when identifying source data.
PADD_NAMES = {1: 'East Coast', 2: 'Midwest', 3: 'Gulf Coast', 4: 'Rocky Mountain', 5: 'West Coast'}
CODES = {'production_kbd': 'MCRFPP{p}1', 'demand_kbd': 'MCRRIP{p}1',
         'imports_kbd': 'MCRIMP{p}1', 'exports_kbd': 'MCREXP{p}1',
         'net_receipts_kbd': 'MCRNRP{p}1', 'adjustments_kbd': 'MCRUA_R{p}0_1',
         'transfers_kbd': 'M_EPC0_TVP_R{p}0_MBBL', 'direct_use_kbd': 'MCRUPP{p}1',
         'stock_kb': 'MCESTP{p}1', 'cdu_capacity_kbd': 'MOCLEP{p}2'}
# Keep this order aligned with SIGNS: the dot product uses +1 for supply
# additions and -1 for disposition to reconstruct net supply.
FLOWS = ['production_kbd', 'demand_kbd', 'imports_kbd', 'exports_kbd',
         'net_receipts_kbd', 'adjustments_kbd', 'transfers_kbd', 'direct_use_kbd']
SIGNS = [1, -1, 1, -1, 1, 1, 1, -1]
# The downloaded transfer series start in January 2022; earlier rows are structural zeros.
SPARSE = {'exports_kbd', 'transfers_kbd', 'net_receipts_kbd'}

# Return one row per PADD/month with flows, month-end stocks, assumption flags,
# and balance diagnostics, starting from the saved long-form EIA observations.
def load_panel(data_dir=HERE / 'data', start='2007-01-01'):
    # Parse timestamps for exact calendar reindexing and month-length calculations.
    raw = pd.read_csv(Path(data_dir) / 'eia_observations.csv', parse_dates=['month'])
    # Raw flow snapshots retain EIA monthly volumes; convert exactly once here.
    raw['component'] = raw.component.replace({c.replace('_kbd', '_kb'): c for c in FLOWS})
    raw = raw[raw.component.isin(FLOWS + ['commercial_stock_kb', 'cdu_capacity_kbd'])].copy()
    raw['component'] = raw.component.replace({'commercial_stock_kb': 'stock_kb'})
    raw['value'] = raw.value.astype(float)
    is_flow = raw.component.isin(FLOWS)
    raw.loc[is_flow, 'value'] /= raw.loc[is_flow, 'month'].dt.days_in_month
    # The pivot requires a single observation per month/PADD/component combination.
    if raw.duplicated(['month', 'padd', 'component']).any():
        raise ValueError('Duplicate EIA observations')
    panels = []
    for p in PADD_NAMES:
        r = raw[raw.padd.eq(p)]
        # Use the latest observed stock month as this district's provisional endpoint.
        end = r[r.component.eq('stock_kb') & r.value.notna()].month.max()
        index = pd.date_range(start, end, freq='MS')
        # Turn component names into columns, then expose absent months on a continuous
        # month-start calendar. Reindexing does not interpolate missing observations.
        wide = r.pivot(index='month', columns='component', values='value').reindex(index)
        for c in SPARSE:
            rows = r[r.component.eq(c)].set_index('month')
            # Missing calendar rows in sparse export/transfer sheets are an explicit
            # no-reported-flow assumption, not interpolation from future values.
            absent = ~wide.index.isin(rows.index)
            # Combine source zero flags with wholly absent sparse-series rows. An explicit
            # unavailable value on an existing row remains missing instead of becoming zero.
            wide[c + '_assumed_zero'] = rows.assumed_zero.reindex(index).astype("boolean").fillna(False).astype(bool) | absent
            wide.loc[absent, c] = 0.0
        wide['padd'] = p
        wide['padd_name'] = PADD_NAMES[p]
        wide.index.name = 'month'
        panels.append(wide.reset_index())
    panel = pd.concat(panels, ignore_index=True)
    core = FLOWS + ['stock_kb', 'cdu_capacity_kbd']
    # Find the last month complete in all five districts, then trim to a shared
    # endpoint so national totals compare the same observation months.
    valid = panel[core].notna().all(axis=1)
    common = panel[valid].groupby('month').padd.nunique()
    end = common[common.eq(5)].index.max()
    panel = panel[panel.month.le(end)].copy()
    # Earlier unresolved gaps still fail validation: trimming the endpoint does not
    # authorize filling withheld values or dropping an incomplete district.
    if panel[core].isna().any().any():
        raise ValueError('Unresolved missing core data; inspect the source snapshots')
    if not panel.groupby('padd').size().eq(len(pd.date_range(start, end, freq='MS'))).all():
        raise ValueError('Incomplete PADD calendar')
    if panel.cdu_capacity_kbd.le(0).any():
        raise ValueError('CDU capacity must be positive')
    panel['utilization_ratio'] = panel.demand_kbd / panel.cdu_capacity_kbd
    # Shift within each district so prior stock never comes from another PADD.
    # The first month has no predecessor and therefore a missing stock change.
    panel['previous_stock_kb'] = panel.groupby('padd').stock_kb.shift(1)
    # Sum signed daily flows. Compare this with stock change divided by month length
    # so both sides of the crude reconciliation are in thousand barrels per day.
    panel['balance_kbd'] = panel[FLOWS].to_numpy() @ SIGNS
    panel['stock_change_kb'] = panel.stock_kb - panel.previous_stock_kb
    # Despite its legacy name, this diagnostic uses four flows: production plus
    # imports minus refinery demand and exports. The remaining terms are omitted.
    panel['five_term_balance_kbd'] = panel.production_kbd + panel.imports_kbd - panel.demand_kbd - panel.exports_kbd
    panel['five_term_residual_kbd'] = panel.stock_change_kb / panel.month.dt.days_in_month - panel.five_term_balance_kbd
    panel['accounting_residual_kbd'] = panel.stock_change_kb / panel.month.dt.days_in_month - panel.balance_kbd
    return panel

# Audit the loaded snapshot: show Gulf Coast source identifiers, each PADD's
# calendar coverage, counts of assumed-zero observations, and latest refinery
# input/capacity. Flags make missing-data assumptions visible before fitting.

panel=load_panel(HERE/'data')
manifest=pd.read_json(HERE/'data/source_manifest.json')
display(manifest[manifest.padd.eq(3)][['component','series_id','title']])
display(panel.groupby('padd').agg(start=('month','min'),end=('month','max'),months=('month','size')))
flags=[c for c in panel if c.endswith('_assumed_zero')]
display(panel.groupby('padd')[flags].sum())
display(panel.sort_values('month').groupby('padd').tail(1)[
    ['padd','padd_name','month','cdu_capacity_kbd','demand_kbd','utilization_ratio']].sort_values('padd'))

# Plot commercial inventories, excluding SPR, at national and regional scales.
# Stocks are stored in thousand barrels; dividing by 1000 gives million barrels.

national_history=panel.groupby('month')[['stock_kb']].sum()
fig,axes=plt.subplots(1,2,figsize=(14,4))
(national_history/1000).rename(columns={'stock_kb':'Commercial crude (target)'}).plot(ax=axes[0])
for p,g in panel.groupby('padd'):
    axes[1].plot(g.month,g.stock_kb/1000,label=f'PADD {p}')
axes[0].set_title('Commercial crude inventory target');axes[1].set_title('Commercial crude by PADD');axes[1].legend()
for ax in axes:ax.set_ylabel('Million barrels')
plt.tight_layout();plt.show()

# Compare residuals from the four-flow shortcut and the full published-flow
# balance. The >2 kbd exception threshold is a diagnostic display rule, not a
# training filter or a rule for forcing the commercial residual to zero.

audit=panel.groupby('padd').agg(
    incomplete_balance_mae_kbd=('five_term_residual_kbd',lambda s:s.abs().mean()),
    complete_balance_mae_kbd=('accounting_residual_kbd',lambda s:s.abs().mean()),
    max_complete_residual_kbd=('accounting_residual_kbd',lambda s:s.abs().max()))
display(audit)
exceptions=panel[panel.accounting_residual_kbd.abs()>2]
display(exceptions[['month','padd','accounting_residual_kbd']])
fig,axes=plt.subplots(1,2,figsize=(14,3))
for p,g in panel.groupby('padd'):
    axes[0].plot(g.month,g.five_term_residual_kbd,label=f'PADD {p}',alpha=.7)
    axes[1].plot(g.month,g.accounting_residual_kbd,label=f'PADD {p}',alpha=.7)
axes[0].set(title='Four-flow commercial balance residual',ylabel='Thousand barrels/day')
axes[1].set(title='Published-flow commercial balance residual',ylabel='Thousand barrels/day')
axes[0].legend(ncol=3,fontsize=8);plt.tight_layout();plt.show()

# Build a trend in years plus 11 calendar-month indicators. January is the
# reference month represented by the regression intercept.
def seasonal_design(months, origin):
    dates = pd.DatetimeIndex(months)
    trend = ((dates.year - origin.year) * 12 + dates.month - origin.month).to_numpy() / 12
    return np.column_stack([trend] + [(dates.month == m).astype(float) for m in range(2, 13)])


# Use only the supplied origin history. The last 60 non-COVID observations may
# span more than 60 calendar months because excluded months are skipped.
def forecast_flows(history, months):
    """Forecast daily flows and crude-input/capacity ratio from origin history."""
    train = history.loc[outside_covid(history.month)].tail(60)
    # For crude, model refinery utilization instead of demand directly, then recover
    # refinery input using utilization times capacity known at the forecast origin.
    targets = [c for c in FLOWS if c != 'demand_kbd'] + ['utilization_ratio']
    model = LinearRegression().fit(
        seasonal_design(train.month, train.month.iloc[0]), train[targets])
    result = pd.DataFrame(model.predict(seasonal_design(months, train.month.iloc[0])),
                          columns=targets, index=pd.DatetimeIndex(months))
    # Clip forecast utilization to 0-100% and freeze capacity at its last known
    # value; future actual capacity must not enter a historical forecast.
    result['utilization_ratio'] = result.utilization_ratio.clip(0, 1)
    result['cdu_capacity_kbd'] = float(history.cdu_capacity_kbd.iloc[-1])
    result['demand_kbd'] = result.utilization_ratio * result.cdu_capacity_kbd
    for c in ['production_kbd', 'imports_kbd', 'exports_kbd', 'direct_use_kbd']:
        result[c] = result[c].clip(lower=0)
    return result[FLOWS + ['cdu_capacity_kbd', 'utilization_ratio']]

# Build target-month inputs from history ending one month earlier. History must
# be sorted, continuous, and contain at least 12 monthly observations.
def feature_row(history, month):
    last = history.iloc[-1]
    # Use prior stock, year-ago stock, and the latest stock change. Sine/cosine encode
    # annual seasonality smoothly across the December/January boundary.
    row = {'stock_lag1': last.stock_kb, 'stock_lag12': history.iloc[-12].stock_kb,
           'change_lag1': last.stock_kb - history.iloc[-2].stock_kb,
           'sin_month': np.sin(2*np.pi*month.month/12), 'cos_month': np.cos(2*np.pi*month.month/12)}
    # Calculate observed changes before filtering; omit any change whose current
    # or previous stock observation falls inside the excluded COVID period.
    changes = history.set_index('month').stock_kb.diff()
    clean = outside_covid(changes.index) & outside_covid(changes.index - pd.DateOffset(months=1))
    changes = changes.loc[clean].tail(60)
    # Average recent eligible changes for the target calendar month; use zero
    # if that month has no usable historical changes.
    seasonal = changes[changes.index.month == month.month].dropna()
    row['seasonal_change'] = float(seasonal.mean()) if len(seasonal) else 0.0
    # Compare the latest stock with earlier same-calendar-month stock levels,
    # excluding the latest observation from its own reference average.
    normal_history = history.iloc[:-1].loc[lambda f: outside_covid(f.month)].tail(60)
    normal = normal_history[normal_history.month.dt.month.eq(last.month.month)].stock_kb
    row['stock_seasonal_gap'] = float(last.stock_kb - normal.mean()) if len(normal) else 0.0
    # Pre-sign prior-month flows so nonnegative regression coefficients retain the
    # supply-versus-demand direction in the constrained model.
    row.update({'lag_' + c: last[c] * sign for c, sign in zip(FLOWS, SIGNS)})
    return row


# Create labeled rows after a 12-month warm-up. Rebuild each feature vector and
# flow-identity benchmark using only history preceding that row's target month.
def supervised(history):
    rows = []
    for i in range(12, len(history)):
        # Keep the target separate: its stock supplies actual_kb/delta_kb only. Actual
        # target-month flows must not become inputs to a prediction of that month.
        past, current = history.iloc[:i], history.iloc[i]
        row = feature_row(past, current.month)
        flow = forecast_flows(past, [current.month]).iloc[0]
        row.update(month=current.month, origin_month=past.month.iloc[-1],
                   actual_kb=current.stock_kb, delta_kb=current.stock_kb-past.stock_kb.iloc[-1],
                   forecast_flow_identity=max(0.0, past.stock_kb.iloc[-1] + flow[FLOWS].to_numpy() @ SIGNS * current.month.days_in_month))
        rows.append(row)
    return pd.DataFrame(rows)

# Inspect a few supervised rows to make the forecast origin/target alignment
# concrete. Each origin must precede its target; actual_kb is an outcome only.

example=supervised(panel[panel.padd.eq(3)].reset_index(drop=True))
display(example[['origin_month','month','stock_lag1','lag_production_kbd','lag_demand_kbd',
                 'forecast_flow_identity','actual_kb']].tail())
assert (example.origin_month<example.month).all()

# Baselines carry stocks forward, repeat year-ago stocks, add a seasonal change,
# or accumulate forecast net flows; none needs a fitted stock estimator.
BASELINES = ['persistence', 'seasonal_naive', 'seasonal_change', 'forecast_flow_identity']
# Candidates share forecast timing and evaluation. Names ending in '_change'
# model stock increments, while constrained_level models inventory levels.
LEARNED = ['constrained_level', 'ridge_change', 'seasonal_ridge_change',
           'polynomial_ridge_change', 'spline_ridge_change', 'random_forest_change',
           'xgboost_change', 'neural_network_change', 'rolling_ridge_change',
           'damped_ridge_change', 'ensemble_change']
LEARNED += ['seasonal_residual_change', 'huber_change']
MODELS = BASELINES + LEARNED
BASE_FEATURES = ['stock_lag1', 'stock_lag12', 'change_lag1'] + ['lag_' + c for c in FLOWS]
FEATURES = BASE_FEATURES + ['sin_month', 'cos_month']

# Select predictor columns explicitly. actual_kb and delta_kb are labels and
# must never enter the estimator's input matrix.
def columns(name):
    if name == 'seasonal_residual_change':
        return FEATURES + ['stock_seasonal_gap', 'seasonal_change']
    if name == 'constrained_level':
        return ['stock_lag1'] + ['lag_' + c for c in FLOWS]
    if name == 'ridge_change':
        return BASE_FEATURES
    return FEATURES


# Return an unfitted candidate. Pipeline transformations are learned inside each
# fit, so test observations cannot influence scaling or spline preparation.
def estimator(name):
    if name == 'seasonal_residual_change':
        # Standardize predictors before L2-regularized regression. Larger alpha shrinks
        # coefficients more strongly; grid search overrides these starting defaults.
        return make_pipeline(StandardScaler(), Ridge(alpha=100))
    if name == 'huber_change':
        return TransformedTargetRegressor(regressor=make_pipeline(StandardScaler(),
            # Huber loss reduces outlier influence; the surrounding target transformer
            # scales the response and reverses that scaling when predicting.
            HuberRegressor(max_iter=1000)), transformer=StandardScaler())
    if name == 'constrained_level':
        # Constrain coefficients, not the intercept, to be nonnegative. Pre-signed demand
        # features therefore contribute negatively without imposing a full flow identity.
        return LinearRegression(positive=True)
    if name in ['ridge_change', 'seasonal_ridge_change', 'rolling_ridge_change', 'damped_ridge_change']:
        # Standardize predictors before L2-regularized regression. Larger alpha shrinks
        # coefficients more strongly; grid search overrides these starting defaults.
        return make_pipeline(StandardScaler(), Ridge(alpha=100))
    if name == 'polynomial_ridge_change':
        # Expand predictors with powers/interactions, rescale, and regularize to control
        # flexibility. The regression supplies the intercept, so no bias column is needed.
        return make_pipeline(StandardScaler(), PolynomialFeatures(2, include_bias=False),
                             StandardScaler(), Ridge(alpha=1000))
    if name == 'spline_ridge_change':
        # Represent each feature using smooth piecewise-polynomial basis functions;
        # linear extrapolation controls behavior outside its observed training range.
        return make_pipeline(SplineTransformer(n_knots=4, degree=2, extrapolation='linear'),
                             StandardScaler(), Ridge(alpha=100))
    if name == 'random_forest_change':
        # Average resampled decision trees. Depth and minimum leaf size limit complexity;
        # the fixed seed makes randomized fitting reproducible.
        return RandomForestRegressor(n_estimators=150, max_depth=4, min_samples_leaf=12,
                                     max_features=0.8, n_jobs=1, random_state=42)
    if name == 'xgboost_change':
        # Boost shallow trees sequentially. Learning rate, regularization, and row/column
        # subsampling control how aggressively the ensemble fits remaining residuals.
        return XGBRegressor(n_estimators=120, max_depth=2, learning_rate=0.03,
            min_child_weight=12, reg_lambda=30, subsample=0.8, colsample_bytree=0.8,
            n_jobs=1, random_state=42, objective='reg:squarederror')
    if name == 'neural_network_change':
        return TransformedTargetRegressor(regressor=make_pipeline(StandardScaler(),
            # Fit a small neural network with LBFGS on scaled inputs and a scaled target.
            # Alpha penalizes weights, and the random seed fixes initialization.
            MLPRegressor(hidden_layer_sizes=(16,), alpha=10, solver='lbfgs',
                         max_iter=3000, random_state=42)), transformer=StandardScaler())
    if name == 'ensemble_change':
        # Average seasonal ridge, random forest, and XGBoost predictions. The tuning
        # grid varies their relative weights.
        return VotingRegressor([(n, estimator(n)) for n in
            ['seasonal_ridge_change', 'random_forest_change', 'xgboost_change']])
    raise ValueError(name)

# Return (estimator, warning_text). Baselines return None because predict computes
# them directly; learned candidates filter training dates before fitting.
def fit(name, train, params=None):
    if name in BASELINES:
        return None, ''
    train = training_rows(train)
    model = estimator(name)
    if params:
        model.set_params(**params)
    if name == 'rolling_ridge_change':
        train = train.tail(60)
    # Choose inventory level or monthly change as response, both in thousand barrels.
    target = 'actual_kb' if name == 'constrained_level' else 'delta_kb'
    # Match the search's numerical settings, including neural-network refits.
    from threadpoolctl import threadpool_limits
    with warnings.catch_warnings(record=True) as caught, threadpool_limits(limits=1):
        warnings.simplefilter('always', ConvergenceWarning)
        # Fit the seasonal-residual candidate to departures from normal seasonal change;
        # predict reconstructs the seasonal component before converting to stock levels.
        values = train[target] - train.seasonal_change if name == 'seasonal_residual_change' else train[target]
        model.fit(train[columns(name)], values)
    return model, '; '.join(str(w.message) for w in caught)


# Return stock levels in kb for every candidate, including change-based models.
# Use identical reconstruction and flooring rules during tuning and forecasting.
def predict(name, model, frame):
    if name == 'persistence':
        return frame.stock_lag1.to_numpy()
    if name == 'seasonal_naive':
        return frame.stock_lag12.to_numpy()
    if name == 'seasonal_change':
        return np.maximum(0, frame.stock_lag1.to_numpy() + frame.seasonal_change.to_numpy())
    if name == 'forecast_flow_identity':
        return frame.forecast_flow_identity.to_numpy()
    prediction = model.predict(frame[columns(name)])
    if name == 'seasonal_residual_change':
        prediction += frame.seasonal_change.to_numpy()
    # Halve the estimated change to shrink the forecast toward unchanged stocks.
    if name == 'damped_ridge_change':
        prediction *= 0.5
    # Add the predicted increment to prior stock to recover an inventory level.
    if name != 'constrained_level':
        prediction += frame.stock_lag1.to_numpy()
    # Prevent a learned model from producing a negative final stock level.
    return np.maximum(prediction, 0)

# Enumerate tunable combinations. Double underscores address nested pipeline
# parameters; an empty grid indicates a candidate with no parameter search.
def parameter_grid(name):
    """Conventional, compact Cartesian grids; structural baselines stay fixed."""
    if name == 'seasonal_residual_change':
        return {'ridge__alpha': [0.1, 1, 10, 100, 1000]}
    if name == 'huber_change':
        return {'regressor__huberregressor__epsilon': [1.35, 1.75],
                'regressor__huberregressor__alpha': [0.01, 1, 10]}
    if name in ['ridge_change', 'seasonal_ridge_change', 'rolling_ridge_change', 'damped_ridge_change']:
        return {'ridge__alpha': [0.01, 0.1, 1, 10, 100, 1000]}
    if name == 'polynomial_ridge_change':
        return {'polynomialfeatures__degree': [1, 2], 'ridge__alpha': [1, 10, 100, 1000]}
    if name == 'spline_ridge_change':
        return {'splinetransformer__n_knots': [3, 5, 7],
                'splinetransformer__degree': [2, 3], 'ridge__alpha': [1, 100, 1000]}
    if name == 'random_forest_change':
        return {'n_estimators': [100, 200], 'max_depth': [3, 5, None],
                'min_samples_leaf': [5, 12]}
    if name == 'xgboost_change':
        return {'n_estimators': [100, 200], 'max_depth': [2, 4],
                'learning_rate': [0.01, 0.05, 0.1], 'reg_lambda': [1, 30]}
    if name == 'neural_network_change':
        return {'regressor__mlpregressor__hidden_layer_sizes': [(16,), (32,), (32, 16)],
                'regressor__mlpregressor__alpha': [0.01, 1, 10]}
    if name == 'ensemble_change':
        return {'weights': [(1, 1, 1), (2, 1, 1), (1, 2, 1), (1, 1, 2)]}
    return {}


# Adapt candidates to sklearn's cloning/search interface. Score reconstructed
# stock levels, including seasonal adjustments, damping, and the zero floor.
class StockRegressor(RegressorMixin, BaseEstimator):
    """Score the actual stock forecast, including damping, floor and rolling fit."""
    def __init__(self, name, model):
        self.name = name
        self.model = model

    def fit(self, X, y):
        from sklearn.base import clone
        # The sliced frame carries labels and dates for eligibility filtering; y is
        # accepted for sklearn compatibility, while the response is selected below.
        eligible = training_rows(X)
        train = eligible.tail(60) if self.name == 'rolling_ridge_change' else eligible
        target = train.actual_kb if self.name == 'constrained_level' else train.delta_kb
        if self.name == 'seasonal_residual_change':
            target = target - train.seasonal_change
        # Clone an unfitted candidate for each fold so fitted state cannot leak across
        # search trials. A trailing underscore conventionally marks fitted attributes.
        self.model_ = clone(self.model).fit(train[columns(self.name)], target)
        return self

    def predict(self, X):
        return predict(self.name, self.model_, X)


# Search only the supplied past window; return parameters, all candidate/fold
# scores, and convergence warnings. The chosen estimator is fitted separately.
def tune(name, train, folds=10, jobs=8, splits=None):
    """Search only the supplied past history; retain every candidate and fold score."""
    grid = parameter_grid(name)
    if not grid:
        return {}, pd.DataFrame(), ''
    cv = splits if splits is not None else validation_splits(train, folds)
    # Maximize negative MAE, equivalent to minimizing stock error. model__ reaches
    # the candidate inside StockRegressor; refit=False avoids an unused final fit.
    search = GridSearchCV(StockRegressor(name, estimator(name)),
        {'model__' + key: value for key, value in grid.items()},
        scoring='neg_mean_absolute_error', cv=cv, n_jobs=jobs,
        refit=False, error_score='raise', return_train_score=True)
    # Threads keep warnings in this process; native numerical threads are capped.
    from joblib import parallel_backend
    from threadpoolctl import threadpool_limits
    with warnings.catch_warnings(record=True) as caught, threadpool_limits(limits=1), parallel_backend('threading'):
        warnings.simplefilter('always', ConvergenceWarning)
        search.fit(train, train.actual_kb)
    results = pd.DataFrame(search.cv_results_)
    results['params'] = results.params.map(lambda p: json.dumps({k.removeprefix('model__'): v for k,v in p.items()}, sort_keys=True))
    # Reverse sklearn's negative-score convention to report positive MAE.
    results['mean_test_mae_kb'] = -results.mean_test_score
    # Remove the wrapper prefix so fit() can set parameters on the candidate itself.
    params = {k.removeprefix('model__'): v for k,v in search.best_params_.items()}
    warning = '; '.join(sorted({str(w.message) for w in caught}))
    return params, results, warning

# Show each candidate's search space before running the experiment; nested
# parameter names identify the relevant estimator or pipeline transformation.

display(pd.DataFrame([{'model':name,'param_grid':parameter_grid(name)} for name in LEARNED]))

# MAE and RMSE are in kb; RMSE weights large misses more. R-squared is unitless
# and may be negative when errors exceed the actual sample's mean-only benchmark.
def score(actual, predicted):
    return {'mae_kb': mean_absolute_error(actual, predicted),
            'rmse_kb': mean_squared_error(actual, predicted)**0.5,
            'r2': r2_score(actual, predicted), 'n': len(actual)}


# Require one row for every PADD in each group before summing additive columns.
# Compute national errors after summation; regional forecast errors can offset.
def aggregate(frame, keys, additive):
    """Reject partial U.S. totals instead of silently summing fewer than five PADDs."""
    if frame.duplicated(keys + ['padd']).any():
        raise ValueError('Duplicate PADD rows in aggregation')
    if not frame.groupby(keys).padd.nunique().eq(5).all():
        raise ValueError('U.S. aggregation requires all five PADDs')
    return frame.groupby(keys, as_index=False)[additive].sum(min_count=5)


# Score each group, expanding either a scalar or tuple group key into the named
# key columns of the resulting metrics table.
def metric_table(predictions, keys):
    return pd.DataFrame([{**dict(zip(keys, k if isinstance(k, tuple) else (k,))),
        **score(g.actual_kb, g.predicted_kb)} for k,g in predictions.groupby(keys)])

# Inspect the candidate roster and example estimator configurations. The count
# assertion guards the expected 17-model comparison; it does not evaluate accuracy.

assert len(MODELS)==17
print('Complete constrained features:',columns('constrained_level'))
print(estimator('neural_network_change'))
print(estimator('ensemble_change'))

Complete constrained features: ['stock_lag1', 'lag_production_kbd', 'lag_demand_kbd', 'lag_imports_kbd', 'lag_exports_kbd', 'lag_net_receipts_kbd', 'lag_adjustments_kbd', 'lag_transfers_kbd', 'lag_direct_use_kbd']
TransformedTargetRegressor(regressor=Pipeline(steps=[('standardscaler',
                                                      StandardScaler()),
                                                     ('mlpregressor',
                                                      MLPRegressor(alpha=10,
                                                                   hidden_layer_sizes=(16,),
                                                                   max_iter=3000,
                                                                   random_state=42,
                                                                   solver='lbfgs'))]),
                           transformer=StandardScaler())
VotingRegressor(estimators=[('seasonal_ridge_change',
                             Pipeline(steps=[('standardscaler',
                                              StandardScaler()),
                                             ('ridge', Ridge(alpha=100))])),
                            ('random_forest_change',
                             RandomForestRegressor(max_depth=4,
                                                   max_features=0.8,
                                                   min_samples_leaf=12,
                                                   n_estimators=150, n_jobs=1,
                                                   random_state=42)),
                            ('xgboost_change',
                             XGBRegressor(base_score=None, booster=None,
                                          callbacks=None,
                                          c...
                                          feature_weights=None, gamma=None,
                                          grow_policy=None,
                                          importance_type=None,
                                          interaction_constraints=None,
                                          learning_rate=0.03, max_bin=None,
                                          max_cat_threshold=None,
                                          max_cat_to_onehot=None,
                                          max_delta_step=None, max_depth=2,
                                          max_leaves=None, min_child_weight=12,
                                          missing=nan,
                                          monotone_constraints=None,
                                          multi_strategy=None, n_estimators=120,
                                          n_jobs=1, num_parallel_tree=None, ...))])

# Choose each PADD's one-step candidate on development data, evaluate the fixed
# choices on the last 24 months, then refit winners for the live outlook.
def evaluate(samples, folds=10, holdout=24, jobs=8):
    """Choose models only on development CV; score untouched final 24 months."""
    predictions, fold_metrics, selections, fitted, warning_rows = [], [], [], {}, []
    best_parameters, searches = {}, []
    for padd, frame in samples.items():
        # Reserve the final holdout rows before any tuning or candidate selection.
        development = frame.iloc[:-holdout]
        if len(development) - folds*12 < 36:
            raise ValueError('Need >=36 initial training months plus 10 annual folds and holdout')
        splits = validation_splits(development, folds, test_size=12)
        for name in MODELS:
            params, results, warning = tune(name, development, folds, jobs, splits)
            best_parameters[padd, name] = params
            if not results.empty:
                searches.append(results.assign(padd=padd, model=name, stage='development',
                    train_end=development.month.max()))
                print(f'PADD {padd} {name}: grid MAE {results.mean_test_mae_kb.min():,.1f}; {params}', flush=True)
            if warning:
                warning_rows.append(dict(padd=padd, model=name, stage='grid_search', warning=warning))
            # Reuse tuning folds for development diagnostics; these are selection scores,
            # not independent nested-CV estimates of future generalization.
            for fold, (tr, te) in enumerate(splits, 1):
                train, test = training_rows(development.iloc[tr]), development.iloc[te]
                model, warning = fit(name, train, params)
                pred = predict(name, model, test)
                fold_metrics.append(dict(padd=padd, model=name, fold=fold,
                    train_end=train.month.max(), test_start=test.month.min(), test_end=test.month.max(),
                    train_mae_kb=mean_absolute_error(train.actual_kb, predict(name, model, train)),
                    **score(test.actual_kb, pred)))
                output = test[['month', 'origin_month', 'actual_kb']].copy()
                output['predicted_kb'], output['padd'], output['model'] = pred, padd, name
                output['split'], output['fold'] = 'cv', fold
                predictions.append(output)
                if warning:
                    warning_rows.append(dict(padd=padd, model=name, stage=f'cv_{fold}', warning=warning))
        current = pd.concat(predictions, ignore_index=True)
        current = current[current.padd.eq(padd)]
        # Pool absolute development errors for this PADD across folds and select
        # the lowest mean; no holdout observations enter this comparison.
        losses = current.assign(error=lambda x: abs(x.actual_kb-x.predicted_kb)).groupby('model').error.mean()
        winner = losses.idxmin()
        selections.append(dict(padd=padd, selected_model=winner, cv_mae_kb=losses[winner]))
        for name in MODELS:
            model, warning = fit(name, development, best_parameters[padd, name])
            # This is one-step evaluation: each row uses actual prior-month inputs, while
            # estimator parameters remain fitted on development data for the entire block.
            test = frame.iloc[-holdout:]
            output = test[['month', 'origin_month', 'actual_kb']].copy()
            output['predicted_kb'] = predict(name, model, test)
            output['padd'], output['model'], output['split'], output['fold'] = padd, name, 'holdout', 0
            predictions.append(output)
            if warning:
                warning_rows.append(dict(padd=padd, model=name, stage='holdout', warning=warning))
        # Refit the already chosen winner on all eligible history for future use.
        # This final estimator is separate from the one that produced holdout scores.
        final, warning = fit(winner, frame, best_parameters[padd, winner])
        fitted[padd] = {'name': winner, 'estimator': final}
        if warning:
            warning_rows.append(dict(padd=padd, model=winner, stage='final', warning=warning))
        print(f'PADD {padd}: selected {winner}; CV MAE {losses[winner]:,.1f} kb', flush=True)
    return (pd.concat(predictions, ignore_index=True), pd.DataFrame(fold_metrics),
            pd.DataFrame(selections), fitted,
            pd.DataFrame(warning_rows, columns=['padd','model','stage','warning']),
            best_parameters, pd.concat(searches, ignore_index=True))

# Produce a recursive path per district plus national sums. After the origin,
# features depend on predicted stocks/flows, not future actual observations.
def forecast(history, fitted, horizon=12):
    all_paths = []
    for p, state in fitted.items():
        past = history[history.padd.eq(p)].sort_values('month').copy()
        origin = past.month.iloc[-1]
        months = pd.date_range(origin + pd.offsets.MonthBegin(), periods=horizon, freq='MS')
        # Forecast the full flow path once at the origin; do not refit flows each step
        # using synthetic observations appended to the stock path.
        flows = forecast_flows(past, months)
        # Track a separate raw flow-identity path initialized at the last observed stock.
        identity_stock = past.stock_kb.iloc[-1]
        for h, month in enumerate(months, 1):
            previous = past.stock_kb.iloc[-1]
            flow = flows.loc[month]
            balance = flow[FLOWS].to_numpy() @ SIGNS
            # Crude balances are daily rates; multiply by target-month days to obtain
            # the monthly stock increment in kb.
            volume = balance * month.days_in_month
            # At longer horizons, past contains earlier predicted rows. This is the recursive
            # feedback through which forecast errors can propagate along the path.
            features = feature_row(past, month)
            features['forecast_flow_identity'] = max(0.0, previous + volume)
            stock = float(predict(state['name'], state['estimator'], pd.DataFrame([features]))[0])
            # Keep the raw flow accumulation unclipped so impossible inventory paths remain
            # visible; it can diverge from the statistical stock forecast.
            identity_stock += volume
            row = dict(month=month, origin_month=origin, horizon=h, padd=p,
                padd_name=PADD_NAMES[p], model=state['name'], stock_kb=stock,
                previous_stock_kb=previous, stock_change_kb=stock-previous,
                identity_stock_unclipped_kb=identity_stock, balance_kbd=balance,
                # Express statistical stock change minus forecast net supply as a daily
                # rate. This discrepancy is not an extra observed EIA adjustment.
                model_reconciliation_kbd=(stock-previous)/month.days_in_month-balance,
                supply_kbd=flow.production_kbd+flow.imports_kbd+flow.net_receipts_kbd+flow.adjustments_kbd+flow.transfers_kbd,
                total_demand_kbd=flow.demand_kbd+flow.exports_kbd+flow.direct_use_kbd, **flow.to_dict())
            all_paths.append(row)
            # Append the predicted row so the next month uses this stock and these flows
            # as its most recent history, without seeing future actual values.
            past = pd.concat([past, pd.DataFrame([row])], ignore_index=True)
    regional = pd.DataFrame(all_paths)
    sums = FLOWS + ['stock_kb', 'previous_stock_kb','stock_change_kb', 'identity_stock_unclipped_kb',
                   'balance_kbd','model_reconciliation_kbd','supply_kbd','total_demand_kbd','cdu_capacity_kbd']
    national = aggregate(regional, ['month','origin_month','horizon'], sums)
    # National utilization is total refinery input divided by total capacity,
    # not a sum or unweighted average of regional utilization ratios.
    national['utilization_ratio'] = national.demand_kbd / national.cdu_capacity_kbd
    return regional, national

# Run controls: use ten expanding validation folds and eight search workers.
# The worker count controls parallel grid trials; numerical threads are capped
# inside fitting to avoid multiplying the number of active CPU threads.

data_dir=HERE/'data'
output_dir=OUTPUT
folds=10
jobs=8

output = Path(output_dir)
# Build chronological supervised samples separately for all five districts,
# then run the full development search, one-step holdout evaluation, and refits.
# This is a computationally expensive cell; progress and warnings are retained.
output.mkdir(parents=True, exist_ok=True)
panel = load_panel(data_dir)
histories = {p: panel[panel.padd.eq(p)].reset_index(drop=True) for p in PADD_NAMES}
samples = {p: supervised(history) for p,history in histories.items()}
predictions, fold_metrics, selection, fitted, fit_warnings, best_parameters, grid_results = evaluate(samples, folds, jobs=jobs)

display(fold_metrics.groupby('model').agg(folds=('fold','count'),test_mae_kb=('mae_kb','mean')))
display(selection)

# For crude, compare the mixture of regional winners with each common model
# family on national development MAE. Keep regional winners for audit, but deploy
# the model with the smallest national error without using holdout scores.

# Label the development-selected regional forecasts as a combined model while
# retaining every original candidate's predictions for comparison.
selected = predictions.merge(selection[['padd','selected_model']], on='padd')
selected = selected[selected.model.eq(selected.selected_model)].copy()
selected['model'] = 'selected_padd_models'
combined = pd.concat([predictions, selected[predictions.columns]], ignore_index=True)
us_predictions = aggregate(combined, ['month','origin_month','model','split','fold'], ['actual_kb','predicted_kb'])
# Choose the aggregation model on development predictions only. A single
# common model family can outperform the mix of regional winners nationally.
one_month_model_scores = metric_table(us_predictions[us_predictions.split.eq('cv')], ['model'])
# Select a common-family model or the mix of PADD winners by national CV MAE.
# This extra selection is also based on development data and can be optimistic.
one_month_model = one_month_model_scores.sort_values(['mae_kb','model']).iloc[0].model
# Preserve regional-only winners before a possible national-model override.
regional_selection = selection.copy()
# When a common family wins nationally, refit that family separately for each
# PADD and update deployment labels while retaining the original comparisons.
if one_month_model != 'selected_padd_models':
    for p, frame in samples.items():
        fitted_estimator, warning = fit(one_month_model, frame, best_parameters[p, one_month_model])
        fitted[p] = {'name': one_month_model, 'estimator': fitted_estimator}
        if warning:
            fit_warnings.loc[len(fit_warnings)] = [p, one_month_model, 'national_model_final', warning]
    selection['selected_model'] = one_month_model
    selection['cv_mae_kb'] = [mean_absolute_error(g.actual_kb, g.predicted_kb)
        for p in selection.padd
        for g in [predictions.query("padd == @p and model == @one_month_model and split == 'cv'")]]
padd_metrics = metric_table(combined, ['padd','model','split'])
us_metrics = metric_table(us_predictions, ['model','split'])
# Generate the next month using the selected one-step estimators. The separate
# year-ahead selection below may choose a different family for the 12-month path.
one_step_regional, one_step_national = forecast(panel, fitted, horizon=1)

display(one_month_model_scores.sort_values('mae_kb'))
print('One-month model:', one_month_model)

# Evaluate every candidate on six complete development-year forecast paths.
# Retune at each historical origin, sum PADD predictions, then choose the common
# year-ahead family using national MAE across all path months.

# Long-horizon model selection is a different task. Six disjoint development
# years end before the final holdout; all candidates use the actual recursive
# forecast function, without seeing any test-year observations.
development_end = samples[1].month.iloc[-25]
recursive_cv = []
for cutoff in recursive_origins(development_end):
    for name in MODELS:
        states = {}
        for p in PADD_NAMES:
            # Restrict fitting to targets observed by this simulated forecast origin.
            train = samples[p][samples[p].month.le(cutoff)]
            # Retune on this origin's history only; later development observations must not
            # inform a model used to forecast an earlier evaluation year.
            params, results, tuning_warning = tune(name, train, folds, jobs)
            if not results.empty:
                grid_results = pd.concat([grid_results, results.assign(padd=p, model=name,
                    stage='recursive_development', train_end=cutoff)], ignore_index=True)
            if tuning_warning:
                fit_warnings.loc[len(fit_warnings)] = [p, name, 'recursive_grid_search', tuning_warning]
            model, warning = fit(name, train, params)
            states[p] = {'name': name, 'estimator': model}
            if warning:
                fit_warnings.loc[len(fit_warnings)] = [p, name, 'recursive_cv', warning]
        path, _ = forecast(panel[panel.month.le(cutoff)], states)
        path = path[['padd','month','origin_month','horizon','stock_kb']].rename(columns={'stock_kb':'predicted_kb'})
        # Attach actual stocks only after forecasting, for scoring. one_to_one rejects
        # duplicate month/PADD observations that could otherwise multiply rows.
        path = path.merge(panel[['padd','month','stock_kb']], on=['padd','month'], validate='one_to_one').rename(columns={'stock_kb':'actual_kb'})
        path['model'] = name
        recursive_cv.append(path)
    print(f'Recursive development origin {cutoff.date()} complete', flush=True)
recursive_cv = pd.concat(recursive_cv, ignore_index=True)
us_recursive_cv = aggregate(recursive_cv, ['month','origin_month','horizon','model'], ['actual_kb','predicted_kb'])
year_ahead_model_scores = metric_table(us_recursive_cv, ['model']).sort_values(['mae_kb','model'])
# Choose the lowest national development-path MAE, with model name as a stable
# tie-breaker; freeze the family before examining recursive holdout results.
year_ahead_model = year_ahead_model_scores.iloc[0].model
# Refit the selected year-ahead family on full eligible history, using parameter
# settings chosen before the holdout; then build the forward 12-month outlook.
horizon_fitted = {}
for p in PADD_NAMES:
    model, warning = fit(year_ahead_model, samples[p], best_parameters[p, year_ahead_model])
    horizon_fitted[p] = {'name': year_ahead_model, 'estimator': model}
    if warning:
        fit_warnings.loc[len(fit_warnings)] = [p, year_ahead_model, 'horizon_final', warning]
regional, national = forecast(panel, horizon_fitted)
print(f'National one-step model: {one_month_model}; 12-month model: {year_ahead_model}', flush=True)

display(year_ahead_model_scores)
print('Twelve-month model:', year_ahead_model)

# Evaluate the fixed year-ahead model, recursive one-step model, and benchmarks
# at 13 monthly origins. Refit on history available at each origin using frozen
# development parameters; overlapping 12-month paths have correlated errors.

# Thirteen rolling 12-month paths lying wholly inside the final holdout.
# They overlap and are correlated; do not interpret them as 13 independent years.
recursive = []
# Move the origin forward one month at a time; each next 12 months remains
# inside the final 24-month evaluation period.
evaluation_models = list(dict.fromkeys([year_ahead_model, one_month_model, 'persistence', 'constrained_level']))
for offset in range(24, 11, -1):
    cutoff = panel.month.max() - pd.DateOffset(months=offset)
    for name in evaluation_models:
        states = {}
        for p, state in fitted.items():
            # Restrict fitting to targets observed by this simulated forecast origin.
            train = samples[p][samples[p].month.le(cutoff)]
            model, warning = fit(name, train, best_parameters[p, name])
            states[p] = {'name': name, 'estimator': model}
            if warning:
                fit_warnings.loc[len(fit_warnings)] = [p, name, 'recursive_holdout', warning]
        path, _ = forecast(panel[panel.month.le(cutoff)], states)
        path['model'] = name
        # Attach actual stocks only after forecasting, for scoring. one_to_one rejects
        # duplicate month/PADD observations that could otherwise multiply rows.
        path = path.merge(panel[['padd','month','stock_kb']], on=['padd','month'], suffixes=('_forecast','_actual'), validate='one_to_one')
        path = path.rename(columns={'stock_kb_forecast':'predicted_kb','stock_kb_actual':'actual_kb'})
        recursive.append(path)
recursive = pd.concat(recursive, ignore_index=True)
us_recursive = aggregate(recursive, ['month','origin_month','horizon','model'], ['actual_kb','predicted_kb'])

display(metric_table(us_recursive, ['model']))

# Combine regional history into national totals and show the next-month
# forecast. The latest table contains five PADD rows plus one U.S. reporting row.

# Build national historical totals from all five districts. Ratios, when needed,
# must be recomputed from their summed numerator and denominator.
us_history = aggregate(panel, ['month'], FLOWS+['stock_kb','cdu_capacity_kbd',
    'balance_kbd','accounting_residual_kbd','five_term_balance_kbd','five_term_residual_kbd'])
us_history['utilization_ratio'] = us_history.demand_kbd / us_history.cdu_capacity_kbd
# Use padd=0 as a reporting label for the U.S. total; it is not a sixth region.
latest_us = one_step_national.assign(padd=0, padd_name='United States (sum of PADDs)',
                                                   model=one_month_model)
latest = pd.concat([one_step_regional, latest_us], ignore_index=True)

display(us_history.tail(3))
display(latest[['padd','padd_name','month','stock_kb']])

# Save result tables and fitted models for reporting and reuse. CSVs carry
# readable diagnostics; joblib preserves estimator state; coefficients expose
# the constrained regression's fitted relationship.

# Collect inputs, predictions, scores, selections, and audits under stable names.
# Each name becomes a CSV stem consumed by later reporting cells.
tables = {'padd_monthly_model':panel,'us_monthly_model':us_history,
    'padd_predictions':combined,'us_predictions':us_predictions,'padd_model_metrics':padd_metrics,
    'us_model_metrics':us_metrics,'fold_metrics':fold_metrics,'selected_models':selection,
    'regional_cv_winners':regional_selection, 'national_model_cv_metrics':one_month_model_scores,
    'recursive_cv_predictions':recursive_cv, 'us_recursive_cv_predictions':us_recursive_cv,
    'year_ahead_model_cv_metrics':year_ahead_model_scores,
    'padd_forecast_12m':regional,'us_forecast_12m':national,'latest_forecast':latest,
    'recursive_holdout_predictions':recursive,'us_recursive_holdout_predictions':us_recursive,
    'recursive_holdout_metrics':metric_table(recursive,['padd','model']),
    'us_recursive_horizon_metrics':metric_table(us_recursive,['horizon','model']), 'fit_warnings':fit_warnings,
    'grid_search_results':grid_results,
    'best_parameters':pd.DataFrame([{'padd':p, 'model':name, 'params':json.dumps(params, sort_keys=True)}
        for (p,name),params in best_parameters.items()])}
# Record target and origin dates plus separate fitting/evaluation eligibility
# flags so readers can inspect exactly which dates the exclusion model affects.
tables['training_sample_audit'] = pd.concat([frame[['month','origin_month']].assign(
    padd=p, training_eligible=frame.index.isin(training_rows(frame).index),
    evaluation_eligible=outside_covid(frame.month)) for p, frame in samples.items()], ignore_index=True)
# Write named result tables without pandas' row index; rerunning replaces files
# with the same names in model_output.
for name, frame in tables.items():
    frame.to_csv(output/f'{name}.csv', index=False)
# Save all final candidate estimators as well as selected models so downstream
# analysis can reproduce predictions without repeating parameter searches.
all_fitted = {}
for p, frame in samples.items():
    for name in MODELS:
        model, warning = fit(name, frame, best_parameters[p, name])
        all_fitted[p, name] = {'name': name, 'estimator': model}
joblib.dump(all_fitted, output/'all_fitted_models.joblib')
joblib.dump(fitted, output/'fitted_models.joblib')
joblib.dump(horizon_fitted, output/'fitted_horizon_models.joblib')
# Export an interpretable constrained-regression coefficient table by PADD.
# These coefficients apply to pre-signed inputs, not unsigned raw demand flows.
coefficients = []
for p, frame in samples.items():
    model, _ = fit('constrained_level', frame)
    coefficients.append({'padd':p,'intercept_kb':model.intercept_,**dict(zip(columns('constrained_level'),model.coef_))})
pd.DataFrame(coefficients).to_csv(output/'constrained_model_coefficients.csv',index=False)
pd.DataFrame(coefficients).to_csv(output/'padd_model_coefficients.csv',index=False)

print('Saved',len(tables),'result tables to',output)

# Persist machine-readable run documentation alongside result tables. Source
# URLs identify provenance; versions and settings describe the actual run,
# including limitations of revised-history and point-only forecasts.

# Record scope, units, actual data dates, selection rules, limitations, and
# software versions alongside outputs for reproducibility and interpretation.
metadata = {'product':'Commercial crude oil; excludes strategic reserves',
    'covid_exclusion': {'start':'2020-03-01', 'end':'2021-03-01', 'inclusive':True,
        'method':'Exclude targets and direct lag inputs touching COVID from stock fitting; exclude COVID observations from flow fits and seasonal averages; preserve the full calendar and historical balance.',
        'evaluation':'Exclude COVID target months from development selection; full recursive development years avoid COVID; final 24 months unchanged.'},
    'stock_units':'thousand barrels at month end', 'flow_units':'thousand barrels per day',
    'equation':'S[t] = S[t-1] + days[t] * (field production[t] + imports[t] + net receipts[t] + adjustments[t] + transfers to crude[t] - refinery input[t] - exports[t] - direct use[t] + residual[t])',
    'data_start':str(panel.month.min().date()),'data_end':str(panel.month.max().date()),
    'cv':f'{folds} expanding folds, 12 eligible test months each; final 24 months excluded from selection',
    'holdout_start':str(samples[1].month.iloc[-24].date()),
    'forecast_timing':'One month after latest observed EIA month; conditional on prior monthly data being available. Current revised history, not a real-time release-vintage backtest.',
    'selection':'GridSearchCV per PADD/model on development data; PADD winners plus common-family models; lowest national development MAE chooses deployment; no holdout tuning',
    'hyperparameter_search':{'method':'GridSearchCV', 'folds':folds, 'scoring':'neg_mean_absolute_error',
        'splitter':'TimeSeriesSplit on non-COVID target dates; 12 eligible observations per development block; default block size inside recursive origins',
        'evaluation':'Development scores reuse tuning folds and are selection scores, not nested CV estimates. Final 24 months excluded from all searches.',
        'recursive':'Each development origin retunes on its past history. Holdout and final refits freeze parameters selected before holdout.'},
    'one_month_model':one_month_model,
    'year_ahead_model':year_ahead_model,
    'horizon_selection':'Lowest national MAE across six disjoint 12-month development paths, six full non-COVID development years; all 17 candidates compared',
    'flow_forecast':'Last 60 non-COVID observations, daily rates, linear trend and month fixed effects; crude run = forecast utilization ratio (clipped 0–1) times CDU capacity fixed at origin',
    'forecast_stock_floor':0, 'uncertainty':'Point forecasts only; no calibrated prediction intervals',
    'recursive_validation':'13 overlapping 12-month paths entirely in final holdout; refit at each origin; correlated errors',
    'versions':{'python':platform.python_version(),'numpy':np.__version__,'pandas':pd.__version__,
                'sklearn':sklearn.__version__, 'xgboost':xgboost.__version__},
    'sources':['https://www.eia.gov/dnav/pet/pet_sum_snd_d_r10_mbbl_m_cur.htm',
               'https://www.eia.gov/dnav/pet/TblDefs/pet_sum_snd_tbldef2.asp']}
(output/'model_metadata.json').write_text(json.dumps(metadata,indent=2))
# Serialize candidate definitions for inspection; repr describes configuration,
# while the joblib files above preserve actual fitted estimator state.
settings = {name: {'features': columns(name), 'estimator': repr(estimator(name)), 'param_grid':parameter_grid(name)}
            for name in LEARNED}
(output/'candidate_models.json').write_text(json.dumps(settings, indent=2))
print(us_metrics[us_metrics.split.eq('holdout')].sort_values('mae_kb').to_string(index=False))

print('Data through:',metadata['data_end'])
print('Fitting warnings:',len(fit_warnings))

folds=tables['fold_metrics']
# Inspect one representative model/PADD because candidates share split dates.
# Here folds is reassigned from the earlier integer setting to the metrics table.
# Assertions require training before testing and development before the holdout.
schedule=folds.query("padd==1 and model=='persistence'")[['fold','train_end','test_start','test_end','n']]
display(schedule)
assert (pd.to_datetime(schedule.train_end)<pd.to_datetime(schedule.test_start)).all()
assert pd.to_datetime(schedule.test_end).max()<pd.Timestamp(metadata['holdout_start'])

# Compare development MAE by candidate and district, then inspect the average
# train/test gap. A large gap suggests overfitting or a changing data regime;
# these reused selection scores are not an independent performance estimate.

display(tables['regional_cv_winners'])
display(tables['national_model_cv_metrics'].sort_values('mae_kb'))
cv=tables['padd_model_metrics'].query("split == 'cv' and model != 'selected_padd_models'")
display(cv.pivot(index='model',columns='padd',values='mae_kb').style.highlight_min(axis=0))
gaps=folds.groupby('model')[['train_mae_kb','mae_kb']].mean().sort_values('mae_kb')
# Subtract in-sample MAE from validation MAE; positive gaps mean worse
# performance on later observations than on fitted training observations.
gaps['test_minus_train_kb']=gaps.mae_kb-gaps.train_mae_kb
display(gaps)
gaps[['train_mae_kb','mae_kb']].plot.barh(figsize=(10,6),title='Mean training and annual-test errors across PADDs')
plt.xlabel('Thousand barrels');plt.show()

# Group validation errors by test-block endpoint to show regime variation.
# These plotted values average regional MAEs, not errors on national stock sums.
# The coefficient CSV provides a separate view of the constrained model.

fig,ax=plt.subplots(figsize=(12,4))
for name in ['persistence','seasonal_change','constrained_level','xgboost_change','rolling_ridge_change']:
    g=folds[folds.model.eq(name)].groupby('test_end').mae_kb.mean()
    ax.plot(pd.to_datetime(g.index),g.values,label=name,marker='o')
ax.set(title='Annual test performance changes across regimes',ylabel='Mean PADD MAE (kb)')
ax.legend(fontsize=8);plt.show()
display(pd.read_csv(OUTPUT/'constrained_model_coefficients.csv'))

# Report final one-step national MAE for the development-selected model against
# persistence. The percentage is 100 * (1 - selected MAE / baseline MAE), so a
# negative value means deterioration. Show settings without selecting on holdout.

scores=tables['us_model_metrics'].query("split=='holdout'").sort_values('mae_kb')
display(scores)
lookup=scores.set_index('model')
selected=lookup.loc[metadata['one_month_model'],'mae_kb']
naive=lookup.loc['persistence','mae_kb']

with pd.option_context('display.max_colwidth',None,'display.max_rows',None):
    display(tables['best_parameters'].query("params != '{}'").reset_index(drop=True))
comparison_names = [metadata['one_month_model'], 'constrained_level', 'persistence']
display(tables['padd_model_metrics'].query("split == 'holdout' and model in @comparison_names"))

# Summarize eligible dates and check that fitting excludes COVID target months.
# Optional comparison CSVs show earlier refits on matching dates when available.

audit = tables['training_sample_audit']
display(audit.groupby('padd').agg(total_target_months=('training_eligible','size'),
    eligible_stock_training_months=('training_eligible','sum'), evaluation_months=('evaluation_eligible','sum')))
assert not audit.loc[audit.training_eligible, 'month'].between(COVID_START, COVID_END).any()
# Display prior-run comparison files only when available; these optional reports
# are read here and are not recomputed by this cell.
comparison_file = OUTPUT/'refit_comparison.csv'
if comparison_file.exists():
    display(pd.read_csv(comparison_file))
addition_file = OUTPUT/'candidate_addition_comparison.csv'
if addition_file.exists():
    display(pd.read_csv(addition_file))

year_comparison = OUTPUT/'year_ahead_refit_comparison.csv'
if year_comparison.exists():
    display(pd.read_csv(year_comparison))

# Compare regional selected-model errors with unchanged stocks on the same
# holdout dates. Normalize MAE by average actual stocks to compare PADD sizes;
# positive percentage improvement means lower error than the benchmark.

# Compare like-for-like regional errors on the same final 24 months.
regional_test = tables['padd_predictions'].query("split == 'holdout'")
selected_names = tables['selected_models'][['padd','selected_model']]
regional_selected = regional_test.merge(selected_names, on='padd')
regional_selected = regional_selected[regional_selected.model.eq(regional_selected.selected_model)]
rows = []
for padd, g in regional_selected.groupby('padd'):
    baseline = regional_test[regional_test.padd.eq(padd) & regional_test.model.eq('persistence')]
    # Align selected and persistence forecasts on the same month within this PADD.
    # Validate unique matches so benchmark improvement is a like-for-like comparison.
    paired = g.merge(baseline[['month','predicted_kb']], on='month', suffixes=('', '_unchanged'), validate='one_to_one')
    # Positive signed error means overestimated stocks; absolute error measures
    # miss size, while mean signed error below measures systematic bias.
    error = paired.predicted_kb - paired.actual_kb
    baseline_mae = (paired.predicted_kb_unchanged - paired.actual_kb).abs().mean()
    # Keep the row label of the largest absolute miss to retrieve its calendar month.
    worst = error.abs().idxmax()
    rows.append({'PADD': padd, 'Average error (kb)': error.abs().mean(),
        'Error / average stocks (%)': 100 * error.abs().mean() / paired.actual_kb.mean(),
        'Improvement over unchanged stocks (%)': 100 * (1 - error.abs().mean() / baseline_mae),
        'Average bias (kb)': error.mean(), 'Worst month': paired.loc[worst, 'month'],
        'Largest miss (kb)': error.abs().max()})
regional_review = pd.DataFrame(rows)
display(regional_review)
# Persist regional diagnostics; charts convert kb to million barrels by /1000.
regional_review.to_csv(OUTPUT/'regional_error_review.csv', index=False)
fig, axes = plt.subplots(5, 1, figsize=(12, 10), sharex=True)
for ax, (padd, g) in zip(axes, regional_selected.groupby('padd')):
    ax.bar(g.month, (g.predicted_kb-g.actual_kb)/1000, width=20)
    ax.axhline(0, color='black', linewidth=.6)
    ax.set_ylabel(f'PADD {padd}\nmillion bbl')
    ax.set_title(g.selected_model.iloc[0], fontsize=9, loc='right')
axes[0].set_title('Monthly inventory forecast errors: positive means stocks were overestimated')
plt.tight_layout(); plt.show()

# Use one model's actual column to plot national stocks once, then overlay
# candidate forecasts and signed selected-model errors. Divide kb by 1000
# for million-barrel axes; a positive error bar means excess predicted stocks.

hold=tables['us_predictions'].query("split=='holdout'")
one_month_model = metadata['one_month_model']
actual=hold[hold.model.eq(one_month_model)]
fig,axes=plt.subplots(2,1,figsize=(12,7),sharex=True)
axes[0].plot(actual.month,actual.actual_kb/1000,color='black',linewidth=2,label='Actual')
for name in [one_month_model,'constrained_level','persistence']:
    g=hold[hold.model.eq(name)]
    axes[0].plot(g.month,g.predicted_kb/1000,label=name,alpha=.8)
axes[0].set(title='Final one-month evaluation: U.S. commercial crude stocks',ylabel='Million barrels');axes[0].legend()
axes[1].bar(actual.month,(actual.predicted_kb-actual.actual_kb)/1000,width=20)
axes[1].axhline(0,color='black');axes[1].set(ylabel='Forecast − actual (million bbl)',title=f'{one_month_model} forecast errors')
plt.tight_layout();plt.show()

# Check recursive development targets precede the holdout, then summarize final
# recursive errors across paths and by lead month. Later horizons repeatedly
# use predicted history, unlike the actual-prior-input one-step evaluation.

display(tables['year_ahead_model_cv_metrics'])
cv_paths=tables['us_recursive_cv_predictions']
display(cv_paths[cv_paths.model.eq('persistence')].groupby('origin_month').agg(first_target=('month','min'),last_target=('month','max')))
assert cv_paths.month.max()<pd.Timestamp(metadata['holdout_start'])
recursive=tables['us_recursive_holdout_predictions'].copy()
recursive['abs_error_kb']=(recursive.predicted_kb-recursive.actual_kb).abs()
recursive_summary=recursive.groupby('model').abs_error_kb.agg(['mean','median','max'])
display(recursive_summary)
horizon=tables['us_recursive_horizon_metrics']
display(horizon.pivot(index='horizon',columns='model',values='mae_kb'))
horizon.pivot(index='horizon',columns='model',values='mae_kb').plot(figsize=(11,4),marker='o',title='Final recursive error by forecast horizon')
plt.ylabel('National MAE (kb)');plt.show()

# Compare recursive holdout MAE with persistence and plot the earliest/latest
# forecast origins as examples. These overlapping paths are correlated evidence,
# so the observed comparison is not a guarantee of future year-ahead skill.

recursive_mae=recursive_summary.loc[metadata['year_ahead_model'],'mean']
persistence_mae=recursive_summary.loc['persistence','mean']

fig,axes=plt.subplots(1,2,figsize=(14,4))
origins=sorted(recursive.origin_month.unique())
for origin,ax in zip([origins[0],origins[-1]],axes):
    g=recursive[recursive.origin_month.eq(origin)]
    a=g[g.model.eq(metadata['year_ahead_model'])]
    ax.plot(a.month,a.actual_kb/1000,color='black',label='Actual')
    comparison_models = list(dict.fromkeys([metadata['year_ahead_model'], metadata['one_month_model'], 'persistence']))
    for name in comparison_models:
        t=g[g.model.eq(name)];ax.plot(t.month,t.predicted_kb/1000,label=name)
    ax.set(title=f'Origin {pd.Timestamp(origin):%Y-%m}',ylabel='Million barrels');ax.legend(fontsize=8)
plt.tight_layout();plt.show()

# Plot the conditional 12-month flow and stock outlook. The statistical stock
# path and raw accumulated flow identity need not coincide: reconciliation
# measures predicted stock change minus forecast net supply, not an observed flow.

outlook=tables['us_forecast_12m']
display(outlook[['month','production_kbd','demand_kbd','imports_kbd','exports_kbd','transfers_kbd',
    'adjustments_kbd','cdu_capacity_kbd','utilization_ratio','balance_kbd','stock_kb','identity_stock_unclipped_kb','model_reconciliation_kbd']])
fig,axes=plt.subplots(2,2,figsize=(14,8))
for c,label in [('supply_kbd','Supply including receipts/adjustments/transfers'),('total_demand_kbd','Refinery input + exports + direct use')]:
    axes[0,0].plot(outlook.month,outlook[c],label=label)
axes[0,0].set(title='U.S. monthly flow outlook',ylabel='Thousand barrels/day');axes[0,0].legend(fontsize=7)
axes[0,1].bar(outlook.month,outlook.balance_kbd,width=20)
axes[0,1].axhline(0,color='black');axes[0,1].set(title='Projected net flow balance',ylabel='Thousand barrels/day')
recent=tables['us_monthly_model'].tail(24)
axes[1,0].plot(recent.month,recent.stock_kb/1000,label='Actual')
axes[1,0].plot(outlook.month,outlook.stock_kb/1000,label=metadata['year_ahead_model'])
axes[1,0].plot(outlook.month,outlook.identity_stock_unclipped_kb/1000,label='Raw flow identity',linestyle='--')
axes[1,0].plot(outlook.month,np.repeat(recent.stock_kb.iloc[-1]/1000,12),label='Persistence baseline',linestyle=':')
axes[1,0].set(title='Conditional stock paths: no calibrated intervals',ylabel='Million barrels');axes[1,0].legend(fontsize=7)
for p,g in tables['padd_forecast_12m'].groupby('padd'):
    axes[1,1].plot(g.month,g.stock_kb/1000,label=f'PADD {p}')
axes[1,1].set(title='Regional commercial stock paths',ylabel='Million barrels');axes[1,1].legend(fontsize=8)
for ax in axes.flat:ax.tick_params(axis='x',rotation=30)
plt.tight_layout();plt.show()
fig,axes=plt.subplots(1,2,figsize=(14,4))
# Convert utilization ratios to percent for display; capacity is plotted
# separately to expose the constant-capacity forecast assumption.
axes[0].plot(tables['us_monthly_model'].tail(36).month,tables['us_monthly_model'].tail(36).utilization_ratio*100,label='Observed')
axes[0].plot(outlook.month,outlook.utilization_ratio*100,label='Forecast')
axes[0].set(title='Crude input / CDU capacity',ylabel='Utilization (%)');axes[0].legend()
for p,g in tables['padd_forecast_12m'].groupby('padd'):
    axes[1].plot(g.month,g.cdu_capacity_kbd,label=f'PADD {p}')
axes[1].set(title='Capacity held constant at forecast origin',ylabel='Thousand barrels/day');axes[1].legend()
plt.tight_layout();plt.show()
# Validate that capacity stays fixed within each PADD's forward path and
# that national demand equals aggregate utilization times aggregate capacity.
assert tables['padd_forecast_12m'].groupby('padd').cdu_capacity_kbd.nunique().eq(1).all()
np.testing.assert_allclose(outlook.demand_kbd,outlook.utilization_ratio*outlook.cdu_capacity_kbd)

# Summarize national forecasts on matched months and origins.
one = tables['us_predictions'].query("split == 'holdout'").copy()
one_name = metadata.get('one_month_model', 'selected_padd_models')
long_name = metadata['year_ahead_model']
long = tables['us_recursive_holdout_predictions'].copy()
history_stocks = tables['us_monthly_model'].set_index('month').stock_kb
base = long[long.model.eq(long_name)].copy()
assert not base.duplicated(['origin_month', 'month']).any()
year_ago = base.copy()
lag_dates = year_ago.month - pd.DateOffset(years=1)
assert (lag_dates <= year_ago.origin_month).all()
year_ago['predicted_kb'] = history_stocks.reindex(lag_dates).to_numpy()
assert year_ago.predicted_kb.notna().all()
year_ago['model'] = 'seasonal_naive'
long = pd.concat([long[~long.model.eq('seasonal_naive')], year_ago], ignore_index=True)

def summary_table(frame):
    rows = []
    reference = frame[frame.model.eq('persistence')].set_index(['origin_month', 'month'])
    for name, group in frame.groupby('model'):
        aligned = group.set_index(['origin_month', 'month']).sort_index()
        pd.testing.assert_index_equal(aligned.index, reference.sort_index().index)
        np.testing.assert_allclose(aligned.actual_kb, reference.sort_index().actual_kb)
        rows.append(dict(model=name, **score(group.actual_kb, group.predicted_kb),
                         bias_kb=(group.predicted_kb-group.actual_kb).mean()))
    result = pd.DataFrame(rows).set_index('model').sort_values('mae_kb')
    for benchmark, label in [('persistence','MAE improvement vs persistence (%)'),
                             ('seasonal_naive','MAE improvement vs year ago (%)')]:
        result[label] = 100*(1-result.mae_kb/result.loc[benchmark,'mae_kb'])
    return result

short_scores, long_scores = summary_table(one), summary_table(long)
report = short_scores.copy()
for column in ['mae_kb','rmse_kb','bias_kb']:
    report[column] /= 1000
display(report.rename(columns={'mae_kb':'MAE (million bbl)', 'rmse_kb':'RMSE (million bbl)',
    'bias_kb':'Bias (million bbl)', 'r2':'R²', 'n':'Forecast observations'}))

report = long_scores.copy()
for column in ['mae_kb','rmse_kb','bias_kb']:
    report[column] /= 1000
display(report.rename(columns={'mae_kb':'MAE (million bbl)', 'rmse_kb':'RMSE (million bbl)',
    'bias_kb':'Bias (million bbl)', 'r2':'R²', 'n':'Forecast observations'}))

selected = one[one.model.eq(one_name)].sort_values('month').copy()
selected['error'] = (selected.predicted_kb-selected.actual_kb)/1000
fig, axes = plt.subplots(1,3,figsize=(16,4))
axes[0].bar(selected.month,selected.error,width=20)
axes[0].set(title='One-month inventory forecast errors',ylabel='Forecast − actual (million bbl)')
axes[1].hist(selected.error,bins=9,edgecolor='white')
axes[1].set(title='Distribution of forecast errors',xlabel='Forecast − actual (million bbl)',ylabel='Months')
axes[2].scatter(selected.predicted_kb/1000,selected.error)
axes[2].set(title='Errors versus forecast stock level',xlabel='Forecast stocks (million bbl)',ylabel='Forecast − actual (million bbl)')
for ax in [axes[0],axes[2]]: ax.axhline(0,color='black',linewidth=.8)
axes[1].axvline(0,color='black',linewidth=.8)
fig.autofmt_xdate(); plt.tight_layout(); plt.show()

fig, ax = plt.subplots(figsize=(11,4))
for name, group in long.groupby('model'):
    errors = group.assign(error=(group.predicted_kb-group.actual_kb).abs()/1000)
    errors.groupby('horizon').error.mean().plot(ax=ax,marker='o',label=name)
ax.set(title='Inventory error by forecast horizon',xlabel='Months ahead',ylabel='MAE (million bbl)')
ax.legend(fontsize=8); plt.tight_layout(); plt.show()

# Check the model's stock scope and accounting identities, then perturb the
# last target observation to verify target-month isolation. Finally require
# the national forecast to equal the sum of the five regional forecasts.

assert not any('spr' in c.lower() for c in panel.columns)
np.testing.assert_allclose(panel.demand_kbd, panel.utilization_ratio * panel.cdu_capacity_kbd)
check_history=panel[panel.padd.eq(3)].iloc[:40].copy()
before=supervised(check_history)
# Perturb the final target month's stocks and flows substantially to test that
# its own inputs cannot see target-month information.
check_history.loc[check_history.index[-1],FLOWS+['stock_kb']]+=100000
changed=supervised(check_history)
# Only target labels may change after perturbing the final observation. All
# predictors, origin dates, and flow-identity forecasts must remain identical.
feature_columns=[c for c in before if c not in ['actual_kb','delta_kb']]
pd.testing.assert_frame_equal(before[feature_columns],changed[feature_columns])
np.testing.assert_allclose(tables['us_forecast_12m'].stock_kb,
    tables['padd_forecast_12m'].groupby('month').stock_kb.sum())
print('Source integrity, stock scope, target-month isolation, and aggregation checks passed.')

Source integrity, stock scope, target-month isolation, and aggregation checks passed.

# Inspect run configuration and saved artifacts, and validate output shape: no
# duplicate forecasts, all five PADDs each month, and a full 12-month outlook.

display(pd.Series(metadata,name='Run configuration'))
assert not tables['padd_predictions'].duplicated(['padd','month','model','split']).any()
assert tables['padd_forecast_12m'].groupby('month').padd.nunique().eq(5).all()
assert len(outlook)==12
assert len(tables['latest_forecast'])==6
assert not ((tables['us_recursive_cv_predictions'].month>=pd.Timestamp(metadata['holdout_start']))).any()
print('Notebook checks passed. Saved result files:')
for p in sorted(OUTPUT.iterdir()):print(p.name)

product                  Commercial crude oil; excludes strategic reserves
covid_exclusion          {'start': '2020-03-01', 'end': '2021-03-01', '...
stock_units                                  thousand barrels at month end
flow_units                                        thousand barrels per day
equation                 S[t] = S[t-1] + days[t] * (field production[t]...
data_start                                                      2007-01-01
data_end                                                        2026-06-01
cv                       10 expanding folds, 12 eligible test months ea...
holdout_start                                                   2024-07-01
forecast_timing          One month after latest observed EIA month; con...
selection                GridSearchCV per PADD/model on development dat...
hyperparameter_search    {'method': 'GridSearchCV', 'folds': 10, 'scori...
one_month_model                                            seasonal_change
year_ahead_model                                         constrained_level
horizon_selection        Lowest national MAE across six disjoint 12-mon...
flow_forecast            Last 60 non-COVID observations, daily rates, l...
forecast_stock_floor                                                     0
uncertainty              Point forecasts only; no calibrated prediction...
recursive_validation     13 overlapping 12-month paths entirely in fina...
versions                 {'python': '3.13.9', 'numpy': '2.3.5', 'pandas...
sources                  [https://www.eia.gov/dnav/pet/pet_sum_snd_d_r1...
Name: Run configuration, dtype: object

Notebook checks passed. Saved result files:
all_fitted_models.joblib
best_parameters.csv
candidate_addition_comparison.csv
candidate_models.json
constrained_model_coefficients.csv
fit_warnings.csv
fitted_horizon_models.joblib
fitted_models.joblib
fold_metrics.csv
grid_search_results.csv
latest_forecast.csv
model_metadata.json
national_model_cv_metrics.csv
padd_forecast_12m.csv
padd_model_coefficients.csv
padd_model_metrics.csv
padd_monthly_model.csv
padd_predictions.csv
recursive_cv_predictions.csv
recursive_holdout_metrics.csv
recursive_holdout_predictions.csv
refit_comparison.csv
regional_cv_winners.csv
regional_error_review.csv
selected_models.csv
training_sample_audit.csv
us_forecast_12m.csv
us_model_metrics.csv
us_monthly_model.csv
us_predictions.csv
us_recursive_cv_predictions.csv
us_recursive_holdout_predictions.csv
us_recursive_horizon_metrics.csv
year_ahead_model_cv_metrics.csv
year_ahead_refit_comparison.csv

	component	series_id	title
20	production_kb	MCRFPP31	Gulf Coast (PADD 3) Field Production of Crude ...
21	demand_kb	MCRRIP31	Gulf Coast (PADD 3) Refinery and Blender Net I...
22	imports_kb	MCRIMP31	Gulf Coast (PADD 3) Imports of Crude Oil (Thou...
23	exports_kb	MCREXP31	Gulf Coast (PADD 3) Exports of Crude Oil (Thou...
24	net_receipts_kb	MCRNRP31	Gulf Coast (PADD 3) Net Receipts by Pipeline, ...
25	adjustments_kb	MCRUA_R30_1	Gulf Coast (PADD 3) Supply Adjustment of Crude...
26	transfers_kb	M_EPC0_TVP_R30_MBBL	Gulf Coast (PADD 3) Transfers to Crude Oil Sup...
27	direct_use_kb	MCRUPP31	Gulf Coast (PADD 3) Product Supplied of Crude ...
28	commercial_stock_kb	MCESTP31	Gulf Coast (PADD 3) Ending Stocks excluding SP...
29	cdu_capacity_kbd	MOCLEP32	Gulf Coast (PADD 3) Operable Crude Oil Distill...

component	exports_kbd_assumed_zero	transfers_kbd_assumed_zero	net_receipts_kbd_assumed_zero
padd
1	21	180	0
2	0	180	0
3	63	180	0
4	98	180	0
5	103	180	43

component	padd	padd_name	month	cdu_capacity_kbd	demand_kbd	utilization_ratio
233	1	East Coast	2026-06-01	928.00	727.00	0.78
467	2	Midwest	2026-06-01	4,283.00	4,291.17	1.00
701	3	Gulf Coast	2026-06-01	9,888.00	9,525.30	0.96
935	4	Rocky Mountain	2026-06-01	653.00	656.33	1.01
1169	5	West Coast	2026-06-01	2,275.00	2,000.47	0.88

	incomplete_balance_mae_kbd	complete_balance_mae_kbd	max_complete_residual_kbd
padd
1	161.80	0.01	0.04
2	825.28	0.02	1.94
3	1,198.28	91.67	1,262.68
4	462.47	0.01	0.19
5	146.96	0.01	0.04

component	month	padd	accounting_residual_kbd
471	2007-04-01	3	-26.33
472	2007-05-01	3	-28.32
475	2007-08-01	3	-5.39
476	2007-09-01	3	-78.87
477	2007-10-01	3	-42.77
...	...	...	...
697	2026-02-01	3	-8.18
698	2026-03-01	3	19.90
699	2026-04-01	3	676.73
700	2026-05-01	3	1,262.68
701	2026-06-01	3	1,106.47

U.S. crude oil supply, demand, and stocks¶

Data refresh — September 16, 2026¶

Data¶

Sources and units¶

Model¶

Forecast inputs¶

Tested models¶

Training¶

Results¶

Regional results¶

Twelve-month outlook¶

Forecast performance¶

One-month national results¶

Twelve-month path results¶

Interpretation¶

Reproduce the results¶

	start	end	months
padd
1	2007-01-01	2026-06-01	234
2	2007-01-01	2026-06-01	234
3	2007-01-01	2026-06-01	234
4	2007-01-01	2026-06-01	234
5	2007-01-01	2026-06-01	234

	origin_month	month	stock_lag1	lag_production_kbd	lag_demand_kbd	forecast_flow_identity	actual_kb
217	2026-01-01	2026-02-01	224,301.00	9,849.23	-8,907.68	243,036.16	249,711.00
218	2026-02-01	2026-03-01	249,711.00	10,211.14	-8,733.11	260,609.30	260,536.00
219	2026-03-01	2026-04-01	260,536.00	10,190.74	-9,276.29	269,918.32	258,758.00
220	2026-04-01	2026-05-01	258,758.00	10,439.80	-9,360.27	261,554.91	244,403.00
221	2026-05-01	2026-06-01	244,403.00	10,221.13	-9,556.10	240,283.49	233,275.00

	model	param_grid
0	constrained_level	{}
1	ridge_change	{'ridge__alpha': [0.01, 0.1, 1, 10, 100, 1000]}
2	seasonal_ridge_change	{'ridge__alpha': [0.01, 0.1, 1, 10, 100, 1000]}
3	polynomial_ridge_change	{'polynomialfeatures__degree': [1, 2], 'ridge_...
4	spline_ridge_change	{'splinetransformer__n_knots': [3, 5, 7], 'spl...
5	random_forest_change	{'n_estimators': [100, 200], 'max_depth': [3, ...
6	xgboost_change	{'n_estimators': [100, 200], 'max_depth': [2, ...
7	neural_network_change	{'regressor__mlpregressor__hidden_layer_sizes'...
8	rolling_ridge_change	{'ridge__alpha': [0.01, 0.1, 1, 10, 100, 1000]}
9	damped_ridge_change	{'ridge__alpha': [0.01, 0.1, 1, 10, 100, 1000]}
10	ensemble_change	{'weights': [(1, 1, 1), (2, 1, 1), (1, 2, 1), ...
11	seasonal_residual_change	{'ridge__alpha': [0.1, 1, 10, 100, 1000]}
12	huber_change	{'regressor__huberregressor__epsilon': [1.35, ...

	fold	train_end	test_start	test_end	n
0	1	2013-05-01	2013-06-01	2014-05-01	12
1	2	2014-05-01	2014-06-01	2015-05-01	12
2	3	2015-05-01	2015-06-01	2016-05-01	12
3	4	2016-05-01	2016-06-01	2017-05-01	12
4	5	2017-05-01	2017-06-01	2018-05-01	12
5	6	2018-05-01	2018-06-01	2019-05-01	12
6	7	2019-05-01	2019-06-01	2021-06-01	12
7	8	2020-02-01	2021-07-01	2022-06-01	12
8	9	2022-06-01	2022-07-01	2023-06-01	12
9	10	2023-06-01	2023-07-01	2024-06-01	12

	padd	selected_model	cv_mae_kb
0	1	neural_network_change	1,046.18
1	2	damped_ridge_change	4,258.24
2	3	seasonal_change	6,222.72
3	4	seasonal_residual_change	738.58
4	5	spline_ridge_change	1,646.08

	model	mae_kb	rmse_kb	r2	n
11	seasonal_change	8,190.01	10,450.37	0.95	120
13	seasonal_residual_change	8,232.03	10,340.39	0.95	120
15	selected_padd_models	8,440.53	10,453.34	0.95	120
17	xgboost_change	9,051.74	11,385.25	0.94	120
8	random_forest_change	9,069.75	11,319.01	0.94	120
16	spline_ridge_change	9,318.61	11,109.33	0.95	120
2	ensemble_change	9,384.75	11,592.29	0.94	120
7	polynomial_ridge_change	9,903.90	12,047.77	0.94	120
4	huber_change	10,014.32	12,339.43	0.93	120
5	neural_network_change	10,044.17	12,263.98	0.93	120
10	rolling_ridge_change	10,098.67	12,533.07	0.93	120
1	damped_ridge_change	10,152.76	12,435.16	0.93	120
14	seasonal_ridge_change	10,162.22	12,429.43	0.93	120
9	ridge_change	10,564.22	12,964.16	0.93	120
6	persistence	10,818.33	13,133.21	0.92	120
0	constrained_level	11,103.92	14,235.30	0.91	120
3	forecast_flow_identity	11,266.31	14,790.30	0.90	120
12	seasonal_naive	39,776.22	50,836.94	-0.15	120

padd	1	2	3	4	5
model
constrained_level	1240.045865	4651.975993	8145.071871	887.980370	1720.386926
damped_ridge_change	1051.416827	4258.244729	7674.048251	785.185145	1728.138425
ensemble_change	1061.771003	4298.513596	7011.194897	773.741179	1756.959525
forecast_flow_identity	2306.896990	5010.608114	9005.547167	811.827747	2013.300947
huber_change	1077.731035	4309.789931	7591.271419	925.121139	1758.239635
neural_network_change	1046.178623	4304.095395	7646.017868	811.739484	1777.984012
persistence	1114.841667	4572.183333	7835.491667	811.408333	1886.983333
polynomial_ridge_change	1060.653663	4260.349419	7551.761366	809.733439	1773.842029
random_forest_change	1074.623520	4401.772969	6683.817144	744.434036	1763.691864
ridge_change	1074.791505	4374.943200	7845.659070	816.848524	1754.866718
rolling_ridge_change	1114.413595	4378.170011	7418.082881	805.171253	1779.949768
seasonal_change	1232.179167	4690.117083	6222.722083	757.695000	1809.377917
seasonal_naive	2099.583333	17976.983333	23205.116667	1372.358333	2907.075000
seasonal_residual_change	1113.166513	4399.248253	6251.620832	738.580768	1659.496021
seasonal_ridge_change	1060.653663	4260.349419	7599.706932	812.964942	1773.842029
spline_ridge_change	1088.705381	4437.362796	7222.400995	797.185309	1646.076458
xgboost_change	1073.722131	4359.164779	6568.209115	743.130029	1732.058887

	train_mae_kb	mae_kb	test_minus_train_kb
model
seasonal_residual_change	2,365.32	2,832.42	467.11
xgboost_change	1,468.34	2,895.26	1,426.92
random_forest_change	1,790.16	2,933.67	1,143.51
seasonal_change	2,635.56	2,942.42	306.86
ensemble_change	2,184.63	2,980.44	795.80
spline_ridge_change	2,251.17	3,038.35	787.18
polynomial_ridge_change	2,457.46	3,091.27	633.80
rolling_ridge_change	2,848.19	3,099.16	250.97
damped_ridge_change	2,529.27	3,099.41	570.14
seasonal_ridge_change	2,477.46	3,101.50	624.04
neural_network_change	2,211.56	3,117.20	905.64
huber_change	2,306.29	3,132.43	826.14
ridge_change	2,534.75	3,173.42	638.67
persistence	2,773.12	3,244.18	471.07
constrained_level	2,451.07	3,329.09	878.02
forecast_flow_identity	3,211.49	3,829.64	618.14
seasonal_naive	7,338.88	9,512.22	2,173.35

	padd	intercept_kb	stock_lag1	lag_production_kbd	lag_demand_kbd	lag_imports_kbd	lag_exports_kbd	lag_net_receipts_kbd	lag_adjustments_kbd	lag_transfers_kbd
0	1	1,600.17	0.78	0.00	0.00	0.81	0.00	1.85	4.21	0.00
1	2	13,336.58	0.92	13.08	12.13	9.12	6.41	10.35	6.09	10.98
2	3	35,906.85	0.90	4.68	4.70	0.26	5.11	1.60	5.88	2.25
3	4	-439.18	0.87	5.23	1.10	5.23	2.57	3.77	5.98	11.97
4	5	22,636.18	0.48	13.95	8.25	6.56	16.32	10.99	8.38	0.00

	model	split	mae_kb	rmse_kb	r2	n
23	seasonal_change	holdout	6,781.31	9,036.61	0.44	24
31	selected_padd_models	holdout	7,250.51	9,222.26	0.42	24
27	seasonal_residual_change	holdout	7,339.84	9,248.33	0.42	24
17	random_forest_change	holdout	7,862.99	9,562.93	0.38	24
33	spline_ridge_change	holdout	8,336.50	10,238.91	0.29	24
5	ensemble_change	holdout	8,574.38	10,287.70	0.28	24
35	xgboost_change	holdout	8,896.36	11,143.55	0.15	24
11	neural_network_change	holdout	9,011.52	10,533.61	0.24	24
29	seasonal_ridge_change	holdout	9,225.54	11,265.22	0.14	24
3	damped_ridge_change	holdout	9,436.51	11,551.45	0.09	24
15	polynomial_ridge_change	holdout	9,509.31	11,410.50	0.11	24
13	persistence	holdout	9,630.54	12,194.98	-0.01	24
21	rolling_ridge_change	holdout	9,952.10	12,379.51	-0.04	24
9	huber_change	holdout	10,015.02	11,610.71	0.08	24
19	ridge_change	holdout	10,045.40	12,095.97	0.00	24
1	constrained_level	holdout	10,755.87	12,128.02	-0.00	24
25	seasonal_naive	holdout	10,990.08	13,923.84	-0.32	24
7	forecast_flow_identity	holdout	11,878.49	14,265.31	-0.39	24

	padd	model	params
0	1	ridge_change	{"ridge__alpha": 10}
1	1	seasonal_ridge_change	{"ridge__alpha": 100}
2	1	polynomial_ridge_change	{"polynomialfeatures__degree": 1, "ridge__alpha": 100}
3	1	spline_ridge_change	{"ridge__alpha": 1000, "splinetransformer__degree": 2, "splinetransformer__n_knots": 3}
4	1	random_forest_change	{"max_depth": 5, "min_samples_leaf": 12, "n_estimators": 200}
5	1	xgboost_change	{"learning_rate": 0.01, "max_depth": 2, "n_estimators": 200, "reg_lambda": 1}
6	1	neural_network_change	{"regressor__mlpregressor__alpha": 10, "regressor__mlpregressor__hidden_layer_sizes": [32, 16]}
7	1	rolling_ridge_change	{"ridge__alpha": 1000}
8	1	damped_ridge_change	{"ridge__alpha": 10}
9	1	ensemble_change	{"weights": [2, 1, 1]}
10	1	seasonal_residual_change	{"ridge__alpha": 100}
11	1	huber_change	{"regressor__huberregressor__alpha": 10, "regressor__huberregressor__epsilon": 1.75}
12	2	ridge_change	{"ridge__alpha": 100}
13	2	seasonal_ridge_change	{"ridge__alpha": 100}
14	2	polynomial_ridge_change	{"polynomialfeatures__degree": 1, "ridge__alpha": 100}
15	2	spline_ridge_change	{"ridge__alpha": 1000, "splinetransformer__degree": 2, "splinetransformer__n_knots": 3}
16	2	random_forest_change	{"max_depth": null, "min_samples_leaf": 12, "n_estimators": 200}
17	2	xgboost_change	{"learning_rate": 0.05, "max_depth": 4, "n_estimators": 100, "reg_lambda": 30}
18	2	neural_network_change	{"regressor__mlpregressor__alpha": 10, "regressor__mlpregressor__hidden_layer_sizes": [32, 16]}
19	2	rolling_ridge_change	{"ridge__alpha": 100}
20	2	damped_ridge_change	{"ridge__alpha": 10}
21	2	ensemble_change	{"weights": [2, 1, 1]}
22	2	seasonal_residual_change	{"ridge__alpha": 10}
23	2	huber_change	{"regressor__huberregressor__alpha": 10, "regressor__huberregressor__epsilon": 1.75}
24	3	ridge_change	{"ridge__alpha": 100}
25	3	seasonal_ridge_change	{"ridge__alpha": 100}
26	3	polynomial_ridge_change	{"polynomialfeatures__degree": 2, "ridge__alpha": 1000}
27	3	spline_ridge_change	{"ridge__alpha": 100, "splinetransformer__degree": 2, "splinetransformer__n_knots": 3}
28	3	random_forest_change	{"max_depth": 5, "min_samples_leaf": 5, "n_estimators": 200}
29	3	xgboost_change	{"learning_rate": 0.1, "max_depth": 4, "n_estimators": 200, "reg_lambda": 30}
30	3	neural_network_change	{"regressor__mlpregressor__alpha": 10, "regressor__mlpregressor__hidden_layer_sizes": [32, 16]}
31	3	rolling_ridge_change	{"ridge__alpha": 10}
32	3	damped_ridge_change	{"ridge__alpha": 100}
33	3	ensemble_change	{"weights": [1, 1, 2]}
34	3	seasonal_residual_change	{"ridge__alpha": 1000}
35	3	huber_change	{"regressor__huberregressor__alpha": 10, "regressor__huberregressor__epsilon": 1.35}
36	4	ridge_change	{"ridge__alpha": 1000}
37	4	seasonal_ridge_change	{"ridge__alpha": 1000}
38	4	polynomial_ridge_change	{"polynomialfeatures__degree": 2, "ridge__alpha": 1000}
39	4	spline_ridge_change	{"ridge__alpha": 1000, "splinetransformer__degree": 2, "splinetransformer__n_knots": 3}

	padd	model	split	mae_kb	rmse_kb	r2	n
1	1	constrained_level	holdout	591.62	776.97	-0.15	24
13	1	persistence	holdout	626.88	792.38	-0.20	24
23	1	seasonal_change	holdout	579.67	772.86	-0.14	24
37	2	constrained_level	holdout	3,333.58	4,041.28	-0.41	24
49	2	persistence	holdout	2,859.67	3,670.33	-0.17	24
59	2	seasonal_change	holdout	2,879.95	3,729.80	-0.20	24
73	3	constrained_level	holdout	7,368.32	9,095.45	-0.05	24
85	3	persistence	holdout	7,646.62	9,277.82	-0.09	24
95	3	seasonal_change	holdout	5,950.43	7,639.34	0.26	24
109	4	constrained_level	holdout	664.08	790.47	0.42	24
121	4	persistence	holdout	681.88	821.33	0.38	24
131	4	seasonal_change	holdout	556.01	815.43	0.39	24
145	5	constrained_level	holdout	1,916.97	2,498.38	-0.39	24
157	5	persistence	holdout	1,964.42	2,519.89	-0.41	24
167	5	seasonal_change	holdout	1,971.54	2,545.46	-0.44	24

	total_target_months	eligible_stock_training_months	evaluation_months
padd
1	222	197	209
2	222	197	209
3	222	197	209
4	222	197	209
5	222	197	209

	region	previous_mae_kb	refit_mae_kb	improvement_pct	test_months
0	U.S.	6,556.73	6,781.31	-3.43	24
1	PADD 1	676.88	579.67	14.36	24
2	PADD 2	3,143.41	2,879.95	8.38	24
3	PADD 3	5,740.99	5,950.43	-3.65	24
4	PADD 4	510.71	556.01	-8.87	24
5	PADD 5	2,086.72	1,971.54	5.52	24

	candidate_set	period	mae_kb
0	Original candidates, COVID excluded	cv	8,190.01
1	Original candidates, COVID excluded	holdout	6,781.31
2	Expanded candidates, COVID excluded	cv	8,190.01
3	Expanded candidates, COVID excluded	holdout	6,781.31

	region	previous_mae_kb	refit_mae_kb	improvement_pct	forecast_observations
0	U.S.	9,111.59	9,055.39	0.62	156
1	PADD 1	875.61	671.05	23.36	156
2	PADD 2	4,033.29	3,004.13	25.52	156
3	PADD 3	7,848.24	7,194.17	8.33	156
4	PADD 4	667.23	840.92	-26.03	156
5	PADD 5	2,204.30	1,382.37	37.29	156

	PADD	Average error (kb)	Error / average stocks (%)	Improvement over unchanged stocks (%)	Average bias (kb)	Worst month	Largest miss (kb)
0	1	579.67	7.12	7.53	-136.53	2024-07-01	2,265.20
1	2	2,879.95	2.72	-0.71	224.34	2026-05-01	9,159.40
2	3	5,950.43	2.49	22.18	599.34	2026-02-01	17,923.00
3	4	556.01	2.33	18.46	-12.15	2026-05-01	2,233.80
4	5	1,971.54	4.23	-0.36	131.87	2024-07-01	5,824.20

	model	mae_kb	rmse_kb	r2	n
0	constrained_level	28,781.81	38,347.05	-0.51	72
1	xgboost_change	28,919.88	39,424.45	-0.60	72
2	damped_ridge_change	29,341.30	37,692.35	-0.46	72
3	random_forest_change	29,822.63	42,305.12	-0.84	72
4	ensemble_change	29,866.45	39,904.73	-0.64	72
5	persistence	29,870.10	36,728.66	-0.39	72
6	seasonal_residual_change	30,349.95	40,776.48	-0.71	72
7	seasonal_change	30,562.80	42,345.68	-0.84	72
8	neural_network_change	30,788.51	41,904.48	-0.81	72
9	seasonal_ridge_change	31,776.08	41,092.61	-0.74	72
10	ridge_change	32,365.60	41,160.37	-0.74	72
11	spline_ridge_change	32,811.55	41,439.32	-0.77	72
12	rolling_ridge_change	33,019.08	43,972.66	-0.99	72
13	polynomial_ridge_change	33,082.74	42,960.75	-0.90	72
14	huber_change	34,169.66	44,236.48	-1.01	72
15	seasonal_naive	40,506.54	49,476.85	-1.52	72
16	forecast_flow_identity	68,966.32	86,352.95	-6.67	72

	first_target	last_target
origin_month
2016-06-01	2016-07-01	2017-06-01
2017-06-01	2017-07-01	2018-06-01
2018-06-01	2018-07-01	2019-06-01
2021-06-01	2021-07-01	2022-06-01
2022-06-01	2022-07-01	2023-06-01
2023-06-01	2023-07-01	2024-06-01

model	constrained_level	persistence	seasonal_change
horizon
1	7,795.34	7,190.69	3,954.18
2	8,871.62	10,117.62	6,913.70
3	9,625.49	10,605.85	8,737.63
4	9,610.92	12,202.69	9,926.27
5	8,631.78	12,647.62	11,138.46
6	7,784.86	11,781.00	12,218.46
7	8,533.26	12,079.77	11,995.90
8	9,323.87	12,631.77	12,896.13
9	10,022.13	11,347.77	14,962.32
10	9,838.57	11,566.15	16,221.69
11	8,981.00	10,312.23	14,879.82
12	9,645.86	8,948.62	11,916.30

	month	production_kbd	demand_kbd	imports_kbd	exports_kbd	transfers_kbd	adjustments_kbd	cdu_capacity_kbd	utilization_ratio	balance_kbd	stock_kb	identity_stock_unclipped_kb	model_reconciliation_kbd
0	2026-07-01	14,241.61	16,880.64	6,514.40	4,375.73	660.18	-34.30	18,027.00	0.94	125.51	402,441.45	408,011.76	-179.69
1	2026-08-01	14,323.07	16,915.37	6,399.69	4,551.10	725.39	1.74	18,027.00	0.94	-16.58	407,121.62	407,497.70	167.56
2	2026-09-01	14,322.92	16,436.70	6,432.55	4,469.70	775.57	-372.50	18,027.00	0.91	252.12	410,311.82	415,061.35	-145.78
3	2026-10-01	14,555.87	15,956.04	6,084.56	4,598.81	800.83	-200.75	18,027.00	0.89	685.65	412,584.02	436,316.35	-612.35
4	2026-11-01	14,599.41	16,652.84	6,340.41	4,624.04	713.67	-124.50	18,027.00	0.92	252.11	418,198.74	443,879.74	-64.96
5	2026-12-01	14,501.79	16,725.74	6,330.31	4,687.29	703.69	-64.41	18,027.00	0.93	58.33	420,626.06	445,688.02	19.97
6	2027-01-01	14,199.84	15,993.16	6,454.74	4,473.16	733.89	-323.94	18,027.00	0.89	598.23	422,217.11	464,233.08	-546.90
7	2027-02-01	14,430.85	15,725.85	6,332.03	4,873.69	781.21	-70.78	18,027.00	0.87	873.76	425,118.40	488,698.47	-770.15
8	2027-03-01	14,617.63	16,309.55	6,204.73	4,755.97	774.71	-63.35	18,027.00	0.90	468.21	430,310.96	503,212.91	-300.71
9	2027-04-01	14,661.55	16,353.71	6,196.55	4,974.10	746.39	-54.83	18,027.00	0.91	221.84	431,980.62	509,868.05	-166.18
10	2027-05-01	14,586.15	16,963.55	6,343.89	4,903.26	821.36	-174.17	18,027.00	0.94	-289.58	433,010.61	500,890.94	322.81
11	2027-06-01	14,661.15	17,218.61	6,343.21	4,766.78	747.89	-153.16	18,027.00	0.96	-386.30	430,462.99	489,302.00	301.38

	MAE (million bbl)	RMSE (million bbl)	R²	Forecast observations	Bias (million bbl)	MAE improvement vs persistence (%)	MAE improvement vs year ago (%)
model
seasonal_change	6.78	9.04	0.44	24	0.81	29.59	38.30
selected_padd_models	7.25	9.22	0.42	24	2.28	24.71	34.03
seasonal_residual_change	7.34	9.25	0.42	24	2.76	23.79	33.21
random_forest_change	7.86	9.56	0.38	24	3.66	18.35	28.45
spline_ridge_change	8.34	10.24	0.29	24	3.22	13.44	24.15
ensemble_change	8.57	10.29	0.28	24	3.21	10.97	21.98
xgboost_change	8.90	11.14	0.15	24	5.58	7.62	19.05
neural_network_change	9.01	10.53	0.24	24	4.12	6.43	18.00
seasonal_ridge_change	9.23	11.27	0.14	24	2.80	4.21	16.06
damped_ridge_change	9.44	11.55	0.09	24	2.51	2.01	14.14
polynomial_ridge_change	9.51	11.41	0.11	24	2.85	1.26	13.47
persistence	9.63	12.19	-0.01	24	1.52	0.00	12.37
rolling_ridge_change	9.95	12.38	-0.04	24	5.40	-3.34	9.44
huber_change	10.02	11.61	0.08	24	5.29	-3.99	8.87
ridge_change	10.05	12.10	0.00	24	2.92	-4.31	8.60
constrained_level	10.76	12.13	-0.00	24	5.85	-11.68	2.13
seasonal_naive	10.99	13.92	-0.32	24	7.55	-14.12	0.00
forecast_flow_identity	11.88	14.27	-0.39	24	-0.18	-23.34	-8.08

	MAE (million bbl)	RMSE (million bbl)	R²	Forecast observations	Bias (million bbl)	MAE improvement vs persistence (%)	MAE improvement vs year ago (%)
model
constrained_level	9.06	11.06	-0.08	156	1.53	17.32	25.41
persistence	10.95	13.47	-0.60	156	2.30	0.00	9.78
seasonal_change	11.31	13.85	-0.69	156	-1.88	-3.29	6.81
seasonal_naive	12.14	15.39	-1.08	156	10.15	-10.84	0.00