diff options
Diffstat (limited to 'R/lossdistrib.c')
| -rw-r--r-- | R/lossdistrib.c | 85 |
1 files changed, 76 insertions, 9 deletions
diff --git a/R/lossdistrib.c b/R/lossdistrib.c index 6348b2c9..c336d254 100644 --- a/R/lossdistrib.c +++ b/R/lossdistrib.c @@ -4,8 +4,10 @@ #include <omp.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 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 daxpy_(int* n, double* da, double* dx, int* incx, double* dy, int* incy);
void lossdistrib(double *p, int *np, double *w, double *S, int *N, int* defaultflag, double *q) {
/* recursive algorithm with first order correction for computing
@@ -309,18 +311,80 @@ void lossdistrib_prepay_joint(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));
+ 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){
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 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{
+ 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 addandmultiply(double *X, double alpha, double *Y, int n) {
+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;
- for(i = 0; i<n; i++){
- Y[i] += alpha*X[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)));
+ }
}
}
@@ -336,7 +400,6 @@ void lossdistrib_prepay_joint_Z(double *dp, double *pp, int *ndp, double *w, double alpha = 1;
double beta = 0;
int one = 1;
- double r;
#pragma omp parallel for private(j)
for(i = 0; i < *nZ; i++){
@@ -348,6 +411,7 @@ void lossdistrib_prepay_joint_Z(double *dp, double *pp, int *ndp, double *w, }
dgemv_("n", &N2, nZ, &alpha, qmat, &N2, wZ, &one, &beta, q, &one);
+
free(dpshocked);
free(ppshocked);
free(qmat);
@@ -375,6 +439,7 @@ void BClossdist(double *SurvProb, int *dim1, int *dim2, double *issuerweights, int t, i, j;
double g;
double *gshocked, *Rshocked, *Sshocked, *Lw, *Rw;
+ int one = 1;
gshocked = Calloc((*dim1), double);
Rshocked = Calloc((*dim1), double);
@@ -395,8 +460,10 @@ void BClossdist(double *SurvProb, int *dim1, int *dim2, double *issuerweights, 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);
+ /* addandmultiply(Lw, w[i], L + t * (*N), *N); */
+ /* addandmultiply(Rw, w[i], R + t * (*N), *N); */
+ daxpy_(N, w + i, Lw, &one, R + t * (*N), &one);
+ daxpy_(N, w + i, Rw, &one, R + t * (*N), &one);
}
}
Free(gshocked);
|