solver/models/

american_put.rs

use super::control::{ControlProblem, StateDerivatives};
use crate::numeric::finite_difference::discretization::{DimensionKind, Transport};
use crate::numeric::finite_difference::pde::PdeProblem;

/// American Put Option
///
/// Models the optimal exercise strategy for an American Put option.
/// This is formulated as an optimal stopping problem, but can be solved via HJB
/// by treating the exercise decision as a control.
///
/// Dynamics follow Geometric Brownian Motion for the underlying asset.
pub struct AmericanPut {
    /// Risk-free interest rate (r).
    pub risk_free_rate: f64,
    /// Volatility of the underlying asset ($\sigma$).
    pub volatility: f64,
    /// Strike price (K).
    pub strike: f64,
    /// Spatial discretization step size (dx) used for finite difference coefficients.
    pub dx: f64,
}

impl ControlProblem<1> for AmericanPut {
    type Control = ();

    fn optimize(&self, _t: f64, _state: &[f64; 1], _derivs: &StateDerivatives<1>) {}

    fn running_reward(&self, _t: f64, _state: &[f64; 1], _control: &()) -> f64 {
        // No running reward; the objective is terminal payoff plus the early
        // exercise obstacle enforced by `apply_constraint`.
        0.0
    }

    fn generator(
        &self,
        _t: f64,
        state: &[f64; 1],
        _control: &(),
        derivs: &StateDerivatives<1>,
    ) -> f64 {
        // Risk-neutral GBM generator L V = r S V_s + 0.5 sigma^2 S^2 V_ss.
        let s = state[0];
        self.risk_free_rate * s * derivs.grad[0]
            + 0.5 * self.volatility.powi(2) * s.powi(2) * derivs.hessian[0]
    }

    fn terminal(&self, state: &[f64; 1]) -> f64 {
        let s = state[0];
        (self.strike - s).max(0.0)
    }

    fn discount_rate(&self, _state: &[f64; 1]) -> f64 {
        self.risk_free_rate
    }

    fn constant_discount_rate(&self) -> Option<f64> {
        Some(self.risk_free_rate)
    }

    fn apply_constraint(&self, state: &[f64; 1], value: f64) -> f64 {
        let s = state[0];
        let intrinsic = (self.strike - s).max(0.0);
        value.max(intrinsic)
    }

    fn next_step(&self, _t: f64, state: &[f64; 1], dt: f64, noise: &[f64; 1]) -> [f64; 1] {
        // Log-Euler for risk-neutral GBM.
        let drift = (self.risk_free_rate - 0.5 * self.volatility.powi(2)) * dt;
        let diffusion = self.volatility * dt.sqrt() * noise[0];
        [state[0] * (drift + diffusion).exp()]
    }

    fn is_diffusion_dimension(&self, _dim: usize) -> bool {
        true
    }
}

impl PdeProblem<1> for AmericanPut {
    fn dimension_kind(&self, _dim: usize) -> DimensionKind {
        DimensionKind::Diffusion
    }

    fn transport(
        &self,
        _t: f64,
        state: &[f64; 1],
        _control: &(),
        _derivs: &StateDerivatives<1>,
    ) -> Transport<1> {
        // L V = r S V_s + 0.5 sigma^2 S^2 V_ss.
        // Upwind the drift and central-difference the diffusion:
        //   plus  = D/h^2 + max(b, 0)/h
        //   minus = D/h^2 + max(-b, 0)/h
        let s = state[0];
        let h = self.dx.max(1e-12);
        let b = self.risk_free_rate * s;
        let d = 0.5 * self.volatility.powi(2) * s.powi(2);
        Transport::new(
            [d / (h * h) + b.max(0.0) / h],
            [d / (h * h) + (-b).max(0.0) / h],
            0.0,
        )
    }
}

impl AmericanPut {
    /// Creates an American put model.
    ///
    /// # Examples
    ///
    /// ```
    /// use solver::models::american_put::AmericanPut;
    /// let put = AmericanPut::new(0.05, 0.2, 100.0, 1.0);
    /// assert_eq!(put.strike, 100.0);
    /// ```
    pub fn new(risk_free_rate: f64, volatility: f64, strike: f64, dx: f64) -> Self {
        Self {
            risk_free_rate,
            volatility,
            strike,
            dx,
        }
    }

    /// Updates the spatial step size used by the finite-difference transport
    /// coefficients.
    ///
    /// # Examples
    ///
    /// ```
    /// use solver::models::american_put::AmericanPut;
    /// let mut put = AmericanPut::new(0.05, 0.2, 100.0, 1.0);
    /// put.set_dx(0.5);
    /// assert_eq!(put.dx, 0.5);
    /// ```
    pub fn set_dx(&mut self, dx: f64) {
        self.dx = dx;
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::models::control::ControlProblem;

    fn default_model() -> AmericanPut {
        AmericanPut::new(0.05, 0.2, 100.0, 1.0)
    }

    #[test]
    fn terminal_is_payoff() {
        let m = default_model();
        assert_eq!(ControlProblem::terminal(&m, &[80.0]), 20.0);
        assert_eq!(ControlProblem::terminal(&m, &[100.0]), 0.0);
        assert_eq!(ControlProblem::terminal(&m, &[120.0]), 0.0);
    }

    #[test]
    fn constraint_enforces_early_exercise() {
        let m = default_model();
        assert_eq!(ControlProblem::apply_constraint(&m, &[80.0], 15.0), 20.0);
        assert_eq!(ControlProblem::apply_constraint(&m, &[80.0], 25.0), 25.0);
    }

    #[test]
    fn discount_rate_is_risk_free() {
        let m = default_model();
        assert_eq!(ControlProblem::discount_rate(&m, &[100.0]), 0.05);
        assert_eq!(ControlProblem::constant_discount_rate(&m), Some(0.05));
    }

    #[test]
    fn control_problem_constraint_matches_model() {
        let m = default_model();
        let v = ControlProblem::apply_constraint(&m, &[80.0], 15.0);
        assert_eq!(v, 20.0);
    }

    #[test]
    fn control_problem_generator_is_black_scholes_operator() {
        let m = default_model();
        let d = StateDerivatives::new([1.0], [0.0]);
        let g = ControlProblem::generator(&m, 0.0, &[100.0], &(), &d);
        // L V = r S V_s (diffusion term zero at V_ss = 0).
        assert!((g - m.risk_free_rate * 100.0).abs() < 1e-12);
    }

    #[test]
    fn control_problem_next_step_is_log_euler() {
        let m = default_model();
        let s = ControlProblem::next_step(&m, 0.0, &[100.0], 0.0, &[0.0]);
        assert!((s[0] - 100.0).abs() < 1e-12);
    }

    #[test]
    fn dimension_kind_is_diffusion() {
        let m = default_model();
        assert_eq!(PdeProblem::dimension_kind(&m, 0), DimensionKind::Diffusion);
    }

    #[test]
    fn transport_is_black_scholes_upwind_stencil() {
        let m = default_model();
        let s: f64 = 100.0;
        let d = 0.5 * m.volatility.powi(2) * s.powi(2);
        let h = m.dx;
        let t = PdeProblem::transport(&m, 0.0, &[s], &(), &StateDerivatives::new([0.0], [0.0]));
        assert!((t.plus[0] - (d / (h * h) + m.risk_free_rate * s / h)).abs() < 1e-12);
        assert!((t.minus[0] - d / (h * h)).abs() < 1e-12);
        assert_eq!(t.source, 0.0);
    }
}