solver/models/
american_put.rs

1use super::traits::{ControlOutput, Gradients, Model};
2
3/// American Put Option
4///
5/// Models the optimal exercise strategy for an American Put option.
6/// This is formulated as an optimal stopping problem, but can be solved via HJB
7/// by treating the exercise decision as a control.
8///
9/// Dynamics follow Geometric Brownian Motion for the underlying asset.
10pub struct AmericanPut {
11    /// Risk-free interest rate (r).
12    pub risk_free_rate: f64,
13    /// Volatility of the underlying asset ($\sigma$).
14    pub volatility: f64,
15    /// Strike price (K).
16    pub strike: f64,
17    /// Spatial discretization step size (dx) used for finite difference coefficients.
18    pub dx: f64,
19}
20
21impl Model<1> for AmericanPut {
22    type Process = ();
23
24    fn process(&self) {}
25
26    fn optimize(&self, state: &[f64; 1], _grads: &Gradients<1>) -> ControlOutput<1> {
27        let s = state[0];
28        let dx = self.dx;
29
30        let mu = self.risk_free_rate * s;
31        let sigma2 = self.volatility.powi(2) * s.powi(2);
32
33        // Central Difference Coefficients
34        // lambda_plus = sigma^2 / (2 dx^2) + mu / (2 dx)
35        // lambda_minus = sigma^2 / (2 dx^2) - mu / (2 dx)
36
37        let diff_term = sigma2 / (2.0 * dx * dx);
38        let drift_term = mu / (2.0 * dx);
39
40        let mut lambda_plus = diff_term + drift_term;
41        let mut lambda_minus = diff_term - drift_term;
42
43        // Upwinding if negative intensity
44        if lambda_minus < 0.0 {
45            // Use forward difference for drift (if mu > 0)
46            // mu * (V(x+dx) - V(x)) / dx
47            // lambda_plus = mu / dx
48            // But we also have diffusion.
49            // Standard upwind for convection-diffusion:
50            // lambda_plus = D/dx^2 + max(0, mu)/dx
51            // lambda_minus = D/dx^2 + max(0, -mu)/dx
52            // Here D = 0.5 sigma^2
53
54            lambda_plus = diff_term + mu.max(0.0) / dx;
55            lambda_minus = diff_term + (-mu).max(0.0) / dx;
56        }
57
58        ControlOutput {
59            lambda_plus: [lambda_plus],
60            lambda_minus: [lambda_minus],
61            flow: 0.0, // No cash flow, just value change
62        }
63    }
64
65    fn terminal(&self, state: &[f64; 1]) -> f64 {
66        let s = state[0];
67        (self.strike - s).max(0.0)
68    }
69
70    fn discount_rate(&self, _state: &[f64; 1]) -> f64 {
71        self.risk_free_rate
72    }
73
74    fn constant_discount_rate(&self) -> Option<f64> {
75        Some(self.risk_free_rate)
76    }
77
78    fn apply_constraint(&self, state: &[f64; 1], value: f64) -> f64 {
79        let s = state[0];
80        let intrinsic = (self.strike - s).max(0.0);
81        value.max(intrinsic)
82    }
83}
84
85impl AmericanPut {
86    pub fn new(risk_free_rate: f64, volatility: f64, strike: f64, dx: f64) -> Self {
87        Self {
88            risk_free_rate,
89            volatility,
90            strike,
91            dx,
92        }
93    }
94
95    // Helper to update dx if grid changes
96    pub fn set_dx(&mut self, dx: f64) {
97        self.dx = dx;
98    }
99}
100
101#[cfg(test)]
102mod tests {
103    use super::super::traits::{Gradients, Model};
104    use super::*;
105
106    fn zero_grads() -> Gradients<1> {
107        Gradients {
108            fwd: [0.0; 1],
109            bwd: [0.0; 1],
110        }
111    }
112
113    fn default_model() -> AmericanPut {
114        AmericanPut::new(0.05, 0.2, 100.0, 1.0)
115    }
116
117    #[test]
118    fn terminal_is_payoff() {
119        let m = default_model();
120        assert_eq!(m.terminal(&[80.0]), 20.0);
121        assert_eq!(m.terminal(&[100.0]), 0.0);
122        assert_eq!(m.terminal(&[120.0]), 0.0);
123    }
124
125    #[test]
126    fn constraint_enforces_early_exercise() {
127        let m = default_model();
128        assert_eq!(m.apply_constraint(&[80.0], 15.0), 20.0);
129        assert_eq!(m.apply_constraint(&[80.0], 25.0), 25.0);
130    }
131
132    #[test]
133    fn discount_rate_is_risk_free() {
134        let m = default_model();
135        assert_eq!(m.discount_rate(&[100.0]), 0.05);
136        assert_eq!(m.constant_discount_rate(), Some(0.05));
137    }
138
139    #[test]
140    fn transport_coefficients_positive() {
141        let m = default_model();
142        let grads = zero_grads();
143        let ctrl = m.optimize(&[100.0], &grads);
144        assert!(ctrl.lambda_plus[0] > 0.0);
145        assert!(ctrl.lambda_minus[0] > 0.0);
146    }
147}
solver/models/american_put.rs

solver/models/
american_put.rs
