tibbleのオプション設定
表示数の変更
tibble（tbl_df）型のデータを表示すると、デフォルトでは先頭の10行が表示されます。 加えて、最大表示数は20行に設定されています。
library(tidyverse) class(diamonds) #> [1] "tbl_df" "tbl" "data.frame" diamonds #> # A tibble: 53,940 x 10 #> carat cut color clarity depth table price x y #> <dbl> <ord> <ord> <ord> <dbl> <dbl> <int> <dbl> <dbl> #> 1 0.23 Ideal E SI2 61.5 55 326 3.95 3.98 #> 2 0.21 Premium E SI1 59.8 61 326 3.89 3.84 #> 3 0.23 Good E VS1 56.9 65 327 4.05 4.07 #> 4 0.290 Premium I VS2 62.4 58 334 4.2 4.23 #> 5 0.31 Good J SI2 63.3 58 335 4.34 4.35 #> 6 0.24 Very G… J VVS2 62.8 57 336 3.94 3.96 #> 7 0.24 Very G… I VVS1 62.3 57 336 3.95 3.98 #> 8 0.26 Very G… H SI1 61.9 55 337 4.07 4.11 #> 9 0.22 Fair E VS2 65.1 61 337 3.87 3.78 #> 10 0.23 Very G… H VS1 59.4 61 338 4 4.05 #> # ... with 5.393e+04 more rows, and 1 more variable: z <dbl>
これらを変更するにはtibble.print_max, tibble.print_minの値を変更すればOKです。
tibble_opt <- list( "tibble.print_max" = 100, "tibble.print_min" = 20 ) options(tibble_opt) getOption("tibble.print_max") #> [1] 100 getOption("tibble.print_min") #> [1] 20 diamonds #> # A tibble: 53,940 x 10 #> carat cut color clarity depth table price x y #> <dbl> <ord> <ord> <ord> <dbl> <dbl> <int> <dbl> <dbl> #> 1 0.23 Ideal E SI2 61.5 55 326 3.95 3.98 #> 2 0.21 Premium E SI1 59.8 61 326 3.89 3.84 #> 3 0.23 Good E VS1 56.9 65 327 4.05 4.07 #> 4 0.290 Premium I VS2 62.4 58 334 4.2 4.23 #> 5 0.31 Good J SI2 63.3 58 335 4.34 4.35 #> 6 0.24 Very G… J VVS2 62.8 57 336 3.94 3.96 #> 7 0.24 Very G… I VVS1 62.3 57 336 3.95 3.98 #> 8 0.26 Very G… H SI1 61.9 55 337 4.07 4.11 #> 9 0.22 Fair E VS2 65.1 61 337 3.87 3.78 #> 10 0.23 Very G… H VS1 59.4 61 338 4 4.05 #> 11 0.3 Good J SI1 64 55 339 4.25 4.28 #> 12 0.23 Ideal J VS1 62.8 56 340 3.93 3.9 #> 13 0.22 Premium F SI1 60.4 61 342 3.88 3.84 #> 14 0.31 Ideal J SI2 62.2 54 344 4.35 4.37 #> 15 0.2 Premium E SI2 60.2 62 345 3.79 3.75 #> 16 0.32 Premium E I1 60.9 58 345 4.38 4.42 #> 17 0.3 Ideal I SI2 62 54 348 4.31 4.34 #> 18 0.3 Good J SI1 63.4 54 351 4.23 4.29 #> 19 0.3 Good J SI1 63.8 56 351 4.23 4.26 #> 20 0.3 Very G… J SI1 62.7 59 351 4.21 4.27 #> # ... with 5.392e+04 more rows, and 1 more variable: z <dbl>
その他のオプション
上記は.Rprofileに仕込んでいる場合も多いのではないかと思い、別に何か有用なオプションがあるかどうか調べてみました。
tibble.width：出力幅をいくつにするか？デフォルトはNULL。tibble.max_extra_cols：表示できないカラム名の最大表示数をいくつにするか？デフォルトは100。pillar.bold：カラム名の表示をボールドフォントにするか？デフォルトはFALSE。pillar.subtle：データの次元やカラムのデータ型の表示に細かい濃淡を付けるか？デフォルトはTRUE。pillar.neg：負の値をハイライトするか？デフォルトはTRUE。pillar.sigfig：表示する有効数字の数をいくつにするか？デフォルトは3。pillar.min_title_chars：カラム名の最低表示文字数をいくつにするか？デフォルトは15。
すべて違いを確認したものの、細かい装飾の類が多く、好みの問題で、どれを設定しておいた方が良さそうというものはないように思いました...
tidy時系列データにおける予測 fasster
時系列データもTidydataとして扱う流れが加速しているようです。
xtsなどで頑張っていた時代からアップデートしていきます。
時系列データをTidy化するtsibbleパッケージに関連して、モデリングに対応するfassterパッケージも絶賛開発されています。
今回はfassterを使って時系列データのモデリングから予測までを行い、その使い心地を試してみます。
バージョン確認
gitのDevelopment cycleにも書かれている通り、早いペースで開発が進められており、今後も大幅な仕様変更も想定されます。 ここではバージョン0.1.0.9000を試しており、ほぼ出たばかりの開発途上のテストバージョンとなっています。 数ヶ月後にはここに書かれていることそのままでは動かない可能性もあり、注意が必要です。
パッケージ読み込み
パッケージはgitからインストールします。
# install.packages("devtools") devtools::install_github("tidyverts/fasster")
gitのREADMEも暫定的な書き方のようなので、READMEをなぞる場合には、tidyverseと合わせて、tsibble関連のパッケージをまとめて読み込んでおいた方が全てスムーズに実行できます。
library(tidyverse) library(tsibble) library(fable) library(fasster) library(tsibbledata)
Tidyデータへの変換
時系列データとしてよく用いられるUKgasをもとに見ていきます。
tsクラスをtbl_tsクラスとして扱うにはas_tsibbleを使うだけです。
UKgas %>% as_tsibble #> # A tsibble: 108 x 2 [1QUARTER] #> index value #> <qtr> <dbl> #> 1 1960 Q1 160. #> 2 1960 Q2 130. #> 3 1960 Q3 84.8 #> 4 1960 Q4 120. #> 5 1961 Q1 160. #> 6 1961 Q2 125. #> 7 1961 Q3 84.8 #> 8 1961 Q4 117. #> 9 1962 Q1 170. #> 10 1962 Q2 141. #> # ... with 98 more rows
1QUARTERとなっているように、4半期ごとに観測されたデータになっています。
モデリング
Local Level + 季節性の単純なものでモデルを当てはめてみます。
fit <- UKgas %>% as_tsibble %>% FASSTER(value ~ poly(1) + seas(4)) fit %>% summary #> FASSTER Model: #> value ~ poly(1) + seas(4) #> #> Estimated variances: #> State noise variances (W): #> poly(1) #> 3.0239e+01 #> seas(4) #> 1.2268e+02 1.2366e-17 2.9014e-19 #> #> Observation noise variance (V): #> 1.0449e+02
formulaに状態空間モデルとして取り入れる項を書くだけで、簡単に観測ノイズと状態ノイズの分散を推定してくれます。
さらに、各要素の値も確認しておきます。
fit %>% components #> # A tsibble: 108 x 3 [1QUARTER] #> index `poly(1)` `seas(4)` #> <qtr> <dbl> <dbl> #> 1 1960 Q1 123. 38.6 #> 2 1960 Q2 123. 6.84 #> 3 1960 Q3 123. -40.0 #> 4 1960 Q4 123. -5.61 #> 5 1961 Q1 124. 37.1 #> 6 1961 Q2 124. 6.65 #> 7 1961 Q3 122. -37.7 #> 8 1961 Q4 122. -3.71 #> 9 1962 Q1 122. 38.2 #> 10 1962 Q2 124. 1.04 #> # ... with 98 more rows
これらの各要素の時系列推移の結果もすぐに可視化できます。
fit %>% autoplot
予測
上記のモデルをもとに、5年先まで予測します。
fit %>% forecast(h = 20) %>% autoplot
まとめ
fassterでTidy化された時系列データに対してモデリングから予測まで行えました。
実際には、dlmModPolyやdlmFilterといったdlmパッケージのお馴染みの関数が中で動いています。
dlm自体は使い慣れないと書きにくい部分もあるので、Tidy化した時系列データに対しては圧倒的に使い心地は良いことがわかりました。
dlmを直接いじることが過去になる日も近そうです。
tidyquantでビットコインの取引データを取得する
2017年は何かと仮想通貨の価格が高騰（or 暴落）が騒がれた一年でした。
仮想通貨が法定通貨と同様に決済手段や金融商品の一つとして定着するようになった際には、今以上に様々な分析が行われると思います。
その前に、tidyverseな世界観でのデータ取得から可視化まではどうやるの？は、世の中的に少し需要があるのではないかと思ったため調べてみました。
今回はデータ取得部分を中心にメモとして書いておきます。
データソース
まず、tidyquant::tq_getで指定できるビットコインに関するコードを探したところ、Quandl にいくつかそれに該当するものが見つかりました。
日本では、日本円とビットコインの取引レート BTC/JPY が興味の対象があると思われるため、大手の仮想通貨取引所の価格と取引高のデータが分析対象としては最適化と思われます。
データ取得
データを取得できる期間が不明なため、bitFlyer、Zaif、coincheck の取引所のデータを取得します。
library(tidyverse) library(tidyquant) codes <- c("BCHARTS/BITFLYERJPY", "BCHARTS/ZAIFJPY", "BCHARTS/COINCHECKJPY") btc_data <- tq_get(x = codes, get = "quandl", from = "2017-01-01", to = "2017-12-21")
取得したデータを確認します。
> btc_data # A tibble: 684 x 9 symbol date open high low close volume.btc. volume.currency. weighted.price <chr> <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> 1 BCHARTS/BITFLYERJPY 2017-07-04 296127 297040 287072 292501 6401.075 1871425686 292361.2 2 BCHARTS/BITFLYERJPY 2017-07-05 292501 295800 286220 292127 9259.952 2689073855 290398.3 3 BCHARTS/BITFLYERJPY 2017-07-06 292252 294710 290650 293850 7480.301 2193005651 293170.8 4 BCHARTS/BITFLYERJPY 2017-07-07 293900 294350 285432 288084 8908.557 2580602677 289676.9 5 BCHARTS/BITFLYERJPY 2017-07-08 288121 291988 285693 291821 5627.901 1625558830 288839.3 6 BCHARTS/BITFLYERJPY 2017-07-09 291943 293333 288500 288855 4854.138 1411796907 290844.0 7 BCHARTS/BITFLYERJPY 2017-07-10 288855 289750 267370 273753 15486.759 4352953461 281075.8 8 BCHARTS/BITFLYERJPY 2017-07-11 273751 277930 265850 266592 14649.179 3972753107 271192.9 9 BCHARTS/BITFLYERJPY 2017-07-12 266770 274839 258888 272350 13183.987 3508777924 266139.4 10 BCHARTS/BITFLYERJPY 2017-07-13 272360 276799 264489 267990 9524.230 2578378059 270717.8 # ... with 674 more rows
データ自体は取得できていそうですが、bitFlyerのデータが7月の初めくらいからしか取得できていなさそうです。
データ可視化
全体を確認するために close の値を可視化します。
btc_data %>% select(symbol, date, close) %>% ggplot(aes(x = date, y = close, colour = symbol)) + geom_line(size = 1) + theme_tq() + scale_color_tq()
どうやら coincheck 以外は、7月以降のデータしか取得できないようです。 加えて、データが欠損（値が0になっている）ところがあるようです。
長期間のデータを使いたい場合は coincheck を選択すれば良さそうです。
欠損値への対応
bitFlyer のデータを使いたい場合は欠損値（この場合は0）を処理する必要があります。 どう欠損値を処理すれば良いかはその時の状況にもよりますが、とりあえず、欠損している部分を削除してしまう場合は、以下のようにすれば良さそうです。
dplyr::if_elseを使う場合は、double型に対応させて0をNA_real_に置き換えてna.omitでその行のデータを取り除きます。
bitflyer_btc_data <- btc_data %>% filter(symbol == "BCHARTS/BITFLYERJPY") %>% mutate_if(is.numeric, funs(if_else(. == 0, NA_real_, .))) %>% drop_na
NA値を削除せずに直前のデータで埋める場合には、tidyr::fillを使うことでできます。
tidyquantを使った分析メモ ポートフォリオのリターン
tidyquantを使った分析メモの続き。前回までに個別銘柄のリターンについて計算したので、ポートフォリオのリターンについて出してみます。
PerformanceAnalytics::Returns.portfolioのラッパー関数であるtq_portfolioを使えば、tidyなデータを入力として簡単にポートフォリオのリターンを計算できます。
データ
TOPIX100 の銘柄の株価データを取得して、対数リターンを計算します。 ここでは、期間を2017年7月から9月末までを指定しています。
> log_returns # A tibble: 6,138 x 4 # Groups: symbol [99] symbol date close log_return <chr> <date> <dbl> <dbl> 1 YJ1605.T 2017-07-03 1087.0 0.000000000 2 YJ1605.T 2017-07-04 1099.0 0.010979067 3 YJ1605.T 2017-07-05 1086.0 -0.011899454 4 YJ1605.T 2017-07-06 1073.0 -0.012042758 5 YJ1605.T 2017-07-07 1069.0 -0.003734832 6 YJ1605.T 2017-07-10 1067.5 -0.001404166 7 YJ1605.T 2017-07-11 1072.0 0.004206597 8 YJ1605.T 2017-07-12 1073.5 0.001398276 9 YJ1605.T 2017-07-13 1079.5 0.005573633 10 YJ1605.T 2017-07-14 1083.0 0.003236997 # ... with 6,128 more rows
ポートフォリオの銘柄
TOPIX100 に含まれる特定の銘柄を使ってポートフォリオのリターンを計算します。 ここでは、「セブン&アイ・ホールディングス」「楽天」「ファナック」「任天堂」「ファーストリテイリング」「ソフトバンク」の6銘柄を使うことにします。
stocks <- c("YJ3382.T", "YJ4755.T", "YJ6954.T", "YJ7974.T", "YJ9983.T", "YJ9984.T")
それぞれの銘柄の対数リターンの時系列推移を可視化します。
log_returns %>% filter(symbol %in% stocks) %>% ggplot(aes(x = date, y = log_return, colour = symbol)) + geom_line() + theme_tq() + scale_color_tq() + facet_wrap(~ symbol, ncol = 3, scales = "free_y")
それぞれの銘柄の期待リターンとボラティリティを計算します。
> log_returns %>% + filter(symbol %in% targets) %>% + group_by(symbol) %>% + summarise(return = mean(log_return), vol = sd(log_return)) # A tibble: 6 x 3 symbol return vol <chr> <dbl> <dbl> 1 YJ4063.T -1.277554e-04 0.01137096 2 YJ4755.T -1.080331e-03 0.01011779 3 YJ6954.T 8.090697e-04 0.01075738 4 YJ7974.T 1.909432e-03 0.02136985 5 YJ9983.T -1.943404e-03 0.01333176 6 YJ9984.T -4.786856e-05 0.01177579
ファナックと任天堂の期待リターンがプラスであることがわかります。
ポートフォリオのリターン
各銘柄のウェイトを決めた上で、tq_portfolioでポートフォリオのリターンを計算します。
> weights <- c(rep(0.15, 4), rep(0.2, 2)) > portfolio_returns <- log_returns %>% + filter(symbol %in% stocks) %>% + tq_portfolio(assets_col = symbol, + returns_col = log_return, + weights = weights) %T>% + print # A tibble: 62 x 2 date portfolio.returns <date> <dbl> 1 2017-07-03 0.0000000000 2 2017-07-04 -0.0076961084 3 2017-07-05 -0.0001438849 4 2017-07-06 -0.0048844521 5 2017-07-07 -0.0031134117 6 2017-07-10 0.0087959661 7 2017-07-11 0.0033212960 8 2017-07-12 -0.0052806985 9 2017-07-13 -0.0014577677 10 2017-07-14 -0.0065329481 # ... with 52 more rows
このポートフォリオの期待リターンは以下のようになります。
> portfolio_returns %>% + summarise(return = mean(portfolio.returns), vol = sd(portfolio.returns)) # A tibble: 1 x 2 return vol <dbl> <dbl> 1 -0.0003031948 0.008354069
各銘柄のボラティリティに比べて、ポートフォリオのボラティリティが小さくなっていることがわかります。
ポートフォリオのリターンの推移を可視化します。
portfolio_returns %>% ggplot(aes(x = date, y = portfolio.returns)) + geom_bar(stat = "identity") + theme_tq() + scale_color_tq() + scale_y_continuous(labels = scales::percent)
ということで、あとはポートフォリオのウェイトを最適化して、 ポートフォリオのリターンやボラティリティを詳細に分析していけば良いことがわかります。
tidyquantを使った分析メモ 期待リターンとボラティリティ
tidyquantを使った分析メモの続き。今回は対数リターンの計算とその可視化までをやってみます。
対数リターンの計算
TOPIX100構成銘柄の株価データを取得したい場合、銘柄コードを用意して、tidyquantのtq_getを使うことで簡単に取得することができるところまで前回確認しました。
> TOPIX100 # A tibble: 6,100 × 8 symbol date open high low close volume <chr> <date> <dbl> <dbl> <dbl> <dbl> <dbl> 1 YJ1605.T 2017-01-04 1182.5 1201.0 1180.0 1190.0 5927100 2 YJ1605.T 2017-01-05 1186.0 1187.0 1148.5 1168.0 8302500 3 YJ1605.T 2017-01-06 1139.0 1158.5 1131.0 1153.5 7737100 4 YJ1605.T 2017-01-10 1141.5 1153.0 1136.0 1136.0 4807700 5 YJ1605.T 2017-01-11 1140.5 1156.5 1133.5 1151.5 5587400 6 YJ1605.T 2017-01-12 1143.0 1156.5 1134.5 1145.5 4831300 7 YJ1605.T 2017-01-13 1130.5 1152.5 1130.5 1150.5 5257000 8 YJ1605.T 2017-01-16 1144.5 1150.0 1127.0 1131.5 3200200 9 YJ1605.T 2017-01-17 1136.0 1138.0 1117.0 1122.0 3625500 10 YJ1605.T 2017-01-18 1116.0 1146.5 1112.0 1143.0 5458900 # ... with 6,090 more rows, and 1 more variables: adjusted <dbl>
株価そのままのデータでは、銘柄ごとの比較などが行いにくいために、対数リターンを計算します。
tidyquantは複数のファイナンス関係のパッケージの便利なラッパーとして機能します。
対数リターンなどの基本的な計算は、元のパッケージの便利な関数を使うことでtidyなデータに対して処理を行えます。
tq_mutateによって、新しい変数として対数リターンを追加します。
適用できる関数の一覧はtq_mutate_fun_optionsによって確認できます。
> tq_mutate_fun_options() $zoo [1] "rollapply" "rollapplyr" [3] "rollmax" "rollmax.default" [5] "rollmaxr" "rollmean" [7] "rollmean.default" "rollmeanr" [9] "rollmedian" "rollmedian.default" [11] "rollmedianr" "rollsum" [13] "rollsum.default" "rollsumr" $xts [1] "apply.daily" "apply.monthly" "apply.quarterly" [4] "apply.weekly" "apply.yearly" "diff.xts" [7] "lag.xts" "period.apply" "period.max" [10] "period.min" "period.prod" "period.sum" [13] "periodicity" "to_period" "to.daily" [16] "to.hourly" "to.minutes" "to.minutes10" ...
対数リターンの計算は、quantmodのperiodReturnを利用します。
log_returns <- TOPIX100 %>% group_by(symbol) %>% tq_mutate( select = close, mutate_fun = periodReturn, period = "daily", type = "log", col_rename = "log_return" ) %>% select(symbol, date, close, log_return)
計算したものを確認します。
> log_returns Source: local data frame [6,100 x 4] Groups: symbol [100] symbol date close log_return <chr> <date> <dbl> <dbl> 1 YJ1605.T 2017-01-04 1190.0 0.000000000 2 YJ1605.T 2017-01-05 1168.0 -0.018660423 3 YJ1605.T 2017-01-06 1153.5 -0.012492086 4 YJ1605.T 2017-01-10 1136.0 -0.015287478 5 YJ1605.T 2017-01-11 1151.5 0.013552120 6 YJ1605.T 2017-01-12 1145.5 -0.005224217 7 YJ1605.T 2017-01-13 1150.5 0.004355408 8 YJ1605.T 2017-01-16 1131.5 -0.016652444 9 YJ1605.T 2017-01-17 1122.0 -0.008431379 10 YJ1605.T 2017-01-18 1143.0 0.018543578 # ... with 6,090 more rows
データの可視化
対数リターンを時系列で可視化します。（銘柄コードが7000番台の銘柄のみ）
log_returns %>% filter(grepl("^YJ7", symbol)) %>% ggplot(aes(x = date, y = log_return, colour = symbol)) + geom_line() + theme_tq() + scale_color_tq() + facet_wrap(~ symbol, ncol = 3, scales = "free_y")
対数リターンの分布を確認します。（銘柄コードが7000番台の銘柄のみ）
log_returns %>% filter(grepl("^YJ7", symbol)) %>% ggplot(aes(x = log_return, fill = symbol)) + geom_histogram(binwidth = 0.005) + theme_tq() + scale_color_tq() + facet_wrap(~ symbol, ncol = 3, scales = "free_y")
期待リターンとボラティリティ
最後に各銘柄の期待リターンとボラティリティを計算し、高い銘柄を確認します。
> return_vol <- log_returns %>% + group_by(symbol) %>% + summarise(return = mean(log_return), vol = sd(log_return)) > > return_vol %>% + arrange(desc(return)) %>% + head(10) # A tibble: 10 × 3 symbol return vol <chr> <dbl> <dbl> 1 YJ6273.T 0.002274630 0.016555097 2 YJ2502.T 0.002011907 0.012576164 3 YJ6758.T 0.002002305 0.013206705 4 YJ6954.T 0.001991144 0.013172032 5 YJ4452.T 0.001708134 0.010844302 6 YJ6861.T 0.001573498 0.009784067 7 YJ4188.T 0.001551808 0.016332976 8 YJ2503.T 0.001353853 0.011863891 9 YJ7741.T 0.001352014 0.011299796 10 YJ8035.T 0.001253836 0.014759141 > > return_vol %>% + arrange(desc(vol)) %>% + head(10) # A tibble: 10 × 3 symbol return vol <chr> <dbl> <dbl> 1 YJ6502.T -2.278774e-03 0.05451714 2 YJ8795.T -8.095544e-05 0.02178426 3 YJ7186.T -2.028825e-03 0.02028579 4 YJ8750.T -1.430680e-04 0.01997049 5 YJ4755.T -7.893335e-04 0.01927645 6 YJ9064.T -7.722044e-04 0.01899148 7 YJ7974.T 8.831755e-04 0.01850704 8 YJ4578.T -1.365315e-04 0.01812934 9 YJ8630.T 8.039943e-06 0.01787958 10 YJ9983.T -3.522327e-03 0.01776665
tidyquantを使った分析メモ データの取得と可視化
R-bloggersでも度々登場して気になっていたtidyquantを使った分析メモになります。
基本的な使い方は、Vignettesが充実しています。 また、日本語での資料としてはtidyquantの使い方があります。 また、Exploratoryと組み合わせて分析した記事もあります。 qiita.com
ここでは、日本の株式銘柄をそれなりに大量取得して、データを可視化するところまでやってみます。
銘柄コード一覧を作成
まず、株価データを取得する際に必要な銘柄コード一覧を作成します。 ここでは、TOPIX100構成銘柄をまとめて取得することを考えます。
銘柄の入れ替え等があった場合に適宜更新を行ってくれると思われる
松井証券さんのページを信用できる情報として、
rvestでスクレイプして作成します。
library(tidyverse) library(tidyquant) library(rvest) url <- "http://www.matsui.co.jp/service/stock/trade/rule_topix100/" html <- read_html(url, encoding = "UTF-8") stock_name_df <- html %>% html_node("table.m-tbl") %>% html_table
htmlの中からm-tblクラスのtable部分をdata.frameとして取得しています。
> stock_name_df 銘柄コード 銘柄名 1 1605 国際石油開発帝石 2 1878 大東建託 3 1925 大和ハウス工業 4 1928 積水ハウス 5 2502 アサヒグループホールディングス 6 2503 キリンホールディングス 7 2802 味の素 8 2914 日本たばこ産業 9 3382 セブン&アイ・ホールディングス 10 3402 東レ 11 3407 旭化成 ...
データ取得で必要なのは銘柄コードのみのため、これだけを切り出せば完了です。
stock_code <- stock_name_df %>% select(`銘柄コード`)
データを確認します。
> stock_code 銘柄コード 1 1605 2 1878 3 1925 4 1928 5 2502 6 2503 7 2802 8 2914 9 3382 10 3402 11 3407 ...
株価データの取得
銘柄コード一覧が作成できたら、それをもとに株価データを取得します。
ここでは、tq_get()のget = "stock.prices.japan"を指定して、ヤフーファイナンスからデータを取得することにします。
この際に銘柄コードは、"YJ9984.T"のようにして問い合わせれば、データを取得できます。 （取得に時間がかかるために、３ヶ月分としています。）
TOPIX100 <- stock_code %>% unlist %>% as.character %>% paste0("YJ", ., ".T") %>% tq_get(get = "stock.prices.japan", from = "2017-01-01", to = "2017-03-31")
データを確認します。
> TOPIX100 # A tibble: 6,100 × 8 symbol date open high low close volume <chr> <date> <dbl> <dbl> <dbl> <dbl> <dbl> 1 YJ1605.T 2017-01-04 1182.5 1201.0 1180.0 1190.0 5927100 2 YJ1605.T 2017-01-05 1186.0 1187.0 1148.5 1168.0 8302500 3 YJ1605.T 2017-01-06 1139.0 1158.5 1131.0 1153.5 7737100 4 YJ1605.T 2017-01-10 1141.5 1153.0 1136.0 1136.0 4807700 5 YJ1605.T 2017-01-11 1140.5 1156.5 1133.5 1151.5 5587400 6 YJ1605.T 2017-01-12 1143.0 1156.5 1134.5 1145.5 4831300 7 YJ1605.T 2017-01-13 1130.5 1152.5 1130.5 1150.5 5257000 8 YJ1605.T 2017-01-16 1144.5 1150.0 1127.0 1131.5 3200200 9 YJ1605.T 2017-01-17 1136.0 1138.0 1117.0 1122.0 3625500 10 YJ1605.T 2017-01-18 1116.0 1146.5 1112.0 1143.0 5458900 # ... with 6,090 more rows, and 1 more variables: adjusted <dbl>
データの可視化
取得したデータを可視化します。全ての銘柄をプロットすると見にくいため、 銘柄コードが7000番台の銘柄のみの株価推移を可視化します。
TOPIX100 %>% filter(grepl("^YJ7", symbol)) %>% ggplot(aes(x = date, y = close, colour = symbol)) + geom_line() + theme_tq() + scale_color_tq() + facet_wrap(~ symbol, ncol = 3, scales = "free_y")
日本国債金利の主成分分析
マイナス金利や長短金利操作付き量的・質的金融緩和導入など、今年も何かと金融政策の変更とそれに伴う金利の変動がニュースを賑わせた一年でした。 久々に何か分析ネタを…ということで、日本国債金利を主成分分析することで、金利の構造の主なファクターが2016年はどのように推移したのかを調べたいと思います。
データはQuandlから持ってきます。
Quandlパッケージを使えば、data.frameやxtsなどの形でデータを使えるようになり、手軽に分析を始めることができます。
まずは、必要なライブラリを読み込みます。
library(Quandl) library(dplyr) library(tidyr) library(ggplot2) library(scales)
Quandlからデータを取ってくるには、データを表すコードと取得開始日を最低限入れます。
ここでは、日本国債金利データをそのままを使うことにします。
金利の変動について分析を行いたい場合は、Quandl関数のtransformでdiffを指定することで取ってくることができます。
irate <- Quandl("MOFJ/INTEREST_RATE_JAPAN", order = "asc", start_date = "2016-01-01")
取ってきたデータを確認します。
> head(irate) Date 1Y 2Y 3Y 4Y 5Y 6Y 1 2016-01-04 -0.042 -0.028 -0.014 0.003 0.026 0.036 2 2016-01-05 -0.043 -0.030 -0.014 0.003 0.026 0.036 3 2016-01-06 -0.044 -0.029 -0.014 0.002 0.022 0.032 4 2016-01-07 -0.045 -0.028 -0.014 -0.003 0.016 0.027 5 2016-01-08 -0.045 -0.031 -0.018 -0.003 0.016 0.022 6 2016-01-12 -0.048 -0.028 -0.018 -0.003 0.017 0.016 7Y 8Y 9Y 10Y 15Y 20Y 25Y 30Y 40Y 1 0.059 0.120 0.186 0.264 0.599 0.986 1.150 1.282 1.402 2 0.059 0.120 0.181 0.259 0.588 0.971 1.135 1.264 1.386 3 0.055 0.115 0.176 0.253 0.583 0.967 1.131 1.260 1.382 4 0.049 0.110 0.171 0.245 0.573 0.958 1.116 1.238 1.363 5 0.039 0.101 0.161 0.230 0.557 0.939 1.096 1.220 1.346 6 0.034 0.097 0.157 0.221 0.548 0.930 1.086 1.207 1.338
1年から40年満期までの金利データが取れました。 早速、データを可視化します。
irate %>% gather(Maturity, Rate, -Date) %>% mutate(Maturity = factor(.$Maturity, levels = paste0(c(1:10, 15, 20, 25, 30, 40), "Y"))) %>% ggplot(aes(x = Date, y = Rate)) + geom_line(aes(colour = Maturity), size = 1) + scale_x_date(labels = date_format("%y/%m"), date_breaks = "3 month") + guides(colour = guide_legend(reverse = TRUE))
１月のマイナス金利導入以降、イールドカーブ全体的に低下傾向となり、２月終盤からは10年満期以下の国債の金利が全てマイナス圏で推移するようになります。 ７月には15年満期の国債の金利もマイナスという状態も見られるようになり、年間を通して一番金利が低下した時期であったことがわかります。 それ以降は全体的に上昇し、アメリカの大統領選挙後には、10年満期の国債の金利もマイナス圏から脱し、金利の上昇が加速しているように見えます。
それでは、このデータに対して主成分分析を行います。 これらの分析については、様々なところで行われているため詳細な説明は省きますが、金利は水準(Level)、傾き(Slope)、曲率(Curveture)の３つのファクターでほぼその構造が説明できると言われています。主成分分析によって、これらのファクターが導出されるのかを確認し、その時系列推移を可視化してみます。
特にひねりを加えるわけでもなく、Rのprcompを金利データに適用します。
irate_pc <- irate %>% select(-Date) %>% prcomp(scale. = TRUE)
主成分分析を適用した結果を見てみます。
> summary(irate_pc) Importance of components: PC1 PC2 PC3 PC4 PC5 Standard deviation 3.654 0.99449 0.75082 0.24082 0.13649 Proportion of Variance 0.890 0.06593 0.03758 0.00387 0.00124 Cumulative Proportion 0.890 0.95589 0.99347 0.99734 0.99858 PC6 PC7 PC8 PC9 PC10 Standard deviation 0.09270 0.07582 0.05226 0.03636 0.03167 Proportion of Variance 0.00057 0.00038 0.00018 0.00009 0.00007 Cumulative Proportion 0.99915 0.99953 0.99972 0.99980 0.99987 PC11 PC12 PC13 PC14 PC15 Standard deviation 0.02533 0.02314 0.01826 0.01654 0.01175 Proportion of Variance 0.00004 0.00004 0.00002 0.00002 0.00001 Cumulative Proportion 0.99991 0.99995 0.99997 0.99999 1.00000
第３主成分までで累積寄与率が99％となり、ほぼこの３つでイールドカーブの構造を説明できていることがわかります。 主成分得点は以下のようになりました。
> head(irate_pc$x) PC1 PC2 PC3 PC4 PC5 [1,] -9.534255 -0.3139355 0.4521807 -0.16670843 -0.2093294 [2,] -9.425388 -0.2323307 0.4851745 -0.16565373 -0.1661761 [3,] -9.329534 -0.2554974 0.5167301 -0.11327598 -0.1747176 [4,] -9.152626 -0.2405698 0.5779134 -0.08359126 -0.2015301 [5,] -8.915227 -0.2231325 0.6567429 -0.03116720 -0.1449082 [6,] -8.801091 -0.1916643 0.7010752 0.05046791 -0.1497915 PC6 PC7 PC8 PC9 PC10 [1,] -0.0451359411 0.057984040 -0.006210947 0.03531776 0.006980638 [2,] -0.0351578833 0.054768411 0.009829469 0.03076594 0.008573947 [3,] -0.0347377299 0.033557155 0.012547679 0.03020654 0.010608943 [4,] -0.0531471021 -0.004868424 0.018739106 0.05401370 -0.006719837 [5,] -0.0005787395 0.023142080 0.005714631 0.02306513 -0.002223937 [6,] -0.0028323655 0.009862543 -0.014780631 -0.01505679 -0.007297795 PC11 PC12 PC13 PC14 PC15 [1,] -0.02586290 -0.03433084 -0.0009471546 -0.016075984 0.02848080 [2,] -0.03508559 -0.03784112 0.0147360705 -0.012732611 0.02126039 [3,] -0.02178643 -0.04119171 0.0143781817 -0.010653086 0.02193993 [4,] -0.02087492 -0.02625881 0.0161335068 -0.015374304 0.01088626 [5,] -0.01278054 0.01291703 0.0168199370 -0.008567272 0.02197490 [6,] -0.02024915 0.05001019 0.0143638963 -0.011216030 0.01494438
そして、各満期に対するファクターローディングスを見てみます。
data.frame(irate_pc$rotation[, 1:3]) %>% mutate(Maturity = factor(rownames(.), levels = rownames(.))) %>% gather(PC, FactorLodings, -Maturity) %>% ggplot(aes(x = Maturity, y = FactorLodings)) + geom_line(aes(colour = PC, group = PC), size = 1)
第１から第３主成分までが金利の水準、曲率、長短期の傾きをそれぞれ表していることがわかります。 （横軸の間隔が満期に対応して一定ではないので、形がちょっと変な感じですが…）
さらに、これらのファクターの時系列推移も可視化して、一年を通してどのように変動したのかを見てみます。
data.frame(Date = irate$Date, irate_pc$x[, 1:3]) %>% gather(PC, Factor, -Date) %>% ggplot(aes(x = Date, y = Factor)) + geom_line(aes(colour = PC), size = 1) + scale_x_date(labels = date_format("%y/%m"), date_breaks = "3 month")
１月のマイナス金利導入以降、金利の水準が大きく下落し、傾きや曲率も変動したことがわかります。 そして、年後半には再び水準が上昇したことがわかります。 年間を通して、金利が下落した後に上昇して戻ってきたと言っても、 水準と傾き、曲率の３つのファクターを見ると、年初と年末ではイールドカーブの構造が異なっていることがわかります。
ただ単純に主成分分析しただけのものなので、ガチなフラットナー/スティープナーといった金利関連の取引に役立つわけではないと思いますが、データの取得から可視化、分析手法適用までRを使って簡単に始められますというものでした。
来年も世界各国の政治イベントで大きな動きがありそうということで、勝手にデータ分析の材料を提供してくれると思うと楽しめそうです。
