Univariate completed

Changed file structure so each type of Hawkes process is in its own
folder.
`Intensity`, `simulation`, `time_change` and `fit` implemented for univariate
processes.
Missing `time_change` and `fit` for multivariate processes.

Author: JoseKling

docs/src/api.md

@@ -107,8 +107,8 @@ HawkesProcess
 ### Univariate
 
 ```@docs
-UnivariateHawkesProcess
 UnmarkedUnivariateHawkesProcess
+UnivariateHawkesProcess
 ```
 
 ### Multivariate

src/PointProcesses.jl

@@ -10,7 +10,7 @@ module PointProcesses
 using DensityInterface: DensityInterface, HasDensity, densityof, logdensityof
 using Distributions: Distributions, UnivariateDistribution, MultivariateDistribution
 using Distributions: Categorical, Exponential, Poisson, Uniform, Dirac
-using Distributions: fit, suffstats, probs, mean, support
+using Distributions: fit, suffstats, probs, mean, support, pdf
 using LinearAlgebra: dot, diagm
 using Random: rand
 using Random: AbstractRNG, default_rng
@@ -73,6 +73,20 @@ include("hawkes/hawkes_process.jl")
 include("hawkes/simulation.jl")
 include("hawkes/intensity.jl")
 include("hawkes/fit.jl")
-include("hawkes/time_change.jl")
+include("hawkes/unmarked_univariate/unmarked_univariate_hawkes.jl")
+include("hawkes/unmarked_univariate/intensity.jl")
+include("hawkes/unmarked_univariate/time_change.jl")
+include("hawkes/unmarked_univariate/fit.jl")
+include("hawkes/unmarked_univariate/simulation.jl")
+include("hawkes/univariate/univariate_hawkes.jl")
+include("hawkes/univariate/intensity.jl")
+include("hawkes/univariate/time_change.jl")
+include("hawkes/univariate/fit.jl")
+include("hawkes/univariate/simulation.jl")
+include("hawkes/multivariate/multivariate_hawkes.jl")
+include("hawkes/multivariate/intensity.jl")
+# include("hawkes/multivariate/time_change.jl")
+# include("hawkes/multivariate/fit.jl")
+include("hawkes/multivariate/simulation.jl")
 
 end

src/hawkes/fit.jl

@@ -1,50 +1,29 @@
-"""
-    StatsAPI.fit(rng, ::Type{UnmarkedUnivariateHawkesProcess{T}}, h::History; step_tol::Float64 = 1e-6, max_iter::Int = 1000) where {T<:Real}
-
-Expectation-Maximization algorithm from [E. Lewis, G. Mohler (2011)](https://arxiv.org/pdf/1801.08273)).
-The relevant calculations are in page 4, equations 6-13.
-
-Let (t₁ < ... < tₙ) be the event times over the interval [0, T). We use the immigrant-descendant representation,
-where immigrants arrive at a constant base rate μ and each each arrival may generate descendants following the
-activation function α exp(-ω(t - tᵢ)).
-
-The algorithm consists in the following steps:
-1. Start with some initial guess for the parameters μ, ψ, and ω. ψ = α ω is the branching factor.
-2. Calculate λ(tᵢ; μ, ψ, ω) (`lambda_ts` in the code) using the procedure in [Ozaki (1979)](https://doi.org/10.1007/bf02480272).
-3. For each tᵢ and each j < i, calculate Dᵢⱼ = P(tᵢ is a descendant of tⱼ) as
-
-    Dᵢⱼ = ψ ω exp(-ω(tᵢ - tⱼ)) / λ(tᵢ; μ, ψ, ω).
-
-    Define D = ∑_{j < i} Dᵢⱼ (expected number of descendants) and div = ∑_{j < i} (tᵢ - tⱼ) Dᵢⱼ. 
-4. Update the parameters as
-        μ = (N - D) / T
-        ψ = D / N
-        ω = D / div
-5. If convergence criterion is met, return updated parameters, otherwise, back to step 2.
-
-Notice that, in the implementation, the process is normalized so the average inter-event time is equal to 1 and, 
-therefore, the interval of the process is transformed from T to N. Also, in equation (8) in the paper,
-
-∑_{i=1:n} pᵢᵢ = ∑_{i=1:n} (1 - ∑_{j < i} Dᵢⱼ) = N - D.
-
-"""
-function StatsAPI.fit(
-    ::Type{<:UnmarkedUnivariateHawkesProcess{T}},
-    h::History;
+#=
+Method for fitting the parameters of marked and unmarked Hawkes
+processes. For unmarked processes, the marks in `history` must
+be either all equal to `nothing` or equal to 1.
+=#
+function fit_hawkes_params(
+    ::Type{R},
+    h::History,
+    div_ψ::Real;
     step_tol::Float64=1e-6,
     max_iter::Int=1000,
     rng::AbstractRNG=default_rng(),
-) where {T<:Real}
+) where {R<:Real}
+    T = float(R)
+    div_ψ = T(div_ψ)
     n = nb_events(h)
-    n == 0 && return HawkesProcess(zero(T), zero(T), zero(T))
+    nT = T(n)
+    n == 0 && return (zero(T), zero(T), zero(T))
 
     tmax = T(duration(h))
     # Normalize times so average inter-event time is 1 (T -> n)
     norm_ts = T.(h.times .* (n / tmax))
 
     # preallocate
-    A = zeros(T, n)            # A[i] = sum_{j<i} exp(-ω (t_i - t_j))
-    S = zeros(T, n)            # S[i] = sum_{j<i} (t_i - t_j) exp(-ω (t_i - t_j))
+    A = zeros(T, n)            # A[i] = sum_{j<i} mⱼ exp(-ω (t_i - t_j))
+    S = zeros(T, n)            # S[i] = sum_{j<i} mⱼ (t_i - t_j) exp(-ω (t_i - t_j))
     lambda_ts = similar(A)     # λ_i
 
     # Λ(T) ≈ n after normalization
@@ -65,9 +44,8 @@ function StatsAPI.fit(
         for i in 2:n
             Δ = norm_ts[i] - norm_ts[i - 1]
             e = exp(-ω * Δ)
-            Ai_1 = A[i - 1]
-            A[i] = e * (one(T) + Ai_1)
-            S[i] = e * (S[i - 1] + Δ * (one(T) + Ai_1))
+            A[i] = update_A(A[i - 1], e, h.marks[i - 1]) # Different updates for real marks and no marks
+            S[i] = (Δ * A[i]) + (e * S[i - 1])
             lambda_ts[i] = μ + (ψ * ω) * A[i]
         end
 
@@ -81,8 +59,8 @@ function StatsAPI.fit(
         end
 
         # Steps 4–5: M-step updates + convergence check
-        new_μ = one(T) - (D / n)
-        new_ψ = D / n
+        new_μ = one(T) - (D / nT)
+        new_ψ = D / div_ψ
         new_ω = D / (div + eps(T))   # small guard to avoid div=0
 
         step = max(abs(μ - new_μ), abs(ψ - new_ψ), abs(ω - new_ω))
@@ -97,13 +75,5 @@ function StatsAPI.fit(
 
     # Unnormalize back to original time scale (T -> tmax):
     # parameters in normalized space (') relate to original by μ0=μ'*(n/tmax), ω0=ω'*(n/tmax), α0=(ψ'ω')*(n/tmax)
-    return HawkesProcess(μ * (n / tmax), ψ * ω * (n / tmax), ω * (n / tmax))
-end
-
-# Type parameter for `HawkesProcess` was NOT explicitly provided
-function StatsAPI.fit(
-    HP::Type{UnmarkedUnivariateHawkesProcess}, h::History{H,M}; kwargs...
-) where {H<:Real,M}
-    T = promote_type(Float64, H)
-    return fit(HP{T}, h; kwargs...)
+    return (μ * (nT / tmax), ψ * ω * (nT / tmax), ω * (nT / tmax))
 end

src/hawkes/hawkes_process.jl

@@ -5,183 +5,13 @@ Common interface for all subtypes of `HawkesProcess`.
 """
 abstract type HawkesProcess <: AbstractPointProcess end
 
-"""
-    UnivariateHawkesProcess{R<:Real,D} <: HawkesProcess
-
-Univariate Hawkes process with exponential decay kernel and mark 
-distribution `D`.
-
-Denote the events of the process by (tᵢ, mᵢ), where tᵢ is the event time
-and mᵢ ∈ M the corresponding mark. The conditional intensity function
-of the Hawkes process is given by
-
-λ(t) = μ + ∑_{tᵢ < t} α mᵢ exp(-ω (t - tᵢ)).
-
-Notice that the mark only affects the jump size.
-
-# Fields
-- `μ::R`: baseline intensity (immigration rate)
-- `α::R`: jump size
-- `ω::R`: decay rate.
-- `mark_dist::D`: distribution of marks
-"""
-struct UnivariateHawkesProcess{T<:Real,D} <: HawkesProcess
-    μ::T
-    α::T
-    ω::T
-    mark_dist::D
-end
-
-function Base.show(io::IO, hp::UnivariateHawkesProcess)
-    return print(io, "UnivariateHawkesProcess($(hp.μ), $(hp.α), $(hp.ω), $(hp.mark_dist))")
-end
-
-"""
-    UnmarkedUnivariateHawkesProcess{R<:Real} <: HawkesProcess
-
-Unmarked univariate Hawkes process with exponential decay kernel.
-
-Denote the events of the process by tᵢ. The conditional intensity function
-of the Hawkes process is given by
-
-λ(t) = μ + ∑_{tᵢ < t} α exp(-ω (t - tᵢ)).
-
-# Fields
-- `μ::R`: baseline intensity (immigration rate)
-- `α::R`: jump size
-- `ω::R`: decay rate.
-
-Alias for UnivariateHawkesProcess{R, Dirac{Nothing}}.
-"""
-const UnmarkedUnivariateHawkesProcess{R<:Real} = UnivariateHawkesProcess{R,Dirac{Nothing}}
-
-function Base.show(io::IO, hp::UnmarkedUnivariateHawkesProcess)
-    return print(io, "UnmarkedUnivariateHawkesProcess$((hp.μ, hp.α, hp.ω))")
-end
-
-"""
-    MultivariateHawkesProcess{R} <: HawkesProcess
-
-Unmarked multivariate temporal Hawkes process with exponential decay
-
-For a process with m = 1, 2, ..., M marginal processes, the conditional intensity
-function of the m-th process is given by
-
-    λₘ(t) = μₘ + ∑_{l=1,...,M} ∑_{tˡᵢ < t} αₘₗ exp(-βₘₗ(t - tˡᵢ)),
-
-where tˡᵢ is the i-th element in the l-th marginal process.
-μ is a vector of length M with elements μₘ. α and ω are M×M
-matrices with elements αₘₗ and ωₘₗ.
-
-The process is represented as a marked process, where each marginal
-process m = 1, 2, ..., M is represented by events (tᵐᵢ, m), tᵐᵢ being
-the i-th element with mark m.
-
-# Fields
-- `μ::Vector{<:Real}`: baseline intensity (immigration rate)
-- `α::Matrix{<:Real}`: Jump size.
-- `ω::Matrix{<:Real}`: Decay rate.
-"""
-struct MultivariateHawkesProcess{T<:Real} <: HawkesProcess
-    μ::T
-    α::Matrix{T}
-    ω::Matrix{T}
-    mark_dist::Categorical{Float64,Vector{Float64}} # To keep consistent with PoissonProcess. Helps in simulation.
-end
-
-function Base.show(io::IO, hp::MultivariateHawkesProcess{T}) where {T<:Real}
-    return print(
-        io,
-        "MultivariateHawkesProcess\nμ = $(T.(hp.μ .* probs(hp.mark_dist)))\nα = $(hp.α)\nω = $(hp.ω)",
-    )
-end
-
-Base.length(mh::MultivariateHawkesProcess) = size(mh.α)[1]
-
-## Constructors
-### UnivariatehawkesProcess
-function HawkesProcess(μ::Real, α::Real, ω::Real, mark_dist; check_args::Bool=true)
-    check_args && check_args_Hawkes(μ, α, ω, mark_dist)
-    return UnivariateHawkesProcess(promote(μ, α, ω)..., mark_dist)
-end
-
-function HawkesProcess(μ::Real, α::Real, ω::Real; check_args::Bool=true)
-    return HawkesProcess(μ, α, ω, Dirac(nothing); check_args=check_args)
-end
-
-### MultivariateHawkesProcess
-function HawkesProcess(
-    μ::Vector{<:Real}, α::Matrix{<:Real}, ω::Matrix{<:Real}; check_args::Bool=true
-)
-    check_args && check_args_Hawkes(μ, α, ω)
-    R = promote_type(eltype(μ), eltype(α), eltype(ω))
-    return MultivariateHawkesProcess(R(sum(μ)), R.(α), R.(ω), Categorical(μ / sum(μ)))
-end
-
-# For M independent marginal processes
-function HawkesProcess(
-    μ::Vector{<:Real}, α::Vector{<:Real}, ω::Vector{<:Real}; check_args::Bool=true
-)
-    check_args && check_args_Hawkes(μ, diagm(α), diagm(ω))
-    R = promote_type(eltype(μ), eltype(α), eltype(ω))
-    return MultivariateHawkesProcess(R(sum(μ)), R.(α), R.(ω), Categorical(μ / sum(μ)))
-end
-
-# Check args
-## Univariate Hawkes
-function check_args_Hawkes(μ::Real, α::Real, ω::Real, mark_dist)
-    if any((μ, α, ω) .< 0)
-        throw(
-            DomainError(
-                "μ = $μ, α = $α, ω = $ω",
-                "HawkesProcess: All parameters must be non-negative.",
-            ),
-        )
-    end
-    mean_α = mark_dist isa Dirac{Nothing} ? α : mean(mark_dist) * α
-    if mean_α >= ω
-        throw(
-            DomainError(
-                "α = $(mean_α), ω = $ω",
-                "HawkesProcess: mᵢα must be, on average, smaller than ω. Stability condition.",
-            ),
-        )
-    end
-    if !isa(mark_dist, Dirac{Nothing}) && minimum(mark_dist) < 0
-        throw(
-            DomainError(
-                "Mark distribution support = $((support(mark_dist).lb, support(mark_dist).ub))",
-                "HawkesProcess: Support of mark distribution must be contained in non-negative numbers",
-            ),
-        )
-    end
-    return nothing
-end
-
-## Multivariate Hawkes
-function check_args_Hawkes(μ::Vector{<:Real}, α::Matrix{<:Real}, ω::Matrix{<:Real})
-    if length(μ) != size(α)[2]
-        throw(DimansionMismatch("α must have size $(length(μ))×$(length(μ))"))
-    end
-    if length(μ) != size(ω)[2]
-        throw(DimansionMismatch("ω must have size $(length(μ))×$(length(μ))"))
-    end
-    if any(sum(α ./ ω; dims=1) .>= 1)
-        throw(
-            DomainError(
-                "α = $α, ω = $ω",
-                "HawkesProcess: Sum of α/β over each row must be smaller than 1. Stability condition.",
-            ),
-        )
-    end
-    if any(μ .< zero(μ)) || any(α .< zero(α)) || any(ω .< zero(ω))
-        throw(
-            DomainError(
-                "μ = $μ, α = $α, ω = $ω",
-                "HawkesProcess: All elements of μ, α and ω must be non-negative.",
-            ),
-        )
-    end
-    ω[α .== 0] .= 1 # Protects against division by 0 in simulation
-    return nothing
-end
+# #=
+# If type parameter for `HawkesProcess` was NOT explicitly provided,
+# use `Float64` as the standard type
+# =#
+# function StatsAPI.fit(
+#     HP::Type{<:HawkesProcess}, h::History{H,M}; kwargs...
+# ) where {H<:Real,M}
+#     T = promote_type(Float64, H)
+#     return fit(HP{T}, h; kwargs...)
+# end