[Python] Why is SARIMAX expensive to apply over high-frequency data in Python? (epochs=5mins)...

I’ve been experimenting with SARIMAX using statsmodels package to model time series data and understand seasonality and non-stationarity. My data is collected with a high frequency of 5-minute epochs, and I’m facing significant challenges:

• 
  Performance Issues: Applying SARIMAX with what I believe is the correct configuration (e.g., daily seasonality with 288 observations per day) is extremely slow. Sometimes, the process crashes due to memory or computational limits.
  
  
• 
  Suitability for High-Frequency Data: I suspect SARIMAX might not be the right choice for high-frequency datasets. Is this true? If so, why is SARIMAX particularly inefficient for such data?
  

Fixes or Alternatives:

Important Note: I am not interested in resampling or downsampling the data because it would cause a loss of valuable information that is critical for my analysis.

Questions:

• Are there any techniques to optimize SARIMAX for high-frequency data without losing information?
• Alternatively, are there better models specifically designed for high-frequency data that still allow me to understand time components like seasonality and trends?

I’m open to suggestions or insights from anyone who has tackled similar challenges with high-frequency data. Thanks in advance for your guidance!

I generated a reproducible example:

# Import additional required libraries
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.statespace.sarimax import SARIMAX
from sklearn.metrics import mean_absolute_error, mean_squared_error
from datetime import datetime, timedelta

# Generate timestamps for 30 consecutive days with 5-minute intervals
start_time = datetime.now().replace(hour=0, minute=0, second=0, microsecond=0)
timestamps = [start_time + timedelta(minutes=5 * i) for i in range(30 * 24 * 12)]

# Generate 'avgcpu' data with seasonality
np.random.seed(42) # For reproducibility
base_pattern = np.sin(np.linspace(0, 2 * np.pi, len(timestamps) // 30)) * 7.5 + 27.5
avgcpu = np.tile(base_pattern, 30) + np.random.normal(0, 2, len(timestamps))
avgcpu = np.clip(avgcpu, 20, 35)

# Add outliers for the 12th day
outlier_start = 12 * 24 * 12 # Start index for the 12th day
outlier_end = outlier_start + 12 # 1-hour duration (12 points)
avgcpu[outlier_startutlier_end] = np.random.uniform(70, 80, outlier_end - outlier_start)

# Add outliers for the 23th day
outlier_start = 23 * 24 * 12 # Start index for the 12th day
outlier_end = outlier_start + 12 # 1-hour duration (12 points)
avgcpu[outlier_startutlier_end] = np.random.uniform(5, 10, outlier_end - outlier_start)

# Create DataFrame
data = pd.DataFrame({
"timestamp": timestamps,
"avgcpu": avgcpu
}).sort_values("timestamp")

# Extract training and testing series
train_series = train_data["avgcpu"]
test_series = test_data["avgcpu"]

# SARIMAX configuration
# Seasonal period is 288 (daily seasonality with 5-minute epochs)
seasonal_period = 288

model = SARIMAX(
train_series,
order=(1, 1, 1), # Simplified initial setup for ARIMA terms
seasonal_order=(1, 1, 1, seasonal_period), # Seasonal ARIMA terms
enforce_stationarity=False,
enforce_invertibility=False
)

# Fit the model
sarimax_result = model.fit(disp=False)

# Forecast for the test set
forecast_result = sarimax_result.get_forecast(steps=test_size, alpha=0.05) # 95% CI
forecast_mean = forecast_result.predicted_mean
forecast_ci = forecast_result.conf_int() # Extracts the confidence intervals

# Calculate metrics
mae = mean_absolute_error(test_series, forecast_mean)
mse = mean_squared_error(test_series, forecast_mean)

# Plot the results
plt.figure(figsize=(14, 7))
plt.plot(data["timestamp"], data["avgcpu"], label="True Data")
plt.plot(test_data["timestamp"], forecast_mean, label="Predicted Mean", color="orange")
plt.fill_between(
test_data["timestamp"],
forecast_ci.iloc[:, 0],
forecast_ci.iloc[:, 1],
color="orange",
alpha=0.3,
label="Prediction Interval"
)
plt.title("SARIMAX Prediction with Test Data")
plt.xlabel("Timestamp")
plt.ylabel("Avg CPU")
plt.ylim(0, 100)
plt.xticks(rotation=45)
plt.legend(title=f"MAE: {mae:.2f}, MSE: {mse:.2f}")
plt.grid()
plt.show()

Continue reading...