CIRモデルの債券オプション
Vasicekモデルと同じ流れでCIRの債券オプションについて
モデルの一部が違うだけでも、CIRモデルとVasicekモデルでこのような違いが出てくるのかといった部分を追っていくのも面白い。ただ、プログラミングするのは少し面倒となる
# zero coupon bond option price under Cox-Ingersoll-Ross 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 # CIRZCBOption <- function(kappa, mu, sigma, r0, t, T, S, K, L) { tau <- S - T g <- sqrt(kappa^2 + 2 * sigma^2) A <- ((2 * g * exp((kappa + g) * tau / 2)) / (2 * g + (kappa + g) * (exp(tau * g) - 1)))^(2 * kappa * mu / sigma^2) B <- 2 * (exp(tau * g) - 1) / (2 * g + (kappa + g) * (exp(tau * g) - 1)) ra <- log(L * A / K) / B l1 <- 2 * g / (sigma^2 * (exp((T - t) * g) - 1)) l2 <- (kappa + g) / sigma^2 d <- 4 * kappa * mu / sigma^2 ncp1 <- 2 * l1^2 * r0 * exp((T - t) * g) / (l1 + l2 + B) ncp2 <- 2 * l1^2 * r0 * exp((T - t) * g) / (l1 + l2) q1 <- 2 * ra * (l1 + l2 + B) q2 <- 2 * ra * (l1 + l2) PT <- CIRZCBond(kappa, mu, sigma, r0, t, T)$price PS <- CIRZCBond(kappa, mu, sigma, r0, t, S)$price M <- K / (L * PS / PT) CP <- L * PS * pchisq(q1, d, ncp1) - K * PT * pchisq(q2, d, ncp2) PP <- CP + K * PT - L * PS data.frame(moneyness = M, callprice = CP, putprice = PP) }
パラメータをセットする
kappa <- 0.2 mu <- 0.015 sigma <- 0.05 r0 <- 0.001 t <- 0 T <- 1/4 S <- 3 K <- 0.98 L <- 1
オプション価格を確認する
> CIRZCBOption(kappa, mu, sigma, r0, t, T, S, K, L) moneyness callprice putprice 1 0.9928814 0.00702689 2.968141e-06
パラメータを調整して、CIRモデルとVasicekモデルを完全に同じ条件で比較することは難しいが、CIRモデルの債券価格の場合と同様にボラティリティの部分を調整することで、確認することができる
オプション満期で債券価格がほぼ同じに評価されるようにパラメータを調整して、債券価格と債券オプション価格を計算してみる
> sigma_v <- sigma * sqrt(mu) > VasicekZCBond(kappa, mu, sigma_v, r0, t, T) maturity price yield 1 0.25 0.9996641 0.001343863 > CIRZCBond(kappa, mu, sigma, r0, t, T) maturity price yield 1 0.25 0.999664 0.001344209 > VasicekZCBOption(kappa, mu, sigma_v, r0, t, T, S, K, L) moneyness callprice putprice 1 0.992794 0.007500676 0.0003899409 > CIRZCBOption(kappa, mu, sigma, r0, t, T, S, K, L) moneyness callprice putprice 1 0.9928814 0.00702689 2.968141e-06
基本的に、CIRモデルとVasicekモデルの債券オプション価格を比較した場合には、Vasicekモデルの短期金利がマイナスになるという可能性がある点で違いが出てくる。CIRモデルは非心カイ自乗分布、Vasicekモデルは正規分布に短期金利が従うことから、CIRモデルはVasicekモデルに比べて、債券価格を低く評価する傾向がある。そこから、債券オプション価格も権利行使価格が額面価格に近づくにつれて、CIRモデルはVasicekモデルよりも低く評価する傾向にある
権利行使価格を変化させながら、オプション価格を計算する
K <- seq(0.97, 1.01, by = 0.001) df1 <- VasicekZCBOption(kappa, mu, sigma_v, r0, t, T, S, K, L) df2 <- CIRZCBOption(kappa, mu, sigma, r0, t, T, S, K, L) df <- cbind(rbind(df1, df2), model = c(rep("Vasicek", nrow(df1)), rep("CIR", nrow(df2))))
ggplot2でプロットする
library(ggplot2) library(ggthemes) ggplot(df, aes(x = moneyness, y = callprice, colour = model)) + geom_line(size = 1.2) + theme_economist() + ggtitle("call option price") + scale_x_continuous(breaks = seq(0.98, 1.02, by = 0.01)) ggplot(df, aes(x = moneyness, y = putprice, colour = model)) + geom_line(size = 1.2) + theme_economist() + ggtitle("put option price") + scale_x_continuous(breaks = seq(0.98, 1.02, by = 0.01))
ATM付近でモデルの違いがオプション価格に現れてくることがわかります。初期金利の設定の仕方によっては変わってくる部分もありますが、細かいところを論ずる部分は教科書等を参考にお願いします…
CIRモデルの債券価格とイールド
Vasicekモデルの方はやったので、CIRモデルの債券価格とイールドについて
CIRモデルの債券価格もVasicekモデルの場合と同様にアフィン型となり、解析的に求めることができる
# zero coupon bond price and yield under Cox-Ingersoll-Ross 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 # CIRZCBond <- function(kappa, mu, sigma, r0, t, T) { tau <- T - t h <- sqrt(kappa^2 + 2 * sigma^2) B <- 2 * (exp(tau * h) - 1) / (2 * h + (kappa + h) * (exp(tau * h) - 1)) A <- ((2 * h * exp((kappa + h) * tau / 2)) / (2 * h + (kappa + h) * (exp(tau * h) - 1)))^(2 * kappa * mu / sigma^2) P <- A * exp(-B * r0) Y <- ifelse(tau != 0, -log(P) / tau, r0) data.frame(maturity = tau, price = P, yield = Y) }
パラメータをセットする
kappa <- 0.2 mu <- 0.015 sigma <- 0.02 r0 <- 0.001 t <- 0 T <- c(1/12, 1/4, 1/2, 1:10, 15, 20, 30)
債券価格とイールドは次のようになる
> CIRZCBond(kappa, mu, sigma, r0, t, T) maturity price yield 1 0.08333333 0.9999070 0.001116021 2 0.25000000 0.9996640 0.001344234 3 0.50000000 0.9991617 0.001677218 4 1.00000000 0.9976916 0.002311055 5 2.00000000 0.9931024 0.003460729 6 3.00000000 0.9866763 0.004471087 7 4.00000000 0.9787842 0.005361026 8 5.00000000 0.9697338 0.006146732 9 6.00000000 0.9597787 0.006842082 10 7.00000000 0.9491268 0.007458982 11 8.00000000 0.9379476 0.008007653 12 9.00000000 0.9263789 0.008496884 13 10.00000000 0.9145324 0.008934236 14 15.00000000 0.8537373 0.010542117 15 20.00000000 0.7940308 0.011531650 16 30.00000000 0.6846598 0.012627772
library(ggplot2) library(ggthemes) r0 <- c(0, 0.01, 0.015, 0.02, 0.03) tl <- length(T) irate <- c(rep(r0[1], tl), rep(r0[2], tl), rep(r0[3], tl), rep(r0[4], tl), rep(r0[5], tl)) y1 <- CIRZCBond(kappa, mu, sigma, r0[1], t, T) y2 <- CIRZCBond(kappa, mu, sigma, r0[2], t, T) y3 <- CIRZCBond(kappa, mu, sigma, r0[3], t, T) y4 <- CIRZCBond(kappa, mu, sigma, r0[4], t, T) y5 <- CIRZCBond(kappa, mu, sigma, r0[5], t, T) df <- data.frame(irate = factor(irate), rbind(y1, y2, y3, y4, y5)) ggplot(data = df, aes(x = maturity, y = yield, colour = irate)) + geom_line(size = 1.2) + theme_economist() + scale_colour_economist()
VasicekモデルとCIRモデルの形を見ての通り、sigma CIR * sqrt(mu) = sigma Vasというようにすれば、ボラティリティの部分がある程度近似されるので、モデルの違いによる価格を比較する場合にはパラメータを調整してセットしてあげればよい
分布の違いによる影響はオプション価格の方が出やすいので、その比較を含めてCIRモデルの債券オプションは次の機会に
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年満期までの金利が完全に水没、水没しかかっているのがこちらでもはっきりわかる。このような状態で金融機関はちゃんと息できているんでしょうか？
