diff options
| -rw-r--r-- | R/lossdistrib.c | 15 | ||||
| -rw-r--r-- | R/tranche_functions.R | 32 |
2 files changed, 33 insertions, 14 deletions
diff --git a/R/lossdistrib.c b/R/lossdistrib.c index 3ad21069..0a2c30ab 100644 --- a/R/lossdistrib.c +++ b/R/lossdistrib.c @@ -14,6 +14,7 @@ extern void openblas_set_num_threads(int); void omp_set_num_threads(int);
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);
@@ -165,8 +166,8 @@ void lossdistrib_Z(double *p, int *np, double *w, double *S, int *N, int *defaul for(j = 0; j < *np; j++){
pshocked[j + (*np) * i] = shockprob(p[j], rho[j], Z[i], 0);
}
- lossdistrib(pshocked + (*np) * i, np, w, S + (*np) * i, N,
- defaultflag, q + (*N) * i);
+ lossdistrib_blas(pshocked + (*np) * i, np, w, S + (*np) * i, N,
+ defaultflag, q + (*N) * i);
}
free(pshocked);
}
@@ -617,7 +618,11 @@ void shockprobvec2(double p, double rho, double* Z, int nZ, double *q){ }
double shockseverity(double S, double Z, double rho, double p){
- return( exp(shockprob(S * p, rho, Z, 1) - shockprob(p, rho, Z, 1)) );
+ 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){
@@ -747,11 +752,11 @@ void BClossdist(double *defaultprob, int *dim1, int *dim2, 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 dim2)
+ 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
+ rho: correlation beta vector (length dim1)
N: number of ticks in the grid
defaultflag: if true, computes the default distribution
outputs:
diff --git a/R/tranche_functions.R b/R/tranche_functions.R index 2b7a6219..f22bdc68 100644 --- a/R/tranche_functions.R +++ b/R/tranche_functions.R @@ -450,17 +450,31 @@ dist.transform <- function(dist.joint){ }
shockprob <- function(p, rho, Z, log.p=F){
- ## computes the shocked default probability as a function of the copula factor
+ ## 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))
+ return(Z<=qnorm(p))
}else{
- pnorm((qnorm(p)-sqrt(rho)*Z)/sqrt(1-rho), log.p=log.p)
+ 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
- Stilde * exp( shockprob(S/Stilde*p, rho, Z, TRUE) - shockprob(p, rho, Z, TRUE))
+ ## 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){
@@ -657,21 +671,21 @@ tranche.pvvec <- function(K, L, R, cs){ }
BClossdist <- function(defaultprob, issuerweights, recov, rho, Z, w,
- N=length(recov)+1, defaultflag=FALSE){
+ 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
}
- LZ <- matrix(0, N, n.int)
- RZ <- matrix(0, N, n.int)
+ 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, 0, issuerweights, S.shocked, N)
+ temp <- lossrecovdist(g.shocked, , issuerweights, S.shocked, N)
LZ[,i] <- temp$L
RZ[,i] <- temp$R
}
|