Time Series Analysis with Python – statsmodels, Prophet & LSTM 2026

Time series data is everywhere — stock prices, sales figures, website traffic, sensor readings, energy consumption. Unlike cross-sectional data, time series observations are ordered and dependent on past values. This guide covers the full time series analysis workflow: decomposition, stationarity, ARIMA, Prophet, and LSTM, so you can confidently tackle any forecasting problem.

Time Series Decomposition

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose, STL

# Generate sample sales data with trend + seasonality + noise
np.random.seed(42)
dates  = pd.date_range('2020-01-01', '2025-12-31', freq='ME')
trend  = np.linspace(100, 200, len(dates))
season = 20 * np.sin(2 * np.pi * np.arange(len(dates)) / 12)
noise  = np.random.normal(0, 5, len(dates))

ts = pd.Series(trend + season + noise, index=dates, name='sales')

# Classical decomposition
decomp = seasonal_decompose(ts, model='additive', period=12)

fig, axes = plt.subplots(4, 1, figsize=(12, 10))
decomp.observed.plot(ax=axes[0], title='Observed')
decomp.trend.plot(ax=axes[1], title='Trend')
decomp.seasonal.plot(ax=axes[2], title='Seasonality')
decomp.resid.plot(ax=axes[3], title='Residuals')
plt.tight_layout(); plt.show()

# STL decomposition (more robust, handles outliers better)
stl = STL(ts, period=12, robust=True)
res = stl.fit()
res.plot(); plt.show()

print(f'Seasonal strength: {max(0, 1 - res.resid.var() / (res.seasonal + res.resid).var()):.3f}')
print(f'Trend strength:    {max(0, 1 - res.resid.var() / (res.trend + res.resid).var()):.3f}')

Stationarity Tests

from statsmodels.tsa.stattools import adfuller, kpss

def check_stationarity(series, name='Series'):
    print(f'=== Stationarity Check: {name} ===')

    # ADF test: H0=non-stationary, reject H0 → stationary
    adf_result = adfuller(series.dropna(), autolag='AIC')
    print(f'ADF p-value: {adf_result[1]:.4f} '
          f'({"✅ stationary" if adf_result[1] < 0.05 else "❌ non-stationary"})')

    # KPSS test: H0=stationary, reject H0 → non-stationary
    kpss_result = kpss(series.dropna(), regression='c', nlags='auto')
    print(f'KPSS p-value: {kpss_result[1]:.4f} '
          f'({"✅ stationary" if kpss_result[1] > 0.05 else "❌ non-stationary"})')
    print()

check_stationarity(ts, 'Raw Sales')

# Make stationary by differencing
ts_diff = ts.diff().dropna()
check_stationarity(ts_diff, 'First Difference')

ARIMA / SARIMA with statsmodels

from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from sklearn.metrics import mean_absolute_error, mean_squared_error

# ACF/PACF to identify p, d, q
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 6))
plot_acf(ts_diff, lags=24, ax=ax1)
plot_pacf(ts_diff, lags=24, ax=ax2)
plt.tight_layout(); plt.show()

# Train/test split (last 12 months = test)
train = ts[:-12]
test  = ts[-12:]

# SARIMA(p,d,q)(P,D,Q,s) — handles seasonal data
model = SARIMAX(train,
                order=(1, 1, 1),
                seasonal_order=(1, 1, 1, 12),
                enforce_stationarity=False,
                enforce_invertibility=False)

result = model.fit(disp=False)
print(result.summary())

# Forecast
forecast = result.get_forecast(steps=12)
mean_fc  = forecast.predicted_mean
ci       = forecast.conf_int()

plt.figure(figsize=(12, 5))
train.plot(label='Train'); test.plot(label='Actual')
mean_fc.plot(label='Forecast', color='red')
plt.fill_between(ci.index, ci.iloc[:, 0], ci.iloc[:, 1],
                 alpha=0.2, color='red', label='95% CI')
plt.legend(); plt.title('SARIMA Forecast'); plt.show()

mae  = mean_absolute_error(test, mean_fc)
rmse = np.sqrt(mean_squared_error(test, mean_fc))
mape = (abs(test - mean_fc) / abs(test)).mean() * 100
print(f'MAE: {mae:.2f} | RMSE: {rmse:.2f} | MAPE: {mape:.1f}%')

Facebook Prophet

pip install prophet

from prophet import Prophet
from prophet.plot import plot_plotly, plot_components_plotly

# Prophet requires 'ds' (date) and 'y' (value) columns
df_prophet = ts.reset_index()
df_prophet.columns = ['ds', 'y']

train_p = df_prophet[:-12]
test_p  = df_prophet[-12:]

m = Prophet(
    yearly_seasonality=True,
    weekly_seasonality=False,
    daily_seasonality=False,
    seasonality_mode='additive',   # or 'multiplicative'
    changepoint_prior_scale=0.05,  # trend flexibility (0.001–0.5)
    seasonality_prior_scale=10,    # seasonality strength
    interval_width=0.95
)

# Add custom seasonality
m.add_seasonality(name='quarterly', period=91.25, fourier_order=5)

# Add regressors (external variables)
# m.add_regressor('holiday_flag')

m.fit(train_p)

# Make future dataframe
future   = m.make_future_dataframe(periods=12, freq='ME')
forecast = m.predict(future)

# Plot
fig1 = m.plot(forecast); plt.title('Prophet Forecast')
fig2 = m.plot_components(forecast)
plt.show()

# Evaluate
pred  = forecast.tail(12)['yhat'].values
actua = test_p['y'].values
mape  = np.mean(np.abs((actua - pred) / actua)) * 100
print(f'Prophet MAPE: {mape:.1f}%')

LSTM for Time Series

import torch
import torch.nn as nn
from sklearn.preprocessing import MinMaxScaler

# Prepare sequences
scaler     = MinMaxScaler()
ts_scaled  = scaler.fit_transform(ts.values.reshape(-1, 1))

SEQ_LEN = 12

def make_sequences(data, seq_len):
    X, y = [], []
    for i in range(len(data) - seq_len):
        X.append(data[i:i + seq_len])
        y.append(data[i + seq_len])
    return np.array(X), np.array(y)

X, y   = make_sequences(ts_scaled, SEQ_LEN)
split  = len(X) - 12
X_tr, X_te = X[:split], X[split:]
y_tr, y_te = y[:split], y[split:]

X_tr = torch.FloatTensor(X_tr)
X_te = torch.FloatTensor(X_te)
y_tr = torch.FloatTensor(y_tr)
y_te = torch.FloatTensor(y_te)

class LSTMForecaster(nn.Module):
    def __init__(self, input_size=1, hidden=64, layers=2, dropout=0.2):
        super().__init__()
        self.lstm = nn.LSTM(input_size, hidden, layers,
                             batch_first=True, dropout=dropout)
        self.fc   = nn.Linear(hidden, 1)

    def forward(self, x):
        out, _ = self.lstm(x)
        return self.fc(out[:, -1, :])

model_lstm = LSTMForecaster()
optimiser  = torch.optim.Adam(model_lstm.parameters(), lr=1e-3)
criterion  = nn.MSELoss()

for epoch in range(200):
    model_lstm.train()
    pred  = model_lstm(X_tr)
    loss  = criterion(pred, y_tr)
    optimiser.zero_grad(); loss.backward(); optimiser.step()
    if (epoch + 1) % 50 == 0:
        print(f'Epoch {epoch+1}: Loss={loss.item():.6f}')

model_lstm.eval()
with torch.no_grad():
    preds = model_lstm(X_te).numpy()
preds = scaler.inverse_transform(preds)
actua = scaler.inverse_transform(y_te.numpy())
mape  = np.mean(np.abs((actua - preds) / actua)) * 100
print(f'LSTM MAPE: {mape:.1f}%')

Choosing the Right Model

Use ARIMA/SARIMA for short, stationary series with clear seasonality and when interpretability matters — it is fast and gives confidence intervals natively. Use Prophet when your series has multiple seasonalities, holidays, or structural change points you want to specify — it handles missing data and outliers gracefully. Use LSTM when you have 2+ years of high-frequency data (daily or finer), multiple related series, or patterns that classical models cannot capture. In practice, always baseline with SARIMA before investing in LSTM — the complexity is only worth it if you have enough data and the MAPE improvement justifies it.

Conclusion

Time series forecasting is part statistics, part domain knowledge. The best model is rarely the most complex one. Always start with decomposition to understand what you are dealing with — trend, seasonality, cycles, noise. Verify stationarity before fitting ARIMA. Use Prophet as your fast, production-ready workhorse for business forecasting. Reserve LSTM for high-data, high-complexity scenarios. And always validate on a holdout set that is temporally separated from training — cross-validation must respect the time ordering of your data.