main.rs

use num_complex::Complex;
use probability::prelude::*;
use rayon::prelude::*;
use serde_derive::Deserialize;
use serde_json::json;
use std::error::Error;
use std::fs::File;
use std::io::prelude::*; //needed for write
use std::io::BufRead;
use std::io::BufReader;

#[derive(Debug, Deserialize)]
#[serde(rename_all = "camelCase")]
struct Parameters {
    lambda: f64,
    q: f64,
    num_u: usize,
    x_min: f64,
    x_max: f64,
    num_x: Option<usize>,
    r_squared: Vec<f64>,
}

//for the bivariate cdf...see https://apps.dtic.mil/dtic/tr/fulltext/u2/a125033.pdf
fn mean_h(h: f64, rho_sq: f64, normal: &Gaussian) -> f64 {
    -normal.density(h) * rho_sq / normal.distribution(h)
}
fn var_h(h: f64, rho_sq: f64, mean: f64) -> f64 {
    1.0 + rho_sq * h * mean - mean.powi(2)
}

//needed to convert from variance/covariance in Merton to Credit Risk Plus
fn biv_gaussian_orig(x1: f64, x2: f64, rho_sq: f64, normal: &Gaussian) -> f64 {
    let mean = mean_h(x1, rho_sq, &normal);
    let variance = var_h(x1, rho_sq, mean);
    normal.distribution(x1) * normal.distribution((x2 - mean) / variance.sqrt())
}

fn biv_gaussian(x1: f64, x2: f64, rho_sq: f64, normal: &Gaussian) -> f64 {
    let c = if x1 > x2 { x1 } else { x2 };
    let d = if x1 > x2 { x2 } else { x1 };
    if c < 0.0 {
        biv_gaussian_orig(c, d, rho_sq, normal)
    } else {
        normal.distribution(d) - biv_gaussian_orig(-c, d, -rho_sq, normal)
    }
}

//https://www.nag.com/doc/techrep/pdf/tr1_16.pdf
fn cov_merton(p: f64, rho_sq: f64) -> f64 {
    let normal = Gaussian::new(0.0, 1.0);
    let x = normal.inverse(p);
    biv_gaussian(x, x, rho_sq, &normal)
}
//converts to variance
fn get_systemic_variance(p: f64, rho_sq: f64) -> f64 {
    cov_merton(p, rho_sq) / p.powi(2) - 1.0
}

fn gamma_mgf(variance: Vec<f64>) -> impl Fn(&[Complex<f64>]) -> Complex<f64> {
    move |u_weights: &[Complex<f64>]| -> Complex<f64> {
        u_weights
            .iter()
            .zip(&variance)
            .map(|(u, v)| -(1.0 - v * u).ln() / v)
            .sum::<Complex<f64>>()
            .exp()
    }
}

fn main() -> Result<(), Box<dyn Error>> {
    let args: Vec<String> = std::env::args().collect();
    let Parameters {
        lambda,
        q,
        num_u,
        x_min,
        x_max,
        r_squared,
        ..
    } = serde_json::from_str(args[1].as_str())?;
    let num_w = r_squared.len();
    let p = 0.05; //just for test. The p has to be constant per r-squared.  This could, for example, be the average PD per industry in the book.
    let systemic_variance = r_squared
        .iter()
        .map(|&r_sq| get_systemic_variance(p, r_sq))
        .collect::<Vec<_>>();

    let liquid_fn = loan_ec::get_liquidity_risk_fn(lambda, q);
    let lgd_fn = |u: &Complex<f64>, l: f64, lgd_v: f64| {
        if lgd_v > 0.0 {
            cf_functions::gamma_cf(&(-u * l), 1.0 / lgd_v, lgd_v)
        } else {
            (-u * l).exp()
        }
    };
    let log_lpm_cf = loan_ec::get_log_lpm_cf(&lgd_fn, &liquid_fn);

    let mut discrete_cf = loan_ec::EconomicCapitalAttributes::new(num_u, num_w);

    let u_domain: Vec<Complex<f64>> = fang_oost::get_u_domain(num_u, x_min, x_max).collect();

    let f = File::open(args[2].as_str())?;
    let file = BufReader::new(&f);
    for line in file.lines() {
        let loan: loan_ec::Loan = serde_json::from_str(&line?)?;
        discrete_cf.process_loan(&loan, &u_domain, &log_lpm_cf);
    }
    let v_mgf = gamma_mgf(systemic_variance);
    let final_cf: Vec<Complex<f64>> = discrete_cf.get_full_cf(&v_mgf);
    if args.len() > 3 {
        let x_domain: Vec<f64> = fang_oost::get_x_domain(1024, x_min, x_max).collect();
        let density: Vec<f64> = fang_oost::get_density(
            x_min,
            x_max,
            fang_oost::get_x_domain(1024, x_min, x_max),
            &final_cf,
        )
        .collect();
        let json_results = json!({"x":x_domain, "density":density});
        let mut file_w = File::create(args[3].as_str())?;
        file_w.write_all(json_results.to_string().as_bytes())?;
    }

    let max_iterations = 100;
    let tolerance = 0.0001;
    let cf_dist_utils::RiskMetric {
        expected_shortfall,
        value_at_risk,
    } = cf_dist_utils::get_expected_shortfall_and_value_at_risk_discrete_cf(
        0.01,
        x_min,
        x_max,
        max_iterations,
        tolerance,
        &final_cf,
    )?;

    println!("This is ES: {}", expected_shortfall);
    println!("This is VaR: {}", value_at_risk);
    Ok(())
}
#[cfg(test)]
mod tests {
    use super::*;
    use approx::*;
    #[test]
    fn test_gamma_cf() {
        let kappa = 2.0;
        let u = Complex::new(0.5, 0.5);
        let theta = 0.5;
        let cf = gamma_mgf(vec![theta]);
        let result = cf(&vec![u]);
        let expected = (1.0 - u * theta).powf(-kappa);
        assert_eq!(result, expected);
    }
    #[test]
    fn cov_merton_compare_r_1() {
        let p = 0.05;
        let rho_sq = 0.5;
        let result = cov_merton(p, rho_sq);
        assert_abs_diff_eq!(result, 0.01218943, epsilon = 0.001);
    }
    #[test]
    fn biv_gaussian_compare_r_1() {
        let rho_sq = 0.7;
        let k = 1.0;
        let h = 0.2;
        let normal = Gaussian::new(0.0, 1.0);
        let result = biv_gaussian(h, k, rho_sq, &normal);
        assert_abs_diff_eq!(result, 0.55818, epsilon = 0.00001);
    }
    #[test]
    fn biv_gaussian_compare_r_2() {
        let rho_sq = -0.7;
        let k = -1.0;
        let h = 0.2;
        let normal = Gaussian::new(0.0, 1.0);
        let result = biv_gaussian_orig(k, h, rho_sq, &normal);
        assert_abs_diff_eq!(result, 0.02108, epsilon = 0.00001);
    }
    #[test]
    fn mean_compare_r_1() {
        let rho_sq = -0.7;
        let normal = Gaussian::new(0.0, 1.0);
        let h = -1.0;
        let result = mean_h(h, rho_sq, &normal);
        assert_abs_diff_eq!(result, 1.0676, epsilon = 0.00001);
    }
    #[test]
    fn var_compare_r_1() {
        let rho_sq = -0.7;
        let h = -1.0;
        let normal = Gaussian::new(0.0, 1.0);
        let mean = mean_h(h, rho_sq, &normal);
        let var = var_h(h, rho_sq, mean);
        assert_abs_diff_eq!(var.sqrt(), 0.77946, epsilon = 0.00001);
    }
    #[test]
    fn cov_compare_r_2() {
        let p = 0.02;
        let rho_sq = (0.6 as f64).powi(2);
        let result = cov_merton(p, rho_sq) - p.powi(2);
        assert_abs_diff_eq!(result, 0.001689042, epsilon = 0.00001);
    }
    #[test]
    fn test_systemic() {
        let p = 0.02;
        let rho_sq = (0.6 as f64).powi(2);
        let cov_merton = cov_merton(p, rho_sq) - p.powi(2);
        let systemic_variance = get_systemic_variance(p, rho_sq);
        let cov_creditriskplus = systemic_variance * p.powi(2);
        assert_abs_diff_eq!(cov_merton, cov_creditriskplus, epsilon = 0.00001);
    }
}