aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--R/lossdistrib.c237
-rw-r--r--R/tranche_functions.R19
2 files changed, 243 insertions, 13 deletions
diff --git a/R/lossdistrib.c b/R/lossdistrib.c
index 22a449e2..36be7be3 100644
--- a/R/lossdistrib.c
+++ b/R/lossdistrib.c
@@ -8,8 +8,9 @@
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);
@@ -22,7 +23,10 @@ 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);
+ 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);
@@ -52,7 +56,8 @@ 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);
-void lossdistrib(double *p, int *np, double *w, double *S, int *N, int *defaultflag, double *q) {
+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
@@ -98,6 +103,55 @@ void lossdistrib(double *p, int *np, double *w, double *S, int *N, int *defaultf
Free(qtemp);
}
+void lossdistrib_blas(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
+ S vector of severities (should be same length as p)
+ N number of ticks in the grid
+ defaultflat if true compute the default distribution
+ q the loss distribution */
+
+ int i, j, d1, d2, M;
+ double lu, d, p1, p2, sum;
+ double *qtemp;
+ int bound;
+ double pbar;
+ int one = 1;
+
+ lu = 1./(*N-1);
+ qtemp = Calloc(*N, 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));
+ pbar = 1-p[i];
+ bound = MIN(M, *N);
+ dscal_(&bound, &pbar, q, &one);
+ bound = MIN(M, *N-d2);
+ daxpy_(&bound, &p1, qtemp, &one, q+d1, &one);
+ daxpy_(&bound, &p2, qtemp, &one, q+d2, &one);
+ M += d2;
+ }
+
+ /* correction for weight loss */
+ if(M > *N){
+ sum = 0;
+ for(j=0; j<MIN(M, *N); j++){
+ sum += q[j];
+ }
+ q[*N-1] += 1-sum;
+ }
+ Free(qtemp);
+}
+
void lossdistrib_Z(double *p, int *np, double *w, double *S, int *N, int *defaultflag,
double *rho, double *Z, double *wZ, int *nZ, double *q){
int i, j;
@@ -159,7 +213,7 @@ void lossdistrib_truncated(double *p, int *np, double *w, double *S, int *N,
Free(q2);
}
-void lossdistrib_joint( double *p, int *np, double *w, double *S, int *N, int *defaultflag, double *q) {
+void lossdistrib_joint(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
@@ -228,6 +282,79 @@ void lossdistrib_joint( double *p, int *np, double *w, double *S, int *N, int *d
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), 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<MIN(My, *N); n++){
+ memcpy(qtemp+n*(*N), q+n*(*N), MIN(Mx, *N) * sizeof(double));
+ }
+
+ bound = MIN(Mx, *N);
+ pbar = 1-p[k];
+ for(n=0; n<MIN(My, *N); n++){
+ dscal_(&bound, &pbar, q+(*N)*n, &one);
+ }
+ bound = MIN(Mx, *N-i-1);
+ for(n=0; n < MIN(My, *N-j-1); n++){
+ daxpy_(&bound, &w1, qtemp+(*N)*n, &one, q+i+(*N)*(j+n), &one);
+ daxpy_(&bound, &w2, qtemp+(*N)*n, &one, q+i+(*N)*(j+1+n), &one);
+ daxpy_(&bound, &w3, qtemp+(*N)*n, &one, q+i+1+(*N)*(j+1+n), &one);
+ daxpy_(&bound, &w4, qtemp+(*N)*n, &one, q+i+1+(*N)*(j+n), &one);
+ }
+ Mx += i+1;
+ My += j+1;
+ }
+ /* correction for weight loss */
+ if(Mx>*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.
@@ -298,6 +425,8 @@ void lossdistrib_prepay_joint(double *dp, double *pp, int *ndp, double *w,
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;
@@ -336,10 +465,9 @@ void lossdistrib_prepay_joint(double *dp, double *pp, int *ndp, double *w,
q[i+1+m+(*N)*(j1+n)] += dpw4 * qtemp[m+(*N)*n];
q[i+m+(*N)*(j2+n)] += ppw1 * qtemp[m+(*N)*n];
- temp = ppw2 * qtemp[m+(*N)*n];
- q[i+m+(*N)*(j2+1+n)] += temp;
+ q[i+m+(*N)*(j2+1+n)] += ppw2 * qtemp[m+(*N)*n];
q[i+1+m+(*N)*(j2+1+n)] += ppw3 * qtemp[m+(*N)*n];
- q[i+1+m+(*N)*(j2+n)] += temp;
+ q[i+1+m+(*N)*(j2+n)] += ppw2 * qtemp[m+(*N)*n];
}
}
}else{
@@ -370,6 +498,89 @@ void lossdistrib_prepay_joint(double *dp, double *pp, int *ndp, double *w,
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), double);
+ q[0] = 1;
+ Mx=1;
+ My=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<MIN(My, *N); n++){
+ memcpy(qtemp+n*(*N), q+n*(*N), MIN(Mx, *N) * sizeof(double));
+ }
+ bound = MIN(Mx, *N);
+ pbar = 1-dp[k]-pp[k];
+ for(n=0; n<MIN(My, *N); n++){
+ dscal_(&bound, &pbar, q+(*N)*n, &one);
+ }
+ if(*defaultflag){
+ ppw1 = alpha1 * alpha1 * pp[k];
+ ppw2 = alpha1 * alpha2 * pp[k];
+ ppw3 = alpha2 * alpha2 * pp[k];
+ for(n=0; n < MIN(My, *N-j2-1); n++){
+ daxpy_(&bound, &dpw1, qtemp+(*N)*n, &one, q+i+(*N)*(j1+n), &one);
+ daxpy_(&bound, &dpw2, qtemp+(*N)*n, &one, q+i+(*N)*(j1+1+n), &one);
+ daxpy_(&bound, &dpw3, qtemp+(*N)*n, &one, q+i+1+(*N)*(j1+1+n), &one);
+ daxpy_(&bound, &dpw4, qtemp+(*N)*n, &one, q+i+1+(*N)*(j1+n), &one);
+ daxpy_(&bound, &ppw1, qtemp+(*N)*n, &one, q+i+(*N)*(j2+n), &one);
+ }
+ }else{
+ bound = MIN(Mx, *N-i-1);
+ ppw1 = pp[k] * (j2+1-y2);
+ ppw2 = pp[k] * (y2-j2);
+ for(n=0; n < MIN(My, *N-j2-1); n++){
+ daxpy_(&bound, &dpw1, qtemp+(*N)*n, &one, q+i+(*N)*(j1+n), &one);
+ daxpy_(&bound, &dpw2, qtemp+(*N)*n, &one, q+i+(*N)*(j1+1+n), &one);
+ daxpy_(&bound, &dpw3, qtemp+(*N)*n, &one, q+i+1+(*N)*(j1+1+n), &one);
+ daxpy_(&bound, &dpw4, qtemp+(*N)*n, &one, q+i+1+(*N)*(j1+n), &one);
+ daxpy_(&bound, &ppw1, qtemp+(*N)*n, &one, q+(*N)*(j2+n), &one);
+ daxpy_(&bound, &ppw2, qtemp+(*N)*n, &one, q+(*N)*(j2+1+n), &one);
+ }
+ }
+ Mx += i + 1;
+ My += j2 + 1;
+ }
+ /* correction for weight loss */
+ if(Mx>*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){
return( pnorm( (qnorm(p, 0, 1, 1, 0) - sqrt(rho) * Z)/sqrt(1 - rho), 0, 1, 1, give_log));
}
@@ -461,8 +672,8 @@ void lossdistrib_prepay_joint_Z(double *dp, double *pp, int *ndp, double *w,
dpshocked[j + (*ndp) * i] = shockprob(dp[j], *rho, Z[i], 0);
ppshocked[j + (*ndp) * i] = shockprob(pp[j], *rho, -Z[i], 0);
}
- lossdistrib_prepay_joint(dpshocked + (*ndp) * i, ppshocked + (*ndp) * i, ndp,
- w, S + (*ndp) * i, N, defaultflag, qmat + N2 * i);
+ 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);
@@ -489,8 +700,8 @@ void lossdistrib_joint_Z(double *dp, int *ndp, double *w,
for(j = 0; j < *ndp; j++){
dpshocked[j + (*ndp) * i] = shockprob(dp[j], *rho, Z[i], 0);
}
- lossdistrib_joint(dpshocked + (*ndp) * i, ndp, w, S + (*ndp) * i, N,
- defaultflag, qmat + N2 * i);
+ 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);
@@ -541,9 +752,9 @@ void BClossdist(double *defaultprob, int *dim1, int *dim2,
Sshocked[j] = shockseverity(1-recov[j], Z[i], *rho, g);
Rshocked[j] = 1 - Sshocked[j];
}
- lossdistrib(gshocked, dim1, issuerweights, Sshocked, N, defaultflag,
+ lossdistrib_blas(gshocked, dim1, issuerweights, Sshocked, N, defaultflag,
Lw + i * (*N));
- lossdistrib(gshocked, dim1, issuerweights, Rshocked, N, defaultflag,
+ lossdistrib_blas(gshocked, dim1, issuerweights, Rshocked, N, defaultflag,
Rw + i * (*N));
}
dgemv_("n", N, n, &alpha, Lw, N, w, &one, &beta, L + t * (*N), &one);
diff --git a/R/tranche_functions.R b/R/tranche_functions.R
index 087f5d75..9ad62feb 100644
--- a/R/tranche_functions.R
+++ b/R/tranche_functions.R
@@ -241,6 +241,15 @@ lossdistC <- function(p, w, S, N, defaultflag=FALSE){
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
+ if(!is.loaded("lossdistrib_blas")){
+ dyn.load(file.path(root.dir, "code", "R", paste0("lossdistrib", .Platform$dynlib.ext)))
+ }
+ .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){
if(!is.loaded("lossdistrib_Z")){
dyn.load(file.path(root.dir, "code", "R", paste0("lossdistrib", .Platform$dynlib.ext)))
@@ -279,6 +288,16 @@ lossdistC.joint <- function(p, w, S, N, defaultflag=FALSE){
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
+ if(!is.loaded("lossdistrib_joint_blas")){
+ dyn.load(file.path(root.dir, "code", "R", paste0("lossdistrib", .Platform$dynlib.ext)))
+ }
+ .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
generated by cgit v1.2.3-70-g09d2 (git 2.52.0) at 2026-01-03 02:07:00 +0000