solver/models/

merton_jump.rs

//! Merton's log-utility portfolio problem with a jump in the risky asset.
//!
//! A validating instance for the neural jump BSDE: a diffusion-plus-jump
//! wealth process whose optimal control and value still have closed forms, but
//! where the jump term genuinely changes the optimal policy (unlike the
//! no-jump [`Merton`], where the policy is constant in
//! the jump parameters).
//!
//! # Dynamics
//!
//! The risky asset follows a geometric jump-diffusion and the wealth `x` is
//! invested with fraction `u` in the risky asset and the rest at the risk-free
//! rate `r`:
//!
//! ```text
//! dS/S  = mu dt + sigma dW + (y - 1) dN
//! dx    = x [ r + u (mu - r) ] dt + x u sigma dW + x u (y - 1) dN
//! ```
//!
//! where `N` is a Poisson process with intensity `lambda` and `y > 0` is a
//! deterministic multiplicative jump, so a jump scales wealth by
//! `1 + u (y - 1)`.
//!
//! # Objective
//!
//! Maximize expected log utility of terminal wealth, `J(u) = E[ ln x_T ]`.
//!
//! # HJB equation
//!
//! ```text
//! 0 = d_t V + sup_u { x (r + u (mu - r)) V_x + 0.5 (x u sigma)^2 V_xx
//!      + lambda [ V(x (1 + u (y - 1))) - V(x) ] }
//! ```
//!
//! # Exact solution
//!
//! The value function is `V(t, x) = ln x + B (T - t)` with
//!
//! ```text
//! B = r + u* (mu - r) - 0.5 (u* sigma)^2 + lambda ln(1 + u* (y - 1)),
//! ```
//!
//! and the optimal fraction `u*` solves the quadratic
//!
//! ```text
//! 0 = (mu - r) - sigma^2 u* + lambda (y - 1) / (1 + u* (y - 1)).
//! ```
//!
//! Setting `lambda = 0` or `y = 1` reduces this to the no-jump Merton policy
//! `u* = (mu - r) / sigma^2`.
use crate::models::merton::Merton;

/// Merton log-utility portfolio with a deterministic multiplicative jump.
///
/// See the module-level documentation for the formulation and exact solution.
#[derive(Clone, Copy, Debug)]
pub struct MertonJump {
    /// Risk-free rate.
    pub risk_free_rate: f64,
    /// Risky-asset drift.
    pub mu: f64,
    /// Risky-asset volatility.
    pub sigma: f64,
    /// Poisson jump intensity.
    pub jump_rate: f64,
    /// Deterministic multiplicative jump (wealth scales by `1 + u (y - 1)`).
    pub jump_multiplier: f64,
}

impl MertonJump {
    /// Creates a Merton jump problem.
    ///
    /// # Panics
    ///
    /// Panics if `sigma` or `jump_multiplier` is not positive, or if
    /// `jump_rate` is negative.
    ///
    /// # Examples
    ///
    /// ```
    /// use solver::models::merton_jump::MertonJump;
    /// let m = MertonJump::new(0.03, 0.08, 0.3, 1.0, 1.1);
    /// assert_eq!(m.jump_multiplier, 1.1);
    /// ```
    pub fn new(
        risk_free_rate: f64,
        mu: f64,
        sigma: f64,
        jump_rate: f64,
        jump_multiplier: f64,
    ) -> Self {
        assert!(sigma > 0.0, "volatility must be positive");
        assert!(jump_rate >= 0.0, "jump rate must be non-negative");
        assert!(jump_multiplier > 0.0, "jump multiplier must be positive");
        Self {
            risk_free_rate,
            mu,
            sigma,
            jump_rate,
            jump_multiplier,
        }
    }

    /// Optimal (constant) portfolio fraction, the admissible root of the
    /// policy quadratic.
    ///
    /// # Examples
    ///
    /// ```
    /// use solver::models::merton_jump::MertonJump;
    /// let m = MertonJump::new(0.03, 0.08, 0.3, 1.0, 1.1);
    /// assert!(m.exact_policy() > 0.0);
    /// ```
    pub fn exact_policy(&self) -> f64 {
        let a = self.mu - self.risk_free_rate;
        let s = self.jump_multiplier - 1.0;
        let sig2 = self.sigma * self.sigma;

        // A zero jump amplitude degenerates to the no-jump Merton policy.
        if s.abs() < 1e-14 {
            return a / sig2;
        }

        // The first-order condition
        //   a - sig2 u + lambda s / (1 + u s) = 0
        // is a quadratic in u after multiplying by (1 + u s).
        let quad_a = -sig2 * s;
        let quad_b = a * s - sig2;
        let quad_c = a + self.jump_rate * s;
        // The discriminant is (a s + sig2)^2 + 4 lambda sig2 s^2 >= 0.
        let disc = quad_b * quad_b - 4.0 * quad_a * quad_c;
        let sqrt_disc = disc.max(0.0).sqrt();
        let r1 = (-quad_b + sqrt_disc) / (2.0 * quad_a);
        let r2 = (-quad_b - sqrt_disc) / (2.0 * quad_a);

        // Exactly one root keeps wealth positive after a jump (1 + u s > 0).
        if 1.0 + r1 * s > 0.0 { r1 } else { r2 }
    }

    /// Closed-form value function `V(t, x)` for remaining horizon `tau = T - t`.
    ///
    /// # Examples
    ///
    /// ```
    /// use solver::models::merton_jump::MertonJump;
    /// let m = MertonJump::new(0.03, 0.08, 0.3, 1.0, 1.1);
    /// let v = m.exact_value(100.0, 1.0);
    /// assert!(v.is_finite());
    /// ```
    pub fn exact_value(&self, wealth: f64, tau: f64) -> f64 {
        let u = self.exact_policy();
        let a = self.mu - self.risk_free_rate;
        let s = self.jump_multiplier - 1.0;
        let growth = self.risk_free_rate + u * a - 0.5 * (u * self.sigma).powi(2)
            + self.jump_rate * (1.0 + u * s).ln();
        wealth.ln() + growth * tau
    }
}

/// The no-jump reduction of [`MertonJump`] with `jump_rate = 0`.
impl From<&MertonJump> for Merton {
    fn from(model: &MertonJump) -> Self {
        Merton::new(model.risk_free_rate, model.mu, model.sigma)
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    fn model() -> MertonJump {
        MertonJump::new(0.03, 0.08, 0.3, 1.0, 1.1)
    }

    #[test]
    fn zero_jump_rate_reduces_to_no_jump_policy() {
        let m = MertonJump::new(0.03, 0.08, 0.3, 0.0, 1.1);
        let no_jump = Merton::new(0.03, 0.08, 0.3);
        assert!((m.exact_policy() - no_jump.exact_policy()).abs() < 1e-12);
    }

    #[test]
    fn unit_multiplier_reduces_to_no_jump_policy() {
        let m = MertonJump::new(0.03, 0.08, 0.3, 1.0, 1.0);
        let no_jump = Merton::new(0.03, 0.08, 0.3);
        assert!((m.exact_policy() - no_jump.exact_policy()).abs() < 1e-12);
    }

    #[test]
    fn policy_satisfies_first_order_condition() {
        let m = model();
        let u = m.exact_policy();
        let a = m.mu - m.risk_free_rate;
        let s = m.jump_multiplier - 1.0;
        let residual = a - m.sigma * m.sigma * u + m.jump_rate * s / (1.0 + u * s);
        assert!(residual.abs() < 1e-10);
    }

    #[test]
    fn jump_changes_the_optimal_policy() {
        let m = model();
        let no_jump = Merton::new(0.03, 0.08, 0.3);
        assert!((m.exact_policy() - no_jump.exact_policy()).abs() > 1e-3);
    }

    #[test]
    fn value_is_finite_and_positive_growth() {
        let m = model();
        let v = m.exact_value(100.0, 1.0);
        assert!(v.is_finite());
        assert!(v > 100.0f64.ln());
    }
}