The main script

The orchestrator is intentionally minimal. It source()s each step in order and nothing else:

source("00_packages.R")
source("01_config.R")
source("02_market_data.R")
source("03_features.R")
source("04_news_sentiment.R")
source("05_assemble.R")
source("06_train_test_split.R")
source("07_models_brent.R")
source("08_gam_gasoline.R")
source("09_counterfactual.R")
source("10_accuracy.R")
source("11_plots.R")
source("12_export.R")
source("13_diagnostic_report.R")
source("14_compact_diagnostic.R")

Configuration (01_config.R)

All parameters are set in one place. This is the entire configuration for the system – dates, tickers, API keys, and output paths:

EIA_API_KEY  <- Sys.getenv("EIA_KEY")
FRED_API_KEY <- Sys.getenv("FRED_API_KEY")

BRENT_TICKER <- "BZ=F"
WTI_TICKER   <- "CL=F"

LOOKBACK_DAYS    <- 540      # 18 months of training data
FORECAST_HORIZON <- 270      # 9 months ahead

DATE_END   <- Sys.Date()
DATE_START <- DATE_END - LOOKBACK_DAYS

Futures curve data

The WTI futures curve gives the market’s forward expectation for oil prices across 12 contract months. The system fetches each contract individually from Yahoo Finance and computes curve metrics:

WTI_CURVE_TICKERS <- c(
  "CLM26.NYM", "CLN26.NYM", "CLQ26.NYM", "CLU26.NYM",
  "CLV26.NYM", "CLX26.NYM", "CLZ26.NYM", "CLF27.NYM",
  "CLG27.NYM", "CLH27.NYM", "CLJ27.NYM", "CLK27.NYM"
)

The curve slope (front-month vs. back-month, annualized) tells you whether the market is in contango (upward-sloping, future prices higher than spot) or backwardation (downward-sloping, future prices below spot). This distinction matters for forecasting: a model that predicts rising prices in a steeply backwardated market is fighting the market’s consensus.

Technical indicators (03_features.R)

Seasonal demand indices

Oil demand has strong seasonal patterns. The system creates continuous indices rather than binary flags:

add_seasonal_demand <- function(df) {
  df |> mutate(
    # Peaks mid-summer (day ~172), zero in winter
    driving_season_idx = pmax(0, sin(2 * pi * (yday(date) - 80) / 365)),
    # Peaks mid-winter, zero in summer
    heating_demand_idx = pmax(0, -cos(2 * pi * yday(date) / 365)),
    # Binary: spring + fall turnaround windows
    refinery_maint = as.integer(month(date) %in% c(3, 4, 9, 10)),
    summer_peak    = as.integer(month(date) %in% 5:8),
    winter_heat    = as.integer(month(date) %in% c(11, 12, 1, 2))
  )
}

News sentiment (04_news_sentiment.R)

The system scrapes RSS feeds from major energy news sources, scores each article using keyword-based topic classification, and computes daily sentiment features across five channels:

A critical design choice: negative news enters the model quickly (0-2 day lag) while positive news enters slowly (5-21 day rolling average). This reflects the empirical asymmetry in how oil markets respond to shocks – supply disruptions spike prices within hours, while diplomatic resolutions take weeks to filter through.

5.1 ARIMA-X

ARIMA (AutoRegressive Integrated Moving Average) with exogenous regressors is the classical time-series workhorse. The model is fit on log-transformed Brent prices with d=1 (first differencing enforced) and the technical indicators as external regressors:

brent_ts <- ts(log(train$brent_close), frequency = 5)

arima_model <- auto.arima(
  brent_ts,
  xreg    = xreg_train,   # scaled feature matrix
  d       = 1,            # force first differencing
  stepwise    = TRUE,
  approximation = TRUE,
  max.p = 5, max.q = 3,
  max.P = 1, max.Q = 1,
  ic = "aicc"
)

Strengths: Captures short-term autocorrelation structure (today’s price predicts tomorrow’s). Ljung-Box test on residuals tells you if the model has captured all available serial correlation.

Weaknesses: Assumes the autocorrelation structure is stable over the entire training window. In regime-shifting markets, this is often false.

5.2 Prophet

Facebook’s Prophet decomposes the series into trend, weekly seasonality, yearly seasonality, and external regressors. It handles missing data and holidays automatically:

m_prophet <- prophet(
  seasonality.mode      = "multiplicative",
  daily.seasonality     = FALSE,
  weekly.seasonality    = TRUE,
  yearly.seasonality    = TRUE,
  changepoint.prior.scale = 0.15  # controls trend flexibility
)

Strengths: Handles structural breaks (changepoints) more gracefully than ARIMA. Multiplicative seasonality is appropriate for percentage-based price moves.

Weaknesses: Not designed for financial time series. The weekly “seasonality” in oil prices is really a day-of-week effect, not a true seasonal cycle.

5.3 LASSO

Least Absolute Shrinkage and Selection Operator regression performs automatic feature selection by penalizing coefficient magnitudes toward zero. Features that do not improve out-of-sample prediction are dropped entirely:

en_spec <- linear_reg(penalty = tune(), mixture = tune()) |>
  set_engine("glmnet")

# Time-series cross-validation (rolling origin)
ts_folds <- rolling_origin(
  train,
  initial    = round(nrow(train) * 0.7),
  assess     = 10,
  cumulative = FALSE,
  skip       = 9
)

Strengths: Handles many correlated features without manual selection. The regularization path reveals which features carry genuine predictive power versus noise.

Weaknesses: Linear model – cannot capture nonlinear interactions unless you engineer them explicitly.

Bias-penalized inverse-RMSE weighting

The three models are combined using inverse-RMSE weights computed from expanding-window cross-validation. In v17, a bias penalty was added: models that consistently under- or over-predict have their effective RMSE inflated, reducing their weight in the ensemble.

# Directional bias: fraction of over-predictions per model
over_pct_vec <- c(
  arima   = mean(xw_all$arima_pred   > xw_all$actual, na.rm = TRUE),
  prophet = mean(xw_all$prophet_pred > xw_all$actual, na.rm = TRUE),
  lasso   = mean(xw_all$lasso_pred   > xw_all$actual, na.rm = TRUE)
)

# Bias penalty: 1.0 (balanced) to 2.0 (completely one-sided)
bias_penalty_vec <- 1 + 2 * abs(over_pct_vec - 0.5)

# Penalized RMSE for weight computation
penalized_rmse_vec <- rmse_vec * bias_penalty_vec

# Inverse-penalized-RMSE weights
wts <- (1 / penalized_rmse_vec) / sum(1 / penalized_rmse_vec)

Why not equal weights?

Equal weighting (1/3 each) is a strong baseline that often beats complex weighting schemes. The bias-penalized approach adds value when models have asymmetric errors – if ARIMA consistently under-predicts by $5 while LASSO is unbiased, giving them equal weight bakes in a -$1.67 bias.

The under-prediction problem

The most persistent diagnostic finding across versions was systematic under-prediction: over_pct=0% (the model under-predicted every single backtest day). The root cause is that the ensemble averages three models trained on 18 months of history, which includes periods at lower prices. In a rising market, this historical anchoring creates a downward drag.

Spot-anchoring correction

The fix blends the raw ensemble forecast with a “no-change from spot” baseline. The blend weight decays exponentially over the forecast horizon using a calibrated half-life:

# Alpha = 1 at day 0 (pure spot), decays toward 0 (pure model)
fc_alpha <- 0.5 ^ (fc_days_out / best_halflife)

ensemble_future <- fc_alpha * spot_current +
                   (1 - fc_alpha) * ensemble_future

The half-life is calibrated by grid search over the expanding-window backtest, choosing the value that minimizes RMSE. A 15-day half-life means the model’s own signal dominates by about day 30.

Curve-consistent blending (v18)

When the ensemble forecasts a direction that contradicts a steep futures curve (e.g., predicting a rally in a -22% backwardated market), the forecast is blended with the curve-implied price path at longer horizons:

if (contradicts_curve) {
  curve_implied <- spot_current * (1 + curve_slope * (fc_days_out / 365))
  curve_weight  <- pmin(0.60, pmax(0, (fc_days_out - 30) / 240 * 0.60))
  ensemble_future <- (1 - curve_weight) * ensemble_future +
                      curve_weight * curve_implied
}

The weight ramps from 0 at day 30 to a maximum of 60% at 9 months. This prevents the model from predicting a 10% rally when the forward curve prices in a 22% decline.

AR(1) residual correction (v18)

If the backtest residuals show strong autocorrelation (ACF lag-1 > 0.3), there is persistent structure the ensemble has not captured. An AR(1) model fit on the residual series projects a decaying correction onto the forecast:

residuals_bt <- test$brent_close - ensemble_test
resid_ar     <- arima(residuals_bt, order = c(1, 0, 0))
resid_forecast <- predict(resid_ar, n.ahead = length(future_dates))$pred

# Cap at +/- 2 SD to prevent runaway extrapolation
max_correction <- 2 * sd(residuals_bt)
resid_forecast <- pmax(-max_correction, pmin(max_correction, resid_forecast))

ensemble_future <- ensemble_future + resid_forecast

GAM with holdout validation

A Generalized Additive Model (GAM) captures the nonlinear relationship between Brent and the spread:

spread_gam <- mgcv::gam(
  spread ~ s(dl_brent_lag, k = 8) +   # nonlinear pass-through
    s(driving_season_idx, k = 5) +     # seasonal cycle
    lag_spread +                        # autoregressive component
    refinery_stress +                   # CA refinery outage intensity
    sent_refinery +                     # refinery sentiment
    sent_ca_regulatory +                # CA regulatory sentiment
    summer_peak + winter_heat,
  data   = spread_df,
  method = "REML"
)

The overfitting problem and OLS fallback (v18)

In v17, the GAM achieved in-sample adjusted R-squared of 0.860 but holdout R-squared of -15.049 – worse than predicting the mean. The GAM was fitting noise in the training data despite reduced basis dimensions (k=8, k=5).

The v18 fix adds automatic holdout validation: after fitting the GAM, hold out the last 15% of training data and check R-squared. If negative, fall back to a simpler OLS model:

gam_val_n    <- min(30, floor(nrow(spread_df) * 0.15))
gam_val_df   <- tail(spread_df, gam_val_n)
gam_train_df <- head(spread_df, nrow(spread_df) - gam_val_n)

gam_val_model <- mgcv::gam(/* same formula */, data = gam_train_df)
gam_val_pred  <- predict(gam_val_model, newdata = gam_val_df)
gam_val_r2    <- 1 - var(gam_val_df$spread - gam_val_pred) /
                     var(gam_val_df$spread)

if (gam_val_r2 < 0) {
  # Fall back to OLS
  spread_gam <- lm(
    spread ~ dl_brent_lag + driving_season_idx + lag_spread +
      refinery_stress + summer_peak + winter_heat + quarter,
    data = spread_df
  )
}

Reading the diagnostic

The improvement log

Every version’s header documents the specific diagnostic finding that motivated each change: