ドル円為替レートの推論値算出

## はじめに  

　ドルを資産運用として考えるときに、期待する利益としてインカムゲインとキャピタルゲインがある。インカムゲインはドル金利により得られ運用リスクはあまりないが、キャピタルゲインはドル購入のタイミングにより収益にもなり得るし、損失にもなり得る。
　円高の時にドルを購入し、円安の時にドルを売却することでキャピタルゲインを得ることができるが、今が円高なのか円安なのかの判断材料を提供したい。
　新聞等で日米金利差によって為替レートが説明されることがある。為替レートはそれだけでなく様々な要因の影響を受けている。
 　今回入手可能な経済データをもとに為替レートの推論値を算出し、指標毎の影響度を比較してみることとする。

　なお、本投稿はスクールで学んだ各手法を用いて分析することを主目的としており、分析結果は必ずしも実務的な正確性や汎用性を保証するものではないことをご了承いただきたい。

## 推論値算出結果
ネットから入手できるデータを用いてドル円為替レートの推論値を算出した。経済収支が為替レートに最も大きく影響を与えるという結果を得た。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/ab5dc5a1-9544-4086-acfd-1c672d4e4923.png)

![01_特徴量寄与度.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/83f443b7-9dbf-ba1a-aed4-4ea25541c460.jpeg)

## 分析手法、および、工夫点
* 使用したモデル：LinearRegression(線形回帰)
　深層学習など他のモデルも考えられるが、今回は特徴量の寄与度を明確にすべく線形回帰を選択した。
* 特徴量の追加：相関係数による選択
　モデル評価後、追加特徴量を「モデル誤差と特徴量の相関係数」を基準に選択し、モデル性能の改善を図った。
* ターゲットリーケージの回避：モデル性能評価に与える問題の回避
　目的変数自身を説明変数に混入してしまいモデル性能を過剰に評価してしまうターゲットリーケージ の回避対応を行った。
* 追加特徴量選択方法の妥当性確認：GridSearchの考え方
　モデルの作成/評価指標算出を追加特徴量候補のすべてに行い、特徴量選択方法の妥当性を確認した。


## 分析するデータ
ネットから入手可能な次のデータを元に円ドルの為替レート推論値を算出する。

### 【目的変数】
#### ドル円為替レート
　みずほ銀行　外国為替相場情報　ヒストリカルデータ　公示データ(2002年～)
　https://www.mizuhobank.co.jp/market/historical/index.html
### 【説明変数】
#### 日本の政策金利
　日本銀行 基準割引率および基準貸付利率　 データコード: IR01'MADR1M
　https://www.stat-search.boj.or.jp/#

####  アメリカの政策金利
FEDERAL RESERVE BANK of ST.LOUIS
目標金利は上限と下限が2008-12から設定されている。それ以前は目標値
上限と下限の中間値を目標値として、2008以前のデータとunionする。

　目標金利　上限　https://fred.stlouisfed.org/series/DFEDTARU
　目標金利　下限　https://fred.stlouisfed.org/series/DFEDTARL
　目標金利　～　2008-12-15　https://fred.stlouisfed.org/series/DFEDTAR
#### 通貨の流通量
マクロでみると通貨レートはそれぞれの通貨の流通量の比になるという論がある

　FEDERAL RESERVE BANK of ST.LOUIS
　https://fred.stlouisfed.org/series/BOGMBASE

　日本銀行　時系列統計データ　マネタリーベース平均残高　 データコード: MD01'MABS1AN11
　https://www.stat-search.boj.or.jp/#　　

#### 経常収支
　貿易収支に加え、サービス収支 その他も含まれる。貿易収支よりも実需をより反映していると考えられる。
 　https://www.stat-search.boj.or.jp/ssi/cgi-bin/famecgi2?cgi=$nme_s060
  　キーワード：「収支」 「アメリカ」 にて抽出

## 実行環境
パソコン：LENOBO IdeaPad L3 15ITL6
開発環境：Google Coraboratory
言語：Python
ライブラリ：Pandas, Seaborn, Matplotlib, Sklearn

## 分析の流れ
1. データの収集
2. データの結合、可視化
3. データの標準化、整理
4. 回帰モデルの作成
5. モデルによる推論値の算出
6. モデルの評価
7. 特徴量追加の検討
8. 特徴量追加でのモデル再作成
9. 補足(特徴量選択方法の妥当性検証)

## 1. データの収集
必要なライブラリのインポート
```python
import matplotlib.pyplot as plt
!pip install japanize-matplotlib
import japanize_matplotlib

import pandas as pd
import seaborn as sns
```
### 1.1 ドル円為替レート
```python
#ドル円レートを読み込む
act_rate = pd.read_csv("https://gist.github.com/GreenGrass48/b0accd52a7088e34b9856bddcde76901/raw/c1afa201cea279079b2aa982645d04e50c611a59/quote.csv", encoding="UTF-8")

act_rate = act_rate.iloc[2:,:2]  #不要行、列削除
act_rate.columns = ['date', 'USD_act'] #カラム名設定

#月末のデータに整理する。 df: act_rate

act_rate['date'] = pd.to_datetime(act_rate['date'])
act_rate = act_rate.loc[act_rate.groupby(act_rate['date'].dt.to_period('M'))['date'].idxmax()].reset_index(drop=True)
act_rate['year_mon'] = act_rate['date'].dt.to_period('M')
act_rate=act_rate.iloc[:,1:]
act_rate= act_rate.reindex(columns=['year_mon','USD_act'])
act_rate['USD_act'] = act_rate['USD_act'].astype(float)

act_rate.head()
```
実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/32cfa853-ea5b-6be6-650b-7ed222878548.png)

### 1.2 日本の政策金利
~~~python
policy_rate_JPN = pd.read_csv('https://gist.github.com/GreenGrass48/0d30a856dec5a891772452e6de5a36cc/raw/b56bcfd8ee2c2376bd30ff76a212cd5fa1b1a4c7/nme_R031.1726314.20250109092503.01.csv', encoding='UTF-8')
policy_rate_JPN = policy_rate_JPN.iloc[1:,:]
policy_rate_JPN['データコード'] = pd.to_datetime(policy_rate_JPN['データコード'], format='%Y/%m').dt.to_period('M')
policy_rate_JPN.columns=['year_mon', 'policy_rate_JPN']
policy_rate_JPN['policy_rate_JPN'] = policy_rate_JPN['policy_rate_JPN'].astype(float)
policy_rate_JPN.head()
~~~
実行結果
![12_日本の政策金利.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/4f60159b-710d-10c1-b771-17d14c4d1561.jpeg)
### 1.3 アメリカの政策金利
~~~python
policy_rate_USA_2008 = pd.read_csv('https://gist.github.com/GreenGrass48/46510d1078b64864f490904d7e1cf5cc/raw/c15ddaf16b64b7641dd0c3ec6bb5c1c6ab5bce26/DFEDTAR.csv', encoding='UTF-8')   # FRB 政策金利 2008-12-15まで
policy_rate_USA_upper = pd.read_csv('https://gist.github.com/GreenGrass48/cb538d646cbfa9e201f380e7228d5937/raw/3b9fe30eb7cdd33ee2d207dd433b25d738a66fd3/DFEDTARU.csv', encoding='UTF-8') # FRB 政策金利 2008-12-16から 範囲上限
policy_rate_USA_lower = pd.read_csv('https://gist.github.com/GreenGrass48/f19ea54db1fb0dc3f82ed7206adf35b9/raw/2db064e300e80447cd8faa6396e3bb1d4e312834/DFEDTARL.csv', encoding='UTF-8') # FRB 政策金利 2008-12-16から 範囲下限
policy_rate_USA_avg = pd.merge(policy_rate_USA_lower, policy_rate_USA_upper, on='observation_date', how='inner')  # 上限と下限の中間値を政策金利とする
policy_rate_USA_avg['policy_rate_USA'] = (policy_rate_USA_avg['DFEDTARL'] + policy_rate_USA_avg['DFEDTARU']) / 2

#FRB政策金利 期間別 をまとめる
policy_rate_USA_2008.columns = ['observation_date', 'policy_rate_USA']
policy_rate_USA = pd.concat([policy_rate_USA_2008, policy_rate_USA_avg]).reset_index(drop=True)

#月末のデータに整理する。
policy_rate_USA['observation_date'] = pd.to_datetime(policy_rate_USA['observation_date'])
policy_rate_USA = policy_rate_USA.loc[policy_rate_USA.groupby(policy_rate_USA['observation_date'].dt.to_period('M'))['observation_date'].idxmax()].reset_index(drop=True)

policy_rate_USA = policy_rate_USA.iloc[:,:2]

# year_mon カラムに整理
policy_rate_USA['observation_date'] = policy_rate_USA['observation_date'] .dt.to_period('M')
policy_rate_USA.columns=['year_mon', 'policy_rate_USA']

policy_rate_USA.head()

~~~
実行結果
![13_アメリカの政策金利.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/4469739b-25bb-2917-8eb3-5b49f70d01ef.jpeg)


### 1.4 円の流通量
~~~python
mb_JPN = pd.read_csv('https://gist.github.com/GreenGrass48/08418eb0720a3e0f1b5bf199b4f8b054/raw/83647f0065c0cddd1a3373f9e3748f47aa05e18c/nme_R031.1741542.20250109134108.02.csv', encoding='UTF-8')
mb_JPN =  mb_JPN.iloc[1:,:]
mb_JPN['データコード'] = pd.to_datetime( mb_JPN['データコード'], format='%Y/%m').dt.to_period('M')
mb_JPN.columns=['year_mon', 'mb_JPN']
mb_JPN['mb_JPN'] = mb_JPN['mb_JPN'].astype(float)
mb_JPN.head()
~~~
実行結果
![14_円の流通量.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/0102d394-bfd2-9e0b-4a34-05002b81b0cb.jpeg)

### 1.5 ドルの流通量
~~~python
mb_USA = pd.read_csv('https://gist.github.com/GreenGrass48/f4792094cda3407cb238e1071e00eada/raw/67970f54e9b2e4edb6361931efeb448d5e69e268/BOGMBASE.csv', encoding='UTF-8')   # FRB マネータリーベース

#月末のデータに整理する。
mb_USA['observation_date'] = pd.to_datetime( mb_USA['observation_date'])
mb_USA =  mb_USA.loc[ mb_USA.groupby( mb_USA['observation_date'].dt.to_period('M'))['observation_date'].idxmax()].reset_index(drop=True)

mb_USA =  mb_USA.iloc[:,:2]

# year_mon カラムに整理
mb_USA['observation_date'] =  mb_USA['observation_date'] .dt.to_period('M')
mb_USA.columns=['year_mon', 'mb_USA']

mb_USA.head()

~~~
実行結果
![15_ドルの流通量.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/21c6d02e-b569-f265-c259-a9efe32bdfa1.jpeg)

### 1.6 経常収支
 経常収支は 貿易収支、サービス収支、その他 が含まれる。金融収支も取得したかったが期間が不足するため取得を断念する。
~~~python
balance_data = pd.read_csv('https://gist.github.com/GreenGrass48/cf68378789e53c26af71b24073ba78d9/raw/bd8599ff1032a0848cc2858c0ac85593c8df36bc/nme_R031.2809015.20250114093622.02.csv', encoding='UTF-8')
balance_data =  balance_data.iloc[1:,:]
balance_data['データコード'] = pd.to_datetime( balance_data['データコード'], format='%Y/%m').dt.to_period('M')
balance_data.columns =['year_mon', '経常収支', '貿易・サービス収支', '貿易収支/ネット', '貿易収支/輸出', '貿易収支/輸入', 'サービス収支/ネット', '第一次所得収支/ネット', '第二次所得収支/ネット', '資本移転等収支/ネット', '金融収支/ネット', '（更新停止）経常収支', '（更新停止）貿易・サービス収支', '（更新停止）貿易収支/尻', '（更新停止）貿易収支/輸出', '（更新停止）貿易収支/輸入', '（更新停止）サービス収支/尻', '（更新停止）サービス収支/受取', '（更新停止）サービス収支/支払', '（更新停止）サービス収支/輸送収支/受取', '（更新停止）サービス収支/輸送収支/支払', '（更新停止）サービス収支/輸送収支/海上輸送/受取', '（更新停止）サービス収支/輸送収支/海上輸送/支払', '（更新停止）サービス収支/輸送収支/航空輸送/受取', '（更新停止）サービス収支/輸送収支/航空輸送/支払', '（更新停止）サービス収支/旅行収支/受取', '（更新停止）サービス収支/旅行収支/支払', '（更新停止）サービス収支/その他サービス収支/受取', '（更新停止）サービス収支/その他サービス収支/支払', '（更新停止）サービス収支/その他サービス収支/通信/受取', '（更新停止）サービス収支/その他サービス収支/通信/支払', '（更新停止）サービス収支/その他サービス収支/建設/受取', '（更新停止）サービス収支/その他サービス収支/建設/支払', '（更新停止）サービス収支/その他サービス収支/保険/受取', '（更新停止）サービス収支/その他サービス収支/保険/支払', '（更新停止）サービス収支/その他サービス収支/金融/受取', '（更新停止）サービス収支/その他サービス収支/金融/支払', '（更新停止）サービス収支/その他サービス収支/情報/受取', '（更新停止）サービス収支/その他サービス収支/情報/支払', '（更新停止）サービス収支/その他サービス収支/特許等使用料/受取', '（更新停止）サービス収支/その他サービス収支/特許等使用料/支払', '（更新停止）サービス収支/その他サービス収支/特許等使用料/工業所有権・鉱業権等使用料/受取', '（更新停止）サービス収支/その他サービス収支/特許等使用料/工業所有権・鉱業権等使用料/支払', '（更新停止）サービス収支/その他サービス収支/特許等使用料/著作権等使用料/受取', '（更新停止）サービス収支/その他サービス収支/特許等使用料/著作権等使用料/支払', '（更新停止）サービス収支/その他サービス収支/その他営利業務/受取', '（更新停止）サービス収支/その他サービス収支/その他営利業務/支払', '（更新停止）サービス収支/その他サービス収支/文化・興行/受取', '（更新停止）サービス収支/その他サービス収支/文化・興行/支払', '（更新停止）サービス収支/その他サービス収支/公的その他サービス/受取', '（更新停止）サービス収支/その他サービス収支/公的その他サービス/支払', '（更新停止）所得収支/尻', '（更新停止）所得収支/受取', '（更新停止）所得収支/支払', '（更新停止）所得収支/雇用者報酬/受取', '（更新停止）所得収支/雇用者報酬/支払', '（更新停止）所得収支/直接投資収益/受取', '（更新停止）所得収支/直接投資収益/支払', '（更新停止）所得収支/証券投資収益/受取', '（更新停止）所得収支/証券投資収益/支払', '（更新停止）所得収支/証券投資収益/配当金/受取', '（更新停止）所得収支/証券投資収益/配当金/支払', '（更新停止）所得収支/証券投資収益/債券利子/受取', '（更新停止）所得収支/証券投資収益/債券利子/支払', '（更新停止）所得収支/その他投資収益/受取', '（更新停止）所得収支/その他投資収益/支払', '（更新停止）経常移転収支/尻', '（更新停止）経常移転収支/受取', '（更新停止）経常移転収支/支払', '（更新停止）経常移転収支/労働者送金/受取', '（更新停止）経常移転収支/労働者送金/支払', '（更新停止）資本収支', '（更新停止）投資収支（資産＋負債）', '（更新停止）投資収支（資産）', '（更新停止）投資収支（負債）', '（更新停止）その他資本収支/尻', '（更新停止）その他資本収支/受取', '（更新停止）その他資本収支/支払', '（更新停止）その他資本収支/資本移転/尻', '（更新停止）その他資本収支/資本移転/受取', '（更新停止）その他資本収支/資本移転/支払', '（更新停止）その他資本収支/その他資産/受取', '（更新停止）その他資本収支/その他資産/支払']
balance_data = balance_data[['year_mon', '経常収支', '（更新停止）経常収支']]
balance_data[['経常収支', '（更新停止）経常収支']] = balance_data[['経常収支', '（更新停止）経常収支']] .astype(float).fillna(0)

balance_data['balance'] = balance_data['経常収支'] + balance_data['（更新停止）経常収支']
balance_data = balance_data[balance_data['balance'] != 0 ][['year_mon', 'balance']].reset_index(drop=True)

balance_data.head()
~~~
実行結果
![16_経常収支.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/d3e11137-1d5e-29f5-5f02-fcee49ae4fd0.jpeg)
## 2.データの結合、可視化
### 2.1. 入手したデータを結合
~~~python
from functools import reduce

dfs = [act_rate, policy_rate_JPN, policy_rate_USA, mb_JPN, mb_USA, balance_data]
df = reduce(lambda left, right: pd.merge(left, right, on='year_mon', how='inner'), dfs)

df.head()
~~~
実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/d53d17cc-631b-3c93-d5d8-9b8a653dff6a.png)
## 2.2. データから特徴量を算出
 政策金利差　　      　policy_rate_diff　: policy_rate_USA - policy_rate_JPN  
 通貨流通量の比  　mb_ratio　  　　        : mb_JPN / mb_USA  
 経常収支        　　　　balance  
~~~python
df['policy_rate_diff'] = df['policy_rate_USA'] - df['policy_rate_JPN']
df['mb_ratio'] = df['mb_JPN'] / df['mb_USA'] / 10  # JPN:億円  USA:10憶ドル  単位を合わせるために 10で除算

df.head()
~~~
実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/f57c4d35-9c08-9009-0c1a-e7613a12cf84.png)
### 目的変数と特徴量データを抽出
目的変数と説明変数を持ったDataFrame df_feaures を作成する
~~~python
features = ['year_mon', 'USD_act', 'policy_rate_diff', 'mb_ratio',  'balance']
df_features = df[features]

df_features.head()
~~~
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/e2743820-9d26-5850-2c10-a477bcb771c8.png)

## 2.3.データの可視化
~~~python
#グラフの日本語化用辞書 グラフ表示する際に日本語表示するための辞書を定義
replace_dict = {'USD_act': '実際値', 'policy_rate_diff': '政策金利差', 'mb_ratio': '通貨流通量の比', 'balance': '経常収支', 'trend':'トレンド', 'pred':'推論値'}

#散布図にプロット
df_features_renamed = df_features.rename(columns=replace_dict)

g = sns.PairGrid(df_features_renamed)
g.map(sns.scatterplot)
~~~
実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/20a29c38-e7fa-a6e6-9ace-c63ee0456d17.png)

## 2.4.相関係数マトリクスを確認
　実際値 [USD_act] と最も相関係数が大きいのは 経常収支 [balance]
~~~python
corrmat = df_features.corr(numeric_only=True)

# USD_actとの相関値に基づいて系列をソート
sorted_indices = corrmat.loc['USD_act'].sort_values(ascending=False).index
sorted_corrmat = corrmat.loc[sorted_indices, sorted_indices]
sorted_corrmat = sorted_corrmat.rename(index=replace_dict, columns=replace_dict)

f, ax = plt.subplots(figsize=(5, 5))

sns.heatmap(sorted_corrmat, annot=True, fmt=".2f", cmap="coolwarm", cbar=True)
~~~
実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/e743a11b-4e94-cf5e-06df-c21240f622a9.png)
## 2.5.時系列でのプロット
実際値[USD_act]とそれぞれの説明変数を時系列でプロットして確認する
#### データ準備
時系列情報を持つカラム'year_mon'はPeriod型なのでプロット可能なTimestamp型のカラムを作成する
~~~python
df_features['year_mon_plt'] = df_features['year_mon'].dt.to_timestamp()    #Period型はsnsではpltできないのでplt用にカラムを作成
~~~
#### プロット関数の定義
実際値と各特徴量とを2軸で表示する関数の定義
~~~python
def show_lineplot_vs_USD_with_2nd_axis(title, df, plt_cols_2nd_axis):
    """
    Parameters:
    title: グラフのタイトル
    df: plotするデータを横持しているDataFrame. Period型のカラム'year_mon_plt'が必要
        USD_act が含まれている必要がある
    plot_cols_2nd_axis: plotするカラム名の辞書
       key  : カラム名
       value: グラフに表示するカラー

    Returns:
    lineplotを表示

    """

    import matplotlib.pyplot as plt

    fig, ax1 = plt.subplots(figsize=(10, 4))

    # 実際値の折れ線グラフをプロット
    line1, = ax1.plot(df_features['year_mon_plt'], df_features['USD_act'], color='blue', label='実際値')
    ax1.set_xlabel('年月')
    ax1.set_ylabel('実際値', color='b')
    ax1.tick_params(axis='y', labelcolor='b')

    # 辞書のカラムのY軸を作成
    for plt_col, color in plt_cols_2nd_axis.items():
        ax2 = ax1.twinx()
        line2, = ax2.plot(df_features['year_mon_plt'], df_features[plt_col], color=color, label=replace_dict[plt_col])
        ax2.set_ylabel(replace_dict[plt_col], color=color)
        ax2.tick_params(axis='y', labelcolor=color)

    # 凡例を1か所にまとめる
    lines = [line1, line2]
    labels = [line.get_label() for line in lines]
    fig.legend(lines, labels, loc='upper right', bbox_to_anchor=(0.25, 0.85), ncol=1)

    plt.title(title)
    plt.show()

~~~

USD_実際値と特徴量を比較するためのグラフ表示を行う関数を定義する
表示に際して特徴量の範囲補正の処理も実施する
~~~python
def show_lineplot_vs_USD_act(title, df, plt_cols_adj):
    """
    Parameters:
    title: グラフのタイトル
    df: plotするデータを横持しているDataFrame. Period型のカラム'year_mon_plt'が必要
        USD_act が含まれている必要がある
    plot_cols_adj: plotするカラム名の辞書
       key  : カラム名
       value: 1: USD_actに対して範囲補正を行う
              0: USD_actに対して範囲補正を行わない

    Returns:
    lineplotを表示

    """
    fig, ax1 = plt.subplots(figsize=(10, 4))
    for plt_col, flg in plt_cols_adj.items():
        if flg == 0:
            sns.lineplot(x=df['year_mon_plt'], y=df[plt_col], label=replace_dict[plt_col])
        else:
            y_adj = df[plt_col]*(df['USD_act'].max()-df['USD_act'].min())/(df[plt_col].max()-df[plt_col].min())
            y_shift = df['USD_act'].mean() - y_adj.mean()
            sns.lineplot(x=df['year_mon_plt'], y = y_adj + y_shift, label=replace_dict[plt_col])
    # X,Y軸表示
    plt.xlabel('年月')
    plt.ylabel('USDレート')
    # タイトルを追加
    plt.title(title)
    #凡例を表示
    plt.legend()

    # グラフを表示
    plt.show()
~~~

#### 経常収支
~~~python
title = 'vs_経常収支'
plt_cols_2nd_axis = {'balance':'green'}

show_lineplot_vs_USD_with_2nd_axis(title, df_features, plt_cols_2nd_axis)
~~~
実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/e541ca8b-c4f0-c7b2-241a-1825e5fe893e.png)


#### 政策金利差 
~~~python
title = 'vs 政策金利差'
# plt_cols_adj = {'USD_act': 0 , 'policy_rate_diff': 1}
plt_cols_2nd_axis = {'policy_rate_diff':'orange'}

show_lineplot_vs_USD_with_2nd_axis(title, df_features, plt_cols_2nd_axis)
~~~
実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/6f0bf82d-2f90-480f-73b1-cfd6a47e3b83.png)

#### 通貨流通量の比 
~~~python
title = 'vs 通貨流通量の比'
plt_cols_2nd_axis = {'mb_ratio':'red'}

show_lineplot_vs_USD_with_2nd_axis(title, df_features, plt_cols_2nd_axis)
~~~
実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/72f51bb9-06fc-132e-db10-2ee75f2eb138.png)


#### 全ての特徴量 
~~~python
title = 'vs 特徴量'
plt_cols_adj = {'USD_act': 0 , 'policy_rate_diff': 1, 'balance': 1, 'mb_ratio': 1}

show_lineplot_vs_USD_act(title, df_features, plt_cols_adj)
~~~
実行結果
特徴量を重ねてみると経常収支だけで推論値が導けそうだが、実際に３つの特徴量で回帰分析を実行してその効果を定量的に評価する。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/9253f389-687e-be96-5e23-0ab991deb4b8.png)


## 3.データの標準化,整理
特徴量の単位がそれぞれ異なるのでデータのスケーリングが必要である。スケーリング方法として標準化と正規化があるが、回帰モデルを適用するので標準化を適用する。
参考url: https://qiita.com/oki_kosuke/items/02ec7bee92c85cf10ac2
#### df_feaures から説明変数 X と目的変数 y を抽出
~~~python
X_columns = ['policy_rate_diff','mb_ratio' ,'balance']
y_columns = ['USD_act']

# 説明変数Xと目的変数yへ分離
X = df_features[X_columns].to_numpy()
y = df_features[y_columns].to_numpy()
~~~
今回は算出した推論値と実績値を比較することが目的なのでテストデータの分離は行わず、現在までのデータすべてを用いてモデルを作成する。
~~~python
X_train = X
y_train = y
~~~

### スケーリングの実行

~~~python
# クラスのインポート
from sklearn.preprocessing import StandardScaler
# クラスのインスタンス化
scaler = StandardScaler()
# X_trainを標準化する変換モデルを生成
scaler.fit(X_train)
~~~
実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/dbdcd63b-7772-d828-6cd5-63f3a26374ce.png)

~~~python
# X_trainのスケール変換
X_train_scaled = scaler.transform(X_train)
# スケール後のdf作成、確認
df_X_train_scaled = pd.DataFrame(X_train_scaled, columns=X_columns)
# スケール後のdf_X_train_scaledは平均 0 標準偏差 1になっている。（＝標準化できている）
df_X_train_scaled.describe()
~~~
実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/ad992f24-c812-9547-f21d-08feb9ce114c.png)

## 4.回帰モデルの作成
~~~python
from sklearn.linear_model import LinearRegression
model = LinearRegression()
# モデルの学習
model.fit(X_train_scaled, y_train)
~~~
実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/cea43d0a-adb1-efd0-17de-98798a50984a.png)

#### 重回帰の計算式の確認 
~~~python
print(model.coef_)
print(model.intercept_)
~~~
実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/0df87983-8bcc-1c88-b5a6-93e316f1065d.png)

### 回帰係数の可視化
~~~python
# seabornによる 回帰係数の可視化
plt.title('特徴量の寄与度')
sns.barplot( x=[replace_dict[key] for key in X_columns], y=model.coef_[0])
~~~
実行結果  実際値との相関係数が最も高かった経常収支の寄与度が想定通り大きい
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/e3eda63b-718a-51ae-225f-8d9e69668388.png)

## 5.モデルによる推論値算出
### 5.1. 推論値算出
~~~python
# 学習させたデータでの推論の実行
pred = model.predict(X_train_scaled)
~~~

### 5.2.推論値と実績値の比較
推論値と実際値の比較するため、推論値を結合
~~~python
#推論値と実際値の比較するため、推論値を結合
df_pred = pd.DataFrame(pred)
df_pred.columns = ['pred']

df_features_pred = pd.concat([df_features, df_pred], axis=1)
# df_features_pred.head()
~~~
推論値と実際値の比較を可視化
### 
~~~python
title = '実際値 vs 推論値'
plt_cols_adj = {'USD_act': 0 , 'pred': 1}

show_lineplot_vs_USD_act(title, df_features_pred, plt_cols_adj)
~~~
実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/e407ff86-4923-cef0-4a60-72dd12ab239e.png)


## 6.モデルの評価  
 算出した推論値がどれくらい当てはまりがあるのか定量的に評価する。  
 単純に可視化した [政策金利差と実際値] と [推論値と実際値]は、見た感じではどちらが当てはまっているのかよくわからないので、定量的な評価を行う。  

 評価方法：平均二乗誤差（Mean Squared Error, MSE）を指標とする。小さい方が良い。

政策金利差は実際値と単位が異なるので単純にMSEの算出ができない。そのためそれぞれのカラムを標準化する
### 6.1.前処理 標準化の実施
~~~python
# クラスのインポート
from sklearn.preprocessing import StandardScaler
# クラスのインスタンス化
scaler = StandardScaler()

mse_columns = ['pred', 'balance', 'USD_act']
mse_pred = df_features_pred[mse_columns].to_numpy()
# X_trainを標準化する変換モデルを生成
scaler.fit(mse_pred)
~~~
実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/460b71e0-0342-f302-20b2-7524caadae24.png)

### 
~~~python
# X_trainのスケール変換
mse_pred_scaled = scaler.transform(mse_pred)
# スケール後のdf作成、確認
df_mse_pred_scaled = pd.DataFrame(mse_pred_scaled, columns=mse_columns)
~~~

### 6.2. MSEを算出
~~~python
# MSEの計算　pred vs USD_act
mse = ((df_mse_pred_scaled['pred'] - df_mse_pred_scaled['USD_act']) ** 2).mean()
print(f'MSE: 推論値     vs 実際値: {mse}')
# MSEの計算　balance vs USD_act
mse = ((df_mse_pred_scaled['balance'] - df_mse_pred_scaled['USD_act']) ** 2).mean()
print(f'MSE: 政策金利差 vs 実際値: {mse}')
~~~
実行結果 
　推論値のMSEは 0.285 で 政策金利差のMSE 0.420 に比べて小さいことがわかる。
　すなわち**推論値は政策金利差より実際値に近い**と判断でき、**回帰分析の有効性が確認**できた。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/201df907-170b-48cf-ce65-a46b28125faf.png)

## 7.特徴量追加の検討
ここまで３つの特徴量にてドル円為替レートの推論値算出を行った。実際の為替レートは市場心理を表す市場トレンドも影響していると考えられる。市場トレンドは短期移動平均と長期移動平均の差で計算し、その差がプラスであれば上昇トレンド、マイナスであれば下降トレンドを示す。これを「トレンド」として追加の特徴量とするかの検討を行う。

トレンドは、差をとる移動平均の期間の組み合わせがパラメータとして存在する。短期と長期のどの組み合わせが最も効果的であるかを、モデル誤差(実勢値と推論値の差)との相関を基準に選択する。

トレンド算出する際、目的変数である現時点の実際値を組み込んでしまうとターゲットリーケージが発生する恐れがある。これを回避するため当月の実際値を前月の実測値に置き換えることで移動平均算出に当月の実測値を含まず前月以前の実際値を使うこととする。

:::note info
ターゲットリーケージ
モデルのトレーニング時に目的変数の情報が直接的または間接的に説明変数に含まれることで、モデルの性能を過剰に評価してしまう問題のこと
:::

### 7.1移動平均データの作成
ターゲットリーケージ回避のために今月データを前月データにシフトする
~~~python
trend = act_rate.copy()
# ターゲットリーケージを回避するために今月データを前月データにシフトする

trend['USD_act'] = trend['USD_act'].shift(1)
trend = trend.dropna()
~~~
期間毎の組み合わせ移動平均データを生成する
~~~python
#期間毎の組み合わせ移動平均データの生成する
num_ma =[2,3,6,9,12]  #期間の単位は月なので 移動平均期間は60日、90日、180日、270日、360日

for n in num_ma:
  trend[f'tmp_{n}ma'] = trend['USD_act'].rolling(window=n).mean()

# print(trend.head(10))
import itertools
# 2つの要素の組み合わせをループで処理
for combo in itertools.combinations(num_ma, 2):
    # print(combo)
    trend[f'trend_{combo[0]}ma_{combo[1]}ma'] = trend[f'tmp_{combo[0]}ma'] - trend[f'tmp_{combo[1]}ma']

trend = trend.dropna()

# 移動平均の差のカラムのみに整理　'tmp'を含むカラムを削除
columns_to_drop = [col for col in trend.columns if 'tmp' in col]
trend = trend.drop(columns=columns_to_drop).reset_index(drop=True)

trend.head()
~~~
実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/079f25c9-23b7-28da-6ca9-d6d3a7b77e8e.png)

### 7.2 モデル誤差との相関算出
~~~python
#モデル誤差の算出
df_features_pred['diff_USD_pred_act'] = df_features_pred['pred'] - df_features_pred['USD_act']
df_features_pred.head()
~~~
実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/e04a5b71-d8b9-9a95-8e5d-a42a4626916b.png)

~~~python
#相関係数算出のためのdfのカラム定義
col_len = len(trend.columns)
plt_cols = trend.columns[2:col_len-1]
cmp_cols_ls = ['year_mon', 'diff_USD_pred_act']+list(plt_cols)
cmp_cols_ls
~~~
実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/aa274738-51c7-23bb-eddc-99bda24704da.png)


~~~python
#相関係数算出のためのdf作成
df_cmp = pd.merge(df_features_pred, trend, on='year_mon', how='inner')[cmp_cols_ls]
df_cmp.head()
~~~
実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/16fc7204-163b-1236-1e60-1a16033a9fd9.png)

~~~python
#相関係数を算出
corr = df_cmp.corr(numeric_only=True)
corr
~~~
実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/0d5f2518-006a-4cdf-7761-1159069c635d.png)

trendtrend_2ma_6ma(60日移動平均と180日移動平均から算出したトレンド)の相関係数の絶対値が最も大きく0.4を超えている。強い相関とは言えないが無相関ではないので
**trend_2ma_6maを新たな特徴量'trend'として追加**する。
### 追加特徴量 trend の作成
~~~python
trend = trend[['year_mon', 'trend_2ma_6ma']].reset_index(drop=True)
trend.columns = ['year_mon', 'trend']
trend.head()
~~~
実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/0e28db4a-c428-9367-f82e-c650b01e6fab.png)

## 8.特徴量追加でのモデル再作成 
### 8.1. データの結合
特徴量データに trendを追加
~~~python
df_features_2 = pd.merge(df_features, trend, on='year_mon', how='inner')
features_2 = ['year_mon', 'USD_act', 'policy_rate_diff', 'mb_ratio',  'balance', 'trend']  #trend追加
df_features_2 = df_features_2[features_2]

df_features_2.head()
~~~
実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/06d37230-d9e7-ae84-4d34-0e57d5550725.png)

### データの可視化
散布図
~~~python
#散布図にプロット
import seaborn as sns
df_features_2_renamed = df_features_2.rename(columns=replace_dict)
g = sns.PairGrid(df_features_2_renamed)
g.map(sns.scatterplot)
~~~
実行結果
![pairgrid_plot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/118b1dcc-10f5-8051-9c86-b99c58f265e6.png)

相関マトリクス
~~~python
import matplotlib.pyplot as plt

corrmat = df_features_2.corr(numeric_only=True)

# USD_actとの相関値に基づいて系列をソート
sorted_indices = corrmat.loc['USD_act'].sort_values(ascending=False).index
sorted_corrmat = corrmat.loc[sorted_indices, sorted_indices]
# インデックスを辞書に基づいて置き換える
sorted_corrmat = sorted_corrmat.rename(index=replace_dict,  columns=replace_dict)


f, ax = plt.subplots(figsize=(5, 5))

sns.heatmap(sorted_corrmat, annot=True, fmt=".2f", cmap="coolwarm", cbar=True)
~~~
実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/4adec992-4601-45fb-1fe6-e477f29b5947.png)

### 8.2.データのスケーリング
df_feaures からXとyを抽出
~~~python
X_columns = ['policy_rate_diff','mb_ratio' ,'balance', 'trend']  #trend追加
y_columns = ['USD_act']
# 説明変数Xと目的変数yへ分離
X = df_features_2[X_columns].to_numpy()
y = df_features_2[y_columns].to_numpy()

X_train = X
y_train = y

# クラスのインスタンス化
scaler = StandardScaler()
# X_trainを標準化する変換モデルを生成
scaler.fit(X_train)
# X_trainのスケール変換
X_train_scaled = scaler.transform(X_train)
# スケール後のdf作成
df_X_train_scaled = pd.DataFrame(X_train_scaled, columns=X_columns)
~~~

### 8.3.回帰モデルの作成

~~~python
from sklearn.linear_model import LinearRegression
model = LinearRegression()
# モデルの学習
model.fit(X_train_scaled, y_train)
# 予測モデルの評価
# 学習データの決定係数
model.score(X_train_scaled, y_train)
# 重回帰の計算式の確認
print(model.coef_)
print(model.intercept_)
~~~
実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/26681d32-8b43-8069-dca4-bc9fb79bdc1e.png)

### 回帰係数の可視化

~~~python
# seabornによる 回帰係数の可視化
plt.title('特徴量の寄与度')
sns.barplot( x=[replace_dict[key] for key in X_columns], y=model.coef_[0])
~~~
実行結果

 
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/55f08160-3594-34b2-1acf-a5a7b3953692.png)
### 8.4.モデルによる推論値算出
推論値算出
~~~python
# 推論値算出
pred = model.predict(X_train_scaled)
~~~
推論値と実績値の比較
~~~python
#推論値と実績値の比較
df_pred = pd.DataFrame(pred)
df_pred.columns = ['pred']

df_features_2_pred = pd.concat([df_features_2, df_pred], axis=1)

#グラフ表示
title = '実際値 vs 推論値(トレンド追加)'
plt_cols_adj = {'USD_act': 0 , 'pred': 1}
df_features_2_pred['year_mon_plt'] = df_features_2_pred['year_mon'].dt.to_timestamp()    #Period型はsnsではpltできないのでplt用にカラムを作成

show_lineplot_vs_USD_act(title, df_features_2_pred, plt_cols_adj)
~~~
実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/49222db9-7e68-2ad2-35f5-6db6e1042ede.png)


### 8.5.モデルの評価
前処理 標準化の実施
~~~python
# クラスのインポート
from sklearn.preprocessing import StandardScaler
# クラスのインスタンス化
scaler = StandardScaler()

mse_columns = ['pred', 'balance', 'USD_act', 'trend']  #trend追加
mse_pred = df_features_2_pred[mse_columns].to_numpy()
# X_trainを標準化する変換モデルを生成
scaler.fit(mse_pred)
# X_trainのスケール変換
mse_pred_scaled = scaler.transform(mse_pred)
# スケール後のdf作成
df_mse_pred_scaled = pd.DataFrame(mse_pred_scaled, columns=mse_columns)
~~~
MSEを算出
~~~python
# 予測値と実績値のmse算出

# MSEの計算　pred vs USD_act
mse = ((df_mse_pred_scaled['pred'] - df_mse_pred_scaled['USD_act']) ** 2).mean()
print(f'MSE: 推論値     vs 実際値: {mse}')
# MSEの計算　balance vs USD_act
mse = ((df_mse_pred_scaled['balance'] - df_mse_pred_scaled['USD_act']) ** 2).mean()
print(f'MSE: 政策金利差 vs 実際値: {mse}')
~~~
実行結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/382b27c3-ec42-4b8b-7680-de8cab3518bf.png)
参考：balanceに対するMSEが、trend 追加前の0.42と異なるのは、以下の理由による。
　　　　trend 追加時に移動平均を算出した。この時 移動平均が算出されない行をdropしたため。
### 8.6.結果
trendを考慮しなかった場合のMSE:0.285 から trendを追加することで MSE:0.203に改善することができた。
**トレンドの影響を考慮することで推論値の改善が図れた**ことが確認できた。


## 補足：トレンド選択方法の妥当性検証
 7.特徴量追加検討で相関係数のもっとも高い移動平均期間組み合わせを選択したが、その妥当性をGridSearch的に総当たりして確認する。
### 線形回帰分析関数 
 総当たりを実行するにあたりMSE算出を関数化する。
~~~python
# クラスのインポート
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression

def execute_linear_regression(df_features, X_columns, y_columns, trend_col):
    """
    Parameters:
    df_features : 目的変数を含んだ分析対象のdf
    X_columns   : 説明変数カラム名
    y_columns   : 目的変数カラム名
    trend_col   : トレンドカラム名

    Returns:
    mse        : 実際値と推論値のMSE
    """

    # 説明変数Xと目的変数yへ分離
    X = df_features[X_columns].to_numpy()
    y = df_features[y_columns].to_numpy()
    X_train = X
    y_train = y

    # クラスのインスタンス化
    scaler = StandardScaler()
    # X_trainを標準化する変換モデルを生成
    scaler.fit(X_train)
    # X_trainのスケール変換
    X_train_scaled = scaler.transform(X_train)
    # スケール後のdf作成
    df_X_train_scaled = pd.DataFrame(X_train_scaled, columns=X_columns)

    #8.3.回帰モデルの作成
    model = LinearRegression()
    # モデルの学習
    model.fit(X_train_scaled, y_train)

    # # 重回帰の計算式の確認
    # print(model.coef_)
    # print(model.intercept_)

    #8.4.モデルによる推論値算出
    # 学習させたデータでの予測の実行
    pred = model.predict(X_train_scaled)
    #予測値と実績値の比較
    df_pred = pd.DataFrame(pred)
    df_pred.columns = ['pred']

    df_features_pred = pd.concat([df_features, df_pred], axis=1)

    # 8.5.モデルの評価
    scaler = StandardScaler()

    mse_columns = ['pred', 'balance', 'USD_act', trend_col]  #trend追加
    mse_pred = df_features_pred[mse_columns].to_numpy()
    # X_trainを標準化する変換モデルを生成
    scaler.fit(mse_pred)

    # X_trainのスケール変換
    mse_pred_scaled = scaler.transform(mse_pred)
    # スケール後のdf作成
    df_mse_pred_scaled = pd.DataFrame(mse_pred_scaled, columns=mse_columns)

    # 予測値と実績値のmse算出

    # MSEの計算　pred vs USD_act
    mse = ((df_mse_pred_scaled['pred'] - df_mse_pred_scaled['USD_act']) ** 2).mean()

    return mse
~~~


### 総当たりの実行
差分を取る移動平均期間の組み合わせを変化させたトレンドデータ各種にて線形回帰を実行してMSEを一覧化する。

~~~python
# トレンドデータ各種を特徴量として一つのDataFrameに統合
features_base = ['year_mon', 'USD_act', 'policy_rate_diff', 'mb_ratio', 'balance']
df_features_3 = df_features_2[features_base]
df_features_3 = pd.merge(df_features_3, df_cmp, on='year_mon', how='inner')

#df_features_3.head()
~~~

~~~python
# GridSearch的総当たりを実行して結果をTableに表示

from prettytable import PrettyTable

# prettyTableのインスタンス作成、見出し設定
table = PrettyTable()
table.field_names = ['トレンド期間', "MSE"]

# 説明変数のトレンドとしてトレンドデータ各種を入れ替えてMSEを算出、一覧化

for i in range(6,len(df_features_3.columns)):
    trend_col = df_features_3.columns[i]

    X_columns = ['policy_rate_diff','mb_ratio' ,'balance', trend_col]
    y_columns = ['USD_act']

    mse = execute_linear_regression(df_features_3, X_columns, y_columns, trend_col)
    table.add_row([trend_col, format(mse,'.3f')])

print(table)
~~~
実行結果
![スクリーンショット 2025-01-27 143431.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782092/5400932d-4d19-b40f-fc62-1c0d140978b8.png)


### 結果
実際に各種トレンドデータを用いて算出したMSEの最小は 0.203のtrend_2ma_6maであった。これより**トレンド検討での移動平均期間組み合わせの選択の妥当性が確認できた。**

## まとめ
### 分析したデータから得られた情報
　ドル円為替レートは昨今、日米のそれぞれの中央銀行の政策に焦点が当たり、金利の差に大きく影響を受けていることがニュースになっている。 実際に収集できるデータを長期にわたり分析してみると、貿易収支を含む経常収支の影響が大きいことがわかった。 実需といわれる部分で実際の貨幣交換によって為替レートが影響を受けるという常識的な結果となった。また、それに次いで市場トレンドの影響が大きく政策金利差の影響は最も小さいという結果となった。データ分析を行うことで新たな発見を得ることができた。
### 分析の手法について
　目的変数に関連ありそうな説明変数のデータを収集･加工してモデルの作成を行うだけでなく、その目的変数に対する知識から新たな特徴量を算出して加えることでモデルの改善を進めることができた。
 　分析技術はもちろんデータに対する知見も重要であることを改めて体験することができた。