Machine Learning – ページ 9

今回は, 取引関係のある企業の株価変化率の時系列を k-Shape という時系列クラスタリング手法で調べてみます。

企業間取引と株価
k-Shape
R で k-Shape

実行環境は以下です。

R: 3.5.1
dtwclust: 5.5.1
dplyr: 0.7.8

企業間取引と株価

総売上に対して10%以上を占める取引先がある場合, 有価証券報告書の「生産、受注及び販売の状況」に開示する義務がある。企業間取引 (BtoB) と株価の関係は, 企業間の取引関係をノードを企業, エッジを取引関係としてネットワーク構造で表すことができる。
今回は, ネットワーク構造の分析ではなく取引関係のある企業の株価変化率 (Price Rate Of Change) の時系列の類似性を調べてみる。
サプライヤー企業の売上高に占める特定顧客に対する売上高の割合が大きいほどカスタマー企業の株価が変化したとき, そのサプライヤー企業の株価への反応も大きいと考えられる。参考書籍では例として, カスタマー企業である任天堂に対してサプライヤー企業であるコネクタやスイッチを製造するホシデン [6804], 電子部品を製造するミツミ電機 (2017年にミネベアと経営統合しミネベアミツミ [6479]) の株価の関係が紹介されている。

k-Shape

k-Shape [Paparrizos and Gravano, 2015] は時系列の形状 (shape) に着目した時系列クラスタリング手法。k-Shape の説明に入る前に時系列の不変性と, よく使われる時系列間の距離について触れる。

Time-Series Invariances

系列間の比較を行うとき多くのタスクで以下の不変性 (Invariances) のいくつかまたは全てを満たしたい。

1. Scaling and translation invariances

系列 x に対して振幅定数 a とオフセット定数 b が異なる系列は似ていると判定してほしい。

$\begin{eqnarray*} \vec{x} = a\vec{x} + b \end{eqnarray*}$

2. Shift invariance

2つの系列が位相のみ異なる場合にも似ていると判定してほしい。

3. Uniform scaling invariance

2つの系列の長さが異なる場合, 短い系列を伸ばす (stretching) または長い系列を縮める (shrinking) ことで比較したい。

4. Occlusion invariance

部分系列が欠測しているときに, その部分は無視して系列を比較したい。

5. Complexity invariance

系列の形状は似ているが複雑さが異なる場合は, アプリケーションの種類によっては系列の類似性を低くまたは高くしたい (e.g. 屋外と屋内で録音した音響信号はノイズに差があるが似ていると判定してほしい)

Time-Series Distance Measures

時系列間の比較の距離尺度として ED (Euclidean Distance) と DTW (Dynamic Time Warping) がよく使われる。

m を系列の長さとして2つの時系列 x = (x1, …, xm), y = (y1, …, ym) の ED は以下となる。

$\begin{eqnarray*} ED(\vec{x}, \vec{y}) = \sqrt{ \sum_{i=1}^m (x_i - y_i)^2} \end{eqnarray*}$

DTW は ED の拡張で, 局所的で非線形なアライメントが可能。

$\begin{eqnarray*} DTW(\vec{x}, \vec{y}) = min \sqrt{ \sum_{i=1}^k w_i \end{eqnarray*}$

k-Shape では Shape-based distance (SBD) という距離を提案している。

k-Shape algorithm

k-Shape クラスタリングは Scaling と Shifting に対する不変性に焦点を当てている。k-Shape の主な特徴として Shape-based distance (SBD) と, 時系列の形状抽出 (Time-Series Shape Extraction) の2つがある。

SBD

相互相関 (Cross-correlation) は信号処理の分野でよく使われる尺度で以下で表せる。計算効率化のために DFT でなく FFT (+α) を用いる。

$\begin{eqnarray*} CC(\vec{x}, \vec{y}) = F^{-1}{F(\vec{x}) ∗ F(\vec{y})} \end{eqnarray*}$

正規化された相互相関 (係数正規化) NCCc は相互相関系列を個々の系列の自己相関の幾何平均で割った値。 NCCc が最大となる位置 ω を検出する。

$\begin{eqnarray*} SBD(\vec{x}, \vec{y}) = 1 - max_w \left(\frac{CC_w(\vec{x}, \vec{y})}{\sqrt{R_0(\vec{x}, \vec{x}) \cdot R_0(\vec{y}, \vec{y})}} \right) \end{eqnarray*}$

SBD は 0 から 2の間の値を取り, 0 に近いほど2つの時系列は類似している。

Shape Extraction

SBD による時系列クラスタリングのための重心 (centroid) ベクトルを以下で求める。詳しい表記は論文を参照。

$\begin{eqnarray*} \vec{\mu}_{k}^{*} = argmax_{\vec{\mu}_{k}} \sum_{\vec{x}_i \in P_k} NCC_c(\vec{x_i}, \vec{\mu}_{k})^2 \\ = argmax_{\vec{\mu}_{k}} \sum_{\vec{x}_i \in P_k} \left(\frac{max_w CC_w(\vec{x_i}, \vec{\mu}_{k})}{\sqrt{R_0(\vec{x_i}, \vec{x_i}) \cdot R_0(\vec{\mu}_{k}}, \vec{\mu}_{k}})}} \right) \\ = argmax_{\vec{\mu}_{k}} \sum_{\vec{x}_i \in P_k} (\vec{x_i} \cdot \vec{\mu}_{k})^2 \\ = argmax_{\vec{\mu}_{k}} \vec{\mu}_{k}^{\mathrm{T}} \cdot \sum_{\vec{x}_i \in P_k} (\vec{x_i} \cdot \vec{x_i}^{\mathrm{T}})^2 \cdot \vec{\mu}_{k} \\ = argmax_{\vec{\mu}_{k}} \vec{\mu}_{k}^{\mathrm{T}} \cdot S \cdot \vec{\mu}_{k} \\ = argmax_{\vec{\mu}_{k}} \frac{\vec{\mu}_{k}^{\mathrm{T}} \cdot Q^{\mathrm{T}} \cdot S \cdot Q \cdot \vec{\mu}_{k}}{\vec{\mu}_{k}^{\mathrm{T}} \cdot \vec{\mu}_{k}} \\ = argmax_{\vec{\mu}_{k}} \frac{\vec{\mu}_{k}^{\mathrm{T}} \cdot M \cdot \vec{\mu}_{k}}{\vec{\mu}_{k}^{\mathrm{T}} \cdot \vec{\mu}_{k}} \end{eqnarray*}$

k-Shape の全体のアルゴリズムは以下となる。

k-Shape は k-means のように反復的な手続きにより各時系列にクラスタを割り当てる。

各時系列を各クラスタの重心ベクトルと比較し, 最も近い重心ベクトルのクラスタに割り当てる
クラスタの重心ベクトルを更新する

クラスタのメンバに変更が発生しないか, 反復回数が最大値に達するまで上記の手順 1, 2 を繰り返す。

論文の実験では, 実行時間の比較やクラスタリングの性能指標である Rand Index での比較, また統計的な解析として複数のデータセットに渡る複数のアルゴリズムの比較のために Friedman test と事後検定 (post-hoc test) のひとつである Nemenyi test で検証を行なっている。

R で k-Shape

R で k-Shape を実行するために dtwclust パッケージを用いる。

今回は, ホシデン, ミネベアミツミに加えてプリント配線板を製造するサプライヤー企業のシライ電子工業 [6658], また任天堂との企業間取引は少ないであろう積水ハウス [1928] とみずほフィナンシャルグループ [8411] の計6銘柄に対して株価変化率の時系列クラスタリングを試してみた。

各銘柄は以下のカラムを持つ。期間は 2014-01-01 ~ 2018-12-31 までの5年間とし, 株価変化率 (Price Rate Of Change) は (当日終値 – n日前終値) / n日前終値の n = 1 で計算した。

> start <- "2014-01-01"
> df_7974 %>%
+     filter(date > as.Date(start))
# A tibble: 1,222 x 10
   date        open  high   low close   volume close_adj change rate_of_change  code
                                 
 1 2014-01-06 14000 14330 13920 14320  1013000     14320    310       0.0221    7974
 2 2014-01-07 14200 14380 14060 14310   887900     14310    -10      -0.000698  7974
 3 2014-01-08 14380 16050 14380 15850  3030500     15850   1540       0.108     7974
 4 2014-01-09 15520 15530 15140 15420  1817400     15420   -430      -0.0271    7974
 5 2014-01-10 15310 16150 15230 16080  2124100     16080    660       0.0428    7974
 6 2014-01-14 15410 15755 15370 15500  1462200     15500   -580      -0.0361    7974
 7 2014-01-15 15750 15880 15265 15360  1186800     15360   -140      -0.00903   7974
 8 2014-01-16 15165 15410 14940 15060  1606600     15060   -300      -0.0195    7974
 9 2014-01-17 15100 15270 14575 14645  1612600     14645   -415      -0.0276    7974
10 2014-01-20 11945 13800 11935 13745 10731500     13745   -900      -0.0615    7974
# … with 1,212 more rows

銘柄ごとに期間内の行数が異なるためもっとも長い銘柄を基準として前営業日の値で欠測補完した。(k-Shape では多少のずれは許容されるが念のため)

align_length <- function(df, criterion) {
  if (nrow(df) == nrow(criterion)) return(df)

  diff <- as.Date(setdiff(criterion$date, df$date), origin = "1970-01-01")

  for (i in diff) {
   position <- df %>%
     filter(date < i) %>%
     nrow()
   row <- df[position,]
   row$date <- as.Date(i, origin = "1970-01-01")
   df <- df %>%
     bind_rows(row) %>%
     arrange(date)
  }

  return(df)
}

stocks <- list(df_1928, df_6479, df_6658, df_6804, df_7974, df_8411)

argmax <- stocks %>%
  map(nrow) %>%
  unlist() %>%
  which.max()

aligned_stocks <- stocks %>%
  map(function(x) {
    x %>%
      align_length(criterion = stocks[argmax][[1]])
  }) %>%
  map(function(x) {
    x %>%
      filter(date > as.Date(start)) %>%
      pull(rate_of_change)
  })

各銘柄の株価と株価変化率。

tsclust() の series に数値ベクトルの list, type に “partitional”, preproc に zscore, distance に “sbd”, centroid = “shape” を指定すると k-Shape となる。

res <- tsclust(series = aligned_stocks,
               type = "partitional",
               preproc = zscore,
               distance = "sbd",
               centroid = "shape")

df_res <- data.frame(cluster = res@cluster,
                     centroid_dist = res@cldist,
                     code = c("1928", "6479", "6658", "6804", "7974", "8411"),
                     name = c("積水ハウス", "ミネベアミツミ", "シライ電子工業", "ホシデン", "任天堂", "みずほ"))

クラスタリングの結果は以下となった。

> df_res %>%
+     arrange(cluster)
  cluster centroid_dist code           name
1       1     0.1897561 1928     積水ハウス
2       1     0.2196533 6479 ミネベアミツミ
3       1     0.1481051 8411         みずほ
4       2     0.3468301 6658 シライ電子工業
5       2     0.2158674 6804       ホシデン
6       2     0.2372485 7974         任天堂

任天堂, ホシデン, シライ電子工業が同じ cluster に割り当てられた。2016年のホシデンの任天堂に対する販売比率は 50.5% であり, 企業間の取引関係が株価の連動にも影響していることが窺える。
一方, ミネベアミツミは別の cluster となったが, ミツミがミネベアと経営統合したのは2017年であり Pockemon Go がリリースされた2016年7月頃の株価の急騰には反応しなかったことも要因のひとつかもしれない。

おわりに

Python だと tslearn に k-Shape が実装されています。
参考書籍は『実践金融データサイエンス』です。

[1] Time Series Classification

タグ: Machine Learning

【R】k-Shape による株価の時系列クラスタリング