Time Series Analysis with Statsmodels, Prophet and LSTM – Complete Guide

Time series analysis is the practice of extracting patterns, structure, and forecasts from data indexed by time — and it requires a fundamentally different toolkit from cross-sectional data analysis. The ordering of observations matters, stationarity assumptions replace i.i.d. assumptions, and models must respect causality (you cannot use future data to predict the past). This guide covers the complete time series analysis workflow: decomposition, stationarity testing, classical statistical models (ARIMA, SARIMA with statsmodels), modern ML-based forecasting (Prophet), and deep learning approaches (LSTM). Each method is illustrated with working code and guidance on when to choose each.

This guide is the analytical companion to our Time Series Forecasting Interview Q&A (which covers the theory and interview questions). The feature engineering for time series — lag features, rolling statistics, calendar features — is covered in our Feature Engineering guide. The LSTM architecture used for sequence modelling is explained in our Neural Network Architectures guide. Visualisation of time series data uses the techniques in our Matplotlib and Seaborn guide.

Time Series Decomposition and Stationarity

Most time series can be decomposed into components: Trend (long-term direction), Seasonality (regular periodic fluctuations — daily, weekly, annual), Cyclical (irregular multi-year fluctuations — economic cycles), and Residual (random noise). Additive decomposition: Y_t = Trend_t + Seasonal_t + Residual_t (appropriate when seasonal amplitude is constant). Multiplicative decomposition: Y_t = Trend_t * Seasonal_t * Residual_t (appropriate when seasonal amplitude grows with the level — common in economic and retail data).

Stationarity — the property that mean, variance, and autocorrelation structure are constant over time — is assumed by most classical time series models (ARIMA, GARCH). A non-stationary series must be made stationary through differencing (subtract consecutive observations) before fitting these models. The Augmented Dickey-Fuller (ADF) test formally tests for a unit root (non-stationarity).

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose, STL
from statsmodels.tsa.stattools import adfuller, kpss
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# Load example time series (monthly retail sales)
df = pd.read_csv('retail_sales.csv', index_col='date', parse_dates=True)
ts = df['sales']

# --- Decomposition ---
# Classical additive decomposition
decomp = seasonal_decompose(ts, model='additive', period=12)
fig, axes = plt.subplots(4, 1, figsize=(12, 10), sharex=True)
for ax, (label, component) in zip(axes, [
    ('Observed',   decomp.observed),
    ('Trend',      decomp.trend),
    ('Seasonal',   decomp.seasonal),
    ('Residual',   decomp.resid),
]):
    ax.plot(component); ax.set_ylabel(label)
plt.suptitle('Time Series Decomposition'); plt.tight_layout(); plt.show()

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

# --- Stationarity testing ---
def test_stationarity(series, name='series'):
    # ADF test: H0 = unit root present (non-stationary)
    adf_stat, adf_p, _, _, crit_vals, _ = adfuller(series.dropna())
    # KPSS test: H0 = series is stationary (reverse of ADF)
    kpss_stat, kpss_p, _, kpss_crit = kpss(series.dropna(), regression='c')

    print('Stationarity tests for: ' + name)
    print('  ADF  stat =', round(adf_stat, 4), '  p =', round(adf_p, 4),
          '  Stationary?', 'Yes' if adf_p < 0.05 else 'No')
    print('  KPSS stat =', round(kpss_stat, 4), '  p =', round(kpss_p, 4),
          '  Stationary?', 'Yes' if kpss_p > 0.05 else 'No')

test_stationarity(ts, 'Raw sales')
ts_diff = ts.diff().dropna()
test_stationarity(ts_diff, 'First-differenced sales')

# --- ACF and PACF plots to identify ARIMA orders ---
fig, axes = plt.subplots(1, 2, figsize=(14, 4))
plot_acf(ts_diff, lags=40, ax=axes[0],  title='ACF (identify MA order q)')
plot_pacf(ts_diff, lags=40, ax=axes[1], title='PACF (identify AR order p)')
plt.tight_layout(); plt.show()

ARIMA and SARIMA with Statsmodels

ARIMA(p, d, q) combines three components: AR(p) — autoregressive: the series depends on its own past p values; I(d) — integrated: d rounds of differencing to achieve stationarity; MA(q) — moving average: the series depends on past q forecast errors. SARIMA(p,d,q)(P,D,Q)[s] adds seasonal counterparts for the seasonal period s. The Box-Jenkins methodology for selecting orders: (1) Make stationary via differencing (determine d); (2) Inspect ACF for MA order q; (3) Inspect PACF for AR order p; (4) Fit, check residuals for white noise (Ljung-Box test).

from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.stats.diagnostic import acorr_ljungbox
import pmdarima as pm    # pip install pmdarima

# --- Auto-ARIMA: automatically selects best (p,d,q)(P,D,Q) ---
auto_model = pm.auto_arima(
    ts,
    seasonal=True, m=12,           # 12-month seasonality
    d=None, D=None,                # let it determine differencing
    stepwise=True,                  # fast search (vs exhaustive)
    information_criterion='aic',    # minimise AIC
    max_p=3, max_q=3,
    max_P=2, max_Q=2,
    trace=True,                     # print search progress
    error_action='ignore',
    suppress_warnings=True
)
print('
Best model:', auto_model.order, auto_model.seasonal_order)

# --- Manual SARIMA fit ---
train = ts[:'2025-12']
test  = ts['2026-01':]

model = SARIMAX(train,
                order=(1, 1, 1),
                seasonal_order=(1, 1, 1, 12),
                trend='c',
                enforce_stationarity=False,
                enforce_invertibility=False)
fit = model.fit(disp=False)
print(fit.summary())

# Residual diagnostics
residuals = fit.resid
lb_test = acorr_ljungbox(residuals, lags=[12], return_df=True)
print('Ljung-Box p-value (should be > 0.05):', lb_test['lb_pvalue'].values[0])

# Forecast and confidence intervals
forecast = fit.get_forecast(steps=len(test))
pred_mean = forecast.predicted_mean
pred_ci   = forecast.conf_int(alpha=0.05)

# Plot forecast
plt.figure(figsize=(12, 5))
plt.plot(train[-24:], label='Train', color='steelblue')
plt.plot(test, label='Actual', color='green')
plt.plot(pred_mean, label='Forecast', color='red', linestyle='--')
plt.fill_between(pred_ci.index, pred_ci.iloc[:,0], pred_ci.iloc[:,1],
                 alpha=0.2, color='red', label='95% CI')
plt.legend(); plt.title('SARIMA Forecast'); plt.tight_layout(); plt.show()

Prophet — Scalable Forecasting for Business Time Series

Prophet (Facebook/Meta, 2017) is designed for business time series with strong seasonal patterns, holidays, and occasional structural breaks. It fits an additive model: y(t) = trend(t) + seasonality(t) + holidays(t) + noise. Trend is modelled as piecewise linear or logistic growth with automatically detected changepoints. Seasonality uses Fourier series decomposition. Strengths: handles missing data gracefully, robust to outliers, automatically detects changepoints, requires minimal parameter tuning, and produces uncertainty intervals. Weaknesses: assumes additive structure, less powerful than ML methods for complex multivariate forecasting.

from prophet import Prophet
from prophet.diagnostics import cross_validation, performance_metrics
import pandas as pd

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

m = Prophet(
    seasonality_mode='multiplicative',  # for series where amplitude scales with level
    yearly_seasonality=True,
    weekly_seasonality=True,
    daily_seasonality=False,
    changepoint_prior_scale=0.05,       # flexibility of trend (0.001=rigid, 0.5=flexible)
    seasonality_prior_scale=10,
    holidays_prior_scale=10,
    interval_width=0.95                 # 95% confidence intervals
)

# Add custom seasonality (e.g., quarterly for business data)
m.add_seasonality(name='quarterly', period=91.25, fourier_order=5)

# Add country holidays
m.add_country_holidays(country_name='IN')

m.fit(df_prophet[df_prophet['ds'] < '2026-01-01'])

# Make future dataframe and forecast
future   = m.make_future_dataframe(periods=12, freq='MS')
forecast = m.predict(future)
m.plot(forecast); plt.show()
m.plot_components(forecast); plt.show()

# Cross-validation: rolling-window backtesting
cv_results = cross_validation(
    m, initial='730 days', period='90 days', horizon='180 days'
)
metrics = performance_metrics(cv_results)
print(metrics[['horizon', 'mae', 'rmse', 'mape']].to_string())

For LSTM-based sequence forecasting (which outperforms ARIMA and Prophet when you have large amounts of training data, multiple input series, or complex non-linear patterns), our Neural Network Architectures guide covers LSTM architecture in detail. For the full set of time series interview questions covering stationarity, ARIMA orders, evaluation metrics (MAPE, RMSE, MASE), and model selection, our Time Series Forecasting Interview Q&A has 35+ questions. For feature engineering specific to time series — lag features, rolling statistics, Fourier features — our Feature Engineering guide covers all standard patterns. For the data pipelines that feed time series models with fresh data, our ETL Pipelines guide covers Airflow scheduling and dbt transformations.