Rの金利期間構造パッケージ ycinterextra
金利の期間構造関連パッケージを使ってみるシリーズ。今回はycinterextraを調査。このパッケージもこれといった情報がないので、普通にパッケージの説明書を読んで、日本の金利データに適用してみる
パッケージの中身は、CRAN Task View: Empirical Financeにも書かれている通り、Nelson-Siegelモデル等の様々なモデルでイールドカーブの補間と補外を行うことができるパッケージだと言うこと
主な関数はycinterとycextraで、その他ユーティリティ的な関数が用意されている
ycinter(yM = NULL, p = NULL, matsin, matsout, method = c("NS", "SV", "SW", "HCSPL"), typeres = c("rates", "prices")) ycextra(yM = NULL, p = NULL, matsin, matsout, method = c("NS", "SV", "SW"), typeres = c("rates", "prices"), UFR, T_UFR = NULL)
ycinterとycextraともに、ゼロクーポン債価格とゼロクーポンイールドを対象としているので、クーポン債に適用したい場合にはその前に処理しておく必要がある
jgbcme <- read.csv("http://www.mof.go.jp/english/jgbs/reference/interest_rate/jgbcme.csv", skip = 1) JGB_rate <- as.numeric(jgbcme[1, -1]) / 100 maturity <- c(1:10, 15, 20, 25, 30, 40)
Nelson-SiegelモデルとSvenssonモデルをデータに適用していく
library(ycinterextra) NSyc <- ycinter(yM = JGB_rate, matsin = maturity, matsout = 1:40, method = "NS", typeres = "rates") SVyc <- ycinter(yM = JGB_rate, matsin = maturity, matsout = 1:40, method = "SV", typeres = "rates")
求まったパラメータの値は次の通り
> coeffs(NSyc) [1] 0.01847065 -0.01415590 -0.03796077 3.00000000 > coeffs(SVyc) [1] 0.02019232 -0.02168073 0.04297268 -0.07961120 2.01990134 3.00000000
3.00000000と出ていることからも上手く最適化できていない可能性が高そうです…
ycplotを使って、当てはまりの診断を見てみる
> ycplot(NSyc) > ycplot(SVyc)
Nelson-Siegelモデル
Svenssonモデル
残差が本当に正規なのというところは疑問が残るところですが、それなりに補間は上手く行っているようです
（method = "HCSPL"を指定すれば、Hermite cubic splineによる補間もできるのだが、上と同じようにしてもエラーが出てなぜか上手く適用できなかった…）
同じように補外の方も1年から60年までの設定で試してみる
ここで、UFRの値は上のパラメータの値から0.02を設定している
NSyc <- ycextra(yM = JGB_rate, matsin = maturity, matsout = 1:60, method = "NS", typeres = "rates", UFR = 0.02) SVyc <- ycextra(yM = JGB_rate, matsin = maturity, matsout = 1:60, method = "SV", typeres = "rates", UFR = 0.02)
求まったパラメータの値は次の通り
> coeffs(NSyc) [1] 0.02000000 -0.01457810 -0.04355021 3.00000000 > coeffs(SVyc) [1] 0.02000000 -0.02162734 0.03623047 -0.07322875 1.90287193 3.00000000
やはり、上手く最適化できていないようです…
こちらもycplotで確認
> ycplot(NSyc) > ycplot(SVyc)
Nelson-Siegelモデル
Svenssonモデル
少々、簡単でしたがycinterextraパッケージの関数を使って、実際データに適用してみました。今の日本の金利データがこれらのモデルに適用するというのが難しい面もありそうですが、パラメータを求めるあたりはあまり上手くいっていないようです。このパッケージでは、これらのモデル以外にも、Smith-Wilson methodも使えるのでお手軽に試してみたいという時は良さそうです
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パッケージを使ったものはほとんど同じという結果
次は、モンテカルロシミュレーションでそれぞれのパスから債券価格を求めてみます
