ARMA, ARIMA, SARIMA¶

Overview¶

• In the previous chapters, we covered stationarity, smoothing, trend, seasonality, and autocorrelation, and built two different models:

MA models: The current value of the series depends linearly on the series' mean and a set of prior (observed) white noise error terms.

AR models: The current value of the series depends linearly on its own previous values and on a stochastic term (an imperfectly predictable term).

• In this chapter we will review these concepts and combine the AR and MA models into three more complicated ones.

In particular, we will cover:

1. Autoregressive Moving Average (ARMA) models.
2. Autoregressive Integrated Moving Average (ARIMA) models.
3. SARIMA models (ARIMA model for data with seasonality).
4. Selecting the best model.

ARMA¶

The ARMA model (also known as the Box-Jenkins approach) combines two models:

• An autoregressive (AR) model of order $p$.
• A moving average (MA) model of order $q$.

• When we have autocorrelation between outcomes and their ancestors, there will be a pattern in the time series.
• This relationship can be modeled using an ARMA model.
• It allows us to predict the future with a confidence level proportional to the strength of the relationship and the proximity to known values (prediction weakens the further out we go).

ℹ️ Note

• ARMA models assume the time series is stationary.
• A good rule of thumb is to have at least 100 observations when fitting an ARMA model.

Example data¶

• In the following, we'll look at the monthly average temperatures between 1907-1972.

• This is how the data looks like.

• These are some statistics.

• This is the run sequence plot.

• We compute the annual mean and plot it on top of the data.

• This gives us an indication that the mean is rather constant over the years.
• We can extract further information abouth the underlying trend and seasonality by performing a seasonal decomposition.
• We can use both the seasonal_decompose and the STL methods.

• The seasonality is well defined.
• There doesn't seem to be a strong, time-varying trend in the data.
  • We can assume the trend is almost constant.

ARMA modeling stages¶

There are three stages in building an ARMA model:

1. Model identification.
2. Model estimation.
3. Model evaluation.

Model identification¶

• Model identification consists in finding the orders $p$ and $q$ of AR and MA components.
• Before performing model identification we need to:
  1. Determine if the time series is stationary.
  2. Determine if the time series has seasonal component.

Determine stationarity¶

• We will use tools we already know (ADF test).
• We can also look at the rolling mean and std.

⚠ Attention!

• Before we continue, let's consider the following result:

• Periodic signals, by their nature, have means and variances that repeat over the period of the cycle.
• This implies that their statistical properties are functions of time within each period.
• For instance, the mean of a periodic signal over one cycle may be constant.
• However, when considering any point in time relative to the cycle, the instantaneous mean of the signal can vary.
• Similarly, the variance can fluctuate within the cycle.

• The ADF test specifically looks for a unit root (more on this later on).
• A unit root indicates that shocks to the time series have a permanent effect, causing drifts in the level of the series.
• A sinusoidal function, by contrast, is inherently mean-reverting within its cycles.
• After a peak a sinusoid reverts to its mean and any "shock" in terms of phase shift or amplitude change does not alter its oscillatory nature.

• It's crucial to note that the ADF test's conclusion of stationarity for a sinusoid does not imply that the sinusoid is stationary.
• The test's conclusion is about the absence of a unit root.
• This does not imply that the mean and variance are constant within the periodic fluctuations.

• The $p$-value indicates that the time series is stationary...
• ... even if it clearly has a periodic component.
• The rolling mean and rolling standard deviation seem globally constant along the time series...
• ... even if they change locally within the period.

Determine seasonality¶

We can determine if seasonality is present by using the following tools:

• Autocorrelation plot.
• Seasonal subseries plot (month plot).
• Fourier Transform.

Let's first look how these plots look like on synthetic data

• Let's look at the real data now.

• Notice that a violinplot can give a very similar information to the month_plot.

• Finally, to obtain the numerical value of the main periodicity we can use the Foruier Transfrom.
• Here, we'll use the function we defined in the first lecture.

Remove the main seasonality¶

• In this case, it is clear that the main seasonality is $L=12$.
• We can remove it with a seasonal differencing.

⚙ Try it yourself

• Try redoing the previous plots on the differenced data!

Identifying $p$ and $q$¶

As we learned in the previous chapter, we will identify the AR order $p$ and the MA order $q$ with:

• Autocorrelation function (ACF) plot.
• Partial autocorrelation function (PACF) plot.

AR($p$)

• The order of the AR model is identified as follows:
  • Plot 95% confidence interval on the PACF (done as default by statsmodels).
  • Choose lag $p$ such that the partial autocorrelation becomes insignificant for $p+1$ and beyond.

• If a process depends on previous values of itself then it is an AR process.
• If it depends on previous errors than it is an MA process.
• An AR process propagates shocks infinitely.
• AR processes will exhibit exponential decay in ACF and a cut-off in PACF.

• It looks like the PACF becomes zero at lag 2.
• However there is a non-zero partial autocorrelation at lag 3.
• The optimal value might be $p=1$, $p=2$, or $p=3$.
• Note that there are high partial autocorrelations at higher lags, especially 12.
  • This is an effect from seasonality and seasonal differencing.
  • It should not be accounted for when choosing $p$.

MA($q$)

• The order of the MA model is identified as follows:
  • Plot 95% confidence interval on the ACF (default in statsmodels).
  • Choose lag $q$ such that ACF becomes statistically zero for $q+1$ and beyond.

• MA models do not propagate shocks infinitely; they die after $q$ lags.
• MA processes will exhibit exponential decay in PACF and a cut-off in ACF.

• Also in this case there are non-zero autocorrelations at lags 1 and 3.
• So, the values to try are $q=1$, $q=2$, or $q=3$.

Model estimation¶

• Once the orders $p$ and $q$ are identified is necessary to estimate parameters $\phi_1, \dots, \phi_p$ of the AR part and parameters $\theta_1, \dots, \theta_q$ of the MA part.
• Estimating the parameters of an ARMA model is a complicated, nonlinear problem.
• Nonlinear least squares and maximum likelihood estimation are common approaches.
• Software packages will fit the ARMA model for us.

• We split the data in two parts:
  • the training set, that will be used to fit the model's parameters.
  • the test set, that will be used later on to evaluate the prediction performance of the model on unseen data.

ARMA Model Validation¶

• How do we know if our ARMA model is good enough?
• We can check the residuals, i.e., what the model was not able to fit.
• The residuals should approximate a Gaussian distribution (aka white noise).
• Otherwise, we might need to select a better model.

🤔 How to test if the residuals look like noise?

• We will use both visual inspection and statistical tests.
• Visual inspection:
  • ACF plot.
  • Histogram.
  • QQ plot.
• Statistical tests:
  • Normality.
  • Autocorrelation.
  • Heteroskedasticity.

Visual inspection¶

ACF plot

• Checks for any autocorrelation in the residuals.
• White noise should show no significant autocorrelation at all lags.

Histogram and QQ-Plot

• Assess the normality of the residuals.
• White noise should ideally follow a normal distribution.

• The plots are conveniently summarized in the function plot_diagnostics() that can be called on the fit model.

Statistical tests¶

Normality: Jarque-Bera and Shapiro-Wilk tests

$H_0$: the residuals are normally distributed.

• The small $p$-values allow us to reject $H_0$.
• Conclusion: the residuals are not normally distributed.
• Is this a contraddiction?

📝 Note

• For reference, let's see what these tests say about data that are actually normally distributed.
• Try executing the cell below multiple times and see how much the results changes each time.
• These tests start to be reliable only for large sample sizes ($N>5000$).

Autocorrelation: Ljung-Box test

$H_0$: the residuals are independently distributed (no autocorrelation).

• There is a $p$-value for each lag.
• Here we just take the mean, but one might also want to look at the at largest lag (pval[-1]).
• It is also not always obvious to select how many lags should be used in the test...

Autocorrelation: Durbin Watson test

• Tests autocorrelation in the residuals.
• We want something between 1-3.
• 2 is ideal (no serial correlation).

Heteroskedasticity test

• Tests for change in variance between residuals.

$H_0$: no heteroskedasticity.

• $H_0$ indicates different things based on the alternative $H_A$:
  • $H_A$: Increasing, $H_0$: the variance is not increasing throughout the series.
  • $H_A$: Decreasing, $H_0$: the variance is not decreasing throughout the series.
  • $H_A$: Two-sided (default), $H_0$: the variance does not increase nor decrease throughout the series.

Summary of our tests

Independence:

• ✅ ACF plot.
• ✅ Ljung-Box test.
• ✅ Durbin Watson test.

Normality:

• ✅ Histogram/Density plot.
• 🤔 QQ-plot
• ❌ Jarque-Bera (reliable for large sample size).
• ❌ Shapiro-Wilk (reliable for large sample size).

Heteroskedasticity

• ❌ Heteroskedasticity test.

• The tests are a bit inconclusive.
• There is no strong evidence that the model is either very good or very bad.
• It is probably wise to try other candidate models, e.g., ARMA(2,0,2), and repeat the tests.

ARMA Model Predictions¶

• Once the model is fit, we can use it to predict the test data.
• The predictions come in the form of a distribution.
• In other words, ARMA performs a probabilistic forecasting.
• The mean (mode) of this distribution correspond to the most likely value and correspond to our forecast.
• The rest of the distribution can be used to compute confidence intervals.

ARIMA Model¶

• ARIMA stands for Auto Regressive Integrated Moving Average.
• ARIMA models have three components:
  • AR model.
  • Integrated component (more on this shortly).
  • MA model.

• The ARIMA model is denoted ARIMA($p, d, q$).
  • $p$ is the order of the AR model.
  • $d$ is the number of times to difference the data.
  • $q$ is the order of the MA model.
  • $p$, $d$, and $q$ are nonnegative integers.

• As we saw previously, differencing nonstationary time series data one or more times can make it stationary.
• That’s the role of the integrated (I) component of ARIMA.
• $d$ is the number of times to perform a lag 1 difference on the data.
  • $d=0$: no differencing.
  • $d=1$: difference once.
  • $d=2$: difference twice.

• The ARMA model is suitable for stationary time series where the mean and variance do not change over time.
• The ARIMA model effectively models non-stationary time series by differencing the data.
• In practice, ARIMA makes the time series stationary before applying the ARMA model.
• Let's see it with an example.

SARIMA¶

• To apply ARMA and ARIMA, we must remove the seasonal component.
• After computing the predictions we had to put the seasonal component back.
• It would be convenient to directly work on data with seasonality.

• SARIMA is an extension of ARIMA that includes seasonal terms.
• The model is specified as SARIMA $(p, d, q) \times (P, D, Q, s)$:
  • Regular ARIMA components $(p, d, q)$.
  • Seasonal components $(P, D, Q, s)$ where:
    • $P$: Seasonal autoregressive order.
    • $D$: Seasonal differencing order.
    • $Q$: Seasonal moving average order.
    • $s$: Number of time steps for a single seasonal period.

How to select the values $s, P, D, Q$?

• $s$:
  • Is the main seasonality in the data.
  • We already know how to find it.

• $P$ and $Q$:
  • A spike at $s$-th lag (and potentially multiples of $s$) should be present in the PACF/ACF plots.
  • For example, if $s = 12$, there could be spikes at $(s*n)^{th}$ lags too.
  • Pick out the lags with largest spikes as candidates for $P$ or $Q$.

• $D$:
  • Is the number of seasonal differencing required to make the time series stationary.
  • Is often determined by trial and error or by examining the seasonally differenced data.

💡 Rule of thumb

• Before selecting $P$ and $Q$, ensure that the series is seasonally stationary by applying seasonal differencing if needed ($D$).

• Look the ACF plot to identify the seasonal moving average order $Q$.

  • Look for significant autocorrelations at seasonal lags (multiples of $s$).
  • If the ACF plot shows a sharp cut-off in the seasonal lags, it suggests the order of the seasonal MA component ($Q$).

• Look at the PACF plot to identify the seasonal autoregressive order $P$.

  • Look for significant spikes at multiples of the seasonality $s$.
  • A sharp cut-off in the PACF at a seasonal lag suggests the order of the seasonal AR terms ($P$) needed.

AutoARIMA¶

• At this point it should be clear that identifying the optimal SARIMA model is difficult.
• It requires careful analysis, trial and errors, and some experience.
• The following cheatsheet summarizes some rules of thumb to select the model.

💡 Cheatsheet for coefficients setting

• An alternative to manual model selection is to use automated procedures.
• Here enters AutoARIMA.
• AutoARIMA requires you to specify the maximum range of values to try.
• Afterwards, it tries to find the best confgiuration among the possible ones.
• See here for a complete list and description of the options available.

AutoARIMA or not AutoARIMA?¶

While being very convenient, like all automated procedures auto_arima comes with drawbacks.

1. auto_arima can be computationally expensive, especially for large datasets and when exploring a wide range of models.

2. Automated model selection lacks the qualitative insights a human might bring to the modeling process.
   • These, include understanding business cycles, external factors, or anomalies in the data.

3. The defaults in auto_arima may not be optimal for all time series data.
   • The range of values to explore should be adjusted properly each time.
   • This might become almost as tricky as doing manual model selection.

4. auto_arima requires a sufficiently long time series to accurately identify patterns and seasonality.

5. The selection of the best model is typically based on statistical criteria such as AIC or BIC.
   • These, might not always align with practical performance metrics such as MSE.

Grid search¶

• Selecting the best ARIMA model is as much an art as it is a science, involving iteration and refinement.
• A common approach is to select a set of candidates for $(p,d,q)\times (P,D,Q,s)$ and fit a model for each possible combination.

• For each fit model:
  • Analyze the residuals using the visual techniques and the statistical tests discussed before.
  • Evaluate the prediction performance. 🤔 How?
  • Evaluate its complexity. 🤔 How?

Prediction performance¶

Use MSE and/or MAPE to evaluate the prediction performance of the model.

Mean Squared Error (MSE)

• MSE is the average of the squared differences between the observed values and the predictions.
• $MSE = \frac{1}{n}\sum_{i=1}^{n}(Y_i - \hat{Y}_i)^2$, where $Y_i$ is the actual observation and $\hat{Y}_i$ is the forecasted value.
• In SARIMA model selection, a model with a lower MSE is preferred, indicating better fit to the data.

Mean Absolute Percentage Error (MAPE)

• MAPE is the average of the absolute percentage errors of forecasts.
• $MAPE = \frac{100\%}{n}\sum_{i=1}^{n}\left|\frac{Y_i - \hat{Y}_i}{Y_i}\right|$, where $Y_i$ is the actual observation and $\hat{Y}_i$ is the forecasted value.
• MAPE expresses errors as a percentage, making it straightforward to understand the magnitude of forecasting errors.
• If you are comparing models that predict different quantities (e.g., dollars vs. units sold), the percentage error allows for a more apples-to-apples comparison.
• Also, when you're more interested in the relative size of the errors than in their absolute size, MAPE is relevant.
• Finally, MAPE is useful when the magnitude of the data varies significantly.

Model complexity¶

Use AIC or BIC to estimate the model's complexity.

Akaike Information Criterion (AIC)

• AIC is a measure of the relative quality of statistical models for a given set of data.
• It deals with the trade-off between the goodness of fit of the model and the complexity of the model.
• $AIC = 2k - 2\ln(\hat{L})$, where $k$ is the number of parameters in the model, and $\hat{L}$ is the maximized value of the likelihood function for the model.
• The model with the lowest AIC value is preferred, as it fits the data well but is not overly complex.

Bayesian Information Criterion (BIC)

• Similar to AIC, the BIC is another criterion for model selection, but it introduces a stronger penalty for models with more parameters.
• $BIC = \ln(n)k - 2\ln(\hat{L})$, where $n$ is the number of observations, $k$ is the number of parameters, and $\hat{L}$ is the maximized likelihood.
• A lower BIC value indicates a better model, preferring simpler models to complex ones, especially as the sample size $n$ increases.

Restricting the search with Exploratory Data Analysis (EDA)¶

• Grid search can be very expensive if done exaustively, especially on limited hardware.
• An Exploratory Data Analysis can help to significantly reduce the number of candidates to try out.

Selecting the candidates for differentiation¶

• Let's start by identifying all the candidates for seasonal and general differencing.
• In this case, we already know that the main seasionality is $s=12$.
• Should we apply first the general or the seasonal differencing?

• If seasonal patterns are dominant and the goal is to remove seasonality before addressing any trend, start with seasonal differencing.
  • This is particularly useful when the seasonal pattern is strong and clear.
• If the trend is the predominant feature, you might start with standard differencing.

• In our case, we go with seasonal differencing first.

Filter-out non-stationary candidates¶

• Among all the differenced time series, keep only those that are stationary (according to ADF).

Select candidates for $p, q, P, Q$¶

• We are going to make a script to extract candidates for the orders of the AR and MA component, both in the general and the seasonal part.
• We will leverage the ACF and PACF functions.
• So far, we looked at the acf_plot and pacf_plot.
• Now we need to use acf and pacf.
• We need to understand how these relate.

Define the grid of values to search¶

Train the models¶

• We defined a function that takes every model configuration and trains a model.
• For each model, we save the MSE, MAPE, AIC and BIC.

Analyze the results¶

• Show the 10 best models according to the performance (MSE, MAPE) and model complexity (AIC, BIC).

• We can check if there are overlaps between the 4 different groups using the set function.

• Show the top model according to each metric.

Compute performance on the test set¶

• Take the best models, compute the predictions and evaluate their performance in terms of MAPE on the w.r.t. test data.

Summary¶

1. Autoregressive Moving Average (ARMA) models.
2. Autoregressive Integrated Moving Average (ARIMA) models.
3. SARIMA models (ARIMA model for data with seasonality).
4. Selecting the best model.

Exercise¶

• Look at sensor data that tracks atmospheric CO2 from continuous air samples at Mauna Loa Observatory in Hawaii. This data includes CO2 samples from MAR 1958 to DEC 1980.

• Determine the presence of main trend and seasonality in the data.
• Determine if the data are stationary.
• Split the data in train (90%) and test (10%).
• Find a set of SARIMAX candidate models by looking at the ACF and PACF.
• Perform a grid search on the model candidates.
• Select the best models, based on performance metrics, model complexity, and normality of the residuals.
• Compare the best model you found with the one from autoarima.