aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rw-r--r--lossdistrib.c30
-rw-r--r--tranche_functions.R81
2 files changed, 66 insertions, 45 deletions
diff --git a/lossdistrib.c b/lossdistrib.c
index 5b44cbda..c8bd0068 100644
--- a/lossdistrib.c
+++ b/lossdistrib.c
@@ -1,15 +1,17 @@
#include <R.h>
#include <Rmath.h>
#include <string.h>
+
#define MIN(x, y) (((x) < (y)) ? (x) : (y))
-void lossdistrib(double *p, int *np, double *w, double *S, int *N, double *q) {
+void lossdistrib(double *p, int *np, double *w, double *S, int *N, int* defaultflag, double *q) {
/* recursive algorithm with first order correction for computing
the loss distribution.
p vector of default probabilities
np length of p
S vector of severities (should be same length as p)
N number of ticks in the grid
+ defaultflat if true compute the default distribution
q the loss distribution */
int i, j, d1, d2, M;
@@ -21,7 +23,7 @@ void lossdistrib(double *p, int *np, double *w, double *S, int *N, double *q) {
q[0] = 1;
M = 1;
for(i=0; i<(*np); i++){
- d = S[i] * w[i]/ lu;
+ d = *defaultflag? w[i]/lu : S[i] * w[i]/ lu;
d1 = floor(d);
d2 = ceil(d);
p1 = p[i] * (d2-d);
@@ -48,7 +50,8 @@ void lossdistrib(double *p, int *np, double *w, double *S, int *N, double *q) {
Free(qtemp);
}
-void lossdistrib_truncated(double *p, int *np, double *w, double *S, int *N, int *T, double *q) {
+void lossdistrib_truncated(double *p, int *np, double *w, double *S, int *N,
+ int *T, int *defaultflag, double *q) {
/* recursive algorithm with first order correction for computing
the loss distribution.
p vector of default probabilities
@@ -56,6 +59,7 @@ void lossdistrib_truncated(double *p, int *np, double *w, double *S, int *N, int
S vector of severities (should be same length as p)
N number of ticks in the grid
T where to truncate the distribution
+ defaultflat if true computes the default distribution
q the loss distribution */
int i, j, d1, d2, M;
@@ -68,7 +72,7 @@ void lossdistrib_truncated(double *p, int *np, double *w, double *S, int *N, int
q[0] = 1;
M = 1;
for(i=0; i<(*np); i++){
- d = S[i] * w[i] / lu;
+ d = *defaultflag? w[i] / lu : S[i] * w[i] / lu;
d1 = floor(d);
d2 = ceil(d);
p1 = p[i] * (d2-d);
@@ -90,7 +94,7 @@ void lossdistrib_truncated(double *p, int *np, double *w, double *S, int *N, int
Free(q2);
}
-void lossdistrib_joint( double *p, int *np, double *w, double *S, int *N, double *q) {
+void lossdistrib_joint( double *p, int *np, double *w, double *S, int *N, int *defaultflag, double *q) {
/* recursive algorithm with first order correction
computes jointly the loss and recovery distribution
p vector of default probabilities
@@ -98,6 +102,7 @@ void lossdistrib_joint( double *p, int *np, double *w, double *S, int *N, double
w vector of issuer weights (length np)
S vector of severities (should be same length as p)
N number of ticks in the grid
+ defaultflag if true computes the default distribution
q the joint probability distribution */
int i, j, k, m, n;
@@ -113,7 +118,7 @@ void lossdistrib_joint( double *p, int *np, double *w, double *S, int *N, double
Mx=1;
My=1;
for(k=0; k<(*np); k++){
- x = S[k] * w[k] / lu;
+ x = *defaultflag? w[k] /lu : S[k] * w[k] / lu;
y = (1-S[k]) * w[k] / lu;
i = floor(x);
j = floor(y);
@@ -212,7 +217,8 @@ void recovdist(double *dp, double *pp, int *n, double *w, double *S, int *N, dou
Free(qtemp);
}
-void lossdistrib_prepay_joint(double *dp, double *pp, int *ndp, double *w, double *S, int *N, double *q) {
+void lossdistrib_prepay_joint(double *dp, double *pp, int *ndp, double *w,
+ double *S, int *N, int *defaultflag, double *q) {
int i, j1, j2, k, m, n;
double lu, x, y1, y2, sum;
double alpha1, alpha2, beta1, beta2;
@@ -226,9 +232,9 @@ void lossdistrib_prepay_joint(double *dp, double *pp, int *ndp, double *w, doubl
Mx=1;
My=1;
for(k=0; k<(*ndp); k++){
- x = S[k] * w[k]/lu;
y1 = (1-S[k]) * w[k]/lu;
y2 = w[k]/lu;
+ x = *defaultflag? y2: y2-y1;
i = floor(x);
j1 = floor(y1);
j2 = floor(y2);
@@ -292,7 +298,8 @@ void addandmultiply(double *X, double alpha, double *Y, int n) {
}
void BClossdist(double *SurvProb, int *dim1, int *dim2, double *issuerweights,
- double *recov, double *Z, double *w, int *n, double *rho, int *N, double *L, double *R) {
+ double *recov, double *Z, double *w, int *n, double *rho, int *N,
+ int *defaultflag, double *L, double *R) {
/*
computes the loss and recovery distribution over time with a flat gaussiancorrelation
inputs:
@@ -303,6 +310,7 @@ void BClossdist(double *SurvProb, int *dim1, int *dim2, double *issuerweights,
w: vector of factor weights (length n)
rho: correlation beta
N: number of ticks in the grid
+ defaultflag: if true, computes the default distribution
outputs:
L: matrix of size (N, dim2)
R: matrix of size (N, dim2)
@@ -328,8 +336,8 @@ void BClossdist(double *SurvProb, int *dim1, int *dim2, double *issuerweights,
/* reset Lw and Rw to 0 */
memset(Lw, 0, *N * sizeof(double));
memset(Rw, 0, *N * sizeof(double));
- lossdistrib(gshocked, dim1, issuerweights, Sshocked, N, Lw);
- lossdistrib(gshocked, dim1, issuerweights, Rshocked, N, Rw);
+ 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);
}
diff --git a/tranche_functions.R b/tranche_functions.R
index 7329c971..479749a8 100644
--- a/tranche_functions.R
+++ b/tranche_functions.R
@@ -35,18 +35,23 @@ lossdistrib.fft <- function(p){
return(1/(n+1)*Re(fft(v)))
}
-lossdistrib2 <- function(p, w, S, N){
- #recursive algorithm with first order correction
- #p vector of default probabilities
- #w vector of weigths
- #S vector of severities
- #N number of ticks in the grid
+lossdistrib2 <- function(p, w, S, N, defaultflag=FALSE){
+ ## recursive algorithm with first order correction
+ ## p vector of default probabilities
+ ## w vector of weigths
+ ## S vector of severities
+ ## N number of ticks in the grid
+ ## defaultflag if true computes the default distribution
n <- length(p)
lu <- 1/(N-1)
q <- rep(0, N)
q[1] <- 1
for(i in 1:n){
- d <- S[i] * w[i] / lu
+ if(defaultflag){
+ d <- w[i] /lu
+ }else{
+ d <- S[i] * w[i] / lu
+ }
d1 <- floor(d)
d2 <- ceiling(d)
p1 <- p[i]*(d2-d)
@@ -127,7 +132,7 @@ recovdist <- function(dp, pp, w, S, N){
return(q)
}
-lossdistrib.joint <- function(p, w, S, N){
+lossdistrib.joint <- function(p, w, S, N, defaultflag=FALSE){
## recursive algorithm with first order correction
## to compute the joint probability distribution of the loss and recovery
## inputs:
@@ -135,6 +140,7 @@ lossdistrib.joint <- function(p, w, S, N){
## w: vector of issuer weights
## S: vector of severities
## N: number of tick sizes on the grid
+ ## defaultflag: if true computes the default distribution
## output:
## q: matrix of joint loss, recovery probability
## colSums(q) is the recovery distribution marginal
@@ -144,7 +150,11 @@ lossdistrib.joint <- function(p, w, S, N){
q <- matrix(0, N, N)
q[1,1] <- 1
for(k in 1:n){
- x <- S[k] * w[k]/lu
+ if(defaultflag){
+ x <- w[k] / lu
+ }else{
+ x <- S[k] * w[k]/lu
+ }
y <- (1-S[k]) * w[k]/lu
i <- floor(x)
j <- floor(y)
@@ -161,7 +171,7 @@ lossdistrib.joint <- function(p, w, S, N){
return(q)
}
-lossdistribprepay.joint <- function(dp, pp, w, S, N){
+lossdistribprepay.joint <- function(dp, pp, w, S, N, defaultflag=FALSE){
## recursive algorithm with first order correction
## to compute the joint probability distribition of the loss and recovery
## inputs:
@@ -170,6 +180,7 @@ lossdistribprepay.joint <- function(dp, pp, w, S, N){
## w: vector of issuer weights
## S: vector of severities
## N: number of tick sizes on the grid
+ ## defaultflag: if true computes the default
## outputs
## q: matrix of joint loss and recovery probability
## colSums(q) is the recovery distribution marginal
@@ -201,13 +212,13 @@ lossdistribprepay.joint <- function(dp, pp, w, S, N){
return(q)
}
-lossdistribC <- function(p, w, S, N){
+lossdistribC <- function(p, w, S, N, defaultflag){
## C version of lossdistrib2, roughly 50 times faster
if(!is.loaded("lossdistrib")){
dyn.load(paste0("lossdistrib", .Platform$dynlib.ext))
}
.C("lossdistrib", as.double(p), as.integer(length(p)),
- as.double(w), as.double(S), as.integer(N), q = double(N))$q
+ as.double(w), as.double(S), as.integer(N), as.logical(defaultflag), q = double(N))$q
}
lossdistribC.truncated <- function(p, w, S, N, T=N){
@@ -229,78 +240,78 @@ recovdistC <- function(dp, pp, w, S, N){
as.double(w), as.double(S), as.integer(N), q = double(N))$q
}
-lossdistribC.joint <- function(p, w, S, N){
+lossdistribC.joint <- function(p, w, S, N, defaultflag=FALSE){
## C version of lossdistrib.joint, roughly 20 times faster
if(!is.loaded("lossdistrib_joint")){
dyn.load(paste0("lossdistrib", .Platform$dynlib.ext))
}
.C("lossdistrib_joint", as.double(p), as.integer(length(p)), as.double(w),
- as.double(S), as.integer(N), q = matrix(0, N, N))$q
+ as.double(S), as.integer(N), as.logical(defaultflag), q = matrix(0, N, N))$q
}
-lossdistribprepayC.joint <- function(dp, pp, w, S, N){
+lossdistribprepayC.joint <- function(dp, pp, w, S, N, defaultflag=FALSE){
## C version of lossdistribprepay.joint
if(!is.loaded("lossdistrib_prepay_joint")){
dyn.load(paste0("lossdistrib", .Platform$dynlib.ext))
}
.C("lossdistrib_prepay_joint", as.double(dp), as.double(pp), as.integer(length(dp)),
- as.double(w), as.double(S), as.integer(N), q=matrix(0, N, N))$q
+ as.double(w), as.double(S), as.integer(N), as.logical(defaultflag), q=matrix(0, N, N))$q
}
-lossrecovdist <- function(defaultprob, prepayprob, w, S, N, useC=TRUE){
+lossrecovdist <- function(defaultprob, prepayprob, w, S, N, defaultflag, useC=TRUE){
if(all(!prepayprob)){
if(useC){
- L <- lossdistribC(defaultprob, w, S, N)
+ L <- lossdistribC(defaultprob, w, S, N, defaultflag)
R <- lossdistribC(defaultprob, w, 1-S, N)
}else{
- L <- lossdistrib2(defaultprob, w, S, N)
+ L <- lossdistrib2(defaultprob, w, S, N, defaultflag)
R <- lossdistrib2(defaultprob, w, 1-S, N)
}
}else{
if(useC){
- L <- lossdistribC(defaultprob, w, S, N)
+ L <- lossdistribC(defaultprob, w, S, N, defaultflag)
R <- recovdistC(defaultprob, prepayprob, w, S, N)
}else{
- L <- lossdistrib2(defaultprob, w, S, N)
+ L <- lossdistrib2(defaultprob, w, S, N, defaultflag)
R <- recovdist(defaultprob, prepayprob, w, S, N)
}
}
return(list(L=L, R=R))
}
-lossrecovdist.term <- function(defaultprob, prepayprob, w, S, N, useC=TRUE){
+lossrecovdist.term <- function(defaultprob, prepayprob, w, S, N, defaultflag=FALSE, useC=TRUE){
## computes the loss and recovery distribution over time
L <- array(0, dim=c(N, ncol(defaultprob)))
R <- array(0, dim=c(N, ncol(defaultprob)))
for(t in 1:ncol(defaultprob)){
- temp <- lossrecovdist(defaultprob[,t], prepayprob[,t], w, S[,t], N, useC)
+ temp <- lossrecovdist(defaultprob[,t], prepayprob[,t], w, S[,t], N, defaultflag, useC)
L[,t] <- temp$L
R[,t] <- temp$R
}
return(list(L=L, R=R))
}
-lossrecovdist.joint.term <- function(defaultprob, prepayprob, w, S, N, useC=TRUE){
+lossrecovdist.joint.term <- function(defaultprob, prepayprob, w, S, N, defaultflag=FALSE, useC=TRUE){
## computes the joint loss and recovery distribution over time
Q <- array(0, dim=c(ncol(defaultprob), N, N))
if(useC){
if(all(!prepayprob)){
for(t in 1:ncol(defaultprob)){
- Q[t,,] <- lossdistribC.joint(defaultprob[,t], w, S[,t], N)
+ Q[t,,] <- lossdistribC.joint(defaultprob[,t], w, S[,t], N, defaultflag)
}
}else{
for(t in 1:ncol(defaultprob)){
- Q[t,,] <- lossdistribprepayC.joint(defaultprob[,t], prepayprob[,t], w, S[,t], N)
+ Q[t,,] <- lossdistribprepayC.joint(defaultprob[,t], prepayprob[,t], w, S[,t], N, defaultflag)
}
}
}else{
if(all(!prepayprob)){
for(t in 1:ncol(defaultprob)){
- Q[t,,] <- lossdistrib.joint(defaultprob[,t], w, S[,t], N)
+ Q[t,,] <- lossdistrib.joint(defaultprob[,t], w, S[,t], N, defaultflag)
}
}else{
for(t in 1:ncol(defaultprob)){
- Q[t,,] <- lossdistribprepay.joint(defaultprob[,t], prepayprob[,t], w, S[,t], N)
+ Q[t,,] <- lossdistribprepay.joint(defaultprob[,t], prepayprob[,t], w, S[,t], N, defaultflag)
}
}
}
@@ -480,7 +491,7 @@ BClossdist <- function(SurvProb, issuerweights, recov, rho, N=length(recov)+1,
}
BClossdistC <- function(SurvProb, issuerweights, recov, rho,
- N=length(issuerweights)+1, n.int=100){
+ N=length(issuerweights)+1, n.int=100, defaultflag){
if(!is.loaded("BClossdist")){
dyn.load(paste0("lossdistrib", .Platform$dynlib.ext))
}
@@ -491,7 +502,7 @@ BClossdistC <- function(SurvProb, issuerweights, recov, rho,
R <- matrix(0, N, dim(SurvProb)[2])
r <- .C("BClossdist", SurvProb, as.integer(dim(SurvProb)[1]), as.integer(dim(SurvProb)[2]),
as.double(issuerweights), as.double(recov), as.double(Z), as.double(w),
- as.integer(n.int), as.double(rho), as.integer(N), L=L, R=R)
+ as.integer(n.int), as.double(rho), as.integer(N), as.logical(defaultflag), L=L, R=R)
return(list(L=r$L,R=r$R))
}
@@ -680,12 +691,14 @@ MFlossrecovdist2 <- function(cl, w, Z, rho, defaultprob, defaultprobmod,
ppshocked <- apply(prepayprobmod, 2, shockprob, rho=rho, Z=-Z[i])
S <- 1 - Rstoch[i,,]
dist <- lossrecovdist.joint.term(dpshocked, ppshocked, issuerweights, S, Ngrid)
+ gc()
return(dist)
}
- ## temp <- sapply(1:length(w), parf)
- clusterExport(cl, list("Rstoch", "defaultprobmod", "prepayprobmod"), envir=environment())
- temp <- parSapply(cl, 1:length(w), parf)
- clusterCall(cl, gc)
+ temp <- sapply(1:length(w), parf)
+ ## clusterExport(cl, list("Rstoch", "defaultprobmod", "prepayprobmod"), envir=environment())
+ ## temp <- parSapply(cl, 1:length(w), parf)
+ ##gc()
+ ## clusterCall(cl, gc)
Q <- array(0, dim=c(ncol(defaultprob), Ngrid, Ngrid))
for(i in 1:length(w)){
Q <- Q + w[i]*array(temp[,i], dim=c(ncol(defaultprob), Ngrid, Ngrid))