概要
feature-engineがきっかけで最近チェックしているPyDATA
この動画を見ながら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を比較用によく使用するが、これにも勝てないことはザラで、原因は経験的に下記
- stationalityがない
white noiseやrandom workな波形は予測が困難で、targetを説明する他のfeaturesが必要となる。
correlogramや1回差分(diff)を取ったそれで簡単に確認可能 - training用とval用のdatasetで傾向が異なる
簡単にunder/over fitしてしまう。
windowの広さやtraining/valの期間の選び方が重要だし、shuffleできないから恣意的に精度を上げることも可能
(Kaggleでup voteもらいやすくなる) - datasetが時間方向に少ない
モデルがパターンを見つけられない。 - datasetがfeature方向に少ない
周期性のはっきりした気温のような予測ならfeatureは過去の気温一択でも良いが、株価では通じないので今回tsfreshに助けてもらう。
ざっとドキュメントを確認すると期待できそうなのは大きく次の2点
- Feature extracting: 時系列関連のmetricsを自動計算
- Feature filtering: その中でtargetに関係しそうなmetricsの選定
ここで言うmetricsがfeaturesになる。
前者については力技でなんとかなるかとも思っていたが実際計算してくれる種類がハンパなく力技でも無理。この多さをもってドメイン知識の少なさをカバーする。
それ以上に期待するのが後者で、そこにはFRESH(FeatuRe Extraction based on Scalable Hypothesis tests)というアルゴリズムを使用しているとのこと。元ペーパのアブストは下記の通り。
“特徴選択における全関連問題とは、強く関連する属性と弱く関連する属性を全て同定することである。この問題は、予知保全や生産ライン最適化などの産業応用における時系列分類や回帰において、各ラベルや回帰対象が複数の時系列やメタ情報と同時に関連するため、特に解決が困難である。本論文では、機械学習パイプラインの初期段階において、分類や回帰タスクに対する重要度に関して利用可能な特徴をフィルタリングし、選択されたが無関係な特徴の割合を制御する、効率的でスケーラブルな時系列向け特徴抽出アルゴリズムを提案する。提案アルゴリズムは、確立された特徴抽出手法と特徴重要度フィルタを組み合わせたものである。このアルゴリズムは計算量が少なく、限られた領域に関する知識しかなくても問題に取り組むことができ、並列化が可能で、スケーラブルであり、よく研究されたノンパラメトリック仮説検定に基づくものである。我々は、UCR時系列分類アーカイブの全ての二値分類問題、および、生産ライン最適化プロジェクトからの時系列と、ダイナミクスの質的変化を基礎とする確率過程のシミュレーションに対して、提案するアルゴリズムをベンチマークする。”
by DeepL
問題
今回の問題は次のような感じ
手元にある対象銘柄の過去ohlcデータのみ使用し、t+1のclose値を予測する。
このままではfeaturesが4つしかない(実際は簡単のためOpenとCloseの2つ)のでtsfreshにかさ増ししてもらう。
悪いことに、今回使用するIntel株の値動きに95%信頼区間でstationalityが無いことを棄却できなかったので精度は期待していない。手順と精度の相対比較のみとなる。
Dataset準備
1. 使用するパケージ
!pip install tsfresh
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は列のひとつにしている。
## 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')
一応確認する。
df_wk.info()
infoのフォーマットが少し変わった気がする。
<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かさ増し
流れは
- 独自フォーマットにdfを整形
- そのdfをwindow単位にrolling
- rollingしたdfをextraction
- yを与え、yを説明しているfeatureをfiltering
Documentsを読んでも1が特にわからなかった。ココに、3typeの期待するDataFrame型が書かれてあるが、そもそもyの意味がわからない。時系列metricsを計算するのにtargetは関係なく、なんでyが必要？
さらにyを入れてrollingした後、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する。
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
featuersがもっとあれば更にconcat()するだけ。
もしかするとidは'0'とか、ひとつに固定して、column_kindで'intl_Open', 'intl_Close'を識別させるのかもしれないが、使い分けがDocumentsだけでは今ひとつわからず無視。
2. そのdfをwindow単位にrolling
rolling前に最小かさ増し量のMinimalFCParametersでdefault_fc_parametersに指定してextractionさせてみる。
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
10個のmetricsが100日間について計算されてるっぽい。一応確認してみる。
cudaエラーはColaboのGPUを有効にすると消える。「並列化が可能」とアブストにあるようにGPUを使って並列計算しているのだろう。
print(Close_df['X'].sum(), Close_df['X'].mean())
print(Open_df['X'].sum(), Open_df['X'].mean())
707.5300000000001 7.075300000000001
711.4 7.114
やりたいのはforecastなので任意の幅のwindowをスライドさせdatasetを準備しなければならない。いつもはpandasのshiftや、feature_engine.timeseries.forecastingのLagFeaturesを使うところだが、ここではtsfreshにやてもらう。
どのsampleもwindowは同じ大きさにしたいのでmax/min_timeshiftに同じ3を指定してみる。
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]
新しいidは(元id, window先頭のDate)のタプルで表記される。
Documentに説明されているように3time stepが1windowになりそうだが、実際は何故か(3+1)time stepになるみたい。それ以外は期待通りの結果
# tupleになったidで指定するときの表記
rolled_df[rolled_df['id'] == ('intl_Close', pd.to_datetime('2005-06-14'))]
3. rollingしたdfをextraction
windowを9+1に広げたrolled_dfからextractionしてもらう。
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
features_df.info()
<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で抜き取り確認する。
rolled_df[rolled_df['id'] == ('intl_Close', pd.to_datetime('2005-03-10'))]['X'].sum()
78.85であっている。
nanが含まれる列を調べる。
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()という関数も用意されているようだが、説明読むと乱暴な処理がされるみたいなので使わないことにする。
このままfilteringするとnanが入っているからダメよ、とエラーが出ることがあるので、impute()する。
impute(features_df)
4. yを与え、yを説明しているfeatureをfiltering
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
777個のmetricsから75個に減っている。features_filteredについてnanを確認する。
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営業日にしたものを使用する。
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()
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で分ける。
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させる。ハイパラ調整は省略する。
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させ、カレンダー情報も追加したもの。
!pip install feature-engine
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を準備する。
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が出来たか確認する。
X_train.columns
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で予測精度を確認する。
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が十分に多ければ精度比較に影響しないと考える。
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化する。
!pip install shap
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))
同様に"2. originalのfeatures"の結果が下図となる。
1.はcwt関連が上位に来ていて、何やらContinuous Wavelet Transform（和文）した結果らしい。
2.は'intl_Close'やそのlag成分が上位に来るかと思ったら、'intl_Low', 'intl_High'が上位だった。この2つについてもlag成分を追加したり、その日の波高値、High - Low なんかも入れてみたい。
以上
2023年のresolution
キャリアアップのため、KaggleでnotebookとdiscussionでExpertをもらうこと！