diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/lossdistrib.c | 760 | ||||
| -rw-r--r-- | R/lossdistrib.h | 46 | ||||
| -rw-r--r-- | R/tranche_functions.R | 886 |
3 files changed, 0 insertions, 1692 deletions
diff --git a/R/lossdistrib.c b/R/lossdistrib.c deleted file mode 100644 index 77ea335c..00000000 --- a/R/lossdistrib.c +++ /dev/null @@ -1,760 +0,0 @@ -#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 void openblas_set_num_threads(int);
-
-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 BClossdist(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);
-}
diff --git a/R/lossdistrib.h b/R/lossdistrib.h deleted file mode 100644 index d4bd974d..00000000 --- a/R/lossdistrib.h +++ /dev/null @@ -1,46 +0,0 @@ -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);
-
-void lossdistrib_Z(double *p, int *np, double *w, double *S, int *N, int *defaultflag,
- double *rho, double *Z, int *nZ, double *q);
-
-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);
-
-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);
-
-void lossdistrib_prepay_joint(double *dp, double *pp, int *ndp, double *w,
- double *S, int *N, int *defaultflag, double *q);
-double dqnorm(double x);
-
-double dshockprob(double p, double rho, double Z);
-
-void shockprobvec2(double p, double rho, double* Z, int nZ, double *q);
-
-double shockseverity(double S, double Z, double rho, double p);
-
-void fitprob(double* Z, double* w, int* nZ, double* rho, double* p0, double* result);
-
-void stochasticrecov(double* R, double* Rtilde, double* Z, double* w, int* nZ,
- double* rho, double* porig, double* pmod, double* q);
-
-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);
-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);
-
-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);
-
-double quantile(double* Z, double* w, int nZ, double p0);
diff --git a/R/tranche_functions.R b/R/tranche_functions.R deleted file mode 100644 index c969c1fd..00000000 --- a/R/tranche_functions.R +++ /dev/null @@ -1,886 +0,0 @@ -## todo:
-## -investigate other ways to interpolate the random severities on the grid
-## I'm thinking that at eah severity that we add to the distribution, round it down
-## and keep track of the missing mass: namely if X_i=S_i w.p p_i, then add
-## X_i=lu*floor(S_i/lu) with probability p_i and propagate
-## X_{i+1}=S_{i+1}+(S_i-lu*floor(S_i/lu)) with the right probability
-## - investigate truncated distributions more (need to compute loss and recov distribution
-## separately, for the 0-10 equity tranche, we need the loss on the 0-0.1 support and
-## recovery with 0.1-1 support, so it's not clear that there is a big gain.
-## - do the correlation adjustments when computing the deltas since it seems to be
-## the market standard
-hostname <- system("hostname", intern=TRUE)
-
-checkSymbol <- function(name){
- if(!is.loaded(name)){
- if(.Platform$OS.type == "unix"){
- root.dir <- "/home/share/CorpCDOs"
- dyn.load(file.path(root.dir, "code", "R", paste0("lossdistrib",
- hostname,
- .Platform$dynlib.ext)))
- }else{
- root.dir <- "//WDSENTINEL/share/CorpCDOs"
- dyn.load(file.path(root.dir, "code", "R", paste0("lossdistrib",
- .Platform$dynlib.ext)))
- }
- }
-}
-lossdistrib <- function(p){
- ## basic recursive algorithm of Andersen, Sidenius and Basu
- 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]
- }
- return(q)
-}
-
-lossdistrib.fft <- function(p){
- ## computes loss distribution using the fft
- ## complexity is of order O(n*m)+O(m*log(m))
- ## where m is the size of the grid and n the number of probabilities.
- ## this is slower than the recursive algorithm
- theta <- 2*pi*1i*(0:n)/(n+1)
- Phi <- 1 - p + p%o%exp(theta)
- v <- apply(Phi, 2, prod)
- return(1/(n+1)*Re(fft(v)))
-}
-
-lossdistrib2 <- function(p, w, S, N, defaultflag=FALSE){
- ## recursive algorithm with first order correction
- ## p vector of default probabilities
- ## w vector of weigths
- ## S vector of severities
- ## N number of ticks in the grid
- ## defaultflag if true computes the default distribution
- n <- length(p)
- lu <- 1/(N-1)
- q <- rep(0, N)
- q[1] <- 1
- for(i in 1:n){
- if(defaultflag){
- d <- w[i] /lu
- }else{
- d <- S[i] * w[i] / lu
- }
- d1 <- floor(d)
- d2 <- ceiling(d)
- p1 <- p[i]*(d2-d)
- p2 <- p[i] - p1
- q1 <- c(rep(0,d1), p1*q[1:(N-d1)])
- q2 <- c(rep(0,d2), p2*q[1:(N-d2)])
- q <- q1 + q2 + (1-p[i])*q
- }
- q[length(q)] <- q[length(q)]+1-sum(q)
- return(q)
-}
-
-lossdistrib2.truncated <- function(p, w, S, N, cutoff=N){
- ## recursive algorithm with first order correction
- ## p vector of default probabilities
- ## w vector of weigths
- ## S vector of severities
- ## N number of ticks in the grid (for best accuracy should
- ## be a multiple of the number of issuers)
- ## cutoff where to stop computing the exact probabilities
- ## (useful for tranche computations)
-
- ## this is actually slower than lossdistrib2. But in C this is
- ## twice as fast.
- ## for high severities, M can become bigger than N, and there is
- ## some probability mass escaping.
- n <- length(p)
- lu <- 1/(N-1)
- q <- rep(0, truncated)
- q[1] <- 1
- M <- 1
- for(i in 1:n){
- d <- S[i] * w[i] / lu
- d1 <- floor(d)
- d2 <- ceiling(d)
- p1 <- p[i]*(d2-d)
- p2 <- p[i] - p1
- q1 <- p1*q[1:min(M, cutoff-d1)]
- q2 <- p2*q[1:min(M, cutoff-d2)]
- q[1:min(M, cutoff)] <- (1-p[i])*q[1:min(M, cutoff)]
- q[(d1+1):min(M+d1, cutoff)] <- q[(d1+1):min(M+d1, cutoff)]+q1
- q[(d2+1):min(M+d2, cutoff)] <- q[(d2+1):min(M+d2, cutoff)]+q2
- M <- M+d2
- }
- return(q)
-}
-
-recovdist <- function(dp, pp, w, S, N){
- ## computes the recovery distribution for a sum of independent variables
- ## R=\sum_{i=1}^n X_i
- ## where X_i = 0 w.p 1-dp[i]-pp[i]
- ## = w[i]*(1-S[i]) w.p dp[i]
- ## = w[i] w.p pp[i]
- ## each non zero value v is interpolated on the grid as
- ## the pair of values floor(v/lu) and ceiling(v/lu) so that
- ## X_i has four non zero values
- n <- length(dp)
- q <- rep(0, N)
- q[1] <- 1
- lu <- 1/(N-1)
- for(i in 1:n){
- d1 <- w[i]*(1-S[i])/lu
- d1l <- floor(d1)
- d1u <- ceiling(d1)
- d2 <- w[i] / lu
- d2l <- floor(d2)
- d2u <- ceiling(d2)
- dp1 <- dp[i] * (d1u-d1)
- dp2 <- dp[i] - dp1
- pp1 <- pp[i] * (d2u - d2)
- pp2 <- pp[i] - pp1
- q1 <- c(rep(0, d1l), dp1 * q[1:(N-d1l)])
- q2 <- c(rep(0, d1u), dp2 * q[1:(N-d1u)])
- q3 <- c(rep(0, d2l), pp1 * q[1:(N-d2l)])
- q4 <- c(rep(0, d2u), pp2 *q[1:(N-d2u)])
- q <- q1+q2+q3+q4+(1-dp[i]-pp[i])*q
- }
- return(q)
-}
-
-lossdist.joint <- function(p, w, S, N, defaultflag=FALSE){
- ## recursive algorithm with first order correction
- ## to compute the joint probability distribution of the loss and recovery
- ## inputs:
- ## p: vector of default probabilities
- ## w: vector of issuer weights
- ## S: vector of severities
- ## N: number of tick sizes on the grid
- ## defaultflag: if true computes the default distribution
- ## output:
- ## q: matrix of joint loss, recovery probability
- ## colSums(q) is the recovery distribution marginal
- ## rowSums(q) is the loss distribution marginal
- n <- length(p)
- lu <- 1/(N-1)
- q <- matrix(0, N, N)
- q[1,1] <- 1
- for(k in 1:n){
- if(defaultflag){
- x <- w[k] / lu
- }else{
- x <- S[k] * w[k]/lu
- }
- y <- (1-S[k]) * w[k]/lu
- i <- floor(x)
- j <- floor(y)
- weights <- c((i+1-x)*(j+1-y), (i+1-x)*(y-j), (x-i)*(y-j), (j+1-y)*(x-i))
- psplit <- p[k] * weights
- qtemp <- matrix(0, N, N)
- qtemp[(i+1):N,(j+1):N] <- qtemp[(i+1):N,(j+1):N] + psplit[1] * q[1:(N-i),1:(N-j)]
- qtemp[(i+1):N,(j+2):N] <- qtemp[(i+1):N,(j+2):N] + psplit[2] * q[1:(N-i), 1:(N-j-1)]
- qtemp[(i+2):N,(j+2):N] <- qtemp[(i+2):N,(j+2):N] + psplit[3] * q[1:(N-i-1), 1:(N-j-1)]
- qtemp[(i+2):N, (j+1):N] <- qtemp[(i+2):N, (j+1):N] + psplit[4] * q[1:(N-i-1), 1:(N-j)]
- q <- qtemp + (1-p[k])*q
- }
- q[length(q)] <- q[length(q)]+1-sum(q)
- return(q)
-}
-
-lossdist.prepay.joint <- function(dp, pp, w, S, N, defaultflag=FALSE){
- ## recursive algorithm with first order correction
- ## to compute the joint probability distribition of the loss and recovery
- ## inputs:
- ## dp: vector of default probabilities
- ## pp: vector of prepay probabilities
- ## w: vector of issuer weights
- ## S: vector of severities
- ## N: number of tick sizes on the grid
- ## defaultflag: if true computes the default
- ## outputs
- ## q: matrix of joint loss and recovery probability
- ## colSums(q) is the recovery distribution marginal
- ## rowSums(q) is the loss distribution marginal
- n <- length(dp)
- lu <- 1/(N-1)
- q <- matrix(0, N, N)
- q[1,1] <- 1
- for(k in 1:n){
- y1 <- (1-S[k]) * w[k]/lu
- y2 <- w[k]/lu
- j1 <- floor(y1)
- j2 <- floor(y2)
- if(defaultflag){
- x <- y2
- i <- j2
- }else{
- x <- y2-y1
- i <- floor(x)
- }
-
- ## weights <- c((i+1-x)*(j1+1-y1), (i+1-x)*(y1-j1), (x-i)*(y1-j1), (j1+1-y1)*(x-i))
- weights1 <- c((i+1-x)*(j1+1-y1), (i+1-x)*(y1-j1), (x-i)*(y1-j1), (j1+1-y1)*(x-i))
- dpsplit <- dp[k] * weights1
-
- 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{
- ppsplit <- pp[k] * c(j2+1-y2, y2-j2)
- }
- qtemp <- matrix(0, N, N)
- qtemp[(i+1):N,(j1+1):N] <- qtemp[(i+1):N,(j1+1):N] + dpsplit[1] * q[1:(N-i),1:(N-j1)]
- qtemp[(i+1):N,(j1+2):N] <- qtemp[(i+1):N,(j1+2):N] + dpsplit[2] * q[1:(N-i), 1:(N-j1-1)]
- qtemp[(i+2):N,(j1+2):N] <- qtemp[(i+2):N,(j1+2):N] + dpsplit[3] * q[1:(N-i-1), 1:(N-j1-1)]
- qtemp[(i+2):N,(j1+1):N] <- qtemp[(i+2):N, (j1+1):N] + dpsplit[4] * q[1:(N-i-1), 1:(N-j1)]
- if(defaultflag){
- qtemp[(i+1):N,(j2+1):N] <- qtemp[(i+1):N,(j2+1):N] + ppsplit[1] * q[1:(N-i),1:(N-j2)]
- 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{
- 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)]
- }
- q <- qtemp + (1-pp[k]-dp[k]) * q
- }
- q[length(q)] <- q[length(q)] + 1 - sum(q)
- return(q)
-}
-
-lossdistC <- function(p, w, S, N, defaultflag=FALSE){
- ## C version of lossdistrib2, roughly 50 times faster
- checkSymbol("lossdistrib")
- .C("lossdistrib", as.double(p), as.integer(length(p)),
- 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
- checkSymbol("lossdistrib_blas")
- .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){
- checkSymbol("lossdistrib_Z")
- .C("lossdistrib_Z", as.double(p), as.integer(length(p)),
- as.double(w), as.double(S), as.integer(N), as.logical(defaultflag),
- as.double(rho), as.double(Z), as.integer(length(Z)),
- q = matrix(0, N, length(Z)))$q
-}
-
-lossdistC.truncated <- function(p, w, S, N, T=N){
- ## C version of lossdistrib2, roughly 50 times faster
- checkSymbol("lossdistrib_truncated")
- .C("lossdistrib_truncated", as.double(p), as.integer(length(p)),
- as.double(w), as.double(S), as.integer(N), as.integer(T), q = double(T))$q
-}
-
-recovdistC <- function(dp, pp, w, S, N){
- ## C version of recovdist
- checkSymbol("recovdist")
- .C("recovdist", as.double(dp), as.double(pp), as.integer(length(dp)),
- as.double(w), as.double(S), as.integer(N), q = double(N))$q
-}
-
-lossdistC.joint <- function(p, w, S, N, defaultflag=FALSE){
- ## C version of lossdistrib.joint, roughly 20 times faster
- checkSymbol("lossdistrib_joint")
- .C("lossdistrib_joint", 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.jointblas <- function(p, w, S, N, defaultflag=FALSE){
- ## C version of lossdistrib.joint, roughly 20 times faster
- checkSymbol("lossdistrib_joint_blas")
- .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
- ## w is of size n.credits
- ## S is of size n.credits by nZ
- ## rho is a double
- ## Z is a vector of length nZ
- ## w is a vector if length wZ
- checkSymbol("lossdistrib_joint_Z")
- r <- .C("lossdistrib_joint_Z", as.double(dp), as.integer(length(dp)),
- as.double(w), as.double(S), as.integer(N), as.logical(defaultflag), as.double(rho),
- as.double(Z), as.double(wZ), as.integer(length(Z)), q = matrix(0, N, N))$q
-}
-
-lossdistC.prepay.joint <- function(dp, pp, w, S, N, defaultflag=FALSE){
- ## C version of lossdist.prepay.joint
- checkSymbol("lossdistrib_prepay_joint")
- r <- .C("lossdistrib_prepay_joint", as.double(dp), as.double(pp), as.integer(length(dp)),
- as.double(w), as.double(S), as.integer(N), as.logical(defaultflag), q=matrix(0, N, N))$q
- return(r)
-}
-
-lossdistC.prepay.jointZ <- function(dp, pp, w, S, N, defaultflag = FALSE, rho, Z, wZ){
- ## N is the size of the grid
- ## dp is of size n.credits
- ## pp is of size n.credits
- ## w is of size n.credits
- ## S is of size n.credits by nZ
- ## rho is a double
- ## Z is a vector of length nZ
- ## w is a vector if length wZ
- checkSymbol("lossdistrib_prepay_joint_Z")
- r <- .C("lossdistrib_prepay_joint_Z", as.double(dp), as.double(pp), as.integer(length(dp)),
- as.double(w), as.double(S), as.integer(N), as.logical(defaultflag), as.double(rho),
- as.double(Z), as.double(wZ), as.integer(length(Z)), q = matrix(0, N, N))$q
-}
-
-lossrecovdist <- function(defaultprob, prepayprob, w, S, N, defaultflag=FALSE, useC=TRUE){
- 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)
- }
- }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)
- }
- }
- return(list(L=L, R=R))
-}
-
-lossrecovdist.term <- function(defaultprob, prepayprob, w, S, N, defaultflag=FALSE, useC=TRUE){
- ## 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)){
- 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)){
- temp <- lossrecovdist(defaultprob[,t], prepayprob[,t], w, S[,t], N, defaultflag, useC)
- L[,t] <- temp$L
- R[,t] <- temp$R
- }
- }
- return(list(L=L, R=R))
-}
-
-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.joint(defaultprob[,t], w, S[,t], N, defaultflag)
- }
- }else{
- for(t in 1:ncol(defaultprob)){
- Q[t,,] <- lossdistC.prepay.joint(defaultprob[,t], prepayprob[,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)
- }
- }
- }
- gc()
- return(Q)
-}
-
-dist.transform <- function(dist.joint){
- ## compute the joint (D, R) distribution
- ## 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(i in 1:Ngrid){
- for(j in 1:Ngrid){
- index <- i+j
- 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] +
- dist.joint[t,i,j]
- }
- }
- }
- distDR[t,,] <- distDR[t,,]/sum(distDR[t,,])
- }
- return( distDR )
-}
-
-shockprob <- function(p, rho, Z, log.p=F){
- ## 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))
- }else{
- 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
- 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){
- dnorm((qnorm(p)-sqrt(rho)*Z)/sqrt(1-rho))*dqnorm(p)/sqrt(1-rho)
-}
-
-dqnorm <- function(x){
- 1/dnorm(qnorm(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
- require(distr)
- if(p0==0){
- return(0)
- }
- if(rho == 1){
- distw <- DiscreteDistribution(Z, w)
- return(pnorm(q(distw)(p0)))
- }
- eps <- 1e-12
- 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)
- phi <- 1
- while ((p-phi*dp)<0 || (p-phi*dp)>1){
- phi <- 0.8*phi
- }
- p <- p - phi*dp
- }
- return(p)
-}
-
-fit.probC <- function(Z, w, rho, p0){
- checkSymbol("fitprob")
- r <- .C("fitprob", as.double(Z), as.double(w), as.integer(length(Z)),
- as.double(rho), as.double(p0), q = double(1))
- return(r$q)
-}
-
-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){
- return(rep(R, length(Z)))
- }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))))
- }
-}
-
-stochasticrecovC <- function(R, Rtilde, Z, w, rho, porig, pmod){
- checkSymbol("stochasticrecov")
- r <- .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)))
- 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) )
-}
-
-tranche.pvvec <- function(K, L, R, cs){
- r <- rep(0, length(K)-1)
- for(i in 1:(length(K)-1)){
- r[i] <- 1/(K[i+1]-K[i]) * tranche.pv(L, R, cs, K[i], K[i+1])
- }
- return( r )
-}
-
-BClossdist <- function(defaultprob, issuerweights, recov, rho, Z, w,
- N=length(recov)+1, defaultflag=FALSE, n.int=500){
- ## do not use if weights are not gaussian, results would be incorrect
- ## since shockseverity is invalid in that case (need to use stochasticrecov)
- 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, , issuerweights, S.shocked, N)
- LZ[,i] <- temp$L
- RZ[,i] <- temp$R
- }
- L[,t] <- LZ%*%w
- R[,t] <- RZ%*%w
- }
- list(L=L, R=R)
-}
-
-BClossdistC <- function(defaultprob, issuerweights, recov, rho, Z, w,
- N=length(issuerweights)+1, defaultflag=FALSE){
- checkSymbol("BClossdist")
- L <- matrix(0, N, dim(defaultprob)[2])
- R <- matrix(0, N, dim(defaultprob)[2])
- rho <- rep(rho, length(issuerweights))
- r <- .C("BClossdist", defaultprob, dim(defaultprob)[1], dim(defaultprob)[2],
- as.double(issuerweights), as.double(recov), as.double(Z), as.double(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.pv <- function(defaultprob, issuerweights, recov, cs, K1, K2, rho1, rho2,
- Z, w, N=length(issuerweights)+1){
- ## computes the protection leg, couponleg, and bond price of a tranche
- ## in the base correlation setting
- if(K1==0){
- if(rho1!=0){
- stop("equity tranche must have 0 lower correlation")
- }
- }
- dK <- K2 - K1
- dist2 <- BClossdistC(defaultprob, issuerweights, recov, rho2, Z, w, N)
- if(rho1!=0){
- dist1 <- BClossdistC(defaultprob, issuerweights, recov, rho1, Z, w, N)
- }
- cl2 <- tranche.cl(dist2$L, dist2$R, cs, 0, K2)
- cl1 <- tranche.cl(dist1$L, dist1$R, cs, 0, K1)
- pl2 <- tranche.pl(dist2$L, cs, 0, K2)
- pl1 <- tranche.pl(dist1$L, cs, 0, K1)
- return(list(pl=(pl2-pl1)/dK, cl=(cl2-cl1)/dK,
- bp=100*(1+(pl2-pl1+cl2-cl1)/dK)))
-}
-
-BCtranche.delta <- function(portfolio, index, coupon, K1, K2, rho1, rho2, Z, w,
- N=length(portolio)+1, tradedate = Sys.Date()){
- ## computes the tranche delta (on current notional) by doing a proportional
- ## blip of all the curves
- ## if K2==1, then computes the delta using the lower attachment only
- ## this makes sense for bottom-up skews
- eps <- 1e-4
- portfolioplus <- portfolio
- portfoliominus <- portfolio
- cs <- couponSchedule(IMMDate(tradedate), index$maturity,"Q", "FIXED", coupon, 0, tradedate,
- IMMDate(tradedate, "prev"))
- for(i in 1:length(portfolio)){
- portfolioplus[[i]]@curve@hazardrates <- portfolioplus[[i]]@curve@hazardrates * (1 + eps)
- portfoliominus[[i]]@curve@hazardrates <- portfoliominus[[i]]@curve@hazardrates * (1- eps)
- }
- dPVindex <- indexpv(portfolioplus, index, tradedate = tradedate, clean = FALSE)$bp -
- indexpv(portfoliominus, index, tradedate = tradedate, clean = FALSE)$bp
- defaultprobplus <- 1 - SPmatrix(portfolioplus, length(cs$dates))
- defaultprobminus <- 1 - SPmatrix(portfoliominus, length(cs$dates))
- if(K2==1){
- K1adj <- adjust.attachments(K1, index$loss, index$factor)
- dPVtranche <- BCtranche.pv(defaultprobplus, issuerweights, recov, cs, 0, K1adj, 0, rho1,
- Z, w, N)$bp -
- BCtranche.pv(defaultprobminus, issuerweights, recov, cs, 0, K1adj, 0, rho1,
- Z, w, N)$bp
- delta <- (1 - dPVtranche/(100*dPVindex) * K1adj)/(1-K1adj)
- }else{
- Kmod <- adjust.attachments(c(K1, K2), index$loss, index$factor)
- dPVtranche <- BCtranche.pv(defaultprobplus, issuerweights, recov, cs, Kmod[1], Kmod[2], rho1, rho2,
- Z, w, N)$bp -
- BCtranche.pv(defaultprobminus, issuerweights, recov, cs, Kmod[1], Kmod[2], rho1, rho2,
- Z, w, N)$bp
- delta <- dPVtranche/(100*dPVindex)
- }
- ## dPVindex <- BCtranche.pv(portfolioplus, index, coupon, 0, 1, 0, 0.5, lu)$bp-
- ## BCtranche.pv(portfoliominus, index, coupon, 0, 1, 0, 0.5, lu)$bp
- return( delta )
-}
-
-BCstrikes <- function(portfolio, index, coupon, K, rho, N=101) {
- ## computes the strikes as a percentage of expected loss
- EL <- c()
- for(i in 2:length(K)){
- EL <- c(EL, -BCtranche.pv(portfolio, index, coupon, K[i-1], K[i], rho[i-1], rho[i], N)$pl)
- }
- Kmodified <- adjust.attachments(K, index$loss, index$factor)
- return(cumsum(EL*diff(Kmodified))/sum(EL*diff(Kmodified)))
-}
-
-delta.factor <- function(K1, K2, index){
- ## compute the factor to convert from delta on current notional to delta on original notional
- ## K1 and K2 original strikes
- factor <- (adjust.attachments(K2, index$loss, index$factor)
- -adjust.attachments(K1, index$loss, index$factor))/(K2-K1)
- return( 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))
- checkSymbol("fitprob")
- 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 )
-}
-
-MFtranche.pv <- function(cl, cs, w, rho, defaultprob, defaultprobmod, issuerweights, recov,
- Kmodified, Ngrid=length(issuerweights)+1){
- ## computes the tranches pv using the modified factor distribution
- ## p is the modified probability so that
- n.credit <- length(issuerweights)
- 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.term(pshocked, 0, issuerweights, S, Ngrid)
- return( tranche.pvvec(Kmodified, dist$L, dist$R, cs))
- }
- clusterExport(cl, list("Rstoch", "p"), envir=environment())
- result <- parSapply(cl, 1:length(w), parf)
- return( list(pv=100*(1+result%*%w), pv.w=result))
-}
|