diff --git a/Cargo.toml b/Cargo.toml index 6cf5ed73..0e4dbcf2 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -97,9 +97,10 @@ plotly = "0.10.0" # https://docs.rs/plotly/latest/plotly/ rand = "0.8.5" # https://docs.rs/rand/latest/rand/ rand_distr = "0.4.3" # https://docs.rs/rand_distr/latest/rand_distr/ rayon = "1.9.0" # https://docs.rs/rayon/latest/rayon/ -rust_decimal = "1.34.3" # https://docs.rs/rust_decimal/latest/rust_decimal/ +rust_decimal = { version = "1.34.3", features = ["maths"] } # https://docs.rs/rust_decimal/latest/rust_decimal/ statrs = "0.17.1" # https://docs.rs/statrs/latest/statrs/ thiserror = "1.0.57" # https://docs.rs/thiserror/latest/thiserror/ +ordered-float = "5.1.0" # https://docs.rs/ordered-float/latest/ordered_float/ # https://docs.rs/ndarray/latest/ndarray/ ndarray = { version = "0.16.1", features = ["rayon"] } @@ -117,4 +118,17 @@ polars = { version = "0.44.0", features = ["docs-selection"] } serde = { version = "1.0.213", features = ["derive"] } # https://docs.rs/crate/pyo3/latest -pyo3 = {version = "0.26.0", features = ["time"] } +pyo3 = {version = "0.26.0", features = ["extension-module", "time"]} + +## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +## PYTHON BINDINGS +## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +# [lib] +# name = "RustQuant" +# crate-type = ["cdylib"] + +# [dependencies.pyo3] +# version = "0.22.0" +# features = ["extension-module"] +# features = ["abi3-py37", "extension-module"] diff --git a/crates/RustQuant_data/Cargo.toml b/crates/RustQuant_data/Cargo.toml index 7554bda3..0a5feef6 100644 --- a/crates/RustQuant_data/Cargo.toml +++ b/crates/RustQuant_data/Cargo.toml @@ -30,8 +30,8 @@ time = { workspace = true } plotly = { workspace = true } argmin = { workspace = true } argmin-math = { workspace = true } -pyo3 = { workspace = true } - +ordered-float = { workspace=true } +pyo3 = { workspace=true } ## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ## RUSTDOC CONFIGURATION diff --git a/crates/RustQuant_data/src/curves.rs b/crates/RustQuant_data/src/curves.rs index bea485b1..0c2b7e0d 100644 --- a/crates/RustQuant_data/src/curves.rs +++ b/crates/RustQuant_data/src/curves.rs @@ -35,7 +35,7 @@ use plotly::{common::Mode, Plot, Scatter}; use pyo3::{pyclass, pymethods, PyResult}; use std::collections::BTreeMap; use std::sync::Arc; -use time::Date; +use ordered_float::OrderedFloat; use RustQuant_math::interpolation::{ExponentialInterpolator, Interpolator, LinearInterpolator}; use RustQuant_stochastics::{CurveModel, NelsonSiegelSvensson}; @@ -84,7 +84,7 @@ pub enum InterpolationMethod { #[pyclass] pub struct Curve { /// The nodes of the curve. - nodes: BTreeMap, + nodes: BTreeMap, f64>, /// The type of the curve. curve_type: CurveType, @@ -93,7 +93,7 @@ pub struct Curve { interpolation_method: InterpolationMethod, /// Interpolator backend. - interpolator: Arc>, + interpolator: Arc>, /// Nelson-Siegel-Svensson curve parameters. nss: Option, @@ -104,12 +104,12 @@ impl Curve { /// Create a new curve. #[new] pub fn new( - dates: Vec, + dates: Vec, rates: Vec, curve_type: CurveType, interpolation_method: InterpolationMethod, ) -> PyResult { - let interpolator: Arc> = match interpolation_method { + let interpolator: Arc> = match interpolation_method { InterpolationMethod::Linear => { Arc::new(LinearInterpolator::new(dates.clone(), rates.clone())?) } @@ -122,10 +122,10 @@ impl Curve { InterpolationMethod::Lagrange => { todo!("Implement LagrangeInterpolator") } - }; + }; Ok(Self { - nodes: dates.into_iter().zip(rates.into_iter()).collect(), + nodes: dates.into_iter().zip(rates).map(|(a, b)| (OrderedFloat(a), b)).collect(), curve_type, interpolation_method, interpolator, @@ -133,38 +133,39 @@ impl Curve { }) } - /// Create a new Curve from a list of nodes. - #[staticmethod] - pub fn from_nodes( - nodes: BTreeMap, - curve_type: CurveType, - interpolation_method: InterpolationMethod, - ) -> PyResult { - let interpolator: Arc> = match interpolation_method { - InterpolationMethod::Linear => Arc::new(LinearInterpolator::new( - nodes.keys().cloned().collect(), - nodes.values().cloned().collect(), - )?), - InterpolationMethod::Exponential => Arc::new(ExponentialInterpolator::new( - nodes.keys().cloned().collect(), - nodes.values().cloned().collect(), - )?), - InterpolationMethod::CubicSpline => { - todo!("Implement CubicSplineInterpolator") - } - InterpolationMethod::Lagrange => { - todo!("Implement LagrangeInterpolator") - } - }; - Ok(Self { - nodes, - curve_type, - interpolation_method, - interpolator, - nss: None, - }) - } + // /// Create a new Curve from a list of nodes. + // #[staticmethod] + // pub fn from_nodes( + // nodes: BTreeMap, + // curve_type: CurveType, + // interpolation_method: InterpolationMethod, + // ) -> PyResult { + // let interpolator: Arc> = match interpolation_method { + // InterpolationMethod::Linear => Arc::new(LinearInterpolator::new( + // nodes.keys().cloned().collect(), + // nodes.values().cloned().collect(), + // )?), + // InterpolationMethod::Exponential => Arc::new(ExponentialInterpolator::new( + // nodes.keys().cloned().collect(), + // nodes.values().cloned().collect(), + // )?), + // InterpolationMethod::CubicSpline => { + // todo!("Implement CubicSplineInterpolator") + // } + // InterpolationMethod::Lagrange => { + // todo!("Implement LagrangeInterpolator") + // } + // }; + + // Ok(Self { + // nodes, + // curve_type, + // interpolation_method, + // interpolator, + // nss: None, + // }) + // } /// Get the interpolation method used by the curve. pub fn interpolation_method(&self) -> InterpolationMethod { @@ -172,32 +173,32 @@ impl Curve { } /// Get a rate from the curve. - pub fn get_rate(&self, date: Date) -> Option { - match self.nodes.get(&date) { + pub fn get_rate(&self, date: f64) -> Option { + match self.nodes.get(&OrderedFloat(date)) { Some(rate) => Some(*rate), None => self.interpolator.interpolate(date).ok(), } } /// Get a rate, and simultaneously add it to the nodes. - pub fn get_rate_and_insert(&mut self, date: Date) -> Option { - match self.nodes.get(&date) { + pub fn get_rate_and_insert(&mut self, date: f64) -> Option { + match self.nodes.get(&OrderedFloat(date)) { Some(rate) => Some(*rate), None => { let rate = self.interpolator.interpolate(date).ok()?; - self.nodes.insert(date, rate); + self.nodes.insert(OrderedFloat(date), rate); Some(rate) } } } /// Get multiple rates from the curve. - pub fn get_rates(&self, dates: Vec) -> Vec> { + pub fn get_rates(&self, dates: Vec) -> Vec> { dates.iter().map(|date| self.get_rate(*date)).collect() } /// Get multiple rates from the curve, and simultaneously add them to the nodes. - pub fn get_rates_and_insert(&mut self, dates: Vec) -> Vec> { + pub fn get_rates_and_insert(&mut self, dates: Vec) -> Vec> { dates .iter() .map(|date| self.get_rate_and_insert(*date)) @@ -205,30 +206,36 @@ impl Curve { } /// Set a rate in the curve. - pub fn set_rate(&mut self, date: Date, rate: f64) { - self.nodes.insert(date, rate); + pub fn set_rate(&mut self, date: f64, rate: f64) { + self.nodes.insert(OrderedFloat(date), rate); } /// Set multiple rates in the curve. - pub fn set_rates(&mut self, rates: Vec<(Date, f64)>) { + pub fn set_rates(&mut self, rates: Vec<(f64, f64)>) { for (date, rate) in rates { self.set_rate(date, rate); } } /// Get the first date in the curve. - pub fn first_date(&self) -> Option<&Date> { - self.nodes.keys().next() + pub fn first_date(&self) -> Option<&f64> { + match self.nodes.keys().next() { + Some(date) => Some(&date.0), + None => None, + } } /// Get the last date in the curve. - pub fn last_date(&self) -> Option<&Date> { - self.nodes.keys().next_back() + pub fn last_date(&self) -> Option<&f64> { + match self.nodes.keys().next_back() { + Some(date) => Some(&date.0), + None => None, + } } /// Get the dates of the curve. - pub fn dates(&self) -> Vec { - self.nodes.keys().cloned().collect() + pub fn dates(&self) -> Vec { + self.nodes.keys().map(|k| k.0).collect()//.cloned().collect() } /// Get the rates of the curve. @@ -257,7 +264,7 @@ impl Curve { } /// Get the bracketing indices for a specific index. - pub fn get_brackets(&self, index: Date) -> (Date, Date) { + pub fn get_brackets(&self, index: f64) -> (f64, f64) { let first = self.first_date().unwrap(); let last = self.last_date().unwrap(); @@ -268,10 +275,10 @@ impl Curve { return (*last, *last); } - let left = self.nodes.range(..index).next_back().unwrap().0; - let right = self.nodes.range(index..).next().unwrap().0; + let left = self.nodes.range(..OrderedFloat(index)).next_back().unwrap().0; + let right = self.nodes.range(OrderedFloat(index)..).next().unwrap().0; - return (*left, *right); + (left.0, right.0) } /// Shift the curve by a constant value. @@ -326,7 +333,7 @@ impl CostFunction for Curve { let y_model = x .into_iter() - .map(|date| curve_function(&nss, *date)) + .map(|date| curve_function(&nss, **date)) .collect::>(); let data = std::iter::zip(y, y_model); diff --git a/crates/RustQuant_math/src/interpolation/b_splines.rs b/crates/RustQuant_math/src/interpolation/b_splines.rs index 73813319..63f8c218 100644 --- a/crates/RustQuant_math/src/interpolation/b_splines.rs +++ b/crates/RustQuant_math/src/interpolation/b_splines.rs @@ -181,32 +181,6 @@ mod tests_b_splines { ); } - #[test] - fn test_b_spline_dates() { - let now = time::OffsetDateTime::now_utc(); - let knots: Vec = vec![ - now, - now + time::Duration::days(1), - now + time::Duration::days(2), - now + time::Duration::days(3), - now + time::Duration::days(4), - now + time::Duration::days(5), - now + time::Duration::days(6), - ]; - let control_points = vec![-1.0, 2.0, 0.0, -1.0]; - - let mut interpolator = BSplineInterpolator::new(knots.clone(), control_points, 2).unwrap(); - let _ = interpolator.fit(); - - assert_approx_equal!( - 1.375, - interpolator - .interpolate(knots[2] + time::Duration::hours(12)) - .unwrap(), - RUSTQUANT_EPSILON - ); - } - #[test] fn test_b_spline_inconsistent_parameters() { let knots = vec![0.0, 1.0, 2.0, 3.0, 4.0]; diff --git a/crates/RustQuant_math/src/interpolation/cubic_spline.rs b/crates/RustQuant_math/src/interpolation/cubic_spline.rs new file mode 100644 index 00000000..da0aea9a --- /dev/null +++ b/crates/RustQuant_math/src/interpolation/cubic_spline.rs @@ -0,0 +1,445 @@ +// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +// RustQuant: A Rust library for quantitative finance tools. +// Copyright (C) 2023 https://github.com/avhz +// Dual licensed under Apache 2.0 and MIT. +// See: +// - LICENSE-APACHE.md +// - LICENSE-MIT.md +// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +//! Module containing functionality for interpolation. + + +use crate::interpolation::{InterpolationIndex, InterpolationValue, Interpolator, IntoValue}; +use RustQuant_error::RustQuantError; + +// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +// STRUCTS & ENUMS +// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +/// Cubic Spline Interpolator. +pub struct CubicSplineInterpolator +where + IndexType: InterpolationIndex, + ValueType: InterpolationValue, +{ + /// X-axis values for the interpolator. + pub xs: Vec, + + /// Y-axis values for the interpolator. + pub ys: Vec, + + /// Whether the interpolator has been fitted. + pub fitted: bool, + + /// The second derivative of the spline at each point. + pub second_derivatives: Vec, + + /// The time steps between each point. + pub time_steps: Vec, +} + +// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +// IMPLEMENTATIONS, FUNCTIONS, AND MACROS +// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +impl CubicSplineInterpolator +where + IndexType: InterpolationIndex, + ValueType: InterpolationValue, +{ + /// Create a new CubicSplineInterpolator. + /// + /// # Errors + /// - `RustQuantError::UnequalLength` if ```xs.length() != ys.length()```. + /// + /// # Panics + /// Panics if NaN is in the index. + pub fn new( + xs: Vec, + ys: Vec, + ) -> Result, RustQuantError> { + if xs.len() != ys.len() { + return Err(RustQuantError::UnequalLength); + } + + let mut tmp: Vec<_> = xs.into_iter().zip(ys).collect(); + + tmp.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap()); + + let (xs, ys): (Vec, Vec) = tmp.into_iter().unzip(); + + Ok(Self { + xs, + ys, + fitted: false, + second_derivatives: vec![], + time_steps: vec![], + }) + } + + fn cholesky_decomposition(&self) -> Vec> { + let n: usize = self.time_steps.len(); + let mut lower_tri_matrix: Vec> = vec![vec![]; n - 1]; + + let mut prev_diag: ValueType = (ValueType::from_f64(2.0).unwrap() * (self.time_steps[0] + self.time_steps[1]).into_value()).square_root(); + lower_tri_matrix[0].push(prev_diag); + let mut prev_sub_diag: ValueType; + + for (i, time) in self.time_steps[1..(self.time_steps.len() - 1)].iter().enumerate() { + prev_sub_diag = time.into_value() / prev_diag; + lower_tri_matrix[i + 1].push(prev_sub_diag); + + prev_diag = (ValueType::from_f64(2.0).unwrap() * (*time + self.time_steps[i + 1]).into_value() - (prev_sub_diag * prev_sub_diag)).square_root(); + lower_tri_matrix[i + 1].push(prev_diag) + } + + lower_tri_matrix + } + + fn backward_substitution(&self, matrix: &[Vec], rhs: &[ValueType]) -> Vec { + + let matrix_len: usize = matrix.len(); + let mut prev_value: ValueType = rhs[matrix_len - 1] / matrix[matrix_len - 1][1]; + let mut result: Vec = vec![prev_value]; + + for i in (1..(matrix_len - 1)).rev() { + prev_value = (rhs[i] - matrix[i + 1][0] * prev_value) / matrix[i][1]; + result.insert(0, prev_value) + } + + result.insert(0, (rhs[0] - matrix[1][0] * prev_value) / matrix[0][0]); + + result + } + + fn forward_substitution(&self, matrix: &[Vec], rhs: &[ValueType]) -> Vec { + + let mut prev_value: ValueType = rhs[0] / matrix[0][0]; + let mut result: Vec = vec![prev_value]; + + for i in 1..matrix.len() { + prev_value = (rhs[i] - matrix[i][0] * prev_value) / matrix[i][1]; + result.push(prev_value) + } + + result + } + + // Compute the spline value at a given point. + #[allow(clippy::too_many_arguments)] + fn spline( + &self, + point: IndexType, + time_step: IndexType::Delta, + x_coordinate_lower: IndexType, + x_coordinate_upper: IndexType, + y_coordinate_lower: ValueType, + y_coordinate_upper: ValueType, + second_derivative_lower: ValueType, + second_derivative_upper: ValueType + ) -> ValueType { + let term_1: ValueType = ( + second_derivative_lower / (ValueType::from_f64(6.0).unwrap() + * time_step.into_value())) + * (x_coordinate_upper - point).into_value() + * (x_coordinate_upper - point).into_value() + * (x_coordinate_upper - point).into_value(); + + let term_2: ValueType = ( + second_derivative_upper / (ValueType::from_f64(6.0).unwrap() + * time_step.into_value())) + * (point - x_coordinate_lower).into_value() + * (point - x_coordinate_lower).into_value() + * (point - x_coordinate_lower).into_value(); + + let term_3: ValueType = ( + (y_coordinate_upper / time_step.into_value()) - (time_step.into_value() + * second_derivative_upper / ValueType::from_f64(6.0).unwrap()) + ) * (point - x_coordinate_lower).into_value(); + + let term_4: ValueType = ( + (y_coordinate_lower / time_step.into_value()) - (time_step.into_value() + * second_derivative_lower / ValueType::from_f64(6.0).unwrap()) + ) * (x_coordinate_upper - point).into_value(); + + term_1 + term_2 + term_3 + term_4 + } +} + +impl Interpolator + for CubicSplineInterpolator +where + IndexType: InterpolationIndex, + ValueType: InterpolationValue, +{ + // The first step in fitting the cubic spline interpolator + // is to compute the 2nd derivative at each point of the spline. + // This is done by solving a system of linear equations Ax = b + // where A is a tridiagonal matrix, b is a vector of values and + // x is the vector of second derivatives at each point. + // + // The system is solved by decomposing A into L * D * L^T + // and subsequently computing the inverse of A. + // + // Once all the second derivatives are computed, we can then + // compute the spline value at a given point. + fn fit(&mut self) -> Result<(), RustQuantError> { + + self.time_steps = self.xs.windows(2) + .map(|x| x[1] - x[0]) + .collect(); + + let mut rhs: Vec = self.xs.iter().zip(self.ys.iter()) + .collect::>() + .windows(3) + .map(|window| { + let x0: &IndexType = window[0].0; + let x1: &IndexType = window[1].0; + let x2: &IndexType = window[2].0; + + let y0: &ValueType = window[0].1; + let y1: &ValueType = window[1].1; + let y2: &ValueType = window[2].1; + + ValueType::from_f64(6.0).unwrap() * ( + (ValueType::one() / (*x2 - *x1).into_value()) * (*y2) + - ( + (ValueType::one() / (*x2 - *x1).into_value()) + + (ValueType::one() / (*x1 - *x0).into_value()) + ) * (*y1) + + (ValueType::one() / (*x1 - *x0).into_value()) * (*y0) + ) + }).collect(); + + let lower_tri_matrix: Vec> = self.cholesky_decomposition(); + rhs = self.forward_substitution(&lower_tri_matrix, &rhs); + self.second_derivatives = self.backward_substitution(&lower_tri_matrix, &rhs); + self.second_derivatives.insert(0, ValueType::zero()); + self.second_derivatives.push(ValueType::zero()); + + self.fitted = true; + + Ok(()) + } + + fn range(&self) -> (IndexType, IndexType) { + (*self.xs.first().unwrap(), *self.xs.last().unwrap()) + } + + fn add_point(&mut self, point: (IndexType, ValueType)) { + let idx = self.xs.partition_point(|&x| x < point.0); + self.xs.insert(idx, point.0); + self.ys.insert(idx, point.1); + self.fitted = false; + } + + fn interpolate(&self, point: IndexType) -> Result { + if !self.fitted { + return Err(RustQuantError::Unfitted); + } + let range = self.range(); + if point.partial_cmp(&range.0).unwrap() == std::cmp::Ordering::Less + || point.partial_cmp(&range.1).unwrap() == std::cmp::Ordering::Greater + { + return Err(RustQuantError::OutsideOfRange); + } + if let Ok(idx) = self + .xs + .binary_search_by(|p| p.partial_cmp(&point).expect("Cannot compare values.")) + { + return Ok(self.ys[idx]); + } + + let mut result: Option = None; + for k in 0..(self.xs.len() - 1) { + if self.xs[k] <= point && point < self.xs[k + 1] { + result = Some( + self.spline( + point, + self.time_steps[k], + self.xs[k], + self.xs[k + 1], + self.ys[k], + self.ys[k + 1], + self.second_derivatives[k], + self.second_derivatives[k + 1], + ) + ); + break; + } + } + + match result { + Some(value) => Ok(value), + None => Err(RustQuantError::OutsideOfRange), + } + } +} + +#[cfg(test)] +mod tests_cubic_spline_interpolation { + use super::*; + use RustQuant_utils::{assert_approx_equal, RUSTQUANT_EPSILON}; + + #[test] + fn test_natural_cubic_interpolation() { + + let xs: Vec = vec![0., 1., 2., 3., 4.]; + let ys: Vec = vec![0., 1., 16., 81., 256.]; + + let mut interpolator: CubicSplineInterpolator = CubicSplineInterpolator::new(xs, ys).unwrap(); + let _ = interpolator.fit(); + + assert_approx_equal!( + 36.660714285714285, + interpolator.interpolate(2.5).unwrap(), + RUSTQUANT_EPSILON + ); + } + + #[test] + fn test_cubic_interpolation_out_of_range() { + let xs: Vec = vec![1., 2., 3., 4., 5.]; + let ys: Vec = vec![1., 2., 3., 4., 5.]; + + let mut interpolator: CubicSplineInterpolator = CubicSplineInterpolator::new(xs, ys).unwrap(); + let _ = interpolator.fit(); + + assert!(interpolator.interpolate(6.).is_err()); + } + + #[test] + fn test_cubic_interpolation_add_range_and_refit() { + let xs: Vec = vec![0., 1., 2., 3., 4.]; + let ys: Vec = vec![0., 1., 16., 81., 256.]; + + let mut interpolator: CubicSplineInterpolator = CubicSplineInterpolator::new(xs, ys).unwrap(); + let _ = interpolator.fit(); + + interpolator.add_point((5.0, 625.0)); + let _ = interpolator.fit(); + + assert_approx_equal!( + 39.97368421052632, + interpolator.interpolate(2.5).unwrap(), + RUSTQUANT_EPSILON + ); + } + + #[test] + fn test_cubic_interpolation_add_range_and_no_refit() { + let xs: Vec = vec![0., 1., 2., 3., 4.]; + let ys: Vec = vec![0., 1., 16., 81., 256.]; + + let mut interpolator: CubicSplineInterpolator = CubicSplineInterpolator::new(xs, ys).unwrap(); + let _ = interpolator.fit(); + + interpolator.add_point((5.0, 625.0)); + + assert!( + matches!( + interpolator.interpolate(2.5), + Err(RustQuantError::Unfitted), + ) + ); + } + + #[test] + fn test_cubic_spline_unfitted() { + + let xs: Vec = vec![0., 1., 2., 3., 4.]; + let ys: Vec = vec![0., 1., 16., 81., 256.]; + + let interpolator: CubicSplineInterpolator = CubicSplineInterpolator::new(xs, ys).unwrap(); + + assert!( + matches!( + interpolator.interpolate(2.5), + Err(RustQuantError::Unfitted), + ) + ); + } +} + +#[cfg(test)] +mod tests_cubic_spline_helper_functions { + use super::*; + + #[test] + fn test_cholesky_decomposition() { + let xs: Vec = vec![1., 2., 3., 4., 5.]; + let ys: Vec = vec![1., 2., 3., 4., 5.]; + + let mut interpolator: CubicSplineInterpolator = CubicSplineInterpolator::new(xs, ys).unwrap(); + interpolator.time_steps = interpolator.xs.windows(2) + .map(|x| (x[1] - x[0])) + .collect(); + let result: Vec> = interpolator.cholesky_decomposition(); + let expected: &[&[f64]] = &[ + &[2.0], + &[0.5, 1.9364916731037085], + &[0.5163977794943222, 1.9321835661585918], + ]; + + assert_eq!(result, expected); + } + + #[test] + fn test_forward_substitution() { + let xs: Vec = vec![1., 2., 3., 4., 5.]; + let ys: Vec = vec![1., 2., 3., 4., 5.]; + + let interpolator: CubicSplineInterpolator = CubicSplineInterpolator::new(xs, ys).unwrap(); + let result: Vec = interpolator.forward_substitution( + &[ + [2.0].to_vec(), + [1.0, 2.0].to_vec(), + [1.0, 2.0].to_vec(), + ], + &[1.0, 1.0, 1.0] + ); + let expected: &[f64] = &[0.5, 0.25, 0.375]; + + assert_eq!(result, expected); + } + + #[test] + fn test_backward_substitution() { + let xs: Vec = vec![1., 2., 3., 4., 5.]; + let ys: Vec = vec![1., 2., 3., 4., 5.]; + + let interpolator: CubicSplineInterpolator = CubicSplineInterpolator::new(xs, ys).unwrap(); + let result: Vec = interpolator.backward_substitution( + &[ + [2.0].to_vec(), + [1.0, 2.0].to_vec(), + [1.0, 2.0].to_vec(), + ], + &[1.0, 1.0, 1.0] + ); + let expected: &[f64] = &[0.375, 0.25, 0.5]; + + assert_eq!(result, expected); + } + + #[test] + fn test_spline() { + let xs: Vec = vec![0., 1., 2., 3., 4.]; + let ys: Vec = vec![0., 1., 16., 81., 256.]; + + let interpolator: CubicSplineInterpolator = CubicSplineInterpolator::new(xs, ys).unwrap(); + let spline_result: f64 = interpolator.spline( + 2.5, + 1.0, + 2.0, + 3.0, + 16.0, + 81.0, + 32.75733333333333, + 156.85714285714286 + ); + + assert_eq!(spline_result, 36.64909523809524); + } +} diff --git a/crates/RustQuant_math/src/interpolation/exponential_interpolator.rs b/crates/RustQuant_math/src/interpolation/exponential_interpolator.rs index bdc90b6d..58d075c9 100644 --- a/crates/RustQuant_math/src/interpolation/exponential_interpolator.rs +++ b/crates/RustQuant_math/src/interpolation/exponential_interpolator.rs @@ -173,7 +173,6 @@ where #[cfg(test)] mod tests_exponential_interpolation { use super::*; - use time::macros::date; use RustQuant_utils::{assert_approx_equal, RUSTQUANT_EPSILON}; #[test] @@ -188,24 +187,4 @@ mod tests_exponential_interpolation { RUSTQUANT_EPSILON ); } - - #[test] - fn test_exponential_interpolation_dates() { - let d_1m = date!(1990 - 06 - 16); - let d_2m = date!(1990 - 07 - 17); - - let r_1m = 0.9870; - let r_2m = 0.9753; - - let dates = vec![d_1m, d_2m]; - let rates = vec![r_1m, r_2m]; - - let interpolator = ExponentialInterpolator::new(dates, rates).unwrap(); - - assert_approx_equal!( - 0.9854824711068088, - interpolator.interpolate(date!(1990 - 06 - 20)).unwrap(), - RUSTQUANT_EPSILON - ); - } } diff --git a/crates/RustQuant_math/src/interpolation/linear_interpolator.rs b/crates/RustQuant_math/src/interpolation/linear_interpolator.rs index f98dec10..d0e9a01f 100644 --- a/crates/RustQuant_math/src/interpolation/linear_interpolator.rs +++ b/crates/RustQuant_math/src/interpolation/linear_interpolator.rs @@ -129,7 +129,6 @@ where #[cfg(test)] mod tests_linear_interpolation { use super::*; - use time::macros::date; use RustQuant_utils::{assert_approx_equal, RUSTQUANT_EPSILON}; #[test] @@ -162,51 +161,4 @@ mod tests_linear_interpolation { assert!(interpolator.interpolate(6.).is_err()); } - - #[test] - fn test_linear_interpolation_dates() { - let now = time::OffsetDateTime::now_utc(); - - let xs = vec![ - now, - now + time::Duration::days(1), - now + time::Duration::days(2), - now + time::Duration::days(3), - now + time::Duration::days(4), - ]; - - let ys = vec![1., 2., 3., 4., 5.]; - - let mut interpolator = LinearInterpolator::new(xs.clone(), ys).unwrap(); - let _ = interpolator.fit(); - - assert_approx_equal!( - 2.5, - interpolator - .interpolate(xs[1] + time::Duration::hours(12)) - .unwrap(), - RUSTQUANT_EPSILON - ); - } - - #[test] - fn test_linear_interpolation_dates_textbook() { - let d_1m = date!(1990 - 06 - 16); - let d_2m = date!(1990 - 07 - 17); - - let r_1m = 0.9870; - let r_2m = 0.9753; - - let dates = vec![d_1m, d_2m]; - let rates = vec![r_1m, r_2m]; - - let interpolator = LinearInterpolator::new(dates, rates).unwrap(); - - let d = date!(1990 - 06 - 20); - assert_approx_equal!( - interpolator.interpolate(d).unwrap(), - 0.9854903225806452, - RUSTQUANT_EPSILON - ); - } } diff --git a/crates/RustQuant_math/src/interpolation/mod.rs b/crates/RustQuant_math/src/interpolation/mod.rs index d3b69d9a..2fa2bc62 100644 --- a/crates/RustQuant_math/src/interpolation/mod.rs +++ b/crates/RustQuant_math/src/interpolation/mod.rs @@ -7,7 +7,8 @@ // - LICENSE-MIT.md // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -use std::ops::{AddAssign, Div, Mul, Sub}; +use std::ops::{Add, AddAssign, Div, Mul, MulAssign, Neg, Sub}; +use rust_decimal::prelude::*; use RustQuant_error::RustQuantError; pub mod linear_interpolator; @@ -19,22 +20,53 @@ pub use exponential_interpolator::*; pub mod b_splines; pub use b_splines::*; +pub mod cubic_spline; +pub use cubic_spline::*; + + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ /// Trait describing requirements to be interpolated. -pub trait InterpolationValue: - num::Num + AddAssign + std::fmt::Debug + Copy + Clone + Sized -{ +pub trait InterpolationValue: num::Num + + num::FromPrimitive + + Neg + + AddAssign + + MulAssign + + SquareRoot + + Copy + + Clone + + Sized + + std::fmt::Display + + std::fmt::Debug + + Send + + Sync {} + +/// Trait to convert a Delta type into a value of `ValueType`. +pub trait IntoValue { + /// Convert `self` into `ValueType`. + fn into_value(self) -> ValueType; +} + +/// Trait to compute the square root of a value. +pub trait SquareRoot { + /// Calculate square root of a value. + fn square_root(self) -> ValueType; } /// Trait describing requirements to be an index of interpolation. pub trait InterpolationIndex: - Sub + PartialOrd + Copy + Clone + Sized + std::fmt::Display + Sub + PartialOrd + Copy + Clone + Sized + std::fmt::Display + Send + Sync { /// Type of the difference of `Self` - `Self` type Delta: Div - + Mul; + + Mul + + Add + + Sub + + IntoValue + + Copy + + Send + + Sync; /// Type of `Delta` / `Delta` type DeltaDiv: InterpolationValue; @@ -66,10 +98,19 @@ where fn add_point(&mut self, point: (IndexType, ValueType)); } -impl InterpolationValue for T where - T: num::Num + AddAssign + std::fmt::Debug + Copy + Clone + Sized -{ -} +impl InterpolationValue for T where T: num::Num + + num::FromPrimitive + + Neg + + AddAssign + + MulAssign + + SquareRoot + + Copy + + Clone + + Sized + + std::fmt::Display + + std::fmt::Debug + + Send + + Sync {} macro_rules! impl_interpolation_index { ($a:ty, $b:ty, $c:ty) => { @@ -80,31 +121,73 @@ macro_rules! impl_interpolation_index { }; } +macro_rules! impl_num_delta_into_value { + ($b:ty, $c:ty) => { + impl IntoValue<$c> for $b { + fn into_value(self) -> $c { + self + } + } + }; +} + +macro_rules! impl_int_square_root { + ($c:ty) => { + impl SquareRoot<$c> for $c { + fn square_root(self) -> $c { + self.isqrt() + } + } + }; +} + +macro_rules! impl_float_square_root { + ($c:ty) => { + impl SquareRoot<$c> for $c { + fn square_root(self) -> $c { + self.sqrt() + } + } + }; +} + +macro_rules! impl_decimal_square_root { + ($c:ty) => { + impl SquareRoot<$c> for $c { + fn square_root(self) -> $c { + self.sqrt().unwrap() + } + } + }; +} + // Implement InterpolationIndex for all signed integer types. impl_interpolation_index!(i8, i8, i8); +impl_num_delta_into_value!(i8, i8); +impl_int_square_root!(i8); impl_interpolation_index!(i16, i16, i16); +impl_num_delta_into_value!(i16, i16); +impl_int_square_root!(i16); impl_interpolation_index!(i32, i32, i32); +impl_num_delta_into_value!(i32, i32); +impl_int_square_root!(i32); impl_interpolation_index!(i64, i64, i64); +impl_num_delta_into_value!(i64, i64); +impl_int_square_root!(i64); impl_interpolation_index!(i128, i128, i128); +impl_num_delta_into_value!(i128, i128); +impl_int_square_root!(i128); impl_interpolation_index!(isize, isize, isize); - -// Implement InterpolationIndex for all unsigned integer types. -impl_interpolation_index!(u8, u8, u8); -impl_interpolation_index!(u16, u16, u16); -impl_interpolation_index!(u32, u32, u32); -impl_interpolation_index!(u64, u64, u64); -impl_interpolation_index!(u128, u128, u128); -impl_interpolation_index!(usize, usize, usize); +impl_num_delta_into_value!(isize, isize); +impl_int_square_root!(isize); // Implement InterpolationIndex for all floating point types. impl_interpolation_index!(f32, f32, f32); +impl_num_delta_into_value!(f32, f32); impl_interpolation_index!(f64, f64, f64); - -// Implement InterpolationIndex for date/time types. -impl_interpolation_index!(time::Date, time::Duration, f64); -impl_interpolation_index!(time::Time, time::Duration, f64); -impl_interpolation_index!(time::OffsetDateTime, time::Duration, f64); -impl_interpolation_index!(time::PrimitiveDateTime, time::Duration, f64); +impl_num_delta_into_value!(f64, f64); +impl_float_square_root!(f32); +impl_float_square_root!(f64); // Implement InterpolationIndex for Decimal type. impl_interpolation_index!( @@ -112,3 +195,5 @@ impl_interpolation_index!( rust_decimal::Decimal, rust_decimal::Decimal ); +impl_num_delta_into_value!(rust_decimal::Decimal, rust_decimal::Decimal); +impl_decimal_square_root!(rust_decimal::Decimal); diff --git a/crates/RustQuant_stochastics/src/curve_model.rs b/crates/RustQuant_stochastics/src/curve_model.rs index fd3f41f2..c34c20b9 100644 --- a/crates/RustQuant_stochastics/src/curve_model.rs +++ b/crates/RustQuant_stochastics/src/curve_model.rs @@ -12,11 +12,11 @@ /// A trait for curve models. pub trait CurveModel { /// Returns the forward rate for a given date. - fn forward_rate(&self, date: time::Date) -> f64; + fn forward_rate(&self, date: f64) -> f64; /// Returns the spot rate for a given date. - fn spot_rate(&self, date: time::Date) -> f64; + fn spot_rate(&self, date: f64) -> f64; /// Returns the discount factor for a given date. - fn discount_factor(&self, date: time::Date) -> f64; + fn discount_factor(&self, date: f64) -> f64; } diff --git a/crates/RustQuant_stochastics/src/nelson_siegel_svensson.rs b/crates/RustQuant_stochastics/src/nelson_siegel_svensson.rs index a3b71210..df5e0219 100644 --- a/crates/RustQuant_stochastics/src/nelson_siegel_svensson.rs +++ b/crates/RustQuant_stochastics/src/nelson_siegel_svensson.rs @@ -65,10 +65,10 @@ impl NelsonSiegelSvensson { impl CurveModel for NelsonSiegelSvensson { /// Returns the forward rate for a given date. - fn forward_rate(&self, date: Date) -> f64 { - assert!(date > today(), "Date must be in the future."); + fn forward_rate(&self, tau: f64) -> f64 { + // assert!(date > today(), "Date must be in the future."); - let tau = DayCountConvention::Actual_365_25.day_count_factor(today(), date); + // let tau = DayCountConvention::Actual_365_25.day_count_factor(today(), date); let term1 = f64::exp(-tau / self.lambda1); let term2 = (tau / self.lambda1) * term1; @@ -78,10 +78,10 @@ impl CurveModel for NelsonSiegelSvensson { } /// Returns the spot rate for a given date. - fn spot_rate(&self, date: Date) -> f64 { - assert!(date > today(), "Date must be in the future."); + fn spot_rate(&self, tau: f64) -> f64 { + // assert!(date > today(), "Date must be in the future."); - let tau = DayCountConvention::Actual_365_25.day_count_factor(today(), date); + // let tau = DayCountConvention::Actual_365_25.day_count_factor(today(), date); let term1 = self.lambda1 * (1. - f64::exp(-tau / self.lambda1)) / tau; let term2 = term1 - f64::exp(-tau / self.lambda1); @@ -91,10 +91,10 @@ impl CurveModel for NelsonSiegelSvensson { self.beta0 + self.beta1 * term1 + self.beta2 * term2 + self.beta3 * term3 } - fn discount_factor(&self, date: Date) -> f64 { - let tau = DayCountConvention::Actual_365_25.day_count_factor(today(), date); + fn discount_factor(&self, tau: f64) -> f64 { + // let tau = DayCountConvention::Actual_365_25.day_count_factor(today(), date); - f64::exp(-self.spot_rate(date) * tau / 100.) + f64::exp(-self.spot_rate(tau) * tau / 100.) } } @@ -105,7 +105,6 @@ impl CurveModel for NelsonSiegelSvensson { #[cfg(test)] mod tests_nelson_siegel_svensson { use super::*; - use time::Duration; #[test] fn test_nelson_siegel_svensson() { @@ -119,8 +118,8 @@ mod tests_nelson_siegel_svensson { }; let dates = (2..365 * 30) - .map(|i| today() + Duration::days(i)) - .collect::>(); + .map(|i| 1.0 / (i as f64)) + .collect::>(); let _forward_curve = dates .iter()