Forecasting Australian Beer Production using ARIMA & SARIMA Models in Python

Understanding the Dataset

The dataset consists of monthly beer production in Australia, indexed by time (Month). Our first step is to load and preprocess the dataset:

#Import Necessary Libraries
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.arima.model import ARIMA
from sklearn.metrics import mean_absolute_error, mean_squared_error

#Load dataset and set Month as index
df = pd.read_csv('monthly-beer-production-in-austr.csv', parse_dates=['Month'], index_col='Month')
print(df.head())

To evaluate model performance, we split the data into 70% training and 30% testing:

# Step 3 - Split the Dataset (70% Train, 30% Test)
train, test = df.iloc[:int(len(df) * 0.7)], df.iloc[int(len(df) * 0.7):]

# Print shapes to confirm split
print(train.shape, test.shape)

Checking Stationarity

A key assumption of ARIMA is that the time series should be stationary (i.e., its statistical properties remain constant over time). We first start by a visual graphical analysis of the Time Series Data before moving on to formal tests of stationarity.

# Step 4 - Plot the training data to check for trends or seasonality
plt.figure(figsize=(10, 5))
plt.plot(train, label="Training Data", color='blue')
plt.xlabel("Year")
plt.ylabel("Monthly Beer Production")
plt.title("Monthly Beer Production Over Time")
plt.legend()
plt.show()

Time Series Plot

We now use the Augmented Dickey-Fuller (ADF) test to check stationarity:

# Step 5 - Perform the ADF Test for Stationarity
result = adfuller(train)
print('ADF Statistic:', result[0])
print('p-value:', result[1])
print('Critical Values:', result[4])
if result[1] <= 0.05:
    print("The time series is stationary")
else:
    print("The time series is non-stationary")

Results: ADF Statistic: -1.43598433169934 p-value: 0.5649084390929975 Critical Values: {‘1%’: -3.451281394993741, ‘5%’: -2.8707595072926293, ‘10%’: -2.571682118921643} The time series is non-stationary

Differencing to Achieve Stationarity

Differencing is a technique that transforms the dataset into a stationary series by subtracting consecutive observations.

# Step 6 - Differencing to Make Time Series Stationary
train_diff = train.diff().dropna()

# Perform ADF Test Again
result_diff = adfuller(train_diff)

# Print ADF test results
print('ADF Statistic (After Differencing):', result_diff[0])
print('p-value:', result_diff[1])
print('Critical Values:', result_diff[4])

# Check if the data is now stationary
if result_diff[1] <= 0.05:
    print("The Differenced Time Series is now Stationary")
else:
    print("The Differenced Time Series is still Non-Stationary")

Results: ADF Statistic (After Differencing): -5.012526073359135 p-value: 2.0960589659990194e-05 Critical Values: {‘1%’: -3.451281394993741, ‘5%’: -2.8707595072926293, ‘10%’: -2.571682118921643} The Differenced Time Series is now Stationary

Visually Confirming Stationarity

# Step 7 - Visualizing Original vs Differenced Data [For Visually Confirming Stationarity]
plt.figure(figsize=(12, 5))

# Original Data
plt.subplot(1, 2, 1)
plt.plot(train, label="Original Data", color="blue")
plt.title("Original Time Series")
plt.legend()

# Differenced Data
plt.subplot(1, 2, 2)
plt.plot(train_diff, label="Differenced Data (1st Order)", color="red")
plt.title("After Differencing (Stationary)")
plt.legend()

plt.tight_layout()
plt.show()

Original vs Differenced Time Series

Identifying ARIMA Parameters (p, d, q)

To determine the ARIMA order parameters (p, d, q), we plot Autocorrelation Function (ACF) and Partial Autocorrelation Function (PACF):

# Step 8 - Plot ACF and PACF to Determine ARIMA Orders (p, d, q)
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

plt.figure(figsize=(12, 5))

# Plot ACF (to determine q)
plt.subplot(1, 2, 1)
plot_acf(train_diff, lags=30, ax=plt.gca())
plt.title("Autocorrelation Function (ACF)")

# Plot PACF (to determine p)
plt.subplot(1, 2, 2)
plot_pacf(train_diff, lags=30, ax=plt.gca(), method='ywm')
plt.title("Partial Autocorrelation Function (PACF)")

plt.tight_layout()
plt.show()

ACF and PACF Plots

Fitting the ARIMA Model with parameters 1, 1, 1

#Step 9 - Fit ARIMA Model
model = ARIMA(train, order=(1, 1, 1))
arima_result = model.fit()

# Print summary of the model
print(arima_result.summary())

Forecasting and Evaluating ARIMA

We forecast the next 143 months (same length as the test set) and evaluate performance using Mean Absolute Error (MAE) and Root Mean Squared Error (RMSE):

# Step 10 - Forecasting
forecast_steps = 143
forecast = arima_result.forecast(steps=forecast_steps)

print(forecast.head(10))

# Step 11 - ARIMA Model Performance

from sklearn.metrics import mean_absolute_error

mae = mean_absolute_error(test, forecast)
print(f"Mean Absolute Error (MAE): {mae:.2f}")

from sklearn.metrics import mean_squared_error
import numpy as np

rmse = np.sqrt(mean_squared_error(test, forecast))
print(f"Root Mean Squared Error (RMSE): {rmse:.2f}")

Results: Mean Absolute Error (MAE): 17.98 Root Mean Squared Error (RMSE): 21.71

Now we visualize the forecast against actual values:

# Step 12 - Plotting the ARIMA Model
plt.figure(figsize=(12, 6))
plt.plot(train.index, train, label="Training Data", color="blue")
plt.plot(test.index, test, label="Actual Test Data", color="yellow")
plt.plot(test.index, forecast, label="Forecasted Data", color="black", linestyle="dashed")
plt.legend()
plt.title("ARIMA Forecast vs Actual")
plt.xlabel("Year")
plt.ylabel("Beer Production")
plt.show()

ARIMA Forecast

As you can see the ARIMA model does not capture seasonality well, we move to SARIMA.

Building the SARIMA Model

SARIMA extends ARIMA by adding seasonal components. The order notation is (p, d, q) x (P, D, Q, s) where s is the seasonal period (12 for monthly data)

# Step 13 - Fitting a SARIMA Model

from statsmodels.tsa.statespace.sarimax import SARIMAX

sarima_model = SARIMAX(train, order=(1,1,1), seasonal_order=(1,1,1,12))  # 12 for monthly data
sarima_result = sarima_model.fit()
print(sarima_result.summary())

Forecasting the SARIMA Model and Evaluating its Performance

# Step 14 - Setting the Steps for Forecast
forecast_steps = 143
sarima_forecast = sarima_result.get_forecast(steps=forecast_steps)

# Extract forecasted values
forecast_mean = sarima_forecast.predicted_mean

# Step 15 - SARIMA Model Performance

mae_sarima = mean_absolute_error(test, forecast_mean)
print(f"Mean Absolute Error (MAE): {mae_sarima:.2f}")

# Compute RMSE
rmse_sarima = np.sqrt(mean_squared_error(test, forecast_mean))
print(f"Root Mean Squared Error (RMSE): {rmse_sarima:.2f}")

Results: Mean Absolute Error (MAE): 12.20 Root Mean Squared Error (RMSE): 14.98

Visualizing SARIMA vs. Actual Data

# Step 16 - Plotting the predicted SARIMA values against the Actual Time Series Values

plt.figure(figsize=(12, 6))

# Plot actual values for the forecasted period
plt.plot(test.index, test, label="Actual", color='blue')

# Plot predicted values for the forecasted period
plt.plot(forecast_mean.index, forecast_mean, label="Predicted", color='red')

plt.xlabel("Time")
plt.ylabel("Beer Production")
plt.title("SARIMA:- Actual vs. Predicted (Forecasted Period Only)")
plt.legend()
plt.show()

SARIMA Forecast

As we can see, the SARIMA Model performed better than the ARIMA Model ,further tuning of parameters would give us an even better fitting model.