rspp 0.1.7

use super::bsm::Option;
extern crate docify;
use crate::tools::tool::jacbin_one_row;
use anyhow::{bail, Result};
use nalgebra::ComplexField;
/// https://www2.zhihu.com/question/327650676
/// https://github.com/RussoMarioDamiano/VanillaSABR/blob/master/SABR.py
/// https://github.com/RussoMarioDamiano/VanillaSABR/blob/master/Implementation.ipynb
use std::f64::consts::E;
use std::f64::INFINITY;
///α 是类似波动率（volatility like)的随机数，用他可以表示在atm时的波动率，α 反映波动率的大小水平，当α 越大，则F的波动率越大；当ß=1和 v=0，也即α 是常数时，SABR模型就是一个Black schloles模型框架
///beta是一个  的常数，他反映波动率曲线的倾斜率，ß 越大，波动率曲线越平；当ß=1时， 服从一个lognormal分布，当ß=0时，  服从一个normal分布，当标的是index，equity时，正常会选取ß=1，当标的是interest的话一般选取ß=0.5；
/// ρ 表示 与 的相关系数，其效果跟ß 一样，ρ 雨打，波动率曲线越平；
/// vv 是 volatility of volatility，他反映曲线的曲度，当 v 越大时，波动率曲线会更弯。
/// IV implied vol

/// 没有使用:ln(vol_atm)=ln(alpha)+(beta-1)*ln(base_foward)  最小二乘求偏导数后 beta 0,1 0.5
/// 没有使用:选定beta,回归得到alpha,使用atm vol逼近得到rho, vv  let alpha = v_atm * self.S.powf(beta - 1.);
impl Option {
    pub fn sabr_iv(&self, alpha: f64, beta: f64, rho: f64, vv: f64) -> f64 {
        let iv = hagan(self.K, self.S, self.T, alpha, beta, rho, vv);
        iv
    }
}

///cross check hagan func with https://github.com/RussoMarioDamiano/VanillaSABR/blob/master/SABR.py
pub fn hagan(K: f64, S: f64, T: f64, alpha: f64, beta: f64, rho: f64, vv: f64) -> f64 {
    let eps = 1e-07;
    let lnfk = (S / K).log(E);
    let b1 = 1. - beta;
    let fkb = (S * K).powf(1. - beta);

    let n1 = b1.powf(2.) * alpha.powf(2.) / (24. * fkb);
    let n2 = 0.25 * (rho * beta * alpha * vv) / (fkb.sqrt());
    let n3 = (2. - 3. * rho.powf(2.)) * vv.powf(2.) / 24.;
    let N = 1. + (n1 + n2 + n3) * T;
    if K == S {
        return alpha * N / fkb.sqrt();
    }
    let d1 = fkb.sqrt();
    let d2 = b1.powf(2.) * lnfk.powf(2.) / 24.;
    let d3 = b1.powf(4.) * lnfk.powf(4.) / 1920.;
    let D = d1 * (1. + d2 + d3);
    let z = (vv / alpha) * fkb.sqrt() * lnfk;
    if z.abs() > eps {
        let mut not_neg = 1. - 2. * rho * z + z.powf(2.);
        if not_neg < 0. {
            not_neg = 0.;
        }
        let _xz = (not_neg).sqrt() + z - rho;
        let x_z = (_xz / (1. - rho)).log(E);
        let IV = alpha * (N / D) * (z / x_z);
        return IV;
    }
    let IV = alpha * (N / D) * z;
    IV
}
pub struct SabrCali {
    pub Ks: Vec<f64>,
    pub ivs: Vec<f64>,
    pub S: f64,
    pub T: f64,
}
impl SabrCali {
    fn cost_func(&self, arg: &Vec<f64>) -> f64 {
        let mut cost = Vec::new();
        let n = self.Ks.len();
        for i in 0..n {
            let K = self.Ks[i];
            let iv = self.ivs[i];
            let pred = hagan(K, self.S, self.T, arg[0], 0.5, arg[1], arg[2]);
            let v = (iv - pred).powf(2.);
            cost.push(v);
        }
        cost.iter().sum()
    }
    pub fn jacbin_one_row(&self, x: &Vec<f64>) -> Vec<f64> {
        let dt = 1e-8;
        let mut jac = Vec::new();
        let y0 = self.cost_func(x);
        for col in 0..x.len() {
            // let djx = if x[col] != 0.0 { x[col] * dt.abs() } else { dt };
            let x1: Vec<f64> = x
                .iter()
                .enumerate()
                // 其他稳着不动,在col方向步进
                .map(|(i, v)| if i != col { *v } else { *v + dt })
                .collect();
            let res = (self.cost_func(&x1) - self.cost_func(x)) / dt;
            jac.push(res);
        }
        jac
    }
    /// alpha rho volvol
    fn fit(&self) -> Vec<f64> {
        let mut arg = vec![0.01, 0.00, 0.10];

        let precise = 1e-3;
        let ncol = arg.len();
        let bounds = vec![(0.0001, INFINITY), (-0.9999, 0.9999), (0.0001, INFINITY)];
        let mut cost = INFINITY;
        for i in 0..5000 {
            let jac = self.jacbin_one_row(&arg);
            let cost = self.cost_func(&arg);
            for col in 0..ncol {
                let limit = bounds[col];
                arg[col] = arg[col] - precise * cost / jac[col];
                if arg[col] < limit.0 {
                    arg[col] = limit.0;
                }
                if arg[col] > limit.1 {
                    arg[col] = limit.1;
                }
            }
            dbg!(cost, &jac);
            dbg!(&arg);
        }
        arg
    }
}
#[test]
#[docify::export]
fn test_sabr() {
    /// calibrate with pysabr
    ///from pysabr import hagan_2002_lognormal_sabr as sabr
    /// res=sabr.lognormal_vol(k, f, t, alpha, beta, rho, volvol)
    // let k = 5.;
    // let f = 5.5;
    // let t = 1.;
    // let alpha = 0.025;
    // let beta = 0.5;
    // let rho = -0.04;
    // let volvol = 0.29;
    // let res = hagan(k, f, t, alpha, beta, rho, volvol);
    // dbg!(res);
    let sabr = SabrCali {
        Ks: vec![1., 2., 3., 4., 5.],
        ivs: vec![0.4, 0.21, 0.18, 0.24, 0.35],
        S: 3.,
        T: 1.,
    };
    sabr.fit();
}