xtsでmatplotのようなプロット
最近はdplyrやtidyrとの相性で、ほぼggplot2を使うようになり、めっきりmatplotを使う機会が無くなってきています。
しかし、時系列データを扱う際にdata.frameよりもxtsを使いたい場合もあり、その時にxtsのままでmatplotのような感じのプロットするにはどうするのがいいのか？と思っていたので調べたメモします。結局、xtsとxtsExtraをセットで使うのが良さそうということがわかりました。
まず、データとして、複数時系列データをQuandlから取ってきます。この時に、Quandl関数でtypeにxtsを指定すれば、xtsクラスでデータを取って来れますが、一旦、ggplot2でプロットするためにdata.frameで取って来ることにします。
library(Quandl) rate_df <- Quandl("MOFJ/INTEREST_RATE_JAPAN", trim_start = "2015-01-01")
このデータをggplot2でプロットします。
library(dplyr) library(tidyr) library(ggplot2) library(scales) rate_df %>% gather(maturity, rate, -Date) %>% ggplot(aes(x = Date, y = rate, colour = maturity)) + geom_line(size = 1) + scale_x_date(labels = date_format("%m/%d"))
その他、色々と付け足していけば、思いどおりのプロットができるでしょう。しかし、本題はデータをdata.frameではなくxtsクラスで扱いたい場合です。まわりくどいですが、Quandlの中でやっている処理と同じようにdata.frameからxtsに変換します。
rate_xts <- xts(rate_df[, -1], order.by = rate_df[, 1])
ここで、そのままplotすれば複数系列を一辺にプロットしてくれると思いきや、警告メッセージとともに残念なプロットがされます。
> plot(rate_xts) Warning message: In plot.xts(rate_xts) : only the univariate series will be plotted
matplotならばどうかというと、xtsプロットの良い部分は完全に失われます。
> matplot(rate_xts, type = "l")
xtsはzooクラスとしても扱えるのだから、plot.zooで良い感じにプロットされるのではないかと思い、やってみるとmatplotとほぼ変わりません…
> plot.zoo(rate_xts, plot.type = "single")
というところまで試してみて、CRANには上がっていないxtsExtraを使うのが一番良さそうということが分かり、使ってみるとxtsプロットの良い部分もあり、matplotのような複数系列のプロットが簡単にできました。
library(devtools) devtools::install_github("joshuaulrich/xtsExtra") library(xtsExtra) plot(rate_xts, col = 1:15)
ただ、デフォルトでplotするだけでは12系列までしかプロットされないようなので、色の指定を系列数分に増やす、または、addSeriesを使って追加する必要があるようです。matplotしたものの軸ラベル等をxtsプロット風に頑張って変更するよりも、断然xtsExtraを使うのが良さそうです。
pforeachで最尤推定した漸近分布を確認する
最尤推定で推定したパラメータの漸近分布をシミュレーションによって確認する
繰り返し推定を行う際にpforeach(hoxo-m/pforeach · GitHub)を使う
CIRモデルのパラメータ推定と関数もパラメータの設定も同じとする
r0 <- 0.01 n <- 1000 delta <- 1/52 kappa <- 0.2 mu <- 0.03 sigma <- 0.08 initparam <- c(kappa, mu, sigma)
パスを発生させ、パラメータ推定を繰り返す。収束しなかったものは取り除く
ggplot2で推定したパラメータの分布をプロットしたいので、最終的にdata.frameで返ってくるようにする
.c = data.frameとして、data.frameのまま結合してもパラメータの名前がついたdata.frameが返って来なかったため、rbindで結果を結合したあとにdata.frameに変換している（これが最適な方法なのかについては…）
library(sde) library(dplyr) library(pforeach) est <- pforeach(i = 1:1000, .c = rbind)({ sde.sim(X0 = r0, N = n, delta = delta, theta = c(kappa * mu, kappa, sigma), model = "CIR") %>% CIRMLE(initparam, delta) }) %>% as.data.frame %>% filter(convergence != 1)
foreachの少し込み入った使い方に面倒くささを感じていたが、pforeachを使えば簡単、簡潔に記述できる
replicateを使う場合には、データを発生させて、パラメータを推定してという一連の処理を書いた関数を書いておいて、その関数を繰り返すということをしなければならないが、pforeachではループの中にザッと処理を書いておくだけで繰り返し推定が行える。返り値もあるので、その後一度に加工して扱える。素晴らしい
推定したパラメータの分布をプロットする
library(ggplot2) ggplot(est, aes(x = kappa)) + geom_histogram(binwidth = 0.1) ggplot(est, aes(x = mu)) + geom_histogram(binwidth = 0.005) ggplot(est, aes(x = sigma)) + geom_histogram(binwidth = 0.0001)
平均回帰パラメータkappaの分布のみここでは確認する。推定された値の平均を計算すると真値よりも大きい値となる。そして、プロットしても上方に推定するバイアスが生じることが確認できる（TRUE = 0.2）
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付近でモデルの違いがオプション価格に現れてくることがわかります。初期金利の設定の仕方によっては変わってくる部分もありますが、細かいところを論ずる部分は教科書等を参考にお願いします…
