Additive Model Time-series Analysis using Python Machine Learning Client for SAP HANA

import matplotlib.pyplot as plt
plt.figure(figsize=(16, 8))
plt.plot(data['datetime'], data['system_demand_actual'])
plt.title('Figure 1. Overview of the dataset')
plt.show()
plt.close()

df_5wk = cc.sql('SELECT * FROM ({}) WHERE "datetime" <= \'2012-03-12 00:00:00\''.format(df.select_statement))
data_5wk = df_5wk.sort('datetime').collect()

plt.figure(figsize=(16, 4))
plt.plot(data_5wk['datetime'], data_5wk['system_demand_actual'])
plt.title('Figure 2: Zoom in view of data of 1st 5 weeks')
plt.show()
plt.close()

Additive models analysis is a new method that treats time-series modeling as a curve-fitting problem with respect to time. In contrast, exponential smoothing and ARIMA try model the dependencies of the current data with the past(inclusive of expected values and errors). Besides, additive model analysis use (partial) Fourier series to model periodic patterns in the input data. By doing so, the problem of non-integer periods is easily resolved, and the problem of multiple periods also become non-critical. Moreover, additive model analysis also allows users to input events(usually cyclical, called holidays in the Python API) that could potentially affect the the shape of the curve for given time-series during a short period of time and break up the regular seasonal pattern. This often leads to better modeling for time-series.

Now let us apply additive model analysis to our dataset of electricity system demand. As a standard procedure in machine learning, we firstly split the dataset into two parts for training and testing purposes, respectively.

df_train = cc.sql('SELECT * FROM ({}) WHERE "datetime" <= \'2016-02-18 00:00:00\''.format(df.select_statement))
df_train.sort('datetime').collect()

70704 rows × 2 columns

Now let us use the training set to fit a model for additive model analysis.

from hana_ml.algorithms.pal.tsa.additive_model_forecast import AdditiveModelForecast # import the python api for additive model analysis from hana_ml
amf = AdditiveModelForecast(seasonality_mode='additive', daily_seasonality='true', weekly_seasonality='true', yearly_seasonality='true')
amf.fit(df_train)

70704 rows × 2 columns

The model fitting procedure takes some time. After it is done, we may apply the trained additive model to the test data, in which case only time-stamps matter irrespective of the demand values.

from hana_ml.algorithms.pal.tsa.exponential_smoothing import AutoExponentialSmoothing
aesm = AutoExponentialSmoothing(forecast_model_name='TESM',
                                forecast_num=df_test.count(),
                                seasonal='additive',
                                optimizer_random_seed=1)
df_train_wid = df_train.add_id('ID', 'datetime')
aesm.fit_predict(data=df_train_wid, key='ID',
                 endog='system_demand_actual')

aesm.forecast_.sort('TIMESTAMP').collect()

df_wid = df.add_id('ID', 'datetime')
df_test_wid = cc.sql('SELECT * FROM ({}) WHERE ID >= {}'.format(df_wid.select_statement, df_train.count()))
aesm_join_res = df_test_wid.join(aesm.forecast_[['TIMESTAMP', 'VALUE']],##the VALUE column stores fitted & predicted demand values
                                 'ID = TIMESTAMP')
aesm_join_res.collect()

aesm_measure = accuracy_measure(data = aesm_join_res[['system_demand_actual', 'VALUE']], 
                                evaluation_metric = ['rmse', 'mad', 'mape'])
aesm_measure.collect()

fit_res = amf.predict(df_train) #obtain the fitting result for the training data
fit_res.sort('datetime').collect()

fit_res_rename = fit_res.rename_columns({'datetime':'date_time'})
fit_join = df_train.join(fit_res_rename, '"datetime" = "date_time"')
fit_join.collect()

residue = cc.sql('SELECT "datetime", YHAT - "system_demand_actual" AS RESID FROM ({})'.format(fit_join.select_statement))
import matplotlib.pyplot as plt
plt.figure(figsize=(16,5))
plt.plot(data['datetime'][0:df_train.count()], residue.collect()['RESID'])
plt.title('Figure 4. Overview of the random component of training data')
plt.show()
plt.close()

abnormal = cc.sql('SELECT * FROM ({}) ORDER BY RESID DESC LIMIT 100'.format(residue.select_statement))
abnormal_dates = abnormal[['datetime']].cast('datetime', 'DATE').drop_duplicates()
abnormal_dates.sort('datetime').collect()

holidays = cc.table('SG_HOLIDAY_TBL').cast('Date', 'DATE')#The Date column in holiday table is of type 'TIMESTAMP', convert to 'DATE'
holidays.head(5).collect()

amf.fit(data = df_train, holiday=holidays)

fit_res = amf.predict(df_train)
fit_res_rename = fit_res.rename_columns({'datetime':'date_time'})
fit_join = df_train.join(fit_res_rename, '"datetime" = "date_time"')
residue = cc.sql('SELECT "datetime", YHAT - "system_demand_actual" AS RESID FROM ({})'.format(fit_join.select_statement))
import matplotlib.pyplot as plt
plt.figure(figsize=(16,5))
plt.plot(data['date'][0:df_train.count()], residue.collect()['RESID'])
plt.title('Figure 4. Overview of random component of training data with holiday info.')
plt.show()
plt.close()

pred_res2 = amf.predict(df_test)
renm_pred2 = pred_res2.rename_columns({'datetime': 'date_time'})
join_pred2 = df_test.join(renm_pred2[['date_time', 'YHAT']], '"datetime"="date_time"')
from hana_ml.algorithms.pal.tsa.accuracy_measure import accuracy_measure
res = accuracy_measure(data = join_pred2[['system_demand_actual', 'YHAT']],
                       evaluation_metric = ['rmse', 'mad', 'mape'])
res.collect()

psg = cc.table('SG_INTER_VISITORS_TBL').cast('month', 'DATE')#'month' column with type 'TIMESTAMP', convert to 'DATE' type for easier processing
psg.collect()

train_num = 408
test_num = psg.count() - train_num
psg_train = psg.head(train_num)
psg_test = psg.sort('month', desc=True).head(test_num).sort('month')

psg_train_data = psg_train.collect()
import matplotlib.pyplot as plt
plt.figure(figsize=(16*(408/455),5))
plt.plot(psg_train_data['month'],
         psg_train_data['total_international_visitor_arrivals'])
plt.show()
plt.close()

event = cc.table('SG_SARS_MONTH_TBL').cast('Month', 'DATE')
event.collect()

It is seen from the training data that,  local variations increases with the growing number of international arrivals, which is typically regards a single a multiplicative seasonality(if there is any).

from hana_ml.algorithms.pal.tsa.additive_model_forecast import AdditiveModelForecast
amf = AdditiveModelForecast(seasonality_mode='multiplicative')
amf.fit(data=psg_train, holiday=event)

psg_fit = amf.predict(psg_train)# obtain fitted value for training data
psg_fit_rename_month = psg_fit.rename_columns({'month': 'Month'})
join_true_and_fit = psg_train.join(psg_fit_rename_month, '"month" = "Month"')
from hana_ml.algorithms.pal.tsa.accuracy_measure import accuracy_measure
fit_res = accuracy_measure(data = join_true_and_fit[['total_international_visitor_arrivals', 'YHAT']], 
                           evaluation_metric = ['rmse', 'mad', 'mape'])
fit_res.collect()

psg_pred = amf.predict(psg_test)
psg_pred_rename_month = psg_pred.rename_columns({'month': 'Month'})
join_true_and_pred = psg_test.join(psg_pred_rename_month, '"month" = "Month"')
psg_pred_acc = accuracy_measure(data = join_true_and_pred[['total_international_visitor_arrivals', 'YHAT']],
                                evaluation_metric = ['rmse', 'mad', 'mape'])
psg_pred_acc.collect()

Since the time-series of interest in this section has larger variations when its level value goes higher, so it is natual that MSE and RMSE have larger values on the test data since the level is expected to be higher within the prediction period. However, mean absolute percentage error is scale-invariant w.r.t. level of time-series, so compared to that of training data, its higher value on the test data strongly suggests that the training data could be over-fitted by the additive time-series model.

For better illustrattion, the real and predicted values are plotted as follows:

pred_data = psg_pred.collect()
psg_test_data = psg_test.collect()
import matplotlib.pyplot as plt
plt.figure(figsize=(16,5))
plt.plot(psg_test_data['month'], psg_test_data['total_international_visitor_arrivals'])
plt.plot(pred_data['month'], pred_data['YHAT'], color='red')
plt.show()
plt.close()

As discussed in the beginning of this section, when building the above additive model, a set of change-points is detected implicitly using a defaulted mechanism, and the detected change-points may not be optimal. So for the input time-series, after a careful examination of training data, we suggest the following set of changepoints:

‘1982-11-01’, ‘1986-01-01’, ‘1991-05-01’, ‘1995-09-01’, ‘1997-10-01’, ‘2000-07-01’, ‘2004-02-01’, ‘2007-12-01’, ‘2009-03-01’.

Let us first mark out the specified change-points in the training data via run-sequence plot.

import pandas as pd
psg_train_data = psg_train.collect()
changepoints = list(pd.to_datetime(pd.Series(['1982-11-01', '1986-01-01', '1991-05-01', '1995-09-01', '1997-10-01', '2000-07-01', '2004-02-01', '2007-12-01', '2009-03-01'])))
event_points = list(pd.to_datetime(pd.Series(['2003-04-01', '2003-05-01', '2003-06-01'])))
cp_data = psg_train_data.loc[[x in changepoints for x in psg_train_data['month']],:]
event_data = psg_train_data.loc[[x in event_points for x in psg_train_data['month']],:]

plt.figure(figsize=(int(16*(408/455)),5))
plt.plot(psg_train_data['month'], psg_train_data['total_international_visitor_arrivals'])
plt.scatter(cp_data['month'], cp_data['total_international_visitor_arrivals'], color='red')
plt.scatter(event_data['month'], event_data['total_international_visitor_arrivals'], color='purple')
plt.legend(['Visitor Arrivals', 'Specified Change-points', 'Event Points'])
plt.show()
plt.close()

Next, we manually specify these points as change-points when building up the additive model for time-series analysis.

amf = AdditiveModelForecast(seasonality_mode='multiplicative',
                            changepoints=['1982-11-01 00:00:00',
                                          '1986-01-01 00:00:00',
                                          '1991-05-01 00:00:00',
                                          '1995-09-01 00:00:00',
                                          '2000-07-01 00:00:00',
                                          '2004-02-01 00:00:00',
                                          '2007-12-01 00:00:00',
                                          '2009-03-01 00:00:00'])
amf.fit(data=psg_train, holiday=event)

Again let us check how well the model fits the training data.

psg_fit2 = amf.predict(psg_train)
psg_fit2_rename_month = psg_fit2.rename_columns({'month': 'Month'})
join_true_and_fit2 = psg_train.join(psg_fit2_rename_month, '"month" = "Month"')
from hana_ml.algorithms.pal.tsa.accuracy_measure import accuracy_measure
fit_res2 = accuracy_measure(data = join_true_and_fit2[['total_international_visitor_arrivals', 'YHAT']], 
                            evaluation_metric = ['rmse', 'mad', 'mape'])
fit_res2.collect()

pred2_data = psg_pred2.collect()
plt.figure(figsize=(16,5))
plt.plot(psg_test_data['month'],
         psg_test_data['total_international_visitor_arrivals'])
plt.plot(pred2_data['month'], pred2_data['YHAT'], color='red')
plt.legend(['Actual Arrivals', 'Predicted Arrivals'])
plt.show()
plt.close()