diff options
Diffstat (limited to 'lossdistrib.c')
| -rw-r--r-- | lossdistrib.c | 374 |
1 files changed, 0 insertions, 374 deletions
diff --git a/lossdistrib.c b/lossdistrib.c deleted file mode 100644 index c873df9e..00000000 --- a/lossdistrib.c +++ /dev/null @@ -1,374 +0,0 @@ -#include <R.h>
-#include <Rmath.h>
-#include <string.h>
-
-#define MIN(x, y) (((x) < (y)) ? (x) : (y))
-
-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
- S vector of severities (should be same length as p)
- N number of ticks in the grid
- defaultflat if true compute the default distribution
- q the loss distribution */
-
- int i, j, d1, d2, M;
- double lu, d, p1, p2, sum;
- double *qtemp;
-
- lu = 1./(*N-1);
- qtemp = Calloc(*N, double);
- q[0] = 1;
- M = 1;
- for(i=0; i<(*np); i++){
- d = (*defaultflag)? w[i]/lu : S[i] * w[i]/ lu;
- d1 = floor(d);
- d2 = ceil(d);
- p1 = p[i] * (d2-d);
- p2 = p[i] - p1;
- memcpy(qtemp, q, MIN(M, *N) * sizeof(double));
- 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_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.
- 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
- defaultflat if true computes the default distribution
- q the loss distribution */
-
- int i, j, d1, d2, M;
- double lu, d, p1, p2, sum;
- double *q1, *q2;
-
- lu = 1./(*N-1);
- q1 = Calloc( *T, double);
- q2 = Calloc( *T, 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), 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 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), 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 temp;
- double *qtemp;
- int Mx, My;
-
- lu = 1./(*N-1);
- qtemp = Calloc((*N) * (*N), double);
- q[0] = 1;
- Mx=1;
- My=1;
- for(k=0; k<(*ndp); k++){
- y1 = (1-S[k]) * w[k]/lu;
- y2 = w[k]/lu;
- x = (*defaultflag)? y2: y2-y1;
- i = floor(x);
- j1 = floor(y1);
- j2 = floor(y2);
- alpha1 = i + 1 - x;
- alpha2 = 1 - alpha1;
- beta1 = j1 + 1 - y1;
- beta2 = 1 - beta1;
- dpw1 = alpha1 * beta1 * dp[k];
- dpw2 = alpha1 * beta2 * dp[k];
- dpw3 = alpha2 * beta2 * dp[k];
- dpw4 = alpha2 * beta1 * dp[k];
-
- /* by default distribution, we mean names fractions of names that disappeared
- either because of default or prepayment */
- for(n=0; n<MIN(My, *N); n++){
- memcpy(qtemp+n*(*N), q+n*(*N), MIN(Mx, *N) * sizeof(double));
- }
- 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];
- temp = ppw2 * qtemp[m+(*N)*n];
- q[i+m+(*N)*(j2+1+n)] += temp;
- q[i+1+m+(*N)*(j2+1+n)] += ppw3 * qtemp[m+(*N)*n];
- q[i+1+m+(*N)*(j2+n)] += temp;
- }
- }
- }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);
-}
-
-double shockprob(double p, double rho, double Z, int give_log){
- return( pnorm( (qnorm(p, 0, 1, 1, 0) -sqrt(rho) * Z)/sqrt(1 - rho), 0, 1, 1, give_log));
-}
-
-double shockseverity(double S, double Z, double rho, double p){
- return( exp(shockprob(S * p, rho, Z, 1) - shockprob(p, rho, Z, 1)) );
-
-}
-
-void addandmultiply(double *X, double alpha, double *Y, int n) {
- int i;
- for(i = 0; i<n; i++){
- Y[i] += alpha*X[i];
- }
-}
-
-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) {
- /*
- computes the loss and recovery distribution over time with a flat gaussiancorrelation
- inputs:
- Survprob: matrix of size dim1 x dim2. dim1 is the number of issuers and dim2 number of time steps
- recov: vector of recoveries (length dim1)
- issuerweights: vector of issuer weights (length dim2)
- Z: vector of factor values (length n)
- w: vector of factor weights (length n)
- rho: correlation beta
- 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;
-
- gshocked = Calloc((*dim1), double);
- Rshocked = Calloc((*dim1), double);
- Sshocked = Calloc((*dim1), double);
- Lw = malloc((*N) * sizeof(double));
- Rw = malloc((*N) * sizeof(double));
-
- for(t=0; t < (*dim2); t++) {
- for(i=0; i < *n; i++){
- for(j=0; j < (*dim1); j++){
- g = 1 - SurvProb[j + (*dim1) * t];
- gshocked[j] = shockprob(g, *rho, Z[i], 0);
- Sshocked[j] = shockseverity(1-recov[j], Z[i], *rho, g);
- Rshocked[j] = 1 - Sshocked[j];
- }
- /* reset Lw and Rw to 0 */
- memset(Lw, 0, *N * sizeof(double));
- memset(Rw, 0, *N * sizeof(double));
- lossdistrib(gshocked, dim1, issuerweights, Sshocked, N, defaultflag, Lw);
- lossdistrib(gshocked, dim1, issuerweights, Rshocked, N, defaultflag, Rw);
- addandmultiply(Lw, w[i], L + t * (*N), *N);
- addandmultiply(Rw, w[i], R + t * (*N), *N);
- }
- }
- Free(gshocked);
- Free(Rshocked);
- Free(Sshocked);
- free(Lw);
- free(Rw);
-}
|