# python time series prediction -- SARIMA

Posted by Twister1004 on Sat, 07 Mar 2020 09:12:27 +0100

SARIMA(p,d,q)(P,D,Q,s) seasonal autoregressive moving average model has seven structural parameters

AR(p) autoregressive model, i.e. use oneself to regress oneself.

The basic assumption is that the current sequence value depends on its previous value with time delay. The maximum lag in the model is recorded as p

To determine the initial p, you need to look at the PACF graph and find the largest significant lag, after which most of the other lags become irrelevant.

MA(q) moving average model is to model the error of time series, and assume that the current error depends on the error with delay. The initial value can be found on the ACF graph.

Combining the above two methods: AR(p)+MA(q)=ARMA(p,q)AR(p)+MA(q)=ARMA(p,q)AR(p)+MA(q)=ARMA(p,q), which is the autoregressive moving average model

Remaining parameters:

I(d)

• order of integration. This is simply the number of nonseasonal differences needed to make the series stationary. In our case, it's just 1 because we used first differences.
Adding this letter to the four gives us the ARIMA
model which can handle non-stationary data with the help of nonseasonal differences. Great, one more letter to go!

S(s)

• this is responsible for seasonality and equals the season period length of the series
With this, we have three parameters: (P,D,Q)

P

• order of autoregression for the seasonal component of the model, which can be derived from PACF. But you need to look at the number of significant lags, which are the multiples of the season period length. For example, if the period equals 24 and we see the 24-th and 48-th lags are significant in the PACF, that means the initial P
should be 2.

Q

• similar logic using the ACF plot instead.

D

• order of seasonal integration. This can be equal to 1 or 0, depending on whether seasonal differeces were applied or not.

OK, let's use SARIMA to model A kind of

# Import package

```import warnings                                  # do not disturbe mode
warnings.filterwarnings('ignore')

import numpy as np                               # vectors and matrices
import pandas as pd                              # tables and data manipulations
import matplotlib.pyplot as plt                  # plots
import seaborn as sns                            # more plots

from dateutil.relativedelta import relativedelta # working with dates with style
from scipy.optimize import minimize              # for function minimization

import statsmodels.formula.api as smf            # statistics and econometrics
import statsmodels.tsa.api as smt
import statsmodels.api as sm
import scipy.stats as scs

from itertools import product                    # some useful functions
from tqdm import tqdm_notebook

# Importing everything from forecasting quality metrics
from sklearn.metrics import r2_score, median_absolute_error, mean_absolute_error
from sklearn.metrics import median_absolute_error, mean_squared_error, mean_squared_log_error
```
```# MAPE
def mean_absolute_percentage_error(y_true, y_pred):
return np.mean(np.abs((y_true - y_pred) / y_true)) * 100

def tsplot(y, lags=None, figsize=(12, 7), style='bmh'):
"""
Plot time series, its ACF and PACF, calculate Dickey–Fuller test

y - timeseries
lags - how many lags to include in ACF, PACF calculation
"""

if not isinstance(y, pd.Series):
y = pd.Series(y)

with plt.style.context(style):
fig = plt.figure(figsize=figsize)
layout = (2, 2)
ts_ax = plt.subplot2grid(layout, (0, 0), colspan=2)
acf_ax = plt.subplot2grid(layout, (1, 0))
pacf_ax = plt.subplot2grid(layout, (1, 1))

y.plot(ax=ts_ax)
ts_ax.set_title('Time Series Analysis Plots\n Dickey-Fuller: p={0:.5f}'.format(p_value))
smt.graphics.plot_acf(y, lags=lags, ax=acf_ax)
smt.graphics.plot_pacf(y, lags=lags, ax=pacf_ax)
plt.tight_layout()
```

```ads = pd.read_csv('ads.csv', index_col=['Time'], parse_dates=['Time'])
``` ```plt.figure(figsize=(18, 6))
plt.grid(True)
plt.show()
``` ```tsplot(ads.Ads, lags=60)
``` ```# The seasonal difference
``` ```ads_diff = ads_diff - ads_diff.shift(1)
``` ## SARIMA parameter

• p is 4 because the last significant lag on the PACF, after which most of the others were not significant.
• d is 1 because of the first order difference
• q according to ACF, it should be about 4
• P may be 2, because the lag of 24 and 48 is important in PACF
• D is 1 because we do the seasonal difference
• Q may be 1, the 24th lag of ACF is obvious, while the 48th lag is not
```# setting initial values and some bounds for them
ps = range(2, 5)
d=1
qs = range(2, 5)
Ps = range(0, 2)
D=1
Qs = range(0, 2)
s = 24 # season length is still 24

# creating list with all the possible combinations of parameters
parameters = product(ps, qs, Ps, Qs)
parameters_list = list(parameters)
len(parameters_list)  # 36
```

```def optimizeSARIMA(parameters_list, d, D, s):
"""Return dataframe with parameters and corresponding AIC

parameters_list - list with (p, q, P, Q) tuples
d - integration order in ARIMA model
D - seasonal integration order
s - length of season
"""

results = []
best_aic = float("inf")

for param in tqdm_notebook(parameters_list):
# we need try-except because on some combinations model fails to converge
try:
seasonal_order=(param, D, param, s)).fit(disp=-1)
except:
continue
aic = model.aic
# saving best model, AIC and parameters
if aic < best_aic:
best_model = model
best_aic = aic
best_param = param
results.append([param, model.aic])

result_table = pd.DataFrame(results)
result_table.columns = ['parameters', 'aic']
# sorting in ascending order, the lower AIC is - the better
result_table = result_table.sort_values(by='aic', ascending=True).reset_index(drop=True)

return result_table
```
```warnings.filterwarnings("ignore")
result_table = optimizeSARIMA(parameters_list, d, D, s)
'''
parameters	    aic
0	(2, 3, 1, 1)	3888.642174
1	(3, 2, 1, 1)	3888.763568
2	(4, 2, 1, 1)	3890.279740
3	(3, 3, 1, 1)	3890.513196
4	(2, 4, 1, 1)	3892.302849
5	(4, 3, 1, 1)	3892.322855
6	(3, 4, 1, 1)	3893.762846
7	(4, 4, 1, 1)	3894.327967
8	(2, 2, 1, 1)	3894.798147
9	(2, 3, 0, 1)	3897.170902
10	(3, 2, 0, 1)	3897.815032
11	(4, 2, 0, 1)	3899.073591
12	(3, 3, 0, 1)	3899.165271
13	(3, 4, 0, 1)	3900.500309
14	(2, 4, 0, 1)	3900.502494
15	(4, 3, 0, 1)	3901.255700
16	(4, 4, 0, 1)	3902.650501
17	(2, 2, 0, 1)	3903.905714
18	(2, 3, 1, 0)	3909.281188
19	(3, 2, 1, 0)	3909.502838
20	(4, 2, 1, 0)	3910.927759
21	(3, 3, 1, 0)	3911.192654
22	(3, 4, 1, 0)	3911.344351
23	(2, 4, 1, 0)	3911.809710
24	(4, 4, 1, 0)	3913.084053
25	(4, 3, 1, 0)	3913.409057
26	(2, 2, 1, 0)	3914.786853
27	(3, 2, 0, 0)	3925.433806
28	(2, 2, 0, 0)	3925.786566
29	(2, 3, 0, 0)	3925.879649
30	(2, 4, 0, 0)	3926.311190
31	(3, 3, 0, 0)	3927.427240
32	(4, 2, 0, 0)	3927.427417
33	(3, 4, 0, 0)	3928.376406
34	(4, 3, 0, 0)	3929.059361
35	(4, 4, 0, 0)	3930.078725
'''
```
```# set the parameters that give the lowest AIC
p, q, P, Q = result_table.parameters

seasonal_order=(P, D, Q, s)).fit(disp=-1)
print(best_model.summary())
'''
Statespace Model Results
==========================================================================================
Dep. Variable:                                Ads   No. Observations:                  216
Model:             SARIMAX(2, 1, 3)x(1, 1, 1, 24)   Log Likelihood               -1936.321
Date:                            Sat, 07 Mar 2020   AIC                           3888.642
Time:                                    15:00:14   BIC                           3914.660
Sample:                                09-13-2017   HQIC                          3899.181
- 09-21-2017
Covariance Type:                              opg
==============================================================================
coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1          0.7913      0.270      2.928      0.003       0.262       1.321
ar.L2         -0.5503      0.306     -1.799      0.072      -1.150       0.049
ma.L1         -0.7316      0.262     -2.793      0.005      -1.245      -0.218
ma.L2          0.5651      0.282      2.005      0.045       0.013       1.118
ma.L3         -0.1811      0.092     -1.964      0.049      -0.362      -0.000
ar.S.L24       0.3312      0.076      4.351      0.000       0.182       0.480
ma.S.L24      -0.7635      0.104     -7.361      0.000      -0.967      -0.560
sigma2      4.574e+07   5.61e-09   8.15e+15      0.000    4.57e+07    4.57e+07
===================================================================================
Ljung-Box (Q):                       43.70   Jarque-Bera (JB):                10.56
Prob(Q):                              0.32   Prob(JB):                         0.01
Heteroskedasticity (H):               0.65   Skew:                            -0.28
Prob(H) (two-sided):                  0.09   Kurtosis:                         4.00
===================================================================================

Warnings:
 Covariance matrix calculated using the outer product of gradients (complex-step).
 Covariance matrix is singular or near-singular, with condition number 1.08e+32. Standard errors may be unstable.
'''
```
```tsplot(best_model.resid[24+1:], lags=60)
``` ```def plotSARIMA(series, model, n_steps):
"""Plots model vs predicted values

series - dataset with timeseries
model - fitted SARIMA model
n_steps - number of steps to predict in the future
"""

data = series.copy()
data.columns = ['actual']
data['sarima_model'] = model.fittedvalues
# making a shift on s+d steps, because these values were unobserved by the model
# due to the differentiating
data['sarima_model'][:s+d] = np.NaN

# forecasting on n_steps forward
forecast = model.predict(start = data.shape, end = data.shape+n_steps)
forecast = data.sarima_model.append(forecast)
# calculate error, again having shifted on s+d steps from the beginning
error = mean_absolute_percentage_error(data['actual'][s+d:], data['sarima_model'][s+d:])

plt.figure(figsize=(15, 7))
plt.title("Mean Absolute Percentage Error: {0:.2f}%".format(error))
plt.plot(forecast, color='r', label="model")
plt.axvspan(data.index[-1], forecast.index[-1], alpha=0.5, color='lightgrey')
plt.plot(data.actual, label="actual")
plt.legend()
plt.grid(True)
```
```plotSARIMA(ads, best_model, 50)
```   290 original articles published, 479 praised, 430000 visitors+