Python初心者が主成分分析を用いてビットコインのファクター分析をやってみた

#はじめに
自己紹介：
金融機関にてEXCELによるデータ分析実務○年
⇒Aidemy（データ分析）を受講
⇒Pythonによるデータ分析を開始

EXCEL大好きマンとしてデータ分析をEXCELで突き通してきたのですが、さすがに限界を感じるのと、Pythonぐらい使えないと時代に取り残されて野垂れ死にしそうだったので、一念発起しスクールを受講しました。
Python戦士Lv.5程度にはなれた…気がするので、腕試しとして興味のある分析をPythonで実施しました。本記事はその際の記録を目的としており、Pythonの技術面では以下の点に触れています。

• investpyのライブラリを使用した金融資産・暗号資産の価格データ取得
• 価格データから収益率データを作成して標準化するフロー
• scikit-learnによる主成分分析（繰り返しによる時系列分析）
• matplotlib.pyplotを利用したグラフ図示（積み上げ棒グラフや折れ線グラフを並べて表示する 等）

前置きが長くなりましたが、以下より本題です↓

#先行研究の紹介
資産配分を決定する上で、各資産クラスの共通変動要因（ファクター）に着目するアプローチが近年注目を集めている。
諸先輩方の数々の先行研究にて、主成分分析により得られたマルチアセット市場におけるファクターの意味付けを定性的・定量的に行った結果、上位3ファクターが以下のように解釈されている¹。

1. 経済成長

経済成長に対して敏感に反応すると考えられるリスク資産、株式価格等の変動に寄与するファクター

2. 実質金利

金利に敏感に反応すると考えられる資産、債券価格等の変動に寄与するファクター

3. インフレーション

物価変動に敏感に反応すると考えられる資産、主にコモディティ等の商品価格の変動に寄与するファクター

#本記事での分析概要と結果
先行研究¹では、以下のような金融資産の収益率に対して主成分分析を行い、前述の3ファクターを抽出している。

• 株式
• 国債
• REIT（不動産投資信託）
• ハイイールド社債
• コモディティ（原油や大豆等の商品）
• 物価連動国債

本記事では、これらの資産クラスに暗号資産を加えて主成分分析を行う。
近年、暗号資産はポートフォリオ分散の観点から、機関投資家でも投資を開始している先が出てきている。従来の投資先である金融資産との比較により、暗号資産の価格変動に寄与するファクターを調査することで、リスクファクターの観点から**暗号資産がポートフォリオ分散に寄与するのか？**という仮説を検証したい。 決して暗号資産に投資して億り人になりたい訳ではない。

結論を先に述べておくと、結果として暗号資産には金融資産とは関係のないファクターが価格変動に寄与していることが分かった。これは、金融資産の価格変動を説明する経済成長・実質金利・インフレとは異なる要因で、暗号資産の価格が変動していることを示しており、ポートフォリオ分散に寄与する可能性を示唆する結果と言える。

以降の章では、以下のような流れで詳細な分析内容を紹介していく。
尚、本分析の実行環境にはGoogle Colaboratoryを利用した。

• 各資産クラスの価格データ取得
• データの前処理
• 主成分分析の実施
• 結果の解釈

#各資産クラスの価格データ取得
各種資産価格データ取得のためにinvestpyを利用する。
Google Colaboratoryにinvestpyはデフォルトで含まれていないため、最初にインストールする必要がある。

cryptopca.ipynb

pip install investpy # データ取得のための準備（最初に実施）

今回の分析に必要なライブラリを設定。

cryptopca.ipynb

import pandas as pd #データ処理(データ形式の処理等)
import numpy as np #価格データ処理(数値計算)
import datetime #日付データ処理
import investpy #資産価格データ取得
import matplotlib.pyplot as plt# 図やグラフの図示
%matplotlib inline

各種資産価格のデータは、暗号資産、指数、ETFのデータより以下を取得した（指数だけで取得できるデータに限界があり取得したい指数にトラックするETFも利用）。
各資産を並べてみると、ビットコインのリスクが異常なほど高く、いかにハイリスクな投資であるかを物語っている。。

資産クラス	分析対象	Investpyデータソース	リスク（年率標準偏差）
暗号資産	ビットコイン	investpy.crypto	188.8%
米国株式	S&P 500	investpy.indices	13.4%
コモディティ	S&P GSCI Commodity	investpy.indices	22.0%
米国REIT	FTSE EPRA/NAREIT USA	investpy.indices	16.7%
米ドル建てハイイールド社債	iShares iBoxx $ High Yield Corporate Bond	investpy.etf	6.6%
米国債	iShares Core US Treasury Bond	investpy.etf	3.8%
米国物価連動国債	iShares TIPS Bond	investpy.etf	4.4%

価格データの取得に向け、最初に期間の設定とデータ格納用の空のデータフレームを作成。本分析では、米国債のデータ開始期間に合わせ、2012年12月～2021年5月の9年弱を分析期間とし、価格は為替の影響を受けないようにUSD建てとした。

cryptopca.ipynb

# 各資産クラスの価格取得（対USD）
# 日付設定
from_date = '01/12/2012' #取得開始日付
to_date = datetime.date.today().strftime('%d/%m/%Y')#取得終了日付(当該コード上では本日日付を利用)
price_df = pd.DataFrame(index=[], columns=[]) # 価格データ格納用の空のデータフレームを作成

investpyでは、データ形式や資産クラスによって関数が異なるため、各資産クラスに応じてコードを微妙に変更する必要がある。
最初に、暗号資産を取得するためにinvestpyの暗号資産レファレンスを参照する。

cryptopca.ipynb

investpy.crypto.get_cryptos()#用意されている全通貨データを確認

暗号資産の価格を取得するコード↓

cryptopca.ipynb

# 暗号資産の価格を取得するコード
#cc_lst = ['Bitcoin' ,'Ethereum' ,'XRP'] #①Bitcoin、②Ethereum、③XRP←他通貨を選択する場合
cc_lst = ['Bitcoin'] #①Bitcoinのみ選択する場合（本分析はこちら）
for crypto in cc_lst:
    cc_df = investpy.crypto.get_crypto_historical_data(crypto, from_date, to_date, as_json=False, order='ascending', interval='Monthly')
    cc_df.reset_index(drop=False, inplace=True)#データフレームのインデックスを再定義
    cc_df = cc_df.drop(['Open', 'High', 'Low', 'Volume', 'Currency'], axis=1) #終値だけ取得
    cc_df_new = cc_df.rename(columns={'Close': crypto})#名前変更
    
    if crypto == cc_lst[0]:
       price_df = price_df.append(cc_df_new, ignore_index=True)#データの追加(最初だけ日付データが入るように設定)
    
    else:
       price_df = pd.concat([price_df, cc_df_new[crypto]], axis=1)#データの追加(価格データのみ抽出)

次に、指数データより金融資産価格データを取得するためにinvestpyの指数レファレンスを参照する。

cryptopca.ipynb

investpy.indices.get_indices(country=None)#用意されている全指数データを確認

指数データより取得できる資産クラスの価格を取得するコード↓

cryptopca.ipynb

# 指数データより取得できる資産クラスの価格を取得
i_lst = ['S&P 500', 'S&P GSCI Commodity TR', 'FTSE EPRA/NAREIT USA'] #④米国株式、⑤コモディティ、⑥米国REIT
for index in i_lst:
    
    if index == i_lst[1]:
       i_country = 'world'#コモディティのみ国指定がworldのため分岐を用意
    else:
       i_country = 'united states'

    i_df = investpy.indices.get_index_historical_data(index, i_country, from_date, to_date, as_json=False, order='ascending', interval='Monthly')
    i_df.reset_index(drop=False, inplace=True)#データフレームのインデックスを再定義
    i_df = i_df.drop(['Open', 'High', 'Low', 'Volume', 'Currency'], axis=1) #終値だけ取得
    i_df_new = i_df.rename(columns={'Close': index})#名前変更
    price_df = pd.concat([price_df, i_df_new[index]], axis=1)##データの追加(価格データのみ抽出)

最後に、ETFデータより金融資産価格データを取得するためにinvestpyのETFレファレンスを参照する。

cryptopca.ipynb

investpy.etfs.get_etfs(country=None)#用意されている全ETFデータを確認

ETFデータより取得できる資産クラスの価格を取得するコード↓

cryptopca.ipynb

# ETFデータより取得できる資産クラスの価格を取得
etf_lst = ['iShares iBoxx $ High Yield Corporate Bond', 'iShares Core US Treasury Bond', 'iShares TIPS Bond'] #⑦ハイイールド社債、⑧米国債、⑨インフレ連動債
etf_country = 'united states'
for etf in etf_lst:
    etf_df = investpy.etfs.get_etf_historical_data(etf, etf_country, from_date, to_date, stock_exchange=None, as_json=False, order='ascending', interval='Monthly')
    etf_df.reset_index(drop=False, inplace=True)#データフレームのインデックスを再定義
    etf_df = etf_df.drop(['Open', 'High', 'Low', 'Volume', 'Currency', 'Exchange'], axis=1) #終値だけ取得
    etf_df_new = etf_df.rename(columns={'Close': etf})#名前変更
    price_df = pd.concat([price_df, etf_df_new[etf]], axis=1)#データの追加(価格データのみ抽出)

想定する全資産クラスの価格データが取得できたことを以下のコードで確認する。

cryptopca.ipynb

display(price_df.tail()) #データの確認

結果：

#データの前処理
最初に、分析のために価格データを収益率データに変換。

cryptopca.ipynb

#価格⇒収益率データに変更
return_df = pd.DataFrame(index=[], columns=[])
for l in list(price_df.columns):
    if l == 'Date':
       return_df[l] = price_df[l]       
    else:
       return_df[l] = price_df[l].pct_change().dropna()

return_df = return_df.dropna(how='any') #欠損値を除外
r_df = return_df.set_index('Date') #日付をindex化したデータに変換(名前を変更)
display(r_df.tail())#データの確認
display(r_df.corr())#相関係数の確認
display(r_df.std())#標準偏差の確認

結果：

続いて、データを標準化する処理を行う（標準化を行わないとビットコインの異常に高い標準偏差に結果が引きずられてしまう）。

cryptopca.ipynb

#データの前処理として標準化を実施
r_df = (r_df - r_df.mean(axis=0)) / r_df.std(axis=0)# 各資産の収益率データについて平均が0、分散・標準偏差が1となるように変換
display(r_df.tail())#データの確認
display(r_df.corr())#相関係数の確認
display(r_df.std())#標準偏差の確認

結果：処理前と比較して相関係数に変化はなく標準偏差が全て1.0になっていることを確認

#主成分分析の実施

主成分分析にはscikit-learnを利用する（尚、Aidemyでは主成分分析の理解のためにPCA関数を自作した）。
9年弱のデータを8年間のローリングで分析を実施し、主成分の時系列変化を追えるようにすることで、結果の解釈をし易いようにしている（極力長期で最低限ローリングできる期間として8年間を採用）。

cryptopca.ipynb

from sklearn.decomposition import PCA #scikit-learnより主成分分析用のライブラリを取得

# 主成分分析のインスタンスを生成
pca = PCA()
c_rate = pd.DataFrame(index=[], columns=[])
eig_vec = pd.DataFrame(index=[], columns=[])

# データから変換モデルを学習し変換する。
t = 8 #分析期間(年間)
for i in range(len(r_df)-t*12):
    X=r_df[i:i+t * 12]
    X_pca = pca.fit_transform(X)
    c_rate = pd.concat([c_rate, pd.DataFrame(pca.explained_variance_ratio_, index=["PC{}".format(x + 1) for x in range(len(X.columns))])], axis=1)# 寄与率の計算
    eig_vec = pd.concat([eig_vec, pd.DataFrame(pca.components_.T, index = X.columns, columns = ["PC{}".format(x + 1) for x in range(len(X.columns))])], axis=1)# 固有ベクトル

c_rate = c_rate.T.set_index(price_df['Date'][(len(r_df)-t*12+1)*-1:-1]) # 寄与度データの転置とIndexの日付化
pc1 = eig_vec['PC1'].T.set_index(price_df['Date'][(len(r_df)-t*12+1)*-1:-1])*-1# PC1の転置とIndexの日付化、符号の反転
pc2 = eig_vec['PC2'].T.set_index(price_df['Date'][(len(r_df)-t*12+1)*-1:-1])# PC2の転置とIndexの日付化
pc3 = eig_vec['PC3'].T.set_index(price_df['Date'][(len(r_df)-t*12+1)*-1:-1])# PC3の転置とIndexの日付化
pc4 = eig_vec['PC4'].T.set_index(price_df['Date'][(len(r_df)-t*12+1)*-1:-1])# PC3の転置とIndexの日付化

#結果の解釈
暗号資産を①「加えない場合」と②「加えた場合」の両方の結果を比較することで、結果の解釈を行う。
前者の解釈については、先行研究の解釈をそのまま利用することとする。

####分析結果を図示する
Pythonのmatplotlib.pyplotを利用して、分析結果を図示してみる。

cryptopca.ipynb

# 寄与度の積み上げ棒グラフ表示(暗号資産含み)
fig, ax = plt.subplots(2, 3, figsize=(15, 10), tight_layout=True, facecolor='whitesmoke')
c_rate.plot(kind='bar', stacked=True, ax=ax[0,0], title='Cumulative Contirbution')

# PCの棒グラフ表示(暗号資産含み)
pc1.plot(kind='line', stacked=False, ax=ax[0,1], title='PC1')
pc2.plot(kind='line', stacked=False, ax=ax[0,2], title='PC2')
pc3.plot(kind='line', stacked=False, ax=ax[1,0], title='PC3')
pc4.plot(kind='line', stacked=False, ax=ax[1,1], title='PC4')
plt.show()

####暗号資産除きのデータを対象とした主成分分析

主成分分析のグラフを確認すると、PC1～PC3の累積寄与度が9割超となり、各資産クラスの価格変動の大半をPC3までで説明できた。また、各主成分の時系列の傾向も以下のように先行研究¹と似たような解釈が可能となった。

PC1（≒経済成長）

株式やREIT、ハイイールド社債といったリスク性資産が相対的に高いプラスで推移する一方で、国債がマイナスで推移していることから、経済成長を示すファクターであることが伺える

　

PC2（≒実質金利）

国債や物価連動国債が相対的に高いプラスで推移する一方で、その他資産がマイナスないしは低位なプラス水準で推移していることから、実質金利を示すファクターであることが伺える

　

PC3（≒インフレーション）

コモディティや物価連動国債が相対的に高いプラスで推移する一方で、その他資産がマイナスないしは低位なプラス水準で推移していることから、インフレーションを示すファクターであることが伺える

####暗号資産含みのデータを対象とした主成分分析

暗号資産であるビットコインを含んだ主成分分析では、PC1～PC4の累積寄与度が9割超となり、各資産クラスの価格変動の大半を説明するためにPC4まで考慮する必要があった。恐らく、この4つの内のどれかが暗号資産の価格変動に寄与するファクターで、それ以外が前述の経済成長・実質金利・インフレーションで解釈されるファクターとなることが想定される。

各主成分の時系列の傾向を確認して、主成分分析の結果を解釈してみた。

PC1（≒経済成長）

株式やREIT、ハイイールド社債といったリスク性資産が相対的に高いプラスで推移する一方で、国債がマイナスで推移していることから、経済成長を示すファクターであることが伺える

　

PC2（≒実質金利）

国債や物価連動国債が相対的に高いプラスで推移する一方で、その他資産がマイナスないしは低位なプラス水準で推移していることから、実質金利を示すファクターであることが伺える

　

PC3（≒???）

暗号資産（ビットコイン）が相対的に高い水準で推移する一方で、その他資産はゼロ近辺で推移していることから、暗号資産の価格変動に寄与する何かしらのファクター（というよりもビットコインそのもの?）であることが伺える

　

PC4（≒インフレーション?）

途中で正負の符号が逆転するものの、コモディティや物価連動国債が相対的に高いプラス（低いマイナス）で推移する一方で、その他資産がマイナス（プラス）ないしはゼロ近辺で推移していることから、恐らくインフレーションを示すファクターであることが伺える

####まとめ
暗号資産除きのデータに対する主成分分析を元にした解釈を前提にすると、上記のように解釈することができ、金融資産に影響を及ぼすPC1，PC2，PC4において、暗号資産（ビットコイン）はいずれもゼロ近辺で推移しており、ほぼ影響を及ぼしていないことが分かる。また、PC3が独立して暗号資産（ビットコイン）のみに影響を与えていることが読み取れた。

以上のように、金融資産の価格変動を説明する経済成長・実質金利・インフレとは異なる要因で、暗号資産の価格が変動していることが示された。本分析により、リスクファクターの観点から、暗号資産がポートフォリオ分散に寄与する可能性を検証し、成立することを確認することができた。

一方で、暗号資産に影響を及ぼすファクターが何であるかの解釈は課題として残った。将来的な分析課題として認識し、ビットコイン以外の暗号資産を加える分析等を実施して、より詳細な分析を行っていきたい~~…多分~~。

#最後に
今回、Pythonによる分析を実施してみて、EXCELやVBA以上に、WEB上にヒントが沢山ころがっているなと思いました。やりたいことが明確でデータさえ用意できるのであれば、~~極端な話、誰でも簡単に？~~Pythonで分析ができる環境が整っている素晴らしい時代になったな～と感じました。。
最後に一連の分析をまとめたコードを残し、本記事の執筆を終えたいと思います。

####一連の分析コードまとめ

分析コード(クリックで展開)

cryptopca.ipynb

pip install investpy # データ取得のための準備（最初に実施）

import pandas as pd #データ処理(データ形式の処理等)
import numpy as np #価格データ処理(数値計算)
import datetime #日付データ処理
import investpy #資産価格データ取得
import matplotlib.pyplot as plt# 図やグラフの図示
%matplotlib inline
import matplotlib.ticker as ticker

# 各資産クラスの価格取得（対USD）
# 日付設定
from_date = '01/12/2012' #取得開始日付
to_date = datetime.date.today().strftime('%d/%m/%Y')#取得終了日付(当該コード上では本日日付を利用)
price_df = pd.DataFrame(index=[], columns=[]) # 価格データ格納用の空のデータフレームを作成

# 暗号資産の価格を取得
#cc_lst = ['Bitcoin' ,'Ethereum' ,'XRP'] #①Bitcoin、②Ethereum、③XRP
#cc_lst = ['Bitcoin' ,'XRP'] #①Bitcoin、②XRP←より長期のデータがある2銘柄を選択
cc_lst = ['Bitcoin'] #①Bitcoinのみ選択
for crypto in cc_lst:
    cc_df = investpy.crypto.get_crypto_historical_data(crypto, from_date, to_date, as_json=False, order='ascending', interval='Monthly')
    cc_df.reset_index(drop=False, inplace=True)
    cc_df = cc_df.drop(['Open', 'High', 'Low', 'Volume', 'Currency'], axis=1) #終値だけ取得
    cc_df_new = cc_df.rename(columns={'Close': crypto})#名前変更
    
    if crypto == cc_lst[0]:
       price_df = price_df.append(cc_df_new, ignore_index=True)#データの追加(最初だけ日付データが入るように設定)
    
    else:
       price_df = pd.concat([price_df, cc_df_new[crypto]], axis=1)#データの追加(価格データのみ抽出)

# 指数データより取得できる資産クラスの価格を取得
i_lst = ['S&P 500', 'S&P GSCI Commodity TR', 'FTSE EPRA/NAREIT USA'] #④米国株式、⑤コモディティ、⑥米国REIT
#i_lst = ['S&P 500', 'FTSE EPRA/NAREIT USA'] #④米国株式、⑤米国REIT
#i_country = 'united states'
for index in i_lst:
    
    if index == i_lst[1]:
       i_country = 'world'
    else:
       i_country = 'united states'

    i_df = investpy.indices.get_index_historical_data(index, i_country, from_date, to_date, as_json=False, order='ascending', interval='Monthly')
    i_df.reset_index(drop=False, inplace=True)
    i_df = i_df.drop(['Open', 'High', 'Low', 'Volume', 'Currency'], axis=1) #終値だけ取得
    i_df_new = i_df.rename(columns={'Close': index})#名前変更
    
    price_df = pd.concat([price_df, i_df_new[index]], axis=1)##データの追加(価格データのみ抽出)

# ETFデータより取得できる資産クラスの価格を取得
etf_lst = ['iShares iBoxx $ High Yield Corporate Bond', 'iShares Core US Treasury Bond', 'iShares TIPS Bond'] #⑦ハイイールド社債、⑧米国債、⑨インフレ連動債
etf_country = 'united states'
for etf in etf_lst:
    etf_df = investpy.etfs.get_etf_historical_data(etf, etf_country, from_date, to_date, stock_exchange=None, as_json=False, order='ascending', interval='Monthly')
    etf_df.reset_index(drop=False, inplace=True)
    etf_df = etf_df.drop(['Open', 'High', 'Low', 'Volume', 'Currency', 'Exchange'], axis=1) #終値だけ取得
    etf_df_new = etf_df.rename(columns={'Close': etf})#名前変更
    price_df = pd.concat([price_df, etf_df_new[etf]], axis=1)#データの追加(価格データのみ抽出)

display(price_df.tail()) #データの確認

# 暗号資産除きのデータ
# 各資産クラスの価格取得（対USD）
# 日付設定
from_date = '01/12/2012' #取得開始日付
to_date = datetime.date.today().strftime('%d/%m/%Y')#取得終了日付(当該コード上では本日日付を利用)
price_finance_df = pd.DataFrame(index=[], columns=[]) # 価格データ格納用の空のデータフレームを作成

# 指数データより取得できる資産クラスの価格を取得
i_lst = ['S&P 500', 'S&P GSCI Commodity TR', 'FTSE EPRA/NAREIT USA'] #④米国株式、⑤コモディティ、⑥米国REIT
#i_lst = ['S&P 500', 'FTSE EPRA/NAREIT USA'] #④米国株式、⑤米国REIT
#i_country = 'united states'
for index in i_lst:
    
    if index == i_lst[1]:
       i_country = 'world'
    else:
       i_country = 'united states'

    i_df = investpy.indices.get_index_historical_data(index, i_country, from_date, to_date, as_json=False, order='ascending', interval='Monthly')
    i_df.reset_index(drop=False, inplace=True)
    i_df = i_df.drop(['Open', 'High', 'Low', 'Volume', 'Currency'], axis=1) #終値だけ取得
    i_df_new = i_df.rename(columns={'Close': index})#名前変更
    
    if index == i_lst[0]:
       price_finance_df = price_finance_df.append(i_df_new, ignore_index=True)#データの追加(最初だけ日付データが入るように設定)
    
    else:
       price_finance_df = pd.concat([price_finance_df, i_df_new[index]], axis=1)##データの追加(価格データのみ抽出)
    
# ETFデータより取得できる資産クラスの価格を取得
etf_lst = ['iShares iBoxx $ High Yield Corporate Bond', 'iShares Core US Treasury Bond', 'iShares TIPS Bond'] #⑦ハイイールド社債、⑧米国債、⑨インフレ連動債
etf_country = 'united states'
for etf in etf_lst:
    etf_df = investpy.etfs.get_etf_historical_data(etf, etf_country, from_date, to_date, stock_exchange=None, as_json=False, order='ascending', interval='Monthly')
    etf_df.reset_index(drop=False, inplace=True)
    etf_df = etf_df.drop(['Open', 'High', 'Low', 'Volume', 'Currency', 'Exchange'], axis=1) #終値だけ取得
    etf_df_new = etf_df.rename(columns={'Close': etf})#名前変更
    price_finance_df = pd.concat([price_finance_df, etf_df_new[etf]], axis=1)#データの追加(価格データのみ抽出)

display(price_finance_df.tail()) #データの確認

#価格⇒収益率データに変更
return_df = pd.DataFrame(index=[], columns=[])
for l in list(price_df.columns):
    if l == 'Date':
       return_df[l] = price_df[l]       
    else:
       return_df[l] = price_df[l].pct_change().dropna()

return_df = return_df.dropna(how='any') #欠損値を除外
r_df = return_df.set_index('Date') #日付をindex化したデータに変換(名前を変更)
#r_df[:-1].tail()#データの確認
print('データの確認')
display(r_df.tail())#データの確認
print()
print('相関係数の確認')
display(r_df.corr())#相関係数の確認
print()
print('標準偏差の確認')
display(r_df.std())#標準偏差の確認
#display(r_df.std()*(12**0.5))#標準偏差の確認

#暗号資産除きのデータ
#価格⇒収益率データに変更
return_finance_df = pd.DataFrame(index=[], columns=[])
for l in list(price_finance_df.columns):
    if l == 'Date':
       return_finance_df[l] = price_finance_df[l]       
    else:
       return_finance_df[l] = price_finance_df[l].pct_change().dropna()

return_finance_df = return_finance_df.dropna(how='any') #欠損値を除外
r_f_df = return_finance_df.set_index('Date') #日付をindex化したデータに変換(名前を変更)
#r_df[:-1].tail()#データの確認
print('データの確認')
display(r_f_df.tail())#データの確認
print()
print('相関係数の確認')
display(r_f_df.corr())#相関係数の確認
print()
print('標準偏差の確認')
display(r_f_df.std())#標準偏差の確認
#display(r_f_df.std()*(12**0.5))#標準偏差の確認

#データの前処理として標準化を実施
r_df = (r_df - r_df.mean(axis=0)) / r_df.std(axis=0)# 各資産の収益率データについて平均が0、分散・標準偏差が1となるように変換
print('データの確認')
display(r_df.tail())#データの確認
print()
print('相関係数の確認')
display(r_df.corr())#相関係数の確認
print()
print('標準偏差の確認')
display(r_df.std())#標準偏差の確認

#暗号資産除きのデータ
#データの前処理として標準化を実施
r_f_df = (r_f_df - r_f_df.mean(axis=0)) / r_f_df.std(axis=0)# 各資産の収益率データについて平均が0、分散・標準偏差が1となるように変換
print('データの確認')
display(r_f_df.tail())#データの確認
print()
print('相関係数の確認')
display(r_f_df.corr())#相関係数の確認
print()
print('標準偏差の確認')
display(r_f_df.std())#標準偏差の確認

from sklearn.decomposition import PCA #scikit-learnより主成分分析用のライブラリを取得

# 主成分分析のインスタンスを生成
pca = PCA()
c_rate = pd.DataFrame(index=[], columns=[])
eig_vec = pd.DataFrame(index=[], columns=[])

# データから変換モデルを学習し変換する。
t = 8 #分析期間(年間)
for i in range(len(r_df)-t*12):
    X=r_df[i:i+t * 12]
    X_pca = pca.fit_transform(X)
    c_rate = pd.concat([c_rate, pd.DataFrame(pca.explained_variance_ratio_, index=["PC{}".format(x + 1) for x in range(len(X.columns))])], axis=1)# 寄与率の計算
    eig_vec = pd.concat([eig_vec, pd.DataFrame(pca.components_.T, index = X.columns, columns = ["PC{}".format(x + 1) for x in range(len(X.columns))])], axis=1)# 固有ベクトル

# 主成分分析のインスタンスを生成
c_rate_f = pd.DataFrame(index=[], columns=[])
eig_vec_f = pd.DataFrame(index=[], columns=[])

# データから変換モデルを学習し変換する。
for i in range(len(r_f_df)-t*12):
    X=r_f_df[i:i+t * 12]
    X_pca = pca.fit_transform(X)
    c_rate_f = pd.concat([c_rate_f, pd.DataFrame(pca.explained_variance_ratio_, index=["PC{}_f".format(x + 1) for x in range(len(X.columns))])], axis=1)# 寄与率の計算
    eig_vec_f = pd.concat([eig_vec_f, pd.DataFrame(pca.components_.T, index = X.columns, columns = ["PC{}_f".format(x + 1) for x in range(len(X.columns))])], axis=1)# 固有ベクトル

c_rate = c_rate.T.set_index(price_df['Date'][(len(r_df)-t*12+1)*-1:-1]) # 寄与度データの転置とIndexの日付化
pc1 = eig_vec['PC1'].T.set_index(price_df['Date'][(len(r_df)-t*12+1)*-1:-1])*-1# PC1の転置とIndexの日付化、符号の反転
pc2 = eig_vec['PC2'].T.set_index(price_df['Date'][(len(r_df)-t*12+1)*-1:-1])# PC2の転置とIndexの日付化
pc3 = eig_vec['PC3'].T.set_index(price_df['Date'][(len(r_df)-t*12+1)*-1:-1])# PC3の転置とIndexの日付化
pc4 = eig_vec['PC4'].T.set_index(price_df['Date'][(len(r_df)-t*12+1)*-1:-1])# PC3の転置とIndexの日付化

c_rate_f = c_rate_f.T.set_index(price_df['Date'][(len(r_f_df)-t*12+1)*-1:-1]) # 寄与度データの転置とIndexの日付化
pc1_f = eig_vec_f['PC1_f'].T.set_index(price_df['Date'][(len(r_f_df)-t*12+1)*-1:-1])*-1# PC1の転置とIndexの日付化、符号の反転
pc2_f = eig_vec_f['PC2_f'].T.set_index(price_df['Date'][(len(r_f_df)-t*12+1)*-1:-1])# PC2の転置とIndexの日付化
pc3_f = eig_vec_f['PC3_f'].T.set_index(price_df['Date'][(len(r_f_df)-t*12+1)*-1:-1])# PC3の転置とIndexの日付化
pc4_f = eig_vec_f['PC4_f'].T.set_index(price_df['Date'][(len(r_f_df)-t*12+1)*-1:-1])# PC3の転置とIndexの日付化

# 寄与度の積み上げ棒グラフ表示
print('暗号資産除きのデータを対象とした主成分分析')
fig, ax = plt.subplots(2, 3, figsize=(15, 10), tight_layout=True, facecolor='whitesmoke')
c_rate_f.plot(kind='bar', stacked=True, ax=ax[0,0], title='Cumulative Contirbution')

# PCの棒グラフ表示
pc1_f.plot(kind='line', stacked=False, ax=ax[0,1], title='PC1')
pc2_f.plot(kind='line', stacked=False, ax=ax[0,2], title='PC2')
pc3_f.plot(kind='line', stacked=False, ax=ax[1,0], title='PC3')
pc4_f.plot(kind='line', stacked=False, ax=ax[1,1], title='PC4')
plt.show()

# 寄与度の積み上げ棒グラフ表示
print()
print('暗号資産含みのデータを対象とした主成分分析')
fig, ax = plt.subplots(2, 3, figsize=(15, 10), tight_layout=True, facecolor='whitesmoke')
c_rate.plot(kind='bar', stacked=True, ax=ax[0,0], title='Cumulative Contirbution')

# PCの棒グラフ表示
pc1.plot(kind='line', stacked=False, ax=ax[0,1], title='PC1')
pc2.plot(kind='line', stacked=False, ax=ax[0,2], title='PC2')
pc3.plot(kind='line', stacked=False, ax=ax[1,0], title='PC3')
pc4.plot(kind='line', stacked=False, ax=ax[1,1], title='PC4')
plt.show()

1. 参考文献：伊藤, 中川 [2018]『マクロファクターの定量化とリスク分析への活用』 ↩ ↩² ↩³