aboutsummaryrefslogtreecommitdiffstats
path: root/R
diff options
context:
space:
mode:
Diffstat (limited to 'R')
-rw-r--r--R/lossdistrib.c85
-rw-r--r--R/tranche_functions.R77
2 files changed, 150 insertions, 12 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);
diff --git a/R/tranche_functions.R b/R/tranche_functions.R
index 6bea5143..306947e4 100644
--- a/R/tranche_functions.R
+++ b/R/tranche_functions.R
@@ -274,12 +274,22 @@ lossdistC.prepay.joint <- function(dp, pp, w, S, N, defaultflag=FALSE){
if(!is.loaded("lossdistrib_prepay_joint")){
dyn.load(file.path(root.dir, "code", "R", paste0("lossdistrib", .Platform$dynlib.ext)))
}
- r <- .C("lossdistrib_prepay_joint", as.double(dp), as.double(pp), as.integer(length(dp)),
- as.double(w), as.double(S), as.integer(N), as.logical(defaultflag), q=matrix(0, N, N))$q
- gc()
+ r <- .C("lossdistrib_prepay_joint_Z", as.double(dp), as.double(pp), as.integer(length(dp)),
+ as.double(w), as.double(S), as.integer(N), as.logical(defaultflag), as.double(rho),
+ as.double(Z), as.double(wZ), as.integer(length(Z)), q=matrix(0, N, N))$q
+
return(r)
}
+lossdistC.prepay.jointZ <- function(dp, pp, w, S, N, defaultflag = FALSE, rho, Z, wZ){
+ if(!is.loaded("lossdistrib_prepay_joint_Z")){
+ dyn.load(file.path(root.dir, "code", "R", paste0("lossdistrib", .Platform$dynlib.ext)))
+ }
+ r <- .C("lossdistrib_prepay_joint_Z", as.double(dp), as.double(pp), as.integer(length(dp)),
+ as.double(w), as.double(S), as.integer(N), as.logical(defaultflag), as.double(rho),
+ as.double(Z), as.double(wZ), as.integer(length(Z)), q = matrix(0, N, N))$q
+}
+
lossrecovdist <- function(defaultprob, prepayprob, w, S, N, defaultflag=FALSE, useC=TRUE){
if(missing(prepayprob)){
if(useC){
@@ -408,6 +418,15 @@ fit.prob <- function(Z, w, rho, p0){
}
}
+fit.probC <- function(Z, w, rho, p0){
+ if(!is.loaded("fitprob")){
+ dyn.load(file.path(root.dir, "code", "R", paste0("lossdistrib", .Platform$dynlib.ext)))
+ }
+ 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)
+}
+
stochasticrecov <- function(R, Rtilde, Z, w, rho, porig, pmod){
## if porig == 0 (probably matured asset) then return orginal recovery
## it shouldn't matter anyway since we never default that asset
@@ -419,6 +438,16 @@ stochasticrecov <- function(R, Rtilde, Z, w, rho, porig, pmod){
}
}
+stochasticrecovC <- function(R, Rtilde, Z, w, rho, porig, pmod){
+ if(!is.loaded("stochasticrecov")){
+ dyn.load(file.path(root.dir, "code", "R", paste0("lossdistrib", .Platform$dynlib.ext)))
+ }
+ r <- .C("stochasticrecov", as.double(R), as.double(Rtilde), as.double(Z),
+ as.double(w), as.integer(length(Z)), as.double(rho), as.double(porig),
+ as.double(pmod), q = double(length(Z)))
+ return(r$q)
+}
+
pos <- function(x){
pmax(x, 0)
}
@@ -650,6 +679,23 @@ MFupdate.prob <- function(Z, w, rho, defaultprob){
return( p )
}
+MFupdate.probC <- function(Z, w, rho, defaultprob){
+ ## update the probabilites based on a non gaussian factor
+ ## distribution so that the pv of the cds stays the same.
+ p <- matrix(0, nrow(defaultprob), ncol(defaultprob))
+ if(!is.loaded("fitprob")){
+ dyn.load(file.path(root.dir, "code", "R", paste0("lossdistrib", .Platform$dynlib.ext)))
+ }
+ for(i in 1:nrow(defaultprob)){
+ for(j in 1:ncol(defaultprob)){
+ p[i,j] <- .C("fitprob", as.double(Z), as.double(w), as.integer(length(Z)),
+ as.double(rho), as.double(defaultprob[i,j]), q = double(1))$q
+ }
+ }
+ return( p )
+}
+
+
MFlossrecovdist <- function(w, Z, rho, defaultprob, defaultprobmod, issuerweights, recov,
Ngrid=2*length(issuerweights)+1, defaultflag=FALSE){
## computes the loss and recovery distribution using the modified factor distribution
@@ -781,6 +827,31 @@ MFlossdist.prepay.joint <- function(cl, w, Z, rho, defaultprob, defaultprobmod,
return( Q )
}
+MFlossdist.prepay.joint2 <- function(w, Z, rho, defaultprob, defaultprobmod,
+ prepayprob, prepayprobmod, issuerweights, recov,
+ Ngrid=2*length(issuerweights)+1, defaultflag=FALSE){
+ ## rowSums is the loss distribution
+ ## colSums is the recovery distribution
+ ## so that recovery is the y axis and L is the x axis
+ ## if we use the persp function, losses is the axes facing us,
+ ## and R is the axis going away from us.
+ n.credit <- length(issuerweights)
+ n.int <- length(w)
+ Rstoch <- array(0, dim=c(n.credit, n.int, ncol(defaultprob)))
+
+ for(t in 1:ncol(defaultprob)){
+ for(i in 1:n.credit){
+ Rstoch[i,,t] <- stochasticrecovC(recov[i], 0, Z, w, rho, defaultprob[i,t], defaultprobmod[i,t])
+ }
+ }
+ Q <- array(0, dim=c(ncol(defaultprob), Ngrid, Ngrid))
+ for(t in 1:ncol(defaultprob)){
+ S <- 1 - Rstoch[,,t]
+ Q[t,,] <- lossdistC.prepay.jointZ(dp[,t], pp[,t], issuerweights, S, Ngrid, defaultflag, rho, Z, w)
+ }
+ return( Q )
+}
+
MFtranche.pv <- function(cl, cs, w, rho, defaultprob, defaultprobmod, issuerweights, recov,
Kmodified, Ngrid=length(issuerweights)+1){
## computes the tranches pv using the modified factor distribution
generated by cgit v1.2.3-70-g09d2 (git 2.52.0) at 2026-01-25 11:21:38 +0000