XGBoostとtsfreshでMultivariate time series forecasting 備忘録

# 概要
feature-engineがきっかけで最近チェックしているPyDATA

https://www.youtube.com/watch?v=9QtL7m3YS9I&list=LL&index=3&t=1631s

この動画を見ながらKaggleの株価予測をsktimeで試してみるか、ってモチベが上がって（よくある）久しぶりのsktimeを触っていたらtsfreshをラップしたapiに気づいたので、先にtsfreshの素性を調べてみるか、と回りくどい流れ。
時系列初心者にはマストな内容の講演動画で超オススメ。ハッ！と思わせるやり取りが質疑応答にあってとても良かった。

今回はおもにtsfreshを使って説明変数を増やすことが目的で、XGBoostは精度確認しSHAPでimportanceを可視化する程度のおまけ。
Daskのように独自の馴染めない世界観がありDocumentsの例もわかりにくく、少し手間取ったので残す。

- 実施期間: 2022年12月
- 環境：Colab
- パケージ：tsfresh, XGBoost, SHAP, feature-engine

## tsfresh
例えば、株価のClose値予測といった周期性がなさそうな問題では、ML modelの精度は悪くなる。
baselineに1time step前の値を予測値としたnaiveなmodelを比較用によく使用するが、これにも勝てないことはザラで、原因は経験的に下記
1. stationalityがない
white noiseやrandom workな波形は予測が困難で、targetを説明する他のfeaturesが必要となる。
correlogramや1回差分(diff)を取ったそれで簡単に確認可能
1. training用とval用のdatasetで傾向が異なる
簡単にunder/over fitしてしまう。
windowの広さやtraining/valの期間の選び方が重要だし、shuffleできないから恣意的に精度を上げることも可能
(Kaggleでup voteもらいやすくなる)
1. datasetが時間方向に少ない
モデルがパターンを見つけられない。
1. datasetがfeature方向に少ない
周期性のはっきりした気温のような予測ならfeatureは過去の気温一択でも良いが、株価では通じないので今回tsfreshに助けてもらう。

ざっとドキュメントを確認すると期待できそうなのは大きく次の2点
1. Feature extracting: 時系列関連のmetricsを自動計算
1. Feature filtering: その中でtargetに関係しそうなmetricsの選定

ここで言うmetricsがfeaturesになる。
前者については力技でなんとかなるかとも思っていたが実際計算してくれる種類がハンパなく力技でも無理。この多さをもってドメイン知識の少なさをカバーする。
それ以上に期待するのが後者で、そこにはFRESH(FeatuRe Extraction based on Scalable Hypothesis tests)というアルゴリズムを使用しているとのこと。[元ペーパ](https://ui.adsabs.harvard.edu/abs/2016arXiv161007717C/abstract)のアブストは下記の通り。

_“特徴選択における全関連問題とは、強く関連する属性と弱く関連する属性を全て同定することである。この問題は、予知保全や生産ライン最適化などの産業応用における時系列分類や回帰において、各ラベルや回帰対象が複数の時系列やメタ情報と同時に関連するため、特に解決が困難である。本論文では、機械学習パイプラインの初期段階において、分類や回帰タスクに対する重要度に関して利用可能な特徴をフィルタリングし、選択されたが無関係な特徴の割合を制御する、効率的でスケーラブルな時系列向け特徴抽出アルゴリズムを提案する。提案アルゴリズムは、確立された特徴抽出手法と特徴重要度フィルタを組み合わせたものである。このアルゴリズムは計算量が少なく、限られた領域に関する知識しかなくても問題に取り組むことができ、並列化が可能で、スケーラブルであり、よく研究されたノンパラメトリック仮説検定に基づくものである。我々は、UCR時系列分類アーカイブの全ての二値分類問題、および、生産ライン最適化プロジェクトからの時系列と、ダイナミクスの質的変化を基礎とする確率過程のシミュレーションに対して、提案するアルゴリズムをベンチマークする。”_
 by DeepL

## 問題
今回の問題は次のような感じ

手元にある対象銘柄の過去ohlcデータのみ使用し、t+1のclose値を予測する。
このままではfeaturesが4つしかない(実際は簡単のためOpenとCloseの2つ)のでtsfreshにかさ増ししてもらう。



悪いことに、今回使用するIntel株の値動きに95%信頼区間でstationalityが無いことを棄却できなかったので精度は期待していない。手順と精度の相対比較のみとなる。


# Dataset準備
## 1. 使用するパケージ

```terminal
!pip install tsfresh
```

```python
import os
from datetime import datetime
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
```

## 2. 元DataFrame作成
KaggleからDLしたintl.us.txtだけをGoogle driveにおいている前提。
本番用に複数銘柄を使用していたのでfor loopになっている。time seriesもののときはDataFrameのindexをdatetime型にしているが（コメントアウト部）、tsfreshはそれを期待していないのでdatetimeは列のひとつにしている。

```python
## For Google Drive mount
from google.colab import drive
drive.mount('/content/drive')

file_path = '/content/drive/My Drive/Colab Notebooks/csv'
file_names = os.listdir(file_path)

ticker_name = []
flag = True
for i, fname in enumerate(file_names):
    print(fname)
    ticker_name.append(fname.split('.')[0])
    df_wk = pd.read_csv(os.path.join(file_path, fname), header=0)
    # df_wk = pd.read_csv(os.path.join(file_path, fname), header=0, index_col='Date', parse_dates=True)
    df_wk = df_wk.drop('OpenInt', axis=1)
    df_wk = df_wk.rename(columns=
                         {'Open': ticker_name[i] + '_Open',
                          'High': ticker_name[i] + '_High',
                          'Low': ticker_name[i] + '_Low',
                          'Close': ticker_name[i] + '_Close',
                          'Volume': ticker_name[i] + '_Volume'})
    df_wk['Date'] = pd.to_datetime(df_wk['Date'], format='%Y-%m-%d')
```

一応確認する。
```python
df_wk.info()
```
infoのフォーマットが少し変わった気がする。
```terminal
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 3192 entries, 0 to 3191
Data columns (total 6 columns):
 #   Column       Non-Null Count  Dtype         
---  ------       --------------  -----         
 0   Date         3192 non-null   datetime64[ns]
 1   intl_Open    3192 non-null   float64       
 2   intl_High    3192 non-null   float64       
 3   intl_Low     3192 non-null   float64       
 4   intl_Close   3192 non-null   float64       
 5   intl_Volume  3192 non-null   int64         
dtypes: datetime64[ns](1), float64(4), int64(1)
memory usage: 149.8 KB
```

# tsfreshでFeaturesかさ増し
流れは
1. 独自フォーマットにdfを整形
1. そのdfをwindow単位にrolling
1. rollingしたdfをextraction
1. yを与え、yを説明しているfeatureをfiltering

Documentsを読んでも1が特にわからなかった。[ココ](https://tsfresh.readthedocs.io/en/latest/text/data_formats.html#data-formats-label)に、3typeの期待するDataFrame型が書かれてあるが、そもそもyの意味がわからない。時系列metricsを計算するのにtargetは関係なく、なんでyが必要？
さらにyを入れてrollingした後、[make_forecasting_frame()](https://tsfresh.readthedocs.io/en/latest/api/tsfresh.utilities.html#tsfresh.utilities.dataframe_functions.make_forecasting_frame)を使って「特異な時系列xを受け取り、時系列予測タスクに使用できるDataFrame dfとターゲットベクトルyを構築します。」らしいが、yを指定するには明らかに引数が少なすぎる。filteringするときに改めてyを渡すことにして、ここではyをガン無視する。


## 1. 独自フォーマットにdfを整形
普通はかさ増しネタとなるfeaturesがintl_Open列やintl_Close列になるが、試行錯誤の結果それらはidで識別させて１列につなげよ、と読んだ。つまり下記の処理となる。
いっぱいあっても読みづらいので100営業日分の下記を使用する。つまりOpen値とClose値から1 time step先（翌営業日）のClose値をforecastする。

```python
Open_df = df_wk.iloc[:100][['Date', 'intl_Open']]
Open_df.rename(columns={'intl_Open': 'X'}, inplace=True)    # 列名は統一する。
Open_df['id'] = 'intl_Open'
Close_df = df_wk.iloc[:100][['Date', 'intl_Close']]
Close_df.rename(columns={'intl_Close': 'X'}, inplace=True)    # 列名は統一する。
Close_df['id'] = 'intl_Close'

df = pd.concat([Open_df, Close_df])
df
```
![Screenshot from 2022-12-31 11-15-19.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/6c1b521c-e0e1-d428-269d-28a5a74c79e7.png)
featuersがもっとあれば更にconcat()するだけ。

もしかするとidは'0'とか、ひとつに固定して、column_kindで'intl_Open', 'intl_Close'を識別させるのかもしれないが、使い分けがDocumentsだけでは今ひとつわからず無視。



## 2. そのdfをwindow単位にrolling
rolling前に最小かさ増し量のMinimalFCParametersでdefault_fc_parametersに指定してextractionさせてみる。

```python
from tsfresh.feature_extraction import extract_features
from tsfresh.feature_extraction import ComprehensiveFCParameters
from tsfresh.feature_extraction import EfficientFCParameters
from tsfresh.feature_extraction import MinimalFCParameters

settings = MinimalFCParameters()
features_df = extract_features(df, default_fc_parameters=settings, column_id='id', column_sort='Date')
features_df
```
![Screenshot from 2022-12-31 11-36-38.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/79c14666-a0ba-e4ef-409e-a08005941362.png)
10個のmetricsが100日間について計算されてるっぽい。一応確認してみる。
cudaエラーはColaboのGPUを有効にすると消える。「並列化が可能」とアブストにあるようにGPUを使って並列計算しているのだろう。

```python
print(Close_df['X'].sum(), Close_df['X'].mean())
print(Open_df['X'].sum(), Open_df['X'].mean())
```

```terminal
707.5300000000001 7.075300000000001
711.4 7.114
```

やりたいのはforecastなので任意の幅のwindowをスライドさせdatasetを準備しなければならない。いつもはpandasのshiftや、feature_engine.timeseries.forecastingの[LagFeatures](https://feature-engine.readthedocs.io/en/1.3.x/api_doc/timeseries/forecasting/LagFeatures.html#feature_engine.timeseries.forecasting.LagFeatures)を使うところだが、ここではtsfreshにやてもらう。

どのsampleもwindowは同じ大きさにしたいのでmax/min_timeshiftに同じ3を指定してみる。

```python
from tsfresh.utilities.dataframe_functions import roll_time_series

rolled_df = roll_time_series(
    df, column_id='id', column_sort='Date', max_timeshift=3, min_timeshift=3
    )
rolled_df.iloc[:20]
```
![Screenshot from 2022-12-31 16-14-44.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/98558edc-1b01-7c73-615a-981604077739.png)

新しいidは(元id, window先頭のDate)のタプルで表記される。
[Documentに説明](https://tsfresh.readthedocs.io/en/latest/api/tsfresh.utilities.html?highlight=roll_time_series#tsfresh.utilities.dataframe_functions.roll_time_series)されているように3time stepが1windowになりそうだが、実際は何故か(3+1)time stepになるみたい。それ以外は期待通りの結果

```python
# tupleになったidで指定するときの表記
rolled_df[rolled_df['id'] == ('intl_Close', pd.to_datetime('2005-06-14'))]
```
![Screenshot from 2022-12-31 16-23-50.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/baa8ab8e-a88e-484a-7989-8996beee3e17.png)

## 3. rollingしたdfをextraction
windowを9+1に広げたrolled_dfからextractionしてもらう。

```python
rolled_df = roll_time_series(
    df, column_id='id', column_sort='Date', max_timeshift=9, min_timeshift=9
    )

# settings = ComprehensiveFCParameters()
# settings = MinimalFCParameters()
settings = EfficientFCParameters()
features_df = extract_features(
    rolled_df, default_fc_parameters=settings, column_id='id', column_sort='Date', column_value='X'
    )
features_df
```

![Screenshot from 2022-12-31 16-42-55.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/c754ef55-e734-33ad-7844-9268bd8305ec.png)

```python
features_df.info()
```
```terminal
<class 'pandas.core.frame.DataFrame'>
MultiIndex: 182 entries, ('intl_Close', Timestamp('2005-03-10 00:00:00')) to ('intl_Open', Timestamp('2005-07-22 00:00:00'))
Columns: 777 entries, X__variance_larger_than_standard_deviation to X__mean_n_absolute_max__number_of_maxima_7
dtypes: float64(777)
memory usage: 1.1+ MB
```

777個のmetricsが得られた。X__sum_valuesで抜き取り確認する。

```python
rolled_df[rolled_df['id'] == ('intl_Close', pd.to_datetime('2005-03-10'))]['X'].sum()
```
78.85であっている。

nanが含まれる列を調べる。

```python
nan_clmn_lst = [(features_df[s].isna().sum(), s) for s in features_df.columns if features_df[s].isna().sum()!=0]
len(nan_clmn_lst)
```
452個の列でほぼすべてのwindowがnanになっている。~~[impute()](https://tsfresh.readthedocs.io/en/latest/api/tsfresh.utilities.html?highlight=impute#tsfresh.utilities.dataframe_functions.impute)という関数も用意されているようだが、説明読むと乱暴な処理がされるみたいなので使わないことにする。~~
このままfilteringするとnanが入っているからダメよ、とエラーが出ることがあるので、impute()する。

```python
impute(features_df)
```

## 4. yを与え、yを説明しているfeatureをfiltering

```python
y_arr = Close_df[Close_df['Date']>='2005-03-10']['X'].to_numpy()
# 使えそうなfeatureを選んでもらう
features_filtered_Open_df = select_features(features_df.loc['intl_Open'], y_arr)
features_filtered_Close_df = select_features(features_df.loc['intl_Close'], y_arr)
features_filtered_Open_df
```
![Screenshot from 2022-12-31 17-36-09.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/0f61d207-5546-eac7-8016-ea007c5baa0c.png)

777個のmetricsから75個に減っている。features_filteredについてnanを確認する。

```python
nan_clmn_lst = [(features_filtered[s].isna().sum(), s) for s in features_filtered.columns if features_filtered[s].isna().sum()!=0]
len(nan_clmn_lst)
nan_clmn_lst
```
[]

# XGBoostでPrediction
## 1. tsfreshのかさ増しfeaturesでpredict
tsfreshは一通りわかったので、100営業日を下図の1000営業日にしたものを使用する。

```python
fig, axes = plt.subplots(figsize = (10, 5), tight_layout=True)
ax = df_wk['intl_Open'].iloc[1000:2000].plot()
df_wk['intl_Close'].iloc[1000:2000].plot(ax=ax)
ax.legend()
```

![Screenshot from 2023-01-02 11-54-00.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/e0127657-155a-a13d-7bba-20cb2f2b435c.png)

```python
Open_df = df_wk.iloc[1000:1999][['Date', 'intl_Open']]
Open_df.rename(columns={'intl_Open': 'X'}, inplace=True)    # 列名は統一する。
Open_df['id'] = 'intl_Open'
Close_df = df_wk.iloc[1000:1999][['Date', 'intl_Close']]
Close_df.rename(columns={'intl_Close': 'X'}, inplace=True)    # 列名は統一する。
Close_df['id'] = 'intl_Close'

df = pd.concat([Open_df, Close_df])

# yもついでに準備
temp_df = df_wk.iloc[1001:2000]    # 1time step先をforecastするから
y_arr = temp_df[temp_df['Date']>='2009-03-17']['intl_Close'].to_numpy()
print(f'y: {len(y_arr)}')
```

intl_Openもintl_Closeも列名はXのだったため、かさ増ししたfeatures_filtered_Open_dfとfeatures_filtered_Close_dfの列名はほぼ重複している。重複ないように列名変更後マージしてそれをtrainとvalに8:2で分ける。

```python
clmn_lst = list(features_filtered_Open_df.columns)
dic= {s: 'Open_' + s for s in clmn_lst}
features_filtered_Open_df2 = features_filtered_Open_df.rename(columns=dic)

clmn_lst = list(features_filtered_Close_df.columns)
dic= {s: 'Close_' + s for s in clmn_lst}
features_filtered_Close_df2 = features_filtered_Close_df.rename(columns=dic)

features_filtered_df = pd.merge(features_filtered_Open_df2, features_filtered_Close_df2, right_index=True, left_index=True)
features_filtered_df

X_train, X_test, y_train, y_test = features_filtered_df.iloc[:800], features_filtered_df.iloc[800:], y_arr[:800], y_arr[800:]
```

XGBoost modelにfitさせる。ハイパラ調整は省略する。

```python
model = XGBRegressor(random_state=123, objective='reg:squarederror')
model.fit(X_train, y_train)

y_pred = model.predict(X_test)
print('MSE: {:.3f}'.format(mean_squared_error(y_test, y_pred, squared=False)))
```
MSE: 0.539

## 2. オリジナルのfeaturesだけでpredict
比較のため時系列予測するときに使う、いつもの最小セットでpredictする。
'intl_Open', 'intl_Close'を同じくwindow幅10time stepでrollingさせ、カレンダー情報も追加したもの。

```python
!pip install feature-engine
```

```python
from feature_engine.timeseries.forecasting import LagFeatures

def rolling_by_LagFeatures(wk_df, clmn_lst, upto=6):
    lag_f = LagFeatures(variables=clmn_lst ,periods=[i for i in range(1,upto)])
    wk_df_tr = lag_f.fit_transform(wk_df)
    wk_df_tr.dropna(axis=0, how='any', inplace=True)
    return wk_df_tr


def make_calendar_df(date_df, clm_name):
    wk_df = pd.DataFrame(columns=['day', 'month', 'year', 'quarter', 'dayofweek', 'weekofyear'])
    wk_df['day'] = date_df[clm_name].dt.day
    wk_df['month'] = date_df[clm_name].dt.month
    wk_df['year'] = date_df[clm_name].dt.year
    wk_df['quarter'] = date_df[clm_name].dt.quarter
    wk_df['dayofweek'] = date_df[clm_name].dt.dayofweek
    wk_df['weekofyear'] = date_df[clm_name].dt.weekofyear
    wk_df = wk_df.reset_index(drop=True)
    return wk_df
```

tsfreshのときと同じ範囲でdatasetを準備する。

```python
y_arr = Close_df[Close_df['Date']>='2009-03-17']['X'].to_numpy()

org_df = df_wk.iloc[1000:1999]
org_date_df = df_wk.iloc[1001:2000]    # 1time step先の日時
org_date_df = make_calendar_df(org_date_df, 'Date')
org_df = org_df.reset_index(drop=True)
org_df = pd.merge(org_df, org_date_df, left_index=True, right_index=True)

X_org = rolling_by_LagFeatures(org_df, ['intl_Open', 'intl_Close'], 10)
X_org = X_org.set_index('Date', drop=True)
X_train, X_test, y_train, y_test = X_org.iloc[:800], X_org.iloc[800:], y_arr[:800], y_arr[800:]
```

intl_Open, intl_Colseについて10time stepのwindowが出来たか確認する。

```python
X_train.columns
```

```terminal
Index(['intl_Open', 'intl_High', 'intl_Low', 'intl_Close', 'intl_Volume',
       'day', 'month', 'year', 'quarter', 'dayofweek', 'weekofyear',
       'intl_Open_lag_1', 'intl_Close_lag_1', 'intl_Open_lag_2',
       'intl_Close_lag_2', 'intl_Open_lag_3', 'intl_Close_lag_3',
       'intl_Open_lag_4', 'intl_Close_lag_4', 'intl_Open_lag_5',
       'intl_Close_lag_5', 'intl_Open_lag_6', 'intl_Close_lag_6',
       'intl_Open_lag_7', 'intl_Close_lag_7', 'intl_Open_lag_8',
       'intl_Close_lag_8', 'intl_Open_lag_9', 'intl_Close_lag_9'],
      dtype='object')
```

同様にこのdatasetで予測精度を確認する。

```python
model_org = XGBRegressor(random_state=123, objective='reg:squarederror')
model_org.fit(X_train, y_train)
y_pred = model_org.predict(X_test)
print('MSE: {:.3f}'.format(mean_squared_error(y_test, y_pred, squared=False)))
```
MSE: 0.418


## 3. baselineの疑似predict
targetのyを1time stepずらしy_predとする。つまり直前のyを予測値にするというもの。
y =      [1, 2, 3, 4, 5]
y_pred = [nan, 1, 2, 3, 4, 5]
従いこの例のように先頭と、後尾のy_predには該当するXがないのでデータ点数は2つ少なくなってしまう。val datasetが十分に多ければ精度比較に影響しないと考える。

```python
y_pred_lst = list(y_test)
y_pred_lst.insert(0,0)
y_pred_lst = y_pred_lst
print(len(y_test), len(y_pred_lst))

print('MSE: {:.3f}'.format(mean_squared_error(y_test[1:], y_pred_lst[1:-1], squared=False)))
```
191 192
MSE: 0.345

## 4. 結果比較

| No. | Dataset | MSE | # features |
|:-:|:-:|:-:|:-:|
| 1  | tsfreshのかさ増しfeatures  | 0.539  | 252 |
| 2  | originalのfeatures | 0.418  | 29 |
| 3  | baseline  | **0.345**  | - |

window size = 10 time step
samples for training = 790

使用した区間(i.e. 1000:2000)を変えると精度は大きく変わるので一概に1.の結果が悪かったとも言えない。多分features数が多い割にwindow幅が小さかったからかもしれない。
window幅、sample数、マルチコ調査とfeature削減(GBDTなので効かないと思うけど)、calendarの追加など行うと精度は必ず上がるはずだけど、2.に勝てるか微妙かも…
続きはKaggleで追い込むこととする。

# SHAPのfeature importance
"1. tsfreshのかさ増しfeatures"で作ったval datasetとXGBoost modelを使ってshapでviz化する。

```terminal
!pip install shap
```

```python
exp = shap.TreeExplainer(model)
explainer = shap.Explainer(model)
shap_exp = explainer(X_test)        # explainerに、.expected_valueはない
shap_values = shap_exp.values

shap.summary_plot(shap_values, X_test, max_display=10, plot_size=(10, 8))
```

![Screenshot from 2023-01-03 14-11-43.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/08b0f401-77c7-c4fe-1a33-d871015185db.png)

同様に"2. originalのfeatures"の結果が下図となる。

![Screenshot from 2023-01-03 14-13-07.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/2769718c-a0d4-6f60-233f-b3195ee7ad7f.png)

1.はcwt関連が上位に来ていて、何やら[Continuous Wavelet Transform](https://pywavelets.readthedocs.io/en/latest/ref/cwt.html)（[和文](https://ja.wikipedia.org/wiki/%E3%82%A6%E3%82%A7%E3%83%BC%E3%83%96%E3%83%AC%E3%83%83%E3%83%88%E5%A4%89%E6%8F%9B)）した結果らしい。
2.は'intl_Close'やそのlag成分が上位に来るかと思ったら、'intl_Low', 'intl_High'が上位だった。この2つについてもlag成分を追加したり、その日の波高値、High - Low なんかも入れてみたい。

以上

### 2023年のresolution
キャリアアップのため、KaggleでnotebookとdiscussionでExpertをもらうこと！