Simple example:How to use foreach and doSNOW packages for parallel computation.

update

************************************************************************************************

I checked whether this example was run collectly or not in Windows XP(32bit) only ! 

************************************************************************************************

In R language, the members at Revolution R provide foreach and doSNOW packages for parallel computation. these packages allow us to compute things in parallel. So, we start to install these packages.

install.packages("foreach")
install.packages("doSNOW")

Created by Pretty R at inside-R.org

In foreach package, you can write the codes which are run not only in parallel but also in sequence. And, these are as following.

library(foreach)
#we get result as list
foreach(i = 1:3) %do% {sqrt(i)}
#we get result as vector with using .combine="c" option
foreach(i = 1:3,.combine = "c") %do% {sqrt(i)}
#if a result is "vector",we can get it as matrix with using .combine="cbind" option
foreach(i = 1:3,.combine = "cbind") %do% {letters[1:4]}
#if you define a function,you can use it as .combine option
#I wrote my function as returning same result that specify .combine="c" 
MyFunc <- function(x,y)c(x,y)
foreach(i = 1:3, .combine = "MyFunc") %do% {
  sqrt(i)
}

Created by Pretty R at inside-R.org

Next, we make clusters by doSNOW package for the purpose of parallel computation.
Because I have dual core machine, I specify two as the number of clusters.

> library(doSNOW)
> getDoParWorkers()
[1] 1
> getDoParName()
NULL
> registerDoSNOW(makeCluster(2, type = "SOCK"))
> getDoParWorkers()
[1] 2
> getDoParName()
[1] "doSNOW"
> getDoParVersion()
[1] "1.0.3"

Created by Pretty R at inside-R.org

Now, We are ready to compute things in parallel. It is easy for us to do that by foreach package. You only have to change "%do%" into "%dopar%". I compared the performance of parallel comutation to single computation as following.

> N <- 10^4
> system.time(foreach(i = 1:N,.combine = "cbind") %do% {
+   sum(rnorm(N))
+ })
   ユーザ   システム       経過
     57.52       0.48      59.60
> system.time(foreach(i = 1:N,.combine = "cbind") %dopar% {
+   sum(rnorm(N))
+ })
   ユーザ   システム       経過
     18.61       0.58      37.74

Created by Pretty R at inside-R.org

(I'm sorry that some terms are written in Japanese!)

You can understand the result of parallel computation is about twice as fast as single computation do !!!


Reference（including PDF）
-http://cran.r-project.org/web/packages/foreach/foreach.pdf
-http://cran.r-project.org/web/packages/foreach/vignettes/foreach.pdf
-http://cran.r-project.org/web/packages/foreach/vignettes/nested.pdf

Valuation of CDO with equal amount

Bank of Japan(BOJ) publish research paper regularly.
And, they issued very interesting paper about valuation of CDO recently.

(The paper is 金融危機時における資産価格変動の相互依存関係：コピュラに基づく評価 新谷　幸平、山田　哲也、吉羽　要直(sorry,japanese only!))

They introduced copula for pricing of CDO,and discussed how different CDO spreads were with using different copula for pricing.
I would like to reproduce their result (especially,P23-Table7)

The condition of calculation is following that

• number of debt（NUM.REFDEBT）:=100
• maturity(MATURITY):=5 year
• recovery rate(RECOVERY.RATE):=40%（constant value）
• probability of default (DEFAULT.PROBABILITY):=5％（in 5 years）
• parameter of nomal copula ρ:=0.15
• parameter of clayton copula α:=0.21

They apporoximated their valuation formula for easy calculation(equation (27))

(They assumed that CDO spread were paid as discounted bond at the begging.)

I simulated valuation of CDO with their method.

The result is following that

copula/tranche	Equity	mezzanine	senior	super senior
normal	1,145.42	62.49	0.52	0.000
t(20)	1,055.28	86.07	2.18	0.004
t(6)	896.74	126.44	8.56	0.044
t(3)	733.31	165.90	23.56	0.191
clayton	857.64	135.73	12.83	0.084

This table reproduce their result(P23-Table7).

And, In senior or super senior,you can understand that the CDO spread which is evaluated by normal copula is lower than the others. It means that normal copula is inadequate in financial crisis.

I show you my programming code（by R language).
If you copy and run my source code, you can duplicate my result easily.
Before you run, please install "copula"package.

library(copula)
#function for CDO spread calculation 
SpreadOfCDO <- function(copula, default.probability, maturity,
  recovery.rate, attachment, detachment, num.path, num.refdebt)
{
  random.copula <- rcopula(copula,num.path)
  num.default <- rowSums(random.copula < default.probability)
  loss.refdebt <- (1-recovery.rate)/num.refdebt*num.default
  loss.tranche <- (pmax(loss.refdebt - attachment,0)-pmax(loss.refdebt - detachment,0))/(detachment-attachment)
  expectation.loss.tranche <- sum(loss.tranche)/num.path
  spread <- -1/maturity*log(1-expectation.loss.tranche)
  return(spread)
}
################  main  ##################
#parameter
NUM.PATH <- 10^3
NUM.REFDEBT <- 100
DEFAULT.PROBABILITY <- 0.05
MATURITY <- 5
RECOVERY.RATE <- 0.4
#copulas which I want to compare.
COPULA <- list(normalCopula(0.15, dim = NUM.REFDEBT),
  tCopula(0.15, dim = NUM.REFDEBT, df = 20),
  tCopula(0.15, dim = NUM.REFDEBT, df = 6),
  tCopula(0.15, dim = NUM.REFDEBT, df = 3),
  claytonCopula(0.21, dim = NUM.REFDEBT)
)
#define for easy programming
SpreadOfCDOWithFixedParameter <- function(copula,attachment, detachment){
  SpreadOfCDO(copula, DEFAULT.PROBABILITY, MATURITY, 
    RECOVERY.RATE, attachment, detachment, NUM.PATH, NUM.REFDEBT)
}
result <- list()
#tranche:equity
result[[1]] <- sapply(COPULA,SpreadOfCDOWithFixedParameter,0.0,0.06)
#tranche:mezzanine
result[[2]] <- sapply(COPULA,SpreadOfCDOWithFixedParameter,0.06,0.18)
#tranche:senior
result[[3]] <- sapply(COPULA,SpreadOfCDOWithFixedParameter,0.18,0.36)
#tranche:super senior
result[[4]] <- sapply(COPULA,SpreadOfCDOWithFixedParameter,0.36,1)
#convert to matrix, and chenge unit to "bp"
result <- 10^4*do.call("cbind", result) 
colnames(result) <- c("equity","mezzanine","senior","super senior")
rownames(result) <- c("normal","t(20)","t(6)","t(3)","clayton")
result

Created by Pretty R at inside-R.org

Principal component analysis to yield curve change

In quantitive finance,it is often said that yield curve change is explained by three factor,
"parallel shift", "twist" and "butterfly".
Because I found that we can get historical yield curve data from FRB's web site, I check whether these proverbial facts are correct or not.Yield curve data can be downloaded to click "Go to download" and "Download File" button. Default data format is csv. If you would like to get data another format, you should click "Build package" button to change format.


I assume that downloaded data is located at  "C:\tmp"

#load data
term.structure <- read.csv("C:\\tmp\\FRB_H15.csv",stringsAsFactors=FALSE)
#use nearest 1000days data only.
term.structure <- tail(term.structure,1000)
#First column is "DATE".I don't need it.
term.structure <- term.structure[,-1]
label.term <- c("1M","3M","6M","1Y","2Y","3Y","5Y","7Y","10Y","20Y","30Y")
colnames(term.structure) <- label.term
#some rows have invalid value. I erase that's row
term.structure　<- subset(term.structure,term.structure$'1M' != "ND")
term.structure <- apply(term.structure,2,as.numeric)
#calculate diff.
term.structure.diff <- diff(term.structure)

Now,I have gotten yield curve change data.

Next,do principal component analysis and plot the result.

term.structure.princomp<- princomp(term.structure.diff)
factor.loadings <- term.structure.princomp$loadings[,1:3]
legend.loadings <- c("First principal component","Second principal component","Third principal component")
par(xaxt="n")
matplot(factor.loadings,type="l",
  lwd=3,lty=1,xlab = "Term", ylab = "Factor loadings")
legend(4,max(factor.loadings),legend=legend.loadings,col=1:3,lty=1,lwd=3)
par(xaxt="s")
axis(1,1:length(label.term),label.term)

Result image is shown like below

These result imply that Each three principal component correspond to "parallel shift", "twist" and "butterfly".

Cumulative Proportion are shown by "summary" function.

> summary(term.structure.princomp)
Importance of components:
                          Comp.1    Comp.2     Comp.3     Comp.4     Comp.5      Comp.6      Comp.7      Comp.8
Standard deviation     0.2028719 0.1381839 0.06938957 0.05234510 0.03430404 0.022611518 0.016081738 0.013068448
Proportion of Variance 0.5862010 0.2719681 0.06857903 0.03902608 0.01676075 0.007282195 0.003683570 0.002432489
Cumulative Proportion  0.5862010 0.8581690 0.92674803 0.96577411 0.98253486 0.989817052 0.993500621 0.995933111

As a result, yield cuve change can be explained by three principal component.