TOPIXcore30銘柄のクラスタリング分析

# はじめに
TOPIXcore30銘柄の株価変化率のデータにクラスタリング分析をかけることにした。
データはYahoo Financeからライブラリ（yfinance）を利用して取得した。

# コード

## データをGoogle Driveに保存するまで
### Google Driveのマウント
```Python
from google.colab import drive
drive.mount('/content/drive')
```
### 必要なライブラリをインポート
```Python
# 必要なライブラリをインポート
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import yfinance as yf
from datetime import datetime, timedelta
```
### TOPIX Core30銘柄の株価を取得してGoogleドライブに保存
同じデータを使って分析を行うことになる場合も多いだろうから、取得したデータをGoogleドライブにcsvで保存しておく。しかし、 TOPIX Core30は結構銘柄入れ替えしてて、構成銘柄をキチンと拾うのが面倒。

```Python
# TOPIX Core30の構成銘柄
tickers = ['3382.T', '4063.T', '4502.T', '4503.T', '4568.T', '6098.T', '6273.T', '6367.T', 
           '6501.T', '6594.T', '6758.T', '6861.T', '6954.T', '6981.T', '7203.T', '7267.T', 
           '7741.T', '7974.T', '8001.T', '8031.T', '8035.T', '8058.T', '8306.T', '8316.T', 
           '8411.T', '8766.T', '9432.T', '9433.T', '9434.T', '9984.T']

# 銘柄と会社名の辞書
company_names = {
    '3382.T': 'Seven&I', '4063.T': 'ShinEtsu', '4502.T': 'Takeda', 
    '4503.T': 'Astellas', '4568.T': 'DaiichiSankyo', '6098.T': 'Recruit', '6273.T': 'SMC', '6367.T': 'Daikin', 
    '6501.T': 'Hitachi', '6594.T': 'Nidec', '6758.T': 'Sony', '6861.T': 'Keyence', '6954.T': 'Fanuc', 
    '6981.T': 'MuRata', '7203.T': 'Toyota', '7267.T': 'Honda', '7741.T': 'HOYA', '7974.T': 'Nintendo', 
    '8001.T': 'ITOCHU', '8031.T': 'Mitsubishi Corp', '8035.T': 'Tokyo Electron', '8058.T': 'Mitsubishi Trading', 
    '8306.T': 'MUFG', '8316.T': 'SMBC', '8411.T': 'Mizuho', '8766.T': 'Tokio Marine', 
    '9432.T': 'NTT', '9433.T': 'KDDI', '9434.T': 'SoftBank', '9984.T': 'SoftBankG'
}

# 株価データを取得
data = yf.download(tickers, period="max")['Adj Close']

# 現在の日付を取得
now = datetime.now()

# 日付を文字列に変換（形式：年-月-日）
date_str = now.strftime("%Y-%m-%d")

# Google DriveのパスにデータをCSVファイルとして保存
data.to_csv(f'/content/drive/MyDrive/stock_data_{date_str}.csv')
```
これで、株価データが今日の日付（2023/5/17）を入れたファイル名"stock_data_20230517.csv"として保存された。

# クラスタリング分析の実行
### 株価データの再読み込みと株価変化率の計算
保存したデータ"stock_data_20230517.csv"を呼び出して、日々の株価変化率を計算する。
```Python
# CSVファイルの再読込
data = pd.read_csv('/content/drive/MyDrive/stock_data_2023-05-17.csv', index_col=0)

### 日次変化率を計算
returns = data.pct_change()

### 欠損値を削除
returns = returns.dropna()
```


### 階層的クラスタリング分析
各種の階層的クラスタリング手法を試みることにした。

```Python
# クラスタリング手法のリスト
methods = ['ward', 'single', 'complete', 'average']

for method in methods:
    g = sns.clustermap(returns.corr(), cmap='coolwarm', linewidths=0.1, 
                       xticklabels=[company_names[ticker] for ticker in returns.columns], 
                       yticklabels=[company_names[ticker] for ticker in returns.columns], 
                       method=method)
    
    plt.setp(g.ax_heatmap.get_xticklabels(), rotation=90)  # x軸のラベルを90度回転
    g.fig.suptitle(f'Clustering using {method} method', fontsize=16, y=1.01)  # タイトルを設定し、少し上にずらす
    plt.show()
```

#### 階層的クラスタリング手法による出力結果（ヒートマップ）
それぞれの手法ごとに描画したところ、以下の図が描けた。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3305624/594dbcec-bbd0-a5e2-a887-61a50e9b9dce.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3305624/72b7f498-f9e5-f2ae-b33b-3731faa8b1dd.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3305624/717785ed-58dd-34d7-b23f-b757785e26a6.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3305624/ec84a0f0-533c-296b-2044-765ad1c80aab.png)


###  非階層的クラスタリング
以下では、非階層的クラスタリング手法のうちk-means法によってみる。クラスター数はシルエット法によって決めることにした。

```Python
# 必要なライブラリの追加インポート
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score, silhouette_samples
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

# クラスタ数の候補
n_clusters_candidates = list(range(2, 11))

# 各クラスタ数でのシルエットスコアを保存するリスト
silhouette_scores = []

# クラスタ数の候補ごとにk-means法とシルエット法を適用
best_n_clusters = None
best_silhouette_score = -1
best_labels = None
best_silhouette_vals = None
for n_clusters in n_clusters_candidates:
    # クラスタリングモデルを作成
    kmeans = KMeans(n_clusters=n_clusters, n_init=10)
    
    # 変動の相関行列に対してクラスタリングを実行
    kmeans.fit(returns.corr())
    
    # クラスタリング結果を取得
    labels = kmeans.labels_
    
    # シルエットスコアを計算
    score = silhouette_score(returns.corr(), labels)
    
    # 各サンプルのシルエットスコアを計算
    silhouette_vals = silhouette_samples(returns.corr(), labels)
    
    # シルエットスコアをリストに追加
    silhouette_scores.append(score)
    
    # シルエットスコアが最高の場合、そのクラスタ数とシルエットスコア、ラベル、シルエット値を保存
    if score > best_silhouette_score:
        best_n_clusters = n_clusters
        best_silhouette_score = score
        best_labels = labels
        best_silhouette_vals = silhouette_vals

# 最適なクラスタ数とその時のシルエットスコアを表示
print(f"The best number of clusters is {best_n_clusters} with a silhouette score of {best_silhouette_score}")

# シルエットスコアのプロット
plt.plot(n_clusters_candidates, silhouette_scores, 'o-')
plt.xlabel('Number of clusters')
plt.ylabel('Silhouette Score')
plt.grid(True)
plt.title('Silhouette scores for different numbers of clusters')
plt.show()

# 最適なクラスタ数でのクラスタリング結果を表示
for i in range(best_n_clusters):
    print(f"Cluster {i+1}:")  # クラスタ番号を1から始める
    print([company_names[ticker] for ticker, label in zip(returns.columns, best_labels) if label == i])

# 銘柄とクラスタラベルのデータフレームを作成
clusters_df = pd.DataFrame({'Ticker': returns.columns, 'Cluster': best_labels + 1}) 

# 銘柄名を追加
clusters_df['Company'] = clusters_df['Ticker'].apply(lambda x: company_names[x])

print(clusters_df)
```

#### 出力結果（シルエット・スコア）
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3305624/c4c3e013-9d12-494e-4032-0e9528e09163.png)

Cluster 1:
['Hitachi', 'Toyota', 'Honda', 'ITOCHU', 'Mitsubishi Corp', 'Mitsubishi Trading', 'MUFG', 'SMBC', 'Mizuho', 'Tokio Marine']
Cluster 2:
['ShinEtsu', 'Fanuc', 'Tokyo Electron']
Cluster 3:
['Recruit', 'SMC', 'Daikin', 'Nidec', 'Sony', 'Keyence', 'MuRata', 'HOYA', 'SoftBankG']
Cluster 4:
['Seven&I', 'Takeda', 'Astellas', 'DaiichiSankyo', 'Nintendo', 'NTT', 'KDDI', 'SoftBank']
    Ticker  Cluster             Company
0   3382.T        4             Seven&I
1   4063.T        2            ShinEtsu
2   4502.T        4              Takeda
3   4503.T        4            Astellas
4   4568.T        4       DaiichiSankyo
5   6098.T        3             Recruit
6   6273.T        3                 SMC
7   6367.T        3              Daikin
8   6501.T        1             Hitachi
9   6594.T        3               Nidec
10  6758.T        3                Sony
11  6861.T        3             Keyence
12  6954.T        2               Fanuc
13  6981.T        3              MuRata
14  7203.T        1              Toyota
15  7267.T        1               Honda
16  7741.T        3                HOYA
17  7974.T        4            Nintendo
18  8001.T        1              ITOCHU
19  8031.T        1     Mitsubishi Corp
20  8035.T        2      Tokyo Electron
21  8058.T        1  Mitsubishi Trading
22  8306.T        1                MUFG
23  8316.T        1                SMBC
24  8411.T        1              Mizuho
25  8766.T        1        Tokio Marine
26  9432.T        4                 NTT
27  9433.T        4                KDDI
28  9434.T        4            SoftBank
29  9984.T        3           SoftBankG

まあ、概ね順当な結果かなとは思うが、Cluster4がセブン＆アイ、製薬会社（武田薬品工業、アステラス製薬、第一三共）、任天堂、通信事業者（NTT,KDDI、ソフトバンク）で構成されているのはやや意外。

もう少し詳しく見るために、シルエット図を描いてみることにした。
```Python
# シルエット図のプロット
fig, ax = plt.subplots(1, 1, figsize=(10, 8))
y_lower, y_upper = 0, 0
yticks = []
for i, cluster in enumerate(np.unique(best_labels)):
    cluster_silhouette_vals = silhouette_vals[best_labels == cluster]
    cluster_silhouette_vals.sort()
    y_upper += len(cluster_silhouette_vals)
    ax.barh(range(y_lower, y_upper), cluster_silhouette_vals, edgecolor='none', height=1)
    ax.text(-0.03, (y_lower + y_upper) / 2, str(i + 1))
    y_lower += len(cluster_silhouette_vals)
    
    # 銘柄名の追加
    for j, ticker in enumerate(returns.columns[best_labels == cluster]):
        company_name = company_names[ticker]  # Tickerを会社名に変換
        ax.text(0.95, y_lower - len(cluster_silhouette_vals) + j, company_name, ha='left', va='center')  # テキストの位置を左に移動

# 平均シルエットスコアを縦線でプロット
avg_score = np.mean(best_silhouette_vals)
ax.axvline(avg_score, linestyle='--', linewidth=2, color='black')
ax.set_yticks([])
ax.set_xlim([-0.1, 1.1])  # x軸の範囲を少し広げて銘柄名を表示できるようにする
ax.set_xlabel('Silhouette coefficient values')
ax.set_ylabel('Cluster labels')
ax.set_title('Silhouette plot for the various clusters', y=1.02);

plt.tight_layout()
plt.show()
```
#### 出力結果（シルエット図）
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3305624/693c2ce4-1c6d-cc1a-024e-b8f753d8c672.png)


#### 主成分分析
まだよくわからないので、主成分分析を使って次元削減して、2次元プロットしてみる。
図に描いたときに銘柄名が重なるのを避けるため、パッケージ（adjustText）を利用したいが、GoogleColabではプリインストールされていないので、ここでインストール。
  
```Python
pip install adjustText
```
そこで、主成分分析。
```Python
# パッケージの追加インポート
from adjustText import adjust_text
from sklearn.decomposition import PCA

# PCAのインスタンスを作成
pca = PCA(n_components=2)

# PCAを適用
pca_results = pca.fit_transform(returns.corr())

# プロット
plt.figure(figsize=(10,8))
texts = []  # テキストオブジェクトを保存するリスト
for i, cluster in enumerate(np.unique(best_labels)):
    plt.scatter(pca_results[best_labels == cluster, 0], pca_results[best_labels == cluster, 1], label=f'Cluster {i+1}')  # クラスタ番号を1から始める
    for j, ticker in enumerate(returns.columns[best_labels == cluster]):
        texts.append(plt.text(pca_results[best_labels == cluster][j, 0],  # テキストオブジェクトを追加
                              pca_results[best_labels == cluster][j, 1], 
                              company_names[ticker]))  # 会社名をプロット

# テキストラベルの位置を自動調整
adjust_text(texts)

plt.xlabel('First Principal Component')
plt.ylabel('Second Principal Component')
plt.title('PCA Scatter Plot for Clusters')
plt.legend()
plt.show()
```

### 出力結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3305624/64af6592-5a77-757e-7a2d-12bf44b95f63.png)

わかったようでわからない図だが、まあこれは主成分分析の常なんで仕方ないか。

# おわりに
何をやってるか結局よくわからなくなったが、Yahoo Financeのデータも利用できたし、クラスタリング分析、さらに主成分分析までできたので、今回はこれで終了。