My Life as a Mock Quant in English

Pricing defaultable discount bond with reduced form model

2012-02-14T23:42:00.004+09:00

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

2011-11-11T08:30:00.003+09:00

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

How much does "Beta" change depending on time?

2011-11-02T07:45:00.000+09:00

You may often use "Beta" to measure the market exposure of your portfolio because it's easy to calculate.

Since I have been wondering how much "Beta" change depending on time, more precisely writing, data-set and the period of return time series, I think that I would like to write about that in this article.

First, I get stock price( I selected Mitsubishi UFJ Financial Group, Inc here) and market index(Nikkei225 ) from yahoo Japan by using RFinanceYJ package. These data period are from September, 2010 to September, 2011.

library(RFinanceYJ)
#download stock prices of nikkei-225 and MUFJ Financial group
mufg   <- quoteStockXtsData("8306.T",   since="2010-09-30",date.end="2011-09-30")$Close
nikkei <- quoteStockXtsData("998407.O", since="2010-09-30",date.end="2011-09-30")$Close

Created by Pretty R at inside-R.org

(you need to install RFinanceYJ Package to run this code)

I convert this price data into return data, and estimate the "Beta" of Mitsubishi UFJ Financial Group, Inc by using rolling linear regression(every regression have 125 sample data).

#convert to return
returns <- merge(mufg   / lag(mufg) - 1, nikkei / lag(nikkei) - 1)
names(returns) <- c("MUFG","NIKKEI")
#estimate beta by rolling regression
size.window <- 125
coefs <- rollapplyr(returns, size.window, function(x)coef(lm(x[,1]~x[,2])), by.column = FALSE)

Created by Pretty R at inside-R.org

As you know, R language provides us a very easy way to analyze, That is it.

So, Let's visualize these result. The result of simple plot is following that.

plot(as.xts(coefs[, 2]))

Created by Pretty R at inside-R.org

This graph shows that how historical beta changes depending on time.
And I create the animation of this result to understand more easily.

As you find, The result of linear regression strongly depends on the period of data set!
If we want to get a more robust result, we should you use something like "robust" regression :D

Here is the code to create this animation.

x <- na.omit(coredata(returns)[ ,2])
y <- na.omit(coredata(returns)[ ,1])
x.max <- c(-max(abs(x)), max(abs(x)))
y.max <- c(-max(abs(y)), max(abs(y)))
x.lab <- names(returns)[2]
y.lab <- names(returns)[1]
Snap <- function(val){
  val.x <- na.omit(coredata(val)[ ,2])
  val.y <- na.omit(coredata(val)[ ,1])
  lm.xy <- lm(val.y~val.x)
  plot(val.x, val.y, xlim = x.max, ylim = y.max, 
       xlab = x.lab, ylab = y.lab)
  abline(lm.xy)
  text(x.max[1], y.max[2], paste("Beta :", round(coef(lm.xy)[2],3)), pos = 4)
  text(x.max[1], y.max[1], as.character(last(index(val))), pos = 4)
}
 
library(animation)
saveGIF({
for(i in 1:(nrow(returns)-size.window)){
  Snap(returns[(i:(i+size.window)),])
}
},interval = 0.01)

Created by Pretty R at inside-R.org

you need to install "animation" package and ImageMagick (another software) to run this code

Enjoy!

Pair trading strategy : how to use "PairTrading" package

2011-10-25T20:14:00.002+09:00

Mr.Ishikawa(my old friend) and I developed "PairTrading" package, and uploaded it on CRAN.

This article shows you how you can use it.

The pair trading is a market neutral trading strategy and gives traders a chance to profit regardless of market conditions. The idea of this strategy is quite simple.
1 : Select two stocks(or any assets) moving similarly
2 : Short out-performing stock, buy under-performing one
3 : If "spread"(price difference between two stocks) converge, close your position.

So, Let's start to explain how to use this package.
(This example is in PDF manual of this package)

0: Install & load package
You can install and load "PairTrading" package via CRAN in the same way as other packages.

> install.packages("PairTrading")
> library(PairTrading)

1: Load sample data
We prepared sample stock price data in our package. You can load it by using "data" command.

> #load sample stock price data
> data(stock.price)

2: Estimate parameters
Next, We extract two stock prices(from 31 Mar, 2008) and estimate parameters.

> #select 2 stocks
> price.pair <- stock.price[,1:2]["2008-03-31::"]
> #Estimate parameters & plot spread
> reg <- EstimateParameters(price.pair, method = lm)

At the moment, we have only normal linear regression method to estimate parameters, but we will develop more sophisticated method in the future.
The estimation result contains the following contents.

> str(reg)
List of 3
 $ spread     :An ‘xts’ object from 2008-03-31 to 2011-08-05 containing:
  Data: num [1:821, 1] 0.26 0.271 0.294 0.275 0.255 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : chr "7203"
  Indexed by objects of class: [Date] TZ: 
  xts Attributes:  
 NULL
 $ hedge.ratio: num 0.285
 $ premium    : num 6.34

The most important thing in this estimation is "spread", then we try to plot it.

> plot(reg$spread)

And, you can check the stationarity of that by using "IsStationary" function.
This function return the result of two types unit root test.

> #check stationarity
> IsStationary(reg$spread, 0.1)
 PP.test adf.test 
    TRUE     TRUE

( augmented Dickey–Fuller test (ADF) and Phillips-Perron test)

3: Estimate parameters for back-test
To run back-test, you have to estimate parameters historically by using "EstimateParametersHistorically"function. This function do something like "rolling regression" to estimate parameters. This point is different from "EstimateParameter" function.

> #estimate parameters for back test
> params <- EstimateParametersHistorically(price.pair, period = 180)
> #create & plot trading signals
> str(params)
An ‘xts’ object from 2008-12-18 to 2011-08-05 containing:
  Data: num [1:642, 1:3] 0.065 0.0577 0.0396 0.0136 0.021 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : chr [1:3] "spread" "hedge.ratio" "premium"
  Indexed by objects of class: [Date] TZ: 
  xts Attributes:  
 NULL

4: Create trading signal

Next, you create trading singal using estimated spread. "Simple" function give a very simple trading strategy(If The spread is more(less) than specified value, you will buy(sell))

> #create & plot trading signals
> signal <- Simple(params$spread, 0.05)
> barplot(signal,col="blue",space = 0, border = "blue",xaxt="n",yaxt="n",xlab="",ylab="")
> par(new=TRUE)
> plot(params$spread)

In this case, The trading signal is drawn as below.

ading signal is
5: Back-test performance

Last, you can check the performance of pair trading by using "Return" function.

> #performance
> return.pairtrading <- Return(price.pair, lag(signal), lag(params$hedge.ratio))
> plot(100 * cumprod(1 + return.pairtrading))

In this case, our strategy seems to work correctly :-)

6: Conclusion and remarks

Pair trading is well-known trading strategy, and I introduced "PairTrading" package in this article.

We would like to modify this package to be more useful and fit in real-market.

If you have any suggestion, please let me know.

And we created a presentation slide to explain the basic concept of pair trading.

It may be useful for you to understand the basic concept of pair trading if you are interested in it.

Introduction to pairtrading

View more presentations from Kohta Ishikawa

Enjoy!