今週も特にありません

進捗どうですか?

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

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

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

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

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

データソース

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

www.quandl.com

日本では、日本円とビットコインの取引レート 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()

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を使うことでできます。

bitFlyer ビットコインを始めるなら安心・安全な取引所で

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

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

bitFlyer ビットコインを始めるなら安心・安全な取引所で

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の東芝ということがわかりました。

bitFlyer ビットコインを始めるなら安心・安全な取引所で

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 bitFlyer ビットコインを始めるなら安心・安全な取引所で

日本国債金利の主成分分析で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を使って簡単に始められますというものでした。

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

bitFlyer ビットコインを始めるなら安心・安全な取引所で

Rにおける最適化アルゴリズムをこれ一つで?

ほぼ使われているところを見ないoptimxパッケージを使ってみたメモ

データ分析を行っている方は、日々パラメータの推定などで最適化に励んでいると思います

Rで目的関数を自分で定義して最適化しようとする場合、optimやnlminb(一変数の場合はoptimize)等を用いて最適化を行います。 optimで最適化する場合、とりあえずデフォルトで実行されるNelder-Mead法を試して、この最適化が上手くいかなかった場合に、BFGS法などのその他の方法を試すという状況がよくあると思います

optimxパッケージは、Rで実行可能な最適化手法をまとめて一通り実行して、その結果を比較できるものです

尤度関数を定義して、そのパラメータを推定する例をやってみます

library(optimx)

loglik <- function(x, param) {
  -sum(dnbinom(x, size = param[1], mu = param[2], log = TRUE))
}

set.seed(71)
data <- rnbinom(500, size = 10, mu = 5)
initparam <- c(10, 5)
optimx(par = initparam, fn = loglik, x = data, control = list(all.methods = TRUE))

いくつかの警告文が出てきますが、最適化結果は以下のようになります

                   p1       p2         value fevals gevals niter convcode kkt1 kkt2 xtimes
BFGS         9.889487 5.038007  1.190725e+03     14      5    NA        0   NA   NA   0.01
CG           9.892441 5.037995  1.190725e+03    305    101    NA        1   NA   NA   0.27
Nelder-Mead  9.892127 5.038049  1.190725e+03     41     NA    NA        0   NA   NA   0.02
L-BFGS-B     9.889511 5.038000  1.190725e+03     13     13    NA        0   NA   NA   0.03
nlm          9.889480 5.037999  1.190725e+03     NA     NA    11        0   NA   NA   0.01
nlminb       9.889512 5.038000  1.190725e+03      8     18     5        0   NA   NA   0.00
spg          9.889768 5.038000  1.190725e+03     24     NA    18        0   NA   NA   0.13
ucminf       9.889507 5.037997  1.190725e+03      9      9    NA        0   NA   NA   0.02
Rcgmin             NA       NA 8.988466e+307     NA     NA    NA     9999   NA   NA   0.00
Rvmmin             NA       NA 8.988466e+307     NA     NA    NA     9999   NA   NA   0.00
newuoa       9.889511 5.038000  1.190725e+03     34     NA    NA        0   NA   NA   0.02
bobyqa       9.889509 5.038000  1.190725e+03     27     NA    NA        0   NA   NA   0.01
nmkb         9.890341 5.037991  1.190725e+03     69     NA    NA        0   NA   NA   0.03
hjkb        10.000000 5.000000  1.190775e+03      1     NA     0     9999   NA   NA   0.02

関数自体の使い方は普通のoptimとほとんど変わりません。 all.methods=TRUEを指定することですべての手法を実行します。 p1、p2が最適化したパラメータの値、valueがその時の目的関数の値で、 convcodeが0以外のものは最適化が正常に終了していないことを表しています

今回の場合では、結果が大きく変わることはないですが、 より複雑な場合だとそれぞれのアルゴリズムを一覧で比較できるところは使いどころかもしれません

より詳しい説明はJournal of statistical softwareにもあります。 CRAN Task View: Optimization and Mathematical Programmingを見ると、知らないパッケージがたくさんあるなー、まだまだRはいろいろとできるのだなーと感じます

bitFlyer ビットコインを始めるなら安心・安全な取引所で

一般化線形モデル周りのRのパッケージと関数まとめ

(GLM、GLMM周りを中心に)統計モデリングのためのパッケージと関数のメモ 使い方の詳しい説明、利用事例等を見つけたら、順次追加予定

lm {stats} Fitting Linear Models

線形モデル

nls {stats} Nonlinear Least Squares

非線形モデル

glm {stats} Fitting Generalized Linear Models

一般化線形モデル。線形回帰、ロジスティクス回帰、ポアソン回帰等の推定が可能。 cv.glm {boot}でクロスバリデーションを実行可能

glm.nb {MASS} Fit a Negative Binomial Generalized Linear Model

負の二項分布を仮定する場合の一般化線形モデル

gam {mgcv} Generalized additive models with integrated smoothness estimation

一般化加法モデル

gamm {mgcv} Generalized Additive Mixed Models

一般化加法混合モデル

plm {plm} Panel Data Estimators

パネルデータに対する線形モデル

polr {MASS} Ordered Logistic or Probit Regression

順序ロジット。プロビットモデルも可。polrは比例オッズモデルから来ている

mlogit {mlogit} Multinomial logit model

多項ロジットモデル。プロビットモデルも利用可能。(Kenneth Train's exercises using the mlogit package for R

mnlogit {mnlogit} Fast estimation of multinomial logit models

多項ロジットモデル。パッケージ作者による詳しい説明あり(Fast Estimation of Multinomial Logit Models: R Package mnlogit

npmlt {mixcat} Mixed effects cumulative link and logistic regression models

多項ロジットモデル

multinom {nnet} Fit Multinomial Log-linear Models

多項ロジットモデル

gmnl {gmnl} Estimate Multinomial Logit Models with Observed and Unobserved Individual Heterogeneity.

多項ロジットモデル。詳しい説明はこちら(http://msarrias.weebly.com/gmnl-package-in-r.html

lmer {lme4} Fit Linear Mixed-Effects Models

線形混合モデル。ランダム効果を複数指定可能。Journal of Statistical Softwareにパッケージの詳しい説明あり(Fitting Linear Mixed-Effects Models Using lme4 | Bates | Journal of Statistical Software

glmer {lme4} Fitting Generalized Linear Mixed-Effects Models

一般化線形混合モデル。ランダム効果を複数指定可

nlmer {lme4} Fitting Nonlinear Mixed-Effects Models

非線形混合モデル

glmer.nb {lme4} Fitting GLMM's for Negative Binomial

負の二項分布を仮定する場合の一般化線形混合モデル

glmmML {glmmML} Generalized Linear Models with random intercept

一般化線形混合モデル。ランダム効果は一つのみ指定可

glmnet {glmnet} fit a GLM with lasso or elasticnet regularization

lasso、ridge、elastic netによる正則化付き一般化線形モデル。線形回帰、ポアソン回帰、(多項)ロジスティクス回帰等の推定が可能。group lassoなども追加されてきている。cv.glmnet {glmnet}でクロスバリデーションを実行可能

glmreg {mpath} fit a GLM with lasso (or elastic net), snet or mnet regularization

正則化付き一般化線形モデル。lasso、scad、mcpなどを用いることができる
cv.glmreg {mpath}でクロスバリデーションを実行可能

zeroinfl {pscl} Zero-inflated Count Data Regression

Zero-inflatedモデル

vglm {VGAM} Fitting Vector Generalized Linear Models

ベクトル一般化線形モデル。Version 1.0-1の時点で「Currently only fixed-effects models are implemented, i.e., no random-effects models.」ということで、ランダム効果は考慮できない。Springerの本をもとにしている(Vector Generalized Linear and Additive Models - With an | Thomas W. Yee | Springer

vgam {VGAM} Fitting Vector Generalized Additive Models

ベクトル一般化加法モデル

bitFlyer ビットコインを始めるなら安心・安全な取引所で