PD and IR

# visit http://reakkt.com for more

# pd <- 0.01  # probability of default

# cr <- 0     # (annual) coupon rate

fa <- 100   # face value

pd <- 0.05  # probability of default

rv <- 0.37   # recovery rate

rf <- 0.00429  # risk free rate

rl <- 1     # recovery time

n  <- 3     # years

# expected (discounted) return
pv <- function(fa,n,cr,rf) {

  -fa+sum(sapply(1:n, function(i) (fa*cr)/(1+rf)^i))+fa/(1+rf)^n

}

## expected (discounted) return at default
# di - default time/year
# rv - recovery rate
# rl - recovery delay/time
de <- function(fa,di,cr,rv,rf,rl=1) {

  rf.amount <- (fa*rv)/(1+rf)^(di+rl)

  if (di==0)

    # default in year 0
    return(-fa+rf.amount) else

    # default after first year
    return(-fa+sum(sapply(1:di, function(i) (fa*cr)/(1+rf)^i))+rf.amount)

}

pv(fa=100,n=3,cr=0.05,rf)

de(fa=100,di=0,cr=0.05,rv=0.4,rf,rl=1) # default in year one; no coupon payments
de(fa=100,di=1,cr=0.05,rv=0.4,rf,rl=1) # default after first coupon payment
de(fa=100,di=2,cr=0.05,rv=0.4,rf,rl=1) # default after second coupon payment
de(fa=100,di=3,cr=0.05,rv=0.4,rf,rl=1) # default after third coupon payment

# k  - initial number of entities
# pd - annual/periodic PD
# n  - years/periods
PortfolioSurvival <- function(k,pd,n) {

  defaulted  <- numeric(n+1)
  portfolio  <- rep(k,n+1)

  for (i in 0:n) {

    if (i==0) defaulted[1] <- ceiling(k*pd) else defaulted[i+1] <- ceiling(portfolio[i]*pd)

    portfolio[i+1] <- ifelse(i==0,k,portfolio[i])-defaulted[i+1]

  }

  list(defaulted=defaulted,
       portfolio=portfolio)

}

# k  - no. of assets
# pd - probability of default
# n  - no. of years/periods
RateGivenPD <- function(k,pd,n,
                        fa,cr,rf,rv,rl) {

  ps <- PortfolioSurvival(k,pd,n)

  portfolio <- ps$portfolio
  defaulted <- ps$defaulted

  cr   <- 0
  retS <- -1

  while(retS<0) {

    cr <- cr+0.001

    ret <- numeric(n+1)

    for (i in 0:n) {

      if (i!=n) {

        ret[i+1] <- de(fa,i,cr,rv,rf,rl)*defaulted[i+1]

      } else {

        ret[i+1] <- de(fa,i,cr,rv,rf,rl)*defaulted[i+1]+pv(fa,n,cr,rf)*portfolio[n+1]

      }

    }

    ret
    k*fa
    retS <- sum(ret)

  }

  return( list(cr=cr,
               ps=ps,
               ret=ret,
               retS=retS) )

}

##

k  <- 1e6   # no. of assets

pd.seq <- seq(0,0.99,0.01)

cr.arr <- numeric(length(pd.seq))

for (pd in pd.seq) cr.arr[which(pd.seq==pd)] <- RateGivenPD(k,pd,n,fa,cr,rf,rv,rl)$cr

plot(pd.seq,cr.arr,type="l",
     log="y",
     xlab="PD",ylab="IR (log)")
# abline(h=rf,col="Green",lty="dashed")

cbind(pd.seq,cr.arr)