A time series is a sequence of observations recorded at successive equally spaced points in time. It is a form of discrete time data. Time series are generally visualised using a line chart. Time series analysis comprises methods for analyzing time series in order to extract meaningful information. Time series forecasting is the use of a model(such as ARIMA) to predict future values based on previously observed values.
A time series is a function T where x is the time variable and T(x) is the value of the time series at a given time x.
There are three components of a time series : trends, seasonality and variance.
Examaples: A time series that is increasing/decreasing has a trend
If the function has some degree of periodicity then we can say that it is seasonal.
Before we can begin with forecasting, we must first learn about some statistical tools. They are: Autocorrelation and Transformations.
Autocorrelation is the correlation of the data with itself. So, instead of measuring the correlation between two random variables, we are measuring the correlation between a random variable against itself. Hence, why it is called auto-correlation.
Correlation is how strongly two variables are related to each other. If the value is 1, the variables are perfectly positively correlated, -1 they are perfectly negatively correlated and 0 there is no correlation.
For time-series, the autocorrelation is the correlation of that time series at two different points in time (also known as lags). In other words, we are measuring the time series against some lagged version of itself.
Where r_k is the autocorrelation at k lags, y bar is the mean, y_t is the value of the variable at time t.
The autocorrelation at k lags is the Pearson coefficient of a time series with itself shifted by k lags.
As stated above, we use autocorrelation to measure the correlation of a time series with a lagged version of itself. This computation allows us to gain some interesting insight into the characteristics of our series:
Seasonality: Let's say we find the correlation at certain lag multiples is in general higher than others. This means we have some seasonal component in our data. For example, if we have daily data and we find that every multiple of 7 lag terms is higher than others, we probably have some weekly seasonality.
Trend: If the correlation for recent lags is higher and slowly decreases as the lags increase, then there is some trend in our data.
Example:
import plotly.express as px
import plotly.graph_objects as go
import pandas as pd
from statsmodels.graphics.tsaplots import plot_acf
import matplotlib.pyplot as plt
import numpy as np
data = pd.read_csv('m.csv', index_col=0)
data.index = pd.to_datetime(data.index)
# Plot the data
fig = px.line(data, x=data.index, y='Sales',
labels=({'Sales': 'Sales', 'Month': 'Date'}))
fig.update_layout(template="simple_white", font=dict(size=18),
title_text='Sales', width=650, title_x=0.5, height=400)
fig.show()
plt.rc("figure", figsize=(8,4))
plot_acf(data['Sales'], lags=48)
plt.ylim(0,1)
plt.xlabel('Lags', fontsize=18)
plt.ylabel('Correlation', fontsize=18)
plt.xticks(fontsize=18)
plt.yticks(fontsize=18)
plt.title('Autocorrelation Plot', fontsize=20)
plt.tight_layout()
plt.show()
The blue shaded region contains points that are statistically irrelevant.
From the above autocorrelation plot we observe that over time the correlation slowly decays. This means that there is a trend in our data. Second, there is a sharp spike in the correlation every 10 to months, implying that there is seasonality in the data. In this case, it is yearly seasonality.
Autocorrelation is an effective first tool in time series analysis in identifying whether a given data set is seasonal or has a trend, if it is not immediately obvious from the graph.
A time series is stationary if it does not exhibit any long term trends or obvious seasonality.
It has:
A major disadvantage of many forecasting models is that they cannot process non stationary data. To make the data stationary, we must apply certain transformations.
A transform is a function that converts a time series into another. This is done for many reasons: to attain a different perspective on the data, or to make it suitable for prediction algorithms.
The purpose of differencing is to remove any trends from the data. This is achieved by plotting the difference between the current and previous function values. Differencing may be done multiple times to achieve trendless data.
from statsmodels.tsa.stattools import adfuller
data["Sales_Diff"] = data["Sales"].diff()
fig2 = px.line(data, x=data.index, y=data["Sales_Diff"], labels={'Sales': "Differenced Sales",'Month' : "Date"})
fig2.update_layout(template="simple_white", font=dict(size=18),
title_text="Differenced Sales", width=650,
title_x=0.5, height=400)
fig2.show()
The data no longer shows any trends or seasonality. Mean seems to be stable, but variance seems to slightly increase with time. Depending on your criteria, you may move forward with a single differencing step or keep going.
The purpose of the logarthim transform is to reduce the amplitude of peaks and troughs, in order to stabilise variance.
data["Sales_Log"] = np.log(data["Sales"])
fig3 = px.line(data, x=data.index, y=data["Sales_Log"], labels={'Sales_Log': "Log Transformed Sales",'Month' : "Date"})
fig3.update_layout(template="simple_white", font=dict(size=18),
title_text="Logarithm Transform", width=650,
title_x=0.5, height=400)
fig3.show()
This equalises the fluctuations to a certain extent, but there is still a trend. Generally this is combined with differencing to achieve the intended result.
data["Sales_Log_Diff"] = data["Sales_Log"].diff()
fig3 = px.line(data, x=data.index, y=data["Sales_Log_Diff"], labels={'Sales_Log_Diff': "Transformed Sales",'Month' : "Date"})
fig3.update_layout(template="simple_white", font=dict(size=18),
title_text="Differenced Logarithm Transform", width=650,
title_x=0.5, height=400)
fig3.show()
The trend has been removed, and the variance has been stabilised to an extent.
The box-cox transformation is much like the logarithm transform, with the same intentions of equalising variance while being slightly more effective.
Where The parameter λ which ranges from -5 to 5 is estimated by seeing which value best transforms the data into a normal distribution. Most packages(like the ones we are using) automatically do this for us.
from scipy.stats import boxcox
data['Sales_box_cox'], lam = boxcox(data['Sales'])
fig4 = px.line(data, x=data.index, y=data["Sales_box_cox"], labels={'Sales_box_cox': "Box-Cox Sales",'Month' : "Date"})
fig4.update_layout(template="simple_white", font=dict(size=18),
title_text="Box-Cox Transform", width=650,
title_x=0.5, height=400)
fig4.show()
Once again , we must difference this to remove the trend.
data['Sales_box_cox_diff'] = data['Sales_box_cox'].diff()
data.dropna(inplace=True)
fig5 = px.line(data, x=data.index, y=data["Sales_box_cox_diff"], labels={'Sales_box_cox_diff': "Differenced Box-Cox Sales",'Month' : "Date"})
fig5.update_layout(template="simple_white", font=dict(size=18),
title_text="Differenced Box-Cox Transform", width=650,
title_x=0.5, height=400)
fig5.show()
The trend is gone, and the variance is moderately stable.
Forecasting is predicting future data points using already available data points. There are many models we can use for this purpose.
Autoregression is when you forecast a time series using some linear weighted combination of the previous values (lags) of that time series. As we are regressing a target value against itself, it is called auto-regression. Mathematically, we can write autoregression as:
Where y is the time series we are forecasting at various time steps, ϕ are the fitted coefficients of the lags for the time series, ε is the error term and p is the number of lagged components included in the model, this is also known as the order.
Where y is the time series we are forecasting at various time steps, ϕ are the fitted coefficients of the lags for the time series, ε is the error term (typically normally distributed) and p is the number of lagged components included in the model, this is also known as the order.
A limitation of this model is that it is viable only for stationary data. Stationary data means that the mean and variance of the data remain constant. Non stationary data must first be transformed using the methods above.
Before fitting and estimating the model, we need to know how many lags (the order), p, to include. One way of doing this is through plotting the autocorrelation function of the time series. This measures how much certain lags directly correlate with each other. Hence, we can deduce which lags are most statistically significant and remove the ones that are not when constructing our model. From the above autocorrelation plot, we can see that the ideal order would be 12.
from statsmodels.tsa.ar_model import AutoReg
from scipy.special import inv_boxcox
train = data.iloc[:-int(len(data) * 0.2)]
test = data.iloc[-int(len(data) * 0.2):]
model = AutoReg(train['Sales_box_cox_diff'], lags=12).fit()
C:\Users\aryad\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
transformed_forecasts = list(model.forecast(steps=len(test)))
boxcox_forecasts = []
for idx in range(len(test)):
if idx == 0:
boxcox_forecast = transformed_forecasts[idx] + train['Sales_box_cox'].iloc[-1]
else:
boxcox_forecast = transformed_forecasts[idx] + boxcox_forecasts[idx-1]
boxcox_forecasts.append(boxcox_forecast)
forecasts = inv_boxcox(boxcox_forecasts, lam)
fig6 = go.Figure()
fig6.add_trace(go.Scatter(x=train.index, y=train['Sales'], name='Train'))
fig6.add_trace(go.Scatter(x=test.index, y=test['Sales'], name='Test'))
fig6.add_trace(go.Scatter(x=test.index, y=forecasts, name='Forecast'))
fig6.update_layout(template="simple_white", font=dict(size=18), title_text='Autoregression',
width=650, title_x=0.5, height=400, xaxis_title='Date',
yaxis_title='Sales')
fig6.show()
C:\Users\aryad\anaconda3\Lib\site-packages\statsmodels\tsa\deterministic.py:308: UserWarning: Only PeriodIndexes, DatetimeIndexes with a frequency set, RangesIndexes, and Index with a unit increment support extending. The index is set will contain the position relative to the data length.
As we can see, the autoregressive model of order 12 [AR(12)] has done a decent job of predicting these values.
The moving average model is regression-like by fitting coefficients, θ, to the previously forecasted errors, ε, also known as white noise error, with the additon of a constant term that is the mean, μ:
This is a MA(q) model, where q is the number of error terms, which is known as the order.
One key requirement of the MA(q) model is that, like autoregression, it needs the data to be stationary. This can be achieved through differencing and stabilising the variance through a Logarithm or Box-Cox transform.
The autocorrelation function (ACF) of an MA(q) process is zero at lag q + 1 and greater. Therefore, we determine the appropriate maximum lag for the estimation by examining the sample autocorrelation function to see where it becomes insignificantly different from zero for all lags beyond a certain lag, which is designated as the maximum lag q. From our autocorrelation plot above, this is 24.
from statsmodels.tsa.arima.model import ARIMA
model2 = ARIMA(train['Sales_box_cox_diff'], order=(0, 1, 24)).fit()
transformed_forecasts = list(model2.forecast(steps=len(test)))
boxcox_forecasts = []
for idx in range(len(test)):
if idx == 0:
boxcox_forecast = transformed_forecasts[idx] + train['Sales_box_cox'].iloc[-1]
else:
boxcox_forecast = transformed_forecasts[idx] + boxcox_forecasts[idx-1]
boxcox_forecasts.append(boxcox_forecast)
forecasts = inv_boxcox(boxcox_forecasts, lam)
fig7 = go.Figure()
fig7.add_trace(go.Scatter(x=train.index, y=train['Sales'], name='Train'))
fig7.add_trace(go.Scatter(x=test.index, y=test['Sales'], name='Test'))
fig7.add_trace(go.Scatter(x=test.index, y=forecasts, name='Forecast'))
fig7.update_layout(template="simple_white", font=dict(size=18), title_text='Moving Average',
width=650, title_x=0.5, height=400, xaxis_title='Date',
yaxis_title='Sales')
fig7.show()
C:\Users\aryad\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used. C:\Users\aryad\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used. C:\Users\aryad\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used. C:\Users\aryad\anaconda3\Lib\site-packages\statsmodels\tsa\statespace\sarimax.py:978: UserWarning: Non-invertible starting MA parameters found. Using zeros as starting parameters. C:\Users\aryad\anaconda3\Lib\site-packages\statsmodels\base\model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
As we can see, the moving average model is not very accurate.
The moving average model and the autoregressive model are used together to build the ARIMA(Autoregressive Integrated Moving Average) model, which is the gold standard for time series forecasting.
Autoregressive Integrated Moving Average better known as ARIMA, is probably the most used time series forecasting model and is a combination of the individual aforementioned models.
The first part of ARIMA - AR stands for the Autoregressive model.
The third part - MA stands for the Moving Average model.
The middle part - I stands for Integrated: This is the number (order d) of differencing required to make the time series stationary.
The AR model is written as AR(p) and the MA model is written as MA(q) where p and q are the orders of the respective models, ARIMA is written as ARIMA(p,d,q) where p and q are the same orders above and d is the number of differencing steps required.
model3 = ARIMA(train['Sales_box_cox_diff'], order=(12, 1, 24)).fit()
transformed_forecasts = list(model3.forecast(steps=len(test)))
boxcox_forecasts = []
for idx in range(len(test)):
if idx == 0:
boxcox_forecast = transformed_forecasts[idx] + train['Sales_box_cox'].iloc[-1]
else:
boxcox_forecast = transformed_forecasts[idx] + boxcox_forecasts[idx-1]
boxcox_forecasts.append(boxcox_forecast)
forecasts = inv_boxcox(boxcox_forecasts, lam)
fig7 = go.Figure()
fig7.add_trace(go.Scatter(x=train.index, y=train['Sales'], name='Train'))
fig7.add_trace(go.Scatter(x=test.index, y=test['Sales'], name='Test'))
fig7.add_trace(go.Scatter(x=test.index, y=forecasts, name='Forecast'))
fig7.update_layout(template="simple_white", font=dict(size=18), title_text='ARIMA',
width=650, title_x=0.5, height=400, xaxis_title='Date',
yaxis_title='Sales')
fig7.show()
C:\Users\aryad\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used. C:\Users\aryad\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used. C:\Users\aryad\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used. C:\Users\aryad\anaconda3\Lib\site-packages\statsmodels\tsa\statespace\sarimax.py:966: UserWarning: Non-stationary starting autoregressive parameters found. Using zeros as starting parameters. C:\Users\aryad\anaconda3\Lib\site-packages\statsmodels\tsa\statespace\sarimax.py:978: UserWarning: Non-invertible starting MA parameters found. Using zeros as starting parameters. C:\Users\aryad\anaconda3\Lib\site-packages\statsmodels\base\model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
The ARIMA model has much better performance than AR or MA individually.