From f11771ed1191a9db78007aa3d7bca002bab0e0f6 Mon Sep 17 00:00:00 2001 From: Guillaume Horel Date: Wed, 8 Oct 2014 10:39:29 -0400 Subject: rename file since there are no tranche functions left --- R/distrib.R | 550 ++++++++++++++++++++++++++++++++++++++++++++++++++ R/tranche_functions.R | 550 -------------------------------------------------- 2 files changed, 550 insertions(+), 550 deletions(-) create mode 100644 R/distrib.R delete mode 100644 R/tranche_functions.R diff --git a/R/distrib.R b/R/distrib.R new file mode 100644 index 0000000..8b10394 --- /dev/null +++ b/R/distrib.R @@ -0,0 +1,550 @@ +## 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 + +#' Gauss-Hermite quadrature weights +#' +#' \code{GHquad} computes the quadrature weights for integrating against a +#' Gaussian distribution. +#' +#' if f is a function, then with(GHquad(100), crossprod(f(Z), w)) +#' +#' @param n Integer, the number of nodes +#' @return A list with two components: +#' \item{Z}{the list of nodes} +#' \item{w}{the corresponding weights} +#' +GHquad <- function(n){ + n <- as.integer(n) + Z <- double(n) + w <- double(n) + result <- .C("GHquad", n, Z=Z, w=w) + result[[1]] <- NULL + return(result) +} + +#' Loss distribution of a portfolio +#' +#' \code{lossdistrib} computes the probability distribution of a sum +#' of independent Bernouilli variables with unequal probabilities. +#' +#' This uses the basic recursive algorithm of Andersen, Sidenius and Basu +#' We compute the probability distribution of S = \sum_{i=1}^n X_i +#' where X_i is Bernouilli(p_i) +#' @param p Numeric vector, the vector of success probabilities +#' @return A vector q such that q[k]=P(S=k) +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) +} + +#' Loss distribution of a portfolio +#' +#' \code{lossdistrib.fft} computes the probability distribution of a sum +#' of independent Bernouilli variables with unequal probabilities. +#' +#' This uses 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. +#' It is slower than the recursive algorithm in practice. +#' We compute the probability distribution of S = \sum_{i=1}^n X_i +#' where X_i is Bernouilli(p_i) +#' @param p Numeric vector, the vector of success probabilities +#' @return A vector such that q[k]=P(S=k) +lossdistrib.fft <- function(p){ + n <- length(p) + 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))) +} + +#' recursive algorithm with first order correction +#' +#' @param p Numeric, vector of default probabilities +#' @param w Numeric, vector of weights +#' @param S Numeric, vector of severities +#' @param N Integer, number of ticks in the grid +#' @param defaultflag Boolean, if True, we compute the default distribution +#' (instead of the loss distribution). +#' @return a Numeric vector of size \code{N} computing the loss (resp. +#' default) distribution if \code{defaultflag} is FALSE (resp. TRUE). +lossdistrib2 <- function(p, w, S, N, defaultflag=FALSE){ + 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 w[i] X_i + ## where X_i = 0 w.p 1 - dp[i] - pp[i] + ## = 1 - S[i] w.p dp[i] + ## = 1 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){ + #S is of size (length(p), length(Z)) + .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.jointblas <- function(dp, pp, w, S, N, defaultflag=FALSE){ + ## C version of lossdist.prepay.joint + r <- .C("lossdistrib_prepay_joint_blas", 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 vector of doubles of size n.credits + ## 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)), output = matrix(0,N,N)) + return(r$output) +} + +lossrecovdist <- function(defaultprob, prepayprob, w, S, N, defaultflag=FALSE, useC=TRUE){ + lossdistrib2 <- if(useC) lossdistC + recovdist <- if(useC) recovdistC + if(missing(prepayprob)){ + L <- lossdistrib2(defaultprob, w, S, N, defaultflag) + R <- lossdistrib2(defaultprob, w, 1-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)) + lossdist.joint <- if(useC) lossdistC.jointblas + lossdist.prepay.joint <- if(useC) lossdistC.prepay.jointblas + 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.jointblas(defaultprob[,t], prepayprob[,t], w, S[,t], N, defaultflag) + } + } + 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) +} + +BClossdist <- function(defaultprob, issuerweights, recov, rho, Z, w, + N=length(recov)+1, defaultflag=FALSE, n.int=500){ + if(missing(Z)){ + quadrature <- GHquad(n.int) + Z <- quadrature$Z + w <- quadrature$w + } + ## 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("BCloss_recov_dist", 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)) +} diff --git a/R/tranche_functions.R b/R/tranche_functions.R deleted file mode 100644 index 8b10394..0000000 --- a/R/tranche_functions.R +++ /dev/null @@ -1,550 +0,0 @@ -## 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 - -#' Gauss-Hermite quadrature weights -#' -#' \code{GHquad} computes the quadrature weights for integrating against a -#' Gaussian distribution. -#' -#' if f is a function, then with(GHquad(100), crossprod(f(Z), w)) -#' -#' @param n Integer, the number of nodes -#' @return A list with two components: -#' \item{Z}{the list of nodes} -#' \item{w}{the corresponding weights} -#' -GHquad <- function(n){ - n <- as.integer(n) - Z <- double(n) - w <- double(n) - result <- .C("GHquad", n, Z=Z, w=w) - result[[1]] <- NULL - return(result) -} - -#' Loss distribution of a portfolio -#' -#' \code{lossdistrib} computes the probability distribution of a sum -#' of independent Bernouilli variables with unequal probabilities. -#' -#' This uses the basic recursive algorithm of Andersen, Sidenius and Basu -#' We compute the probability distribution of S = \sum_{i=1}^n X_i -#' where X_i is Bernouilli(p_i) -#' @param p Numeric vector, the vector of success probabilities -#' @return A vector q such that q[k]=P(S=k) -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) -} - -#' Loss distribution of a portfolio -#' -#' \code{lossdistrib.fft} computes the probability distribution of a sum -#' of independent Bernouilli variables with unequal probabilities. -#' -#' This uses 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. -#' It is slower than the recursive algorithm in practice. -#' We compute the probability distribution of S = \sum_{i=1}^n X_i -#' where X_i is Bernouilli(p_i) -#' @param p Numeric vector, the vector of success probabilities -#' @return A vector such that q[k]=P(S=k) -lossdistrib.fft <- function(p){ - n <- length(p) - 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))) -} - -#' recursive algorithm with first order correction -#' -#' @param p Numeric, vector of default probabilities -#' @param w Numeric, vector of weights -#' @param S Numeric, vector of severities -#' @param N Integer, number of ticks in the grid -#' @param defaultflag Boolean, if True, we compute the default distribution -#' (instead of the loss distribution). -#' @return a Numeric vector of size \code{N} computing the loss (resp. -#' default) distribution if \code{defaultflag} is FALSE (resp. TRUE). -lossdistrib2 <- function(p, w, S, N, defaultflag=FALSE){ - 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 w[i] X_i - ## where X_i = 0 w.p 1 - dp[i] - pp[i] - ## = 1 - S[i] w.p dp[i] - ## = 1 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){ - #S is of size (length(p), length(Z)) - .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.jointblas <- function(dp, pp, w, S, N, defaultflag=FALSE){ - ## C version of lossdist.prepay.joint - r <- .C("lossdistrib_prepay_joint_blas", 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 vector of doubles of size n.credits - ## 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)), output = matrix(0,N,N)) - return(r$output) -} - -lossrecovdist <- function(defaultprob, prepayprob, w, S, N, defaultflag=FALSE, useC=TRUE){ - lossdistrib2 <- if(useC) lossdistC - recovdist <- if(useC) recovdistC - if(missing(prepayprob)){ - L <- lossdistrib2(defaultprob, w, S, N, defaultflag) - R <- lossdistrib2(defaultprob, w, 1-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)) - lossdist.joint <- if(useC) lossdistC.jointblas - lossdist.prepay.joint <- if(useC) lossdistC.prepay.jointblas - 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.jointblas(defaultprob[,t], prepayprob[,t], w, S[,t], N, defaultflag) - } - } - 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) -} - -BClossdist <- function(defaultprob, issuerweights, recov, rho, Z, w, - N=length(recov)+1, defaultflag=FALSE, n.int=500){ - if(missing(Z)){ - quadrature <- GHquad(n.int) - Z <- quadrature$Z - w <- quadrature$w - } - ## 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("BCloss_recov_dist", 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)) -} -- cgit v1.2.3-70-g09d2