Attempt to fix failing XIRR edge case (near -1)

Cargo.toml

@@ -33,6 +33,7 @@ pyo3 = "0.20"
 numpy = "0.20"
 time = { version = "0.3", features = ["parsing", "macros"] }
 ndarray = "0.15"
+# num-complex = "0.4"
 
 [dev-dependencies]
 assert_approx_eq = "1.1"

python/pyxirr/_pyxirr.pyi

@@ -327,7 +327,7 @@ def rate(  # type: ignore[misc]
     fv: _Amount = 0,
     *,
     pmt_at_beginning: bool = False,
-    guess: _Guess = 0.1,
+    guess: _Guess = None,
 ) -> Optional[float]:
     ...
 
@@ -340,7 +340,7 @@ def rate(
     fv: _ScalarOrArrayLike[_Amount] = 0,
     *,
     pmt_at_beginning: _ScalarOrArrayLike[bool] = False,
-    guess: _Guess = 0.1,
+    guess: _Guess = None,
 ) -> List[Optional[float]]:
     ...
 
@@ -472,7 +472,7 @@ def cumprinc(
 def irr(
     amounts: _AmountArray,
     *,
-    guess: _Guess = 0.1,
+    guess: _Guess = None,
     silent: bool = False,
 ) -> Optional[float]:
     ...
@@ -493,7 +493,7 @@ def xirr(
     dates: _DateLikeArray,
     amounts: _AmountArray,
     *,
-    guess: _Guess = 0.1,
+    guess: _Guess = None,
     silent: bool = False,
     day_count: _DayCount = DayCount.ACT_365F,
 ) -> Optional[float]:
@@ -504,7 +504,7 @@ def xirr(
 def xirr(
     dates: _CashFlow,
     *,
-    guess: _Guess = 0.1,
+    guess: _Guess = None,
     silent: bool = False,
     day_count: _DayCount = DayCount.ACT_365F,
 ) -> Optional[float]:

src/core/mod.rs

@@ -4,6 +4,7 @@ mod models;
 mod optimize;
 pub mod periodic;
 mod scheduled;
+mod utils;
 
 pub use models::{DateLike, InvalidPaymentsError};
 pub use periodic::*;

src/core/optimize.rs

@@ -149,3 +149,84 @@ where
         .flat_map(|x| x.windows(2).map(|pair| brentq(f, pair[0], pair[1], 100)))
         .filter(|r| r.is_finite() && f(*r).abs() < 1e-3)
 }
+
+// use std::f64::consts::PI;
+//
+// use num_complex::Complex;
+
+// pub fn durand_kerner(coefficients: &[f64]) -> Vec<f64> {
+//     // https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/durand_kerner_roots.cpp#L109
+//
+//     // numerical errors less when the first coefficient is "1"
+//     // hence, we normalize the first coefficient
+//     let coefficients: Vec<_> = coefficients.iter().map(|x| x / coefficients[0]).collect();
+//     let degree = coefficients.len() - 1;
+//     let accuracy = 1e-10;
+//
+//     let mut roots: Vec<_> = (0..degree)
+//         .into_iter()
+//         .map(|i| Complex::<f64>::new(PI * (i as f64 / degree as f64), 0.0))
+//         .collect();
+//
+//     let mut prev_delta = f64::INFINITY;
+//
+//     for _ in 0..MAX_ITERATIONS {
+//         let mut tol_condition = 0.0f64;
+//
+//         for n in 0..degree {
+//             let numerator = polyval(&coefficients, roots[n]);
+//
+//             let mut denominator = Complex::new(1.0, 0.0);
+//             for i in 0..degree {
+//                 if i != n {
+//                     denominator *= roots[n] - roots[i];
+//                 }
+//             }
+//
+//             let delta = numerator / denominator;
+//
+//             if !delta.norm().is_finite() {
+//                 break;
+//             }
+//
+//             roots[n] -= delta;
+//
+//             tol_condition = tol_condition.max(delta.norm())
+//         }
+//
+//         if (prev_delta - tol_condition).abs() <= accuracy || tol_condition < accuracy {
+//             break;
+//         }
+//
+//         prev_delta = tol_condition
+//     }
+//
+//     roots.into_iter().map(|x| x.norm()).collect()
+// }
+
+// valuate a polynomial at specific values.
+// fn polyval(coefficients: &[f64], x: Complex<f64>) -> Complex<f64> {
+//     let degree = coefficients.len() - 1;
+//     coefficients.iter().enumerate().map(|(i, c)| c * x.powf((degree - i) as f64)).sum()
+// }
+
+// #[cfg(test)]
+// mod tests {
+//     use super::*;
+//     use assert_approx_eq::assert_approx_eq;
+//     use rstest::rstest;
+//
+//     #[rstest]
+//     fn test_durand_kerner() {
+//         let cf = &[-1e6, 5000., -3.];
+//         let roots = durand_kerner(cf);
+//
+//         dbg!(&roots);
+//         for root in roots {
+//             let guess = root - 1.;
+//             dbg!(guess);
+//             let rate = crate::core::irr(cf, Some(guess)).unwrap();
+//             assert_approx_eq!(crate::core::npv(rate, cf, None), 0.0);
+//         }
+//     }
+// }

src/core/periodic.rs

@@ -4,7 +4,8 @@ use ndarray::{ArrayD, ArrayViewD};
 
 use super::{
     models::{validate, InvalidPaymentsError},
-    optimize::{brentq_grid_search, newton_raphson, newton_raphson_with_default_deriv},
+    optimize::{brentq, brentq_grid_search, newton_raphson, newton_raphson_with_default_deriv},
+    utils,
 };
 use crate::{broadcast_together, broadcasting::BroadcastingError};
 
@@ -443,7 +444,21 @@ pub fn irr(values: &[f64], guess: Option<f64>) -> Result<f64, InvalidPaymentsErr
     let df = |rate| self::npv_deriv(rate, values);
     let is_good_rate = |rate: f64| rate.is_finite() && f(rate).abs() < 1e-3;
 
-    let rate = newton_raphson(guess.unwrap_or(0.1), &f, &df);
+    let guess = match guess {
+        Some(g) => g,
+        None => {
+            let (outflows, inflows) = utils::sum_negatives_positives(values);
+            inflows / -outflows - 1.0
+        }
+    };
+
+    let rate = newton_raphson(guess, &f, &df);
+
+    if is_good_rate(rate) {
+        return Ok(rate);
+    }
+
+    let rate = brentq(&f, -0.999999999999999, 100., 100);
 
     if is_good_rate(rate) {
         return Ok(rate);

src/core/private_equity.rs

@@ -1,12 +1,13 @@
 // https://www.insead.edu/sites/default/files/assets/dept/centres/gpei/docs/Measuring_PE_Fund-Performance-2019.pdf
 
+use super::utils;
 use super::InvalidPaymentsError;
 
 type Result<T> = std::result::Result<T, InvalidPaymentsError>;
 
 #[doc = include_str!("../../docs/_inline/pe/dpi.md")]
 pub fn dpi(amounts: &[f64]) -> Result<f64> {
-    let (cs, ds) = sum_negatives_positives(amounts);
+    let (cs, ds) = utils::sum_negatives_positives(amounts);
     check_zero_contributions(cs)?;
     Ok(ds / -cs)
 }
@@ -23,13 +24,13 @@ pub fn dpi_2(contributions: &[f64], distributions: &[f64]) -> Result<f64> {
 pub fn rvpi(contributions: &[f64], nav: f64) -> Result<f64> {
     let cs: f64 = contributions.iter().sum();
     check_zero_contributions(cs)?;
-    let sign = series_signum(contributions);
+    let sign = utils::series_signum(contributions);
     Ok(nav / (sign * cs))
 }
 
 #[doc = include_str!("../../docs/_inline/pe/tvpi.md")]
 pub fn tvpi(amounts: &[f64], nav: f64) -> Result<f64> {
-    let (cs, ds) = sum_negatives_positives(amounts);
+    let (cs, ds) = utils::sum_negatives_positives(amounts);
     check_zero_contributions(cs)?;
     Ok((ds + nav) / -cs)
 }
@@ -62,7 +63,7 @@ pub fn moic_2(contributions: &[f64], distributions: &[f64], nav: f64) -> Result<
 pub fn ks_pme_flows(amounts: &[f64], index: &[f64]) -> Result<Vec<f64>> {
     check_input_len(amounts, index)?;
 
-    Ok(pairwise_mul(amounts, &index_performance(index)))
+    Ok(utils::pairwise_mul(amounts, &index_performance(index)))
 }
 
 #[doc = include_str!("../../docs/_inline/pe/ks_pme_flows.md")]
@@ -75,8 +76,8 @@ pub fn ks_pme_flows_2(
     check_input_len(distributions, index)?;
 
     let px = index_performance(index);
-    let c = pairwise_mul(contributions, &px);
-    let d = pairwise_mul(distributions, &px);
+    let c = utils::pairwise_mul(contributions, &px);
+    let d = utils::pairwise_mul(distributions, &px);
 
     Ok((c, d))
 }
@@ -154,7 +155,7 @@ pub fn pme_plus_flows_2(
     nav: f64,
 ) -> Result<Vec<f64>> {
     let lambda = pme_plus_lambda_2(contributions, distributions, index, nav)?;
-    Ok(scale(distributions, lambda))
+    Ok(utils::scale(distributions, lambda))
 }
 
 #[doc = include_str!("../../docs/_inline/pe/pme_plus_lambda.md")]
@@ -176,8 +177,8 @@ pub fn pme_plus_lambda_2(
     check_input_len(distributions, index)?;
 
     let px = index_performance(index);
-    let ds = sum_pairwise_mul(distributions, &px);
-    let cs = sum_pairwise_mul(contributions, &px);
+    let ds = utils::sum_pairwise_mul(distributions, &px);
+    let cs = utils::sum_pairwise_mul(contributions, &px);
 
     Ok((cs - nav) / ds)
 }
@@ -212,7 +213,7 @@ pub fn pme_plus_2(
 #[doc = include_str!("../../docs/_inline/pe/ln_pme_nav.md")]
 pub fn ln_pme_nav(amounts: &[f64], index: &[f64]) -> Result<f64> {
     check_input_len(amounts, index)?;
-    Ok(-sum_pairwise_mul(amounts, &index_performance(index)))
+    Ok(-utils::sum_pairwise_mul(amounts, &index_performance(index)))
 }
 
 #[doc = include_str!("../../docs/_inline/pe/ln_pme_nav.md")]
@@ -302,33 +303,6 @@ fn index_performance(index: &[f64]) -> Vec<f64> {
     index.iter().map(|p| last / p).collect()
 }
 
-fn scale(values: &[f64], factor: f64) -> Vec<f64> {
-    values.iter().map(|v| v * factor).collect()
-}
-
-fn sum_pairwise_mul(a: &[f64], b: &[f64]) -> f64 {
-    a.iter().zip(b).map(|(x, y)| x * y).sum()
-}
-
-fn pairwise_mul(a: &[f64], b: &[f64]) -> Vec<f64> {
-    a.iter().zip(b).map(|(x, y)| x * y).collect()
-}
-
-fn series_signum(a: &[f64]) -> f64 {
-    // returns -1. if any item is negative, otherwise +1.
-    a.iter().any(|x| x.is_sign_negative()).then_some(-1.).unwrap_or(1.)
-}
-
-fn sum_negatives_positives(values: &[f64]) -> (f64, f64) {
-    values.iter().fold((0., 0.), |acc, x| {
-        if x.is_sign_negative() {
-            (acc.0 + x, acc.1)
-        } else {
-            (acc.0, acc.1 + x)
-        }
-    })
-}
-
 #[cfg(test)]
 mod tests {
     use super::*;

src/core/scheduled/xirr.rs

@@ -1,7 +1,7 @@
 use super::{year_fraction, DayCount};
 use crate::core::{
     models::{validate, DateLike, InvalidPaymentsError},
-    optimize::{brentq_grid_search, newton_raphson},
+    optimize::{brentq, newton_raphson},
 };
 
 pub fn xirr(
@@ -30,13 +30,24 @@ pub fn xirr(
         return Ok(rate);
     }
 
-    // strategy: closest to zero
-    // let breakpoints: &[f64] = &[0.0, 0.25, -0.25, 0.5, -0.5, 1.0, -0.9, -0.99999999999999, 1e9];
-    // strategy: pessimistic
-    let breakpoints: &[f64] = &[-0.99999999999999, -0.75, -0.5, -0.25, 0., 0.25, 0.5, 1.0, 1e6];
-    let rate = brentq_grid_search(&[breakpoints], &f).next();
+    let rate = brentq(&f, -0.999999999999999, 100., 100);
 
-    Ok(rate.unwrap_or(f64::NAN))
+    if is_good_rate(rate) {
+        return Ok(rate);
+    }
+
+    let mut step = 0.01;
+    let mut guess = -0.99999999999999;
+    while guess < 1.0 {
+        let rate = newton_raphson(guess, &f, &df);
+        if is_good_rate(rate) {
+            return Ok(rate);
+        }
+        guess += step;
+        step = (step * 1.1).min(0.1);
+    }
+
+    Ok(f64::NAN)
 }
 
 /// Calculate the net present value of a series of payments at irregular intervals.
@@ -53,13 +64,17 @@ pub fn xnpv(
 }
 
 pub fn sign_changes(v: &[f64]) -> i32 {
-    v.windows(2).map(|v| (v[0].signum() != v[1].signum()) as i32).sum()
+    v.windows(2)
+        .map(|p| (p[0].is_finite() && p[1].is_finite() && p[0].signum() != p[1].signum()) as i32)
+        .sum()
 }
 
 pub fn zero_crossing_points(v: &[f64]) -> Vec<usize> {
     v.windows(2)
         .enumerate()
-        .filter_map(|(i, p)| (p[0].signum() != p[1].signum()).then_some(i))
+        .filter_map(|(i, p)| {
+            (p[0].is_finite() && p[1].is_finite() && p[0].signum() != p[1].signum()).then_some(i)
+        })
         .collect()
 }
 
@@ -94,6 +109,7 @@ mod tests {
         assert_eq!(sign_changes(&[1., -2., 3.]), 2);
         assert_eq!(sign_changes(&[-1., 2., -3.]), 2);
         assert_eq!(sign_changes(&[-1., -2., -3.]), 0);
+        assert_eq!(sign_changes(&[1., f64::NAN, 3.]), 0);
     }
 
     #[rstest]
@@ -102,9 +118,10 @@ mod tests {
         assert_eq!(zero_crossing_points(&[1., -2., -3.]), vec![0]);
         assert_eq!(zero_crossing_points(&[1., -2., 3.]), vec![0, 1]);
         assert_eq!(zero_crossing_points(&[-1., -2., 3.]), vec![1]);
+        assert_eq!(zero_crossing_points(&[1., f64::NAN, 3.]), vec![]);
 
-        assert_eq!(zero_crossing_points(
-            &[7., 6., -3., -4., -7., 8., 3., -6., 7., 8.]),
+        assert_eq!(
+            zero_crossing_points(&[7., 6., -3., -4., -7., 8., 3., -6., 7., 8.]),
             vec![1, 4, 6, 7],
         );
     }