aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--R/lossdistrib.c107
1 files changed, 53 insertions, 54 deletions
diff --git a/R/lossdistrib.c b/R/lossdistrib.c
index 6d40b711..22a449e2 100644
--- a/R/lossdistrib.c
+++ b/R/lossdistrib.c
@@ -6,7 +6,7 @@
#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);
+ 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);
@@ -16,17 +16,18 @@ void lossdistrib(double *p, int *np, double *w, double *S, int *N, int *defaultf
double shockprob(double p, double rho, double Z, int give_log);
void lossdistrib_Z(double *p, int *np, double *w, double *S, int *N, int *defaultflag,
- double *rho, double *Z, double *wZ, int *nZ, double *q);
+ double *rho, double *Z, double *wZ, int *nZ, double *q);
void lossdistrib_truncated(double *p, int *np, double *w, double *S, int *N,
- int *T, int *defaultflag, double *q);
+ int *T, int *defaultflag, double *q);
-void lossdistrib_joint( double *p, int *np, double *w, double *S, int *N, int *defaultflag, double *q);
+void lossdistrib_joint(double *p, int *np, double *w, double *S, int *N,
+ int *defaultflag, double *q);
void recovdist(double *dp, double *pp, int *n, 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);
+ double *S, int *N, int *defaultflag, double *q);
double dqnorm(double x);
double dshockprob(double p, double rho, double Z);
@@ -37,19 +38,19 @@ double shockseverity(double S, double Z, double rho, double p);
void fitprob(double* Z, double* w, int* nZ, double* rho, double* p0, double* result);
-void stochasticrecov(double* R, double* Rtilde, double* Z, double* w, int* nZ, double* rho,
- double* porig, double* pmod, double* q);
+void stochasticrecov(double* R, double* Rtilde, double* Z, double* w, int* nZ,
+ double* rho, double* porig, double* pmod, double* q);
void lossdistrib_prepay_joint_Z(double *dp, double *pp, int *ndp, double *w,
- double *S, int *N, int *defaultflag, double *rho,
- double *Z, double *wZ, int *nZ, double *q);
+ double *S, int *N, int *defaultflag, double *rho,
+ double *Z, double *wZ, int *nZ, double *q);
void lossdistrib_joint_Z(double *dp, int *ndp, double *w,
- double *S, int *N, int *defaultflag, double *rho,
- double *Z, double *wZ, int *nZ, double *q);
+ double *S, int *N, int *defaultflag, double *rho,
+ double *Z, double *wZ, int *nZ, double *q);
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 *recov, double *Z, double *w, int *n, double *rho, int *N,
+ int *defaultflag, double *L, double *R);
void lossdistrib(double *p, int *np, double *w, double *S, int *N, int *defaultflag, double *q) {
/* recursive algorithm with first order correction for computing
@@ -87,7 +88,7 @@ void lossdistrib(double *p, int *np, double *w, double *S, int *N, int *defaultf
}
/* correction for weight loss */
- if(M>*N){
+ if(M > *N){
sum = 0;
for(j=0; j<MIN(M, *N); j++){
sum += q[j];
@@ -98,7 +99,7 @@ void lossdistrib(double *p, int *np, double *w, double *S, int *N, int *defaultf
}
void lossdistrib_Z(double *p, int *np, double *w, double *S, int *N, int *defaultflag,
- double *rho, double *Z, double *wZ, int *nZ, double *q){
+ double *rho, double *Z, double *wZ, int *nZ, double *q){
int i, j;
double* pshocked = malloc(sizeof(double) * (*np) * (*nZ));
@@ -113,15 +114,17 @@ 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) {
+ int *T, int *defaultflag, double *q) {
/* recursive algorithm with first order correction for computing
the loss distribution.
+ input:
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
T where to truncate the distribution
- defaultflat if true computes the default distribution
+ defaultflag if true computes the default distribution
+ output:
q the loss distribution */
int i, j, d1, d2, M;
@@ -394,13 +397,6 @@ 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;
@@ -427,8 +423,8 @@ void fitprob(double* Z, double* w, int* nZ, double* rho, double* p0, double* res
free(q);
}
-void stochasticrecov(double* R, double* Rtilde, double* Z, double* w, int* nZ, double* rho,
- double* porig, double* pmod, double* q){
+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;
if(*porig==0){
@@ -440,14 +436,15 @@ void stochasticrecov(double* R, double* Rtilde, double* Z, double* w, int* nZ, d
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)));
+ q[i] = fabs(1 - (1 - *Rtilde) * exp( shockprob(ptilde, *rho, Z[i], 1) -
+ shockprob(*pmod, *rho, Z[i], 1)));
}
}
}
void lossdistrib_prepay_joint_Z(double *dp, double *pp, int *ndp, double *w,
- double *S, int *N, int *defaultflag, double *rho,
- double *Z, double *wZ, int *nZ, double *q) {
+ double *S, int *N, int *defaultflag, double *rho,
+ double *Z, double *wZ, int *nZ, double *q) {
int i, j;
double* dpshocked = malloc(sizeof(double) * (*ndp) * (*nZ));
double* ppshocked = malloc(sizeof(double) * (*ndp) * (*nZ));
@@ -464,8 +461,8 @@ void lossdistrib_prepay_joint_Z(double *dp, double *pp, int *ndp, double *w,
dpshocked[j + (*ndp) * i] = shockprob(dp[j], *rho, Z[i], 0);
ppshocked[j + (*ndp) * i] = shockprob(pp[j], *rho, -Z[i], 0);
}
- lossdistrib_prepay_joint(dpshocked + (*ndp) * i, ppshocked + (*ndp) * i, ndp, w,
- S + (*ndp) * i, N, defaultflag, qmat + N2 * i);
+ lossdistrib_prepay_joint(dpshocked + (*ndp) * i, ppshocked + (*ndp) * i, ndp,
+ w, S + (*ndp) * i, N, defaultflag, qmat + N2 * i);
}
dgemv_("n", &N2, nZ, &alpha, qmat, &N2, wZ, &one, &beta, q, &one);
@@ -476,8 +473,8 @@ void lossdistrib_prepay_joint_Z(double *dp, double *pp, int *ndp, double *w,
}
void lossdistrib_joint_Z(double *dp, int *ndp, double *w,
- double *S, int *N, int *defaultflag, double *rho,
- double *Z, double *wZ, int *nZ, double *q) {
+ double *S, int *N, int *defaultflag, double *rho,
+ double *Z, double *wZ, int *nZ, double *q) {
int i, j;
double* dpshocked = malloc(sizeof(double) * (*ndp) * (*nZ));
int N2 = (*N) * (*N);
@@ -493,7 +490,7 @@ void lossdistrib_joint_Z(double *dp, int *ndp, double *w,
dpshocked[j + (*ndp) * i] = shockprob(dp[j], *rho, Z[i], 0);
}
lossdistrib_joint(dpshocked + (*ndp) * i, ndp, w, S + (*ndp) * i, N,
- defaultflag, qmat + N2 * i);
+ defaultflag, qmat + N2 * i);
}
dgemv_("n", &N2, nZ, &alpha, qmat, &N2, wZ, &one, &beta, q, &one);
@@ -502,13 +499,16 @@ void lossdistrib_joint_Z(double *dp, int *ndp, double *w,
free(qmat);
}
-void BClossdist(double *defaultprob, 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 BClossdist(double *defaultprob, int *dim1, int *dim2,
+ double *issuerweights, 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
+ computes the loss and recovery distribution over time with a flat gaussian
+ correlation
inputs:
- Survprob: matrix of size dim1 x dim2. dim1 is the number of issuers and dim2 number of time steps
+ Survprob: matrix of size dim1 x dim2. dim1 is the number of issuers
+ and dim2 number of time steps
recov: vector of recoveries (length dim1)
issuerweights: vector of issuer weights (length dim2)
Z: vector of factor values (length n)
@@ -524,31 +524,30 @@ void BClossdist(double *defaultprob, int *dim1, int *dim2, double *issuerweights
double g;
double *gshocked, *Rshocked, *Sshocked, *Lw, *Rw;
int one = 1;
+ double alpha = 1;
+ double beta = 0;
gshocked = Calloc((*dim1), double);
Rshocked = Calloc((*dim1), double);
Sshocked = Calloc((*dim1), double);
- Lw = malloc((*N) * sizeof(double));
- Rw = malloc((*N) * sizeof(double));
+ Lw = malloc((*N) * (*n) * sizeof(double));
+ Rw = malloc((*N) * (*n) * sizeof(double));
for(t=0; t < (*dim2); t++) {
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);
- Rshocked[j] = 1 - Sshocked[j];
+ g = defaultprob[j + (*dim1) * t];
+ gshocked[j] = shockprob(g, *rho, Z[i], 0);
+ Sshocked[j] = shockseverity(1-recov[j], Z[i], *rho, g);
+ Rshocked[j] = 1 - Sshocked[j];
}
- /* reset Lw and Rw to 0 */
- memset(Lw, 0, *N * sizeof(double));
- 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); */
- daxpy_(N, w + i, Lw, &one, L + t * (*N), &one);
- daxpy_(N, w + i, Rw, &one, R + t * (*N), &one);
+ lossdistrib(gshocked, dim1, issuerweights, Sshocked, N, defaultflag,
+ Lw + i * (*N));
+ lossdistrib(gshocked, dim1, issuerweights, Rshocked, N, defaultflag,
+ Rw + i * (*N));
}
+ dgemv_("n", N, n, &alpha, Lw, N, w, &one, &beta, L + t * (*N), &one);
+ dgemv_("n", N, n, &alpha, Rw, N, w, &one, &beta, L + t * (*N), &one);
}
Free(gshocked);
Free(Rshocked);
generated by cgit v1.2.3-70-g09d2 (git 2.52.0) at 2026-01-10 07:14:11 +0000