Time series Concepts and modelling

In this notebook, we shall study and experiment with the statistical tools used for time series data. We will begin with examining characteristics of a time series followed by modelling it. As an example, we take "node_memory_Active_bytes" which are the memory bytes that have been recently used by a node. It reflects the memory pressure on the node. This notebook can be used as a guideline to analyze, understand, and model any time series data.

Contents

• Load data

• Concepts

  • Standard time series

    • White noise
    • Random walk

  • Stationarity

    • Dicky Fuller Test

  • Seasonality

    • Naive decomposition - Additive
    • Multiplicative
    • Synthetic seasonal data

  • Auto Correlation

    • ACF
    • PACF

• Modelling

  • Manual

    • Linear Regression
    • Exponential Smoothing
    • AR
    • MA
    • ARMA
    • ARIMA
    • SARIMAX

  • Automatic

    • Auto ARIMA

Load Data

Concepts

Standard time series

Before defining statistical properties, let's understand two standard time series: white noise and random walk. These time series are often used to understand and compare others.

White noise

• White noise has constant mean, constant variance, and no autocorrelation (no relationship between past and present of time series)

Random walk

• A random walk is modelled using this equation: Pt = Pt-1 + $ϵ$t where $ϵ$t is white noise.
• Random walks are not predictable

Stationarity

Intuitively, a stationary time series does not change the way it changes. In other words, the statistical properties of the process generating the time series does not change over time. For real world problems, we use the weak form or covariance stationarity definition:

• Constant $\mu$
• Constant $\sigma$
• Cov (xn, xn+k) = Cov (xm, xm+k)

Ideally we would want the time series to be stationary for modelling. In the real world it is rarely stationary. We shall now look at the Dickey Fuller Test which determines if a series is stationary or not.

Dickey Fuller Test

• This test determines if the time series is stationary or not.
• If p value > 0.05 then process is not stationary.
• If p < 0.05, null hypothesis is rejected and the process is stationary.

For this series, the p-value is 0 which means that the null hypothesis can be rejected and the series is stationary. The test statistic of -15.6815 is less than 1% critical value (-3.4304) meaning that the chances of this result being a fluke is less than 1%.

As we can see from the above, we know by definition that that white noise is stationary and the random walk is not. Since the results of the Dickey Fuller test verified our expected results, we know that we can at least have some confidence in this statistical test when it indicates to us that our real world time series data set is stationary.

Seasonality

• There are certain trends in some time series that appear cyclically.
• We can decompose the time series data into trend, seasonal, residual
• Trend: general pattern present in the time series
• Seasonal: cyclical effect in the time series
• Residual: Error of prediction

Naive Decomposition

One way to decompose a time series in order to check for seasonality is through naive decomposition that assumes an additive or multiplicative nature between the three components: trend, seasonal, residual.

Additive

• Observed = trend + seasonal + residual

Multiplicative

• Observed = trend seasonal residual

Both additive and multiplicative models give a similar output. There is a downward trend in memory usage at the begining of the week and it remains low throughout the week increasing slightly towards the end. The seasonal values show an oscillating behavior, however we can see that within a day, it starts high and reduces towards the end. The residuals however do not seem exactly random (like white noise) indicating that the decomposition is not perfect. For a better understanding, let's try these methods on a synthetic seasonal time series.

No surprises here, but we see that the time series being perfectly decomposed into the trend and the seasonal variations. The residuals look like white noise as well. We can use such decomposition for simpifying and understanding all time series data.

Auto-correlation

Auto correlation refers to correlation of a time series with a lagged version of itself.

ACF

Auto correlation function for how the number of lags we are interested in.

• X axis represents the number of lags.
• Y axis represents the correlation value of the time series with a lagged version of itself
• Blue area depicts significance. In this plot, all lines are higher than significance which means that all the autocorrelation coefficients are significant.
• Interpretation: There is time dependence in the data since we see high correlation. In other words, we can even use 50 points behind to predict what comes next.
• This curve is also indicative of some non-stationarity in the dataset. Although the dickey-fuller test shows stationarity, ACF plot with many significant coefficients may require integrations before modelling (explained in ARIMA section).

• We see that most coefficients are not significant
• Randomly, some of them cross the blue area
• Hence, white noise has no auto correlation or no time dependence

PACF (Partial Auto correlation)

• Accounts for direct effect of xlagged on xt
• For example, ACF for a lag of 3, will also include effect of xt-3 on xt-2, xt-2 on xt-1, and xt-1 on xt
• PACF for a lag of 3, only accounts for effect of xt-3 on xt
• The first value is same for both as there are no indirect influences in the first computation.

Great! now we have a basic understanding of statistical properties of the time series data. Next, we'll look at how to model a time series data.

Modelling

Before we dive into various time series models, here are some general points to keep in mind:

• Start with a simple parsimonious model and add more parameters and complexity as required.
• If coefficients of added variables are significantly different from zero keep them otherwise prefer simpler models to favour generalizability.
• To see if we should select the more complex model, we can compare log likelihood, information criterion like AIC, and BIC.
• Residuals: If the model fits the data well then the residuals should resemble white noise (or in other words, there is no trend left to capture).

Linear Regression Model

Linear regression is a linear model, a model that assumes a linear relationship between the input variables (x) and the single output variable (y). y can be calculated from a linear combination of the input variables (x). If we try to draw a relatonshop between these two variables, we get a straight line that can be mathematically expressed as :

• y = mx + b

Where b is the intercept and m is the slope of the line.

• yt = β0 + β1xt + εt

The coefficients β0 and β1 denote the intercept and the slope of the line respectively. The intercept β0 represents the predicted value of y when x=0. The slope β1 represents the average predicted change in y resulting from a one unit increase in x.

So basically, the linear regression algorithm gives us the most optimal value for the intercept and the slope in a two dimensional environment. The x and y values

As discussed previously, we saw that Linear Regression model looks for the best value for the intercept and result in a line that best fits the data. Let's see how we can retrieve those parameters.

This means that for every one unit of change in time, the change in the value is about -4126.27. However, this information doesn't help a lot in knowing how the time series data would change with time. This method can be mostly be useful in cases where the data is linear. In this scenario we have some trends and seasonality as well. We will now move towards Exponential Smoothing models to see how they are different and better than linear models.

Exponential Smoothing Model

Forecasts produced using exponential smoothing methods are weighted averages of past observations, with the weights decaying exponentially as the observations get older. In other words, the more recent the observation the higher the associated weight.

Simple Exponential Smoothing

Simple exponential smoothing method is most suitable for forecasting data with no clear trend or seasonal pattern. It requires a single parameter, called alpha (α), also called the smoothing factor or smoothing coefficient.

• Alpha (smoothing_level) controls the rate at which the influence of the observations at prior time steps decay exponentially.
• Alpha is often set to a value between 0 and 1.
• Large values mean that the model pays attention mainly to the most recent past observations, whereas smaller values mean more of the history is taken into account when making a prediction.

• As we mentioned, that SES has only one component, let's see how it is represented in mathematical form :

  • Forecast : pt = lt
  • Level : lt = α yt + (1-α)lt-1

  Where p is the forecast, l is the level and α is the smoothing factor.

The above model is a result of the autofit model by statsmodels. As mentioned above, the depends on one parameter alpha also known as the smoothing factor or smoothing coefficient. According to autofit, the best alpha for this data would be 0.995. That means the model pays attention mainly to the most recent past observations.

Optimize Alpha

We could also try other values for alpha as follows and will chose the minimum AIC i.e. Akaike’s Information Criterion which is defined as :

AIC = -2log(L) + 2k

where L is the likelihood of the model and k is the total number of parameters and initial states that have been estimated (including the residual variance).

The measures on the training set (training sample) are not really suitable as basis for model selection. It is because in the training sample it is always possible to overfit, and the richer the model, the better the fit will be. Meanwhile, information criteria like AIC or BIC take this into account and penalize for the model complexity accordingly. Therefore, they generally are suitable for model selection. For this model let's try to minimize the value of AIC in order to find the best model.

The above implementation of Simple Exponential Smoothing allows us to compare the autofit and manual grid search for best alpha and best model. Both the alphas are pretty comparable and gives us the best model. However, this model performs better where there is no seasonality and trend. So we will move to Holt's model that covers the trend in data but no seasonality.

Holt's Double Exponential Smoothing Model

Holt's Double Exponential Smoothing method is similar to Simple Exponential Smoothing. It calculates the level component to measure the level in the Forecast.In addition to the alpha parameter for controlling smoothing factor for the level, an additional smoothing factor is added to control the decay of the influence of the change in trend called beta (β).

• Forecast Equation : pt+h|t = lt + hbt

• Level Equation : lt = αyt + (1-α)lt-1

• Trend Equation : bt = β(lt - lt-1) + (1-β)bt-1

  Where p is the forecast, l is the level, α is the smoothing factor and β is the smoothing slope.

Just like we did for Simple Exponential Smoothing model, for this model also let's try to minimize the value of AIC in order to find the best model. The only difference here would be that we have two parameters : alpha and beta.

The above implementation of Simple Exponential Smoothing allows us to compare the autofit and manual grid search for best alpha, beta and best model. Both the alphas and betas are very different from each other. However, this model performs better where there is no seasonality. So we will move to other approaches that capture seasonality as well.

Auto regressive (AR) Model

The next model we will explore is the AR model. Linear regression models the linear relationship of predictors whereas ARIMA models the linear relationship of past observation. This makes ARIMA naturally adaptive to new data [source]. The following equation shows AR(1) or auto regression model that learns from just one past observation.

xt = c + $\varphi$xt-1 + $ϵ$t

• $\varphi$ is a numerical constant from -1 to 1 (otherwise it'll blow up)
• $ϵ$t represents unpredictable shocks.

In real life, looking at just 1 lag may not give reasonable results. To determine the optimal number of lags, we look at PACF plots and incrementally increase the number. If we have coefficients close to 1, that means they can be good candidates for the model. Note: AR model stationary process better. For non-stationary time series, it is often converted to stationary process first.

Interpretation:

• ar.L1.value is the value of $\varphi$ in the AR(1) equation.
• const is the value of c in AR(1) equation.
• p < 0.05 suggests that these coefficients are significant.
• Since the ar.L1.value is not 0, this variable adds to the model.
• Great! let's increase the order and see if the more complex model is better.

Alright, we have a more complex model, but is it adding anything new? The log likelihood ratio (LLR) test determines if the new model adds value to the modelling significantly. It returns a p value, if p<0.05 then the complex model is better.

So the second model is significantly better than the first.

The above code suggests that going from order 3 to order 4 does not add any significant value in the model. We can also see that by looking at the summaries of the results. Side note: These are rules of thumb but as we see in the pacf plot, some higher lags than 3 are also significant, it could be that a higher lag model may fit the series better.

At the end the more complex model should have both insignificant $\varphi$ values and insignificant log value. When we reach that stage, in this case at the lag 3, we stop.

Side note: If you want to compare different time series, make sure to normalize values. One way to do this is to divide the series by the first number.

Analyzing the residuals

• This is a very important step in analysing the model performance.
• Remember that residuals should be similar to white noise reflecting that there is no trend remaining to be captured.
• We will now plot the residuals, check for stationarity, and check ACF plots.

Interpretation

• The plot shows that the residuals look random but are not exactly like white noise.
• The Dickey Fuller test suggests that the residuals are stationary.
• The ACF plots show coefficients that are significant further indicating that the residuals are not perfectly white noise and there could be better predictors.

Moving average model

• Better suited if the data has unexpected shocks
• Accounts for the past residuals
• Absorbs shocks (accounts from past errors) and correct

rt = c + $\theta$$ϵ$t-1 + $ϵ$t

• For finding the optimal number of lags for this model, we rely more on ACF than PACF because it learns from residuals, so the direct effect of a past observation on the present day is not relevant.
• Mathematical fact: If MA coefficients are close to 1 then it is basically an approximation of AR model.
• Both AR, MA models are bad with non-stationary data.

Most the coefficients in the lags are significant implying that MA(n+1) will always be better than MA(n). This type of data suggests that some type of a combination of AR and MA is required to correctly model the data. For reference, let's train a MA model of order 1.

model_ma_1 = ARMA(ts, order=(0, 1)) ## p, q where p is AR and q is MA param

model_ma_results_1 = ARMA(ts, order=(0, 1)).fit() model_ma_results_1.summary()

The results show that the coefficient is significant.

ARMA (AR and MA models together)

• The ARMA model considers both lagged value and lagged errors.
• Following equation represents ARMA(1,1):

rt = c + $\varphi$rt-1 + $\theta$$ϵ$t-1 + $ϵ$t

• For ARMA model, we select complex model and move to simpler models.
• For determining p, q take hint from PACF and ACF respectively.
• The coefficients of the complex model should be significantly better than 0 and the either the AIC should be lower or the LLR test should give significant p-values.
• LLR test only works for nested models (one model should have all the lag terms of the other). While comparing ARMA(1,3) and ARMA(3,1), we have to look at AIC. Lower values of AIC are prefered over higher values.

We found 3 lags for AR component and MA can be anything.

ARMA(4,3) and ARMA(3,3) have low AIC, with all significant coefficients. Both of these models are good candidates for fitting the series.

Analyzing the residuals

• This plot has are several significant values in the PACF.
• We can train higher order ARMA models but usually we have to go for ARIMA models.

ARIMA models

• Auto-regressive integrated moving average model.
• Since ARMA models do not work well with non-stationary data, we introduce integrations to make the data the stationary.
• The model builds on the difference between the time periods (integerated values) instead of the values itself.
• For example, converting prices to returns in stock market data would be similar to 1 degree of integration.
• Intuitively, modelling returns instead of prices may be better since returns may resemble stationary data.
• ARIMA has p, d, q hyperparameters.
• d -> number of degree of changes to ensure stationarity.
• One way of finding d is to compute delta of the data unitll stationarity is found.
• ARIMA method uses AIC criteria to compare models with different orders.
• The following equation represents ARIMA(1,1,1).
• Start from a simple model and increase orders by looking at residuals.

$\Delta$Pt = c + $\varphi$$\Delta$Pt-1 + $\theta$$ϵ$t-1 + $ϵ$t

We see that the third lag is significant. We can try all models from ARIMA(0,0,0) to ARIMA(3,1,3).

Great, we see that the first 4 coefficients are now not significant. At this point, as an additional exercise, we could try models with lag 5,6, and 7 as well to further improve performance. Next, we are going to look at how to select the value for "d", or the number of integrations required for modelling the data. Remember that integrations are done to introduce stationarity and should be avoided if the data is already stationary. In this case, the adfuller test says that the data is stationary but ACF plot shows many significant correlations. Hence we try d as 1. Let's manually compute 1 degree of integration and check stationarity.

Cool, the ACF plot has less significant correlations and the adfuller test also indicates that after one integration the data is stationarity. For exercise, and a rigorous search for the best model, we could check for models with d=2 as well. However, we prefer simpler models over complex models if the improvements are not significant. Let's fit ARIMA(3,1,1) and evaluate the results.

Cool, we see 9 insignificant coefficients indicating that this model is a good choice.

ARIMAX (Exogenous factors)

In addition to past observations, sometimes external factors also affect the future value of a time series

$\Delta$Pt = c + $\beta$X + $\varphi$$\Delta$Pt-1 + $\theta$$ϵ$t-1 + $ϵ$t

• X can be anything (data should be available for each period).
• E.g. S&P prices for predicting a particular stock price.
• $\beta$ is a fitted coefficient.

SARIMAX (Seasonality)

Some time series have a seasonal aspect to them. For example, consider the number of views for a christmas song every month. Every year in December we can expect a spike in the series.

• SARIMAX model includes a lag term for the last period.
• Hyper parameters of SARIMAX: (p,d,q)(P,D,Q,s).
• (p,d,q) is the ARIMA order.
• 's' signifies how far behind do we observe the seasonal affect; s=1 means no seasonality.
• In our number of songs example, the value of s will be 12.
• P, Q specifies how many seasonal AR and MA terms we will have so for s=12, P=2, Q=1, the equation will have yt-12, yt-24, and $ϵ$t-12
• D specifies seasonal integrations.
• The total number of coefficients for the model will be: P + Q + p + q coeffs.
• In our current example we don't see any seasonality but let's train a SARIMAX model to understand how it works.

Automated

• In the final section of this notebook, we show the use of auto_arima function of pmdarima library that automatically finds hyperparameters for a SARIMAX model.
• However, we need to proceed with caution and do our due diligence before selecting a new model as automated values may choose a complicated model but in real life we may prefer a simpler one.
• The function does dickey fuller test to find d; and uses one information criteria to select between models (default being "aic").
• The function does not look at whether all coefficients are significant or not, we have to consider that on our own.

• Here we see the auto_arima selects an ARIMA model with order 2,1,3. In our manual analysis using statsmodel package, the model with 2,1,3 could not converge. We can make changes to the default parameters to the function to give it more information.
• An important information criteria is oob which selects model based on the score performance of the model on a validation set. Let's try that.

With updated parameters, especially the information criterion as "oob" we see the function selects an inferior model in terms of AIC and significance of the coefficients. However this model performs well on the validation set. We need to conduct further experiments with a test set to find the generalizability of these models.

Conclusion

Awesome! Apologies for a really long notebook, but we have learnt how to understand characteristics of a time series and then how to model it. Next, we are going to look at how we can forecast using these models.

Referrences

• A lot of the content in this notebook has been adapted from this youtube playlist.