summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/lossdistrib.c1682
1 files changed, 841 insertions, 841 deletions
diff --git a/src/lossdistrib.c b/src/lossdistrib.c
index e5f2c53..61f3c84 100644
--- a/src/lossdistrib.c
+++ b/src/lossdistrib.c
@@ -1,841 +1,841 @@
-#include <R.h>
-#include <Rmath.h>
-#include <string.h>
-#include <omp.h>
-#include "lossdistrib.h"
-
-#define MIN(x, y) (((x) < (y)) ? (x) : (y))
-
-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 int dstev_(char* JOBZ, int* n, double* D, double* E, double* Z, int* ldz, double* WORK, int* INFO);
-
-extern void openblas_set_num_threads(int);
-
-void GHquad(int* n, double* Z, double* w) {
- // Setup for eigenvalue computations
- char JOBZ = 'V'; // Compute eigenvalues & vectors
- int INFO;
- int i;
- // Initialize array for workspace
- double * WORK = malloc(sizeof(double)*(2*(*n)-2));
-
- // Initialize array for eigenvectors
- double * V = malloc(sizeof(double)*(*n)*(*n));
-
- for(i = 0; i<(*n)-1; i++){
- w[i] = sqrt((i+1.)/2);
- }
-
- // Run eigen decomposition
- dstev_(&JOBZ, n, Z, w, V, n, WORK, &INFO);
-
- for (i=0; i<(*n); i++) {
- w[i] = V[i*(*n)] * V[i*(*n)];
- Z[i] *= sqrt(2);
- }
-
- // Deallocate temporary arrays
- free(WORK);
- free(V);
-}
-
-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
- np length of p
- w issuer weights
- S vector of severities (should be same length as p)
- N number of ticks in the grid
- defaultflag if true compute the default distribution
- q the loss distribution */
-
- int i, j, d1, d2, M;
- double lu, d, p1, p2, sum;
- double *qtemp;
-
- lu = 1./(*N-1);
- qtemp = calloc(*N, sizeof(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));
- for(j=0; j < MIN(M, *N); j++){
- q[j] = (1-p[i]) * q[j];
- }
- for(j=0; j < MIN(M, *N-d2); j++){
- q[d1+j] += p1 * qtemp[j];
- q[d2+j] += p2 * qtemp[j];
- };
- 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_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
- w: issuer weights
- S: vector of severities (should be same length as p)
- N: number of ticks in the grid
- defaultflag: 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;
- openblas_set_num_threads(1);
- lu = 1./(*N-1);
- qtemp = calloc(*N, sizeof(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, int *nZ, double *q){
- int i, j;
- double* pshocked = malloc(sizeof(double) * (*np) * (*nZ));
-
- #pragma omp parallel for private(j)
- for(i = 0; i < *nZ; i++){
- for(j = 0; j < *np; j++){
- pshocked[j + (*np) * i] = shockprob(p[j], rho[j], Z[i], 0);
- }
- lossdistrib_blas(pshocked + (*np) * i, np, w, S + (*np) * i, N,
- defaultflag, q + (*N) * i);
- }
- free(pshocked);
-}
-
-void lossdistrib_truncated(double *p, int *np, double *w, double *S, int *N,
- int *T, int *defaultflag, double *q) {
- /* recursive algorithm with first order correction for computing
- the loss distribution.
- input:
- 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
- T where to truncate the distribution
- defaultflag if true computes the default distribution
- output:
- q the loss distribution */
-
- int i, j, d1, d2, M;
- double lu, d, p1, p2;
- double *q1, *q2;
-
- lu = 1./(*N-1);
- q1 = calloc( *T, sizeof(double));
- q2 = calloc( *T, sizeof(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;
- for(j=0; j < MIN(M, *T); j++){
- q1[j] = p1 * q[j];
- q2[j] = p2 * q[j];
- q[j] = (1-p[i]) * q[j];
- }
- for(j=0; j < MIN(M, *T-d1); j++){
- q[d1+j] += q1[j];
- };
- for(j=0; j < MIN(M, *T-d2); j++){
- q[d2+j] += q2[j];
- };
- M += d2;
- }
- free(q1);
- free(q2);
-}
-
-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
- 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;
- double alpha1, alpha2, beta1, beta2;
- double w1, w2, w3, w4;
- double *qtemp;
-
- lu = 1./(*N-1);
- qtemp = calloc( (*N) * (*N), sizeof(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;
- w2 = alpha1 * beta2;
- w3 = alpha2 * beta2;
- w4 = alpha2 * beta1;
-
- for(n=0; n<MIN(My, *N); n++){
- memcpy(qtemp+n*(*N), q+n*(*N), MIN(Mx, *N) * sizeof(double));
- }
- for(n=0; n<MIN(My, *N); n++){
- for(m=0; m<MIN(Mx, *N); m++){
- q[m+(*N)*n] = (1-p[k])* q[m+(*N)*n];
- }
- }
- for(n=0; n < MIN(My, *N-j-1); n++){
- for(m=0; m < MIN(Mx, *N-i-1); m++){
- q[i+m+(*N)*(j+n)] += w1 * p[k] * qtemp[m+(*N)*n];
- q[i+m+(*N)*(j+1+n)] += w2 * p[k] * qtemp[m+(*N)*n];
- q[i+1+m+(*N)*(j+1+n)] += w3 * p[k] * qtemp[m+(*N)*n];
- q[i+1+m+(*N)*(j+n)] += w4 * p[k] *qtemp[m+(*N)*n];
- }
- }
- 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 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), sizeof(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.
- dp vector of default probabilities
- pp vector of prepay probabilities
- n length of p
- S vector of severities (should be same length as p)
- w vector of weights
- N number of ticks in the grid
- q the loss distribution */
-
- int i, j, d1l, d1u, d2l, d2u;
- int M;
- double lu, d1, d2, dp1, dp2, pp1, pp2, sum;
- double *qtemp;
-
- lu = 1./(*N - 1);
- qtemp = calloc( (*N), sizeof(double));
- q[0] = 1;
- M=1;
- for(i=0; i<(*n); i++){
- d1 = w[i] * (1-S[i]) /lu;
- d2 = w[i]/lu;
- d1l = floor(d1);
- d1u = d1l + 1;
- d2l = floor(d2);
- d2u = d2l + 1;
- dp1 = dp[i] * (d1u - d1);
- dp2 = dp[i] - dp1;
- pp1 = pp[i] * (d2u - d2);
- pp2 = pp[i] - pp1;
- memcpy(qtemp, q, MIN(M, *N) * sizeof(double));
- for(j = 0; j< MIN(M, *N); j++){
- q[j] = (1-dp[i]-pp[i]) * q[j];
- }
- for(j=0; j < MIN(M, *N-d2u); j++){
- q[d1l+j] += dp1 * qtemp[j];
- q[d1u+j] += dp2 * qtemp[j];
- q[d2l+j] += pp1 * qtemp[j];
- q[d2u+j] += pp2 * qtemp[j];
- };
- M += d2u;
- }
- /* 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_prepay_joint(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;
-
- lu = 1./(*N-1);
- qtemp = calloc((*N) * (*N), sizeof(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));
- }
- for(n=0; n<MIN(My, *N); n++){
- for(m=0; m<MIN(Mx, *N); m++){
- q[m+(*N)*n] = (1-dp[k]-pp[k]) * q[m+(*N)*n];
- }
- }
- 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++){
- for(m=0; m < MIN(Mx, *N-i-1); m++){
- q[i+m+(*N)*(j1+n)] += dpw1 * qtemp[m+(*N)*n];
- q[i+m+(*N)*(j1+1+n)] += dpw2 * qtemp[m+(*N)*n];
- q[i+1+m+(*N)*(j1+1+n)] += dpw3 * qtemp[m+(*N)*n];
- q[i+1+m+(*N)*(j1+n)] += dpw4 * qtemp[m+(*N)*n];
-
- q[i+m+(*N)*(j2+n)] += ppw1 * qtemp[m+(*N)*n];
- 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)] += ppw2 * qtemp[m+(*N)*n];
- }
- }
- }else{
- for(n=0; n < MIN(My, *N-j2-1); n++){
- for(m=0; m < MIN(Mx, *N-i-1); m++){
- q[i+m+(*N)*(j1+n)] += dpw1 * qtemp[m+(*N)*n];
- q[i+m+(*N)*(j1+1+n)] += dpw2 * qtemp[m+(*N)*n];
- q[i+1+m+(*N)*(j1+1+n)] += dpw3 * qtemp[m+(*N)*n];
- q[i+1+m+(*N)*(j1+n)] += dpw4 * qtemp[m+(*N)*n];
- q[m+(*N)*(j2+n)] += pp[k]*(j2+1-y2) * qtemp[m+(*N)*n];
- q[m+(*N)*(j2+1+n)] += pp[k]*(y2-j2) * qtemp[m+(*N)*n];
- }
- }
- }
- 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);
-}
-
-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), sizeof(double));
- 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;
- 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);
- }
- bound = MIN(Mx, *N-i-1);
- 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);
- daxpy_(&bound, &ppw2, qtemp+(*N)*n, &one, q+i+(*N)*(j2+1+n), &one);
- daxpy_(&bound, &ppw3, qtemp+(*N)*n, &one, q+i+1+(*N)*(j2+1+n), &one);
- daxpy_(&bound, &ppw2, qtemp+(*N)*n, &one, q+i+1+(*N)*(j2+n), &one);
- }
- }else{
- 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){
- if(rho==1){
- return((double)(Z<=qnorm(p, 0, 1, 1, 0)));
- }else{
- return( pnorm( (qnorm(p, 0, 1, 1, 0) - sqrt(rho) * Z)/sqrt(1 - rho), 0, 1, 1, give_log));
- }
-}
-
-double dqnorm(double x){
- return 1/dnorm(qnorm(x, 0, 1, 1, 0), 0, 1, 0);
-}
-
-double dshockprob(double p, double rho, double Z){
- return( dnorm((qnorm(p, 0, 1, 1, 0) - sqrt(rho) * Z)/sqrt(1-rho), 0, 1, 0) * dqnorm(p)/sqrt(1-rho) );
-}
-
-void shockprobvec2(double p, double rho, double* Z, int nZ, double *q){
- /* return a two column vectors with shockprob in the first column
- and dshockprob in the second column*/
- int i;
- #pragma omp parallel for
- for(i = 0; i < nZ; i++){
- q[i] = shockprob(p, rho, Z[i], 0);
- q[i + nZ] = dshockprob(p, rho, Z[i]);
- }
-}
-
-double shockseverity(double S, double Z, double rho, double p){
- 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){
- double cumw;
- int i;
- cumw = 0;
- for(i=0; i<nZ; i++){
- cumw += w[i];
- if(cumw > p0){
- break;
- }
- }
- return( Z[i] );
-}
-
-void fitprob(double* Z, double* w, int* nZ, double* rho, double* p0, double* result){
- double eps = 1e-12;
- int one = 1;
- double *q = malloc( 2 * (*nZ) * sizeof(double));
- double dp, p, phi;
-
- if(*p0 == 0){
- *result = 0;
- }else if(*rho == 1){
- *result = pnorm(quantile(Z, w, *nZ, *p0), 0, 1, 1, 0);
- }else{
- shockprobvec2(*p0, *rho, Z, *nZ, q);
- dp = (ddot_(nZ, q, &one, w, &one) - *p0)/ddot_(nZ, q + *nZ, &one, w, &one);
- p = *p0;
- while(fabs(dp) > eps){
- phi = 1;
- while( (p - phi * dp) < 0 || (p - phi * dp) > 1){
- phi *= 0.8;
- }
- p -= phi * dp;
- shockprobvec2(p, *rho, Z, *nZ, q);
- dp = (ddot_(nZ, q, &one, w, &one) - *p0)/ddot_(nZ, q + *nZ, &one, w, &one);
- }
- *result = p;
- }
- free(q);
-}
-
-void stochasticrecov(double* R, double* Rtilde, double* Z, double* w, int* nZ,
- double* rho, double* porig, double* pmod, double* q){
- double ptemp, ptilde;
- int i;
- if(*porig==0){
- for(i = 0; i < *nZ; i++){
- q[i] = *R;
- }
- }else{
- ptemp = (1 - *R) / (1 - *Rtilde) * *porig;
- fitprob(Z, w, nZ, rho, &ptemp, &ptilde);
- #pragma omp parallel for
- for(i = 0; i < *nZ; i++){
- q[i] = fabs(1 - (1 - *Rtilde) * exp( shockprob(ptilde, *rho, Z[i], 1) -
- shockprob(*pmod, *rho, Z[i], 1)));
- }
- }
-}
-
-void lossdistrib_prepay_joint_Z(double *dp, double *pp, int *ndp, double *w,
- double *S, int *N, int *defaultflag, double *rho,
- double *Z, double *wZ, int *nZ, double *q) {
- int i, j;
- double* dpshocked = malloc(sizeof(double) * (*ndp) * (*nZ));
- double* ppshocked = malloc(sizeof(double) * (*ndp) * (*nZ));
- int N2 = (*N) * (*N);
- double* qmat = malloc(sizeof(double) * N2 * (*nZ));
-
- double alpha = 1;
- double beta = 0;
- int one = 1;
-
-#pragma omp parallel for private(j)
- for(i = 0; i < *nZ; i++){
- for(j = 0; j < *ndp; j++){
- dpshocked[j + (*ndp) * i] = shockprob(dp[j], rho[j], Z[i], 0);
- ppshocked[j + (*ndp) * i] = shockprob(pp[j], rho[j], -Z[i], 0);
- }
- 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);
-
- free(dpshocked);
- free(ppshocked);
- free(qmat);
-}
-
-void lossdistrib_joint_Z(double *dp, int *ndp, double *w,
- double *S, int *N, int *defaultflag, double *rho,
- double *Z, double *wZ, int *nZ, double *q) {
- int i, j;
- double* dpshocked = malloc(sizeof(double) * (*ndp) * (*nZ));
- int N2 = (*N) * (*N);
- double* qmat = malloc(sizeof(double) * N2 * (*nZ));
-
- double alpha = 1;
- double beta = 0;
- int one = 1;
-
-#pragma omp parallel for private(j)
- for(i = 0; i < *nZ; i++){
- for(j = 0; j < *ndp; j++){
- dpshocked[j + (*ndp) * i] = shockprob(dp[j], rho[j], Z[i], 0);
- }
- 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);
-
- free(dpshocked);
- free(qmat);
-}
-
-void BCloss_recov_dist(double *defaultprob, int *dim1, int *dim2,
- double *issuerweights, double *recov, double *Z, double *w,
- int *n, double *rho, int *N, int *defaultflag,
- double *L, double *R) {
- /*
- computes the loss and recovery distribution over time with a flat gaussian
- correlation
- 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 dim1)
- recov: vector of recoveries (length dim1)
- Z: vector of factor values (length n)
- w: vector of factor weights (length n)
- rho: correlation beta vector (length dim1)
- N: number of ticks in the grid
- defaultflag: if true, computes the default distribution
- outputs:
- L: matrix of size (N, dim2)
- R: matrix of size (N, dim2)
- */
- int t, i, j;
- double g;
- double *gshocked, *Rshocked, *Sshocked, *Lw, *Rw;
- int one = 1;
- double alpha = 1;
- double beta = 0;
-
- gshocked = malloc((*dim1) * (*n) * sizeof(double));
- Rshocked = malloc((*dim1) * (*n) * sizeof(double));
- Sshocked = malloc((*dim1) * (*n) * sizeof(double));
- Lw = malloc((*N) * (*n) * sizeof(double));
- Rw = malloc((*N) * (*n) * sizeof(double));
-
-
- for(t=0; t < (*dim2); t++) {
- memset(Lw, 0, (*N) * (*n) * sizeof(double));
- memset(Rw, 0, (*N) * (*n) * sizeof(double));
- #pragma omp parallel for private(j, g)
- for(i=0; i < *n; i++){
- for(j=0; j < (*dim1); j++){
- g = defaultprob[j + (*dim1) * t];
- gshocked[j+(*dim1)*i] = shockprob(g, rho[j], Z[i], 0);
- Sshocked[j+(*dim1)*i] = shockseverity(1-recov[j], Z[i], rho[j], g);
- Rshocked[j+(*dim1)*i] = 1 - Sshocked[j+(*dim1)*i];
- }
- lossdistrib_blas(gshocked + (*dim1) * i, dim1, issuerweights, Sshocked + (*dim1)*i, N, defaultflag,
- Lw + i * (*N));
- lossdistrib_blas(gshocked + (*dim1) * i, dim1, issuerweights, Rshocked + (*dim1)*i, N, defaultflag,
- Rw + i * (*N));
- }
- dgemv_("n", N, n, &alpha, Lw, N, w, &one, &beta, L + t * (*N), &one);
- dgemv_("n", N, n, &alpha, Rw, N, w, &one, &beta, R + t * (*N), &one);
- }
- free(gshocked);
- free(Rshocked);
- free(Sshocked);
- free(Lw);
- free(Rw);
-}
-
-
-void BCloss_dist(double *defaultprob, int *dim1, int *dim2,
- double *issuerweights, double *recov, double *Z, double *w,
- int *n, double *rho, int *N, int *defaultflag,
- double *L) {
- /*
- computes the loss distribution over time with a flat gaussian
- correlation
- 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 dim1)
- recov: vector of recoveries (length dim1)
- Z: vector of factor values (length n)
- w: vector of factor weights (length n)
- rho: correlation beta vector (length dim1)
- N: number of ticks in the grid
- defaultflag: if true, computes the default distribution
- outputs:
- L: matrix of size (N, dim2)
- */
- int t, i, j;
- double g;
- double *gshocked, *Sshocked, *Lw;
- int one = 1;
- double alpha = 1;
- double beta = 0;
-
- gshocked = malloc((*dim1) * (*n) * sizeof(double));
- Sshocked = malloc((*dim1) * (*n) * sizeof(double));
- Lw = malloc((*N) * (*n) * sizeof(double));
-
- for(t=0; t < (*dim2); t++) {
- memset(Lw, 0, (*N) * (*n) * sizeof(double));
- #pragma omp parallel for private(j, g)
- for(i=0; i < *n; i++){
- for(j=0; j < (*dim1); j++){
- g = defaultprob[j + (*dim1) * t];
- gshocked[j+(*dim1)*i] = shockprob(g, rho[j], Z[i], 0);
- Sshocked[j+(*dim1)*i] = shockseverity(1-recov[j], Z[i], rho[j], g);
- }
- lossdistrib_blas(gshocked + (*dim1) * i, dim1, issuerweights, Sshocked + (*dim1)*i, N, defaultflag,
- Lw + i * (*N));
- }
- dgemv_("n", N, n, &alpha, Lw, N, w, &one, &beta, L + t * (*N), &one);
- }
- free(gshocked);
- free(Sshocked);
- free(Lw);
-}
+#include <R.h>
+#include <Rmath.h>
+#include <string.h>
+#include <omp.h>
+#include "lossdistrib.h"
+
+#define MIN(x, y) (((x) < (y)) ? (x) : (y))
+
+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 int dstev_(char* JOBZ, int* n, double* D, double* E, double* Z, int* ldz, double* WORK, int* INFO);
+
+extern void openblas_set_num_threads(int);
+
+void GHquad(int* n, double* Z, double* w) {
+ // Setup for eigenvalue computations
+ char JOBZ = 'V'; // Compute eigenvalues & vectors
+ int INFO;
+ int i;
+ // Initialize array for workspace
+ double * WORK = malloc(sizeof(double)*(2*(*n)-2));
+
+ // Initialize array for eigenvectors
+ double * V = malloc(sizeof(double)*(*n)*(*n));
+
+ for(i = 0; i<(*n)-1; i++){
+ w[i] = sqrt((i+1.)/2);
+ }
+
+ // Run eigen decomposition
+ dstev_(&JOBZ, n, Z, w, V, n, WORK, &INFO);
+
+ for (i=0; i<(*n); i++) {
+ w[i] = V[i*(*n)] * V[i*(*n)];
+ Z[i] *= sqrt(2);
+ }
+
+ // Deallocate temporary arrays
+ free(WORK);
+ free(V);
+}
+
+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
+ np length of p
+ w issuer weights
+ S vector of severities (should be same length as p)
+ N number of ticks in the grid
+ defaultflag if true compute the default distribution
+ q the loss distribution */
+
+ int i, j, d1, d2, M;
+ double lu, d, p1, p2, sum;
+ double *qtemp;
+
+ lu = 1./(*N-1);
+ qtemp = calloc(*N, sizeof(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));
+ for(j=0; j < MIN(M, *N); j++){
+ q[j] = (1-p[i]) * q[j];
+ }
+ for(j=0; j < MIN(M, *N-d2); j++){
+ q[d1+j] += p1 * qtemp[j];
+ q[d2+j] += p2 * qtemp[j];
+ }
+ 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_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
+ w: issuer weights
+ S: vector of severities (should be same length as p)
+ N: number of ticks in the grid
+ defaultflag: 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;
+ openblas_set_num_threads(1);
+ lu = 1./(*N-1);
+ qtemp = calloc(*N, sizeof(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, int *nZ, double *q){
+ int i, j;
+ double* pshocked = malloc(sizeof(double) * (*np) * (*nZ));
+
+#pragma omp parallel for private(j)
+ for(i = 0; i < *nZ; i++){
+ for(j = 0; j < *np; j++){
+ pshocked[j + (*np) * i] = shockprob(p[j], rho[j], Z[i], 0);
+ }
+ lossdistrib_blas(pshocked + (*np) * i, np, w, S + (*np) * i, N,
+ defaultflag, q + (*N) * i);
+ }
+ free(pshocked);
+}
+
+void lossdistrib_truncated(double *p, int *np, double *w, double *S, int *N,
+ int *T, int *defaultflag, double *q) {
+ /* recursive algorithm with first order correction for computing
+ the loss distribution.
+ input:
+ 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
+ T where to truncate the distribution
+ defaultflag if true computes the default distribution
+ output:
+ q the loss distribution */
+
+ int i, j, d1, d2, M;
+ double lu, d, p1, p2;
+ double *q1, *q2;
+
+ lu = 1./(*N-1);
+ q1 = calloc( *T, sizeof(double));
+ q2 = calloc( *T, sizeof(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;
+ for(j=0; j < MIN(M, *T); j++){
+ q1[j] = p1 * q[j];
+ q2[j] = p2 * q[j];
+ q[j] = (1-p[i]) * q[j];
+ }
+ for(j=0; j < MIN(M, *T-d1); j++){
+ q[d1+j] += q1[j];
+ };
+ for(j=0; j < MIN(M, *T-d2); j++){
+ q[d2+j] += q2[j];
+ };
+ M += d2;
+ }
+ free(q1);
+ free(q2);
+}
+
+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
+ 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;
+ double alpha1, alpha2, beta1, beta2;
+ double w1, w2, w3, w4;
+ double *qtemp;
+
+ lu = 1./(*N-1);
+ qtemp = calloc( (*N) * (*N), sizeof(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;
+ w2 = alpha1 * beta2;
+ w3 = alpha2 * beta2;
+ w4 = alpha2 * beta1;
+
+ for(n=0; n<MIN(My, *N); n++){
+ memcpy(qtemp+n*(*N), q+n*(*N), MIN(Mx, *N) * sizeof(double));
+ }
+ for(n=0; n<MIN(My, *N); n++){
+ for(m=0; m<MIN(Mx, *N); m++){
+ q[m+(*N)*n] = (1-p[k])* q[m+(*N)*n];
+ }
+ }
+ for(n=0; n < MIN(My, *N-j-1); n++){
+ for(m=0; m < MIN(Mx, *N-i-1); m++){
+ q[i+m+(*N)*(j+n)] += w1 * p[k] * qtemp[m+(*N)*n];
+ q[i+m+(*N)*(j+1+n)] += w2 * p[k] * qtemp[m+(*N)*n];
+ q[i+1+m+(*N)*(j+1+n)] += w3 * p[k] * qtemp[m+(*N)*n];
+ q[i+1+m+(*N)*(j+n)] += w4 * p[k] *qtemp[m+(*N)*n];
+ }
+ }
+ 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 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), sizeof(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.
+ dp vector of default probabilities
+ pp vector of prepay probabilities
+ n length of p
+ S vector of severities (should be same length as p)
+ w vector of weights
+ N number of ticks in the grid
+ q the loss distribution */
+
+ int i, j, d1l, d1u, d2l, d2u;
+ int M;
+ double lu, d1, d2, dp1, dp2, pp1, pp2, sum;
+ double *qtemp;
+
+ lu = 1./(*N - 1);
+ qtemp = calloc( (*N), sizeof(double));
+ q[0] = 1;
+ M=1;
+ for(i=0; i<(*n); i++){
+ d1 = w[i] * (1-S[i]) /lu;
+ d2 = w[i]/lu;
+ d1l = floor(d1);
+ d1u = d1l + 1;
+ d2l = floor(d2);
+ d2u = d2l + 1;
+ dp1 = dp[i] * (d1u - d1);
+ dp2 = dp[i] - dp1;
+ pp1 = pp[i] * (d2u - d2);
+ pp2 = pp[i] - pp1;
+ memcpy(qtemp, q, MIN(M, *N) * sizeof(double));
+ for(j = 0; j< MIN(M, *N); j++){
+ q[j] = (1-dp[i]-pp[i]) * q[j];
+ }
+ for(j=0; j < MIN(M, *N-d2u); j++){
+ q[d1l+j] += dp1 * qtemp[j];
+ q[d1u+j] += dp2 * qtemp[j];
+ q[d2l+j] += pp1 * qtemp[j];
+ q[d2u+j] += pp2 * qtemp[j];
+ };
+ M += d2u;
+ }
+ /* 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_prepay_joint(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;
+
+ lu = 1./(*N-1);
+ qtemp = calloc((*N) * (*N), sizeof(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));
+ }
+ for(n=0; n<MIN(My, *N); n++){
+ for(m=0; m<MIN(Mx, *N); m++){
+ q[m+(*N)*n] = (1-dp[k]-pp[k]) * q[m+(*N)*n];
+ }
+ }
+ 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++){
+ for(m=0; m < MIN(Mx, *N-i-1); m++){
+ q[i+m+(*N)*(j1+n)] += dpw1 * qtemp[m+(*N)*n];
+ q[i+m+(*N)*(j1+1+n)] += dpw2 * qtemp[m+(*N)*n];
+ q[i+1+m+(*N)*(j1+1+n)] += dpw3 * qtemp[m+(*N)*n];
+ q[i+1+m+(*N)*(j1+n)] += dpw4 * qtemp[m+(*N)*n];
+
+ q[i+m+(*N)*(j2+n)] += ppw1 * qtemp[m+(*N)*n];
+ 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)] += ppw2 * qtemp[m+(*N)*n];
+ }
+ }
+ }else{
+ for(n=0; n < MIN(My, *N-j2-1); n++){
+ for(m=0; m < MIN(Mx, *N-i-1); m++){
+ q[i+m+(*N)*(j1+n)] += dpw1 * qtemp[m+(*N)*n];
+ q[i+m+(*N)*(j1+1+n)] += dpw2 * qtemp[m+(*N)*n];
+ q[i+1+m+(*N)*(j1+1+n)] += dpw3 * qtemp[m+(*N)*n];
+ q[i+1+m+(*N)*(j1+n)] += dpw4 * qtemp[m+(*N)*n];
+ q[m+(*N)*(j2+n)] += pp[k]*(j2+1-y2) * qtemp[m+(*N)*n];
+ q[m+(*N)*(j2+1+n)] += pp[k]*(y2-j2) * qtemp[m+(*N)*n];
+ }
+ }
+ }
+ 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);
+}
+
+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), sizeof(double));
+ 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;
+ 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);
+ }
+ bound = MIN(Mx, *N-i-1);
+ 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);
+ daxpy_(&bound, &ppw2, qtemp+(*N)*n, &one, q+i+(*N)*(j2+1+n), &one);
+ daxpy_(&bound, &ppw3, qtemp+(*N)*n, &one, q+i+1+(*N)*(j2+1+n), &one);
+ daxpy_(&bound, &ppw2, qtemp+(*N)*n, &one, q+i+1+(*N)*(j2+n), &one);
+ }
+ }else{
+ 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){
+ if(rho==1){
+ return((double)(Z<=qnorm(p, 0, 1, 1, 0)));
+ }else{
+ return( pnorm( (qnorm(p, 0, 1, 1, 0) - sqrt(rho) * Z)/sqrt(1 - rho), 0, 1, 1, give_log));
+ }
+}
+
+double dqnorm(double x){
+ return 1/dnorm(qnorm(x, 0, 1, 1, 0), 0, 1, 0);
+}
+
+double dshockprob(double p, double rho, double Z){
+ return( dnorm((qnorm(p, 0, 1, 1, 0) - sqrt(rho) * Z)/sqrt(1-rho), 0, 1, 0) * dqnorm(p)/sqrt(1-rho) );
+}
+
+void shockprobvec2(double p, double rho, double* Z, int nZ, double *q){
+ /* return a two column vectors with shockprob in the first column
+ and dshockprob in the second column*/
+ int i;
+ #pragma omp parallel for
+ for(i = 0; i < nZ; i++){
+ q[i] = shockprob(p, rho, Z[i], 0);
+ q[i + nZ] = dshockprob(p, rho, Z[i]);
+ }
+}
+
+double shockseverity(double S, double Z, double rho, double p){
+ 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){
+ double cumw;
+ int i;
+ cumw = 0;
+ for(i=0; i<nZ; i++){
+ cumw += w[i];
+ if(cumw > p0){
+ break;
+ }
+ }
+ return( Z[i] );
+}
+
+void fitprob(double* Z, double* w, int* nZ, double* rho, double* p0, double* result){
+ double eps = 1e-12;
+ int one = 1;
+ double *q = malloc( 2 * (*nZ) * sizeof(double));
+ double dp, p, phi;
+
+ if(*p0 == 0){
+ *result = 0;
+ }else if(*rho == 1){
+ *result = pnorm(quantile(Z, w, *nZ, *p0), 0, 1, 1, 0);
+ }else{
+ shockprobvec2(*p0, *rho, Z, *nZ, q);
+ dp = (ddot_(nZ, q, &one, w, &one) - *p0)/ddot_(nZ, q + *nZ, &one, w, &one);
+ p = *p0;
+ while(fabs(dp) > eps){
+ phi = 1;
+ while( (p - phi * dp) < 0 || (p - phi * dp) > 1){
+ phi *= 0.8;
+ }
+ p -= phi * dp;
+ shockprobvec2(p, *rho, Z, *nZ, q);
+ dp = (ddot_(nZ, q, &one, w, &one) - *p0)/ddot_(nZ, q + *nZ, &one, w, &one);
+ }
+ *result = p;
+ }
+ free(q);
+}
+
+void stochasticrecov(double* R, double* Rtilde, double* Z, double* w, int* nZ,
+ double* rho, double* porig, double* pmod, double* q){
+ double ptemp, ptilde;
+ int i;
+ if(*porig==0){
+ for(i = 0; i < *nZ; i++){
+ q[i] = *R;
+ }
+ }else{
+ ptemp = (1 - *R) / (1 - *Rtilde) * *porig;
+ fitprob(Z, w, nZ, rho, &ptemp, &ptilde);
+#pragma omp parallel for
+ for(i = 0; i < *nZ; i++){
+ q[i] = fabs(1 - (1 - *Rtilde) * exp( shockprob(ptilde, *rho, Z[i], 1) -
+ shockprob(*pmod, *rho, Z[i], 1)));
+ }
+ }
+}
+
+void lossdistrib_prepay_joint_Z(double *dp, double *pp, int *ndp, double *w,
+ double *S, int *N, int *defaultflag, double *rho,
+ double *Z, double *wZ, int *nZ, double *q) {
+ int i, j;
+ double* dpshocked = malloc(sizeof(double) * (*ndp) * (*nZ));
+ double* ppshocked = malloc(sizeof(double) * (*ndp) * (*nZ));
+ int N2 = (*N) * (*N);
+ double* qmat = malloc(sizeof(double) * N2 * (*nZ));
+
+ double alpha = 1;
+ double beta = 0;
+ int one = 1;
+
+ #pragma omp parallel for private(j)
+ for(i = 0; i < *nZ; i++){
+ for(j = 0; j < *ndp; j++){
+ dpshocked[j + (*ndp) * i] = shockprob(dp[j], rho[j], Z[i], 0);
+ ppshocked[j + (*ndp) * i] = shockprob(pp[j], rho[j], -Z[i], 0);
+ }
+ 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);
+
+ free(dpshocked);
+ free(ppshocked);
+ free(qmat);
+}
+
+void lossdistrib_joint_Z(double *dp, int *ndp, double *w,
+ double *S, int *N, int *defaultflag, double *rho,
+ double *Z, double *wZ, int *nZ, double *q) {
+ int i, j;
+ double* dpshocked = malloc(sizeof(double) * (*ndp) * (*nZ));
+ int N2 = (*N) * (*N);
+ double* qmat = malloc(sizeof(double) * N2 * (*nZ));
+
+ double alpha = 1;
+ double beta = 0;
+ int one = 1;
+
+#pragma omp parallel for private(j)
+ for(i = 0; i < *nZ; i++){
+ for(j = 0; j < *ndp; j++){
+ dpshocked[j + (*ndp) * i] = shockprob(dp[j], rho[j], Z[i], 0);
+ }
+ 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);
+
+ free(dpshocked);
+ free(qmat);
+}
+
+void BCloss_recov_dist(double *defaultprob, int *dim1, int *dim2,
+ double *issuerweights, double *recov, double *Z, double *w,
+ int *n, double *rho, int *N, int *defaultflag,
+ double *L, double *R) {
+ /*
+ computes the loss and recovery distribution over time with a flat gaussian
+ correlation
+ 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 dim1)
+ recov: vector of recoveries (length dim1)
+ Z: vector of factor values (length n)
+ w: vector of factor weights (length n)
+ rho: correlation beta vector (length dim1)
+ N: number of ticks in the grid
+ defaultflag: if true, computes the default distribution
+ outputs:
+ L: matrix of size (N, dim2)
+ R: matrix of size (N, dim2)
+ */
+ int t, i, j;
+ double g;
+ double *gshocked, *Rshocked, *Sshocked, *Lw, *Rw;
+ int one = 1;
+ double alpha = 1;
+ double beta = 0;
+
+ gshocked = malloc((*dim1) * (*n) * sizeof(double));
+ Rshocked = malloc((*dim1) * (*n) * sizeof(double));
+ Sshocked = malloc((*dim1) * (*n) * sizeof(double));
+ Lw = malloc((*N) * (*n) * sizeof(double));
+ Rw = malloc((*N) * (*n) * sizeof(double));
+
+
+ for(t=0; t < (*dim2); t++) {
+ memset(Lw, 0, (*N) * (*n) * sizeof(double));
+ memset(Rw, 0, (*N) * (*n) * sizeof(double));
+ #pragma omp parallel for private(j, g)
+ for(i=0; i < *n; i++){
+ for(j=0; j < (*dim1); j++){
+ g = defaultprob[j + (*dim1) * t];
+ gshocked[j+(*dim1)*i] = shockprob(g, rho[j], Z[i], 0);
+ Sshocked[j+(*dim1)*i] = shockseverity(1-recov[j], Z[i], rho[j], g);
+ Rshocked[j+(*dim1)*i] = 1 - Sshocked[j+(*dim1)*i];
+ }
+ lossdistrib_blas(gshocked + (*dim1) * i, dim1, issuerweights, Sshocked + (*dim1)*i, N, defaultflag,
+ Lw + i * (*N));
+ lossdistrib_blas(gshocked + (*dim1) * i, dim1, issuerweights, Rshocked + (*dim1)*i, N, defaultflag,
+ Rw + i * (*N));
+ }
+ dgemv_("n", N, n, &alpha, Lw, N, w, &one, &beta, L + t * (*N), &one);
+ dgemv_("n", N, n, &alpha, Rw, N, w, &one, &beta, R + t * (*N), &one);
+ }
+ free(gshocked);
+ free(Rshocked);
+ free(Sshocked);
+ free(Lw);
+ free(Rw);
+}
+
+
+void BCloss_dist(double *defaultprob, int *dim1, int *dim2,
+ double *issuerweights, double *recov, double *Z, double *w,
+ int *n, double *rho, int *N, int *defaultflag,
+ double *L) {
+ /*
+ computes the loss distribution over time with a flat gaussian
+ correlation
+ 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 dim1)
+ recov: vector of recoveries (length dim1)
+ Z: vector of factor values (length n)
+ w: vector of factor weights (length n)
+ rho: correlation beta vector (length dim1)
+ N: number of ticks in the grid
+ defaultflag: if true, computes the default distribution
+ outputs:
+ L: matrix of size (N, dim2)
+ */
+ int t, i, j;
+ double g;
+ double *gshocked, *Sshocked, *Lw;
+ int one = 1;
+ double alpha = 1;
+ double beta = 0;
+
+ gshocked = malloc((*dim1) * (*n) * sizeof(double));
+ Sshocked = malloc((*dim1) * (*n) * sizeof(double));
+ Lw = malloc((*N) * (*n) * sizeof(double));
+
+ for(t=0; t < (*dim2); t++) {
+ memset(Lw, 0, (*N) * (*n) * sizeof(double));
+ #pragma omp parallel for private(j, g)
+ for(i=0; i < *n; i++){
+ for(j=0; j < (*dim1); j++){
+ g = defaultprob[j + (*dim1) * t];
+ gshocked[j+(*dim1)*i] = shockprob(g, rho[j], Z[i], 0);
+ Sshocked[j+(*dim1)*i] = shockseverity(1-recov[j], Z[i], rho[j], g);
+ }
+ lossdistrib_blas(gshocked + (*dim1) * i, dim1, issuerweights, Sshocked + (*dim1)*i, N, defaultflag,
+ Lw + i * (*N));
+ }
+ dgemv_("n", N, n, &alpha, Lw, N, w, &one, &beta, L + t * (*N), &one);
+ }
+ free(gshocked);
+ free(Sshocked);
+ free(Lw);
+}
generated by cgit v1.2.3-70-g09d2 (git 2.52.0) at 2025-12-21 21:13:10 +0000