Time Series ML: Forecasting, Anomaly Detection, and Feature Engineering

Time series data is everywhere: server metrics, sales figures, sensor readings, financial prices, user activity counts. Standard ML models applied to time series without modification often fail because they ignore temporal structure. This guide covers the practical ML engineering for time series work.

Why Time Series Is Different

Standard ML assumptions that break:

1. Independent and identically distributed (i.i.d.): Samples are not independent - yesterday's value predicts today's.
2. Random train/test split: Splitting randomly allows future data into training (data leakage). You must always train on the past, test on the future.
3. Stationarity: Many time series have trends and seasonality that shift the distribution over time.

Feature Engineering for Time Series

This is where most of the value comes from. Raw time series values are rarely the best features:

pythonCopy
import pandas as pd
import numpy as np

def create_time_features(df: pd.DataFrame, timestamp_col: str, value_col: str) -> pd.DataFrame:
    df = df.copy()
    df[timestamp_col] = pd.to_datetime(df[timestamp_col])
    df = df.sort_values(timestamp_col)

    # Calendar features (encode seasonality)
    df['hour'] = df[timestamp_col].dt.hour
    df['day_of_week'] = df[timestamp_col].dt.dayofweek
    df['day_of_month'] = df[timestamp_col].dt.day
    df['month'] = df[timestamp_col].dt.month
    df['is_weekend'] = (df['day_of_week'] >= 5).astype(int)
    df['quarter'] = df[timestamp_col].dt.quarter

    # Lag features (past values as features)
    for lag in [1, 2, 3, 7, 14, 28]:
        df[f'lag_{lag}'] = df[value_col].shift(lag)

    # Rolling statistics (local trend and variance)
    for window in [7, 14, 30]:
        df[f'rolling_mean_{window}'] = df[value_col].shift(1).rolling(window).mean()
        df[f'rolling_std_{window}'] = df[value_col].shift(1).rolling(window).std()
        df[f'rolling_min_{window}'] = df[value_col].shift(1).rolling(window).min()
        df[f'rolling_max_{window}'] = df[value_col].shift(1).rolling(window).max()

    # Difference features (captures rate of change)
    df['diff_1'] = df[value_col].diff(1)
    df['diff_7'] = df[value_col].diff(7)   # week-over-week change
    df['pct_change_7'] = df[value_col].pct_change(7)

    # Expanding features (cumulative statistics  -  careful of leakage)
    df['expanding_mean'] = df[value_col].shift(1).expanding().mean()

    return df.dropna()

Critical: Always shift lag features by at least 1 (shift(1)) to ensure you are not including the current value in its own feature computation.

Temporal Cross-Validation

Never use k-fold on time series. Use time series CV with a walk-forward approach:

pythonCopy
from sklearn.model_selection import TimeSeriesSplit
import numpy as np

tscv = TimeSeriesSplit(n_splits=5, gap=0, test_size=30)  # 30 days per fold

scores = []
for train_idx, val_idx in tscv.split(X):
    X_train, X_val = X[train_idx], X[val_idx]
    y_train, y_val = y[train_idx], y[val_idx]

    model.fit(X_train, y_train)
    preds = model.predict(X_val)
    mae = np.mean(np.abs(preds - y_val))
    scores.append(mae)

print(f"MAE per fold: {scores}")
print(f"Mean MAE: {np.mean(scores):.2f} ± {np.std(scores):.2f}")

Forecasting Models

Gradient Boosted Trees (best for tabular time series):

pythonCopy
import lightgbm as lgb
import pandas as pd

# GBTs with the lag/rolling features above work extremely well
# They have won most major time series forecasting competitions (M5, etc.)

model = lgb.LGBMRegressor(
    n_estimators=1000,
    learning_rate=0.05,
    num_leaves=63,
    feature_fraction=0.8,
    bagging_fraction=0.8,
    bagging_freq=5,
    min_child_samples=20
)
model.fit(
    X_train, y_train,
    eval_set=[(X_val, y_val)],
    callbacks=[lgb.early_stopping(50), lgb.log_evaluation(100)]
)

Prophet (for interpretable trend + seasonality decomposition):

pythonCopy
from prophet import Prophet

# Prophet expects 'ds' (datestamp) and 'y' (value)
df_prophet = df.rename(columns={'date': 'ds', 'sales': 'y'})

model = Prophet(
    changepoint_prior_scale=0.1,  # flexibility of trend changes
    seasonality_prior_scale=10,
    yearly_seasonality=True,
    weekly_seasonality=True,
    daily_seasonality=False  # for daily data
)

# Add custom seasonality if needed (e.g., monthly patterns)
model.add_seasonality(name='monthly', period=30.5, fourier_order=5)

model.fit(df_prophet)
future = model.make_future_dataframe(periods=30)  # forecast 30 days
forecast = model.predict(future)

# Forecast components: trend + weekly seasonality + yearly seasonality
model.plot_components(forecast)

Anomaly Detection

Statistical approach (Z-score on residuals):

pythonCopy
def detect_anomalies_zscore(series: pd.Series, window: int = 30,
                              threshold: float = 3.0) -> pd.Series:
    """
    Flag points that deviate from the rolling mean by more than threshold standard deviations.
    """
    rolling_mean = series.rolling(window=window, center=True).mean()
    rolling_std = series.rolling(window=window, center=True).std()

    z_scores = (series - rolling_mean) / (rolling_std + 1e-8)
    return z_scores.abs() > threshold

anomalies = detect_anomalies_zscore(df['request_count'], window=24, threshold=3.0)
print(f"Anomalies detected: {anomalies.sum()}")

Isolation Forest (handles multivariate anomalies):

pythonCopy
from sklearn.ensemble import IsolationForest
import numpy as np

# Use multiple correlated features for better anomaly detection
X_multi = df[['cpu_usage', 'memory_usage', 'request_count', 'error_rate']].values

iso_forest = IsolationForest(
    contamination=0.05,  # expected fraction of anomalies
    n_estimators=200,
    random_state=42
)
anomaly_labels = iso_forest.fit_predict(X_multi)
# -1 = anomaly, 1 = normal

anomaly_scores = iso_forest.score_samples(X_multi)  # lower = more anomalous
df['anomaly_score'] = anomaly_scores
df['is_anomaly'] = anomaly_labels == -1

Evaluation Metrics for Forecasting

Metric	Formula	When to use
MAE	mean(	y - ŷ
RMSE	sqrt(mean((y - ŷ)²))	When large errors matter more
MAPE	mean(	y - ŷ
SMAPE	mean(2	y - ŷ
MASE	MAE / MAE(naive)	Scales to the naive seasonal forecast

MASE (Mean Absolute Scaled Error) is the most robust for comparing across time series of different scales - it compares your model to a naive seasonal baseline (predict what you observed one season ago).

Handling Non-Stationarity

Most real-world time series have trends and seasonality. Test and remove them:

pythonCopy
from statsmodels.tsa.stattools import adfuller

def test_stationarity(series: pd.Series) -> dict:
    """Augmented Dickey-Fuller test. Null hypothesis: series has unit root (non-stationary)."""
    result = adfuller(series.dropna())
    return {
        "adf_statistic": result[0],
        "p_value": result[1],
        "is_stationary": result[1] < 0.05
    }

# If not stationary, difference the series
df['sales_diff'] = df['sales'].diff(1)
df['sales_diff_seasonal'] = df['sales'].diff(7)  # remove weekly seasonality

# After modeling on differenced series, un-difference predictions
last_actual = df['sales'].iloc[-1]
predictions_original = last_actual + np.cumsum(predictions_diff)

Common Mistakes

Using random train/test splits on time-series data. A random split allows future data to leak into the training set because random samples from later time periods end up in training. This produces severely optimistic evaluation metrics that collapse in production where the model always predicts the future from the past. Always split chronologically: train on data before a cutoff, test on data after it.

Treating seasonality as noise. Seasonality is a deterministic, predictable signal and is often the dominant component of a time-series. Models that do not explicitly represent seasonal patterns (e.g., a plain ARIMA without seasonal terms) will leave systematic error on the table that a simple seasonal naive baseline would have captured. Decompose your series before choosing a model family.

Not checking for autocorrelation in residuals. After fitting any time-series model, plot the ACF (autocorrelation function) of residuals. If the residuals are autocorrelated, the model has left predictable signal in the errors, which means a better model exists. White-noise residuals are the minimum bar for a time-series model that has extracted the available signal.

What to Practice Next

• Take any time-series dataset, split it chronologically at the 80% mark, train a last-value (naive) baseline and a simple model (ARIMA or linear trend), and compute MAE for both on the test set.
• Decompose a seasonal time series into trend, seasonal, and residual components using statsmodels.tsa.seasonal_decompose; confirm the residual component looks like white noise.
• Plot the ACF and PACF of a raw time series and its first difference; use them to identify candidate ARIMA orders before fitting.