diff options
Diffstat (limited to 'R/lossdistrib.c')
| -rw-r--r-- | R/lossdistrib.c | 35 |
1 files changed, 28 insertions, 7 deletions
diff --git a/R/lossdistrib.c b/R/lossdistrib.c index 25097060..a8241b63 100644 --- a/R/lossdistrib.c +++ b/R/lossdistrib.c @@ -57,6 +57,8 @@ void BClossdist(double *SurvProb, int *dim1, int *dim2, double *issuerweights, double *recov, double *Z, double *w, int *n, double *rho, int *N,
int *defaultflag, double *L, double *R);
+double quantile(double* Z, double* w, int nZ, double p0);
+
void lossdistrib(double *p, int *np, double *w, double *S, int *N, int *defaultflag,
double *q) {
/* recursive algorithm with first order correction for computing
@@ -161,7 +163,7 @@ void lossdistrib_Z(double *p, int *np, double *w, double *S, int *N, int *defaul #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, Z[i], 0);
+ pshocked[j + (*np) * i] = shockprob(p[j], rho[j], Z[i], 0);
}
lossdistrib(pshocked + (*np) * i, np, w, S + (*np) * i, N,
defaultflag, q + (*N) * i);
@@ -584,7 +586,11 @@ void lossdistrib_prepay_joint_blas(double *dp, double *pp, int *ndp, double *w, }
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));
+ 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){
@@ -610,6 +616,19 @@ double shockseverity(double S, double Z, double rho, double p){ 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;
@@ -618,6 +637,8 @@ void fitprob(double* Z, double* w, int* nZ, double* rho, double* p0, double* res 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);
@@ -671,8 +692,8 @@ void lossdistrib_prepay_joint_Z(double *dp, double *pp, int *ndp, double *w, #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, Z[i], 0);
- ppshocked[j + (*ndp) * i] = shockprob(pp[j], *rho, -Z[i], 0);
+ 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);
@@ -700,7 +721,7 @@ void lossdistrib_joint_Z(double *dp, int *ndp, double *w, #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, Z[i], 0);
+ 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);
@@ -752,8 +773,8 @@ void BClossdist(double *defaultprob, int *dim1, int *dim2, for(i=0; i < *n; i++){
for(j=0; j < (*dim1); j++){
g = defaultprob[j + (*dim1) * t];
- gshocked[j] = shockprob(g, *rho, Z[i], 0);
- Sshocked[j] = shockseverity(1-recov[j], Z[i], *rho, g);
+ gshocked[j] = shockprob(g, rho[j], Z[i], 0);
+ Sshocked[j] = shockseverity(1-recov[j], Z[i], rho[j], g);
Rshocked[j] = 1 - Sshocked[j];
}
lossdistrib(gshocked, dim1, issuerweights, Sshocked, N, defaultflag,
|