Auto ARIMA hyperparameter search
In this blog, I will present the method for automatised search of the key parameters for (S)ARIMA forecasting models.
Introduction
This developed method for auto hyperparameter search for forecasting, could take a lot of time, but it can pay off because when you find the best parameter of your forecasting model, than it will save the time and raise the precision of your predictions. Also, manual trying could take you the most time, but this method could be helpful in some situation, having in mind that understanding why we should use exact parameters in forecasting model could be of crucial importance for better understanding the business problem. Nevertheless, this approach could fasten your exploration time and help you to get better understanding of your problem.
The hyperparameter we will tune in forecasting model of (S)ARIMA are seasonality parameter (S), autoregressive parameter (AR), differencing parameter (I) and moving average (MA). You can read more about this parameters in my previous blog if needed (here).
Example with CO2 data
We will check and plot trend and seasonality of the timeseries about CO2 emission in last 50 years and we see increasing trend and clear seasonal pattern due to the heating periods in the winter (full code with data available here).
import pandas as pd
import numpy as np
from statsmodels.tsa.seasonal import STL
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller
import statsmodels.api as smplt.rc('figure',figsize=(10,6))
plt.rc('font',size=13)
stl = STL(co2, seasonal=7)
res = stl.fit()
fig = res.plot()
We should check the stationarity of the CO2 emissions through the years, after differencing with created functions that plots rolling mean (for 12 data points) and standard deviation of our differenced time series.
def test_stationarity(timeseries, rolling=12):
#Determing rolling statistics
rolmean = timeseries.rolling(rolling).mean()
rolstd = timeseries.rolling(rolling).std()
#Plot rolling statistics:
plt.figure(figsize=(10, 5))
orig = plt.plot(timeseries, color='blue',label='Original')
mean = plt.plot(rolmean, color='red', label='Rolling Mean')
std = plt.plot(rolstd, color='black', label = 'Rolling Std')
plt.title('Power consumption Old data')
plt.xlabel('Time - periods(30s)')
plt.ylabel('Power consumption in Watts')
plt.legend(loc='best')
plt.title('Rolling Mean & Standard Deviation')
plt.show()
#Perform Dickey-Fuller test:
print ('Results of Dickey-Fuller Test:')
dftest = adfuller(timeseries, autolag='AIC')
dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
for key,value in dftest[4].items():
dfoutput['Critical Value (%s)'%key] = value
print (dfoutput)
if dfoutput['p-value'] < 0.05:
print('The time series is stationary at 95% level of confidence')
else:
print('The time series is not stationary at 95% level of confidence')co2_diff = co2 - co2.shift(1)
co2_diff = co2_diff.dropna()
test_stationarity(co2_diff, rolling=12)
Before we go to the auto arima search, we will check the correlogram for 40 time spots to see how are data are correlated across lags and from this we could confirm high AR parameter as well as strong seasonality pattern (partial correlation which captures the correlation of residuals after regression of all other variables).
fig, axes = plt.subplots(1, 2, figsize=(15,4))
fig = sm.graphics.tsa.plot_acf(co2,
lags=40,
ax=axes[0])
fig = sm.graphics.tsa.plot_pacf(co2,
lags=40,
ax=axes[1])
Auto tuning
With the next part of the code, we have created several combinations of p (AR), d(I) and q(MA) parameter. In order to relax the search, we defined the seasonal (s) parameter to 12 as we are sure there is yearly pattern. For evaluating the goodness of the model we use Akaike information criteria (aic). This criteria is crucial for statistical model of some data, by calculating the difference of the k — the number of estimated parameters in the model and L be the maximum value of the likelihood function for the model. This means that the model with best AIC will have the optimum value of likelihood and the smallest possible number of the parameters, thus enabling simplicity and efficiency of the model (parsimony) which leads to the optimal performance. The AIC value of the model is the following:
AIC=2k−2ln(L)
import itertools
import warnings
warnings.filterwarnings('ignore')
# Grid Search
p = d = q = range(0,3) # p, d, and q can be either 0, 1, or 2
pdq = list(itertools.product(p,d,q)) # gets all possible combinations of p, d, and q
p2 = d2 = q2 = range(0, 2) # second set of p's, d's, and q's for seasonal parameters
pdq2 = list(itertools.product(p2,d2,q2)) # similar to code above but for seasonal parameters
s = 12 # here I use twelve but the number here is representative of the periodicity of the seasonal cycle
pdqs2 = [(c[0], c[1], c[2], s) for c in pdq2]
combs = {}
aics = []
# Grid Search Continued
for combination in pdq:
for seasonal_combination in pdqs2:
try:
model = sm.tsa.statespace.SARIMAX(co2, order=combination, seasonal_order=seasonal_combination,
enforce_stationarity=True,
enforce_invertibility=True)
model = model.fit(disp=False)
combs.update({model.aic : [combination, seasonal_combination]})
aics.append(model.aic)
except:
continue
best_aic = min(aics)Now we can see what hyperparameters for forecasting will give the best predictions.
print('best aic is: ', round(best_aic, 3))
print(40*'==')
print ('ARIMA parameters: ', '\n', 'AR: ', combs[best_aic][0][0], '\n', 'I: ',combs[best_aic][0][1], '\n', 'MA: ',combs[best_aic][0][2])
print('Seasonal parameters:', combs[best_aic][1])
Congratulations! With this hyperparameters, we can do a train test split and make a predictions for a test set and than compare the predictions and test results.
Forecasting with selected hyperparameters
After creation of train — test split function, we make actual forecast:
def train_test_split(timeseries, lags_for_prediction=12):
split=len(timeseries)-lags_for_prediction
train=timeseries[:split]
test=timeseries[split:]
return train, testtrain_series, test_series = train_test_split(co2, 12)
The good way of checking the goodness of the model is plotting the forecasted versus test data. In real application, the real forecast will go more into the future, but this step assures us that we have a good model.
def forecasting (p,d,q, season, lags_for_forecast, train_series):
model = sm.tsa.statespace.SARIMAX(train_series, order=(p,d,q), seasonal_order=(p,d,q,season),
simple_differencing=0, #if True time series provided as endog is literally differenced and an ARMA model is fit to the resulting new time series
enforce_stationarity=True,
enforce_invertibility=False)
fitted = model.fit(disp=-1)
# Forecast
forecast = fitted.forecast(lags_for_forecast)
# Plot
plt.figure(figsize=(12,5), dpi=100)
plt.plot(train_series, color='blue', label='train')
plt.plot(test_series, color='green', label='test', alpha=0.6)
plt.plot(forecast, color='red', label='forecast')
plt.title('Forecast vs Actuals')
plt.legend(loc='upper left', fontsize=8)
plt.show()
RSS=np.sqrt(sum(forecast.values-test_series.values.reshape(-1))**2)/lags_for_forecast
print("\n", '\033[1m' +'Root Squared Error (RSS) of SARIMAX model(p,d,q)(p,d,q,s)' + '\033[0m',(p,d,q),season,':', round(RSS, 3),"\n")
print(fitted.summary())forecasting (2,1,1, 12, 12, train_series)
As we can see, the auto tuning works quite good. The summary option of statsmodels gives us more in depth view about statistical significance of our parameters. The very small root mean square error (0.313) gives us additional assurance that we have reliable and good model. That would be it. In some cases, old traditional ARIMA could beat the complexed LSTM :)
