How to Implement SARIMA in Python

Key Insights

SARIMA models excel at forecasting time series data with seasonal patterns by combining non-seasonal ARIMA components (p,d,q) with seasonal components (P,D,Q,s), where ’s’ represents the seasonal period
Parameter selection is critical—use ACF/PACF plots for manual identification or leverage automated grid search with AIC/BIC criteria to find optimal model configurations
Always validate SARIMA models with proper train/test splits and walk-forward validation rather than relying solely on in-sample fit statistics, as overfitting is a common pitfall

Introduction to SARIMA Models

SARIMA (Seasonal AutoRegressive Integrated Moving Average) models are the go-to solution for time series forecasting when your data exhibits both trend and seasonal patterns. Unlike basic ARIMA models, SARIMA explicitly accounts for seasonality—those recurring patterns that appear at regular intervals like weekly sales spikes, quarterly revenue cycles, or monthly temperature fluctuations.

The SARIMA notation looks intimidating at first: SARIMA(p,d,q)(P,D,Q,s). Here’s what each parameter means:

p: Number of autoregressive terms (non-seasonal)
d: Number of differencing operations needed for stationarity (non-seasonal)
q: Number of moving average terms (non-seasonal)
P: Number of seasonal autoregressive terms
D: Number of seasonal differencing operations
Q: Number of seasonal moving average terms
s: Length of the seasonal cycle (12 for monthly data with yearly seasonality, 7 for daily data with weekly patterns)

Use SARIMA when you have time series data with clear seasonal patterns that repeat at known intervals. It outperforms simple ARIMA when seasonality is present and provides more accurate forecasts than naive seasonal decomposition methods.

Data Preparation and Exploration

Before building any SARIMA model, you need to understand your data’s characteristics. Let’s work with a practical example using retail sales data.

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

# Load your time series data
df = pd.read_csv('retail_sales.csv', parse_dates=['date'], index_col='date')
ts = df['sales']

# Visualize the time series
plt.figure(figsize=(12, 6))
plt.plot(ts)
plt.title('Monthly Retail Sales')
plt.xlabel('Date')
plt.ylabel('Sales')
plt.grid(True)
plt.show()

Check for stationarity using the Augmented Dickey-Fuller test. SARIMA requires stationary data, which means constant mean and variance over time.

def check_stationarity(timeseries):
    result = adfuller(timeseries.dropna())
    print('ADF Statistic:', result[0])
    print('p-value:', result[1])
    print('Critical Values:')
    for key, value in result[4].items():
        print(f'\t{key}: {value}')
    
    if result[1] <= 0.05:
        print("Data is stationary")
    else:
        print("Data is non-stationary")

check_stationarity(ts)

Visualize ACF and PACF plots to identify patterns and potential parameter values:

fig, axes = plt.subplots(1, 2, figsize=(14, 4))
plot_acf(ts.dropna(), lags=40, ax=axes[0])
plot_pacf(ts.dropna(), lags=40, ax=axes[1])
plt.tight_layout()
plt.show()

# Decompose to visualize trend and seasonality
decomposition = seasonal_decompose(ts, model='multiplicative', period=12)
fig = decomposition.plot()
fig.set_size_inches(12, 8)
plt.show()

Determining SARIMA Parameters

Parameter selection makes or breaks your SARIMA model. You have two approaches: manual selection using ACF/PACF plots or automated grid search.

For manual selection, examine the ACF and PACF plots:

p: Significant lags in PACF before cutoff
q: Significant lags in ACF before cutoff
P and Q: Significant lags at seasonal intervals

Here’s an automated grid search approach that’s more reliable:

import itertools
from statsmodels.tsa.statespace.sarimax import SARIMAX
import warnings
warnings.filterwarnings('ignore')

# Define parameter ranges
p = d = q = range(0, 3)
pdq = list(itertools.product(p, d, q))
seasonal_pdq = [(x[0], x[1], x[2], 12) for x in pdq]

best_aic = np.inf
best_params = None
best_seasonal_params = None

for param in pdq:
    for param_seasonal in seasonal_pdq:
        try:
            model = SARIMAX(ts,
                          order=param,
                          seasonal_order=param_seasonal,
                          enforce_stationarity=False,
                          enforce_invertibility=False)
            results = model.fit(disp=False)
            
            if results.aic < best_aic:
                best_aic = results.aic
                best_params = param
                best_seasonal_params = param_seasonal
                
        except:
            continue

print(f'Best SARIMA{best_params}x{best_seasonal_params} - AIC:{best_aic}')

Alternatively, use pmdarima for automated parameter selection:

from pmdarima import auto_arima

auto_model = auto_arima(ts, 
                       start_p=0, start_q=0,
                       max_p=3, max_q=3,
                       seasonal=True,
                       start_P=0, start_Q=0,
                       max_P=3, max_Q=3,
                       m=12,  # seasonal period
                       d=None, D=None,  # let auto_arima determine
                       trace=True,
                       error_action='ignore',
                       suppress_warnings=True,
                       stepwise=True)

print(auto_model.summary())

Building and Fitting the SARIMA Model

Once you’ve identified optimal parameters, build and fit your model using statsmodels:

from statsmodels.tsa.statespace.sarimax import SARIMAX

# Split data into train and test sets
train_size = int(len(ts) * 0.8)
train, test = ts[:train_size], ts[train_size:]

# Create and fit the SARIMA model
model = SARIMAX(train,
               order=(1, 1, 1),
               seasonal_order=(1, 1, 1, 12),
               enforce_stationarity=False,
               enforce_invertibility=False)

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

Examine diagnostic plots to validate model assumptions:

# Diagnostic plots
fitted_model.plot_diagnostics(figsize=(12, 8))
plt.tight_layout()
plt.show()

Look for:

Residuals should resemble white noise (no patterns)
Q-Q plot should show residuals along the diagonal line
Correlogram should show no significant autocorrelation
Ljung-Box test p-values should be above 0.05

Forecasting and Model Validation

Generate forecasts with confidence intervals:

# Forecast
forecast_steps = len(test)
forecast = fitted_model.forecast(steps=forecast_steps)
forecast_ci = fitted_model.get_forecast(steps=forecast_steps).conf_int()

# Plot results
plt.figure(figsize=(12, 6))
plt.plot(train.index, train, label='Training Data')
plt.plot(test.index, test, label='Actual', color='green')
plt.plot(test.index, forecast, label='Forecast', color='red')
plt.fill_between(test.index,
                forecast_ci.iloc[:, 0],
                forecast_ci.iloc[:, 1],
                color='pink', alpha=0.3)
plt.legend()
plt.title('SARIMA Forecast vs Actual')
plt.show()

Calculate accuracy metrics:

from sklearn.metrics import mean_squared_error, mean_absolute_error

rmse = np.sqrt(mean_squared_error(test, forecast))
mae = mean_absolute_error(test, forecast)
mape = np.mean(np.abs((test - forecast) / test)) * 100

print(f'RMSE: {rmse:.2f}')
print(f'MAE: {mae:.2f}')
print(f'MAPE: {mape:.2f}%')

Implement walk-forward validation for robust evaluation:

def walk_forward_validation(data, train_size, order, seasonal_order):
    predictions = []
    actuals = []
    
    for i in range(train_size, len(data)):
        train = data[:i]
        test_point = data[i]
        
        model = SARIMAX(train, order=order, seasonal_order=seasonal_order)
        fitted = model.fit(disp=False)
        forecast = fitted.forecast(steps=1)
        
        predictions.append(forecast.iloc[0])
        actuals.append(test_point)
    
    return np.array(actuals), np.array(predictions)

actuals, preds = walk_forward_validation(ts, train_size, (1,1,1), (1,1,1,12))
wf_rmse = np.sqrt(mean_squared_error(actuals, preds))
print(f'Walk-Forward RMSE: {wf_rmse:.2f}')

Practical Tips and Common Pitfalls

Handle missing data properly before fitting SARIMA models. Don’t just drop missing values—interpolate or use forward-fill depending on your domain:

# Linear interpolation for missing values
ts_clean = ts.interpolate(method='linear')

# Or forward fill for short gaps
ts_clean = ts.fillna(method='ffill', limit=2)

Watch for overfitting. Complex models with high parameter values fit training data perfectly but fail on new data. Always validate on held-out test sets and prefer simpler models when AIC/BIC values are close.

SARIMA assumptions matter. If your data has multiple seasonal patterns (daily and weekly), consider TBATS or Prophet instead. If variance changes over time (heteroscedasticity), apply log transformation:

ts_log = np.log(ts)
# Fit model on log-transformed data
# Transform predictions back: np.exp(forecast)

Computational cost increases rapidly with parameter complexity and data size. For datasets with thousands of observations, consider using method='nm' in the fit function or sampling your data for initial parameter exploration.

Compare multiple models before settling on one:

models = {
    'SARIMA(1,1,1)(1,1,1,12)': ((1,1,1), (1,1,1,12)),
    'SARIMA(2,1,2)(1,1,1,12)': ((2,1,2), (1,1,1,12)),
}

for name, params in models.items():
    model = SARIMAX(train, order=params[0], seasonal_order=params[1])
    fitted = model.fit(disp=False)
    forecast = fitted.forecast(steps=len(test))
    rmse = np.sqrt(mean_squared_error(test, forecast))
    print(f'{name} - RMSE: {rmse:.2f}, AIC: {fitted.aic:.2f}')

SARIMA works best with regular, clean seasonal patterns. If your data has structural breaks, regime changes, or irregular seasonality, consider ensemble methods or more flexible models like machine learning approaches with engineered features.