solver/models/

merton.rs

//! Merton's log-utility portfolio problem.
//!
//! A validating instance of [`ControlProblem`] with a known closed-form
//! solution. The state is wealth `x`; the control `u in [0, 1]` is the
//! fraction of wealth allocated to a risky asset.
//!
//! # Dynamics
//!
//! ```text
//! dS/S = mu dt + sigma dW        (risky asset)
//! r     = risk-free rate          (riskless asset)
//! dx    = x [ r + u (mu - r) ] dt + x u sigma dW
//! ```
//!
//! # Objective
//!
//! Maximize expected log utility of terminal wealth:
//!
//! ```text
//! 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 }
//! ```
//!
//! # Exact solution
//!
//! The optimal portfolio fraction is constant:
//!
//! ```text
//! u* = (mu - r) / sigma^2
//! ```
//!
//! and the value function is
//!
//! ```text
//! V(t, x) = ln x + [ r + 0.5 (mu - r)^2 / sigma^2 ] (T - t)
//! ```
//!
//! where `T - t` is the remaining horizon.
use crate::models::control::{ControlProblem, StateDerivatives};
use crate::numeric::finite_difference::discretization::{DimensionKind, Transport};
use crate::numeric::finite_difference::pde::PdeProblem;

/// Merton log-utility portfolio problem.
///
/// See the module-level documentation for the formulation and exact solution.
#[derive(Clone, Copy, Debug)]
pub struct Merton {
    /// Risk-free rate.
    pub risk_free_rate: f64,
    /// Risky-asset drift.
    pub mu: f64,
    /// Risky-asset volatility.
    pub sigma: f64,
}

impl Merton {
    /// Creates a Merton problem.
    ///
    /// # Panics
    ///
    /// Panics if `sigma` is not positive.
    ///
    /// # Examples
    ///
    /// ```
    /// use solver::models::merton::Merton;
    /// let m = Merton::new(0.03, 0.08, 0.2);
    /// assert!((m.sigma - 0.2).abs() < 1e-12);
    /// ```
    pub fn new(risk_free_rate: f64, mu: f64, sigma: f64) -> Self {
        assert!(sigma > 0.0, "volatility must be positive");
        Self {
            risk_free_rate,
            mu,
            sigma,
        }
    }

    /// Optimal (constant) portfolio fraction.
    ///
    /// # Examples
    ///
    /// ```
    /// use solver::models::merton::Merton;
    /// let m = Merton::new(0.03, 0.08, 0.2);
    /// assert!((m.exact_policy() - 0.05 / 0.04).abs() < 1e-12);
    /// ```
    pub fn exact_policy(&self) -> f64 {
        (self.mu - self.risk_free_rate) / (self.sigma * self.sigma)
    }

    /// Closed-form value function `V(t, x)` for remaining horizon `tau = T - t`.
    ///
    /// # Examples
    ///
    /// ```
    /// use solver::models::merton::Merton;
    /// let m = Merton::new(0.03, 0.08, 0.2);
    /// 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 growth = self.risk_free_rate + 0.5 * (self.mu - self.risk_free_rate) * u;
        wealth.ln() + growth * tau
    }
}

impl ControlProblem<1> for Merton {
    type Control = f64;

    fn optimize(&self, _t: f64, state: &[f64; 1], _derivs: &StateDerivatives<1>) -> Self::Control {
        // The optimal fraction is independent of wealth and time. Clamp to a
        // non-negative admissible value so the policy stays well defined even
        // for extrapolated states.
        let _ = state;
        self.exact_policy().max(0.0)
    }

    fn running_reward(&self, _t: f64, _state: &[f64; 1], _control: &Self::Control) -> f64 {
        // No running reward: log-utility has a Mayer (terminal-only) objective.
        0.0
    }

    fn generator(
        &self,
        _t: f64,
        state: &[f64; 1],
        control: &Self::Control,
        derivs: &StateDerivatives<1>,
    ) -> f64 {
        let x = state[0];
        let vx = derivs.grad[0];
        let vxx = derivs.hessian[0];
        let u = *control;

        x * (self.risk_free_rate + u * (self.mu - self.risk_free_rate)) * vx
            + 0.5 * (x * u * self.sigma).powi(2) * vxx
    }

    fn terminal(&self, state: &[f64; 1]) -> f64 {
        state[0].ln()
    }

    fn next_step(&self, _t: f64, state: &[f64; 1], dt: f64, noise: &[f64; 1]) -> [f64; 1] {
        let x = state[0];
        let u = self.exact_policy().max(0.0);
        let drift = x * (self.risk_free_rate + u * (self.mu - self.risk_free_rate)) * dt;
        let diffusion = x * u * self.sigma * dt.sqrt() * noise[0];
        [x + drift + diffusion]
    }

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

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

    fn transport(
        &self,
        _t: f64,
        state: &[f64; 1],
        control: &Self::Control,
        _derivs: &StateDerivatives<1>,
    ) -> Transport<1> {
        // dx = x [r + u (mu - r)] dt + x u sigma dW.
        let x = state[0];
        let u = *control;
        let b = x * (self.risk_free_rate + u * (self.mu - self.risk_free_rate));
        let d = 0.5 * (x * u * self.sigma).powi(2);
        let h = x.max(1e-12);
        Transport::new(
            [d / (h * h) + b.max(0.0) / h],
            [d / (h * h) + (-b).max(0.0) / h],
            0.0,
        )
    }
}

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

    fn model() -> Merton {
        Merton::new(0.03, 0.08, 0.2)
    }

    #[test]
    fn exact_policy_matches_closed_form() {
        let m = model();
        assert!((m.exact_policy() - 0.05 / 0.04).abs() < 1e-12);
    }

    #[test]
    fn optimize_returns_exact_policy() {
        let m = model();
        let d = StateDerivatives::new([0.0], [0.0]);
        let u = m.optimize(0.0, &[100.0], &d);
        assert!((u - m.exact_policy()).abs() < 1e-12);
    }

    #[test]
    fn terminal_is_log_wealth() {
        let m = model();
        assert!((m.terminal(&[100.0]) - 100.0f64.ln()).abs() < 1e-12);
    }

    #[test]
    fn driver_matches_hjb_form_at_exact_derivatives() {
        // For V = ln x + c tau, V_x = 1/x and V_xx = -1/x^2.
        let m = model();
        let x = 100.0;
        let u = m.exact_policy();
        let d = StateDerivatives::new([1.0 / x], [-1.0 / (x * x)]);

        let expected =
            u * (m.mu - m.risk_free_rate) + m.risk_free_rate - 0.5 * u * u * m.sigma * m.sigma;
        assert!((m.driver(0.0, &[x], &u, &d) - expected).abs() < 1e-12);
    }

    #[test]
    fn next_step_preserves_positive_wealth_for_small_dt() {
        let m = model();
        let s = m.next_step(0.0, &[100.0], 0.001, &[0.0]);
        assert!(s[0] > 0.0);
    }

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