ggplot2で日本国債の時系列金利水没マップ
WBSのコメンテーターでおなじみのみずほ総研の高田さんが「世界の金利「水没」マップ、金融機関はどう生き残るのか」というレポートを1月27日に出されていた
ここ数日は日本国債の短期金利がマイナスという状態を脱して来ているが、時系列での金利水没具合に興味があったのでggplot2で可視化してみた
パッケージをひと通り読み込んでおく
library(dplyr) library(tidyr) library(ggplot2) library(ggthemes)
金利データは財務省の金利情報から取得する。日本語版は和暦なので、英語版ページ（Interest Rate : Ministry of Finance）から持ってきた方が良い。データは現在の月のデータと2010年からのデータの二つを使う
jgbcme <- read.csv("http://www.mof.go.jp/english/jgbs/reference/interest_rate/jgbcme.csv") jgbcme_2010 <- read.csv("http://www.mof.go.jp/english/jgbs/reference/interest_rate/historical/jgbcme_2010.csv") JGB_rate <- rbind(jgbcme_2010, jgbcme) JGB_rate$Date <- as.Date(JGB_rate$Date)
データは2010年の1月4日から2015年の2月5日まで
> head(JGB_rate) Date X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 X15 X20 X25 X30 X40 1 2010-01-04 0.131 0.147 0.238 0.383 0.495 0.662 0.828 1.009 1.184 1.322 1.810 2.110 2.236 2.266 2.311 2 2010-01-05 0.136 0.148 0.234 0.376 0.491 0.660 0.834 1.015 1.188 1.329 1.819 2.119 2.244 2.275 2.319 3 2010-01-06 0.135 0.154 0.237 0.378 0.502 0.675 0.851 1.029 1.198 1.337 1.828 2.128 2.252 2.280 2.327 4 2010-01-07 0.137 0.154 0.243 0.388 0.506 0.689 0.872 1.045 1.208 1.345 1.847 2.144 2.269 2.294 2.340 5 2010-01-08 0.136 0.164 0.253 0.398 0.521 0.707 0.892 1.064 1.222 1.359 1.855 2.147 2.270 2.297 2.346 6 2010-01-12 0.131 0.164 0.254 0.396 0.519 0.697 0.884 1.057 1.220 1.355 1.854 2.144 2.266 2.293 2.346 > tail(JGB_rate) Date X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 X15 X20 X25 X30 X40 1245 2015-01-29 0.005 0.005 0.021 0.035 0.051 0.057 0.099 0.171 0.236 0.293 0.666 1.064 1.230 1.313 1.428 1246 2015-01-30 0.009 0.010 0.023 0.032 0.046 0.045 0.081 0.156 0.225 0.288 0.656 1.054 1.216 1.299 1.415 1247 2015-02-02 0.010 0.010 0.029 0.045 0.056 0.055 0.093 0.166 0.226 0.293 0.646 1.045 1.206 1.295 1.406 1248 2015-02-03 0.020 0.020 0.044 0.071 0.096 0.104 0.158 0.232 0.297 0.363 0.709 1.102 1.266 1.362 1.473 1249 2015-02-04 0.019 0.030 0.050 0.081 0.100 0.109 0.163 0.245 0.314 0.382 0.761 1.167 1.338 1.437 1.541 1250 2015-02-05 0.019 0.043 0.065 0.091 0.105 0.106 0.158 0.238 0.299 0.362 0.732 1.129 1.275 1.372 1.483
gatherしてから、geom_tileでプロット
JGB_rate_gather <- JGB_rate %>% gather(maturity, yield, -Date) ggplot(JGB_rate_gather, aes(x = Date, y = maturity, fill = yield)) + geom_tile() + theme_tufte() + scale_fill_gradientn(colours=c("black", "blue", "white"), limits = c(-0.05, 1), na.value = "white")
ggplot2が賢いためか、勝手に営業日（Business Day）以外の部分もプロットしてしまうようだ。そこで、 営業日スケールに直せるbdscaleパッケージを使う
library(bdscale) library(scales) ggplot(JGB_rate_gather, aes(x = Date, y = maturity, fill = yield)) + geom_tile() + theme_tufte() + scale_fill_gradientn(colours=c("black", "blue", "white"), limits = c(-0.05, 1), na.value = "white") + scale_x_bd(business.dates = JGB_rate$Date, max.major.breaks = 10, labels = date_format("%Y"))
2014年末からは青黒い部分が広がっていて、過去に比べ金利低下傾向が一段と強いことがわかる。このグラフではどの満期までが水没しかかっているかは確認しやすいが、実際の金利がどのように推移しているのかはわかりにくい
そこで、時系列プロットにgeom_ribbbonで単純に水没域を重ねてみる
ggplot(JGB_rate_gather, aes(x = Date, y = yield, group = maturity)) + geom_line(aes(colour = maturity)) + theme_tufte() + geom_ribbon(fill = "black", alpha = 0.05, ymin = -1, ymax = 0) + geom_ribbon(fill = "blue", alpha = 0.05, ymin = 0, ymax = 0.5) + geom_ribbon(fill = "lightblue", alpha = 0.05, ymin = 0.5, ymax = 1)
10年満期までの金利が完全に水没、水没しかかっているのがこちらでもはっきりわかる。このような状態で金融機関はちゃんと息できているんでしょうか？
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()
