aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--calibrate_tranches.R28
-rw-r--r--cds_functions_generic.R4
-rw-r--r--lossdistrib.c132
-rw-r--r--tranche_functions.R33
4 files changed, 159 insertions, 38 deletions
diff --git a/calibrate_tranches.R b/calibrate_tranches.R
index d0a52022..6fb49934 100644
--- a/calibrate_tranches.R
+++ b/calibrate_tranches.R
@@ -90,6 +90,7 @@ for(i in 2:length(Kmodified)){
rhovec <- c(rhovec, rho)
}
+rhovec <- c(0, rhovec)
deltas <- c()
for(i in 2:5){
deltas <- c(deltas, BCtranche.delta(hy17portfolio.tweaked, hy17, 0.05, K[i-1], K[i], rhovec[i-1], rhovec[i], Ngrid))
@@ -199,6 +200,28 @@ MFlossdistrib <- function(cl, cs, w, rho, defaultprob, p, issuerweights,
return( list(L=L, R=R) )
}
+MFlossdistrib2 <- function(cl, cs, w, rho, defaultprob, p, issuerweights,
+ Kmodified, Ngrid=2*length(issuerweights)+1, useC=TRUE){
+ n.credit <- length(issuerweights)
+ Rstoch <- array(0, dim=c(n.int, n.credit, ncol(defaultprob)))
+ for(t in 1:ncol(defaultprob)){
+ for(i in 1:n.credit){
+ Rstoch[,i,t] <- stochasticrecov(recov[i], 0, Z, w, rho, defaultprob[i,t], p[i,t])
+ }
+ }
+ parf <- function(i){
+ pshocked <- apply(p, 2, shockprob, rho=rho, Z=Z[i])
+ S <- 1 - Rstoch[i,,]
+ dist <- lossrecovdist.joint.term(pshocked, 0, issuerweights, S, Ngrid, useC)
+ return(dist)
+ }
+ Q <- array(0, dim=c(ncol(defaultprob), Ngrid, Ngrid))
+ for(i in 1:length(w)){
+ Q <- Q + w[i] * parf(i)
+ }
+ return( Q )
+}
+
## computes MFdeltas
newportf <- hy17portfolio.tweaked
eps <- 1e-4
@@ -251,6 +274,7 @@ l <- matrix(0, 17, 99)
r <- matrix(0, 17, 99)
d1 <- seq(0.01, 0.99, 0.01)
Ngrid <- 5*96+1
+
MFdist <- MFlossdistrib(cl, cs, w.mod, rho, defaultprob, p, issuerweights, Kmodified, Ngrid)
for(i in 1:17){
Lfun <- splinefun(c(0, cumsum(MFdist$L[,i])),c(0, seq(0, 1, length=Ngrid)), "monoH.FC")
@@ -260,3 +284,7 @@ for(i in 1:17){
r[i, j] <- Rfun(d1[j])
}
}
+
+## joint generation of scenarios for loss and recovery distribution
+
+system.time(MFdist2 <- MFlossdistrib2(cl, cs, w.mod, rho, defaultprob, p, issuerweights, Kmodified, Ngrid, FALSE))
diff --git a/cds_functions_generic.R b/cds_functions_generic.R
index 9183b3e1..962d3c7d 100644
--- a/cds_functions_generic.R
+++ b/cds_functions_generic.R
@@ -497,7 +497,7 @@ indexpv <- function(portfolio, index, epsilon=0){
pl[i] <- defaultleg(cs, portfolio[[i]]@curve, portfolio[[i]]@recovery)
}
}
- return( list(cl=mean(cl), pl=mean(pl), bp=1+pl-cl))
+ return( list(cl=mean(cl), pl=mean(pl), bp=1+mean(cl-pl)))
}
indexduration <- function(portfolio, index){
@@ -522,7 +522,7 @@ tweakcurves <- function(portfolio, index){
## computes the tweaking factor
epsilon <- 0
f <- function(epsilon, ...){
- abs(indexpv(portfolio, index, epsilon)-index$indexref)
+ abs(indexpv(portfolio, index, epsilon)$bp-index$indexref)
}
epsilon <- optimize(f, c(-0.5, 0.5), portfolio, index, tol=1e-6)$minimum
portfolio.new <- portfolio
diff --git a/lossdistrib.c b/lossdistrib.c
index d9af0f89..3939c239 100644
--- a/lossdistrib.c
+++ b/lossdistrib.c
@@ -1,6 +1,6 @@
#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) {
@@ -12,45 +12,50 @@ void lossdistrib(double *p, int *np, double *w, double *S, int *N, double *q) {
N number of ticks in the grid
q the loss distribution */
- int i, j, d1, d2;
+ int i, j, d1, d2, M;
double lu, d, p1, p2, sum;
- double *q1, *q2;
+ /* double *q1, *q2; */
lu = 1./(*N-1);
- q1 = malloc( *N * sizeof(double));
- q2 = malloc( *N * sizeof(double));
+ /* qtemp = calloc( *N, sizeof(double)); */
+ /*q2 = calloc( *N, sizeof(double));*/
q[0] = 1;
+ M = 1;
for(i=0; i<(*np); i++){
d = S[i] * w[i]/ lu;
d1 = floor(d);
d2 = ceil(d);
p1 = p[i] * (d2-d);
p2 = p[i] - p1;
-
- for(j=0; j<d1; j++){
- q1[j]=0;
- };
- for(j=0; j<*N-d1; j++){
- q1[d1+j] = p1*q[j];
- };
- for(j=0; j<d2; j++){
- q2[j]=0;
- };
- for(j=0; j<*N-d2; j++){
- q2[d2+j] = p2*q[j];
- };
- for(j=0; j<*N; j++){
- q[j] = q1[j] + q2[j] + (1-p[i]) * q[j];
+ /*memcpy(q, qtemp, M * sizeof(double));*/
+ /* for(j=0; j < M; j++){ */
+ /* q1[j] = p1 * q[j]; */
+ /* q2[j] = p2 * q[j]; */
+ /* q[j] = (1-p[i]) * q[j]; */
+ /* } */
+ for(j=0; j < MIN(M, *N); j++){
+ if(d1==0){
+ q[d2+j] += p2 * q[j];
+ q[j] = (1-p2) * q[j];
+ }else{
+ q[d1+j] += p1 * q[j];
+ q[d2+j] += p2 * q[j];
+ q[j] = (1-p[i]) * q[j];
+ }
};
+ M += d2;
}
+
/* correction for weight loss */
- /* sum = 0; */
- /* for(j=0; j<*N; j++){ */
- /* sum += q[j]; */
- /* } */
- /* q[*N-1] += 1-sum; */
- free(q1);
- free(q2);
+ if(M>*N){
+ sum = 0;
+ for(j=0; j<MIN(M, *N); j++){
+ sum += q[j];
+ }
+ q[*N-1] += 1-sum;
+ }
+ /* free(q1); */
+ /* free(q2); */
}
void lossdistrib_truncated(double *p, int *np, double *w, double *S, int *N, int *T, double *q) {
@@ -73,7 +78,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 = S[i] * w[i] / lu;
d1 = floor(d);
d2 = ceil(d);
p1 = p[i] * (d2-d);
@@ -120,7 +125,7 @@ void lossdistrib_joint( double *p, int *np, double *w, double *S, int *N, double
double w1, w2, w3, w4;
double *qtemp;
- lu = 1/(*N-1);
+ lu = 1./(*N-1);
qtemp = malloc( (*N) * (*N) * sizeof(double));
q[0] = 1;
for(k=0; k<(*np); k++){
@@ -157,6 +162,75 @@ void lossdistrib_joint( double *p, int *np, double *w, double *S, int *N, double
free(qtemp);
}
+/* void lossdistrib_joint_fast( double *p, int *np, double *w, double *S, int *N, double *q) { */
+/* /\* recursive algorithm with first order correction */
+/* computes jointly the loss and recovery distribution */
+/* p vector of default probabilities */
+/* np length of p */
+/* w vector of issuer weights (length np) */
+/* S vector of severities (should be same length as p) */
+/* N number of ticks in the grid */
+/* q the joint probability distribution *\/ */
+
+/* int i, j, k, l; */
+/* double lu, x, y, sum; */
+/* double alpha1, alpha2, beta1, beta2; */
+/* double w1, w2, w3, w4; */
+/* double *q1, *q2, *q3, *q4; */
+/* div_t index; */
+/* int temp; */
+
+/* lu = 1./(*N-1); */
+/* q1 = calloc( (*N) * (*N), sizeof(double)); */
+/* q2 = calloc( (*N) * (*N), sizeof(double)); */
+/* q3 = calloc( (*N) * (*N), sizeof(double)); */
+/* q4 = calloc( (*N) * (*N), sizeof(double)); */
+/* q[0] = 1; */
+/* Mx=1; */
+/* My=1; */
+/* for(k=0; k<(*np); k++){ */
+/* x = S[k] * w[k] / lu; */
+/* y = (1-S[k]) * w[k] / lu; */
+/* i = floor(x); */
+/* j = floor(y); */
+/* alpha1 = i + 1 - x; */
+/* alpha2 = 1 - alpha1; */
+/* beta1 = j + 1 - y; */
+/* beta2 = 1 - beta1; */
+/* w1 = alpha1 * beta1; */
+/* w2 = alpha1 * beta2; */
+/* w3 = alpha2 * beta2; */
+/* w4 = alpha2 * beta1; */
+/* /\* init qtemp to 0 *\/ */
+/* for(l=0; l < Mx * My; l++){ */
+/* index = div(l, Mx); */
+/* temp = index.quot * (*N) + index.rem; */
+/* q1[temp] = w1 * p[k] * q[l]; */
+/* q2[temp] = w2 * p[k] * q[l]; */
+/* q3[temp] = w3 * p[k] * q[l]; */
+/* q4[temp] = w4 * p[k] * q[l]; */
+/* q[temp] = (1-p[k]) * q[temp] */
+/* } */
+/* for(l=0; l < MIN(Mx * My, (*N)*(*N)); l++){ */
+/* q[ */
+/* aux2(q, qtemp, w1 * p[k], *N, i, j); */
+/* aux2(q, qtemp, w2 * p[k], *N, i, j+1); */
+/* aux2(q, qtemp, w3 * p[k], *N, i+1, j+1); */
+/* aux2(q, qtemp, w4 * p[k], *N, i+1, j); */
+/* for(l=0; l < (*N) * (*N); l++){ */
+/* q[l] = qtemp[l] + (1-p[k]) * q[l]; */
+/* } */
+/* } */
+/* /\* correction for weight loss *\/ */
+/* sum = 0; */
+/* for(k=0; k < (*N)*(*N); k++){ */
+/* sum+=q[k]; */
+/* } */
+/* q[(*N)*(*N)-1]+= 1 - sum; */
+/* free(qtemp); */
+/* } */
+
+
inline void aux(double *q, double *qtemp, double p, int N, int d) {
/* helper function for the lossdistrib functions */
int j;
diff --git a/tranche_functions.R b/tranche_functions.R
index b553332f..89ad5e77 100644
--- a/tranche_functions.R
+++ b/tranche_functions.R
@@ -212,7 +212,7 @@ lossdistribC <- function(p, w, S, N){
lossdistribC.truncated <- function(p, w, S, N, T=N){
## C version of lossdistrib2, roughly 50 times faster
dyn.load("lossdistrib.dll")
- .C("lossdistrib_fast", as.double(p), as.integer(length(p)),
+ .C("lossdistrib_truncated", as.double(p), as.integer(length(p)),
as.double(w), as.double(S), as.integer(N), as.integer(T), q = double(T))$q
}
@@ -267,6 +267,25 @@ lossrecovdist.term <- function(defaultprob, prepayprob, w, S, N, useC=TRUE){
return(list(L=L, R=R))
}
+lossrecovdist.joint.term <- function(defaultprob, prepayprob, w, S, N, 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)
+ }
+ }
+ }else{
+ if(all(!prepayprob)){
+ for(t in 1:ncol(defaultprob)){
+ Q[t,,] <- lossdistrib.joint(defaultprob[,t], w, S[,t], N)
+ }
+ }
+ }
+ return(Q)
+}
+
shockprob <- function(p, rho, Z, log.p=F){
## computes the shocked default probability as a function of the copula factor
pnorm((qnorm(p)-sqrt(rho)*Z)/sqrt(1-rho), log.p=log.p)
@@ -391,7 +410,7 @@ tranche.pl.scenarios <- function(l, cs, K1, K2, scaled=FALSE){
}
tranche.pv <- function(L, R, cs, K1, K2, Ngrid=nrow(L)){
- return( tranche.pl(L, cs, K1, K2, Ngrid, TRUE) + tranche.cl(L, R, cs, K1, K2, Ngrid, TRUE))
+ return( tranche.pl(L, cs, K1, K2, Ngrid) + tranche.cl(L, R, cs, K1, K2, Ngrid))
}
tranche.pv.scenarios <- function(l, r, cs, K1, K2){
@@ -490,15 +509,15 @@ BCtranche.delta <- function(portfolio, index, coupon, K1, K2, rho1, rho2, N=leng
portfolioplus[[i]]@curve@hazardrates <- portfolioplus[[i]]@curve@hazardrates * (1 + eps)
portfoliominus[[i]]@curve@hazardrates <- portfoliominus[[i]]@curve@hazardrates * (1- eps)
}
- dPVindex <- indexpv(portfolioplus, index) - indexpv(portfoliominus, index)
+ dPVindex <- indexpv(portfolioplus, index)$bp - indexpv(portfoliominus, index)$bp
if(K2==1){
- dPVtranche <- BCtranche.pv(portfolioplus, index, coupon, 0, K1, 0, rho1, lu)$bp -
- BCtranche.pv(portfoliominus, index, coupon, 0, K1, 0, rho1, lu)$bp
+ dPVtranche <- BCtranche.pv(portfolioplus, index, coupon, 0, K1, 0, rho1, N)$bp -
+ BCtranche.pv(portfoliominus, index, coupon, 0, K1, 0, rho1, N)$bp
K1adj <- adjust.attachments(K1, index$loss, index$factor)
delta <- (1 - dPVtranche/(100*dPVindex) * K1adj)/(1-K1adj)
}else{
- dPVtranche <- BCtranche.pv(portfolioplus, index, coupon, K1, K2, rho1, rho2, lu)$bp -
- BCtranche.pv(portfoliominus, index, coupon, K1, K2, rho1, rho2, lu)$bp
+ dPVtranche <- BCtranche.pv(portfolioplus, index, coupon, K1, K2, rho1, rho2, N)$bp -
+ BCtranche.pv(portfoliominus, index, coupon, K1, K2, rho1, rho2, N)$bp
delta <- dPVtranche/(100*dPVindex)
}
## dPVindex <- BCtranche.pv(portfolioplus, index, coupon, 0, 1, 0, 0.5, lu)$bp-
generated by cgit v1.2.3-70-g09d2 (git 2.52.0) at 2026-01-16 21:42:40 +0000