Skip to main content

rustyqlib/core/interpolation/
bilinear.rs

1//! Bilinear interpolation on a rectangular grid — the baseline 2-D
2//! scheme for strike x expiry volatility grids.
3
4use super::linear::bracket;
5use crate::core::errors::RustyQLibError;
6
7/// A rectangular grid `z[i][j]` sampled at `(xs[i], ys[j])`, both axes
8/// strictly increasing. Queries outside the grid are clamped (flat
9/// extrapolation), matching the wing behavior of the 1-D smiles.
10#[derive(Debug, Clone)]
11pub struct BilinearGrid {
12    xs: Vec<f64>,
13    ys: Vec<f64>,
14    z: Vec<Vec<f64>>,
15}
16
17impl BilinearGrid {
18    pub fn new(xs: &[f64], ys: &[f64], z: &[Vec<f64>]) -> Result<Self, RustyQLibError> {
19        if xs.len() < 2 || ys.len() < 2 {
20            return Err(RustyQLibError::invalid_input("bilinear", "need at least a 2 x 2 grid"));
21        }
22        if xs.windows(2).any(|w| w[1] <= w[0]) || ys.windows(2).any(|w| w[1] <= w[0]) {
23            return Err(RustyQLibError::invalid_input("bilinear", "grid axes must be strictly increasing"));
24        }
25        if z.len() != xs.len() || z.iter().any(|row| row.len() != ys.len()) {
26            return Err(RustyQLibError::invalid_input("bilinear", "z must be an xs.len() x ys.len() matrix"));
27        }
28        Ok(BilinearGrid { xs: xs.to_vec(), ys: ys.to_vec(), z: z.to_vec() })
29    }
30
31    pub fn eval(&self, x: f64, y: f64) -> f64 {
32        let (i, wx) = bracket(&self.xs, x);
33        let (j, wy) = bracket(&self.ys, y);
34        let z00 = self.z[i - 1][j - 1];
35        let z01 = self.z[i - 1][j];
36        let z10 = self.z[i][j - 1];
37        let z11 = self.z[i][j];
38        z00 * (1.0 - wx) * (1.0 - wy)
39            + z01 * (1.0 - wx) * wy
40            + z10 * wx * (1.0 - wy)
41            + z11 * wx * wy
42    }
43}
44
45#[cfg(test)]
46mod tests {
47    use super::*;
48
49    #[test]
50    fn reproduces_a_bilinear_function_exactly() {
51        // f = 2 + 3x + 4y + 5xy is in the bilinear span
52        let f = |x: f64, y: f64| 2.0 + 3.0 * x + 4.0 * y + 5.0 * x * y;
53        let xs = [0.0, 0.7, 1.5, 3.0];
54        let ys = [-1.0, 0.5, 2.0];
55        let z: Vec<Vec<f64>> =
56            xs.iter().map(|&x| ys.iter().map(|&y| f(x, y)).collect()).collect();
57        let grid = BilinearGrid::new(&xs, &ys, &z).unwrap();
58        for i in 0..=30 {
59            for j in 0..=30 {
60                let (x, y) = (i as f64 * 0.1, -1.0 + j as f64 * 0.1);
61                assert!((grid.eval(x, y) - f(x, y)).abs() < 1e-12, "({x}, {y})");
62            }
63        }
64    }
65
66    #[test]
67    fn clamps_outside_the_grid() {
68        let grid = BilinearGrid::new(
69            &[0.0, 1.0],
70            &[0.0, 1.0],
71            &[vec![1.0, 2.0], vec![3.0, 4.0]],
72        )
73        .unwrap();
74        assert_eq!(grid.eval(-5.0, -5.0), 1.0);
75        assert_eq!(grid.eval(9.0, 9.0), 4.0);
76        assert_eq!(grid.eval(0.5, -3.0), 2.0); // clamped in y only
77    }
78}