今週も特にありません

進捗どうですか?

Rの金利の期間構造まわりのパッケージ(3)

パッケージを使って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")

f:id:masaqol:20150303005923p:plain

次の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))

f:id:masaqol:20150303011230p:plain

その次は、そのままコピペで大丈夫

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")

f:id:masaqol:20150303011652p:plain

この状況では、期間が短い部分では負の方向に、中期では正の方向に効いていて、その後、長期のフォワードレート4.2%の方向に向かってゆるやかにゼロに近づいていく

各時間におけるゼロクーポン債の価格もプロットしてみる

plot(1:100, Curve$P(1:100), xlab = "Term", ylab = "Price")

f:id:masaqol:20150304005843p:plain

というところで、ひと通りR-bloggersに書かれていることは再現できた。上記の例がシンプルな分、分かったような分からないような…もう少し細かい部分は、もう一つの関数fFitSmithWilsonYieldCurveToInstrumentsを使いながらまたの機会に