CIRモデルのパラメータ推定
CIRモデルについてもVasicekモデルと同様にシンプルに最尤推定でモデルのパラメータを推定できる関数を用意しておく
# maximum likelihood estimation for CIR model # # rate : interest rate # param : model parameter # delta : time step # initparam : initial parameter # CIRloglik <- function(rate, param, delta) { n <- length(rate) c <- 4 * param[1] / (param[3] ^ 2 * (1 - exp(-param[1] * delta))) df <- 4 * param[1] * param[2] / (param[3] ^ 2) ncp <- c * rate[1:(n - 1)] * exp(-param[1] * delta) sum(dchisq(c * rate[2:n], df = df, ncp = ncp, log = TRUE) + log(c)) } CIRMLE <- function(rate, initparam, delta) { opt <- optim(initparam, CIRloglik, rate = rate, delta = delta, control = list(fnscale = -1)) c(kappa = opt$par[1], mu = opt$par[2], sigma = opt$par[3], loglik = opt$value, convergence = opt$convergence) }
パラメータを設定して、パスを発生させる
library(sde) r0 <- 0.01 n <- 500 delta <- 1/52 kappa <- 0.2 mu <- 0.03 sigma <- 0.08 initparam <- c(kappa, mu, sigma) set.seed(12345) rate <- sde.sim(X0 = r0, N = n, delta = delta, theta = c(kappa * mu, kappa, sigma), model = "CIR")
推定結果を確認する
> CIRMLE(rate, initparam, delta) kappa mu sigma loglik convergence 2.789964e-01 6.105776e-02 8.106622e-02 2.376850e+03 0.000000e+00
初期値の与え方にもよるが、Vasicekモデルと比べて多少パラメータの推定値にバイアスが出ることになります。 特に平均回帰パラメータkappaについては上方に推定する傾向にあります
Rからアメリカ雇用統計データ収集
FXをやっている人なら誰しも気になる経済統計といえば、アメリカの雇用統計。 【米国】雇用統計 - 経済指標 - Yahoo!ファイナンスでも直近の数年間のデータを詳しく見ることができますが、簡単にRから直接データを収集して、分析したいという人はblsAPIパッケージを使ってみるのが良さそうです
こちらmikeasilva/blsAPI · GitHubにもサンプルがありますが、発表時に一番注目される非農業部門雇用者数と失業率のデータを使った例がないようなので、そのデータを取得してから可視化するところまでをメモ
使い方は簡単で、パッケージに唯一含まれているblsAPIという関数に欲しいデータのIDコードを与えてやればデータを取って来れます。その他、list形式でパラメータを指定してやれば、いろいろと自由にデータ収集可能です
library(rjson) library(blsAPI) library(dplyr) payload <- list( 'seriesid' = c('CES0000000001', 'LNS14000000'), 'startyear' = 2007, 'endyear' = 2015) response <- blsAPI(payload) %>% fromJSON()
これだけで、2007年から2015年までの雇用者数と失業率のデータが取って来れます。 responseの中に2系列入っており、年月についても整形してdata.frameに変換してやります
series1 <- response$Results$series[[1]]$data series2 <- response$Results$series[[2]]$data date_str <- 0 value1 <- value2 <- 0 for (i in 1:length(series1)) { date_str[i] <- paste(series1[[i]]$year, sub('M', '', series1[[i]]$period), "01", sep = "-") value1[i] <- series1[[i]]$value value2[i] <- series2[[i]]$value } date <- as.Date(date_str[order(date_str)])[-1] data1 <- diff(as.numeric(value1[order(date_str)])) data2 <- as.numeric(value2[order(date_str)])[-1] df <- data.frame(date, data1, data2)
データを確認します
> head(df) date data1 data2 1 2007-02-01 88 4.5 2 2007-03-01 188 4.4 3 2007-04-01 78 4.5 4 2007-05-01 145 4.4 5 2007-06-01 71 4.6 6 2007-07-01 -34 4.7 > tail(df) date data1 data2 95 2014-12-01 329 5.6 96 2015-01-01 201 5.7 97 2015-02-01 266 5.5 98 2015-03-01 119 5.5 99 2015-04-01 221 5.4 100 2015-05-01 280 5.5
前月比で雇用者数の変化を棒グラフにする
library(ggplot2) library(scales) df %>% ggplot(aes(x = date)) + geom_bar(aes(y = data1), stat = "identity", fill = "blue") + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + scale_x_date(labels = date_format("%m/%y"), breaks = date_breaks("6 months")) + scale_y_continuous(breaks = seq(-1000, 1000, by = 200)) + ylab("1-Month Net Change")
失業率もプロットする
df %>% ggplot(aes(x = date)) + geom_line(aes(y = data2), size = 1, colour = "red") + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + scale_x_date(labels = date_format("%m/%y"), breaks = date_breaks("6 months")) + scale_y_continuous(breaks = seq(0, 12, by = 1)) + ylab("Unemployment Rate")
Vasicekモデルのパラメータ推定
SMFI5パッケージでもVasicekモデルのパラメータ推定はできるが、参考程度にパラメータを推定する際、シンプルに最尤推定できるものを用意しておいて、簡単に使い回せるようにしておきたい
# maximum likelihood estimation for Vasicek model # # rate : interest rate # param : Vasicek model parameter # delta : time step # initparam : initial parameter # Vasicekloglik <- function(rate, param, delta) { n <- length(rate) m <- param[2] + (rate[1:(n-1)] - param[2]) * exp(-param[1] * delta) v <- param[3]^2 * (1 - exp(-2 * param[1] * delta)) / (2 * param[1]) sum(dnorm(rate[2:n], mean = m, sd = sqrt(v), log = TRUE)) } VasicekMLE <- function(rate, initparam, delta) { opt <- optim(initparam, Vasicekloglik, rate = rate, delta = delta, control = list(fnscale = -1)) c(kappa = opt$par[1], mu = opt$par[2], sigma = opt$par[3], loglik = opt$value, convergence = opt$convergence) }
optimのhessian = TRUEを足したり、２つの関数を１つにしてしまったりということもできる
sdeパッケージを使ってパスを発生させ、パラメータ推定する
library(sde) r0 <- 0.01 n <- 500 delta <- 1/52 kappa <- 0.2 mu <- 0.03 sigma <- 0.01 initparam <- c(kappa, mu, sigma) set.seed(123) rate <- sde.sim(X0 = r0, N = n, delta = delta, theta = c(kappa * mu, kappa, sigma), model = "OU")
結果を確認する
> VasicekMLE(rate, initparam, delta) kappa mu sigma loglik convergence 5.670395e-01 3.114274e-02 9.740629e-03 2.596687e+03 0.000000e+00
最適化の初期値に真の値を入れていることもあり、良い値に推定できています
金利の平均回帰速度kappaはだいぶ真の値よりも上方に推定されていますが、推定バイアスを改善させましょうというような論文がいくつかありますので、この辺をじっくり読むというのも大変面白いです
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パッケージの調査をまたの機会に
