Pricing defaultable discount bond with reduced form model

[This article was first published on My Life as a Mock Quant in English, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

I often use R language to write “prototype” program. As you know, It has very high productivity and smart grammar. In this article, I would like to show you how to write the program to evaluate the price of defaultable bond by “reduced-form model”.

Before write a program, we need to understand how to price these bond.
Under risk neutral measure, we can evaluate the price of defaultable (discount) bond as below.



In this equation, I set variables as below.

  • v(t,T) : The price of defaultable bond  with maturity T at time t 
  • r(t) : Short rate at time t
  • delta_(t) : Recovery rate at time t
  • tau : default time
  • 1{} : Indicator function
And the second term in this equation means usual pay back of principal, on the other hand, the first term means the value you can get when this bond will default.

For montecarlo simulation, we need to model short rate and default intensity process, I modeled default intensity process and short rate as CIR model.

CIR model assume that some stochastic process obey below Stochastic differential equation.

Each variables have following meaning.

  • x(t) : The value of stochastic variable obeying CIR process at time t.
  • kappa : The speed of adjustment
  • theta : The mean of  stochastic variable
  • sigma : The volatility of  stochastic variable
  • W_t : Standard Brownian motion

And next, I set the simulation condition as below.

——————————– Simulation condition ——————————-
  • The number of path : 500
  • The number of grid per year : 250
  • Recovery rate : 0.7
  • Maturity of bond : 2
  • The correlation between interest rate and default intensity : 0.3
——————————– Interest rate parameter —————————-
  • kappa : 0.6
  • theta : 0.05
  • sigma : 0.05
  • initial value : 0.05
——————————– Default intensity parameter ————————-
  • kappa : 0.6
  • theta : 0.2
  • sigma : 0.5
  • initial value : 0.2
———————————————————————————
After these settings, we can write simulation program. The entire of program is folloing that.
(Sorry for the comment written in Japanese, To tell the truth, this article is posted in other site in Japanese…)
(If you have any question, don’t hesitate to ask me :))
###########################################################
#パスの数
N <- 500
#一年あたりの刻み数
M <- 250
#割引債の満期
T <- 2
#金利モデルとデフォルト強度過程の相関係数
rho <- 0.3
###########################################################
#金利モデルのパラメータ
k.ir <- 0.6
theta.ir <- 0.05
vola.ir <- 0.05
init.ir <- 0.05
#デフォルト強度過程のパラメータ
k.default <- 0.6
theta.default <- 0.2
vola.default <- 0.5
init.default <- 0.2
###########################################################
#単位時間
dt <- 1 / M
#損失率(=1-回収率)
L <- function(t){0.3}
#相関付きの2変量標準正規分布を生成するための関数
#(外部パッケージを使用したくなかったので作成)
CorrelatedTwoVariableNormalRand <- function(n, rho)
{
result <- matrix(0, nrow = n, ncol = 2)
result[, 1] <- rnorm(n)
result[, 2] <- rho * result[,1] + sqrt(1 - rho^2) * rnorm(n)
result
}
#CIRモデルに従うパス生成機の生成機
CIRFactory <- function(k, theta, vola, dt)
{
#CIR model
function(x, dw){
x.sqrt <- ifelse(x >= 0, sqrt(x), 0)
x + k * (theta - x) * dt + vola * x.sqrt * sqrt(dt) * dw
}
}
#1パスでの価格を評価
EvaluateOnePath <- function(mma, index.default, dt)
{
size <- length(mma)
if(index.default > size){
#デフォルトなし(満期時点のMMA)
mma[size]
}else{
#デフォルト時点のMMA × 回収率(1 - 損失率)
mma[index.default] * (1 - L(index.default * dt))
}
}
#与えられたモデルから1パスを作る関数
CreateOnePath <- function(model, rand.numbers, init)
{
Reduce(function(x,y)model(x, y), rand.numbers, init, accumulate = TRUE)
}
#デフォルト時点(インデックス)の算出
DefaultIndex <- function(path.default, dt)
{
#指数分布に従う乱数(閾値用)
threshold <- rexp(1)
#デフォルト強度の累積
path.default.cum <- cumsum(path.default) * dt
#デフォルト時点(インデックスベース)の算出
#デフォルト強度の累積値がある閾値(指数分布に従う乱数)を超えた場合に倒産と判定
default.points <- which(path.default.cum >= threshold)
#デフォルトしてる場合minで箇所特定、してない時「満期+1時点」のインデックス返す
ifelse(length(default.points) > 0, min(default.points), length(path.default) + 1)
}
###########################################################
#金利&デフォルトモデル生成
model.ir <- CIRFactory(k.ir, theta.ir, vola.ir, dt)
model.default <- CIRFactory(k.default, theta.default, vola.default, dt)
#相関付きの乱数の生成
random.numbers <- CorrelatedTwoVariableNormalRand(N*T*M, rho)
#金利モデルを使ってパス&MMA(マネーマーケットアカウント)生成
random.number <- matrix(random.numbers[, 1], nrow = N, ncol = T * M)
path.ir <- t(apply(random.number, 1,function(x)CreateOnePath(model.ir, x, init.ir)))
path.ir <- path.ir[, -1]
mma <- exp(- t(apply(path.ir, 1, cumsum) * dt))
#デフォルトモデルを使ってデフォルト強度過程の生成
random.number <- matrix(random.numbers[, 2], nrow = N, ncol = T * M)
path.default <- t(apply(random.number, 1, function(x)CreateOnePath(model.default, x, init.default)))
path.default <- path.default[, -1]
#デフォルト時点の計算
index.default <- apply(path.default, 1, function(x)DefaultIndex(x, dt))
#各パスの評価+平均=信用リスク込みの割引債価格(モンテカルロ)
price <- sapply(1:N, function(i)EvaluateOnePath(mma[i, ], index.default[i], dt))
result <- mean(price)
result


The result of this program is about 0.830±0.005(It depends on the seed of random number generator).
Enjoy !

To leave a comment for the author, please follow the link and comment on their blog: My Life as a Mock Quant in English.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)