今週も特にありません

進捗どうですか?

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の値を変更すれば良い。

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。

すべて違いを確認したものの、細かい装飾の類が多く、好みの問題でどれを設定しておいた方が良さそうというものはなかった...

時系列 Tidydata の予測 fasster を試してみる

時系列データもTidydataとして扱う流れが加速しているようです。 xtsなどで頑張っていた時代からアップデートしていきましょう。

時系列データをTidy化するtsibbleパッケージに関連して、モデリングに対応するfassterパッケージも絶賛開発されています。 今回はfassterを使って時系列データのモデリングから予測までを行い、その使い心地を試してみます。

バージョン確認

gitのDevelopment cycleにも書かれている通り、早いペースで開発が進められており、今後も大幅な仕様変更も想定されます。 ここではバージョン0.1.0.9000を試しており、ほぼ出たばかりの開発途上のテストバージョンとなっています。 数ヶ月後にはここに書かれていることそのままでは動かない可能性もありますのでご注意ください。

github.com

パッケージ読み込み

パッケージは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

f:id:masaqol:20180716122118j:plain

予測

上記のモデルをもとに、5年先まで予測します。

fit %>% 
  forecast(h = 20) %>%
  autoplot

f:id:masaqol:20180716122600j:plain

まとめ

fassterでTidy化された時系列データに対してモデリングから予測まで行えました。

実際には、dlmModPolydlmFilterといったdlmパッケージのお馴染みの関数が中で動いています。 dlm自体は使い慣れないと書きにくい部分もあるので、Tidy化した時系列データに対しては圧倒的に使い心地は良いことがわかりました。

dlmを直接いじることが過去になる日も近そうです。

tidyquantでビットコインの取引データを取得する

2017年は何かと仮想通貨の価格が高騰(or 暴落)が騒がれた一年でした。

仮想通貨が法定通貨と同様に決済手段金融商品の一つとして定着するようになった際には、今以上に様々な分析が行われることでしょう。

その前に、tidyverseな世界観でのデータ取得から可視化まではどうやるの?は、世の中的に少し需要があるのではないかと思ったため調べてみました。

今回はデータ取得部分を中心にメモとして書いておきます。

データソース

まず、tidyquant::tq_getで指定できるビットコインに関するコードを探したところ、Quandl にいくつかそれに該当するものが見つかりました。

www.quandl.com

日本では、日本円とビットコインの取引レート BTC/JPY が興味の対象があると思われるため、大手の仮想通貨取引所の価格と取引高のデータが分析対象としては最適そうです。

データ取得

データを取得できる期間が不明なため、bitFlyerZaifcoincheck の取引所のデータを取得します。

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()

f:id:masaqol:20171222181725p:plain

どうやら 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を使った分析メモ(3)

tidyquantを使った分析メモ

前回までに個別銘柄のリターンについて計算したので、ポートフォリオのリターンについて出してみます。 PerformanceAnalytics::Returns.portfolioのラッパー関数であるtq_portfolioを使えば、tidyなデータを入力として簡単にポートフォリオのリターンを計算できます。

github.com

masaqol.hatenablog.com

masaqol.hatenablog.com

データ

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")

f:id:masaqol:20171001162821p:plain

それぞれの銘柄の期待リターンとボラティリティを計算します。

> 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)

f:id:masaqol:20171001162832p:plain

ということで、あとはポートフォリオのウェイトを最適化して、 ポートフォリオのリターンやボラティリティを詳細に分析していけば良いことがわかります。

tidyquantを使った分析メモ(2)

tidyquantを使った分析メモ

今回は対数リターンの計算とその可視化までをやってみます。

github.com

masaqol.hatenablog.com

対数リターンの計算

TOPIX100構成銘柄の株価データを取得したい場合、銘柄コードを用意して、tidyquanttq_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" 
...

対数リターンの計算は、quantmodperiodReturnを利用します。

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")

f:id:masaqol:20170722111530p:plain

対数リターンの分布を確認します。(銘柄コードが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")

f:id:masaqol:20170722111649p:plain

期待リターンとボラティリティ

最後に各銘柄の期待リターンとボラティリティを計算し、高い銘柄を確認します。

> 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

期待リターントップは、6273のSMC。ボラティリティトップは6502の東芝ということがわかりました。

tidyquantを使った分析メモ

R-bloggersでも度々登場して気になっていたtidyquantを使った分析メモ

github.com

基本的な使い方は、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"のようにして問い合わせれば、データを取得できます。 (取得に時間がかかるために、3ヶ月分としています。)

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")

f:id:masaqol:20170521123132p:plain

日本国債金利の主成分分析で2016年振り返り

マイナス金利や長短金利操作付き量的・質的金融緩和導入など、今年も何かと金融政策の変更とそれに伴う金利の変動がニュースを賑わせた一年でした。 久々に何か分析ネタを…ということで、日本国債金利を主成分分析することで、金利の構造の主なファクターが2016年はどのように推移したのかを調べたいと思います。

データはQuandlから持ってきます。 Quandlパッケージを使えば、data.framextsなどの形でデータを使えるようになり、手軽に分析を始めることができます。

まずは、必要なライブラリを読み込みます。

library(Quandl)
library(dplyr)
library(tidyr)
library(ggplot2)
library(scales)

Quandlからデータを取ってくるには、データを表すコードと取得開始日を最低限入れます。 ここでは、日本国債金利データをそのままを使うことにします。 金利の変動について分析を行いたい場合は、Quandl関数のtransformdiffを指定することで取ってくることができます。

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))

f:id:masaqol:20161231123705p:plain

1月のマイナス金利導入以降、イールドカーブ全体的に低下傾向となり、2月終盤からは10年満期以下の国債金利が全てマイナス圏で推移するようになります。 7月には15年満期の国債金利もマイナスという状態も見られるようになり、年間を通して一番金利が低下した時期であったことがわかります。 それ以降は全体的に上昇し、アメリカの大統領選挙後には、10年満期の国債金利もマイナス圏から脱し、金利の上昇が加速しているように見えます。

それでは、このデータに対して主成分分析を行います。 これらの分析については、様々なところで行われているため詳細な説明は省きますが、金利は水準(Level)、傾き(Slope)、曲率(Curveture)の3つのファクターでほぼその構造が説明できると言われています。主成分分析によって、これらのファクターが導出されるのかを確認し、その時系列推移を可視化してみます。

特にひねりを加えるわけでもなく、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

第3主成分までで累積寄与率が99%となり、ほぼこの3つでイールドカーブの構造を説明できていることがわかります。 主成分得点は以下のようになりました。

> 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) 

f:id:masaqol:20161231124810p:plain

第1から第3主成分までが金利の水準、曲率、長短期の傾きをそれぞれ表していることがわかります。 (横軸の間隔が満期に対応して一定ではないので、形がちょっと変な感じですが…)

さらに、これらのファクターの時系列推移も可視化して、一年を通してどのように変動したのかを見てみます。

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")

f:id:masaqol:20161231125156p:plain

1月のマイナス金利導入以降、金利の水準が大きく下落し、傾きや曲率も変動したことがわかります。 そして、年後半には再び水準が上昇したことがわかります。 年間を通して、金利が下落した後に上昇して戻ってきたと言っても、 水準と傾き、曲率の3つのファクターを見ると、年初と年末ではイールドカーブの構造が異なっていることがわかります。

ただ単純に主成分分析しただけのものなので、ガチなフラットナー/スティープナーといった金利関連の取引に役立つわけではないと思いますが、データの取得から可視化、分析手法適用までRを使って簡単に始められますというものでした。

来年も世界各国の政治イベントで大きな動きがありそうということで、勝手にデータ分析の材料を提供してくれると思うと楽しめそうです。