solver/models/
merton_jump_lognormal.rs

1//! Merton's log-utility portfolio with log-normal jumps in the risky asset.
2//!
3//! The jump-size multiplier `Y` is log-normal, `ln Y ~ N(m, delta^2)`, so the
4//! jump distribution is a continuum rather than a single deterministic size.
5//! As in [`MertonJump`] the jump genuinely
6//! changes the optimal control, but here the policy has no closed form: it is
7//! the root of a transcendental first-order condition.
8//!
9//! # Dynamics
10//!
11//! ```text
12//! dS/S  = mu dt + sigma dW + (Y - 1) dN,   ln Y ~ N(m, delta^2)
13//! dx    = x [ r + u (mu - r) ] dt + x u sigma dW + x u (Y - 1) dN
14//! ```
15//!
16//! where `N` is a Poisson process with intensity `lambda`. Log utility requires
17//! wealth to stay positive, so the admissible control is `0 <= u <= 1`.
18//!
19//! # Exact solution
20//!
21//! The value function is `V(t, x) = ln x + B (T - t)` with
22//!
23//! ```text
24//! B = r + u* (mu - r) - 0.5 (u* sigma)^2 + lambda E[ ln(1 + u* (Y - 1)) ],
25//! ```
26//!
27//! and `u*` solves the transcendental first-order condition
28//!
29//! ```text
30//! 0 = (mu - r) - sigma^2 u* + lambda E[ (Y - 1) / (1 + u* (Y - 1)) ].
31//! ```
32//!
33//! The expectations are computed with 32-point Gauss-Hermite quadrature on the
34//! standard normal `z = (ln Y - m) / delta`, and `u*` is bracketed on
35//! `[0, 1]` and refined by bisection. With `delta = 0` the jump degenerates to
36//! the deterministic multiplier `y = e^m`, reducing this model to
37//! [`MertonJump`].
38use crate::models::merton_jump::MertonJump;
39
40/// 32-point Gauss-Hermite nodes for the standard normal `z ~ N(0, 1)`.
41const GH_NODES: [f64; 32] = [
42    -1.007_742_267_422_946_5e1,
43    -9.064_399_210_702_405,
44    -8.219_728_765_382_245,
45    -7.460_755_754_121_519,
46    -6.755_930_830_540_705,
47    -6.088_964_309_076_987,
48    -5.450_033_273_623_428,
49    -4.832_604_613_244_489,
50    -4.232_021_109_995_41,
51    -3.644_781_249_880_833,
52    -3.068_135_169_013_121_6,
53    -2.499_840_415_187_395_4,
54    -1.938_004_905_925_717_4,
55    -1.380_980_199_272_144_2,
56    -8.272_849_037_797_652e-1,
57    -2.755_464_192_302_758_4e-1,
58    2.755_464_192_302_758_4e-1,
59    8.272_849_037_797_652e-1,
60    1.380_980_199_272_144_2,
61    1.938_004_905_925_717_4,
62    2.499_840_415_187_395_4,
63    3.068_135_169_013_121_6,
64    3.644_781_249_880_833,
65    4.232_021_109_995_41,
66    4.832_604_613_244_489,
67    5.450_033_273_623_428,
68    6.088_964_309_076_987,
69    6.755_930_830_540_705,
70    7.460_755_754_121_519,
71    8.219_728_765_382_245,
72    9.064_399_210_702_405,
73    1.007_742_267_422_946_5e1,
74];
75
76/// Gauss-Hermite weights, normalized so `sum_i w_i f(z_i) ~ E[f(Z)]`.
77const GH_WEIGHTS: [f64; 32] = [
78    4.124_607_489_018_283_5e-23,
79    5.208_449_591_960_913e-19,
80    6.755_290_223_670_204e-16,
81    2.378_064_855_777_808_6e-13,
82    3.347_501_239_801_196e-11,
83    2.312_518_412_074_239_3e-9,
84    8.881_290_713_105_89e-8,
85    2.059_622_103_953_427_5e-6,
86    3.055_980_306_089_633_6e-5,
87    3.025_570_258_170_628_5e-4,
88    2.062_051_051_307_880_8e-3,
89    9.903_461_702_320_593e-3,
90    3.410_984_772_609_207e-2,
91    8.534_480_827_208_074e-2,
92    1.565_389_937_575_985e-1,
93    2.117_055_698_804_792_8e-1,
94    2.117_055_698_804_792_8e-1,
95    1.565_389_937_575_985e-1,
96    8.534_480_827_208_074e-2,
97    3.410_984_772_609_207e-2,
98    9.903_461_702_320_593e-3,
99    2.062_051_051_307_880_8e-3,
100    3.025_570_258_170_628_5e-4,
101    3.055_980_306_089_633_6e-5,
102    2.059_622_103_953_427_5e-6,
103    8.881_290_713_105_89e-8,
104    2.312_518_412_074_239_3e-9,
105    3.347_501_239_801_196e-11,
106    2.378_064_855_777_808_6e-13,
107    6.755_290_223_670_204e-16,
108    5.208_449_591_960_913e-19,
109    4.124_607_489_018_283_5e-23,
110];
111
112/// Merton log-utility portfolio with log-normal jumps.
113///
114/// See the module-level documentation for the formulation and numerical
115/// solution.
116#[derive(Clone, Copy, Debug)]
117pub struct MertonJumpLogNormal {
118    /// Risk-free rate.
119    pub risk_free_rate: f64,
120    /// Risky-asset drift.
121    pub mu: f64,
122    /// Risky-asset volatility.
123    pub sigma: f64,
124    /// Poisson jump intensity.
125    pub jump_rate: f64,
126    /// Mean of the log jump size, `E[ln Y]`.
127    pub jump_mean: f64,
128    /// Volatility of the log jump size, `std(ln Y)`.
129    pub jump_vol: f64,
130}
131
132impl MertonJumpLogNormal {
133    /// Creates a Merton log-normal jump problem.
134    ///
135    /// # Panics
136    ///
137    /// Panics if `sigma` or `jump_vol` is negative, if `jump_rate` is negative,
138    /// or if the first-order condition does not bracket an interior optimum on
139    /// `[0, 1]` (i.e. the chosen parameters drive the policy to a corner).
140    ///
141    /// # Examples
142    ///
143    /// ```
144    /// use solver::models::merton_jump_lognormal::MertonJumpLogNormal;
145    /// let m = MertonJumpLogNormal::new(0.03, 0.08, 0.3, 1.0, -0.02, 0.05);
146    /// assert!(m.jump_vol >= 0.0);
147    /// ```
148    pub fn new(
149        risk_free_rate: f64,
150        mu: f64,
151        sigma: f64,
152        jump_rate: f64,
153        jump_mean: f64,
154        jump_vol: f64,
155    ) -> Self {
156        assert!(sigma > 0.0, "volatility must be positive");
157        assert!(jump_rate >= 0.0, "jump rate must be non-negative");
158        assert!(jump_vol >= 0.0, "jump volatility must be non-negative");
159        let model = Self {
160            risk_free_rate,
161            mu,
162            sigma,
163            jump_rate,
164            jump_mean,
165            jump_vol,
166        };
167        let f0 = model.first_order_condition(0.0);
168        let f1 = model.first_order_condition(1.0);
169        assert!(
170            f0 > 0.0 && f1 < 0.0,
171            "parameters do not yield an interior optimum on [0, 1]"
172        );
173        model
174    }
175
176    /// `E[(Y - 1) / (1 + u (Y - 1))]` via Gauss-Hermite quadrature.
177    fn expected_jump_fraction(&self, u: f64) -> f64 {
178        let mut acc = 0.0;
179        for i in 0..32 {
180            let y = (self.jump_mean + self.jump_vol * GH_NODES[i]).exp();
181            acc += GH_WEIGHTS[i] * (y - 1.0) / (1.0 + u * (y - 1.0));
182        }
183        acc
184    }
185
186    /// `E[ln(1 + u (Y - 1))]` via Gauss-Hermite quadrature.
187    fn expected_log_jump(&self, u: f64) -> f64 {
188        let mut acc = 0.0;
189        for i in 0..32 {
190            let y = (self.jump_mean + self.jump_vol * GH_NODES[i]).exp();
191            acc += GH_WEIGHTS[i] * (1.0 + u * (y - 1.0)).ln();
192        }
193        acc
194    }
195
196    /// First-order condition of the log-utility objective in `u`.
197    fn first_order_condition(&self, u: f64) -> f64 {
198        let a = self.mu - self.risk_free_rate;
199        a - self.sigma * self.sigma * u + self.jump_rate * self.expected_jump_fraction(u)
200    }
201
202    /// Optimal portfolio fraction, bracketed on `[0, 1]` and refined by
203    /// bisection.
204    ///
205    /// # Examples
206    ///
207    /// ```
208    /// use solver::models::merton_jump_lognormal::MertonJumpLogNormal;
209    /// let m = MertonJumpLogNormal::new(0.03, 0.08, 0.3, 1.0, -0.02, 0.05);
210    /// assert!(m.exact_policy() > 0.0 && m.exact_policy() < 1.0);
211    /// ```
212    pub fn exact_policy(&self) -> f64 {
213        let mut lo = 0.0;
214        let mut hi = 1.0;
215        let mut f_lo = self.first_order_condition(lo);
216        for _ in 0..80 {
217            let mid = 0.5 * (lo + hi);
218            let f_mid = self.first_order_condition(mid);
219            if f_lo * f_mid <= 0.0 {
220                hi = mid;
221            } else {
222                lo = mid;
223                f_lo = f_mid;
224            }
225        }
226        0.5 * (lo + hi)
227    }
228
229    /// Closed-form value function `V(t, x)` for remaining horizon `tau = T - t`.
230    ///
231    /// # Examples
232    ///
233    /// ```
234    /// use solver::models::merton_jump_lognormal::MertonJumpLogNormal;
235    /// let m = MertonJumpLogNormal::new(0.03, 0.08, 0.3, 1.0, -0.02, 0.05);
236    /// let v = m.exact_value(100.0, 1.0);
237    /// assert!(v.is_finite());
238    /// ```
239    pub fn exact_value(&self, wealth: f64, tau: f64) -> f64 {
240        let u = self.exact_policy();
241        let a = self.mu - self.risk_free_rate;
242        let growth = self.risk_free_rate + u * a - 0.5 * (u * self.sigma).powi(2)
243            + self.jump_rate * self.expected_log_jump(u);
244        wealth.ln() + growth * tau
245    }
246}
247
248/// The zero-volatility reduction of [`MertonJumpLogNormal`] to the
249/// deterministic-jump [`MertonJump`] with multiplier `y = e^m`.
250impl From<&MertonJumpLogNormal> for MertonJump {
251    fn from(model: &MertonJumpLogNormal) -> Self {
252        MertonJump::new(
253            model.risk_free_rate,
254            model.mu,
255            model.sigma,
256            model.jump_rate,
257            model.jump_mean.exp(),
258        )
259    }
260}
261
262#[cfg(test)]
263mod tests {
264    use super::*;
265
266    fn model() -> MertonJumpLogNormal {
267        MertonJumpLogNormal::new(0.03, 0.08, 0.3, 1.0, -0.02, 0.05)
268    }
269
270    #[test]
271    fn policy_is_interior() {
272        let m = model();
273        let u = m.exact_policy();
274        assert!(u > 0.0 && u < 1.0);
275    }
276
277    #[test]
278    fn policy_satisfies_first_order_condition() {
279        let m = model();
280        let u = m.exact_policy();
281        assert!(m.first_order_condition(u).abs() < 1e-10);
282    }
283
284    #[test]
285    fn zero_volatility_reduces_to_deterministic_jump() {
286        let m = MertonJumpLogNormal::new(0.03, 0.08, 0.3, 1.0, -0.02, 0.0);
287        let det = MertonJump::new(0.03, 0.08, 0.3, 1.0, (-0.02f64).exp());
288        assert!((m.exact_policy() - det.exact_policy()).abs() < 1e-10);
289        assert!((m.exact_value(100.0, 1.0) - det.exact_value(100.0, 1.0)).abs() < 1e-10);
290    }
291
292    #[test]
293    fn jump_reduces_exposure_relative_to_no_jump() {
294        let m = model();
295        let no_jump = crate::models::merton::Merton::new(0.03, 0.08, 0.3);
296        assert!(m.exact_policy() < no_jump.exact_policy());
297    }
298}
solver/models/merton_jump_lognormal.rs

solver/models/
merton_jump_lognormal.rs
