atelier_quant 0.0.12

//! Excitation kernel abstraction for Hawkes processes.
//!
//! The [`ExcitationKernel`] trait decouples the Hawkes process simulation
//! and forecasting machinery from a specific kernel shape.  Any kernel
//! that can express its contribution through a recursive *state* variable
//! — initialised once from history in O(n), then maintained in O(1) —
//! can be plugged into [`HawkesProcess<K>`](super::HawkesProcess).
//!
//! ## Provided implementation
//!
//! [`ExponentialKernel`] — the classical $\phi(t) = \alpha e^{-\beta t}$
//! kernel, whose state is a single `f64` (the recursive auxiliary $A$).

/// Trait abstracting the excitation kernel of a univariate Hawkes process.
///
/// # Associated type
///
/// - **`State`**: the recursive summary sufficient for O(1) intensity
///   evaluation.  For the exponential kernel this is a single `f64`
///   (the auxiliary $A$); for a sum-of-exponentials kernel it would be
///   a `Vec<f64>`, etc.
///
/// # Contract
///
/// Implementors must satisfy the invariant:
///
/// ```text
/// init_state(t, history)            ≡ brute-force kernel sum at t
/// decay_state(s, dt) ; excite(s)    ≡ init_state(t + dt, history ∪ {t + dt})
/// ```
pub trait ExcitationKernel: Clone {
    /// Recursive state sufficient for O(1) intensity evaluation.
    type State: Clone;

    // ── point-wise evaluation (for tests / brute-force) ──────────

    /// Evaluate the kernel at lag $\delta$: $\phi(\delta)$.
    fn evaluate(&self, delta: f64) -> f64;

    /// Integrated kernel from 0 to $s$: $\int_0^s \phi(u)\,du$.
    ///
    /// Used by the compensator and the conditional-mean quadrature.
    fn integrate(&self, s: f64) -> f64;

    // ── recursive state management ───────────────────────────────

    /// Initialise the recursive state from a full event history.
    ///
    /// $$\text{State}(t) = \sum_{t_i < t} (\text{kernel contribution of } t_i)$$
    ///
    /// Cost: O(n) in the history length.
    fn init_state(&self, t: f64, history: &[f64]) -> Self::State;

    /// Decay the state across a time gap $\Delta t$ (no new event).
    fn decay_state(&self, state: &Self::State, dt: f64) -> Self::State;

    /// Excite the state upon acceptance of a new event at the current time.
    fn excite_state(&self, state: &Self::State) -> Self::State;

    /// Intensity contribution from the kernel at the given state:
    /// $\alpha \cdot \text{state\_summary}$.
    fn intensity_contribution(&self, state: &Self::State) -> f64;

    /// Integrated intensity contribution starting from `state` over a
    /// future interval of length $s$:
    ///
    /// $$\int_0^s \alpha \cdot g(\text{state}, u)\,du$$
    ///
    /// where $g$ captures how the *existing* state decays into the
    /// future.  For the exponential kernel this is:
    ///
    /// $$\frac{\alpha}{\beta} \cdot A \cdot (1 - e^{-\beta s})$$
    ///
    /// This is the key ingredient for the deterministic conditional-mean
    /// forecast: the survival function of the next gap is
    ///
    /// $$S(s) = \exp\!\bigl(-\mu s - \text{integrated\_intensity}(s)\bigr)$$
    fn integrated_intensity_contribution(&self, state: &Self::State, s: f64) -> f64;
}

// ═══════════════════════════════════════════════════════════════════
//  Exponential kernel:  φ(t) = α · exp(−β·t)
// ═══════════════════════════════════════════════════════════════════

/// Classical exponential excitation kernel.
///
/// $$\phi(t) = \alpha \, e^{-\beta t}$$
///
/// The recursive state is a single `f64`:
///
/// $$A(t) = \sum_{t_i < t} e^{-\beta(t - t_i)}$$
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct ExponentialKernel {
    /// Excitation factor $\alpha$ (jump size per event). Must be ≥ 0.
    pub alpha: f64,
    /// Exponential decay rate $\beta$. Must be > 0.
    pub beta: f64,
}

impl ExcitationKernel for ExponentialKernel {
    type State = f64;

    #[inline]
    fn evaluate(&self, delta: f64) -> f64 {
        self.alpha * (-self.beta * delta).exp()
    }

    #[inline]
    fn integrate(&self, s: f64) -> f64 {
        // ∫₀ˢ α exp(−β u) du = (α/β)(1 − exp(−β s))
        (self.alpha / self.beta) * (1.0 - (-self.beta * s).exp())
    }

    fn init_state(&self, t: f64, history: &[f64]) -> f64 {
        let mut a = 0.0_f64;
        for &ti in history {
            if ti < t {
                a += (-self.beta * (t - ti)).exp();
            }
        }
        a
    }

    #[inline]
    fn decay_state(&self, state: &f64, dt: f64) -> f64 {
        state * (-self.beta * dt).exp()
    }

    #[inline]
    fn excite_state(&self, state: &f64) -> f64 {
        state + 1.0
    }

    #[inline]
    fn intensity_contribution(&self, state: &f64) -> f64 {
        self.alpha * state
    }

    #[inline]
    fn integrated_intensity_contribution(&self, state: &f64, s: f64) -> f64 {
        // ∫₀ˢ α · A · exp(−β u) du = (α/β) · A · (1 − exp(−β s))
        (self.alpha / self.beta) * state * (1.0 - (-self.beta * s).exp())
    }
}