Pricing defaultable discount bond with reduced form model

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


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

Changing world, Changing JGB term structure

Writing the article "How much does "Beta" change depending on time?", I learned how to create an animation by using R language. Then, I would like to continue do that in this article.

In this article, I visualize time series of JGB term structure Ministry of Finance Japan publishes because I'm japanese!

You can download these data from here. You can get daily yield curve data, but I visualized these as monthly data for simplicity.

You can easily understand how JGB term structure behave and Japanese yield gradually go down.

The source code to create it is below.

(To run, you need to install xts, animation package, and ImageMagick)

library(xts)
#Source of JGB curve
source.jgb <- NULL
source.jgb[[length(source.jgb) + 1]] <- "http://www.mof.go.jp/jgbs/reference/interest_rate/data/jgbcm_1974-1979.csv"
source.jgb[[length(source.jgb) + 1]] <- "http://www.mof.go.jp/jgbs/reference/interest_rate/data/jgbcm_1980-1989.csv"
source.jgb[[length(source.jgb) + 1]] <- "http://www.mof.go.jp/jgbs/reference/interest_rate/data/jgbcm_1990-1999.csv"
source.jgb[[length(source.jgb) + 1]] <- "http://www.mof.go.jp/jgbs/reference/interest_rate/data/jgbcm_2000-2009.csv"
source.jgb[[length(source.jgb) + 1]] <- "http://www.mof.go.jp/jgbs/reference/interest_rate/data/jgbcm_2010.csv"
source.jgb[[length(source.jgb) + 1]] <- "http://www.mof.go.jp/jgbs/reference/interest_rate/jgbcm.csv"
#From Japanese era To ChristianEra
ToChristianEra <- function(x)
{
  era  <- substr(x, 1, 1)
  year <- as.numeric(substr(x, 2, nchar(x)))
  if(era == "H"){
    year <- year + 1988
  }else if(era == "S"){
    year <- year + 1925
  }
  as.character(year)
}
#Down load yield curve and convert to xts object
GetJGBYield <- function(source.url)
{
  jgb <- read.csv(source.url, stringsAsFactors = FALSE)
  #Extract date only
  jgb.day  <- strsplit(jgb[, 1], "\\.")
  #stop warning
  warn.old <- getOption("warn")
  options(warn = -1)
  #From Japanese era To ChristianEra
  jgb.day <- lapply(jgb.day, function(x)c(ToChristianEra(x[1]), x[2:length(x)]))
  #From date string to date object
  jgb[,  1] <- as.Date(sapply(jgb.day, function(x)Reduce(function(y,z)paste(y,z, sep="-"),x)))
  #Convert data from string to numeric
  jgb[, -1] <- apply(jgb[, -1], 2, as.numeric)
  options(warn = warn.old)
  as.xts(read.zoo(jgb))
}
#Down load JBG yield
jgb.list <- lapply(source.jgb, GetJGBYield)
#convert one xts object
jgb.xts <- Reduce(rbind, jgb.list)
#Interpolate(nearest value)
coredata(jgb.xts) <- na.locf(t(na.locf(t(coredata(jgb.xts)))))
#to monthly
jgb.xts <- jgb.xts[endpoints(jgb.xts, on="months",k = 1)]
 
#Label for x-axis
label.term <- paste(gsub("X", "", colnames(jgb.xts)), "Y", sep="")
#The range of y 
y.max <- c(min(jgb.xts), max(jgb.xts))
#plot one image
Snap <- function(val){
  term.structure <- coredata(val)
  index.date     <- index(val)
  par(xaxt="n")
  plot(t(term.structure),type="l",lwd=3, col = 2, xlab = "Term", ylab = "Rate", ylim = y.max)
  par(xaxt="s")
  axis(1, 1:length(label.term), label.term)
  text(0.5, y.max[2], as.character(index.date), pos = 4)
}
#save as animation
library(animation)
saveGIF({
  for(i in 1:nrow(jgb.xts)){Snap(jgb.xts[i])}
},interval = 0.005)

Created by Pretty R at inside-R.org