summaryrefslogtreecommitdiffstats
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/lossdistrib.c119
-rw-r--r--src/lossdistrib.h20
2 files changed, 122 insertions, 17 deletions
diff --git a/src/lossdistrib.c b/src/lossdistrib.c
index 51e273b..9330cfb 100644
--- a/src/lossdistrib.c
+++ b/src/lossdistrib.c
@@ -6,15 +6,6 @@
#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
@@ -170,7 +161,7 @@ void lossdistrib_truncated(double *p, int *np, double *w, double *S, int *N,
output:
q the loss distribution */
- int i, j, d1, d2, M;
+ int i, d1, d2, M;
double lu, d, p1, p2;
double *qtemp;
int one = 1;
@@ -200,20 +191,27 @@ void lossdistrib_truncated(double *p, int *np, double *w, double *S, int *N,
free(qtemp);
}
+static inline void posK(int T, double K, double lu, double* val){
+ int i = 0;
+ for(i = 0; i < T; i++){
+ val[i] = K-lu*i;
+ }
+}
+
void exp_trunc(double *p, int *np, double *w, double *S, int *N, double *K,
double *r) {
double lu;
double *qtemp;
- double lambda;
+
lu = 1./(*N+1);
int T = (int) floor((*K) * (*N))+1;
int zero = 0;
+ int one = 1;
qtemp = calloc( T, sizeof(double));
- int i;
lossdistrib_truncated(p, np, w, S, N, &T, &zero, qtemp);
- for(i = 0; i < T; i++){
- *r += (*K - lu*i) * qtemp[i];
- }
+ double* val = malloc(T * sizeof(double));
+ posK(T, *K, lu, val);
+ *r = ddot_(&T, val, &one, qtemp, &one);
free(qtemp);
}
@@ -678,7 +676,7 @@ void stochasticrecov(double* R, double* Rtilde, double* Z, double* w, int* nZ,
}else{
ptemp = (1 - *R) / (1 - *Rtilde) * *porig;
fitprob(Z, w, nZ, rho, &ptemp, &ptilde);
-#pragma omp parallel for
+ #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)));
@@ -804,6 +802,95 @@ void BCloss_recov_dist(double *defaultprob, int *dim1, int *dim2,
free(Rw);
}
+void BCloss_recov_trunc(double *defaultprob, int *dim1, int *dim2,
+ double *issuerweights, double *recov, double *Z, double *w,
+ int *n, double *rho, int *N, double * K, int *defaultflag,
+ double *ELt, double *ERt) {
+ /*
+ computes EL_t = E[(K-L_t)^+] and ER_t = E[(K-(1-R_t))^+]
+ the the put options on loss and recovery 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
+ K: the strike
+ defaultflag: if true, computes the default distribution
+ outputs:
+ ELt: vector of length dim2
+ ERt: vector of length dim2
+ */
+ int t, i, j;
+ double g, Ktilde;
+ int one = 1;
+ int T = (int) floor((*K) * (*N))+1;
+ int Ttilde;
+ double lu = 1./(*N+1);
+ double* valL = malloc( T * sizeof(double));
+ posK(T, *K, lu, valL);
+ double* EL = malloc( (*n) * sizeof(double));
+ double* ER = malloc( (*n) * sizeof(double));
+ double* Lw = malloc(T * (*n) * sizeof(double));
+ double alpha = 1;
+ double beta = 0;
+
+ for(t=0; t < (*dim2); t++) {
+ memset(Lw, 0, T * (*n) * sizeof(double));
+ #pragma omp parallel for private(j, g, Ktilde, Ttilde)
+ for(i=0; i < *n; i++){
+ double* Rw = NULL;
+ double* Rshocked = malloc((*dim1) * sizeof(double));
+ double* Sshocked = malloc((*dim1) * sizeof(double));
+ double* gshocked = malloc((*dim1) * sizeof(double));
+ double* gshockedbar = malloc((*dim1) * sizeof(double));
+ double* valR = NULL;
+
+ for(j=0; j < (*dim1); j++){
+ g = defaultprob[j + (*dim1) * t];
+ gshocked[j] = shockprob(g, rho[j], Z[i], 0);
+ Sshocked[j] = shockseverity(1-recov[j], Z[i], rho[j], g);
+ gshockedbar[j] = 1 - gshocked[j];
+ Rshocked[j] = 1 - Sshocked[j];
+ }
+
+ lossdistrib_truncated(gshocked, dim1, issuerweights, Sshocked,
+ N, &T, defaultflag, Lw + i * T);
+ ER[i] = 0;
+ Ktilde = *K - ddot_(dim1, issuerweights, &one, Sshocked, &one);
+ if(Ktilde > 0){
+ Ttilde = (int) floor(Ktilde * (*N))+1;
+ Rw = calloc(Ttilde, sizeof(double));
+ lossdistrib_truncated(gshockedbar, dim1, issuerweights, Rshocked,
+ N, &Ttilde, defaultflag, Rw);
+ valR = malloc(Ttilde * sizeof(double));
+ posK(Ttilde, Ktilde, lu, valR);
+ ER[i] = ddot_(&Ttilde, Rw, &one, valR, &one);
+ }
+ if(Rw != NULL){
+ free(Rw);
+ }
+ if(valR != NULL){
+ free(valR);
+ }
+ free(Rshocked);
+ free(Sshocked);
+ free(gshocked);
+ free(gshockedbar);
+ }
+ dgemv_("t", &T, n, &alpha, Lw, &T, valL, &one, &beta, EL, &one);
+ ELt[t] = ddot_(n, EL, &one, w, &one);
+ ERt[t] = ddot_(n, ER, &one, w, &one);
+ }
+ free(Lw);
+ free(EL);
+ free(ER);
+ free(valL);
+}
void BCloss_dist(double *defaultprob, int *dim1, int *dim2,
double *issuerweights, double *recov, double *Z, double *w,
diff --git a/src/lossdistrib.h b/src/lossdistrib.h
index 785b4e8..d3b8f2a 100644
--- a/src/lossdistrib.h
+++ b/src/lossdistrib.h
@@ -1,3 +1,11 @@
+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 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);
@@ -9,8 +17,13 @@ void lossdistrib_Z(double *p, int *np, double *w, double *S, int *N, int *defaul
void lossdistrib_truncated(double *p, int *np, double *w, double *S, int *N,
int *T, int *defaultflag, double *q);
+static inline void posK(int T, double K, double lu, double* val);
+
+void exp_trunc(double *p, int *np, double *w, double *S, int *N, double *K,
+ double *r);
+
void lossdistrib_joint(double *p, int *np, double *w, double *S, int *N,
- int *defaultflag, double *q);
+ int *defaultflag, double *q);
void lossdistrib_joint_blas(double *p, int *np, double *w, double *S, int *N,
int *defaultflag, double *q);
@@ -43,6 +56,11 @@ 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);
+void BCloss_recov_trunc(double *defaultprob, int *dim1, int *dim2,
+ double *issuerweights, double *recov, double *Z, double *w,
+ int *n, double *rho, int *N, double * K, int *defaultflag,
+ double *ELt, double *ERt);
+
double quantile(double* Z, double* w, int nZ, double p0);
void GHquad(int *n, double* Z, double* w);
generated by cgit v1.2.3-70-g09d2 (git 2.52.0) at 2025-12-21 14:46:18 +0000