From 2ca4f88a4fb14026b502ee7e54ebdc12d2b4b8b9 Mon Sep 17 00:00:00 2001 From: Guillaume Horel Date: Thu, 24 Apr 2014 13:43:37 -0400 Subject: initial import of lossdistrib package --- DESCRIPTION | 9 + NAMESPACE | 2 + R/tranche_functions.R | 879 ++++++++++++++++++++++++++++++++++++++++++++++++++ src/Makevars | 2 + src/Makevars.win | 2 + src/lossdistrib.c | 760 +++++++++++++++++++++++++++++++++++++++++++ src/lossdistrib.h | 46 +++ 7 files changed, 1700 insertions(+) create mode 100644 DESCRIPTION create mode 100644 NAMESPACE create mode 100644 R/tranche_functions.R create mode 100644 src/Makevars create mode 100644 src/Makevars.win create mode 100644 src/lossdistrib.c create mode 100644 src/lossdistrib.h diff --git a/DESCRIPTION b/DESCRIPTION new file mode 100644 index 0000000..166cebd --- /dev/null +++ b/DESCRIPTION @@ -0,0 +1,9 @@ +Package: lossdistrib +Type: Package +Title: What the package does (short line) +Version: 1.0 +Date: 2014-04-24 +Author: Who wrote it +Maintainer: Who to complain to +Description: More about what it does (maybe more than one line) +License: What license is it under? diff --git a/NAMESPACE b/NAMESPACE new file mode 100644 index 0000000..f2e6748 --- /dev/null +++ b/NAMESPACE @@ -0,0 +1,2 @@ +exportPattern("^[[:alpha:]]+") +useDynLib("lossdistrib") diff --git a/R/tranche_functions.R b/R/tranche_functions.R new file mode 100644 index 0000000..2e7f0f3 --- /dev/null +++ b/R/tranche_functions.R @@ -0,0 +1,879 @@ +library(statmod) + +## todo: +## -investigate other ways to interpolate the random severities on the grid +## I'm thinking that at eah severity that we add to the distribution, round it down +## and keep track of the missing mass: namely if X_i=S_i w.p p_i, then add +## X_i=lu*floor(S_i/lu) with probability p_i and propagate +## X_{i+1}=S_{i+1}+(S_i-lu*floor(S_i/lu)) with the right probability +## - investigate truncated distributions more (need to compute loss and recov distribution +## separately, for the 0-10 equity tranche, we need the loss on the 0-0.1 support and +## recovery with 0.1-1 support, so it's not clear that there is a big gain. +## - do the correlation adjustments when computing the deltas since it seems to be +## the market standard +hostname <- system("hostname", intern=TRUE) + +checkSymbol <- function(name){ + if(!is.loaded(name)){ + if(.Platform$OS.type == "unix"){ + root.dir <- "/home/share/CorpCDOs" + dyn.load(file.path(root.dir, "code", "R", paste0("lossdistrib", + hostname, + .Platform$dynlib.ext))) + }else{ + root.dir <- "//WDSENTINEL/share/CorpCDOs" + dyn.load(file.path(root.dir, "code", "R", paste0("lossdistrib", + .Platform$dynlib.ext))) + } + } +} +lossdistrib <- function(p){ + ## basic recursive algorithm of Andersen, Sidenius and Basu + n <- length(p) + q <- rep(0, (n+1)) + q[1] <- 1 + for(i in 1:n){ + q[-1] <- p[i]*q[-(n+1)]+(1-p[i])*q[-1] + q[1] <- (1-p[i])*q[1] + } + return(q) +} + +lossdistrib.fft <- function(p){ + ## computes loss distribution using the fft + ## complexity is of order O(n*m)+O(m*log(m)) + ## where m is the size of the grid and n the number of probabilities. + ## this is slower than the recursive algorithm + theta <- 2*pi*1i*(0:n)/(n+1) + Phi <- 1 - p + p%o%exp(theta) + v <- apply(Phi, 2, prod) + return(1/(n+1)*Re(fft(v))) +} + +lossdistrib2 <- function(p, w, S, N, defaultflag=FALSE){ + ## recursive algorithm with first order correction + ## p vector of default probabilities + ## w vector of weigths + ## S vector of severities + ## N number of ticks in the grid + ## defaultflag if true computes the default distribution + n <- length(p) + lu <- 1/(N-1) + q <- rep(0, N) + q[1] <- 1 + for(i in 1:n){ + if(defaultflag){ + d <- w[i] /lu + }else{ + d <- S[i] * w[i] / lu + } + d1 <- floor(d) + d2 <- ceiling(d) + p1 <- p[i]*(d2-d) + p2 <- p[i] - p1 + q1 <- c(rep(0,d1), p1*q[1:(N-d1)]) + q2 <- c(rep(0,d2), p2*q[1:(N-d2)]) + q <- q1 + q2 + (1-p[i])*q + } + q[length(q)] <- q[length(q)]+1-sum(q) + return(q) +} + +lossdistrib2.truncated <- function(p, w, S, N, cutoff=N){ + ## recursive algorithm with first order correction + ## p vector of default probabilities + ## w vector of weigths + ## S vector of severities + ## N number of ticks in the grid (for best accuracy should + ## be a multiple of the number of issuers) + ## cutoff where to stop computing the exact probabilities + ## (useful for tranche computations) + + ## this is actually slower than lossdistrib2. But in C this is + ## twice as fast. + ## for high severities, M can become bigger than N, and there is + ## some probability mass escaping. + n <- length(p) + lu <- 1/(N-1) + q <- rep(0, truncated) + q[1] <- 1 + M <- 1 + for(i in 1:n){ + d <- S[i] * w[i] / lu + d1 <- floor(d) + d2 <- ceiling(d) + p1 <- p[i]*(d2-d) + p2 <- p[i] - p1 + q1 <- p1*q[1:min(M, cutoff-d1)] + q2 <- p2*q[1:min(M, cutoff-d2)] + q[1:min(M, cutoff)] <- (1-p[i])*q[1:min(M, cutoff)] + q[(d1+1):min(M+d1, cutoff)] <- q[(d1+1):min(M+d1, cutoff)]+q1 + q[(d2+1):min(M+d2, cutoff)] <- q[(d2+1):min(M+d2, cutoff)]+q2 + M <- M+d2 + } + return(q) +} + +recovdist <- function(dp, pp, w, S, N){ + ## computes the recovery distribution for a sum of independent variables + ## R=\sum_{i=1}^n X_i + ## where X_i = 0 w.p 1-dp[i]-pp[i] + ## = w[i]*(1-S[i]) w.p dp[i] + ## = w[i] w.p pp[i] + ## each non zero value v is interpolated on the grid as + ## the pair of values floor(v/lu) and ceiling(v/lu) so that + ## X_i has four non zero values + n <- length(dp) + q <- rep(0, N) + q[1] <- 1 + lu <- 1/(N-1) + for(i in 1:n){ + d1 <- w[i]*(1-S[i])/lu + d1l <- floor(d1) + d1u <- ceiling(d1) + d2 <- w[i] / lu + d2l <- floor(d2) + d2u <- ceiling(d2) + dp1 <- dp[i] * (d1u-d1) + dp2 <- dp[i] - dp1 + pp1 <- pp[i] * (d2u - d2) + pp2 <- pp[i] - pp1 + q1 <- c(rep(0, d1l), dp1 * q[1:(N-d1l)]) + q2 <- c(rep(0, d1u), dp2 * q[1:(N-d1u)]) + q3 <- c(rep(0, d2l), pp1 * q[1:(N-d2l)]) + q4 <- c(rep(0, d2u), pp2 *q[1:(N-d2u)]) + q <- q1+q2+q3+q4+(1-dp[i]-pp[i])*q + } + return(q) +} + +lossdist.joint <- function(p, w, S, N, defaultflag=FALSE){ + ## recursive algorithm with first order correction + ## to compute the joint probability distribution of the loss and recovery + ## inputs: + ## p: vector of default probabilities + ## w: vector of issuer weights + ## S: vector of severities + ## N: number of tick sizes on the grid + ## defaultflag: if true computes the default distribution + ## output: + ## q: matrix of joint loss, recovery probability + ## colSums(q) is the recovery distribution marginal + ## rowSums(q) is the loss distribution marginal + n <- length(p) + lu <- 1/(N-1) + q <- matrix(0, N, N) + q[1,1] <- 1 + for(k in 1:n){ + if(defaultflag){ + x <- w[k] / lu + }else{ + x <- S[k] * w[k]/lu + } + y <- (1-S[k]) * w[k]/lu + i <- floor(x) + j <- floor(y) + weights <- c((i+1-x)*(j+1-y), (i+1-x)*(y-j), (x-i)*(y-j), (j+1-y)*(x-i)) + psplit <- p[k] * weights + qtemp <- matrix(0, N, N) + qtemp[(i+1):N,(j+1):N] <- qtemp[(i+1):N,(j+1):N] + psplit[1] * q[1:(N-i),1:(N-j)] + qtemp[(i+1):N,(j+2):N] <- qtemp[(i+1):N,(j+2):N] + psplit[2] * q[1:(N-i), 1:(N-j-1)] + qtemp[(i+2):N,(j+2):N] <- qtemp[(i+2):N,(j+2):N] + psplit[3] * q[1:(N-i-1), 1:(N-j-1)] + qtemp[(i+2):N, (j+1):N] <- qtemp[(i+2):N, (j+1):N] + psplit[4] * q[1:(N-i-1), 1:(N-j)] + q <- qtemp + (1-p[k])*q + } + q[length(q)] <- q[length(q)]+1-sum(q) + return(q) +} + +lossdist.prepay.joint <- function(dp, pp, w, S, N, defaultflag=FALSE){ + ## recursive algorithm with first order correction + ## to compute the joint probability distribition of the loss and recovery + ## inputs: + ## dp: vector of default probabilities + ## pp: vector of prepay probabilities + ## w: vector of issuer weights + ## S: vector of severities + ## N: number of tick sizes on the grid + ## defaultflag: if true computes the default + ## outputs + ## q: matrix of joint loss and recovery probability + ## colSums(q) is the recovery distribution marginal + ## rowSums(q) is the loss distribution marginal + n <- length(dp) + lu <- 1/(N-1) + q <- matrix(0, N, N) + q[1,1] <- 1 + for(k in 1:n){ + y1 <- (1-S[k]) * w[k]/lu + y2 <- w[k]/lu + j1 <- floor(y1) + j2 <- floor(y2) + if(defaultflag){ + x <- y2 + i <- j2 + }else{ + x <- y2-y1 + i <- floor(x) + } + + ## weights <- c((i+1-x)*(j1+1-y1), (i+1-x)*(y1-j1), (x-i)*(y1-j1), (j1+1-y1)*(x-i)) + weights1 <- c((i+1-x)*(j1+1-y1), (i+1-x)*(y1-j1), (x-i)*(y1-j1), (j1+1-y1)*(x-i)) + dpsplit <- dp[k] * weights1 + + if(defaultflag){ + weights2 <- c((i+1-x)*(j2+1-y2), (i+1-x)*(y2-j2), (x-i)*(y2-j2), (j2+1-y2)*(x-i)) + ppsplit <- pp[k] * weights2 + }else{ + ppsplit <- pp[k] * c(j2+1-y2, y2-j2) + } + qtemp <- matrix(0, N, N) + qtemp[(i+1):N,(j1+1):N] <- qtemp[(i+1):N,(j1+1):N] + dpsplit[1] * q[1:(N-i),1:(N-j1)] + qtemp[(i+1):N,(j1+2):N] <- qtemp[(i+1):N,(j1+2):N] + dpsplit[2] * q[1:(N-i), 1:(N-j1-1)] + qtemp[(i+2):N,(j1+2):N] <- qtemp[(i+2):N,(j1+2):N] + dpsplit[3] * q[1:(N-i-1), 1:(N-j1-1)] + qtemp[(i+2):N,(j1+1):N] <- qtemp[(i+2):N, (j1+1):N] + dpsplit[4] * q[1:(N-i-1), 1:(N-j1)] + if(defaultflag){ + qtemp[(i+1):N,(j2+1):N] <- qtemp[(i+1):N,(j2+1):N] + ppsplit[1] * q[1:(N-i),1:(N-j2)] + qtemp[(i+1):N,(j2+2):N] <- qtemp[(i+1):N,(j2+2):N] + ppsplit[2] * q[1:(N-i), 1:(N-j2-1)] + qtemp[(i+2):N,(j2+2):N] <- qtemp[(i+2):N,(j2+2):N] + ppsplit[3] * q[1:(N-i-1), 1:(N-j2-1)] + qtemp[(i+2):N,(j2+1):N] <- qtemp[(i+2):N, (j2+1):N] + ppsplit[4] * q[1:(N-i-1), 1:(N-j2)] + }else{ + qtemp[, (j2+1):N] <- qtemp[,(j2+1):N]+ppsplit[1]*q[,1:(N-j2)] + qtemp[, (j2+2):N] <- qtemp[,(j2+2):N]+ppsplit[2]*q[,1:(N-j2-1)] + } + q <- qtemp + (1-pp[k]-dp[k]) * q + } + q[length(q)] <- q[length(q)] + 1 - sum(q) + return(q) +} + +lossdistC <- function(p, w, S, N, defaultflag=FALSE){ + ## C version of lossdistrib2, roughly 50 times faster + .C("lossdistrib", as.double(p), as.integer(length(p)), + as.double(w), as.double(S), as.integer(N), as.logical(defaultflag), q = double(N))$q +} + +lossdistCblas <- function(p, w, S, N, defaultflag=FALSE){ + ## C version of lossdistrib2, roughly 50 times faster + .C("lossdistrib_blas", as.double(p), as.integer(length(p)), + as.double(w), as.double(S), as.integer(N), as.logical(defaultflag), q = double(N))$q +} + +lossdistCZ <- function(p, w, S, N, defaultflag=FALSE, rho, Z, wZ){ + .C("lossdistrib_Z", as.double(p), as.integer(length(p)), + as.double(w), as.double(S), as.integer(N), as.logical(defaultflag), + as.double(rho), as.double(Z), as.integer(length(Z)), + q = matrix(0, N, length(Z)))$q +} + +lossdistC.truncated <- function(p, w, S, N, T=N){ + ## C version of lossdistrib2, roughly 50 times faster + .C("lossdistrib_truncated", as.double(p), as.integer(length(p)), + as.double(w), as.double(S), as.integer(N), as.integer(T), q = double(T))$q +} + +recovdistC <- function(dp, pp, w, S, N){ + ## C version of recovdist + .C("recovdist", as.double(dp), as.double(pp), as.integer(length(dp)), + as.double(w), as.double(S), as.integer(N), q = double(N))$q +} + +lossdistC.joint <- function(p, w, S, N, defaultflag=FALSE){ + ## C version of lossdistrib.joint, roughly 20 times faster + .C("lossdistrib_joint", as.double(p), as.integer(length(p)), as.double(w), + as.double(S), as.integer(N), as.logical(defaultflag), q = matrix(0, N, N))$q +} + +lossdistC.jointblas <- function(p, w, S, N, defaultflag=FALSE){ + ## C version of lossdistrib.joint, roughly 20 times faster + .C("lossdistrib_joint_blas", as.double(p), as.integer(length(p)), as.double(w), + as.double(S), as.integer(N), as.logical(defaultflag), q = matrix(0, N, N))$q +} + + +lossdistC.jointZ <- function(dp, w, S, N, defaultflag = FALSE, rho, Z, wZ){ + ## N is the size of the grid + ## dp is of size n.credits + ## w is of size n.credits + ## S is of size n.credits by nZ + ## rho is a double + ## Z is a vector of length nZ + ## w is a vector if length wZ + r <- .C("lossdistrib_joint_Z", as.double(dp), as.integer(length(dp)), + as.double(w), as.double(S), as.integer(N), as.logical(defaultflag), as.double(rho), + as.double(Z), as.double(wZ), as.integer(length(Z)), q = matrix(0, N, N))$q +} + +lossdistC.prepay.joint <- function(dp, pp, w, S, N, defaultflag=FALSE){ + ## C version of lossdist.prepay.joint + r <- .C("lossdistrib_prepay_joint", as.double(dp), as.double(pp), as.integer(length(dp)), + as.double(w), as.double(S), as.integer(N), as.logical(defaultflag), q=matrix(0, N, N))$q + return(r) +} + +lossdistC.prepay.jointZ <- function(dp, pp, w, S, N, defaultflag = FALSE, rho, Z, wZ){ + ## N is the size of the grid + ## dp is of size n.credits + ## pp is of size n.credits + ## w is of size n.credits + ## S is of size n.credits by nZ + ## rho is a double + ## Z is a vector of length nZ + ## w is a vector if length wZ + r <- .C("lossdistrib_prepay_joint_Z", as.double(dp), as.double(pp), as.integer(length(dp)), + as.double(w), as.double(S), as.integer(N), as.logical(defaultflag), as.double(rho), + as.double(Z), as.double(wZ), as.integer(length(Z)), q = matrix(0, N, N))$q +} + +lossrecovdist <- function(defaultprob, prepayprob, w, S, N, defaultflag=FALSE, useC=TRUE){ + if(missing(prepayprob)){ + if(useC){ + L <- lossdistC(defaultprob, w, S, N, defaultflag) + R <- lossdistC(defaultprob, w, 1-S, N) + }else{ + L <- lossdistrib2(defaultprob, w, S, N, defaultflag) + R <- lossdistrib2(defaultprob, w, 1-S, N) + } + }else{ + if(useC){ + L <- lossdistC(defaultprob+defaultflag*prepayprob, w, S, N, defaultflag) + R <- recovdistC(defaultprob, prepayprob, w, S, N) + }else{ + L <- lossdistrib2(defaultprob+defaultflag*prepayprob, w, S, N, defaultflag) + R <- recovdist(defaultprob, prepayprob, w, S, N) + } + } + return(list(L=L, R=R)) +} + +lossrecovdist.term <- function(defaultprob, prepayprob, w, S, N, defaultflag=FALSE, useC=TRUE){ + ## computes the loss and recovery distribution over time + L <- array(0, dim=c(N, ncol(defaultprob))) + R <- array(0, dim=c(N, ncol(defaultprob))) + if(missing(prepayprob)){ + for(t in 1:ncol(defaultprob)){ + temp <- lossrecovdist(defaultprob[,t], , w, S[,t], N, defaultflag, useC) + L[,t] <- temp$L + R[,t] <- temp$R + } + }else{ + for(t in 1:ncol(defaultprob)){ + temp <- lossrecovdist(defaultprob[,t], prepayprob[,t], w, S[,t], N, defaultflag, useC) + L[,t] <- temp$L + R[,t] <- temp$R + } + } + return(list(L=L, R=R)) +} + +lossrecovdist.joint.term <- function(defaultprob, prepayprob, w, S, N, defaultflag=FALSE, useC=TRUE){ + ## computes the joint loss and recovery distribution over time + Q <- array(0, dim=c(ncol(defaultprob), N, N)) + if(useC){ + if(missing(prepayprob)){ + for(t in 1:ncol(defaultprob)){ + Q[t,,] <- lossdistC.joint(defaultprob[,t], w, S[,t], N, defaultflag) + } + }else{ + for(t in 1:ncol(defaultprob)){ + Q[t,,] <- lossdistC.prepay.joint(defaultprob[,t], prepayprob[,t], w, S[,t], N, defaultflag) + } + } + }else{ + if(missing(prepayprob)){ + for(t in 1:ncol(defaultprob)){ + Q[t,,] <- lossdist.joint(defaultprob[,t], w, S[,t], N, defaultflag) + } + }else{ + for(t in 1:ncol(defaultprob)){ + Q[t,,] <- lossdist.prepay.joint(defaultprob[,t], prepayprob[,t], w, S[,t], N, defaultflag) + } + } + } + gc() + return(Q) +} + +dist.transform <- function(dist.joint){ + ## compute the joint (D, R) distribution + ## from the (L, R) distribution using D = L+R + distDR <- array(0, dim=dim(dist.joint)) + Ngrid <- dim(dist.joint)[2] + for(t in 1:dim(dist.joint)[1]){ + for(i in 1:Ngrid){ + for(j in 1:Ngrid){ + index <- i+j + if(index <= Ngrid){ + distDR[t,index,j] <- distDR[t,index,j] + dist.joint[t,i,j] + }else{ + distDR[t,Ngrid,j] <- distDR[t,Ngrid,j] + + dist.joint[t,i,j] + } + } + } + distDR[t,,] <- distDR[t,,]/sum(distDR[t,,]) + } + return( distDR ) +} + +shockprob <- function(p, rho, Z, log.p=F){ + ## computes the shocked default probability as a function of the copula factor + ## function is vectorized provided the below caveats: + ## p and rho are vectors of same length n, Z is a scalar, returns vector of length n + ## p and rho are scalars, Z is a vector of length n, returns vector of length n + if(length(p)==1){ + if(rho==1){ + return(Z<=qnorm(p)) + }else{ + return(pnorm((qnorm(p)-sqrt(rho)*Z)/sqrt(1-rho), log.p=log.p)) + } + }else{ + result <- double(length(p)) + result[rho==1] <- Z<=qnorm(p[rho==1]) + result[rho<1] <- pnorm((qnorm(p[rho<1])-sqrt(rho[rho<1])*Z)/sqrt(1-rho[rho<1]), log.p=log.p) + return( result ) + } +} + +shockseverity <- function(S, Stilde=1, Z, rho, p){ + ## computes the severity as a function of the copula factor Z + result <- double(length(S)) + result[p==0] <- 0 + result[p!=0] <- Stilde * exp( shockprob(S[p!=0]/Stilde*p[p!=0], rho[p!=0], Z, TRUE) - + shockprob(p[p!=0], rho[p!=0], Z, TRUE)) + return(result) +} + +dshockprob <- function(p,rho,Z){ + dnorm((qnorm(p)-sqrt(rho)*Z)/sqrt(1-rho))*dqnorm(p)/sqrt(1-rho) +} + +dqnorm <- function(x){ + 1/dnorm(qnorm(x)) +} + +fit.prob <- function(Z, w, rho, p0){ + ## if the weights are not perfectly gaussian, find the probability p such + ## E_w(shockprob(p, rho, Z)) = p0 + require(distr) + if(p0==0){ + return(0) + } + if(rho == 1){ + distw <- DiscreteDistribution(Z, w) + return(pnorm(q(distw)(p0))) + } + eps <- 1e-12 + dp <- (crossprod(shockprob(p0,rho,Z),w)-p0)/crossprod(dshockprob(p0,rho,Z),w) + p <- p0 + while(abs(dp) > eps){ + dp <- (crossprod(shockprob(p,rho,Z),w)-p0)/crossprod(dshockprob(p,rho,Z),w) + phi <- 1 + while ((p-phi*dp)<0 || (p-phi*dp)>1){ + phi <- 0.8*phi + } + p <- p - phi*dp + } + return(p) +} + +fit.probC <- function(Z, w, rho, p0){ + r <- .C("fitprob", as.double(Z), as.double(w), as.integer(length(Z)), + as.double(rho), as.double(p0), q = double(1)) + return(r$q) +} + +stochasticrecov <- function(R, Rtilde, Z, w, rho, porig, pmod){ + ## if porig == 0 (probably matured asset) then return orginal recovery + ## it shouldn't matter anyway since we never default that asset + if(porig == 0){ + return(rep(R, length(Z))) + }else{ + ptilde <- fit.prob(Z, w, rho, (1-R)/(1-Rtilde) * porig) + return(abs(1-(1-Rtilde) * exp(shockprob(ptilde, rho, Z, TRUE) - shockprob(pmod, rho, Z, TRUE)))) + } +} + +stochasticrecovC <- function(R, Rtilde, Z, w, rho, porig, pmod){ + r <- .C("stochasticrecov", as.double(R), as.double(Rtilde), as.double(Z), + as.double(w), as.integer(length(Z)), as.double(rho), as.double(porig), + as.double(pmod), q = double(length(Z))) + return(r$q) +} + +pos <- function(x){ + pmax(x, 0) +} + +trancheloss <- function(L, K1, K2){ + pos(L - K1) - pos(L - K2) +} + +trancherecov <- function(R, K1, K2){ + pos(R - 1 + K2) - pos(R - 1 +K1) +} + +tranche.cl <- function(L, R, cs, K1, K2, Ngrid=nrow(L), scaled=FALSE){ + ## computes the couponleg of a tranche + ## if scaled is TRUE, scale it by the size of the tranche (K2-K1) + ## can make use of the fact that the loss and recov distribution are + ## truncated (in that case nrow(L) != Ngrid + if(K1==K2){ + return( 0 ) + }else{ + support <- seq(0, 1, length=Ngrid)[1:nrow(L)] + size <- K2 - K1 - crossprod(trancheloss(support, K1, K2), L) - + crossprod(trancherecov(support, K1, K2), R) + sizeadj <- as.numeric(0.5 * (size + c(K2-K1, size[-length(size)]))) + if(scaled){ + return( 1/(K2-K1) * crossprod(sizeadj * cs$coupons, cs$df) ) + }else{ + return( crossprod(sizeadj * cs$coupons, cs$df) ) + } + } +} + +tranche.cl.scenarios <- function(l, r, cs, K1, K2, scaled=FALSE){ + ## computes the couponleg of a tranche for one scenario + ## if scaled is TRUE, scale it by the size of the tranche (K2-K1) + ## can make use of the fact that the loss and recov distribution are + ## truncated (in that case nrow(L) != Ngrid + if(K1==K2){ + return( 0 ) + }else{ + size <- K2 - K1 - trancheloss(l, K1, K2) - trancherecov(r, K1, K2) + sizeadj <- as.numeric(0.5 * (size + c(K2-K1, size[-length(size)]))) + if(scaled){ + return( 1/(K2-K1) * crossprod(sizeadj * cs$coupons, cs$df) ) + }else{ + return( crossprod(sizeadj * cs$coupons, cs$df) ) + } + } +} + +funded.tranche.pv <- function(L, R, cs, K1, K2, scaled = FALSE){ + if(K1==K2){ + return(0) + }else{ + size <- K2 - K1 -trancheloss(L, K1, K2) - trancherecov(R, K1, K2) + sizeadj <- as.numeric(0.5 * (size + c(K2-K1, size[-length(size)]))) + interest <- crossprod(sizeadj * cs$coupons, cs$df) + principal <- diff(c(0, trancherecov(R, K1, K2))) + principal[length(principal)] <- principal[length(principal)] + size[length(size)] + principal <- crossprod(cs$df, principal) + if(scaled){ + pv <- (interest + principal)/(K2-K1) + }else{ + pv <- (interest + principal) + } + return(pv) + } +} + +tranche.pl <- function(L, cs, K1, K2, Ngrid=nrow(L), scaled=FALSE){ + ## computes the protection leg of a tranche + ## if scaled + if(K1==K2){ + return(0) + }else{ + support <- seq(0, 1, length=Ngrid)[1:nrow(L)] + cf <- K2 - K1 - crossprod(trancheloss(support, K1, K2), L) + cf <- c(K2 - K1, cf) + if(scaled){ + return( 1/(K2-K1) * crossprod(diff(cf), cs$df)) + }else{ + return( crossprod(diff(cf), cs$df)) + } + } +} + +tranche.pl.scenarios <- function(l, cs, K1, K2, scaled=FALSE){ + ## computes the protection leg of a tranche + ## if scaled + if(K1==K2){ + return(0) + }else{ + cf <- K2 - K1 - trancheloss(l, K1, K2) + cf <- c(K2 - K1, cf) + if(scaled){ + return( 1/(K2-K1) * as.numeric(crossprod(diff(cf), cs$df))) + }else{ + return( as.numeric(crossprod(diff(cf), cs$df))) + } + } +} + +tranche.pv <- function(L, R, cs, K1, K2, Ngrid=nrow(L)){ + return( tranche.pl(L, cs, K1, K2, Ngrid) + tranche.cl(L, R, cs, K1, K2, Ngrid)) +} + +tranche.pv.scenarios <- function(l, r, cs, K1, K2){ + return( tranche.pl.scenarios(l, cs, K1, K2, TRUE) + + tranche.cl.scenarios(l, r, cs, K1, K2, TRUE)) +} + + +adjust.attachments <- function(K, losstodate, factor){ + ## computes the attachments adjusted for losses + ## on current notional + return( pmin(pmax((K-losstodate)/factor, 0),1) ) +} + +tranche.pvvec <- function(K, L, R, cs){ + r <- rep(0, length(K)-1) + for(i in 1:(length(K)-1)){ + r[i] <- 1/(K[i+1]-K[i]) * tranche.pv(L, R, cs, K[i], K[i+1]) + } + return( r ) +} + +BClossdist <- function(defaultprob, issuerweights, recov, rho, Z, w, + N=length(recov)+1, defaultflag=FALSE, n.int=500){ + if(missing(Z)){ + quadrature <- gauss.quad.prob(n.int, "normal") + Z <- quadrature$nodes + w <- quadrature$weights + } + ## do not use if weights are not gaussian, results would be incorrect + ## since shockseverity is invalid in that case (need to use stochasticrecov) + LZ <- matrix(0, N, length(Z)) + RZ <- matrix(0, N, length(Z)) + L <- matrix(0, N, ncol(defaultprob)) + R <- matrix(0, N, ncol(defaultprob)) + for(t in 1:ncol(defaultprob)){ + for(i in 1:length(Z)){ + g.shocked <- shockprob(defaultprob[,t], rho, Z[i]) + S.shocked <- shockseverity(1-recov, 1, Z[i], rho, defaultprob[,t]) + temp <- lossrecovdist(g.shocked, , issuerweights, S.shocked, N) + LZ[,i] <- temp$L + RZ[,i] <- temp$R + } + L[,t] <- LZ%*%w + R[,t] <- RZ%*%w + } + list(L=L, R=R) +} + +BClossdistC <- function(defaultprob, issuerweights, recov, rho, Z, w, + N=length(issuerweights)+1, defaultflag=FALSE){ + L <- matrix(0, N, dim(defaultprob)[2]) + R <- matrix(0, N, dim(defaultprob)[2]) + rho <- rep(rho, length(issuerweights)) + r <- .C("BClossdist", defaultprob, dim(defaultprob)[1], dim(defaultprob)[2], + as.double(issuerweights), as.double(recov), as.double(Z), as.double(w), + as.integer(length(Z)), as.double(rho), as.integer(N), as.logical(defaultflag), L=L, R=R) + return(list(L=r$L,R=r$R)) +} + +BCtranche.pv <- function(defaultprob, issuerweights, recov, cs, K1, K2, rho1, rho2, + Z, w, N=length(issuerweights)+1){ + ## computes the protection leg, couponleg, and bond price of a tranche + ## in the base correlation setting + if(K1==0){ + if(rho1!=0){ + stop("equity tranche must have 0 lower correlation") + } + } + dK <- K2 - K1 + dist2 <- BClossdistC(defaultprob, issuerweights, recov, rho2, Z, w, N) + if(rho1!=0){ + dist1 <- BClossdistC(defaultprob, issuerweights, recov, rho1, Z, w, N) + } + cl2 <- tranche.cl(dist2$L, dist2$R, cs, 0, K2) + cl1 <- tranche.cl(dist1$L, dist1$R, cs, 0, K1) + pl2 <- tranche.pl(dist2$L, cs, 0, K2) + pl1 <- tranche.pl(dist1$L, cs, 0, K1) + return(list(pl=(pl2-pl1)/dK, cl=(cl2-cl1)/dK, + bp=100*(1+(pl2-pl1+cl2-cl1)/dK))) +} + +BCtranche.delta <- function(portfolio, index, coupon, K1, K2, rho1, rho2, Z, w, + N=length(portolio)+1, tradedate = Sys.Date()){ + ## computes the tranche delta (on current notional) by doing a proportional + ## blip of all the curves + ## if K2==1, then computes the delta using the lower attachment only + ## this makes sense for bottom-up skews + eps <- 1e-4 + portfolioplus <- portfolio + portfoliominus <- portfolio + cs <- couponSchedule(IMMDate(tradedate), index$maturity,"Q", "FIXED", coupon, 0, tradedate, + IMMDate(tradedate, "prev")) + for(i in 1:length(portfolio)){ + portfolioplus[[i]]@curve@hazardrates <- portfolioplus[[i]]@curve@hazardrates * (1 + eps) + portfoliominus[[i]]@curve@hazardrates <- portfoliominus[[i]]@curve@hazardrates * (1- eps) + } + dPVindex <- indexpv(portfolioplus, index, tradedate = tradedate, clean = FALSE)$bp - + indexpv(portfoliominus, index, tradedate = tradedate, clean = FALSE)$bp + defaultprobplus <- 1 - SPmatrix(portfolioplus, length(cs$dates)) + defaultprobminus <- 1 - SPmatrix(portfoliominus, length(cs$dates)) + if(K2==1){ + K1adj <- adjust.attachments(K1, index$loss, index$factor) + dPVtranche <- BCtranche.pv(defaultprobplus, issuerweights, recov, cs, 0, K1adj, 0, rho1, + Z, w, N)$bp - + BCtranche.pv(defaultprobminus, issuerweights, recov, cs, 0, K1adj, 0, rho1, + Z, w, N)$bp + delta <- (1 - dPVtranche/(100*dPVindex) * K1adj)/(1-K1adj) + }else{ + Kmod <- adjust.attachments(c(K1, K2), index$loss, index$factor) + dPVtranche <- BCtranche.pv(defaultprobplus, issuerweights, recov, cs, Kmod[1], Kmod[2], rho1, rho2, + Z, w, N)$bp - + BCtranche.pv(defaultprobminus, issuerweights, recov, cs, Kmod[1], Kmod[2], rho1, rho2, + Z, w, N)$bp + delta <- dPVtranche/(100*dPVindex) + } + ## dPVindex <- BCtranche.pv(portfolioplus, index, coupon, 0, 1, 0, 0.5, lu)$bp- + ## BCtranche.pv(portfoliominus, index, coupon, 0, 1, 0, 0.5, lu)$bp + return( delta ) +} + +BCstrikes <- function(portfolio, index, coupon, K, rho, N=101) { + ## computes the strikes as a percentage of expected loss + EL <- c() + for(i in 2:length(K)){ + EL <- c(EL, -BCtranche.pv(portfolio, index, coupon, K[i-1], K[i], rho[i-1], rho[i], N)$pl) + } + Kmodified <- adjust.attachments(K, index$loss, index$factor) + return(cumsum(EL*diff(Kmodified))/sum(EL*diff(Kmodified))) +} + +delta.factor <- function(K1, K2, index){ + ## compute the factor to convert from delta on current notional to delta on original notional + ## K1 and K2 original strikes + factor <- (adjust.attachments(K2, index$loss, index$factor) + -adjust.attachments(K1, index$loss, index$factor))/(K2-K1) + return( factor ) +} + +MFupdate.prob <- function(Z, w, rho, defaultprob){ + ## update the probabilites based on a non gaussian factor + ## distribution so that the pv of the cds stays the same. + p <- matrix(0, nrow(defaultprob), ncol(defaultprob)) + for(i in 1:nrow(defaultprob)){ + for(j in 1:ncol(defaultprob)){ + p[i,j] <- fit.prob(Z, w, rho[i], defaultprob[i,j]) + } + } + return( p ) +} + +MFupdate.probC <- function(Z, w, rho, defaultprob){ + ## update the probabilities based on a non gaussian factor + ## distribution so that the pv of the cds stays the same. + p <- matrix(0, nrow(defaultprob), ncol(defaultprob)) + for(i in 1:nrow(defaultprob)){ + for(j in 1:ncol(defaultprob)){ + p[i,j] <- .C("fitprob", as.double(Z), as.double(w), as.integer(length(Z)), + as.double(rho[i]), as.double(defaultprob[i,j]), q = double(1))$q + } + } + return( p ) +} + +MFlossrecovdist.prepay <- function(w, Z, rho, defaultprob, defaultprobmod, prepayprob, prepayprobmod, + issuerweights, recov, Ngrid=2*length(issuerweights)+1, defaultflag=FALSE){ + ## computes the loss and recovery distribution using the modified factor distribution + n.credit <- length(issuerweights) + n.int <- length(w) + Rstoch <- array(0, dim=c(n.int, n.credit, ncol(defaultprob))) + for(t in 1:ncol(defaultprob)){ + for(i in 1:n.credit){ + Rstoch[,i,t] <- stochasticrecov(recov[i], 0, Z, w, rho, defaultprob[i,t], defaultprobmod[i,t]) + } + } + parf <- function(i){ + dpshocked <- apply(defaultprobmod, 2, shockprob, rho=rho, Z=Z[i]) + ppshocked <- apply(prepayprobmod, 2, shockprob, rho=rho, Z=-Z[i]) + S <- 1 - Rstoch[i,,] + dist <- lossrecovdist.term(dpshocked, ppshocked, issuerweights, S, Ngrid, defaultflag) + } + L <- matrix(0, Ngrid, ncol(defaultprob)) + R <- matrix(0, Ngrid, ncol(defaultprob)) + for(i in 1:length(w)){ + dist <- parf(i) + L <- L + dist$L * w[i] + R <- R + dist$R * w[i] + } + return( list(L=L, R=R) ) +} + +MFlossdist.joint <- function(cl, w, Z, rho, defaultprob, defaultprobmod, issuerweights, recov, + Ngrid=2*length(issuerweights)+1, defaultflag=FALSE){ + ## rowSums(Q) is the default/loss distribution depending if + ## defaultflag is TRUE or FALSE (default setting is FALSE) + ## colSums(Q) is the recovery distribution + ## so that recovery is the y axis and L/D is the x axis + ## if we use the persp function, losses is the axes facing us, + ## and R is the axis going away from us. + n.credit <- length(issuerweights) + n.int <- lenth(w) + Rstoch <- array(0, dim=c(n.int, n.credit, ncol(defaultprob))) + for(t in 1:ncol(defaultprob)){ + for(i in 1:n.credit){ + Rstoch[,i,t] <- stochasticrecov(recov[i], 0, Z, w, rho, defaultprob[i,t], defaultprobmod[i,t]) + } + } + parf <- function(i){ + pshocked <- apply(defaultprobmod, 2, shockprob, rho=rho, Z=Z[i]) + S <- 1 - Rstoch[i,,] + dist <- lossrecovdist.joint.term(pshocked, 0, issuerweights, S, Ngrid, defaultflag) + gc() + return(dist) + } + temp <- parSapply(cl, 1:length(w), parf) + clusterCall(cl, gc) + Q <- array(0, dim=c(ncol(defaultprob), Ngrid, Ngrid)) + for(i in 1:length(w)){ + Q <- Q + w[i]*array(temp[,i], dim=c(ncol(defaultprob), Ngrid, Ngrid)) + } + return( Q ) +} + +MFlossdist.prepay.joint <- function(w, Z, rho, defaultprob, defaultprobmod, + prepayprob, prepayprobmod, issuerweights, recov, + Ngrid=2*length(issuerweights)+1, defaultflag=FALSE){ + ## rowSums is the loss distribution + ## colSums is the recovery distribution + ## so that recovery is the y axis and L is the x axis + ## if we use the persp function, losses is the axes facing us, + ## and R is the axis going away from us. + n.credit <- length(issuerweights) + n.int <- length(w) + Rstoch <- array(0, dim=c(n.credit, n.int, ncol(defaultprob))) + + for(t in 1:ncol(defaultprob)){ + for(i in 1:n.credit){ + Rstoch[i,,t] <- stochasticrecovC(recov[i], 0, Z, w, rho[i], + defaultprob[i,t], defaultprobmod[i,t]) + } + } + + Q <- array(0, dim=c(ncol(defaultprob), Ngrid, Ngrid)) + for(t in 1:ncol(defaultprob)){ + S <- 1 - Rstoch[,,t] + Q[t,,] <- lossdistC.prepay.jointZ(defaultprobmod[,t], prepayprobmod[,t], issuerweights, + S, Ngrid, defaultflag, rho, Z, w) + } + return( Q ) +} + +MFtranche.pv <- function(cl, cs, w, rho, defaultprob, defaultprobmod, issuerweights, recov, + Kmodified, Ngrid=length(issuerweights)+1){ + ## computes the tranches pv using the modified factor distribution + ## p is the modified probability so that + n.credit <- length(issuerweights) + Rstoch <- array(0, dim=c(n.int, n.credit, ncol(defaultprob))) + for(t in 1:ncol(defaultprob)){ + for(i in 1:n.credit){ + Rstoch[,i,t] <- stochasticrecov(recov[i], 0, Z, w, rho, defaultprob[i,t], defaultprobmod[i,t]) + } + } + parf <- function(i){ + pshocked <- apply(defaultprobmod, 2, shockprob, rho=rho, Z=Z[i]) + S <- 1 - Rstoch[i,,] + dist <- lossrecovdist.term(pshocked, 0, issuerweights, S, Ngrid) + return( tranche.pvvec(Kmodified, dist$L, dist$R, cs)) + } + clusterExport(cl, list("Rstoch", "p"), envir=environment()) + result <- parSapply(cl, 1:length(w), parf) + return( list(pv=100*(1+result%*%w), pv.w=result)) +} diff --git a/src/Makevars b/src/Makevars new file mode 100644 index 0000000..f6790ff --- /dev/null +++ b/src/Makevars @@ -0,0 +1,2 @@ +PKG_CFLAGS=$(SHLIB_OPENMP_CFLAGS) +PKG_LIBS=$(SHLIB_OPENMP_CFLAGS) diff --git a/src/Makevars.win b/src/Makevars.win new file mode 100644 index 0000000..6eb170c --- /dev/null +++ b/src/Makevars.win @@ -0,0 +1,2 @@ +PKG_CFLAGS=$(SHLIB_OPENMP_CFLAGS) +PKG_LIBS=$(SHLIB_OPENMP_CFLAGS) -LOpenBLAS/lib -lopenblas diff --git a/src/lossdistrib.c b/src/lossdistrib.c new file mode 100644 index 0000000..77ea335 --- /dev/null +++ b/src/lossdistrib.c @@ -0,0 +1,760 @@ +#include +#include +#include +#include +#include "lossdistrib.h" + +#define MIN(x, y) (((x) < (y)) ? (x) : (y)) + +extern int dgemv_(char* trans, int *m, int *n, double* alpha, double* A, int* lda, + double* x, int* incx, double* beta, double* y, int* incy); +extern double ddot_(int* n, double* dx, int* incx, double* dy, int* incy); +extern int dscal_(int* n, double* da, double* dx, int* incx); +extern int daxpy_(int* n, double* da, double* dx, int* incx, double* dy, int* incy); +extern void openblas_set_num_threads(int); + +void lossdistrib(double *p, int *np, double *w, double *S, int *N, int *defaultflag, + double *q) { + /* recursive algorithm with first order correction for computing + the loss distribution. + p vector of default probabilities + np length of p + w issuer weights + S vector of severities (should be same length as p) + N number of ticks in the grid + defaultflag if true compute the default distribution + q the loss distribution */ + + int i, j, d1, d2, M; + double lu, d, p1, p2, sum; + double *qtemp; + + lu = 1./(*N-1); + qtemp = calloc(*N, sizeof(double)); + q[0] = 1; + M = 1; + for(i=0; i<(*np); i++){ + d = (*defaultflag)? w[i]/lu : S[i] * w[i]/ lu; + d1 = floor(d); + d2 = ceil(d); + p1 = p[i] * (d2-d); + p2 = p[i] - p1; + memcpy(qtemp, q, MIN(M, *N) * sizeof(double)); + for(j=0; j < MIN(M, *N); j++){ + q[j] = (1-p[i]) * q[j]; + } + for(j=0; j < MIN(M, *N-d2); j++){ + q[d1+j] += p1 * qtemp[j]; + q[d2+j] += p2 * qtemp[j]; + }; + M+=d2; + } + + /* correction for weight loss */ + if(M > *N){ + sum = 0; + for(j=0; j *N){ + sum = 0; + for(j=0; j*N || My>*N){ + sum = 0; + for(m=0; m < MIN(Mx, *N); m++){ + for(n=0; n < MIN(My, *N); n++){ + sum += q[m+n*(*N)]; + } + } + q[MIN(*N, Mx)*MIN(My,*N)-1] += 1 - sum; + } + free(qtemp); +} + +void lossdistrib_joint_blas(double *p, int *np, double *w, double *S, int *N, int *defaultflag, double *q) { + /* recursive algorithm with first order correction + computes jointly the loss and recovery distribution + p vector of default probabilities + np length of p + w vector of issuer weights (length np) + S vector of severities (should be same length as p) + N number of ticks in the grid + defaultflag if true computes the default distribution + q the joint probability distribution */ + + int i, j, k, m, n; + int Mx, My; + double lu, x, y, sum, pbar; + double alpha1, alpha2, beta1, beta2; + double w1, w2, w3, w4; + double *qtemp; + int bound; + int one = 1; + + /* only use one thread, performance is horrible if use multiple threads */ + openblas_set_num_threads(1); + + lu = 1./(*N-1); + qtemp = calloc( (*N) * (*N), sizeof(double)); + q[0] = 1; + Mx=1; + My=1; + for(k=0; k<(*np); k++){ + x = (*defaultflag)? w[k] /lu : S[k] * w[k] / lu; + y = (1-S[k]) * w[k] / lu; + i = floor(x); + j = floor(y); + alpha1 = i + 1 - x; + alpha2 = 1 - alpha1; + beta1 = j + 1 - y; + beta2 = 1 - beta1; + w1 = alpha1 * beta1 * p[k]; + w2 = alpha1 * beta2 * p[k]; + w3 = alpha2 * beta2 * p[k]; + w4 = alpha2 * beta1 * p[k]; + + for(n=0; n*N || My>*N){ + sum = 0; + for(m=0; m < MIN(Mx, *N); m++){ + for(n=0; n < MIN(My, *N); n++){ + sum += q[m+n*(*N)]; + } + } + q[MIN(*N, Mx)*MIN(My,*N)-1] += 1 - sum; + } + free(qtemp); +} + +void recovdist(double *dp, double *pp, int *n, double *w, double *S, int *N, double *q) { + /* recursive algorithm with first order correction for computing + the recovery distribution in case of prepayment. + dp vector of default probabilities + pp vector of prepay probabilities + n length of p + S vector of severities (should be same length as p) + w vector of weights + N number of ticks in the grid + q the loss distribution */ + + int i, j, d1l, d1u, d2l, d2u; + int M; + double lu, d1, d2, dp1, dp2, pp1, pp2, sum; + double *qtemp; + + lu = 1./(*N - 1); + qtemp = calloc( (*N), sizeof(double)); + q[0] = 1; + M=1; + for(i=0; i<(*n); i++){ + d1 = w[i] * (1-S[i]) /lu; + d2 = w[i]/lu; + d1l = floor(d1); + d1u = d1l + 1; + d2l = floor(d2); + d2u = d2l + 1; + dp1 = dp[i] * (d1u - d1); + dp2 = dp[i] - dp1; + pp1 = pp[i] * (d2u - d2); + pp2 = pp[i] - pp1; + memcpy(qtemp, q, MIN(M, *N) * sizeof(double)); + for(j = 0; j< MIN(M, *N); j++){ + q[j] = (1-dp[i]-pp[i]) * q[j]; + } + for(j=0; j < MIN(M, *N-d2u); j++){ + q[d1l+j] += dp1 * qtemp[j]; + q[d1u+j] += dp2 * qtemp[j]; + q[d2l+j] += pp1 * qtemp[j]; + q[d2u+j] += pp2 * qtemp[j]; + }; + M += d2u; + } + /* correction for weight loss */ + if(M>*N){ + sum = 0; + for(j=0; j*N || My>*N){ + sum = 0; + for(m=0; m < MIN(Mx, *N); m++){ + for(n=0; n < MIN(My, *N); n++){ + sum += q[m+n*(*N)]; + } + } + q[MIN(*N, Mx)*MIN(My,*N)-1] += 1 - sum; + } + free(qtemp); +} + +void lossdistrib_prepay_joint_blas(double *dp, double *pp, int *ndp, double *w, + double *S, int *N, int *defaultflag, double *q) { + int i, j1, j2, k, m, n; + double lu, x, y1, y2, sum; + double alpha1, alpha2, beta1, beta2; + double dpw1, dpw2, dpw3, dpw4; + double ppw1, ppw2, ppw3; + double *qtemp; + int Mx, My, bound; + double pbar; + int one = 1; + + lu = 1./(*N-1); + qtemp = calloc((*N) * (*N), sizeof(double)); + q[0] = 1; + Mx=1; + My=1; + + /* only use one thread, performance is horrible if use multiple threads */ + openblas_set_num_threads(1); + for(k=0; k<(*ndp); k++){ + y1 = (1-S[k]) * w[k]/lu; + y2 = w[k]/lu; + x = (*defaultflag)? y2: y2-y1; + i = floor(x); + j1 = floor(y1); + j2 = floor(y2); + alpha1 = i + 1 - x; + alpha2 = 1 - alpha1; + beta1 = j1 + 1 - y1; + beta2 = 1 - beta1; + dpw1 = alpha1 * beta1 * dp[k]; + dpw2 = alpha1 * beta2 * dp[k]; + dpw3 = alpha2 * beta2 * dp[k]; + dpw4 = alpha2 * beta1 * dp[k]; + + /* by default distribution, we mean names fractions of names that disappeared + either because of default or prepayment */ + for(n=0; n*N || My>*N){ + sum = 0; + for(m=0; m < MIN(Mx, *N); m++){ + for(n=0; n < MIN(My, *N); n++){ + sum += q[m+n*(*N)]; + } + } + q[MIN(*N, Mx)*MIN(My,*N)-1] += 1 - sum; + } + free(qtemp); +} + +double shockprob(double p, double rho, double Z, int give_log){ + if(rho==1){ + return((double)(Z<=qnorm(p, 0, 1, 1, 0))); + }else{ + return( pnorm( (qnorm(p, 0, 1, 1, 0) - sqrt(rho) * Z)/sqrt(1 - rho), 0, 1, 1, give_log)); + } +} + +double dqnorm(double x){ + return 1/dnorm(qnorm(x, 0, 1, 1, 0), 0, 1, 0); +} + +double dshockprob(double p, double rho, double Z){ + return( dnorm((qnorm(p, 0, 1, 1, 0) - sqrt(rho) * Z)/sqrt(1-rho), 0, 1, 0) * dqnorm(p)/sqrt(1-rho) ); +} + +void shockprobvec2(double p, double rho, double* Z, int nZ, double *q){ + /* return a two column vectors with shockprob in the first column + and dshockprob in the second column*/ + int i; + #pragma omp parallel for + for(i = 0; i < nZ; i++){ + q[i] = shockprob(p, rho, Z[i], 0); + q[i + nZ] = dshockprob(p, rho, Z[i]); + } +} + +double shockseverity(double S, double Z, double rho, double p){ + if(p==0){ + return 0; + }else{ + return( exp(shockprob(S * p, rho, Z, 1) - shockprob(p, rho, Z, 1)) ); + } +} + +double quantile(double* Z, double* w, int nZ, double p0){ + double cumw; + int i; + cumw = 0; + for(i=0; i p0){ + break; + } + } + return( Z[i] ); +} + +void fitprob(double* Z, double* w, int* nZ, double* rho, double* p0, double* result){ + double eps = 1e-12; + int one = 1; + double *q = malloc( 2 * (*nZ) * sizeof(double)); + double dp, p, phi; + + if(*p0 == 0){ + *result = 0; + }else if(*rho == 1){ + *result = pnorm(quantile(Z, w, *nZ, *p0), 0, 1, 1, 0); + }else{ + shockprobvec2(*p0, *rho, Z, *nZ, q); + dp = (ddot_(nZ, q, &one, w, &one) - *p0)/ddot_(nZ, q + *nZ, &one, w, &one); + p = *p0; + while(fabs(dp) > eps){ + phi = 1; + while( (p - phi * dp) < 0 || (p - phi * dp) > 1){ + phi *= 0.8; + } + p -= phi * dp; + shockprobvec2(p, *rho, Z, *nZ, q); + dp = (ddot_(nZ, q, &one, w, &one) - *p0)/ddot_(nZ, q + *nZ, &one, w, &one); + } + *result = p; + } + free(q); +} + +void stochasticrecov(double* R, double* Rtilde, double* Z, double* w, int* nZ, + double* rho, double* porig, double* pmod, double* q){ + double ptemp, ptilde; + int i; + if(*porig==0){ + for(i = 0; i < *nZ; i++){ + q[i] = *R; + } + }else{ + ptemp = (1 - *R) / (1 - *Rtilde) * *porig; + fitprob(Z, w, nZ, rho, &ptemp, &ptilde); + #pragma omp parallel for + for(i = 0; i < *nZ; i++){ + q[i] = fabs(1 - (1 - *Rtilde) * exp( shockprob(ptilde, *rho, Z[i], 1) - + shockprob(*pmod, *rho, Z[i], 1))); + } + } +} + +void lossdistrib_prepay_joint_Z(double *dp, double *pp, int *ndp, double *w, + double *S, int *N, int *defaultflag, double *rho, + double *Z, double *wZ, int *nZ, double *q) { + int i, j; + double* dpshocked = malloc(sizeof(double) * (*ndp) * (*nZ)); + double* ppshocked = malloc(sizeof(double) * (*ndp) * (*nZ)); + int N2 = (*N) * (*N); + double* qmat = malloc(sizeof(double) * N2 * (*nZ)); + + double alpha = 1; + double beta = 0; + int one = 1; + +#pragma omp parallel for private(j) + for(i = 0; i < *nZ; i++){ + for(j = 0; j < *ndp; j++){ + dpshocked[j + (*ndp) * i] = shockprob(dp[j], rho[j], Z[i], 0); + ppshocked[j + (*ndp) * i] = shockprob(pp[j], rho[j], -Z[i], 0); + } + lossdistrib_prepay_joint_blas(dpshocked + (*ndp) * i, ppshocked + (*ndp) * i, ndp, + w, S + (*ndp) * i, N, defaultflag, qmat + N2 * i); + } + + dgemv_("n", &N2, nZ, &alpha, qmat, &N2, wZ, &one, &beta, q, &one); + + free(dpshocked); + free(ppshocked); + free(qmat); +} + +void lossdistrib_joint_Z(double *dp, int *ndp, double *w, + double *S, int *N, int *defaultflag, double *rho, + double *Z, double *wZ, int *nZ, double *q) { + int i, j; + double* dpshocked = malloc(sizeof(double) * (*ndp) * (*nZ)); + int N2 = (*N) * (*N); + double* qmat = malloc(sizeof(double) * N2 * (*nZ)); + + double alpha = 1; + double beta = 0; + int one = 1; + +#pragma omp parallel for private(j) + for(i = 0; i < *nZ; i++){ + for(j = 0; j < *ndp; j++){ + dpshocked[j + (*ndp) * i] = shockprob(dp[j], rho[j], Z[i], 0); + } + lossdistrib_joint_blas(dpshocked + (*ndp) * i, ndp, w, S + (*ndp) * i, N, + defaultflag, qmat + N2 * i); + } + + dgemv_("n", &N2, nZ, &alpha, qmat, &N2, wZ, &one, &beta, q, &one); + + free(dpshocked); + free(qmat); +} + +void BClossdist(double *defaultprob, int *dim1, int *dim2, + double *issuerweights, double *recov, double *Z, double *w, + int *n, double *rho, int *N, int *defaultflag, + double *L, double *R) { + /* + computes the loss and recovery distribution over time with a flat gaussian + correlation + inputs: + defaultprob: matrix of size dim1 x dim2. dim1 is the number of issuers + and dim2 number of time steps + issuerweights: vector of issuer weights (length dim1) + recov: vector of recoveries (length dim1) + Z: vector of factor values (length n) + w: vector of factor weights (length n) + rho: correlation beta vector (length dim1) + N: number of ticks in the grid + defaultflag: if true, computes the default distribution + outputs: + L: matrix of size (N, dim2) + R: matrix of size (N, dim2) + */ + int t, i, j; + double g; + double *gshocked, *Rshocked, *Sshocked, *Lw, *Rw; + int one = 1; + double alpha = 1; + double beta = 0; + + gshocked = malloc((*dim1) * (*n) * sizeof(double)); + Rshocked = malloc((*dim1) * (*n) * sizeof(double)); + Sshocked = malloc((*dim1) * (*n) * sizeof(double)); + Lw = malloc((*N) * (*n) * sizeof(double)); + Rw = malloc((*N) * (*n) * sizeof(double)); + + + for(t=0; t < (*dim2); t++) { + memset(Lw, 0, (*N) * (*n) * sizeof(double)); + memset(Rw, 0, (*N) * (*n) * sizeof(double)); + #pragma omp parallel for private(j, g) + for(i=0; i < *n; i++){ + for(j=0; j < (*dim1); j++){ + g = defaultprob[j + (*dim1) * t]; + gshocked[j+(*dim1)*i] = shockprob(g, rho[j], Z[i], 0); + Sshocked[j+(*dim1)*i] = shockseverity(1-recov[j], Z[i], rho[j], g); + Rshocked[j+(*dim1)*i] = 1 - Sshocked[j+(*dim1)*i]; + } + lossdistrib_blas(gshocked + (*dim1) * i, dim1, issuerweights, Sshocked + (*dim1)*i, N, defaultflag, + Lw + i * (*N)); + lossdistrib_blas(gshocked + (*dim1) * i, dim1, issuerweights, Rshocked + (*dim1)*i, N, defaultflag, + Rw + i * (*N)); + } + dgemv_("n", N, n, &alpha, Lw, N, w, &one, &beta, L + t * (*N), &one); + dgemv_("n", N, n, &alpha, Rw, N, w, &one, &beta, R + t * (*N), &one); + } + free(gshocked); + free(Rshocked); + free(Sshocked); + free(Lw); + free(Rw); +} diff --git a/src/lossdistrib.h b/src/lossdistrib.h new file mode 100644 index 0000000..d4bd974 --- /dev/null +++ b/src/lossdistrib.h @@ -0,0 +1,46 @@ +void lossdistrib(double *p, int *np, double *w, double *S, int *N, int *defaultflag, double *q); +void lossdistrib_blas(double *p, int *np, double *w, double *S, int *N, int *defaultflag, double *q); + +double shockprob(double p, double rho, double Z, int give_log); + +void lossdistrib_Z(double *p, int *np, double *w, double *S, int *N, int *defaultflag, + double *rho, double *Z, int *nZ, double *q); + +void lossdistrib_truncated(double *p, int *np, double *w, double *S, int *N, + int *T, int *defaultflag, double *q); + +void lossdistrib_joint(double *p, int *np, double *w, double *S, int *N, + int *defaultflag, double *q); + +void lossdistrib_joint_blas(double *p, int *np, double *w, double *S, int *N, + int *defaultflag, double *q); + +void recovdist(double *dp, double *pp, int *n, double *w, double *S, int *N, double *q); + +void lossdistrib_prepay_joint(double *dp, double *pp, int *ndp, double *w, + double *S, int *N, int *defaultflag, double *q); +double dqnorm(double x); + +double dshockprob(double p, double rho, double Z); + +void shockprobvec2(double p, double rho, double* Z, int nZ, double *q); + +double shockseverity(double S, double Z, double rho, double p); + +void fitprob(double* Z, double* w, int* nZ, double* rho, double* p0, double* result); + +void stochasticrecov(double* R, double* Rtilde, double* Z, double* w, int* nZ, + double* rho, double* porig, double* pmod, double* q); + +void lossdistrib_prepay_joint_Z(double *dp, double *pp, int *ndp, double *w, + double *S, int *N, int *defaultflag, double *rho, + double *Z, double *wZ, int *nZ, double *q); +void lossdistrib_joint_Z(double *dp, int *ndp, double *w, + double *S, int *N, int *defaultflag, double *rho, + double *Z, double *wZ, int *nZ, double *q); + +void BClossdist(double *SurvProb, int *dim1, int *dim2, double *issuerweights, + double *recov, double *Z, double *w, int *n, double *rho, int *N, + int *defaultflag, double *L, double *R); + +double quantile(double* Z, double* w, int nZ, double p0); -- cgit v1.2.3-70-g09d2