diff options
| -rw-r--r-- | src/lossdistrib.cpp | 796 | ||||
| -rw-r--r-- | src/lossdistrib.hpp | 75 |
2 files changed, 871 insertions, 0 deletions
diff --git a/src/lossdistrib.cpp b/src/lossdistrib.cpp new file mode 100644 index 0000000..27841be --- /dev/null +++ b/src/lossdistrib.cpp @@ -0,0 +1,796 @@ +#include <omp.h> +#include "lossdistrib.hpp" +#include <Rcpp.h> +#include <cstring> + +#define USE_BLAS + +using namespace Rcpp; + +// [[Rcpp::plugins(cpp11, openmp)]] + +// [[Rcpp::export]] +List GHquadCpp(int n) { + // Setup for eigenvalue computations + char JOBZ = 'V'; // Compute eigenvalues & vectors + int INFO; + // Initialize array for workspace + double * WORK = new double[2*n-2]; + + // Initialize array for eigenvectors + double * V = new double[n*n]; + + std::vector<double> w(n); + std::vector<double> Z(n); + for(int i = 0; i<n-1; i++) { + w[i] = sqrt((i+1.)/2); + } + + // Run eigen decomposition + dstev_(&JOBZ, &n, Z.data(), w.data(), V, &n, WORK, &INFO); + + for (int i=0; i<n; i++) { + w[i] = V[i*n] * V[i*n]; + Z[i] *= sqrt(2); + } + + // Deallocate temporary arrays + delete[] WORK; + delete[] V; + return List::create(Named("Z") = Z, + Named("w") = w); +} + +void lossdistrib(const std::vector<double>& p, const std::vector<double>& w, + const double S[], const int N, const int T, double q[], + bool defaultflag = false) { + /* recursive algorithm with first order correction for computing + the loss distribution. + p: vector of default probabilities + w: issuer weights + S: vector of severities (should be same length as p) + N: number of ticks in the grid + defaultflag: if true compute the default distribution */ + + int d1, d2; + double d, p1, p2; + double* qtemp = new double[T]; + q[0] = 1; + int bound; + double pbar; + int one = 1; + openblas_set_num_threads(1); + double lu = 1./(N-1); + + int M = 1; + for(size_t i = 0; i < p.size(); i++){ + d = defaultflag? w[i]/lu : S[i] * w[i]/ lu; + d1 = floor(d); + d2 = ceil(d); + p1 = p[i] * (d2-d); + p2 = p[i] - p1; + std::memcpy(qtemp, q, std::min(M, T) * sizeof(double)); + pbar = 1-p[i]; + bound = std::min(M, T); +#ifdef USE_BLAS + dscal_(&bound, &pbar, q, &one); +#else + for(int j=0; j < bound; j++) { + q[j] *= pbar * q[j]; + } +#endif + bound = std::min(M, T-d2); +#ifdef USE_BLAS + daxpy_(&bound, &p1, qtemp, &one, q+d1, &one); + daxpy_(&bound, &p2, qtemp, &one, q+d2, &one); +#else + for(int j=0; j < bound; j++) { + q[d1+j] += p1 * qtemp[j]; + q[d2+j] += p2 * qtemp[j]; + } +#endif + M += d2; + } + /* correction for weight loss */ + if(M > N && T==N){ + double sum = 0; + for(int j=0; j < std::min(M, N); j++){ + sum += q[j]; + } + q[N-1] += 1-sum; + } + delete[] qtemp; +} + +// [[Rcpp::export]] +inline NumericVector lossdistrib(const NumericVector& p, const NumericVector& w, + const NumericVector& S, const int N, + bool defaultflag = false) { + double* q = new double[N]; + lossdistrib(as<std::vector<double>>(p), as<std::vector<double>>(w), + S.begin(), N, N, q, defaultflag); + NumericVector r(N, q); + delete[] q; + return r; +} + +inline std::vector<double> lossdistrib(const NumericVector& p, const NumericVector& w, + const NumericVector& S, const int N, const int T, + bool defaultflag = false) { + std::vector<double> q(T); + lossdistrib(as<std::vector<double>>(p), as<std::vector<double>>(w), + S.begin(), N, T, q.data(), defaultflag); + return q; +} + +// [[Rcpp::export]] +NumericMatrix lossdistrib_Z(const std::vector<double>& p, const std::vector<double>& w, + const NumericMatrix S, int N, const std::vector<double>& rho, + const std::vector<double>& Z, bool defaultflag=false){ + + double* q = new double[N*Z.size()]; + int p_size = p.size(); + #pragma omp parallel for + for(size_t i = 0; i < Z.size(); i++) { + std::vector<double> pshocked(p.size()); + for(size_t j = 0; j < p_size; j++){ + pshocked[j] = shockprob(p[j], rho[j], Z[i], 0); + } + lossdistrib(pshocked, w, S(_,i).begin(), N, N, q+N*i, defaultflag); + } + NumericMatrix qmat(N, Z.size(), q); + delete[] q; + return qmat; +} + +static inline double* posK(int T, double K, double lu){ + double* val = new double[T]; + for(int i = 0; i < T; i++){ + val[i] = K-lu*i; + } + return val; +} + +// [[Rcpp::export]] +double exp_trunc(const NumericVector& p, const NumericVector& w, const NumericVector& S, int N, double K) { + double lu = 1./(N+1); + int T = (int) floor(K * N)+1; + int one = 1; + double* q = new double[T]; + lossdistrib(as<std::vector<double>>(p), as<std::vector<double>>(w), + S.begin(), N, T, q); + double* val = posK(T, K, lu); + int r = ddot_(&T, val, &one, q, &one); + delete[] q; + delete[] val; + return r; +} + +void lossdistrib_joint(const std::vector<double>& p, const std::vector<double>& w, + const double S[], int N, matrix& q, bool defaultflag=false) { + /* recursive algorithm with first order correction + computes jointly the loss and recovery distribution + p vector of default probabilities + w vector of issuer weights (same length as p) + S vector of severities (should be same length as p) + N number of ticks in the grid + defaultflag if true computes the default distribution + returns the joint probability distribution */ + + int i, j; + double x, y; + double alpha1, alpha2, beta1, beta2; + double w1, w2, w3, w4; + matrix qtemp(N, N); + double lu = 1./(N-1); + q(0,0) = 1; + int Mx=1, My=1; + int one = 1; + for(size_t k = 0; k < p.size(); k++) { + x = defaultflag? w[k] /lu : 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 * p[k]; + w2 = alpha1 * beta2 * p[k]; + w3 = alpha2 * beta2 * p[k]; + w4 = alpha2 * beta1 * p[k]; + + for(int n=0; n<std::min(My, N); n++) { + std::memcpy(qtemp(n), q(n), std::min(Mx, N) * sizeof(double)); + } + int bound = std::min(Mx, N); + double pbar = 1-p[k]; + for(int n=0; n < std::min(My, N); n++) { +#ifdef USE_BLAS + dscal_(&bound, &pbar, q(n), &one); +#else + for(int m=0; m < bound; m++){ + q(m, n) *= pbar; + } +#endif + } + bound = std::min(Mx, N-i-1); + for(int n=0; n < std::min(My, N-j-1); n++) { +#ifdef USE_BLAS + daxpy_(&bound, &w1, qtemp(n), &one, &q(i,j+n), &one); + daxpy_(&bound, &w2, qtemp(n), &one, &q(i,j+1+n), &one); + daxpy_(&bound, &w3, qtemp(n), &one, &q(i+1,j+1+n), &one); + daxpy_(&bound, &w4, qtemp(n), &one, &q(i+1,j+n), &one); +#else + double temp; + for(int m = 0; m < bound; m++) { + temp = qtemp(m,n); + q(i+m,j+n) += w1 * temp; + q(i+m,j+1+n) += w2 * temp; + q(i+1+m,j+1+n] += w3 * temp; + q(i+1+m,j+n) += w4 * temp; + } +#endif + } + Mx += i+1; + My += j+1; + } + /* correction for weight loss */ + if(Mx>N || My>N){ + double sum = 0; + for(int n=0; n < std::min(My, N); n++){ + for(int m=0; n < std::min(Mx, N); m++){ + sum += q(m,n); + } + } + q[std::min(N, Mx) * std::min(N, My)-1] += 1 - sum; + } +} + +// [[Rcpp::export]] +NumericMatrix lossdistrib_joint(const std::vector<double>& p, const std::vector<double>& w, + const NumericVector S, int N, bool defaultflag=false) { + matrix q(N,N); + lossdistrib_joint(p, w, S.begin(), N, q, defaultflag); + NumericMatrix M(N, N, q.data()); + return M; +} + +// [[Rcpp::export]] +NumericVector recovdist(const NumericVector& dp, const NumericVector& pp, + const NumericVector& w, const NumericVector& S, int N) { + /* recursive algorithm with first order correction for computing + the recovery distribution in case of prepayment. + dp vector of default probabilities + pp vector of prepay probabilities + S vector of severities (should be same length as p) + w vector of weights + N number of ticks in the grid + returns the loss distribution */ + + int d1l, d1u, d2l, d2u; + double d1, d2, dp1, dp2, pp1, pp2; + + double lu = 1./(N - 1); + NumericVector qtemp(N); + NumericVector q(N); + q[0] = 1; + int M = 1; + for(size_t i=0; i<dp.size(); i++) { + d1 = w[i] * (1-S[i]) /lu; + d2 = w[i]/lu; + d1l = floor(d1); + d1u = d1l + 1; + d2l = floor(d2); + d2u = d2l + 1; + dp1 = dp[i] * (d1u - d1); + dp2 = dp[i] - dp1; + pp1 = pp[i] * (d2u - d2); + pp2 = pp[i] - pp1; + std::copy_n(q.begin(), std::min(M, N), qtemp.begin()); + for(int j = 0; j< std::min(M, N); j++) { + q[j] = (1-dp[i]-pp[i]) * q[j]; + } + for(int j=0; j < std::min(M, N-d2u); j++) { + q[d1l+j] += dp1 * qtemp[j]; + q[d1u+j] += dp2 * qtemp[j]; + q[d2l+j] += pp1 * qtemp[j]; + q[d2u+j] += pp2 * qtemp[j]; + }; + M += d2u; + } + /* correction for weight loss */ + if(M > N) { + double sum = 0; + for(int j=0; j < std::min(M, N); j++) { + sum += q[j]; + } + q[N-1] += 1-sum; + } + return q; +} + +NumericVector lossdistrib_prepay_joint(const NumericVector& dp, const NumericVector& pp, + const NumericVector& w, const NumericVector& S, + int N, bool defaultflag=false) { + int i, j1, j2; + double x, y1, y2; + double alpha1, alpha2, beta1, beta2; + double dpw1, dpw2, dpw3, dpw4; + double ppw1, ppw2, ppw3; + + double lu = 1./(N-1); + double* qtemp = new double[N*N]; + double* q = new double[N*N]; + q[0] = 1; + + int Mx = 1; + int My = 1; + for(size_t k = 0; k < dp.size(); k++) { + 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); + alpha1 = i + 1 - x; + alpha2 = 1 - alpha1; + beta1 = j1 + 1 - y1; + beta2 = 1 - beta1; + dpw1 = alpha1 * beta1 * dp[k]; + dpw2 = alpha1 * beta2 * dp[k]; + dpw3 = alpha2 * beta2 * dp[k]; + dpw4 = alpha2 * beta1 * dp[k]; + + /* by default distribution, we mean names fractions of names that disappeared + either because of default or prepayment */ + for(int n=0; n<std::min(My, N); n++){ + std::memcpy(qtemp+n*N, q+n*N, std::min(Mx, N) * sizeof(double)); + } + int bound = std::min(Mx, N); + double pbar = 1-dp[k]-pp[k]; + int one; + for(int n = 0; n < std::min(My, N); n++) { +#ifdef USE_BLAS + dscal_(&bound, &pbar, q+N*n, &one); +#else + for(int m=0; m < bound; m++) { + q[m+N*n] *= pbar; + } +#endif + } + double temp; + double* begin; + if(defaultflag) { + ppw1 = alpha1 * alpha1 * pp[k]; + ppw2 = alpha1 * alpha2 * pp[k]; + ppw3 = alpha2 * alpha2 * pp[k]; + bound = std::min(Mx, N-i-1); + for(int n=0; n < std::min(My, N-j2-1); n++) { + begin = qtemp+N*n; +#ifdef USE_BLAS + + daxpy_(&bound, &dpw1, begin, &one, q+i+N*(j1+n), &one); + daxpy_(&bound, &dpw2, begin, &one, q+i+N*(j1+1+n), &one); + daxpy_(&bound, &dpw3, begin, &one, q+i+1+N*(j1+1+n), &one); + daxpy_(&bound, &dpw4, begin, &one, q+i+1+N*(j1+n), &one); + + daxpy_(&bound, &ppw1, begin, &one, q+i+N*(j2+n), &one); + daxpy_(&bound, &ppw2, begin, &one, q+i+N*(j2+1+n), &one); + daxpy_(&bound, &ppw3, begin, &one, q+i+1+N*(j2+1+n), &one); + daxpy_(&bound, &ppw2, begin, &one, q+i+1+N*(j2+n), &one); +#else + for(int m=0; m < bound; m++) { + temp = *(begin + m); + q[i+m+N*(j1+n)] += dpw1 * temp; + q[i+1+m+N*(j1+n)] += dpw4 * temp; + q[i+m+N*(j1+1+n)] += dpw2 * temp; + q[i+1+m+N*(j1+1+n)] += dpw3 * temp; + + q[i+m+N*(j2+n)] += ppw1 * temp; + q[i+1+m+N*(j2+n)] += ppw2 * temp; + q[i+m+N*(j2+1+n)] += ppw2 * temp; + q[i+1+m+N*(j2+1+n)] += ppw3 * temp; + + } +#endif + } + } else { + ppw1 = pp[k] * (j2+1-y2); + ppw2 = pp[k] * (y2-j2); + bound = std::min(Mx, N-i-1); + for(int n = 0; n < std::min(My, N-j2-1); n++) { + begin = qtemp + N*n; +#ifdef USE_BLAS + + daxpy_(&bound, &dpw1, begin, &one, q+i+N*(j1+n), &one); + daxpy_(&bound, &dpw2, begin, &one, q+i+N*(j1+1+n), &one); + daxpy_(&bound, &dpw3, begin, &one, q+i+1+N*(j1+1+n), &one); + daxpy_(&bound, &dpw4, begin, &one, q+i+1+N*(j1+n), &one); + daxpy_(&bound, &ppw1, begin, &one, q+N*(j2+n), &one); + daxpy_(&bound, &ppw2, begin, &one, q+N*(j2+1+n), &one); +#else + for(int m = 0; m < bound; m++) { + temp = *(begin+m); + + q[i+m+N*(j1+n)] += dpw1 * temp; + q[i+1+m+N*(j1+n)] += dpw4 * temp; + q[i+m+N*(j1+1+n)] += dpw2 * temp; + q[i+1+m+N*(j1+1+n)] += dpw3 * temp; + + q[m+N*(j2+n)] += ppw1 * temp; + q[m+N*(j2+1+n)] += ppw2 * temp; + } +#endif + } + } + Mx += i + 1; + My += j2 + 1; + } + /* correction for weight loss */ + if(Mx > N || My > N){ + double sum = 0; + for(int m = 0; m < std::min(Mx, N); m++) { + for(int n=0; n < std::min(My, N); n++) { + sum += q[m+n*N]; + } + } + q[std::min(N, Mx)*std::min(My,N)-1] += 1 - sum; + } + free(qtemp); +} + +double shockprob(double p, double rho, double Z, bool give_log=false) { + if(rho==1){ + return((double)(Z<=R::qnorm(p, 0, 1, 1, 0))); + }else{ + return( R::pnorm( (R::qnorm(p, 0, 1, 1, 0) - sqrt(rho) * Z)/sqrt(1 - rho), 0, 1, 1, give_log)); + } +} + +static inline double dqnorm(double x){ + return 1/R::dnorm(R::qnorm(x, 0, 1, 1, 0), 0, 1, 0); +} + +double dshockprob(double p, double rho, double Z){ + return( R::dnorm((R::qnorm(p, 0, 1, 1, 0) - sqrt(rho) * Z)/sqrt(1-rho), 0, 1, 0) * dqnorm(p)/sqrt(1-rho) ); +} + +NumericMatrix shockprobvec2(double p, double rho, NumericVector Z){ + /* return a two column matrix with shockprob in the first column + and dshockprob in the second column*/ + NumericMatrix q(Z.size(), 2); + #pragma omp parallel for + for(size_t i = 0; i < Z.size(); i++){ + q[i] = shockprob(p, rho, Z[i]); + q[i + Z.size()] = dshockprob(p, rho, Z[i]); + } + return q; +} + +double shockseverity(double S, double Z, double rho, double p){ + if(p==0){ + return 0; + }else{ + return( exp(shockprob(S * p, rho, Z, true) - shockprob(p, rho, Z, true)) ); + } +} + +double quantile(const NumericVector& Z, const NumericVector& w, double p0){ + double cumw = 0; + size_t i; + for(i=0; i < Z.size(); i++) { + cumw += w[i]; + if(cumw > p0) { + break; + } + } + return( Z[i] ); +} + +// [[Rcpp::export]] +double fitprob(const NumericVector& Z, const NumericVector& w, double rho, double p0){ + if(p0 == 0){ + return 0.; + }else if(rho == 1){ + return R::pnorm(quantile(Z, w, p0), 0, 1, 1, 0); + }else{ + int one = 1; + double eps = 1e-12; + NumericMatrix q = shockprobvec2(p0, rho, Z); + int nZ = Z.size(); + double dp = (ddot_(&nZ, q(_,0).begin(), &one, w.begin(), &one) - p0) / + ddot_(&nZ, q(_,1).begin(), &one, w.begin(), &one); + double p = p0; + double phi = 1; + while(fabs(dp) > eps){ + while( (p - phi * dp) < 0 || (p - phi * dp) > 1){ + phi *= 0.8; + } + p -= phi * dp; + q = shockprobvec2(p, rho, Z); + dp = (ddot_(&nZ, q(_,0).begin(), &one, w.begin(), &one) - p0) / + ddot_(&nZ, q(_,1).begin(), &one, w.begin(), &one); + } + return p; + } +} + +// [[Rcpp::export]] +std::vector<double> stochasticrecov(double R, double Rtilde, + const NumericVector& Z, const NumericVector& w, + double rho, double porig, double pmod){ + if(porig==0){ + std::vector<double> q(Z.size(), R); + return q; + }else{ + double ptemp = (1 - R) / (1 - Rtilde) * porig; + double ptilde = fitprob(Z, w, rho, ptemp); + std::vector<double> q(Z.size()); + #pragma omp parallel for + for(size_t i = 0; i < Z.size(); i++){ + q[i] = fabs(1 - (1 - Rtilde) * exp( shockprob(ptilde, rho, Z[i], true) - + shockprob(pmod, rho, Z[i], true))); + } + return 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) { +// int i, j; +// double* dpshocked = malloc(sizeof(double) * (*ndp) * (*nZ)); +// double* ppshocked = malloc(sizeof(double) * (*ndp) * (*nZ)); +// int N2 = (*N) * (*N); +// double* qmat = malloc(sizeof(double) * N2 * (*nZ)); + +// double alpha = 1; +// double beta = 0; +// int one = 1; + +// #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[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); +// } + +// dgemv_("n", &N2, nZ, &alpha, qmat, &N2, wZ, &one, &beta, q, &one); + +// free(dpshocked); +// free(ppshocked); +// free(qmat); +// } + +NumericMatrix lossdistrib_joint_Z(const std::vector<double>& dp, const std::vector<double>& w, + const NumericMatrix S, int N, bool defaultflag, const std::vector<double>& rho, + const std::vector<double>& Z, const std::vector<double>& wZ) { + + int m = Z.size(), ndp = dp.size(); + NumericMatrix q(N, N); + int N2 = N * N; + double* qmat = new double[N2 * m]; + double alpha = 1; + double beta = 0; + int one = 1; + + #pragma omp parallel for + for(size_t i = 0; i < m; i++) { + std::vector<double> dpshocked(ndp); + matrix qZ(N, N, qmat+i); + for(size_t j =0; j < ndp; j++) { + dpshocked[j] = shockprob(dp[j], rho[j], Z[i], 0); + } + lossdistrib_joint(dpshocked, w, S(_,i).begin(), N, qZ, defaultflag); + } + dgemv_((char*)"n", &N2, &m, &alpha, qmat, &N2, wZ.data(), &one, &beta, q.begin(), &one); + delete[] qmat; + return q; +} + +// [[Rcpp::export]] +List BCloss_recov_dist(NumericMatrix defaultprob, + const std::vector<double>& issuerweights, NumericVector recov, + NumericVector Z, NumericVector w, + NumericVector rho, int N, bool defaultflag = false) { + /* + computes the loss and recovery distribution over time with a flat gaussian + correlation + inputs: + defaultprob: matrix of size dim1 x dim2. dim1 is the number of issuers + and dim2 number of time steps + issuerweights: vector of issuer weights (length dim1) + recov: vector of recoveries (length dim1) + Z: vector of factor values (length n) + w: vector of factor weights (length n) + rho: correlation beta vector (length dim1) + 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) + */ + int one = 1; + double alpha = 1; + double beta = 0; + + NumericMatrix L(N, defaultprob.ncol()); + NumericMatrix R(N, defaultprob.ncol()); + double* Lw = new double[N*Z.size()]; + double* Rw = new double[N*Z.size()]; + for(int t = 0; t < defaultprob.ncol(); t++) { + #pragma omp parallel + { + std::vector<double> gshocked(defaultprob.nrow()); + double* Rshocked = new double[defaultprob.nrow()]; + double* Sshocked = new double[defaultprob.nrow()]; + double g; + #pragma omp for + for(int i=0; i < Z.size(); i++) { + for(int j = 0; j < defaultprob.nrow(); j++) { + g = defaultprob(j, t);; + 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, issuerweights, Sshocked, N, N, + Lw+i*N, defaultflag); + lossdistrib(gshocked, issuerweights, Rshocked, N, N, + Rw+i*N, defaultflag); + } + delete[] Rshocked; + delete[] Sshocked; + } + int n = Z.size(); + dgemv_((char*)"n", &N, &n, &alpha, Lw, &N, w.begin(), &one, &beta, L(_,t).begin(), &one); + dgemv_((char*)"n", &N, &n, &alpha, Rw, &N, w.begin(), &one, &beta, R(_,t).begin(), &one); + } + delete[] Lw; + delete[] Rw; + return List::create(Named("L")=L, Named("R")=R); +} + +// void BCloss_recov_trunc(double *defaultprob, int *dim1, int *dim2, +// double *issuerweights, double *recov, double *Z, double *w, +// int *n, double *rho, int *N, double * K, int *defaultflag, +// double *ELt, double *ERt) { +// /* +// computes EL_t = E[(K-L_t)^+] and ER_t = E[(K-(1-R_t))^+] +// the the put options on loss and recovery over time +// with a flat gaussian correlation. +// inputs: +// defaultprob: matrix of size dim1 x dim2. dim1 is the number of issuers +// and dim2 number of time steps +// issuerweights: vector of issuer weights (length dim1) +// recov: vector of recoveries (length dim1) +// Z: vector of factor values (length n) +// w: vector of factor weights (length n) +// rho: correlation beta vector (length dim1) +// N: number of ticks in the grid +// K: the strike +// defaultflag: if true, computes the default distribution +// outputs: +// ELt: vector of length dim2 +// ERt: vector of length dim2 +// */ +// int t, i, j; +// double g, Ktilde; +// int one = 1; +// int T = (int) floor((*K) * (*N))+1; +// int Ttilde; +// double lu = 1./(*N+1); +// double* valL = malloc( T * sizeof(double)); +// posK(T, *K, lu, valL); +// double* EL = malloc( (*n) * sizeof(double)); +// double* ER = malloc( (*n) * sizeof(double)); +// double* Lw = malloc(T * (*n) * sizeof(double)); +// double alpha = 1; +// double beta = 0; + +// for(t=0; t < (*dim2); t++) { +// memset(Lw, 0, T * (*n) * sizeof(double)); +// #pragma omp parallel for private(j, g, Ktilde, Ttilde) +// for(i=0; i < *n; i++){ +// double* Rw = NULL; +// double* Rshocked = malloc((*dim1) * sizeof(double)); +// double* Sshocked = malloc((*dim1) * sizeof(double)); +// double* gshocked = malloc((*dim1) * sizeof(double)); +// double* gshockedbar = malloc((*dim1) * sizeof(double)); +// double* valR = NULL; + +// for(j=0; j < (*dim1); j++){ +// g = defaultprob[j + (*dim1) * t]; +// gshocked[j] = shockprob(g, rho[j], Z[i], 0); +// Sshocked[j] = shockseverity(1-recov[j], Z[i], rho[j], g); +// gshockedbar[j] = 1 - gshocked[j]; +// Rshocked[j] = 1 - Sshocked[j]; +// } + +// lossdistrib_truncated(gshocked, dim1, issuerweights, Sshocked, +// N, &T, defaultflag, Lw + i * T); +// ER[i] = 0; +// Ktilde = *K - ddot_(dim1, issuerweights, &one, Sshocked, &one); +// if(Ktilde > 0){ +// Ttilde = (int) floor(Ktilde * (*N))+1; +// Rw = calloc(Ttilde, sizeof(double)); +// lossdistrib_truncated(gshockedbar, dim1, issuerweights, Rshocked, +// N, &Ttilde, defaultflag, Rw); +// valR = malloc(Ttilde * sizeof(double)); +// posK(Ttilde, Ktilde, lu, valR); +// ER[i] = ddot_(&Ttilde, Rw, &one, valR, &one); +// } +// if(Rw != NULL){ +// free(Rw); +// } +// if(valR != NULL){ +// free(valR); +// } +// free(Rshocked); +// free(Sshocked); +// free(gshocked); +// free(gshockedbar); +// } +// dgemv_("t", &T, n, &alpha, Lw, &T, valL, &one, &beta, EL, &one); +// ELt[t] = ddot_(n, EL, &one, w, &one); +// ERt[t] = ddot_(n, ER, &one, w, &one); +// } +// free(Lw); +// free(EL); +// free(ER); +// free(valL); +// } + +// void BCloss_dist(double *defaultprob, int *dim1, int *dim2, +// double *issuerweights, double *recov, double *Z, double *w, +// int *n, double *rho, int *N, int *defaultflag, +// double *L) { +// /* +// computes the loss distribution over time with a flat gaussian +// correlation +// inputs: +// defaultprob: matrix of size dim1 x dim2. dim1 is the number of issuers +// and dim2 number of time steps +// issuerweights: vector of issuer weights (length dim1) +// recov: vector of recoveries (length dim1) +// Z: vector of factor values (length n) +// w: vector of factor weights (length n) +// rho: correlation beta vector (length dim1) +// N: number of ticks in the grid +// defaultflag: if true, computes the default distribution +// outputs: +// L: matrix of size (N, dim2) +// */ +// int t, i, j; +// double g; +// double *gshocked, *Sshocked, *Lw; +// int one = 1; +// double alpha = 1; +// double beta = 0; + +// gshocked = malloc((*dim1) * (*n) * sizeof(double)); +// Sshocked = malloc((*dim1) * (*n) * sizeof(double)); +// Lw = malloc((*N) * (*n) * sizeof(double)); + +// for(t=0; t < (*dim2); t++) { +// memset(Lw, 0, (*N) * (*n) * sizeof(double)); +// #pragma omp parallel for private(j, g) +// for(i=0; i < *n; i++){ +// for(j=0; j < (*dim1); j++){ +// g = defaultprob[j + (*dim1) * t]; +// gshocked[j+(*dim1)*i] = shockprob(g, rho[j], Z[i], 0); +// Sshocked[j+(*dim1)*i] = shockseverity(1-recov[j], Z[i], rho[j], g); +// } +// lossdistrib_blas(gshocked + (*dim1) * i, dim1, issuerweights, Sshocked + (*dim1)*i, N, defaultflag, +// Lw + i * (*N)); +// } +// dgemv_("n", N, n, &alpha, Lw, N, w, &one, &beta, L + t * (*N), &one); +// } +// free(gshocked); +// free(Sshocked); +// free(Lw); +// } +// diff --git a/src/lossdistrib.hpp b/src/lossdistrib.hpp new file mode 100644 index 0000000..0bde76a --- /dev/null +++ b/src/lossdistrib.hpp @@ -0,0 +1,75 @@ +#include <algorithm> + +extern "C" { + int dstev_(char* JOBZ, int* n, double* D, double* E, double* Z, int* ldz, double* WORK, int* INFO); + int dscal_(int* n, double* da, double* dx, int* incx); + int daxpy_(int* n, double* da, double* dx, int* incx, double* dy, int* incy); + void openblas_set_num_threads(int); + int dgemv_(char* trans, int *m, int *n, double* alpha, double* A, int* lda, + const double* x, int* incx, double* beta, double* y, int* incy); + double ddot_(int* n, const double* dx, int* incx, const double* dy, int* incy); +} + +double shockprob(double p, double rho, double Z, bool give_log); +double dshockprob(double p, double rho, double Z); + +//simple class which wraps a dynamically allocated array +//but allows it to be externally allocated +struct array { + array(int n) : n(n), extern_alloc(false) { + arr = new double[n]; + } + array(int n, double* const ptr) :n(n), extern_alloc(true) { + arr = new(ptr) double[n]; + } + ~array() { + if(!extern_alloc) { + delete[] arr; + } + } + double* begin() { + return arr; + } + double* begin() const { + return arr; + } + double* end() { + return arr+n; + } + double* end() const { + return arr+n; + } +private: + const int n; + double* arr; + const bool extern_alloc; +}; + +struct matrix { + matrix(int m, int n) : m(m), n(n), mat(m*n) {}; + matrix(int m, int n, double x) : m(m), n(n), mat(m*n) { + std::fill(mat.begin(), mat.end(), x); + } + matrix(int m, int n, double* ptr) : m(m), n(n), mat(m*n, ptr) {}; + double& operator()(int i, int j) { + return *(mat.begin()+j*m+i); + } + double operator()(int i, int j) const { + return *(mat.begin()+j*m+i); + } + double* operator()(int i) { + return mat.begin()+i*m; + } + double* data() { + return mat.begin(); + } + double& operator[](int i) { + return *(mat.begin()+i); + } + double operator[](int i) const { + return *(mat.begin()+i); + } +private: + int m, n; + array mat; +}; |