<jupyter_start><jupyter_text>Contents1. [Background](1)2. [Setup the Environment](2)3. [Prepare the Data](3)4. [Algorithm Comparison](4)   1. Background Time series forecast is a very common problem in many real-world applications. A wide spectrum of algorithms have been proposed to solve this problem. However, it is usually difficult to benchmark different algorithms and compare their performance due to the various implementation of the algorithms. This notebook tries to show an example how to benchmark different time series forecast algorithms by only using the Glutonts library.   2. Setup the Environment  2.1 install R for extenal r forecast algorithms<jupyter_code># intall the missing lib for R
!sudo yum install libXt-1.1.4-6.1.9.amzn1.x86_64 -y
# install python r interface
!conda install -c r rpy2==2.9.4 --yes
# install forecast R packages
!R -e 'install.packages(c("forecast", "nnfor"), repos="https://cloud.r-project.org")'<jupyter_output><empty_output><jupyter_text>2.2 Install Python Pacakages<jupyter_code># install Prophet python packages
!conda install -c plotly plotly==3.10.0 --yes
!conda install -c conda-forge fbprophet=0.5=py36he1b5a44_0 --yes
!conda install -c anaconda ujson=1.35=py36h14c3975_0 --yes
# install gluonts
!pip install gluonts==0.4.2
# install mxnet
!pip install mxnet==1.4.1<jupyter_output><empty_output><jupyter_text>2.3 import libraries<jupyter_code>%matplotlib inline

import matplotlib.pyplot as plt
import numpy as np
import os
import pandas as pd
import random
import sys
import zipfile
import mxnet as mx

from urllib.request import urlretrieve

from gluonts.model.trivial.mean import MeanPredictor
from gluonts.model.seasonal_naive import SeasonalNaivePredictor
from gluonts.model.r_forecast import RForecastPredictor
from gluonts.model.prophet import ProphetPredictor
from gluonts.model.deepar import DeepAREstimator
from gluonts.trainer import Trainer
from gluonts.dataset.common import ListDataset
from itertools import islice
from gluonts.evaluation.backtest import make_evaluation_predictions
from gluonts.evaluation import Evaluator
# fix some plot issues caused by Prophet model
# pls refer to https://darektidwell.com/typeerror-float-argument-must-be-a-string-or-a-number-not-period-facebook-prophet-and-pandas/
pd.plotting.register_matplotlib_converters()
# set random seeds for reproducibility
np.random.seed(42)
random.seed(42)
mx.random.seed(42)<jupyter_output><empty_output><jupyter_text>3. Prepare the Data Import electricity dataset, we need to download the original data set of from the UCI data set repository.  3.1 Download data<jupyter_code>DATA_HOST = "https://archive.ics.uci.edu"
DATA_PATH = "/ml/machine-learning-databases/00321/"
ARCHIVE_NAME = "LD2011_2014.txt.zip"
FILE_NAME = ARCHIVE_NAME[:-4]
def progress_report_hook(count, block_size, total_size):
    mb = int(count * block_size // 1e6)
    if count % 500 == 0:
        sys.stdout.write("\r{} MB downloaded".format(mb))
        sys.stdout.flush()

if not os.path.isfile(FILE_NAME):
    print("downloading dataset (258MB), can take a few minutes depending on your connection")
    urlretrieve(DATA_HOST + DATA_PATH + ARCHIVE_NAME, ARCHIVE_NAME, reporthook=progress_report_hook)

    print("\nextracting data archive")
    zip_ref = zipfile.ZipFile(ARCHIVE_NAME, 'r')
    zip_ref.extractall("./")
    zip_ref.close()
else:
    print("File found skipping download")
data = pd.read_csv(FILE_NAME, sep=";", index_col=0, parse_dates=True, decimal=',')
num_timeseries = data.shape[1]
data_kw = data.resample('2H').sum() / 8
timeseries = []
for i in range(num_timeseries):
    timeseries.append(np.trim_zeros(data_kw.iloc[:,i], trim='f'))
fig, axs = plt.subplots(5, 2, figsize=(20, 20), sharex=True)
axx = axs.ravel()
for i in range(0, 10):
    timeseries[i].loc["2014-01-01":"2014-01-14"].plot(ax=axx[i])
    axx[i].set_xlabel("date")    
    axx[i].set_ylabel("kW consumption")   
    axx[i].grid(which='minor', axis='x')<jupyter_output><empty_output><jupyter_text>3.2 Train Test Data Split<jupyter_code>def split_train_test_data(timeseries,
                          start_dataset,
                          end_training,
                          num_test_windows):
    # create training data.
    training_data = [
    {
        "start": str(start_dataset),
        "target": ts[start_dataset:end_training - 1].tolist(),  # We use -1, because pandas indexing includes the upper bound
        "feat_static_cat": [id]
    }
    for id, ts in enumerate(timeseries)
    ]
    
    # create testing data.
    test_data = [
        {
            "start": str(start_dataset),
            "target": ts[start_dataset:end_training + k * prediction_length].tolist(),
            "feat_static_cat": [id]
        }
        for k in range(1, num_test_windows + 1)
        for id, ts in enumerate(timeseries)
    ]
    
    return training_data, test_data
# we use 2 hour frequency for the time series
freq = '2H'

# we predict for 1 day
prediction_length = 1 * 12

# we also use 7 days as context length, this is the number of state updates accomplished before making predictions
context_length = 7 * 12

# The moving window for forecast
num_test_windows = 1

# training/test Split
start_dataset = pd.Timestamp("2014-01-01 00:00:00", freq=freq)
end_training = pd.Timestamp("2014-09-01 00:00:00", freq=freq)

# number of time series selected
n_timeseries = 50
training_data, test_data = split_train_test_data(timeseries[:n_timeseries],
                          start_dataset,
                          end_training,
                          num_test_windows)<jupyter_output><empty_output><jupyter_text>4. Algorithm Comparison<jupyter_code>def train_and_test(training_data, 
                   test_data,
                   freq,
                   num_test_windows,
                   model,
                   train_per_ts=False,
                   require_train=False
                   ):
    forecasts = []
    tss = []
    # if training per time series is required.
    if train_per_ts:
        # iterate over the timeseries
        count = 0
        for training_ts in training_data:
            # get the training time series
            training_ts = ListDataset(
                      [training_ts],
                      freq = freq
            )
            # get the related testing time series
            test_tss = test_data[count*num_test_windows: (count+1)*num_test_windows]

            test_tss = ListDataset(
                test_tss,
                freq = freq
            )
            if require_train:
                predictor = model.train(training_data=training_ts)
            else:
                predictor = model
            
            forecast_it, ts_it = make_evaluation_predictions(test_tss, predictor=predictor, num_samples=100)
            forecasts.extend(list(forecast_it))
            tss.extend(list(ts_it))
            count += 1
    else:
        training_data = ListDataset(
                      training_data,
                      freq = freq
            )
        test_data = ListDataset(
            test_data,
            freq = freq
        )
        if require_train:
            predictor = model.train(training_data=training_data)
        else:
            predictor = model
        
        forecast_it, ts_it = make_evaluation_predictions(test_data, predictor=predictor, num_samples=100)
        forecasts.extend(list(forecast_it))
        tss.extend(list(ts_it))
        
    return forecasts, tss<jupyter_output><empty_output><jupyter_text>4.1 Training and Testing   Mean<jupyter_code>%%time
mean= MeanPredictor(freq=freq, prediction_length=prediction_length,
                                  context_length=context_length)
forecasts_mean, tss_mean = train_and_test(training_data, 
               test_data,
               freq,
               num_test_windows,
               mean,
               train_per_ts=True,
               require_train=False
               )<jupyter_output><empty_output><jupyter_text>Seaonal Naive<jupyter_code>%%time
seasonal = SeasonalNaivePredictor(freq=freq,
                                  prediction_length=prediction_length,
                                  season_length=context_length)
forecasts_seasonal, tss_seasonal = train_and_test(training_data, 
               test_data,
               freq,
               num_test_windows,
               seasonal,
               train_per_ts=True,
               require_train=False
               )<jupyter_output><empty_output><jupyter_text>ARIMA<jupyter_code>%%time
arima = RForecastPredictor(freq=freq,
                           prediction_length=prediction_length,
                           method_name='arima')
forecasts_arima, tss_arima = train_and_test(training_data, 
               test_data,
               freq,
               num_test_windows,
               arima,
               train_per_ts=True,
               require_train=False
               )<jupyter_output><empty_output><jupyter_text>Prophet<jupyter_code>%%time
prophet = ProphetPredictor(freq, prediction_length)
forecasts_prophet, tss_prophet = train_and_test(training_data, 
               test_data,
               freq,
               num_test_windows,
               prophet,
               train_per_ts=True,
               require_train=False
               )<jupyter_output><empty_output><jupyter_text>DeepAR<jupyter_code>%%time
deepar = DeepAREstimator(freq=freq,
                         use_feat_static_cat=True,
                         cardinality=[n_timeseries],
                        prediction_length=prediction_length,
                        trainer=Trainer(epochs=100),
                        num_cells=40)
forecasts_deepar, tss_deepar = train_and_test(training_data, 
               test_data,
               freq,
               num_test_windows,
               deepar,
               train_per_ts=False,
               require_train=True
               )<jupyter_output><empty_output><jupyter_text>4.2 Results  Evaluation<jupyter_code>evaluator = Evaluator(quantiles=[0.5], seasonality=None)
agg_metrics_mean, item_metrics_mean = evaluator(iter(tss_mean), iter(forecasts_mean), num_series=len(forecasts_mean))
print(agg_metrics_mean)
evaluator = Evaluator(quantiles=[0.5], seasonality=None)
agg_metrics_seasonal, item_metrics_seasonal = evaluator(iter(tss_seasonal), 
                                                        iter(forecasts_seasonal), 
                                                        num_series=len(forecasts_seasonal))
print(agg_metrics_seasonal)
evaluator = Evaluator(quantiles=[0.5], seasonality=None)
agg_metrics_arima, item_metrics_arima = evaluator(iter(tss_arima), 
                                                        iter(forecasts_arima), 
                                                        num_series=len(forecasts_arima))
print(agg_metrics_arima)
evaluator = Evaluator(quantiles=[0.5], seasonality=None)
agg_metrics_prophet, item_metrics_prophet = evaluator(iter(tss_prophet), 
                                                        iter(forecasts_prophet), 
                                                        num_series=len(forecasts_prophet))
print(agg_metrics_prophet)
evaluator = Evaluator(quantiles=[0.5], seasonality=None)
agg_metrics_deepar, item_metrics_deepar = evaluator(iter(tss_deepar), 
                                                        iter(forecasts_deepar), 
                                                        num_series=len(forecasts_deepar))
print(agg_metrics_deepar)
df_metrics = pd.concat(
    [
        pd.DataFrame.from_dict(agg_metrics_deepar, orient='index').rename(columns={0: "DeepAR"}),
     pd.DataFrame.from_dict(agg_metrics_prophet, orient='index').rename(columns={0: "Prophet"}),
     pd.DataFrame.from_dict(agg_metrics_arima, orient='index').rename(columns={0: "ARIMA"}),
    pd.DataFrame.from_dict(agg_metrics_seasonal, orient='index').rename(columns={0: "Seasonal naive"}),
    pd.DataFrame.from_dict(agg_metrics_mean, orient='index').rename(columns={0: "Mean"})], axis=1
)
df_metrics.loc[["MASE", "RMSE", "sMAPE"]]<jupyter_output><empty_output><jupyter_text>Plot Example Forecast<jupyter_code>def plot_forecasts(tss, forecasts, past_length, start, stop, step, title):
    for target, forecast in islice(zip(tss, forecasts), start, stop, step):
        ax = target[-past_length:].plot(figsize=(12, 5), linewidth=2)
        forecast.plot(color='g')
        plt.title(title)
        plt.grid(which='both')
        plt.legend(["observations", "median prediction", "90% confidence interval", "50% confidence interval"])
        plt.show()
start, stop, step = 10, 11, 1
plot_forecasts(tss_mean, forecasts_mean, past_length=100, start=start, stop=stop, step=step, title="mean")
plot_forecasts(tss_seasonal, forecasts_seasonal, past_length=100, start=start, stop=stop, step=step, title="seasonal")
plot_forecasts(tss_arima, forecasts_arima, past_length=100, start=start, stop=stop, step=step, title="arima")
plot_forecasts(tss_prophet, forecasts_prophet, past_length=100, start=start, stop=stop, step=step, title="prophet")
plot_forecasts(tss_deepar, forecasts_deepar, past_length=100, start=start, stop=stop, step=step, title="deepar")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09456801331517628, 'fraction_numerical': 0.017627477681948858, 'mean_word_length': 2.7532651902328222, 'pattern_counts': {'":': 7, '<': 53, '<?xml version=': 0, '>': 53, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 24, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '249248', 'license': 'MIT-0', 'n_tokens_mistral': 4215, 'n_tokens_neox': 3981, 'n_words': 815, 'path': 'amazon-sagemaker-forecast-algorithms-benchmark-using-gluonts.ipynb', 'repo_name': 'pinak-p/amazon-sagemaker-forecast-algorithms-benchmark-using-gluonts', 'size': '13218'}