Hi. can I get some help with this math problem?

[luxiant]

Hi. I'm now trying a project to analyze option data. During the project, I faced with a problem to calculate implied volatility.

Here is the part of my code.

    pub fn calculate_implied_volatility_adjusted(snd_params: &HashMap<String, f64>, option_type: &str, underlying: f64, strike: f64, day_count_convention: f64, risk_free: f64, option_price: f64) -> Option<f64> {
        match option_type {
            "C" => {
                let f = |x: AD| -> AD {
                    let component_one = (underlying/strike).ln();
                    let component_two = 0.5*f64::powf(x.x(), 2.0);
                    let component_three = (x.x())*(day_count_convention.sqrt());
                    let d_one = (component_one+((risk_free+component_two)*day_count_convention))/component_three;
                    let d_two = (component_one+((risk_free-component_two)*day_count_convention))/component_three;
                    // calculate stats of d1, and d2, according to skew normal distribution instead of normal distribution
                    let adjusted_nd_one = snd_cumulative_distribution(d_one, &snd_params);
                    let adjusted_nd_two = snd_cumulative_distribution(d_two, &snd_params);
                    let discount = (-risk_free*day_count_convention).exp();
                    AD0((underlying*adjusted_nd_one)-(strike*discount*adjusted_nd_two)-option_price)
                };

                let adjusted_vol_option = bisection(f, (0f64, 10f64), 100, 0.0001);
                match adjusted_vol_option {
                    Ok(x) => {
                        return Some(x);
                    },
                    Err(_) => {
                        return None;
                    }
                }
            },
            "P" => {
                let f = |x: AD| -> AD {
                    let component_one = (underlying/strike).ln();
                    let component_two = 0.5*f64::powf(x.x(), 2.0);
                    let component_three = (x.x())*(day_count_convention.sqrt());
                    let d_one = (component_one+((risk_free+component_two)*day_count_convention))/component_three;
                    let d_two = (component_one+((risk_free-component_two)*day_count_convention))/component_three;
                    // calculate stats of d1, and d2, according to skew normal distribution instead of normal distribution
                    let adjusted_nd_one = snd_cumulative_distribution(-d_one, &snd_params);
                    let adjusted_nd_two = snd_cumulative_distribution(-d_two, &snd_params);
                    let discount = (-risk_free*day_count_convention).exp();
                    AD0((strike*discount*adjusted_nd_two)-(underlying*adjusted_nd_one)-option_price)
                };

                let adjusted_vol_option = bisection(f, (0f64, 10f64), 100, 0.0001);
                match adjusted_vol_option {
                    Ok(x) => {
                        return Some(x);
                    },
                    Err(_) => {
                        return None;
                    }
                }
            },
            &_ => todo!()
        }
    }

But this does not derive any solution. When I check how AD function f is doing, it just calculates 0.0, 5.0, 10.0, 0.0, 5.0, 10.0..... repeatedly until the end of the iteration.

What am I missing?

[Axect]

Hi @luxiant,

The internal implementation of the Bisection method in Peroxide is outlined below:

I(a, b) => {
    let f_ad = ADFn::new(|x| self.f(x));
    let x = 0.5 * (a + b);
    let fa = f_ad.call_stable(a);
    let fx = f_ad.call_stable(x);
    let fb = f_ad.call_stable(b);
    if (a - b).abs() <= self.tol {
        self.find = RootBool::Find;
        self.root = x;
    } else {
        if fa * fx < 0f64 {
            self.curr = I(a, x);
        } else if fx * fb < 0f64 {
            self.curr = I(x, b);
        } else if fx == 0f64 {
            self.find = RootBool::Find;
            self.root = x;
        } else {
            self.find = RootBool::Error;
        }
    }
}

This algorithm is based on the assumption that the root of f should be within an open interval (a, b). This means there should exist a c such that f(c) = 0 and a < c < b.

If this assumption is not met, the Bisection method may not converge. Consider the following example:

use peroxide::fuga::*;

fn main() {
    let h = |x: AD| -> AD {
        let y = if x.x() < 5.0 { -5.0 } else { 5.0 };
        println!("y = {:>2}", y);
        AD0(y)
    };
    let sol = bisection(h, (0.0, 10.0), 5, 0.0001);
    match sol {
        Ok(x) => println!("{}", x),
        Err(_) => println!("No solution found"),
    }
}

The output of this function is as follows:

y = -5 # <- h(0)
y =  5 # <- h(5)
y =  5 # <- h(10)
y = -5
y = -5
y =  5
y = -5
y = -5
y =  5
y = -5
y = -5
y =  5
y = -5
y = -5
y =  5
No solution found

The first three values correspond to h(0), h(5), and h(10) (i.e., a, (a + b) / 2, b). Since h(0) * h(5) < 0, the next interval becomes (0, 5). The following three values correspond to h(0), h(2.5), and h(5). As h(2.5) * h(5) < 0, the next interval is (2.5, 3.75, 5). However, the iteration continues without finding a solution because there isn't one within the interval.

It seems a similar issue may be occurring in your case. Therefore, it might be necessary to verify that the solution lies within an open interval.

[luxiant]

Thanks. Now I found that this part is not a root cause of the error.