Vasicekモデルの債券オプション
RでVasicekモデルの債券価格とイールドで、金利モデル関連のRパッケージがない…と書いたが、SMFI5というパッケージで一応、VasicekモデルやCIRモデルの債券価格を求めたりできるようだ
「Statistical Methods for Financial Engineering」の5章に関してまとめパッケージとなっている。本を読んでいないので、なぜ5章だけなのかは分からない
このパッケージでも債券オプションについては入っていないので、前回に続き、Vasicekモデルにおける債券オプション価格について
# zero coupon bond option price under Vasicek model # # Args: # kappa : speed of reversion # mu : long term mean level # sigma : instaneous volatility # r0 : current short rate # t : current time # T : bond maturity # S : option maturity # K : strike price # L : face value # # Return: # data.frame of moneyness, call and put option prices # VasicekZCBOption <- function(kappa, mu, sigma, r0, t, T, S, K, L) { sp <- sigma / kappa * (1 - exp(-kappa * (S - T))) * sqrt((1 - exp(-2 * kappa * (T - t))) / (2 * kappa)) PT <- VasicekZCBond(kappa, mu, sigma, r0, t, T)$price PS <- VasicekZCBond(kappa, mu, sigma, r0, t, S)$price d1 <- log(L * PS / (K * PT)) / sp + sp / 2 d2 <- d1 - sp M <- K * PT / (L * PS) CP <- L * PS * pnorm(d1) - K * PT * pnorm(d2) PP <- K * PT * pnorm(-d2) - L * PS * pnorm(-d1) data.frame(moneyness = M, callprice = CP, putprice = PP) }
関数の中でVasicekモデルの債券価格関数を使っているので、読み込んでおく。前回と同様のモデルのパラメータをセット
kappa <- 0.2 mu <- 0.015 sigma <- 0.003 r0 <- -0.001 t <- 0 T <- 1/12 S <- 3 K <- 0.98 L <- 1 df <- VasicekZCBOption(kappa, mu, sigma, r0, t, T, S, K, L)
オプション価格は次のようになる
> df moneyness callprice putprice 1 0.9888111 0.01108998 4.579769e-13
プットコールパリティが成立することを確認
> L * VasicekZCBond(kappa, mu, sigma, r0, t, S)$price + df$putprice [1] 0.9911608 > K * VasicekZCBond(kappa, mu, sigma, r0, t, T)$price + df$callprice [1] 0.9911608
マネネスが１となるあたりでオプション価格変化をプロット
library(dplyr) library(tidyr) library(ggplot2) library(ggthemes) K <- seq(0.975, 1.0075, by = 0.001) df <- VasicekZCBOption(kappa, mu, sigma, r0, t, T, S, K, L) df %>% gather(type, price, -moneyness) %>% ggplot(aes(x = moneyness, y = price, colour = type)) + geom_line(size = 1.2) + theme_economist() + ggtitle("option price") + scale_x_continuous(breaks = seq(0.98, 1.02, by = 0.01))
エコノミスト誌風な色が付けられるscale_colour_economistもあったが２本のみだと良い感じではなかったので使用しなかった
ゼロクーポン債に対するヨーロピアンオプションで、1ファクターVasicekモデルでの評価なんて実際には使われていないだろうけれども、これを基本にしてスワップションなどのより複雑なデリバティブを考えていくにはよい勉強になる
Vasicekモデルの債券価格とイールド
金利のモデル、デリバティブあたりをまとめた便利なRのパッケージってない…？ 少なくとも、有名どころのパッケージを聞かないので、多くの人に安定的に使われているものはなさそう
未だ使い慣れないggplot2の練習も含めて、Vasicekモデルにおける債券価格とイールドについて
# zero coupon bond price and yield under Vasicek model # # Args: # kappa : speed of reversion # mu : long term mean level # sigma : instaneous volatility # r0 : current short rate # t : current time # T : bond maturity # # Return: # data.frame of maturity, bond price and yield # VasicekZCBond <- function(kappa, mu, sigma, r0, t, T) { tau <- T - t B <- (1 - exp(-kappa * tau)) / kappa A <- exp((mu - sigma^2 / (2 * kappa^2)) * (B - tau) - sigma ^ 2 / (4 * kappa) * B^2) P <- A * exp(-B * r0) Y <- ifelse(tau != 0, -log(P) / tau, r0) data.frame(maturity = tau, price = P, yield = Y) }
現在の日本やドイツのイールドカーブにフィットするようなパラメータは大体の教科書に出てくるパラメータとはだいぶかけ離れている
Vasicekモデルなので、お遊びで現在のスポットレートをマイナスにセットしても面白いかもしれない
kappa <- 0.2 mu <- 0.015 sigma <- 0.003 r0 <- -0.001 t <- 0 T <- c(1/12, 1/4, 1/2, 1:10, 15, 20, 30) df <- VasicekZCBond(kappa, mu, sigma, r0, t, T)
債券価格とイールドは次のようになる
> df maturity price yield 1 0.08333333 1.0000723 -0.0008674146 2 0.25000000 1.0001517 -0.0006066745 3 0.50000000 1.0001132 -0.0002263613 4 1.00000000 0.9995030 0.0004971657 5 2.00000000 0.9963899 0.0018083087 6 3.00000000 0.9911608 0.0029594929 7 4.00000000 0.9842342 0.0039728389 8 5.00000000 0.9759579 0.0048671608 9 6.00000000 0.9666189 0.0056584930 10 7.00000000 0.9564529 0.0063605269 11 8.00000000 0.9456528 0.0069849753 12 9.00000000 0.9343755 0.0075418761 13 10.00000000 0.9227486 0.0080398472 14 15.00000000 0.8623588 0.0098722601 15 20.00000000 0.8024879 0.0110019245 16 30.00000000 0.6923496 0.0122554754
イールドカーブをプロットする
library(ggplot2) library(ggthemes) ggplot(df, aes(x = maturity, y = yield)) + geom_line(size = 1.2) + theme_economist() + ggtitle("normal yield curve") + xlab("time to maturity")
ここでは、Vasicekモデルで表現できるハンプ、逆イールドについてもプロット
kappa <- 0.15 mu <- 0.0012 sigma <- 0.003 r0 <- 0.001 df <- VasicekZCBond(kappa, mu, sigma, r0, t, T) ggplot(df, aes(x = maturity, y = yield)) + geom_line(size = 1.2) + theme_economist() + ggtitle("humped yield curve") + xlab("time to maturity")
kappa <- 0.15 mu <- 0.008 sigma <- 0.003 r0 <- 0.009 df <- VasicekZCBond(kappa, mu, sigma, r0, t, T) ggplot(df, aes(x = maturity, y = yield)) + geom_line(size = 1.2) + theme_economist() + ggtitle("inverted yield curve") + xlab("time to maturity")
Vasicekモデルは実務に耐え得るモデルではないものの、いくつかの短期金利モデルを比較するような論文には登場することは多いので、これぐらいはパッケージ化されていても良い気もしますが…
rvestでスクレイプした時の文字化けへの対処
「rvestで日本語のページをスクレイプしても自動でエンコーディングを推論しているから問題ないよ」的な記事を見かけた
日銀のページはcharset=Shift_JISだが、charset=utf-8のページでも文字化けするケースがあったので、以下を試した
対象はBloombergの国債金利ページ。普通にスクレイプして、２つ目のテーブルだけ取り出すと文字化け
> library(rvest) > url <- "http://www.bloomberg.co.jp/markets/rates.html" > HTML <- html(url) > HTML %>% html_table() %>% .[[2]] a<U+009C><U+009F>e<U+0096><U+0093> a<U+0082>￣a<U+0083><U+00BC>a<U+0083><U+009D>a<U+0083>3 a<U+0084><U+009F>e<U+0082><U+0084>a<U+0097>\\ 1 1a1´ 0.1 01/15/2016 2 2a1´ 0.1 01/15/2017 3 3a1´ 1.5 12/20/2017 4 4a1´ 1.3 12/20/2018 5 5a1´ 0.1 12/20/2019 6 6a1´ 1.2 12/20/2020 7 7a1´ 0.9 03/20/2022 8 8a1´ 0.7 12/20/2022 9 9a1´ 0.6 12/20/2023 10 10a1´ 0.3 12/20/2024 11 15a1´ 2.1 12/20/2029 12 20a1´ 1.2 12/20/2034 13 30a1´ 1.5 12/20/2044 a<U+00BE>!a<U+00A0><U+00BC> a<U+0088>ca<U+009B><U+009E>a<U+0082><U+008A> a<U+00BE>!a<U+00A0><U+00BC>a<U+00A4><U+0089>a<U+008C><U+0096>c<U+008E><U+0087> 1 100.12 -0.02 -0.006 2 100.28 -0.04 0.006 3 104.44 -0.03 -0.020 4 105.09 0.00 -0.078 5 100.40 0.02 -0.074 6 106.95 0.02 -0.082 7 105.75 0.09 -0.278 8 104.42 0.13 -0.284 9 103.61 0.19 -0.302 10 100.54 0.24 -0.276 11 121.89 0.52 -0.386 12 105.09 0.90 -0.138 13 108.98 1.10 -0.428 a<U+0088>ca<U+009B><U+009E>a<U+0082><U+008A>a<U+00A4><U+0089>a<U+008C><U+0096>c<U+008E><U+0087> a<U+009B>´a<U+0096>°a<U+0099><U+0082>e<U+0096><U+0093> 1 0.004 00:00 2 -0.004 00:00 3 0.001 00:00 4 0.016 00:00 5 0.015 00:00 6 0.011 00:00 7 0.036 00:00 8 0.033 00:00 9 0.033 00:00 10 0.028 00:00 11 0.022 00:00 12 0.008 00:00 13 0.017 00:00
readLinesを使って、ページを読めば解消できた（ただし、警告は出た…）
> HTML <- html(readLines(url, encoding = "UTF-8")) 警告メッセージ: 1: In if (grepl("^http", x)) { : 条件が長さが 2 以上なので、最初の 1 つだけが使われます 2: In if (grepl("<|>", x)) { : 条件が長さが 2 以上なので、最初の 1 つだけが使われます > HTML %>% html_table() %>% .[[2]] 期間 クーポン 償還日 価格 利回り 価格変化率 利回り変化率 更新時間 1 1年 0.1 01/15/2016 100.12 -0.02 -0.006 0.005 21:00 2 2年 0.1 01/15/2017 100.28 -0.04 0.006 -0.004 21:00 3 3年 1.5 12/20/2017 104.44 -0.03 -0.020 0.001 21:00 4 4年 1.3 12/20/2018 105.09 0.00 -0.078 0.016 21:00 5 5年 0.1 12/20/2019 100.40 0.02 -0.074 0.015 21:00 6 6年 1.2 12/20/2020 106.95 0.02 -0.082 0.011 21:00 7 7年 0.9 03/20/2022 105.75 0.09 -0.278 0.036 21:00 8 8年 0.7 12/20/2022 104.42 0.13 -0.284 0.033 21:00 9 9年 0.6 12/20/2023 103.61 0.19 -0.302 0.033 21:00 10 10年 0.3 12/20/2024 100.54 0.24 -0.276 0.028 21:00 11 15年 2.1 12/20/2029 121.89 0.52 -0.386 0.022 21:00 12 20年 1.2 12/20/2034 105.09 0.90 -0.138 0.008 21:00 13 30年 1.5 12/20/2044 108.98 1.10 -0.428 0.017 21:00
html関数の中にもencodingのオプションがあるが、何か上手くいかなかった。説明に書かれているようにstringiパッケージを使うべきなのかもしれない。最適な方法かどうかはわからないけれども、こちらは大変簡単なのでメモしておく
Rの金利期間構造パッケージ termstrc
ファイナンス系パッケージ探訪。fBonds、YieldCurveパッケージに続きtermstrcパッケージを使ってみたメモ
Journal of Statistical SoftwareからZero-Coupon Yield Curve Estimation with the Package termstrcが出ており、パッケージの詳しい説明が書かれているので、本気で使いたい場合は一読が必要
fBonds、YieldCurveと比較して、クーポン債にも適用できる、Cubicスプラインによる推定ができる、adjusted Svenssonモデルも適用できる、Rcpp使っている、その他、有用な関数が入っている、ということがtermstrcの特徴になっている
まずは、パッケージを読み込んで、中身を確認
library(termstrc)
パッケージマニュアルには、関数がかなり細かく説明されているが、逆に何を使ったらよいかわからなくなっている
基本的な使い方としては、データがゼロクーポン債（zeroyields）なのかクーポン債（couponbonds）なのかを指定して、estim_nssでモデルのパラメータを推定する
> data(gobbonds) > class(govbonds) [1] "couponbonds" > str(govbonds$GERMANY) List of 8 $ ISIN : chr [1:52] "DE0001141414" "DE0001137131" "DE0001141422" "DE0001137149" ... $ MATURITYDATE: Date[1:52], format: "2008-02-15" "2008-03-14" ... $ ISSUEDATE : Date[1:52], format: "2002-08-14" "2006-03-08" ... $ COUPONRATE : num [1:52] 0.0425 0.03 0.03 0.0325 0.0413 ... $ PRICE : num [1:52] 100 99.9 99.8 99.8 100.1 ... $ ACCRUED : num [1:52] 4.09 2.66 2.43 2.07 2.39 ... $ CASHFLOWS :List of 3 ..$ ISIN: chr [1:384] "DE0001141414" "DE0001137131" "DE0001141422" "DE0001137149" ... ..$ CF : num [1:384] 104 103 103 103 104 ... ..$ DATE: Date[1:384], format: "2008-02-15" "2008-03-14" ... $ TODAY : Date[1:1], format: "2008-01-30" > ns_res <- estim_nss(govbonds, c("GERMANY", "AUSTRIA", "FRANCE"), + matrange = c(0, 30), method = "ns", tauconstr = list(c(0.2, 5, 0.1), c(0.2, 5, 0.1), c(0.2, 5, 0.1))) [1] "Searching startparameters for GERMANY" beta0 beta1 beta2 tau1 5.132836 -1.274357 -3.208435 2.700100 [1] "Searching startparameters for AUSTRIA" beta0 beta1 beta2 tau1 5.050193 -1.327244 -2.629411 2.500100 [1] "Searching startparameters for FRANCE" beta0 beta1 beta2 tau1 5.108886 -1.217795 -3.068065 2.500100 > summary(ns_res) --------------------------------------------------- Goodness of fit: --------------------------------------------------- GERMANY AUSTRIA FRANCE RMSE-Prices 0.3582276 0.1801092 0.2214637 AABSE-Prices 0.1992019 0.1224709 0.1182047 RMSE-Yields (in %) 0.0847062 0.0185987 0.0392355 AABSE-Yields (in %) 0.0498615 0.0155659 0.0275024 --------------------------------------------------- Startparameters: --------------------------------------------------- beta0 beta1 beta2 tau1 GERMANY 5.13284 -1.27436 -3.20844 2.70010 AUSTRIA 5.05019 -1.32724 -2.62941 2.50010 FRANCE 5.10889 -1.21779 -3.06807 2.50010 --------------------------------------------------- Convergence information: --------------------------------------------------- optim() convergence info GERMANY 0 AUSTRIA 0 FRANCE 0 optim() solver message GERMANY NULL AUSTRIA NULL FRANCE NULL
クーポン債の場合には、データセットのように関数が適用できる形にデータを直したりするのが面倒かもしれない
パラメータを推定して、その時系列推移を見るという場合に結構使えるかもしれない。 ustycパッケージを使って、米国債イールドデータに適用してみる
library(ustyc) maturities <- c(1/12, 1/4, 1/2, 1, 2, 3, 5, 7, 10, 20, 30) yc <- getYieldCurve(year = "2014") yields <- as.matrix(yc$df[-12]) dates <- rownames(yields) datazeroyields <- zeroyields(maturities, yields, dates)
一応、これだけで、3Dプロットが出力できる
> datazeroyields This is a data set of zero-coupon yields. Maturities range from 0.0833333333333333 to 30 years. There are 250 observations between 2014-01-02 and 2014-12-31 . > class(datazeroyields) [1] "zeroyields" > plot(datazeroyields)
Nelson-Siegelモデルを当てはめる
ns_res <- estim_nss(datazeroyields, "ns", tauconstr = c(0.2, 6, 0.1))
まとめて、出力を見てみる
> ns_res --------------------------------------------------- Estimated Nelson/Siegel parameters: --------------------------------------------------- Number of oberservations: 250 [[1]] beta_0 beta_1 beta_2 tau_1 Min. :2.826 Min. :-4.350 Min. :-5.656 Min. :1.024 1st Qu.:3.358 1st Qu.:-3.853 1st Qu.:-4.890 1st Qu.:1.284 Median :3.698 Median :-3.634 Median :-4.343 Median :1.381 Mean :3.674 Mean :-3.601 Mean :-4.392 Mean :1.388 3rd Qu.:3.953 3rd Qu.:-3.313 3rd Qu.:-3.988 3rd Qu.:1.485 Max. :4.476 Max. :-2.801 Max. :-2.575 Max. :1.668 > summary(ns_res) --------------------------------------------------- Goodness of fit: --------------------------------------------------- [,1] RMSE-Yields (in %) 0.0334253 AABSE-Yields (in %) 0.0277030 --------------------------------------------------- Convergence information from optim (): --------------------------------------------------- [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 0 0 0 0 0 0 0 0 0 0 [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [1,] 0 0 0 0 0 0 0 0 0 … > head(param(ns_res)[[1]]) beta_0 beta_1 beta_2 tau_1 [1,] 4.471780 -4.345616 -5.655850 1.545545 [2,] 4.476129 -4.349701 -5.596621 1.545772 [3,] 4.447632 -4.340263 -5.538249 1.561728 [4,] 4.421710 -4.323720 -5.432554 1.570774 [5,] 4.420168 -4.326042 -5.409182 1.511038 [6,] 4.380350 -4.294240 -5.321058 1.522402
Nelson-Siegelモデルならば、パッケージの使用例に書かれているtauconstrを指定すればよいが、Svensson、adjusted Svenssonだとoptimのところで収束しない、もしくは、パラメータの推定値がおかしな場合があるので、変えた方がよい
パラメータの時系列推移と各ファクターの寄与をプロット
> plot(param(ns_res)) > fcontrib(param(ns_res), index = 1, m = 1:30, method = "ns")
一年を通して、長期金利が低下傾向だったので、このような感じになる
満期が長くなるにつれて、長期水準以外のファクターの影響は小さくなっていく
関数が適用できる形にデータを持っていけば、いろいろなモデルに当てはめられる、良い感じの出力やプロットが得られるところがtermstrcのメリット。簡単に使えるシンプルさだったら、YieldCurveパッケージ
ggplot2を使いたい、もっと自由にいろいろデータを加工、適用したいという時にはtermstrcは少し面倒。このようなパッケージを参考に、自分で関数を作りながらやったほうが勉強にもなるし、分析はしやすい…
Rの金利期間構造パッケージ fBonds
ほとんど使ったことがなかったRのファイナンス系のパッケージ（CRAN Task View: Empirical Finance）を探訪する
前回で利用したYieldCurveパッケージ以外にも金利の期間構造に関するパッケージが複数あることを知ったので、それを調査したメモ
今回、調査したのは、fBondsとYieldCurveパッケージ
まずは、パッケージの読み込み
library(fBonds) library(YieldCurve)
fBondsパッケージは、xtsと同様にRのcoreパッケージということのようだが、Nelson-SiegelモデルとSvenssonモデルへ当てはめる関数の二つのみで、coreパッケージというにしてはあっさりしている
Nelson-Siegelモデルをイールドカーブに当てはめる関数で比較
データは、財務省の国債金利情報から持ってきた
maturity <- c(1:10, 15, 20, 25, 30, 40) yield <- c(-0.013, -0.020, -0.014, -0.008, 0.016, 0.017, 0.073, 0.148, 0.216, 0.284, 0.589, 0.938, 1.063, 1.150, 1.260)
マイナス金利というモデルが作られた時には想定していない事態かもしれないが、このまま推定
par1 <- fBonds::NelsonSiegel(yield, maturity, doplot = TRUE) par2 <- YieldCurve::Nelson.Siegel(yield, maturity)
同じNelson-Siegelモデルへの当てはめでも、最適化に使う関数が違って、fBondsの方はnlminbを使い、YieldCurveの方はoptimizeでパラメータlambdaのみ一旦推定して、lmで線形モデルに当てはめ、残りのパラメータを推定している
このデータから推定したパラメータは、以下の通り
> par1$par beta0 beta1 beta2 tau1 1.805390 -1.659496 -2.991222 4.242192 > par2 beta_0 beta_1 beta_2 lambda [1,] 1.795689 -1.642322 -3.009073 0.239118
fBondsはフォワードレート、YieldCurveはスポットレートの形でモデルのパラメータ定義がされているので、ここでのパラメータの関係は、lambda = 1 / tau1となっている
比較のために、fBondsで推定した結果を変換して、格納しておく
> par1$par.t <- c(par1$par[-4], 1 / par1$par[4]) > names(par1$par.t) <- c("beta_0", "beta_1", "beta_2", "lambda") > par1$par.t beta_0 beta_1 beta_2 lambda 1.8053901 -1.6594964 -2.9912217 0.2357272
次は、この推定したパラメータからモデルイールドの値を計算
fBondsパッケージには、推定したパラメータから、Nelson-Siegelモデルのイールドを計算する関数は定義されていない
YieldCurveパッケージには、NSratesという関数が入っているが、時系列オブジェクトを要求してくるので、今回のような一日分のデータから、モデルイールドを求めたい場合には、一度xts等に変換が必要で面倒
ということで、Nelson-Siegelモデルのイールドを計算する関数を書いて比較
NSModelRates <- function(par, tau) { par[1] + par[2] * (1 - exp(-par[4] * tau)) / (par[4] * tau) + par[3] * ((1 - exp(-par[4] * tau)) / (par[4] * tau) - exp(-par[4] * tau)) }
それぞれのパッケージの関数で推定したパラメータを入れる
> NSModelRates(par1$par.t, maturity) [1] 0.025227873 -0.035976844 -0.053874588 -0.040757689 -0.005934127 [6] 0.043589340 0.102575961 0.167146462 0.234459168 0.302462379 [11] 0.615565858 0.854583921 1.026647584 1.150845555 1.312439110 > NSModelRates(par2, maturity) [1] 0.027722244 -0.036202635 -0.055344088 -0.042543069 -0.007486335 [6] 0.042562722 0.102197788 0.167431439 0.235356597 0.303885415 [11] 0.618089435 0.856425926 1.027194182 1.150083261 1.309627120
fBonds::NelsonSiegelを実行した際に出てくるプロットにYieldCurveで推定されたイールドを重ねる
lines(maturity, NSModelRates(par2, maturity), lty = 2, col = 2)
青がfBonds、赤がYieldCurveで推定したもの
両方のパッケージを使って推定したイールドが重なることが確認できる
イールドカーブの変動は、水準、傾き、曲率でほとんどが説明できると言われており、Nelson-Siegelモデルでは、それぞれの係数がそれらに対応している。これらはファクターローディングと呼ばれる
簡単な関数なので、自分で書いても良いが面倒な場合は、YieldCurveの説明には書かれていないが名前空間には定義されているので、それを使うことができる
> YieldCurve::<tab> YieldCurve::.beta1Forward YieldCurve::.beta1Spot YieldCurve::.beta2Forward YieldCurve::.beta2Spot YieldCurve::.factorBeta1 YieldCurve::.factorBeta2 YieldCurve::.NS.estimator YieldCurve::.NSS.estimator YieldCurve::Nelson.Siegel YieldCurve::NSrates YieldCurve::Srates YieldCurve::Svensson
YieldCurve::.factorBeta1とYieldCurve::.factorBeta2を使えば、各満期におけるファクターローディングが得られる
beta0.fl <- rep(1, length(maturity)) beta1.fl <- YieldCurve::.factorBeta1(par2[4], maturity) beta2.fl <- YieldCurve::.factorBeta2(par2[4], maturity) matplot(maturity, cbind(beta0.fl, beta1.fl, beta2.fl), type = "l", lwd = 2, col = c(1, 2, 4)) legend(30, 0.8, c(expression(beta[0]), expression(beta[1]), expression(beta[2])), lty = 1:3, lwd = 2, col = c(1, 2, 4)) grid()
以上をまとめると、fBondsパッケージは使う必要はあまり感じない…
YieldCurveパッケージは、もう少しいろいろと充実していただけると…
Rから米国債イールドデータ収集
ジェネリック国債のデータを確認しようとしたら、ustycなるRのパッケージを見つけたので、簡単な使い方をメモ
ustycは、米国財務省が提供しているDaily Treasury Yield Curve Ratesのデータをまとめて取ってくるという関数が入ったパッケージ。ustycはもちろん、US Treasury Yield Curveの略
（Bloombergのサイトから、GJGB3M, GJGB6M, …, GJGB10, …, GJGB30と各満期の日本のジェネリック国債データを確認したかったが、サイトからでは確認できなくなった？）
使い方はパッケージを読み込んで、関数を実行して多少待つだけ
library(ustyc) yc <- getYieldCurve()
特定の年、月だけ取りたい場合は、getYieldCurve(year = "2014", month = "12")のようにして指定する
yc$df以外は必要ないので、これだけCSVファイルなどに書き出しておく
write.csv(yc$df, file = "TreasuryYield.csv")
利用する時に、xtsに変換してやる
library(xts) yc <- read.csv("TreasuryYield.csv") yc.xts <- xts(yc[-1], order.by = as.Date(yc[, 1]))
BC_30YEARDISPLAYが何を指すのかがよくわからなかったが…
> tail(yc.xts) BC_1MONTH BC_3MONTH BC_6MONTH BC_1YEAR BC_2YEAR BC_3YEAR BC_5YEAR 2014-12-24 0.01 0.01 0.14 0.26 0.73 1.18 1.76 2014-12-26 0.01 0.01 0.10 0.26 0.73 1.19 1.75 2014-12-29 0.01 0.03 0.12 0.25 0.72 1.14 1.72 2014-12-30 0.03 0.03 0.12 0.23 0.69 1.11 1.68 2014-12-31 0.03 0.04 0.12 0.25 0.67 1.10 1.65 2015-01-02 0.02 0.02 0.11 0.25 0.66 1.07 1.61 BC_7YEAR BC_10YEAR BC_20YEAR BC_30YEAR BC_30YEARDISPLAY 2014-12-24 2.09 2.27 2.56 2.83 2.83 2014-12-26 2.07 2.25 2.54 2.81 2.81 2014-12-29 2.02 2.22 2.51 2.78 2.78 2014-12-30 2.00 2.20 2.49 2.76 2.76 2014-12-31 1.97 2.17 2.47 2.75 2.75 2015-01-02 1.92 2.12 2.41 2.69 2.69 > plot(yc.xts[, 9], main = "BC_10YEAR", ylab = "yield (%)")
YieldCurveパッケージを使ったNelson-Siegelモデル、Svenssonモデルへの当てはめ
library(YieldCurve) maturity <- c(1/12, 1/4, 1/2, 1, 2, 3, 5, 7, 10, 20, 30) obs.rate <- yc.xts["2015-01-02", -12] NS.par <- Nelson.Siegel(obs.rate, maturity) NS.rate <- NSrates(NS.par, maturity) SV.par <- Svensson(obs.rate, maturity) SV.rate <- Srates(SV.par, maturity, "Spot") plot(maturity, obs.rate, main = "Fitting yield curve (2015-01-02)", type = "o", ylab = "yield (%)") lines(maturity, NS.rate, col = 2) lines(maturity, SV.rate, col = 4) legend("bottomright", legend = c("observed yield curve", "Nelson-Siegel model", "Svensson model"), col = c(1, 2, 4), lty = 1) grid()
rvestを使って、投資家の傾向を探る
前回に続き、rvestを使ってみた記録になります。
ヤフーファイナンスの個別銘柄の詳細ページから「この銘柄を見た人はこんな銘柄も見ています」の情報を取得して、何か投資家の傾向が見られないかを探ってみます。
今回も時価総額上位100銘柄を対象とします。前回のプログラムをひと通り流して、時価総額上位100銘柄の証券コードを取得しておきます。
証券コードから個別銘柄の詳細ページをスクレイプします。
detail.HTML <- list() for (i in 1:100) { detail.URL <- paste0("http://stocks.finance.yahoo.co.jp/stocks/detail/?code=", stock.code[i] ,".t") detail.HTML[[i]] <- html(detail.URL) }
分析中に何度もスクレイプすると良くないので、ページの情報はdetail.HTMLに一旦入れておきます。
「この銘柄を見た人はこんな銘柄も見ています」は、html_nodes("ol.affirank li a")で指定すれば、上の方にある1位から5位と下の方にある6位から15位の銘柄が取り出せます。「もっと見る」の部分だけ不要なので、取り除きます。
library(dplyr) affinity.list <- list() for (i in 1:100) { affinity.list[[i]] <- detail.HTML[[i]] %>% html_nodes("ol.affirank li a") %>% html_text() %>% .[-which(. == "もっと見る")] }
これ以下の実行列は、2014年12月30日の市場終了時点
トヨタ自動車とファナックの「この銘柄を見た人はこんな銘柄も見ています」
> stock.code[1] [1] "7203" > affinity.list[[1]] [1] "トヨタグループ株式ファンド" [2] "トヨタグループ世界債券ファンド(年2回)" [3] "デンソー" [4] "ホンダ" [5] "DCトヨタグループ株式ファンド" [6] "豊田織" [7] "日産自" [8] "富士重" [9] "アイシン精" [10] "日野自" [11] "ソフバンク" [12] "マツダ" [13] "トヨタ紡" [14] "DIAMストラテジックJ-REITファンド" [15] "ファインシ" > stock.code[14] [1] "6954" > affinity.list[[14]] [1] "キーエンス" [2] "村田製" [3] "京セラ" [4] "ＳＭＣ" [5] "Ｆリテイリ" [6] "日電産" [7] "信越化" [8] "ローム" [9] "シマノ" [10] "東エレク" [11] "いちよし証" [12] "ナブテスコ" [13] "安川電" [14] "グリコ" [15] "ヒロセ電"
関連ファンドやグループ企業、そして、同業他社が多く見られている
「この銘柄を見た人はこんな銘柄も見ています」によく登場する銘柄
> do.call("rbind", affinity.list) %>% + table() %>% + sort(dec = TRUE) %>% + head(50) ＴＤＫ ダイキン工 京セラ 信越化 村田製 10 10 10 10 10 武田薬 ファナック ローム エーザイ コマツ 10 9 9 8 8 サントリ食 日電産 Ｆリテイリ ダイハツ工 ＫＤＤＩ 8 8 7 7 6 アステラ薬 オムロン カルビー キーエンス テルモ 6 6 6 6 6 トヨタ ブリヂスト ホンダ マツダ 塩野義薬 6 6 6 6 6 三井物 住友電設 大塚ＨＬＤ 大日本住友 第一三共 6 6 6 6 6 東芝 日産自 日東電 日野自 豊田通商 6 6 6 6 6 ＮＥＣ ＮＴＴ ＮＴＴ都市 アイシン精 カカクコム 5 5 5 5 5 キヤノン スズキ ソニーＦＨ ライオン 科研薬 5 5 5 5 5 三井住友Ｆ 三菱ＵＦＪ 三菱自 三菱商 小野薬 5 5 5 5 5
順位に関係する情報は落ちてしまいますが、アソシエーション分析してみます。
library(arules) library(arulesViz) affinity.tran <- as(affinity.list, "transactions")
変換されたものを確認します。
> affinity.tran transactions in sparse format with 100 transactions (rows) and 801 items (columns)
> affinity.rules <- apriori(affinity.tran, parameter = list(support = 0.03, confidence = 0.8)) > affinity.rules set of 1102 rules
1102ルール見つけることができました。supportが0.04以上のものだけに絞ります。
subset.rules <- subset(sort(affinity.rules, by = "support"), support >= 0.04)
一部だけ確認します。
> inspect(subset.rules[1:30]) lhs rhs support confidence lift 1 {小野薬} => {塩野義薬} 0.05 1.0000000 16.666667 2 {塩野義薬} => {小野薬} 0.05 0.8333333 16.666667 3 {富士重} => {日産自} 0.05 1.0000000 16.666667 4 {日産自} => {富士重} 0.05 0.8333333 16.666667 5 {スズキ} => {ダイハツ工} 0.05 1.0000000 14.285714 6 {日野自} => {マツダ} 0.05 0.8333333 13.888889 7 {マツダ} => {日野自} 0.05 0.8333333 13.888889 8 {日野自} => {ダイハツ工} 0.05 0.8333333 11.904762 9 {第一三共} => {塩野義薬} 0.05 0.8333333 13.888889 10 {塩野義薬} => {第一三共} 0.05 0.8333333 13.888889 11 {第一三共} => {アステラ薬} 0.05 0.8333333 13.888889 12 {アステラ薬} => {第一三共} 0.05 0.8333333 13.888889 13 {第一三共} => {エーザイ} 0.05 0.8333333 10.416667 14 {塩野義薬} => {アステラ薬} 0.05 0.8333333 13.888889 15 {アステラ薬} => {塩野義薬} 0.05 0.8333333 13.888889 16 {塩野義薬} => {エーザイ} 0.05 0.8333333 10.416667 17 {アステラ薬} => {エーザイ} 0.05 0.8333333 10.416667 18 {マツダ} => {ダイハツ工} 0.05 0.8333333 11.904762 19 {ＴＤＫ, ローム} => {村田製} 0.05 0.8333333 8.333333 20 {ローム, 村田製} => {ＴＤＫ} 0.05 0.8333333 8.333333 21 {ＴＤＫ, 村田製} => {ローム} 0.05 0.8333333 9.259259 22 {ＴＤＫ, ローム} => {京セラ} 0.05 0.8333333 8.333333 23 {ローム, 京セラ} => {ＴＤＫ} 0.05 0.8333333 8.333333 24 {ＴＤＫ, ローム} => {信越化} 0.05 0.8333333 8.333333 25 {ローム, 信越化} => {ＴＤＫ} 0.05 0.8333333 8.333333 26 {ＴＤＫ, 信越化} => {ローム} 0.05 0.8333333 9.259259 27 {ローム, 村田製} => {京セラ} 0.05 0.8333333 8.333333 28 {ローム, 京セラ} => {村田製} 0.05 0.8333333 8.333333 29 {ローム, 村田製} => {信越化} 0.05 0.8333333 8.333333 30 {ローム, 信越化} => {村田製} 0.05 0.8333333 8.333333
最後にプロットします。
plot(subset.rules, method = "graph", control = list(type="items"))
医薬品、輸送用機器、電気機器・半導体関連が大きなクラスタのようになりました。業種間をつなげるような銘柄は見当たらず…明らかに配当利回りの高い銘柄どうしを比較しているわけでもなさそうです
