Prediction Intervals with PyTorch: Quantile Regression Using Neural Networks
Part 3 of the Quantile Regression Series
In Part 1, we implemented quantile regression using CatBoost. In Part 2, we took a closer look at quantile regression and pinball loss. Now let’s do it in PyTorch.
We’ll train a small neural network, visualize its prediction intervals, and compare it with CatBoost.
We’ll use the same Yahoo Finance SPY example as Part 1, predicting P₁₀, P₅₀, and P₉₀ of the next trading session’s maximum upside from its open. We’ll keep the same five features, 75/25 chronological split, and CatBoost settings so the comparison is on identical test dates.
Prediction target: Yₜ = 100 × (Highₜ₊₁ / Openₜ₊₁ − 1). The P₁₀–P₉₀ band is a nominal 80% interval for tomorrow’s maximum upside, not the day’s entire high–low range or closing price.
Download real SPY data from YFinance
Install once, if needed, in Jupyter: %pip install -U yfinance torch catboost scikit-learn pandas numpy matplotlib
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import torch
from torch import nn
from sklearn.preprocessing import StandardScaler
from catboost import CatBoostRegressor
from pathlib import Path
from IPython.display import display
import yfinance as yf
SEED = 42
np.random.seed(SEED)
torch.manual_seed(SEED)
torch.set_num_threads(2)
TAUS = (0.10, 0.50, 0.90)
NAMES = ("P10", "P50", "P90")
FEATURES = ["ret1", "vol20", "range_today", "upside_mean5", "volume_vs_20d"]
Path("figures").mkdir(exist_ok=True)
# Exclude the still-forming NY trading day, if any.
end_exclusive = pd.Timestamp.now(tz="America/New_York").strftime("%Y-%m-%d")
raw = yf.download(
"SPY", start="2015-01-01", end=end_exclusive,
auto_adjust=True, progress=False, threads=False,
)
if raw.empty:
raise RuntimeError("SPY download failed. Check internet access and rerun.")
if isinstance(raw.columns, pd.MultiIndex):
raw.columns = raw.columns.get_level_values(0)
raw.columns = [str(c).title() for c in raw.columns]
z = raw[["Open", "High", "Low", "Close", "Volume"]].sort_index().dropna()
z = z[~z.index.duplicated(keep="last")]
print(f"Downloaded {len(z):,} SPY trading days: {z.index[0].date()} to {z.index[-1].date()}")
display(z.tail(3))
2. Same features and target as Part 1
Every feature is available after today’s close. Tomorrow’s high and open appear only in the target, never in the inputs.
d = z.copy()
d["upside_today"] = 100 * (d["High"] / d["Open"] - 1)
d["target"] = d["upside_today"].shift(-1)
d["ret1"] = 100 * d["Close"].pct_change()
d["vol20"] = d["ret1"].rolling(20).std()
d["range_today"] = 100 * (d["High"] / d["Low"] - 1)
d["upside_mean5"] = d["upside_today"].rolling(5).mean()
d["volume_vs_20d"] = d["Volume"] / d["Volume"].rolling(20).mean() - 1
d = d.replace([np.inf, -np.inf], np.nan).dropna(subset=FEATURES + ["target"])
# Same 75/25 chronological split as the original CatBoost quickstart.
split = int(0.75 * len(d))
train = d.iloc[:split - 1].copy() # one-row gap: train label cannot use first test session
test = d.iloc[split:].copy()
X_train, y_train = train[FEATURES], train["target"]
X_test, y_test = test[FEATURES], test["target"]
assert train.index.max() < test.index.min()
print(f"Train: {len(train):,} rows | Test: {len(test):,} rows")
print(f"Test signal dates: {test.index.min().date()} to {test.index.max().date()}")Train: 2,219 rows | Test: 740 rows
Test signal dates: 2023-11-27 to 2026-11-063. PyTorch: three outputs + pinball loss
The neural network predicts P₁₀, P₅₀, and P₉₀ together. Pinball loss is the main conceptual change from standard neural network regression. Each quantile uses a different τ, and PyTorch handles the gradients through backpropagation.
Just as with CatBoost, we could also train three separate neural networks, one for each quantile, rather than a single network with three outputs. The joint approach allows the quantiles to share learned features, while separate models learn each quantile independently. Neither approach guarantees better performance, and both can suffer from quantile crossing.
Get Trade Mamba’s stories in your inbox
Join Medium for free to get updates from this writer.
Inputs are standardized using training data only. We use a chronological train/test split, without using information from the test period to fit the network.
scaler = StandardScaler()
x_tr = torch.tensor(scaler.fit_transform(X_train), dtype=torch.float32)
x_te = torch.tensor(scaler.transform(X_test), dtype=torch.float32)
y_tr = torch.tensor(y_train.to_numpy(dtype=np.float32)[:, None])
taus = torch.tensor(TAUS, dtype=torch.float32).view(1, -1)
class QuantileMLP(nn.Module):
def __init__(self, n_features):
super().__init__()
self.net = nn.Sequential(
nn.Linear(n_features, 48), nn.ReLU(),
nn.Linear(48, 48), nn.ReLU(),
nn.Linear(48, 3),
)
def forward(self, x):
return self.net(x)
def pinball_loss(pred, y):
error = y - pred # y: (N,1), pred: (N,3); broadcasting over quantiles
return torch.maximum(taus * error, (taus - 1) * error).mean()
model = QuantileMLP(n_features=len(FEATURES))
optimizer = torch.optim.Adam(model.parameters(), lr=0.003, weight_decay=1e-4)
loss_history = []
for epoch in range(500):
model.train()
optimizer.zero_grad()
loss = pinball_loss(model(x_tr), y_tr)
loss.backward()
optimizer.step()
loss_history.append(loss.item())
model.eval()
with torch.no_grad():
q_torch = model(x_te).cpu().numpy()
print(f"PyTorch final training mean pinball loss: {loss_history[-1]:.4f}")
l()PyTorch final training mean pinball loss: 0.09464. CatBoost: original comparison models
Use Part 1’s settings: three separately fitted Quantile models and one joint MultiQuantile model, all trained on the same rows as PyTorch.
q_separate = []
for tau in TAUS:
cb = CatBoostRegressor(
loss_function=f"Quantile:alpha={tau}",
iterations=300, depth=5, learning_rate=0.04,
random_seed=SEED, verbose=False, allow_writing_files=False,
)
cb.fit(X_train, y_train)
q_separate.append(cb.predict(X_test))
q_separate = np.column_stack(q_separate)
cb_joint = CatBoostRegressor(
loss_function="MultiQuantile:alpha=0.1,0.5,0.9",
iterations=300, depth=5, learning_rate=0.04,
random_seed=SEED, task_type="CPU",
verbose=False, allow_writing_files=False,
)
cb_joint.fit(X_train, y_train)
q_joint = np.asarray(cb_joint.predict(X_test)).reshape(-1, 3)q_separate = []
for tau in TAUS:
cb = CatBoostRegressor(
loss_function=f"Quantile:alpha={tau}",
iterations=300, depth=5, learning_rate=0.04,
random_seed=SEED, verbose=False, allow_writing_files=False,
)
cb.fit(X_train, y_train)
q_separate.append(cb.predict(X_test))
q_separate = np.column_stack(q_separate)
cb_joint = CatBoostRegressor(
loss_function="MultiQuantile:alpha=0.1,0.5,0.9",
iterations=300, depth=5, learning_rate=0.04,
random_seed=SEED, task_type="CPU",
verbose=False, allow_writing_files=False,
)
cb_joint.fit(X_train, y_train)
q_joint = np.asarray(cb_joint.predict(X_test)).reshape(-1, 3)5. Compare pinball loss and interval coverage
For calibrated quantiles, about 10% of observations should fall below P₁₀, 90% below P₉₀, and 80% within P₁₀–P₉₀ on average. We also report mean interval width and quantile-crossing frequency. Lower pinball loss is better; coverage closer to 80% is desirable but does not, by itself, imply useful intervals.
def evaluate(y_true, predictions):
y = np.asarray(y_true, dtype=float).reshape(-1, 1)
q = np.asarray(predictions, dtype=float)
t = np.array(TAUS)[None, :]
err = y - q
return {
"mean_pinball": np.maximum(t * err, (t - 1) * err).mean(),
"P10_actual_below_%": 100 * np.mean(y[:, 0] < q[:, 0]),
"P90_actual_below_%": 100 * np.mean(y[:, 0] < q[:, 2]),
"P10-P90_coverage_%": 100 * np.mean((y[:, 0] >= q[:, 0]) & (y[:, 0] <= q[:, 2])),
"interval_width_pct_points": np.mean(q[:, 2] - q[:, 0]),
"crossing_%": 100 * np.mean((q[:, 0] > q[:, 1]) | (q[:, 1] > q[:, 2])),
}
results = pd.DataFrame({
"PyTorch MLP": evaluate(y_test, q_torch),
"CatBoost separate": evaluate(y_test, q_separate),
"CatBoost MultiQuantile": evaluate(y_test, q_joint),
}).T.sort_values("mean_pinball")
display(results.round(3))
6. Visualize predictions on real SPY test dates
The shaded band is the predicted P₁₀–P₉₀ interval for next-session maximum upside. Each point is indexed by the date we made the prediction (at the close), not by the following session’s date.
view = predictions.tail(100)
fig, ax = plt.subplots(figsize=(12, 4.5))
ax.plot(view.index, view["actual_max_upside_%"], linewidth=1,
label="Actual next-session max upside")
ax.plot(view.index, view["PyTorch_P50"], linewidth=1.5, label="PyTorch P₅₀")
ax.fill_between(view.index,
view["PyTorch_P10"].to_numpy(),
view["PyTorch_P90"].to_numpy(),
alpha=0.23, label="PyTorch P₁₀–P₉₀ (nominal 80%)")
ax.set(title="SPY: PyTorch quantile forecasts on held-out data",
ylabel="Next-session max upside (%)", xlabel="Signal date")
ax.legend()
fig.tight_layout()
fig.savefig("figures/spy_pytorch_prediction_intervals.png", dpi=160)
plt.show()
Takeaways
The results are quite close across all three models.
CatBoost MultiQuantile achieved the lowest mean pinball loss (0.095), followed by PyTorch (0.096) and separate CatBoost models (0.097).
PyTorch came closest to the nominal 80% prediction interval coverage, achieving 78.8%, compared with 76.6% for CatBoost MultiQuantile. However, its intervals were also slightly wider (0.897 vs. 0.844 percentage points).
None of the models was perfectly calibrated. For example, approximately 12.6% of actual outcomes fell below PyTorch’s P₁₀ predictions, rather than the expected 10%.
Quantile crossing was rare: 0.135% for PyTorch, 0.405% for separate CatBoost models, and none for CatBoost MultiQuantile.
The main takeaway is that a simple PyTorch MLP performed competitively with CatBoost, but didn’t clearly outperform it. Jointly predicting multiple quantiles worked well with both approaches, and CatBoost remained slightly better on pinball loss.
Of course, these results come from one dataset and test period, so we shouldn’t draw broader conclusions about which approach is better. But in my experience, CatBoost is king when it comes to tabular data.
