Amazon SageMaker と GluonTS で株価の予測をしてみよう

pip install --upgrade mxnet==1.7
# pip install --upgrade mxnet-cu101 GPUの場合
pip install gluonts

from datetime import datetime, timedelta
import numpy as np
import pandas as pd


START_DATE = datetime(year=2018, month=1, day=4)
END_DATE = datetime(year=2019, month=12, day=30)
TIMEDELTA = END_DATE - START_DATE


# TIMEDELTAの日数で銘柄コードと土日を除いたBusinessDayのグリッドを作成
keys, dates = np.meshgrid(
    [stock_price['Local Code'].unique()],
    [(START_DATE + timedelta(days=i)) for i in range(TIMEDELTA.days + 1) if (START_DATE + timedelta(days=i)).weekday() not in [5,6]]
)

business_day = pd.DataFrame({
    "Local Code": keys.ravel(), 
    "EndOfDayQuote Date": dates.ravel(),
})

# 作成したグリッドにstock_priceをjoinする
stock_dataset = pd.merge(
    business_day, 
    stock_price, 
    on=["Local Code", "EndOfDayQuote Date"], 
    how="left").sort_values(
    ["Local Code", "EndOfDayQuote Date"]).reset_index(drop=True)

from gluonts.dataset.common import ListDataset
from gluonts.dataset.field_names import FieldName


# 予測期間を20営業日
prediction_length = 20

# データをindexが銘柄コード、カラムが日付、終値が値となるようにPivot展開
stock_dataset = stock_dataset.pivot(index='Local Code', columns='EndOfDayQuote Date', values='EndOfDayQuote Close')

# 学習は最後の20営業日を除いたデータ、テストは全て含まれたデータ
train_target_values = [ts[:-prediction_length] for ts in stock_dataset.values]
test_target_values = stock_dataset.copy().values

# データの開始日付を銘柄数分作成
dates = [pd.Timestamp("2018-01-04", freq='1B') for _ in range(stock_dataset.shape[0])]

# 学習、テストListDatasetを作成
train_ds = ListDataset([
    {
        FieldName.TARGET: target,
        FieldName.START: start,
        FieldName.ITEM_ID: code,
    }
    for (target, start, code) in zip(train_target_values, dates, stock_dataset.index)
], freq="1B")

test_ds = ListDataset([
    {
        FieldName.TARGET: target,
        FieldName.START: start,
        FieldName.ITEM_ID: code,
    }
    for (target, start, code) in zip(test_target_values, dates, stock_dataset.index)
], freq="1B")

from gluonts.model.deepar import DeepAREstimator
from gluonts.mx.distribution.student_t import StudentTOutput
from gluonts.mx.trainer import Trainer


estimator = DeepAREstimator(
    freq='1B', # 必須
    prediction_length=prediction_length, # 必須
    trainer = Trainer(batch_size=32,
                      clip_gradient=10.0, 
                      ctx='cpu',
                      epochs=100, 
                      hybridize=True, 
                      init="xavier", 
                      learning_rate=0.001, 
                      learning_rate_decay_factor=0.5, 
                      minimum_learning_rate=5e-05, 
                      num_batches_per_epoch=50, 
                      patience=10,
                      weight_decay=1e-08), 
    context_length = prediction_length, 
    num_layers = 2,
    num_cells = 40, 
    cell_type = 'lstm', 
    dropout_rate = 0.1, 
    use_feat_dynamic_real = False, 
    use_feat_static_cat = False, 
    use_feat_static_real = False, 
    cardinality = None, 
    embedding_dimension = None, 
    distr_output = StudentTOutput(), 
    scaling = True, 
    lags_seq = None, 
    time_features = None, 
    num_parallel_samples = 100, 
)

predictor = estimator.train(train_ds)

from gluonts.evaluation.backtest import make_evaluation_predictions

forecast_it, ts_it = make_evaluation_predictions(
    dataset=test_ds,                                             
    predictor=predictor,
    num_samples=100
)

print("Obtaining time series conditioning values ...")
tss = list(tqdm(ts_it, total=len(test_ds)))

print("Obtaining time series predictions ...")
forecasts = list(tqdm(forecast_it, total=len(test_ds)))

from gluonts.evaluation import Evaluator

evaluator = Evaluator(quantiles=[0.5])
agg_metrics, item_metrics = evaluator(iter(tss), iter(forecasts), num_series=len(test_ds))

print(json.dumps(agg_metrics, indent=4))

item_metrics.head()

import matplotlib.pyplot as plt
import os


plot_log_path = "./plots/"
directory = os.path.dirname(plot_log_path)
if not os.path.exists(directory):
    os.makedirs(directory)


def plot_prob_forecasts(ts_entry, forecast_entry, path, sample_ide, inline=True):
    plot_length = 40
    prediction_intervals = (50, 90)
    legend = ["observations", "median prediction"] + [f"{k}% prediction interval" for k in prediction_intervals][::-1]

    _, ax = plt.subplots(1, 1, figsize=(10, 7))
    ts_entry[-plot_length:].plot(ax=ax)
    forecast_entry.plot(prediction_intervals=prediction_intervals, color='g')
    ax.axvline(ts_entry.index[-prediction_length], color='r')
    plt.legend(legend, loc="upper left")
    plt.title(forecast_entry.item_id)
    if inline:
        plt.show()
        plt.clf()
    else:
        plt.savefig('{}forecast_{}.png'.format(path, forecast_entry.item_id))
        plt.close()


print("Plotting time series predictions ...")
for i in tqdm(range(5)):
    ts_entry = tss[i]
    forecast_entry = forecasts[i]
    plot_prob_forecasts(ts_entry, forecast_entry, plot_log_path, i, inline=True)