aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--R/build_scenarios.R20
-rw-r--r--R/lossdistrib.c35
-rw-r--r--R/script_calibrate_tranches.R5
-rw-r--r--R/tranche_functions.R142
4 files changed, 110 insertions, 92 deletions
diff --git a/R/build_scenarios.R b/R/build_scenarios.R
index 7a250a28..4f0ebee5 100644
--- a/R/build_scenarios.R
+++ b/R/build_scenarios.R
@@ -78,7 +78,6 @@ calibration <- read.table(file.path(root.dir, "Scenarios", "Calibration",
Z <- calibration$Z
w <- calibration$w
-rho <- 0.45
Ngrid <- 201
## dealnames <- list.files(file.path(root.dir, "Scenarios", paste0("Portfolios_", calibration.date)),
## pattern="*.RData")
@@ -100,9 +99,9 @@ for(deal.name in dealnames){
dp <- A$DP
pp <- A$PP
- dpmod <- MFupdate.probC(Z, w, rho, dp)
- ppmod <- MFupdate.probC(-Z, w, rho, pp)
- dist.joint <- MFlossdist.prepay.joint(w, Z, rho, dp, dpmod, pp, ppmod,
+ dpmod <- MFupdate.probC(Z, w, deal.portfolio$beta, dp)
+ ppmod <- MFupdate.probC(-Z, w, deal.portfolio$beta, pp)
+ dist.joint <- MFlossdist.prepay.joint(w, Z, deal.portfolio$beta, dp, dpmod, pp, ppmod,
deal.weights, 1-S, Ngrid)
distDR <- dist.transform(dist.joint)
@@ -158,9 +157,11 @@ for(deal.name in dealnames){
deal.datesmonthlylagged <- getdealschedule(deal.data, "1 month", workdate, 92)
cdrmonthly <- matrix(0, n.scenarios, length(deal.datesmonthly))
recoverymonthly <- matrix(0, n.scenarios, length(deal.datesmonthly))
+ scenariosrmonthly <- matrix(0, n.scenarios, length(deal.datesmonthly))
for(i in 1:n.scenarios){
cdrmonthly[i,] <- approx(deal.dates, cdr[i,], deal.datesmonthly, rule=2)$y
recoverymonthly[i,] <- approx(deal.dates, intexrecov[i,], deal.datesmonthly, rule=2)$y
+ scenariosrmonthly[i,] <- approx(deal.dates, intexrecov[i,], deal.datesmonthly, rule=2)$y
}
## compute reinvestment price
@@ -169,6 +170,17 @@ for(deal.name in dealnames){
forwards <- DC$forwards
reinvprices <- compute.reinvprices(df, forwards, cdrmonthly, recoverymonthly, 0.025, 0.07, 84, 3)
+ ## prev <- scenariosrmonthly[,1]
+ ## loanprices <- c()
+ ## bondprices <- c()
+ ## for(t in 2:dim(cdrmonthly)[2]){
+ ## w <- scenariosrmonthly[,t]-scenariosrmonthly[,t-1]+
+ ## prev*cdrmonthly[,t]/100 * intexrecov[,t] + yearFrac(deal.datesmonthly[t-1], deal.datesmonthly[t])
+ ## prev <- w
+ ## loanprices <- c(loanprices, crossprod(w, reinvprices$loan)/sum(w))
+ ## bondprices <- c(bondprices, crossprod(w, reinvprices$bond)/sum(w))
+ ## }
+ ## browser()
save.dir <- file.path(root.dir, "Scenarios", paste("Intex curves", workdate, sep="_"), "csv")
if(!file.exists(save.dir)){
dir.create(save.dir, recursive = T)
diff --git a/R/lossdistrib.c b/R/lossdistrib.c
index 25097060..a8241b63 100644
--- a/R/lossdistrib.c
+++ b/R/lossdistrib.c
@@ -57,6 +57,8 @@ 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);
+double quantile(double* Z, double* w, int nZ, double p0);
+
void lossdistrib(double *p, int *np, double *w, double *S, int *N, int *defaultflag,
double *q) {
/* recursive algorithm with first order correction for computing
@@ -161,7 +163,7 @@ void lossdistrib_Z(double *p, int *np, double *w, double *S, int *N, int *defaul
#pragma omp parallel for private(j)
for(i = 0; i < *nZ; i++){
for(j = 0; j < *np; j++){
- pshocked[j + (*np) * i] = shockprob(p[j], *rho, Z[i], 0);
+ pshocked[j + (*np) * i] = shockprob(p[j], rho[j], Z[i], 0);
}
lossdistrib(pshocked + (*np) * i, np, w, S + (*np) * i, N,
defaultflag, q + (*N) * i);
@@ -584,7 +586,11 @@ void lossdistrib_prepay_joint_blas(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));
+ if(rho==1){
+ return((double)(Z<=qnorm(p, 0, 1, 1, 0)));
+ }else{
+ return( pnorm( (qnorm(p, 0, 1, 1, 0) - sqrt(rho) * Z)/sqrt(1 - rho), 0, 1, 1, give_log));
+ }
}
double dqnorm(double x){
@@ -610,6 +616,19 @@ double shockseverity(double S, double Z, double rho, double p){
return( exp(shockprob(S * p, rho, Z, 1) - shockprob(p, rho, Z, 1)) );
}
+double quantile(double* Z, double* w, int nZ, double p0){
+ double cumw;
+ int i;
+ cumw = 0;
+ for(i=0; i<nZ; i++){
+ cumw += w[i];
+ if(cumw > p0){
+ break;
+ }
+ }
+ return( Z[i] );
+}
+
void fitprob(double* Z, double* w, int* nZ, double* rho, double* p0, double* result){
double eps = 1e-12;
int one = 1;
@@ -618,6 +637,8 @@ void fitprob(double* Z, double* w, int* nZ, double* rho, double* p0, double* res
if(*p0 == 0){
*result = 0;
+ }else if(*rho == 1){
+ *result = pnorm(quantile(Z, w, *nZ, *p0), 0, 1, 1, 0);
}else{
shockprobvec2(*p0, *rho, Z, *nZ, q);
dp = (ddot_(nZ, q, &one, w, &one) - *p0)/ddot_(nZ, q + *nZ, &one, w, &one);
@@ -671,8 +692,8 @@ void lossdistrib_prepay_joint_Z(double *dp, double *pp, int *ndp, double *w,
#pragma omp parallel for private(j)
for(i = 0; i < *nZ; i++){
for(j = 0; j < *ndp; j++){
- dpshocked[j + (*ndp) * i] = shockprob(dp[j], *rho, Z[i], 0);
- ppshocked[j + (*ndp) * i] = shockprob(pp[j], *rho, -Z[i], 0);
+ dpshocked[j + (*ndp) * i] = shockprob(dp[j], rho[j], Z[i], 0);
+ ppshocked[j + (*ndp) * i] = shockprob(pp[j], rho[j], -Z[i], 0);
}
lossdistrib_prepay_joint_blas(dpshocked + (*ndp) * i, ppshocked + (*ndp) * i, ndp,
w, S + (*ndp) * i, N, defaultflag, qmat + N2 * i);
@@ -700,7 +721,7 @@ void lossdistrib_joint_Z(double *dp, int *ndp, double *w,
#pragma omp parallel for private(j)
for(i = 0; i < *nZ; i++){
for(j = 0; j < *ndp; j++){
- dpshocked[j + (*ndp) * i] = shockprob(dp[j], *rho, Z[i], 0);
+ dpshocked[j + (*ndp) * i] = shockprob(dp[j], rho[j], Z[i], 0);
}
lossdistrib_joint_blas(dpshocked + (*ndp) * i, ndp, w, S + (*ndp) * i, N,
defaultflag, qmat + N2 * i);
@@ -752,8 +773,8 @@ void BClossdist(double *defaultprob, int *dim1, int *dim2,
for(i=0; i < *n; i++){
for(j=0; j < (*dim1); j++){
g = defaultprob[j + (*dim1) * t];
- gshocked[j] = shockprob(g, *rho, Z[i], 0);
- Sshocked[j] = shockseverity(1-recov[j], Z[i], *rho, g);
+ gshocked[j] = shockprob(g, rho[j], Z[i], 0);
+ Sshocked[j] = shockseverity(1-recov[j], Z[i], rho[j], g);
Rshocked[j] = 1 - Sshocked[j];
}
lossdistrib(gshocked, dim1, issuerweights, Sshocked, N, defaultflag,
diff --git a/R/script_calibrate_tranches.R b/R/script_calibrate_tranches.R
index 6ad98263..ab86218b 100644
--- a/R/script_calibrate_tranches.R
+++ b/R/script_calibrate_tranches.R
@@ -1,5 +1,6 @@
#!/usr/bin/Rscript
require(methods)
+options(warn=2)
args <- commandArgs(trailingOnly=TRUE)
if(.Platform$OS.type == "unix"){
@@ -85,7 +86,7 @@ Z <- quadrature$nodes
w.mod <- w
defaultprob <- 1 - SurvProb
p <- defaultprob
-rho <- 0.45
+rho <- rep(0.45, n.credit)
result <- matrix(0, 4, n.int)
for(l in 1:100){
@@ -111,7 +112,7 @@ for(l in 1:100){
err <- 0
for(i in 1:n.credit){
for(j in 1:ncol(p)){
- err <- err + abs(crossprod(shockprob(p[i,j], rho, Z), program$weight) - defaultprob[i,j])
+ err <- err + abs(crossprod(shockprob(p[i,j], rho[i], Z), program$weight) - defaultprob[i,j])
}
}
errvec <- c(errvec, err)
diff --git a/R/tranche_functions.R b/R/tranche_functions.R
index f08dd6c9..2b7a6219 100644
--- a/R/tranche_functions.R
+++ b/R/tranche_functions.R
@@ -451,7 +451,11 @@ dist.transform <- function(dist.joint){
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)
+ if(rho==1){
+ return(Z<=qnorm(p))
+ }else{
+ pnorm((qnorm(p)-sqrt(rho)*Z)/sqrt(1-rho), log.p=log.p)
+ }
}
shockseverity <- function(S, Stilde=1, Z, rho, p){
@@ -470,22 +474,26 @@ dqnorm <- function(x){
fit.prob <- function(Z, w, rho, p0){
## if the weights are not perfectly gaussian, find the probability p such
## E_w(shockprob(p, rho, Z)) = p0
+ require(distr)
if(p0==0){
return(0)
- }else{
- eps <- 1e-12
- dp <- (crossprod(shockprob(p0,rho,Z),w)-p0)/crossprod(dshockprob(p0,rho,Z),w)
- p <- p0
- while(abs(dp) > eps){
- dp <- (crossprod(shockprob(p,rho,Z),w)-p0)/crossprod(dshockprob(p,rho,Z),w)
- phi <- 1
- while ((p-phi*dp)<0 || (p-phi*dp)>1){
- phi <- 0.8*phi
- }
- p <- p - phi*dp
+ }
+ if(rho == 1){
+ distw <- DiscreteDistribution(Z, w)
+ return(pnorm(q(distw)(p0)))
+ }
+ eps <- 1e-12
+ dp <- (crossprod(shockprob(p0,rho,Z),w)-p0)/crossprod(dshockprob(p0,rho,Z),w)
+ p <- p0
+ while(abs(dp) > eps){
+ dp <- (crossprod(shockprob(p,rho,Z),w)-p0)/crossprod(dshockprob(p,rho,Z),w)
+ phi <- 1
+ while ((p-phi*dp)<0 || (p-phi*dp)>1){
+ phi <- 0.8*phi
}
- return(p)
+ p <- p - phi*dp
}
+ return(p)
}
fit.probC <- function(Z, w, rho, p0){
@@ -648,27 +656,29 @@ tranche.pvvec <- function(K, L, R, cs){
return( r )
}
-BClossdist <- function(defaultprob, issuerweights, recov, rho, N=length(recov)+1,
- n.int=100){
+BClossdist <- function(defaultprob, issuerweights, recov, rho, Z, w,
+ N=length(recov)+1, defaultflag=FALSE){
+ if(missing(Z)){
quadrature <- gauss.quad.prob(n.int, "normal")
Z <- quadrature$nodes
w <- quadrature$weights
- LZ <- matrix(0, N, n.int)
- RZ <- matrix(0, N, n.int)
- 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, 0, issuerweights, S.shocked, N)
- LZ[,i] <- temp$L
- RZ[,i] <- temp$R
- }
- L[,t] <- LZ%*%w
- R[,t] <- RZ%*%w
+ }
+ LZ <- matrix(0, N, n.int)
+ RZ <- matrix(0, N, n.int)
+ 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, 0, issuerweights, S.shocked, N)
+ LZ[,i] <- temp$L
+ RZ[,i] <- temp$R
}
- list(L=L, R=R)
+ L[,t] <- LZ%*%w
+ R[,t] <- RZ%*%w
+ }
+ list(L=L, R=R)
}
BClossdistC <- function(defaultprob, issuerweights, recov, rho, Z, w,
@@ -764,7 +774,7 @@ MFupdate.prob <- function(Z, w, rho, defaultprob){
p <- matrix(0, nrow(defaultprob), ncol(defaultprob))
for(i in 1:nrow(defaultprob)){
for(j in 1:ncol(defaultprob)){
- p[i,j] <- fit.prob(Z, w, rho, defaultprob[i,j])
+ p[i,j] <- fit.prob(Z, w, rho[i], defaultprob[i,j])
}
}
return( p )
@@ -782,39 +792,12 @@ MFupdate.probC <- function(Z, w, rho, defaultprob){
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
+ as.double(rho[i]), 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
- n.credit <- length(issuerweights)
- n.int <- length(w)
- 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], defaultprobmod[i,t])
- }
- }
- parf <- function(i){
- pshocked <- apply(defaultprobmod, 2, shockprob, rho=rho, Z=Z[i])
- S <- 1 - Rstoch[i,,]
- dist <- lossrecovdist.term(pshocked, 0, issuerweights, S, Ngrid, defaultflag)
- }
- L <- matrix(0, Ngrid, ncol(defaultprob))
- R <- matrix(0, Ngrid, ncol(defaultprob))
- for(i in 1:length(w)){
- dist <- parf(i)
- L <- L + dist$L * w[i]
- R <- R + dist$R * w[i]
- }
- return( list(L=L, R=R) )
-}
-
MFlossrecovdist.prepay <- function(w, Z, rho, defaultprob, defaultprobmod, prepayprob, prepayprobmod,
issuerweights, recov, Ngrid=2*length(issuerweights)+1, defaultflag=FALSE){
## computes the loss and recovery distribution using the modified factor distribution
@@ -875,30 +858,31 @@ MFlossdist.joint <- function(cl, w, Z, rho, defaultprob, defaultprobmod, issuerw
}
MFlossdist.prepay.joint <- 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)))
+ 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])
- }
+ for(t in 1:ncol(defaultprob)){
+ for(i in 1:n.credit){
+ Rstoch[i,,t] <- stochasticrecovC(recov[i], 0, Z, w, rho[i],
+ 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(defaultprobmod[,t], prepayprobmod[,t], issuerweights,
- S, Ngrid, defaultflag, rho, Z, w)
+ Q <- array(0, dim=c(ncol(defaultprob), Ngrid, Ngrid))
+ for(t in 1:ncol(defaultprob)){
+ S <- 1 - Rstoch[,,t]
+ Q[t,,] <- lossdistC.prepay.jointZ(defaultprobmod[,t], prepayprobmod[,t], issuerweights,
+ S, Ngrid, defaultflag, rho, Z, w)
}
- return( Q )
+ return( Q )
}
MFtranche.pv <- function(cl, cs, w, rho, defaultprob, defaultprobmod, issuerweights, recov,
generated by cgit v1.2.3-70-g09d2 (git 2.52.0) at 2026-01-06 03:20:28 +0000