Rの金利期間構造パッケージ 続SmithWilsonYieldCurve
前回に続き、SmithWilsonYieldCurveパッケージを使ったメモ
ノルウェーの金融監督庁（正式日本語表記わからず）から出ていたA Technical Note on the Smith-Wilson MethodのWorked examplesをやってみる
想定している状況を定義して、fFitSmithWilsonYieldCurveToInstrumentsに適用する
library(SmithWilsonYieldCurve) alpha <- 0.1 ufr <- log(1 + 0.042) Type <- rep("SWAP", 4) Tenor <- c(1, 2, 3, 5) Rate <- c(0.01, 0.02, 0.026, 0.034) Frequency <- rep(1, 4) dfInstruments <- data.frame(Type, Tenor, Rate, Frequency) Curve <- fFitSmithWilsonYieldCurveToInstruments(dfInstruments, ufr, alpha)
結果を確認する
> Curve $P function (t) { fBase(t) + t(KernelWeights) %*% fCompoundKernel(t) } <environment: 0x0000000010ddc478> $xi [,1] [1,] 57.790688 [2,] -33.507208 [3,] 11.396473 [4,] -5.466968 $K function (t) { CashflowMatrix %*% fKernel(t, TimesVector) } <environment: 0x0000000010ddc478> attr(,"class") [1] "SmithWilsonYieldCurve" "YieldCurve"
xiの値が書かれているものと違う…おそらく、書き間違え
不安なので、計算過程を確認する。まず、Wilson関数を計算する
t <- 1:5 u <- 1:5 g <- expand.grid(t, u) w <- sapply(1:nrow(g), function(i) fWilson(g$Var1[i], g$Var2[i], ufr, alpha)) W <- matrix(w, length(t), length(u))
結果を確認する
> W [,1] [,2] [,3] [,4] [,5] [1,] 0.00862561 0.01590149 0.02188057 0.02674724 0.03066104 [2,] 0.01590149 0.02982485 0.04139268 0.05081327 0.05839446 [3,] 0.02188057 0.04139268 0.05813003 0.07188306 0.08296295 [4,] 0.02674724 0.05081327 0.07188306 0.08970176 0.10417950 [5,] 0.03066104 0.05839446 0.08296295 0.10417950 0.12189849
続けて、キャッシュ・フローを書く
C <- rbind(c(1.01, 0, 0, 0, 0), c(0.02, 1.02, 0, 0, 0), c(0.026, 0.026, 1.026, 0, 0), c(0.034, 0.034, 0.034, 0.034, 1.034))
順々に計算する
> C %*% W %*% t(C) [,1] [,2] [,3] [,4] [1,] 0.008798985 0.01655595 0.02331804 0.03453269 [2,] 0.016555948 0.03168201 0.04499267 0.06705478 [3,] 0.023318045 0.04499267 0.06461534 0.09733744 [4,] 0.034532685 0.06705478 0.09733744 0.15049244 > solve(C %*% W %*% t(C)) [,1] [,2] [,3] [,4] [1,] 10658.6389 -10190.359 3652.987 -267.9968 [2,] -10190.3594 14337.601 -7987.511 1116.2003 [3,] 3652.9867 -7987.511 6252.300 -1323.1861 [4,] -267.9968 1116.200 -1323.186 426.6236 > mu <- exp(-ufr * 1:5) > 1 - C %*% mu [,1] [1,] 0.03071017 [2,] 0.04137547 [3,] 0.04423345 [4,] 0.03541536 > xi <- solve(C %*% W %*% t(C)) %*% (1 - C %*% mu) > xi [,1] [1,] 57.790688 [2,] -33.507208 [3,] 11.396473 [4,] -5.466968
パッケージを使って計算したxiの値と同じになることが確認できた
さらに、その下の部分についても計算してみる
> t <- 4 > w4 <- sapply(1:5, function(i) fWilson(t, i, ufr, alpha)) > w4 [1] 0.02674724 0.05081327 0.07188306 0.08970176 0.10417950 > K4 <- t(w4) %*% t(C) > K4 [,1] [,2] [,3] [,4] [1,] 0.02701471 0.05236448 0.07576859 0.1158525 > KF <- K4 %*% xi > KF [,1] [1,] 0.03674387 > exp(-ufr * 4) [1] 0.8482603 > P <- exp(-ufr * 4) + KF > P [,1] [1,] 0.8850041 > P ^ (-1/4) - 1 [,1] [1,] 0.03101189
一致することが確認できた
Example2もキャッシュ・フローの頻度が4回になるように書くだけで確認できる
> alpha <- 0.1 > ufr <- log(1 + 0.042) > Type <- rep("SWAP", 4) > Tenor <- c(1, 2, 3, 5) > Rate <- c(0.01, 0.02, 0.026, 0.034) > Frequency <- rep(4, 4) > dfInstruments <- data.frame(Type, Tenor, Rate, Frequency) > Curve <- fFitSmithWilsonYieldCurveToInstruments(dfInstruments, ufr, alpha) > Curve$xi [,1] [1,] 58.629220 [2,] -34.081520 [3,] 11.818684 [4,] -5.744844 > plot(Curve)
ycinterextraというパッケージではNelson-SiegelやSvenssonと同様にSmith-Wilsonも指定すれば使えるようなので、ycinterextraパッケージの調査をまたの機会に
Rの金利期間構造パッケージ SmithWilsonYieldCurve
パッケージを使ってSmith-Wilsonメソッドをやってみたメモ
パッケージはSmithWilsonYieldCurve。一応、作成者による説明があります
しかし、TeX記法でそのままズラズラ書いてあるし、何も読まずプログラムをコピペするだけだと動かず…ということで、これをひと通り再現してみる
マニュアルには関数の説明が詳しく書いてあるが、fFitSmithWilsonYieldCurveとfFitSmithWilsonYieldCurveToInstrumentsしかpublicとして使えないので、Method Descriptionの部分のプロットをするためにはWilson関数を書く必要がある
fWilson <- function(t, u, ufr, alpha) { mintu <- min(t, u) maxtu <- max(t, u) exp(-ufr * (t + u)) * (alpha * mintu - 0.5 * exp(-alpha * maxtu) * (exp(alpha * mintu) - exp(-alpha * mintu))) }
ggplot2でプロットする
library(ggplot2) t <- 0:100 u <- c(1, 5, 10, 25, 50) f <- 0.042 alpha <- 0.01 g <- expand.grid(t, u) w <- sapply(1:nrow(g), function(i) fWilson(g$Var1[i], g$Var2[i], ufr = f, alpha = alpha)) df <- data.frame(g$Var1, as.factor(g$Var2), w) colnames(df) <- c("t", "u", "value") ggplot(df, aes(x = t, y = value, group = u)) + geom_line(aes(colour = u), size = 1) + labs(colour = "Value of u")
次のSimple Exampleでは、基本の割引関数のレートが4.2％で、価格が0.88の5年物と0.37の20年物のゼロクーポン債、それぞれの利回りが2.56％と4.97％に相当するものが存在すると仮定している
library(SmithWilsonYieldCurve) m1 <- 0.88 m2 <- 0.37 y1 <- 0.0256 y2 <- 0.0497 u1 <- 5 u2 <- 20 f <- 0.042 alpha <- 0.1 C <- diag(2) m <- c(m1, m2) u <- c(u1, u2) Curve <- fFitSmithWilsonYieldCurve(TimesVector = u, CashflowMatrix = C, MarketValueVector = m, ufr = f, alpha = alpha)
結果を確認する
> Curve $P function (t) { fBase(t) + t(KernelWeights) %*% fCompoundKernel(t) } <environment: 0x0000000010ccc728> $xi [,1] [1,] 2.524894 [2,] -1.568533 $K function (t) { CashflowMatrix %*% fKernel(t, TimesVector) } <environment: 0x0000000010ccc728> attr(,"class") [1] "SmithWilsonYieldCurve" "YieldCurve"
カーネルに適用される重みが2.524、-1.569であるとわかる
この状況をプロットしたものは次のようにして再現可能
Term <- 1:100 Rate <- -log(Curve$P(Term)) / Term plot(Term, Rate, ylim = c(0, 0.055)) abline(h = c(f, y1, y2), col = c(2, 4, 5)) legend("bottomright", c("Curve", expression(r[1]), expression(r[2]), "f"), fill = c(1, 4, 5, 2))
その次は、そのままコピペで大丈夫
par(bg = "white") plot(1:100, -log(1 + t(Curve$xi) %*% Curve$K(1:100)/exp(-f * (1:100)))/1:100, xlab = "Term", ylab = "Rate", col = "green")
この状況では、期間が短い部分では負の方向に、中期では正の方向に効いていて、その後、長期のフォワードレート4.2%の方向に向かってゆるやかにゼロに近づいていく
各時間におけるゼロクーポン債の価格もプロットしてみる
plot(1:100, Curve$P(1:100), xlab = "Term", ylab = "Price")
というところで、ひと通りR-bloggersに書かれていることは再現できた。上記の例がシンプルな分、分かったような分からないような…もう少し細かい部分は、もう一つの関数fFitSmithWilsonYieldCurveToInstrumentsを使いながらまたの機会に
モンテカルロシミュレーションでVasicekモデルの債券価格
Vasicekモデルでパスの発生方法を比較の続き、モンテカルロシミュレーションから債券価格を求めてみる
パラメータなどは前回と同じで、それぞれ1000回シミュレーション
nsim <- 1000 p1 <- replicate(nsim, Vasicek_euler(kappa, mu, sigma, r0, dt, n)) p2 <- replicate(nsim, Vasicek_exact(kappa, mu, sigma, r0, dt, n)) p3 <- replicate(nsim, sde.sim(X0 = r0, N = n, delta = dt, theta = c(kappa * mu, kappa, sigma), model = "OU")[-(n+1)]) p4 <- replicate(nsim, sde.sim(X0 = r0, N = n, delta = dt, method = "cdist", rcdist = rcOU, theta = c(kappa * mu, kappa, sigma))[-(n+1)]) p5 <- replicate(nsim, sde.sim(X0 = r0, N = n, delta = dt, drift = expression(0.2 * (0.015 - x)), sigma = expression(0.003))[-(n+1)]) p6 <- replicate(nsim, simulate(model, true.param = list(kappa = kappa, mu = mu, sigma = sigma), sampling = sampling)@data@original.data[-(n+1)])
replicateはシミュレーションしたパスをmatrixの形で返してくる。ここから、債券価格を求める関数を作る
MCZCBond <- function(r, dt) { rsum <- colSums(r) sbond <- exp(-rsum * dt) mean(sbond) }
それぞれの方法について債券価格を計算する
> MCZCBond(p1, dt) [1] 0.9162368 > MCZCBond(p2, dt) [1] 0.9140108 > MCZCBond(p3, dt) [1] 0.9145012 > MCZCBond(p4, dt) [1] 0.9144614 > MCZCBond(p5, dt) [1] 0.9147488 > MCZCBond(p6, dt) [1] 0.9155063
明示的に求まるVasicekモデルの債券価格も計算する
> t <- 0 > T <- dt * n > VasicekZCBond(kappa, mu, sigma, r0, t, T)$price [1] 0.9148043
今回はsdeパッケージを使ったドリフト項とディフュージョン項を指定する方法が最も差が小さいという結果に
何回かシミュレーションをやってみると、必ずどの方法が明示的に求まる債券価格に近いということはなく、パスの発生方法の違いによって生じる誤差よりもシミュレーションした時に生じる誤差の方が大きいと言えそうです。モンテカルロシミュレーションに限らず、Vasicekモデルは基本を確認するには大変役立ちます
Vasicekモデルでパスの発生方法を比較
Vasicekモデルの債券価格と債券オプションについて書いたので、次はパスの発生方法を比較してみる
はじめはEuler–Maruyama method - Wikipediaに書かれているPythonコードをそのままRに書き直したもの。愚直にforループを回してみる
Vasicek_euler <- function(kappa, mu, sigma, r0, dt, n) { r <- numeric(n) r[1] <- r0 w <- rnorm(n) for (i in 2:n) { r[i] <- r[i-1] + kappa * (mu - r[i-1]) * dt + sigma * sqrt(dt) * w[i-1] } r }
週次データで10年までを想定してパラメータをセット
kappa <- 0.2 mu <- 0.015 sigma <- 0.003 r0 <- 0.001 dt <- 1/52 n <- 520
パスを発生させる
> set.seed(123) > p1 <- Vasicek_euler(kappa, mu, sigma, r0, dt, n) > head(p1) [1] 0.0010000000 0.0008206742 0.0007794504 0.0014826067 0.0015639300 0.0016693941
複数パスを発生させて、プロット
library(dplyr) library(tidyr) library(ggplot2) time <- 1:n * dt paths <- replicate(30, Vasicek_euler(kappa, mu, sigma, r0, dt, n)) cbind(time, paths) %>% as.data.frame() %>% gather(path, rate, -time) %>% ggplot(aes(x = time, y = rate, group = path)) + geom_line(aes(colour = path)) + ggtitle("Short Rate Paths of Vasicek Model") + theme(legend.position = "none")
時間経過とともに長期の平均回帰水準0.015付近で推移していくようになる。今ではものすごくおかしなことではないのかもしれないが、Vasicekモデルでは金利がマイナスの状態になることもある
VasicekモデルはOU過程であり、その平均と分散も明確にわかるので、より厳密に書くことができる
Vasicek_exact <- function(kappa, mu, sigma, r0, dt, n) { r <- numeric(n) r[1] <- r0 w <- rnorm(n) e <- exp(-kappa * dt) s <- sigma * sqrt((1 - e^2) / (2 * kappa)) for (i in 2:n) { r[i] <- e * r[i-1] + mu * (1 - e) + s * w[i-1] } r }
こちらも同じパラメータでパスを発生させて、前者と比較する
> set.seed(123) > p2 <- Vasicek_exact(kappa, mu, sigma, r0, dt, n) > head(p2) [1] 0.0010000000 0.0008210185 0.0007798725 0.0014816771 0.0015628477 0.0016681135 > head(p1 - p2) [1] 0.000000e+00 -3.442722e-07 -4.220666e-07 9.296463e-07 1.082250e-06 1.280618e-06
多少ではあるが、誤差が出ることがわかる
パッケージを使ったものとも比較してみる。sdeパッケージのsde.simは異なる方法でパスを発生させることできる。初期値からnステップ発生するので、n+1番目は抜いておく
library(sde) set.seed(123) p3 <- sde.sim(X0 = r0, N = n, delta = dt, theta = c(kappa * mu, kappa, sigma), model = "OU")[-(n+1)] set.seed(123) p4 <- sde.sim(X0 = r0, N = n, delta = dt, method = "cdist", rcdist = rcOU, theta = c(kappa * mu, kappa, sigma))[-(n+1)] set.seed(123) p5 <- sde.sim(X0 = r0, N = n, delta = dt, drift = expression(0.2 * (0.015 - x)), sigma = expression(0.003))[-(n+1)]
それぞれ、モデルを直接指定する方法、条件付き分布から発生させる方法、ドリフト項とディフュージョン項を指定する方法
こちらも一部を確認
> head(p3) [1] 0.0010000000 0.0008210185 0.0007798725 0.0014816771 0.0015628477 0.0016681135 > head(p4) [1] 0.0010000000 0.0008210185 0.0007798725 0.0014816771 0.0015628477 0.0016681135 > head(p5) [1] 0.0010000000 0.0008210190 0.0007798732 0.0014816757 0.0015628461 0.0016681116
当然だが、自作関数のものとほとんど同じ値になることを確認できる
yuimaパッケージはさらに複雑なものもシミュレーションすることができるようになっている。ここでは、一番簡単そうな方法を一つだけやってみる
library(yuima) set.seed(123) sampling <- setSampling(n = n, delta = dt) model <- setModel(drift = "kappa * (mu - x)", diffusion = "sigma", state.var = "x", time.var = "t", solve.var = "x", xinit = r0) p6 <- simulate(model, true.param = list(kappa = kappa, mu = mu, sigma = sigma), sampling = sampling)@data@original.data[-(n+1)]
一部を確認
> head(p6) [1] 0.0010000000 0.0008206742 0.0007794504 0.0014826067 0.0015639300 0.0016693941
結局、どれくらいの違いがあるのか？
> sum((p1 - p2)^2) [1] 5.358932e-08 > sum((p1 - p3)^2) [1] 5.358932e-08 > sum((p1 - p4)^2) [1] 5.358932e-08 > sum((p1 - p5)^2) [1] 5.373786e-08 > sum((p1 - p6)^2) [1] 1.982772e-32 > sum((p2 - p3)^2) [1] 3.662136e-30 > sum((p2 - p4)^2) [1] 3.662136e-30 > sum((p2 - p5)^2) [1] 1.03119e-13 > sum((p2 - p6)^2) [1] 5.358932e-08 > sum((p3 - p4)^2) [1] 0 > sum((p3 - p5)^2) [1] 1.03119e-13 > sum((p3 - p6)^2) [1] 5.358932e-08 > sum((p4 - p5)^2) [1] 1.03119e-13 > sum((p4 - p6)^2) [1] 5.358932e-08 > sum((p5 - p6)^2) [1] 5.373786e-08
多少の違いはあるようですが、sdeパッケージを使った3番目と4番目は全く同じ。2番目のものとsdeパッケージを使ったものとはほとんど同じ。1番はじめのEuler法のものとyuimaパッケージを使ったものはほとんど同じという結果
次は、モンテカルロシミュレーションでそれぞれのパスから債券価格を求めてみます
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モデルは実務に耐え得るモデルではないものの、いくつかの短期金利モデルを比較するような論文には登場することは多いので、これぐらいはパッケージ化されていても良い気もしますが…
