Time Series Forecasting in 2026: From ARIMA to Transformers and Foundation Models
Time series forecasting — predicting future values of a sequence based on historical patterns — is one of the most consequential applications of data science. Accurate forecasts power inventory management at Amazon, grid stability at energy utilities, cash flow planning at banks, and demand prediction at logistics companies. A 10% improvement in forecast accuracy for a major retailer can translate to hundreds of millions of dollars in reduced inventory costs and improved service levels.
The field has transformed dramatically: classical statistical methods (ARIMA, Exponential Smoothing) that dominated for decades have been supplemented and often surpassed by machine learning approaches (gradient boosting, LSTMs) and, most recently, transformer-based architectures and foundation models pretrained on millions of time series. This guide covers the full spectrum of modern time series forecasting: classical methods, ML approaches, deep learning architectures, and the emerging class of time series foundation models.
Classical Statistical Methods
ARIMA and SARIMA
ARIMA (AutoRegressive Integrated Moving Average) models time series as a function of past values (AR), past errors (MA), and their differences (I for integration, to achieve stationarity). SARIMA extends ARIMA with seasonal components. These models are interpretable, have well-understood statistical properties, and work well for univariate series with stable patterns.
import pandas as pd
import numpy as np
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.stattools import adfuller
from pmdarima import auto_arima
import matplotlib.pyplot as plt
# Load and prepare data
df = pd.read_csv('monthly_sales.csv', parse_dates=['date'], index_col='date')
df = df.asfreq('MS') # Monthly frequency
# Test for stationarity
adf_result = adfuller(df['sales'].dropna())
print(f'ADF Statistic: {adf_result[0]:.4f}')
print(f'p-value: {adf_result[1]:.4f}') # p < 0.05 = stationary
# Auto-select SARIMA parameters
model = auto_arima(
df['sales'],
seasonal=True,
m=12, # Monthly seasonality period
information_criterion='aic',
trace=True,
error_action='ignore',
suppress_warnings=True,
stepwise=True, # Faster search
)
print(model.summary())
# Forecast 12 months ahead
forecast, conf_int = model.predict(n_periods=12, return_conf_int=True)
future_dates = pd.date_range(df.index[-1] + pd.DateOffset(months=1), periods=12, freq='MS')
plt.figure(figsize=(12, 5))
plt.plot(df.index, df['sales'], label='Historical')
plt.plot(future_dates, forecast, label='Forecast', color='red')
plt.fill_between(future_dates, conf_int[:, 0], conf_int[:, 1], alpha=0.3, color='red')
plt.legend()
plt.title('SARIMA Forecast')
plt.show()Prophet
Facebook's Prophet models time series as the sum of trend, seasonality, and holiday effects. It handles missing data, outliers, and multiple seasonality periods (daily, weekly, yearly) with minimal tuning. Prophet is particularly effective for business time series with strong seasonal patterns and known holiday effects.
from prophet import Prophet
import pandas as pd
# Prophet requires columns 'ds' (date) and 'y' (value)
df_prophet = df.reset_index().rename(columns={'date': 'ds', 'sales': 'y'})
# Configure model
m = Prophet(
changepoint_prior_scale=0.1, # Flexibility of trend changepoints
seasonality_prior_scale=10, # Strength of seasonality
holidays_prior_scale=10,
yearly_seasonality=True,
weekly_seasonality=False, # Monthly data, no weekly pattern
daily_seasonality=False,
)
# Add custom seasonality
m.add_seasonality(name='quarterly', period=91.25, fourier_order=8)
# Add country holidays
m.add_country_holidays(country_name='US')
m.fit(df_prophet)
# Create future dataframe
future = m.make_future_dataframe(periods=12, freq='MS')
forecast = m.predict(future)
# Plot components
fig = m.plot_components(forecast)Machine Learning for Time Series
Feature Engineering for Tabular ML
Gradient boosting models (LightGBM, XGBoost) can achieve state-of-the-art time series forecasting when trained on carefully engineered features. The key is transforming the forecasting problem into a supervised learning problem: for each point in time t, create a row with lag features (values at t-1, t-2, ..., t-n), rolling window features (mean, std, min, max over past windows), and calendar features (day of week, month, holiday indicators).
import pandas as pd
import numpy as np
import lightgbm as lgb
from sklearn.model_selection import TimeSeriesSplit
def create_time_features(df, target_col, lags=[1,2,3,6,12,24],
windows=[7, 14, 30, 90]):
df = df.copy()
# Lag features
for lag in lags:
df[f'lag_{lag}'] = df[target_col].shift(lag)
# Rolling statistics
for window in windows:
df[f'rolling_mean_{window}'] = df[target_col].shift(1).rolling(window).mean()
df[f'rolling_std_{window}'] = df[target_col].shift(1).rolling(window).std()
df[f'rolling_min_{window}'] = df[target_col].shift(1).rolling(window).min()
df[f'rolling_max_{window}'] = df[target_col].shift(1).rolling(window).max()
# Calendar features
df['month'] = df.index.month
df['quarter'] = df.index.quarter
df['year'] = df.index.year
df['day_of_week'] = df.index.dayofweek
df['is_month_end'] = df.index.is_month_end.astype(int)
# Trend feature
df['time_idx'] = np.arange(len(df))
return df.dropna()
df_features = create_time_features(df, 'sales')
# Time series cross-validation
tscv = TimeSeriesSplit(n_splits=5, test_size=12) # 12 months test
for fold, (train_idx, val_idx) in enumerate(tscv.split(df_features)):
X_train = df_features.iloc[train_idx].drop(columns=['sales'])
y_train = df_features.iloc[train_idx]['sales']
X_val = df_features.iloc[val_idx].drop(columns=['sales'])
y_val = df_features.iloc[val_idx]['sales']
model = lgb.LGBMRegressor(
n_estimators=1000, learning_rate=0.05,
num_leaves=64, min_child_samples=20,
early_stopping_rounds=50, verbose=-1
)
model.fit(X_train, y_train,
eval_set=[(X_val, y_val)])Deep Learning for Time Series
N-BEATS and N-HiTS
N-BEATS (Neural Basis Expansion Analysis for Time Series) and its successor N-HiTS are interpretable neural forecasting models that decompose the forecast into trend and seasonality components using basis function expansion. They achieved state-of-the-art results on M4 competition benchmarks and are effective for univariate time series.
from neuralforecast import NeuralForecast
from neuralforecast.models import NHITS, NBEATS, PatchTST
from neuralforecast.losses.pytorch import MAE, SMAPE
import pandas as pd
# Prepare data in NeuralForecast format
# Requires columns: unique_id, ds, y
nf = NeuralForecast(
models=[
NHITS(
h=12, # forecast horizon (12 months)
input_size=24, # lookback window
loss=MAE(),
scaler_type='standard',
n_freq_downsample=[2, 1, 1],
max_steps=1000,
early_stop_patience_steps=10,
),
NBEATS(
h=12,
input_size=24,
stack_types=['trend', 'seasonality'],
n_blocks=[3, 3],
loss=MAE(),
max_steps=1000,
),
],
freq='MS', # Monthly frequency
)
nf.fit(df_panel) # df_panel has unique_id, ds, y columns
forecast = nf.predict() # Returns forecasts for all seriesTemporal Fusion Transformer (TFT)
TFT is a powerful architecture for multi-horizon forecasting with multiple covariates. It incorporates: Variable Selection Networks that learn which input variables are relevant for each step; Gated Residual Networks for non-linear feature processing; LSTM encoders for processing sequential information; Multi-head attention for capturing long-range dependencies; and quantile outputs for probabilistic forecasting.
from pytorch_forecasting import TemporalFusionTransformer, TimeSeriesDataSet
from pytorch_forecasting.metrics import QuantileLoss
from lightning.pytorch import Trainer
# Prepare dataset
max_encoder_length = 36
max_prediction_length = 12
training_cutoff = df['time_idx'].max() - max_prediction_length
training = TimeSeriesDataSet(
df[df.time_idx <= training_cutoff],
time_idx='time_idx',
target='sales',
group_ids=['store_id', 'product_id'],
min_encoder_length=max_encoder_length // 2,
max_encoder_length=max_encoder_length,
min_prediction_length=1,
max_prediction_length=max_prediction_length,
static_categoricals=['store_id', 'product_id', 'category'],
static_reals=['avg_store_size'],
time_varying_known_categoricals=['holiday'],
time_varying_known_reals=['price', 'time_idx'],
time_varying_unknown_reals=['sales', 'volume'],
target_normalizer=GroupNormalizer(groups=['store_id', 'product_id']),
)
train_dataloader = training.to_dataloader(train=True, batch_size=128)
tft = TemporalFusionTransformer.from_dataset(
training,
learning_rate=0.03,
hidden_size=64,
attention_head_size=4,
dropout=0.1,
loss=QuantileLoss(),
)
trainer = Trainer(max_epochs=30, gradient_clip_val=0.1, accelerator='gpu')
trainer.fit(tft, train_dataloaders=train_dataloader)Time Series Foundation Models
The most significant development in time series forecasting in 2024-2026 is the emergence of foundation models pretrained on massive collections of time series data. Analogous to how LLMs are pretrained on text and then adapted for downstream tasks, time series foundation models are pretrained on millions of diverse time series and can forecast new series zero-shot (without any fine-tuning) or with minimal examples.
TimesFM (Google DeepMind)
TimesFM is pretrained on 100 billion time points from Google Trends, Wikipedia page views, and synthetic data. It uses a decoder-only transformer architecture similar to GPT, treating time series as sequences of numerical patches. TimesFM achieves zero-shot accuracy comparable to supervised models fine-tuned on the target dataset.
Moirai (Salesforce Research)
Moirai is an open-source time series foundation model trained on LOTSA data — a 27-billion time point dataset spanning diverse domains (web traffic, energy consumption, retail sales, financial data). Moirai uses a flexible input tokenization that handles multiple frequencies without frequency-specific fine-tuning.
from uni2ts.model.moirai import MoiraiModule, MoiraiForecast
from uni2ts.model.moirai import MoiraiPipeline
# Zero-shot forecasting with Moirai
pipeline = MoiraiPipeline.from_pretrained(
"Salesforce/moirai-1.0-R-large",
prediction_length=12,
context_length=200,
patch_size="auto",
num_samples=100, # For probabilistic forecasting
target_dim=1,
feat_dynamic_real_dim=0,
)
forecast = pipeline(test_data, past_dynamic_reals=None)
print(f"Forecast shape: {forecast.samples.shape}") # (100 samples, 12 timesteps)Chronos (Amazon)
Amazon's Chronos tokenizes time series values into discrete bins and trains a T5 transformer to predict future token sequences. Uniquely, Chronos is trained primarily on synthetic data generated from Gaussian processes and other stochastic processes, augmented with real-world data. This synthetic pretraining reduces the need for large real-world datasets while achieving strong zero-shot accuracy.
from chronos import ChronosPipeline
import torch
pipeline = ChronosPipeline.from_pretrained(
"amazon/chronos-t5-large",
device_map="cuda",
torch_dtype=torch.bfloat16,
)
# Forecast using historical data
context = torch.tensor(monthly_sales_history)
forecast = pipeline.predict(
context=context.unsqueeze(0), # Add batch dimension
prediction_length=12,
num_samples=20,
)
# Compute quantiles for prediction intervals
low, median, high = np.quantile(forecast[0].numpy(), [0.1, 0.5, 0.9], axis=0)
print(f"Median forecast: {median}")
print(f"90% interval: [{low}, {high}]")Evaluation and Practical Considerations
Metrics for Forecasting Accuracy
MAE (Mean Absolute Error): Average absolute difference between forecast and actual. Easy to interpret; same units as the target. Robust to outliers compared to MSE.
RMSE (Root Mean Squared Error): Square root of mean squared error. Penalizes large errors more than MAE; more sensitive to outliers.
MAPE (Mean Absolute Percentage Error): Mean absolute percentage deviation. Scale-independent, enabling comparison across series with different magnitudes. Problematic when actual values are near zero (percentage error explodes).
SMAPE (Symmetric MAPE): Addresses MAPE's asymmetry; bounded between 0% and 200%. Standard in the M competitions.
MASE (Mean Absolute Scaled Error): Scales MAE by the MAE of a naive seasonal model. MASE > 1 means the model is worse than the seasonal naive baseline; the target for any reasonable model is MASE < 1.
Backtesting with Expanding and Rolling Windows
import numpy as np
def time_series_cv(series, model_fn, h=12, n_windows=5, step_size=6):
"""Expanding window cross-validation for time series"""
n = len(series)
min_train_size = n - (n_windows - 1) * step_size - h
results = []
for i in range(n_windows):
train_end = min_train_size + i * step_size
val_end = train_end + h
if val_end > n:
break
train = series[:train_end]
actual = series[train_end:val_end]
model = model_fn()
model.fit(train)
forecast = model.predict(h)
mae = np.mean(np.abs(actual - forecast))
smape = np.mean(2 * np.abs(actual - forecast) / (np.abs(actual) + np.abs(forecast)))
results.append({'window': i, 'mae': mae, 'smape': smape})
return pd.DataFrame(results)
# Evaluate SARIMA
def sarima_model():
return SARIMAX(order=(1,1,1), seasonal_order=(1,1,1,12))
cv_results = time_series_cv(df['sales'].values, sarima_model, h=12, n_windows=5)
print(cv_results.mean())Choosing the Right Model
Classical (ARIMA, Prophet): When you have a single or small number of series, need interpretable models with confidence intervals, have sufficient history (2+ full seasonal cycles), and want to incorporate domain knowledge (holidays, events). Fast to train and evaluate.
ML (LightGBM, XGBoost): When you have many related series (retail demand planning, store-level forecasting), rich exogenous features (pricing, promotions, weather), and tabular data that is well-represented. Often top performers in Kaggle competitions.
Deep Learning (TFT, N-HiTS): When you have many related series and want to model cross-series relationships; when you have many covariates; when you need probabilistic forecasts; when you have enough data to train neural networks (thousands of time steps per series).
Foundation Models (Chronos, Moirai, TimesFM): When you have few historical observations (cold start); when you need zero-shot forecasting across diverse series without per-series training; when you want a strong baseline with minimal effort. As these models mature, they are becoming the default starting point for new forecasting tasks.
The state of the art in 2026 is often a hybrid: use a foundation model for zero-shot baseline, fine-tune on domain-specific data if available, and blend with domain-specific statistical models for special events (promotions, holidays) that may not be captured by the historical patterns the model was trained on.
Comments
Post a Comment