summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--R/distrib.R80
-rw-r--r--src/lossdistrib.c119
-rw-r--r--src/lossdistrib.h20
3 files changed, 174 insertions, 45 deletions
diff --git a/R/distrib.R b/R/distrib.R
index b2c028a..536ed20 100644
--- a/R/distrib.R
+++ b/R/distrib.R
@@ -122,7 +122,7 @@ lossdistrib2.truncated <- function(p, w, S, N, cutoff=N){
## some probability mass escaping.
n <- length(p)
lu <- 1/(N-1)
- q <- rep(0, truncated)
+ q <- rep(0, cutoff)
q[1] <- 1
M <- 1
for(i in 1:n){
@@ -508,6 +508,7 @@ fit.prob <- function(Z, w, rho, p0){
}
fit.probC <- function(Z, w, rho, p0){
+ stopifnot(length(Z)==length(w))
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)
@@ -533,38 +534,61 @@ stochasticrecovC <- function(R, Rtilde, Z, w, rho, porig, pmod){
BClossdist <- function(defaultprob, issuerweights, recov, rho, Z, w,
N=length(recov)+1, defaultflag=FALSE, n.int=500){
- if(missing(Z)){
- quadrature <- GHquad(n.int)
- Z <- quadrature$Z
- w <- quadrature$w
- }
- ## 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
+
+ if(missing(Z)){
+ quadrature <- GHquad(n.int)
+ Z <- quadrature$Z
+ w <- quadrature$w
}
- L[,t] <- LZ%*%w
- R[,t] <- RZ%*%w
- }
- list(L=L, R=R)
+ stopifnot(length(Z)==length(w),
+ nrow(defaultprob) == length(issuerweights))
+
+ ## 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){
- L <- matrix(0, N, dim(defaultprob)[2])
- R <- matrix(0, N, dim(defaultprob)[2])
+ stopifnot(length(Z)==length(w),
+ nrow(defaultprob)==length(issuerweights))
+ L <- matrix(0, N, ncol(defaultprob))
+ R <- matrix(0, N, ncol(defaultprob))
rho <- rep(rho, length(issuerweights))
- r <- .C("BCloss_recov_dist", 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)
+ r <- .C("BCloss_recov_dist", defaultprob, as.integer(nrow(defaultprob)),
+ as.integer(ncol(defaultprob)), 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))
}
+
+BCER <- function(defaultprob, issuerweights, recov, K, rho, Z, w,
+ N=length(issuwerweights)+1, defaultflag=FALSE){
+ stopifnot(length(Z)==length(w),
+ nrow(defaultprob)==length(issuerweights))
+
+ rho <- rep(rho, length(issuerweights))
+ ELt <- numeric(ncol(defaultprob))
+ ERt <- numeric(ncol(defaultprob))
+ r <- .C("BCloss_recov_trunc", defaultprob, as.integer(nrow(defaultprob)),
+ as.integer(ncol(defaultprob)),
+ as.double(issuerweights), as.double(recov), as.double(Z), as.double(w),
+ as.integer(length(Z)), as.double(rho), as.integer(N), as.double(K),
+ as.logical(defaultflag), ELt=ELt, ERt=ERt)
+ return(list(ELt=r$ELt, ERt=r$ERt))
+}
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 16:57:33 +0000