Time Series Forecasting

Modelling the Electric production time series data

Photo by Leeloo Thefirst

This project deals with prediction of consumption of electricity in the coming future. Time series analysis is important to know trend seasonality and other characteristics of the data for efficient modelling.

The dataset is downloaded from kaggle here.

About Dataset

The dataset has 2 columns which represents date and consumption percentage.

Import Libraries and load the dataset

# Import libraries
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.tsa.stattools import adfuller,acf, pacf
from statsmodels.tsa.arima_model import ARIMA
import statsmodels.api as sm 
from pylab import rcParams
# set plot size
rcParams['figure.figsize'] = 10, 6
# read data
folder_path = '/content/drive/MyDrive/ETE_projects/time_series/'
df = pd.read_csv(folder_path+ 'Electric_Production.csv')
df.head()

Visualize the data

def visualization(data):
  print(data.info())
  print('-------------------------------------------------------')
  print(data.shape)
  print('-------------------------------------------------------')
  print(data.isnull().sum())

visualization(df)

Change the date to datetime.

df['DATE'] = pd.to_datetime(df.DATE)
df = df.set_index(df.DATE)
df.drop('DATE', axis=1, inplace=True)
print(df.dtypes)
df.head()

df.tail(3)

df.describe()

Plot Trend and Seasonality

## Time series Characteristics
def characteristics(data, x, y, title):
  ax = sns.lineplot(data = data, lw=1)
  ax.lines[0].set_linestyle('solid')
  ax.set_title(title)

  df_temp = data.copy()
  df_temp['Year'] = pd.DatetimeIndex(df_temp.index).year
  df_temp['Month'] = pd.DatetimeIndex(df_temp.index).month
  # Stacked Line Plot
  plt.figure(figsize=(8,8))
  plt.title('Seasonality of Time Series')
  sns.set(style="ticks", rc={"lines.linewidth": 0.7})
  sns.pointplot(x='Month', y= 'IPG2211A2N', hue='Year', data = df_temp,palette='mako', linestyles='-' )

characteristics(df, 'DATE', 'IPG2211A2N','Trend of Series')

• We observe upward trend the graph.
• There is also slight increase in the variance over time.(Length of the spikes)

• We also observe that the data is seasonal and there is dip during April-May(4–5) and rise during July-August(7–8) and November- December(11–12)
• The data is seasonal yearly as well as half-yearly.

Check for Stationarity

A time series is said to be stationary if the mean, variance and autocorrelation does not vary over time.

Trend, seasonality and other characteristics varying according to time will affect modelling as underlying pattern of the data is not learnt by the model due to non-stationarity. Thus stationarity needs to be achieved.

The following are the methods for stationarity test:

• Rolling Window Analysis
• Dickey Fuller Test

Rolling Window Analysis

• It refers to calculating the values based on previous ‘window size’ of values. It analyses if the variation is time dependent and if yes for what size of the window.
• This helps use to know the nature of coefficients of the data with respect to time. Our aim is to achieve parameter stability of the data.
• As we have found that the data is seasonal in the above graph, we set the window size = 12.

Augmented Dickey Fuller Test

On an explanatory perspective,

• Null hypothesis H0 implies that unit root=1, the time series is not stationary.
• Alternate hypothesis H1 implies that unit root <1, and the series is stationary.

We us test statistic value, critical values at 1%, 5% and 10% confidence intervals and p value for evaluating the test results.

Get Pooja’s stories in your inbox

Join Medium for free to get updates from this writer.

Remember me for faster sign in

When the test statistic < = critical values, p value < 0.05 indicates that null hypothesis can be rejected and data is stationary.

Stationary Test

sns.set_style('whitegrid')
col = 'IPG2211A2N'
def stationarity_test(data, window, title, col_name):
  data['r_mean'] = data[col_name].rolling(window = window).mean()
  data['r_std'] = data[col_name].rolling(window=window).std()
  # Plot rolling statistic
  ax = sns.lineplot(data = data, lw=1, palette=['navy', 'darkgreen', 'red'])
  ax.lines[0].set_linestyle('solid')
  ax.set_title(title)
  # Dickey Fuller test
  df_test = adfuller(data[col_name])
  df_output = pd.Series(df_test[0:4], index=['Test Statistic',
                                             'p-value',
                                             'Lags used for t-statistic',
                                             'No_of_observation_used'])
  for key, value in df_test[4].items():
    df_output['Critical value(%s)'%key] = value
  print(df_output)

stationarity_test(df, 12, 'Original Data',col)

• We observe that standard deviation is some what constant, but mean is not constant.

• T values corresponding to ADF test. Since critical value -2.25>-3.44,-2.86,-2.57 (t-values at 1%,5%and 10% confidence intervals), null hypothesis cannot be rejected.
• Hence there is non stationarity in the data.
• Also p-value of 0.18>0.05(if we take 5% significance level or 95% confidence interval), null hypothesis is accepted.
• Hence data is non stationary (that means it has relation with time)

Converting Non-Stationary Data into Stationary Data

Differencing, mathematical transformations like log, weighted moving average or combination of these methods can be used to make data stationary.

Differencing

It refers to substracting the value with its previous value(positive lag) in the time series. Lag is specified for differencing. Eg:

Differenced data is subjected to stationarity test.

def differencing(data,order):
  df_diff = data.diff(periods=order)
  # Data transformation may add na values
  df_diff.dropna(inplace=True)
  return df_diff

df_diff = differencing(df, 1)
stationarity_test(df_diff,12, 'Stationary Test for First Order Differenced Data',col )

• Mean is constant.
• There is variation and upward trend in standard deviation.

• Test statistic Value is less than (1%, 5%, 10%) critical.
• p-value is lesser than 0.05.
• Thus it is stationary series.
• There is variance along with time

Let’s experiment with other transformations.

Log Transformation

df_log = np.log(df)
stationarity_test(df_log,12, 'Stationary Test for Log Data',col )

• Log transformed series is stationary with 95% confidence only. Test statistic is lesser than 5%(1–0.05 = 0.95) but greater than 1% critical.

Differencing Log Data

df_log_diff = differencing(df_log, 1)
stationarity_test(df_log_diff,12,'Stationary Test for Log Differenced Data',col)

• Log transformed and differenced series is stationary with test statistic is lesser than 1%, 5% and 10% critical.
• p value is less than 0.05.
• mean is constant and standard deviation is slightly varying with no significant trend.

We use log differenced data for modelling. As we have found that the data series is seasonal, we use SARIMAX for modelling.

SARIMAX

In order to model with SARIMAX we need to find (p,d,q) and seasonal component (P,D,Q) of the series.

• p represents order of auto regression (AR), q represents the order of moving average (MA) and d represents order of differencing.

ACF and PACF Plots

• PACF graph indicates the correlation of the y and lag y at t after removing the effect of previous lags where t represents the lag order.
• ACF graph indicates the correlation of the lag y with its own lags at order t.

from statsmodels.graphics.tsaplots import plot_acf,plot_pacf

fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(df_log_diff[col].iloc[13:],lags=20,ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(df_log_diff[col].iloc[13:],lags=20,ax=ax2)

q=3 as we find significant lag at 3 and not beyond in ACF. p=4 as initially 4 lags are significant in PACF graph.

Finding P, Q D

We decompose the time series data into components to visualize the characteristics. Multiplicative model is used as the seasonal variation increases over time.

Seasonal Decomposition using Multiplicative model

col = 'IPG2211A2N'
sns.set_style('whitegrid')
from statsmodels.tsa.seasonal import seasonal_decompose
result_m = seasonal_decompose(df[col], model='multiplicative',extrapolate_trend='freq')
result_m.plot()
plt.show() 

Check for stationarity of seasonal Component

from statsmodels.tsa.stattools import adfuller
import numpy as np
def check_stationarity(ts):
    dftest = adfuller(ts)
    adf = dftest[0]
    pvalue = dftest[1]
    critical_value = dftest[4]['5%']
    if (pvalue < 0.05) and (adf < critical_value):
        print('The series is stationary')
    else:
        print('The series is NOT stationary')

seasonal_a = result_m.seasonal
check_stationarity(seasonal_a)

The series is stationary

Plot ACF and PACF for seasonal component.

fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(seasonal_a,lags=6,ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(seasonal_a,lags=6,ax=ax2) 

• All the lags are outside confidence interval in ACF. So we consider Q=6 after which it is cyclic in nature.
• P=3 is chosen considering lags at 1,2 and 6 in PACF graph.
• D=0 as we are not differencing the seasonal components.

Modelling using SARIMAX

Train test split

df_train = df_log[:278]
df_test = df_log[279:]
# Copy of data
df_b = df_log.copy(deep=True)

Fit the model

model_b = sm.tsa.statespace.SARIMAX(df_train['IPG2211A2N'],order=(4,1,3), seasonal_order = (3, 0, 6,12))
results_b = model_b.fit()

Get prediction for the test data

pred_b= results_b.predict(start = df_test.index[0], end = df_test.index[-1])
df_test['pred_b'] = pred_b.values

Plot the result

sns.lineplot(data = df_b, x = df_b.index, y = col, color='aqua', lw=1.1) 
sns.lineplot(data = df_test, x = df_test.index, y = 'pred_b', color='darkorange', lw=1.1)

Evaluation

• Mean Squared Error(MSE): It is the average of the sum of the squares of the distance between actual values and predicted values. Lower values indicates the forecast is closer to actual.
• Root mean squared error(RMSE): It is square root of MSE. It is the average of difference between the predicted and actual values. It represents the accuracy of the prediction.
• Both MSE and RMSE is used to measure accuracy.
• Mean absolute error: It is average of the sum of absolute difference between original and predicted values.

mse = mean_squared_error(df_test[col], df_test['pred_b'])
print('mse:',mse)
rmse = np.sqrt(mean_squared_error(df_test[col], df_test['pred_b']))
print('rmse:', rmse)

mae = np.mean(np.abs(results_b.resid))
print('MAE: %.3f' % mae)

Residual Diagnostic Plots

Residuals refers to the difference obtained by substracting predicted value from actual value.

results_b.plot_diagnostics()
plt.show() 

• Standardized residual plot: This plot does not necessarily have a pattern which indicates that the model is functioning well.
• Histogram and density plot: This is the density plot of the residuals. The normal distribution depicted in green and the residual density plot depicted in orange has small differences. If there is no difference it indicates the model is best. Here the model is better.
• Q_Q Plot: The distribution of the residuals appears to be normal with few data points at the ends.
• Correlogram Plot: The ACF plot of the residuals indicates that there is no significant correlation outside the confidence interval, indicating that most of the information is captured by the model. Thus it is a good model.

Summary Statistics

results_b.summary()

• The null hypothesis for Ljung-Box test states that there is no correlation in the residuals.
• Prob(Q) is 0.92>0.05 indicates that null hypothesis is accepted. The residuals are independent and not correlated.
• The null hypothesis for Jarque-Bera test states that the residuals are guassian normally distributed.
• Prob(JB) is 0.00<0.05 indicating that the null hypothesis is rejected. Thus the residuals are not normally distributed.

Thus residuals are not correlated and residuals are not normally distributed.

Prediction with Future Data

Create the Dataframe for future data.

futureDate = pd.DataFrame(pd.date_range(start='2018-02-01', end='2021-01-01', freq='MS'), columns=['Dates'])
futureDate.set_index('Dates', inplace=True)
future_df=pd.concat([df_log,futureDate]) 

Predict and plot the prediction.

future_df['forecast_b'] = results_b.predict(start='2018-02-01', end='2021-01-01', dynamic= True) 
sns.lineplot(data = future_df, x = future_df.index, y = col, color='royalblue', lw=1)
sns.lineplot(data = future_df, x = future_df.index, y = future_df['forecast_b'], color='darkorange', lw=1)

Conclusion

The following steps and techniques are applied for this time series forecasting:

• Visualize the data
• Check for stationarity using rolling statistics and augmented dickey fuller test
• Convert non stationary to stationary data by applying tranaformation and differencing techniques.
• find p, d, q using acf, pacf plots
• Seasonal Decomposition using multiplicative model
• Find P D Q values of the seasonal components
• Modelling using SARIMAX
• Evaluation using MSE, RMSE, MAE.
• Interpret Residual Diagnostic Plots
• Interpret Summary Statistics
• Prediction with future data.

References:

https://www.kaggle.com/code/nholloway/seasonality-and-sarimax/notebook

https://towardsdatascience.com/understanding-the-seasonal-order-of-the-sarima-model-ebef613e40fa

https://www.kaggle.com/code/satishgunjal/tutorial-time-series-analysis-and-forecasting/notebook

https://www.jadsmkbdatalab.nl/forecasting-with-sarimax-models/

https://www.mathworks.com/help/econ/rolling-window-estimation-of-state-space-models.html

https://online.stat.psu.edu/