diff options
| -rw-r--r-- | R/distrib.R | 92 |
1 files changed, 46 insertions, 46 deletions
diff --git a/R/distrib.R b/R/distrib.R index ae5e51f..7dfafdd 100644 --- a/R/distrib.R +++ b/R/distrib.R @@ -49,9 +49,9 @@ lossdistrib <- function(p){ 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] + 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) } @@ -86,9 +86,9 @@ convolve <- function(dist1, dist2) { #' @export lossdistrib.fft <- function(p) { ## haven't tested when p is not a power of 2. - if(length(p) == 1){ - c(1-p, p) - }else { + if(length(p) == 1) { + c(1 - p, p) + } else { convolve(lossdistrib.fft(p[1:(length(p)/2)]), lossdistrib.fft(p[-(1:(length(p)/2))])) } @@ -121,7 +121,7 @@ lossdistrib2 <- function(p, w, S, N, defaultflag = FALSE){ for(i in 1:n){ if(defaultflag){ d <- w[i] /lu - }else{ + } else { d <- S[i] * w[i] / lu } d1 <- floor(d) @@ -238,10 +238,10 @@ lossdist.joint <- function(p, w, S, N, defaultflag = FALSE){ lu <- 1/(N-1) q <- matrix(0, N, N) q[1,1] <- 1 - for(k in 1:n){ - if(defaultflag){ + for(k in 1:n) { + if(defaultflag) { x <- w[k] / lu - }else{ + } else { x <- S[k] * w[k]/lu } y <- (1-S[k]) * w[k]/lu @@ -299,7 +299,7 @@ lossdist.prepay.joint <- function(dp, pp, w, S, N, defaultflag = FALSE){ 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{ + } else { ppsplit <- pp[k] * c(j2+1-y2, y2-j2) } qtemp <- matrix(0, N, N) @@ -312,7 +312,7 @@ lossdist.prepay.joint <- function(dp, pp, w, S, N, defaultflag = FALSE){ 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{ + } 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)] } @@ -410,7 +410,7 @@ lossrecovdist <- function(defaultprob, prepayprob, w, S, N, defaultflag=FALSE, u if(missing(prepayprob)){ L <- lossdistrib2(defaultprob, w, S, N, defaultflag) R <- lossdistrib2(defaultprob, w, 1-S, N) - }else{ + } else { L <- lossdistrib2(defaultprob+defaultflag*prepayprob, w, S, N, defaultflag) R <- recovdist(defaultprob, prepayprob, w, S, N) } @@ -422,14 +422,14 @@ lossrecovdist.term <- function(defaultprob, prepayprob, w, S, N, defaultflag=FAL ## 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)){ + 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)){ + } 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 @@ -462,11 +462,11 @@ dist.transform <- function(dist.joint) { ## 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(t in 1:dim(dist.joint)[1]) { for(i in 1:Ngrid){ for(j in 1:Ngrid){ index <- i+j - if(index <= Ngrid){ + 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] + @@ -485,13 +485,13 @@ shockprob <- function(p, rho, Z, log.p=F) { ## 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{ + 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{ + } 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) @@ -504,8 +504,8 @@ 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)) + 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) } @@ -515,7 +515,7 @@ dshockprob <- function(p,rho,Z) { } #' @export -dqnorm <- function(x){ +dqnorm <- function(x) { 1/dnorm(qnorm(x)) } @@ -523,23 +523,23 @@ dqnorm <- function(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 - if(p0==0){ + if(p0 == 0) { return(0) } - if(rho == 1){ + if(rho == 1) { distw <- distr::DiscreteDistribution(Z, w) return(distr::pnorm(distr::q(distw)(p0))) } eps <- 1e-12 - dp <- (crossprod(shockprob(p0,rho,Z),w)-p0)/crossprod(dshockprob(p0,rho,Z),w) + 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) + 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 + while ((p - phi * dp) < 0 || (p - phi * dp) > 1) { + phi <- 0.8 * phi } - p <- p - phi*dp + p <- p - phi * dp } return(p) } @@ -556,16 +556,16 @@ fit.probC <- function(Z, w, rho, p0) { 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){ + if(porig == 0) { return(rep(R, length(Z))) - }else{ + } 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)))) } } #' @export -stochasticrecovC <- function(R, Rtilde, Z, w, rho, porig, pmod){ +stochasticrecovC <- function(R, Rtilde, Z, w, rho, porig, pmod) { r <- .C(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))) @@ -574,14 +574,14 @@ stochasticrecovC <- function(R, Rtilde, Z, w, rho, porig, pmod){ #' @export BClossdist <- function(defaultprob, issuerweights, recov, rho, Z, w, - N=length(recov)+1, defaultflag=FALSE, n.int=500){ + N = length(recov)+1, defaultflag = FALSE, n.int = 500) { - if(missing(Z)){ + if(missing(Z)) { quadrature <- GHquad(n.int) Z <- quadrature$Z w <- quadrature$w } - stopifnot(length(Z)==length(w), + stopifnot(length(Z) == length(w), nrow(defaultprob) == length(issuerweights)) ## do not use if weights are not gaussian, results would be incorrect @@ -606,11 +606,11 @@ BClossdist <- function(defaultprob, issuerweights, recov, rho, Z, w, #' @export BClossdistC <- function(defaultprob, issuerweights, recov, rho, Z, w, - N=length(issuerweights)+1, defaultflag=FALSE){ - if(is.null(dim(defaultprob))){ + N = length(issuerweights)+1, defaultflag = FALSE) { + if(is.null(dim(defaultprob))) { dim(defaultprob) <- c(length(defaultprob),1) } - stopifnot(length(Z)==length(w), + stopifnot(length(Z) == length(w), nrow(defaultprob)==length(issuerweights), nrow(defaultprob)==length(recov)) L <- matrix(0, N, ncol(defaultprob)) @@ -625,9 +625,9 @@ BClossdistC <- function(defaultprob, issuerweights, recov, rho, Z, w, #' @export BCER <- function(defaultprob, issuerweights, recov, K, rho, Z, w, - N=length(issuerweights)+1, defaultflag=FALSE){ - stopifnot(length(Z)==length(w), - nrow(defaultprob)==length(issuerweights)) + N = length(issuerweights)+1, defaultflag = FALSE) { + stopifnot(length(Z) == length(w), + nrow(defaultprob) == length(issuerweights)) rho <- rep(rho, length(issuerweights)) ELt <- numeric(ncol(defaultprob)) |