From 1e2aea18da610f7f225a744c77081f6ced9379fe Mon Sep 17 00:00:00 2001 From: Guillaume Horel Date: Fri, 3 Oct 2014 15:01:19 -0400 Subject: big cleanup --- R/tranche_functions.R | 588 ++------------------------------------------------ 1 file changed, 13 insertions(+), 575 deletions(-) (limited to 'R/tranche_functions.R') diff --git a/R/tranche_functions.R b/R/tranche_functions.R index 97fbc83..8b10394 100644 --- a/R/tranche_functions.R +++ b/R/tranche_functions.R @@ -355,22 +355,14 @@ lossdistC.prepay.jointZ <- function(dp, pp, w, S, N, defaultflag = FALSE, rho, Z } lossrecovdist <- function(defaultprob, prepayprob, w, S, N, defaultflag=FALSE, useC=TRUE){ + lossdistrib2 <- if(useC) lossdistC + recovdist <- if(useC) recovdistC 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) - } + 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) - } + L <- lossdistrib2(defaultprob+defaultflag*prepayprob, w, S, N, defaultflag) + R <- recovdist(defaultprob, prepayprob, w, S, N) } return(list(L=L, R=R)) } @@ -398,28 +390,17 @@ lossrecovdist.term <- function(defaultprob, prepayprob, w, S, N, defaultflag=FAL 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.jointblas(defaultprob[,t], w, S[,t], N, defaultflag) - } - }else{ - for(t in 1:ncol(defaultprob)){ - Q[t,,] <- lossdistC.prepay.jointblas(defaultprob[,t], prepayprob[,t], w, S[,t], N, defaultflag) - } + 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{ - 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) - } + for(t in 1:ncol(defaultprob)){ + Q[t,,] <- lossdist.prepay.jointblas(defaultprob[,t], prepayprob[,t], w, S[,t], N, defaultflag) } } - gc() return(Q) } @@ -530,123 +511,6 @@ stochasticrecovC <- function(R, Rtilde, Z, w, rho, porig, pmod){ 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) ) -} - BClossdist <- function(defaultprob, issuerweights, recov, rho, Z, w, N=length(recov)+1, defaultflag=FALSE, n.int=500){ if(missing(Z)){ @@ -684,429 +548,3 @@ BClossdistC <- function(defaultprob, issuerweights, recov, rho, Z, 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.legs <- function(index, K, rho, complement=FALSE){ - ## computes the protection leg and couponleg of a 0-K tranche - ## if complement==TRUE, computes the protection leg and coupon leg of a K-1 tranche - if((K==0 && !complement) || (K==1 && complement)){ - return(list(cl=0, pl=0)) - }else if((K==1 && !complement) || (K==0 && complement)){ - return(BCindex.pv(index)) - }else{ - dist <- BClossdistC(index$defaultprob, index$issuerweights, index$recov, rho, index$Z, index$w, index$N) - if(complement){ - return(list(cl=tranche.cl(dist$L, dist$R, index$cs, K, 1), - pl=tranche.pl(dist$L, index$cs, K, 1))) - }else{ - return(list(cl=tranche.cl(dist$L, dist$R, index$cs, 0, K), - pl=tranche.pl(dist$L, index$cs, 0, K))) - } - } -} - -BCtranche.pv <- function(index, protection=FALSE, complement=FALSE){ - ## computes the protection leg, couponleg, and bond price of a tranche - ## in the base correlation setting - ## if complement=FALSE compute the pvs starting from 0 (bottom-up skew) - ## if complement=TRUE compute the pvs starting from 1 (top-down skew) - pl <- rep(0, length(index$rho)) - cl <- rep(0, length(index$rho)) - for(i in seq_along(index$rho)){ - temp <- BCtranche.legs(index, index$K[i], index$rho[i], complement) - pl[i] <- temp$pl - cl[i] <- temp$cl - } - dK <- diff(index$K) - plvec <- diff(pl)/dK - clvec <- diff(cl)/dK*index$tranches$running - if(complement){ - plvec <- -plvec - clvec <- -clvec - } - if(protection){ - bp <- -plvec-clvec - }else{ - bp <- 1+plvec+clvec - } - return(list(pl=plvec, cl=clvec, bp=bp)) -} - -MFtranche.pv <- function(index, dist, protection=FALSE){ - n.tranches <- length(index$K)-1 - pl <- rep(0, n.tranches) - cl <- rep(0, n.tranches) - for(i in seq.int(n.tranches)){ - pl[i] <- tranche.pl(dist$L, index$cs, index$K[i], index$K[i+1]) - cl[i] <- tranche.cl(dist$L, dist$R, index$cs, index$K[i], index$K[i+1]) - } - dK <- diff(index$K) - plvec <- pl/dK - clvec <- cl/dK*index$tranches$running - if(protection){ - bp <- -plvec-clvec - }else{ - bp <- 1+plvec+clvec - } - return(list(pl=plvec, cl=clvec, bp=bp)) -} - -adjust.skew <- function(index1, index2, method="ATM"){ - #index1 is the index for which we already have computed the skew - #index2 is the index we're mapping to - # if method="ATM", do simple at the money mapping - # method="TLP", do tranche loss proportion mapping - # method="PM", do probability matching - - K1 <- index1$K[-c(1,length(index1$K))] - K2 <- index2$K[-c(1,length(index2$K))] - aux <- function(x, index1, el1, skew, index2, el2, K2){ - return(abs(ELtrunc(index1, x, skew(x))/el1- - ELtrunc(index2, K2, skew(x))/el2)) - } - aux2 <- function(x, index1, skew, index2, K2){ - return(abs(log(Probtrunc(index1, x, skew(x)))- - log(Probtrunc(index2, K2, skew(x))))) - } - if(method %in% c("ATM", "TLP")){ - el1 <- EL(index1) - el2 <- EL(index2) - } - skew <- function(x){ - #we cap the correlation at 0.99 and 0.01 - f <- splinefun(K1, index1$rho[-c(1, length(index1$rho))], "natural") - return(pmax(pmin(f(x), 0.99), 0.01)) - } - if(method=="ATM"){ - K1eq <- el1/el2 * K2 - }else if(method == "TLP"){ - K1eq <- c() - m <- max(K2) + 0.3 - for(K2val in K2){ - prog <- optimize(aux, interval=c(0,m), - index1=index1, el1=el1, skew=skew, - index2=index2, el2=el2, K2=K2val) - K1eq <- c(K1eq, prog$minimum) - } - }else if (method=="PM"){ - K1eq <- c() - m <- max(K2) + 0.25 - for(K2val in K2){ - prog <- optimize(aux2, interval=c(0, m), - index1=index1, skew=skew, index2=index2, K2=K2val) - K1eq <- c(K1eq, prog$minimum) - } - } - return(c(NA, skew(K1eq), NA)) -} - -theta.adjust.skew <- function(index, shortened=4, method="ATM"){ - #ajust the correlation skew by doing ATM mapping on the expected loss - indexshort <- index - N <- nrow(index$cs)-shortened - indexshort$defaultprob <- indexshort$defaultprob[,1:N] - indexshort$cs <- indexshort$cs[1:N,] - return(adjust.skew(index, indexshort, method)) -} - -BCtranche.theta <- function(index, shortened=4, complement=FALSE, method="ATM"){ - temp <- BCtranche.pv(index, complement=complement) - rho.adj <- theta.adjust.skew(index, shortened, method) - if(any(rho.adj[-c(1, length(rho.adj))]<=0)){ - print("probable inverted skew: no adjustment") - }else{ - index$rho <- rho.adj - } - N <- nrow(index$cs) - shortened - index$cs <- index$cs[1:N,] - index$defaultprob <- index$defaultprob[,1:N] - temp2 <- BCtranche.pv(index, complement=complement) - temp3 <- BCtranche.delta(index, complement=complement) - return(data.frame(theta=temp2$bp-temp$bp+index$tranches$running, - fw.delta=temp3$delta)) -} - -BCtranche.delta <- function(index, complement=FALSE){ - ## computes the tranche delta (on current notional) by doing a proportional - ## blip of all the curves - ## if complement is False, then computes deltas bottom-up - ## if complement is True, then computes deltas top-down - eps <- 1e-4 - index$N <- 301 ## for gamma computations we need all the precision we can get - ## we build a lit of 4 indices with various shocks - index.list <- lapply(c(0, eps, -eps, 2*eps), function(x){ - if(x==0){ - return(index) - }else{ - newindex <- index - newindex$portfolio <- tweakportfolio(newindex$portfolio, x) - newindex$defaultprob <- 1 - SPmatrix(newindex$portfolio, length(index$cs$dates)) - return(newindex) - } - }) - bp <- matrix(0, length(index$K)-1, length(index.list)) - indexbp <- rep(0, length(index.list)) - for(j in seq_along(index.list)){ - indexbp[j] <- BCindex.pv(index.list[[j]])$bp - bp[,j] <- BCtranche.pv(index.list[[j]], complement=complement)$bp - } - - deltas <- (bp[,2]-bp[,3])/(indexbp[2]-indexbp[3])*tranche.factor(index)/index$factor - deltasplus <- (bp[,4]-bp[,1])/(indexbp[4]-indexbp[1])*tranche.factor(index)/index$factor - gammas <- (deltasplus-deltas)/(indexbp[2]-indexbp[1])/100 - - return( data.frame(delta=deltas, gamma=gammas) ) -} - -MFtranche.delta <- function(index){ - ## computes the tranche delta (on current notional) by doing a proportional - ## blip of all the curves - ## if complement is False, then computes deltas bottom-up - ## if complement is True, then computes deltas top-down - eps <- 1e-4 - index$Ngrid <- 301 ## for gamma computations we need all the precision we can get - ## we build a lit of 4 indices with various shocks - index.list <- lapply(c(0, eps, -eps, 2*eps), function(x){ - if(x==0){ - return(index) - }else{ - newindex <- index - newindex$portfolio <- tweakportfolio(newindex$portfolio, x) - newindex$defaultprob <- 1 - SPmatrix(newindex$portfolio, length(index$cs$dates)) - return(newindex) - } - }) - bp <- matrix(0, length(index$K)-1, length(index.list)) - indexbp <- rep(0, length(index.list)) - for(j in seq_along(index.list)){ - indexbp[j] <- BCindex.pv(index.list[[j]])$bp - dist <- MFdist(index.list[[j]]) - bp[,j] <- BCtranche.pv(index.list[[j]], dist) - } - - deltas <- (bp[,2]-bp[,3])/(indexbp[2]-indexbp[3])*tranche.factor(index)/index$factor - deltasplus <- (bp[,4]-bp[,1])/(indexbp[4]-indexbp[1])*tranche.factor(index)/index$factor - gammas <- (deltasplus-deltas)/(indexbp[2]-indexbp[1])/100 - - return( list(deltas=deltas, gammas=gammas) ) -} - -BCtranche.corr01 <- function(index, eps=0.01, complement=FALSE){ - ##does a parallel shift of the skew and computes the change in pv - before <- BCtranche.pv(index, complement=complement) - index$rho[-1] <- index$rho[-1]+eps - after <- BCtranche.pv(index, complement=complement) - return(after$bp-before$bp) -} - -EL <- function(index, discounted=TRUE, shortened=0){ - ## computes the expected loss of a portfolio (time discounted if discounted is TRUE) - ## given the default curves and recovery - ## should be very close to the protection leg of the portfolio of cds - ## index should be a list with issuerweights, recov, defaultprob and cs parameters - ## shortened: number of quarters to shorten the maturity by - Ncol <- ncol(index$defaultprob)-shortened - ELvec <- as.numeric(crossprod(index$issuerweights * (1-index$recov), index$defaultprob[,1:Ncol])) - if(!discounted){ - return( ELvec[length(ELvec)] ) - }else{ - return( sum(index$cs$df[1:Ncol]*diff(c(0, ELvec))) ) - } -} - -BCindex.pv <- function(index, discounted=TRUE, shortened=0){ - Ncol <- ncol(index$defaultprob)-shortened - ELvec <- as.numeric(crossprod(index$issuerweights * (1-index$recov), index$defaultprob[,1:Ncol])) - size <- 1-as.numeric(crossprod(index$issuerweights, index$defaultprob[,1:Ncol])) - sizeadj <- 0.5*(c(1, size[-length(size)])+size) - if(!discounted){ - pl <- -ELvec[length(ELvec)] - cl <- as.numeric(crossprod(index$cs$coupons[1:Ncol], sizeadj)) - bp <- 1+cl+pl - }else{ - pl <- -sum(index$cs$df[1:Ncol]* diff(c(0, ELvec))) - cl <- as.numeric(crossprod(index$cs$coupons[1:Ncol], sizeadj * index$cs$df[1:Ncol] )) - bp <- 1+cl+pl - } - bp <- 1+cl*index$quotes$spread+pl - return(list(pl=pl, cl=cl, bp=bp)) -} - -ELtrunc <- function(index, K, rho){ - ## computes the expected loss of a portfolio below strike K - ## could be written faster by using a truncated version of lossdist - ## index should be a list with issuerweights, recov, defaultprob and cs parameters - dist <- BClossdistC(index$defaultprob, index$issuerweights, index$recov, rho, index$Z, index$w, index$N) - return( -tranche.pl(dist$L, index$cs, 0, K)) -} - -Probtrunc <- function(index, K, rho){ - dist <- BClossdistC(index$defaultprob, index$issuerweights, index$recov, rho, index$Z, index$w, index$N) - p <- cumsum(dist$L[,ncol(dist$L)]) - support <- seq(0, 1, length=index$N) - probfun <- splinefun(support, p, method="hyman") - return(probfun(K)) -} - - -BCstrikes <- function(index, K, rho) { - ## computes the strikes as a percentage of expected loss - ## Kmodified is the current attachment points (adjusted for losses) - el <- EL(index) - ELvec <- c() - for(i in 2:length(K)){ - ELvec <- c(ELvec, ELtrunc(index, K[i], rho[i])) - } - return( ELvec/el ) -} - -tranche.factor <- function(index){ - ## compute the factor to convert from delta on current notional to delta on original notional - ## K1 and K2 original strikes - return( diff(index$K)/diff(index$K.orig)*index$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 ) -} - -MFrecovery <- function(index, defaultprobmod){ - n.credit <- length(index$issuerweights) - n.int <- length(index$Z) - Rstoch <- array(0, dim=c(n.credit, n.int, ncol(index$defaultprob))) - rho <- rep(0.45, n.credit) - for(t in 1:ncol(index$defaultprob)){ - for(i in 1:n.credit){ - Rstoch[i,,t] <- stochasticrecovC(index$recov[i], 0, index$Z, index$w.mod, - rho[i], index$defaultprob[i,t], defaultprobmod[i,t]) - } - } - return( Rstoch ) -} - -MFlossdist <- function(index){ - n.credit <- length(index$issuerweights) - rho <- rep(0.45, n.credit) - defaultprobmod <- MFupdate.probC(index$Z, index$w.mod, rho, index$defaultprob) - n.credit <- length(index$issuerweights) - n.int <- length(index$Z) - Rstoch <- MFrecovery(index, defaultprobmod) - Lw <- matrix(0, index$N, n.int) - Rw <- matrix(0, index$N, n.int) - L <- matrix(0, index$N, ncol(index$defaultprob)) - R <- matrix(0, index$N, ncol(index$defaultprob)) - for(t in 1:ncol(index$defaultprob)){ - S <- 1 - Rstoch[,,t] - Lw <- lossdistCZ(defaultprobmod[,t], index$issuerweights, S, index$N, 0, rho, index$Z) - Rw <- lossdistCZ(defaultprobmod[,t], index$issuerweights, 1-S, Ngrid, 0, rho, index$Z) - L[,t] <- Lw%*%index$w.mod - R[,t] <- Rw%*%index$w.mod - } - return(list(L=L, R=R)) -} -- cgit v1.2.3-70-g09d2