Model Complexity と頑健性

### モデルの複雑さと、データに軽いノイズを乗っけたときのモデル精度の悪化を検証

モデル構築に用いたデータは、過去起こり得た事象(population)の一つのサンプルであったという発想を持ちがちな私として、ただ一つの過去に見えたサンプルが多少ぶれうる事象であったことも気にかけたくなるもの。「サンプルが多少ぶれることを考えること」＝「目に見えたサンプルににノイズを足す」とすれば、自分の作るモデルの頑健性(データへの依存度)を確認する一つの方法として使えるのはないかと思い、作成

ざっくり言うと
テストデータ、学習データに分けるののではなく、これだ！！というモデルを作り、本当はちょっとデータ違い得るです。というときにモデルは精度がどれだけ悪化するか
という現象を確認したかったというのが動機

* 予測対象：USDJPYの翌週リターン
* 予測の用いるデータ：過去4週分のリターン
* 予測方法：PolynomialFeatures & LinearRegression
* 複雑さとして、PolynomialFeaturesに用いるdegreeを1～6と設定
* 足すノイズ：元データの0.3倍のvolatilityを持つ正規分布乱数

モデルシンプルで、ノイズもシンプルなものなら、数式を用いて予測精度の劣化は計算できると思うが、Simulationベースにすることで、多様なモデル、i.i.dでも正規分布でもないノイズにでも対応できるのではないかと

Import など

```python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# import statsmodels.api as sm
# from statsmodels.tsa.arima_model import ARIMA
import pandas_datareader.data as web
import yfinance as yf
from numpy.random import *

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression

import seaborn as sns
import warnings
import sys

warnings.filterwarnings('ignore')

plt.style.use('seaborn-darkgrid')
plt.rcParams['axes.xmargin'] = 0.01
plt.rcParams['axes.ymargin'] = 0.01
```

yfinanceから、'USDJPY'を取得し、週次リターン作成

```python
ReadDF = yf.download('JPY=X', start="1995-01-01", end="2019-10-30")
ReadDF.index = pd.to_datetime(ReadDF.index)
IndexValueReadDF_rsmpl = ReadDF.resample('W').last()['Adj Close']
ReadDF = IndexValueReadDF_rsmpl / IndexValueReadDF_rsmpl.shift(1) - 1
```

説明変数(X)、被説明変数(y)データ作成

```python
ret_df = pd.DataFrame()
ret_df[mkt] = ReadDF

test = pd.DataFrame(ret_df[mkt])
for i in range(1, 5):
    test['i_' + str(i)] = test['USDJPY'].shift(1 * i)

test = test.dropna(axis=0)
X = test.ix[:, 1:]
y = test[mkt]
```

モデルの作成(実サンプル&ノイズを乗っけたもの)

```python
m_ = 3  # 0.3倍のvolatilityを持たせるための掛け目
output_degree = {}

for k in range(1, 7):
    polynomial_features = PolynomialFeatures(degree=k, include_bias=False)
    linear_regression = LinearRegression()
    pipeline = Pipeline([("polynomial_features", polynomial_features),
                         ("linear_regression", linear_regression)])

    pipeline.fit(X, y)

    k_sample = pd.DataFrame()

    for l in range(0, 300):

        # 正規乱数発生
        eps_0 = pd.DataFrame(randn(X.shape[0], X.shape[1]))
        eps_1 = eps_0.apply(lambda x: x * list(X.std()), axis=1)

        eps_1.columns = X.columns
        eps_1.index = X.index

        # 元データに足す
        X_r = X + m_/10 * eps_1

        signal = pd.DataFrame()
        signal[mkt] = np.sign(pd.DataFrame(pipeline.predict(X_r)))[0]

        signal.index = y.index

        k_sample['s_' + str(l)] = (pd.DataFrame(signal[mkt]) * pd.DataFrame(y)).ix[:, 0]

    signal_IS = pd.DataFrame()
    # signal_IS[mkt] = np.sign(pd.DataFrame(pipeline.predict(X.ix[:, 0][:, np.newaxis])))[0]

    signal_IS[mkt] = np.sign(pd.DataFrame(pipeline.predict(X)))[0]

    signal_IS.index = y.index

    k_sample['IS'] = (pd.DataFrame(signal_IS[mkt])*pd.DataFrame(y)).ix[:, 0]
    k_sample[mkt] = pd.DataFrame(y).ix[:, 0]

    output_degree['degree_' + str(k)] = k_sample
```

アウトプットヒストグラム(degree N毎)

* グレー線：元データ(USDJPY)のリターン/リスク
* 赤線：元サンプルへモデル適用した戦略のリターン/リスク
* 青線：ノイズを乗っけたサンプルへモデル適用した戦略の平均リターン/リスク
* ヒストグラム：青線の元データ
* グレー線：

```python

Performance_sim = pd.DataFrame()

fig = plt.figure(figsize=(15, 7), dpi=80)
for k in range(1, 7):
    ax = fig.add_subplot(2, 3, k)

    for_stats = output_degree['degree_' + str(k)]
    Performance = pd.DataFrame(for_stats.mean()*50 / (for_stats.std()*np.sqrt(50))).T

    Performance_tmp = Performance.ix[:, 1:].T

    ax.hist(Performance.drop(['IS', mkt], axis=1), bins=30, color="dodgerblue", alpha=0.8)
    ax.axvline(x=float(Performance.drop(['IS', mkt], axis=1).mean(axis=1)), color="b")
    ax.axvline(x=float(Performance['IS']), color="tomato")
    ax.axvline(x=float(Performance[mkt]), color="gray")
    ax.set_ylim([0, 40])
    ax.set_xlim([-0.3, 2.5])
    ax.set_title('degree-N polynomial: ' + str(m_))
fig.show()

```

####　結果
![Model_Complexity.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/509334/ed7477cc-fb64-9bdc-d524-dffa9234fcec.png)

* 同然インサンプルでも、週次リターンの簡単なモデルではUSDJPYは予測できない→外れてもパフォーマンスそんなかわらない(もともと精度高くないので)
* インサンプルで複雑なモデルは当然当てはまるものを見つけられるが、激しく元データに依存

シンプルなデータだと当然な結果だが、自分がモデルを作った際の検証には役立ちそう

### 追記1
上では、モデルの複雑さをPolynomialFeaturesのdegree=kとしたが、degree=3、モデルをRandomForestRegressorとし、モデルの複雑さをmax_depth=kとする

```python
    pipeline = Pipeline([("polynomial_features", polynomial_features),
                         ("rf_regression", RandomForestRegressor(max_depth=k))])
```

####　結果
![Model_Complexity_RF.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/509334/5ccb7c39-7c94-7de1-b181-581cd42bc185.png)


### 追記2
もう一度、モデルの複雑さをPolynomialFeaturesのdegree=kとする例で、説明変数を追加し、キリキリもモデルを作りたくなる気持ちをシステマチックに制御するという意味で、回帰部分にLassoを使って説明変数の選択を実施

```python
    pipeline = Pipeline([("polynomial_features", polynomial_features),
                         ("linear_regression", LassoCV(cv=5))])
```
![Model_Complexity_LassoCV.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/509334/c2ca181a-84f9-3260-f30d-fec546ac2b96.png)

* そもそも、次週リターンを予測するのは難しく、予測できたというのはデータに過度に依存(Over-fittingの可能性高い)しおた結果であったというのがつかめる。