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”.Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
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
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
- kappa : 0.6
- theta : 0.05
- sigma : 0.05
- initial value : 0.05
- 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 :))
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
########################################################### | |
#パスの数 | |
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.