diff --git a/docs/src/api.md b/docs/src/api.md index e55431b..ce3fa8b 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -99,6 +99,18 @@ MultivariatePoissonProcessPrior HawkesProcess ``` +### Univariate + +```@docs +UnivariateHawkesProcess +``` + +### Multivariate + +```@docs +MultivariateHawkesProcess +``` + ## Goodness-of-fit tests ```@docs diff --git a/src/PointProcesses.jl b/src/PointProcesses.jl index e91597b..e62804a 100644 --- a/src/PointProcesses.jl +++ b/src/PointProcesses.jl @@ -10,8 +10,8 @@ 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 -using LinearAlgebra: dot +using Distributions: fit, suffstats, probs, mean, support, pdf +using LinearAlgebra: dot, diagm using Random: rand using Random: AbstractRNG, default_rng, Xoshiro using StatsAPI: StatsAPI, fit, HypothesisTest, pvalue @@ -47,10 +47,15 @@ export simulate_ogata, simulate ## Models +### Poisson processes export PoissonProcess export UnivariatePoissonProcess export MultivariatePoissonProcess, MultivariatePoissonProcessPrior + +### Hawkes processes export HawkesProcess +export UnivariateHawkesProcess, UnmarkedUnivariateHawkesProcess +export MultivariateHawkesProcess ## Goodness of fit tests tests @@ -71,6 +76,16 @@ include("poisson/fit.jl") include("poisson/simulation.jl") include("hawkes/hawkes_process.jl") +include("hawkes/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/simulation.jl") include("HypothesisTests/point_process_tests.jl") include("HypothesisTests/statistic.jl") diff --git a/src/hawkes/hawkes_process.jl b/src/hawkes/hawkes_process.jl index a628c17..633c5f5 100644 --- a/src/hawkes/hawkes_process.jl +++ b/src/hawkes/hawkes_process.jl @@ -1,238 +1,42 @@ """ - HawkesProcess{T<:Real} - -Univariate Hawkes process with exponential decay kernel. - -A Hawkes process is a self-exciting point process where each event increases the probability -of future events. The conditional intensity function is given by: - - λ(t) = μ + α ∑_{tᵢ < t} exp(-ω(t - tᵢ)) - -where the sum is over all previous event times tᵢ. - -# Fields - -- `μ::T`: baseline intensity (immigration rate) -- `α::T`: jump size (immediate increase in intensity after an event) -- `ω::T`: decay rate (how quickly the excitement fades) - -Conditions: -- μ, α, ω >= 0 -- ψ = α/ω < 1 → Stability condition. ψ is the expected number of events each event generates - -Following the notation from [E. Lewis, G. Mohler (2011)](https://arxiv.org/pdf/1801.08273). -""" -struct HawkesProcess{T<:Real} <: AbstractPointProcess - μ::T - α::T - ω::T - - function HawkesProcess(μ::T1, α::T2, ω::T3) where {T1,T2,T3} - any((μ, α, ω) .< 0) && - throw(DomainError((μ, α, ω), "All parameters must be non-negative.")) - (α > 0 && α >= ω) && - throw(DomainError((α, ω), "Parameter ω must be strictly smaller than α")) - T = promote_type(T1, T2, T3) - (μ_T, α_T, ω_T) = convert.(T, (μ, α, ω)) - new{T}(μ_T, α_T, ω_T) - end -end - -function simulate(rng::AbstractRNG, hp::HawkesProcess, tmin, tmax) - sim = simulate_poisson_times(rng, hp.μ, tmin, tmax) # Simulate Poisson process with base rate - sim_desc = generate_descendants(rng, sim, tmax, hp.α, hp.ω) # Recursively generates descendants from first events - append!(sim, sim_desc) - sort!(sim) - return History(; times=sim, tmin=tmin, tmax=tmax, check=false) -end + HawkesProcess <: AbstractPointProcess +Common interface for all subtypes of `HawkesProcess`. """ - StatsAPI.fit(rng, ::Type{HawkesProcess{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{HawkesProcess{T}}, - h::History; - step_tol::Float64=1e-6, - max_iter::Int=1000, - rng::AbstractRNG=default_rng(), -) where {T<:Real} - n = nb_events(h) - n == 0 && return HawkesProcess(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= step_tol) && (n_iters < max_iter) - # compute A, S, and λ - 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)) - lambda_ts[i] = μ + (ψ * ω) * A[i] - end - - # E-step aggregates - D = zero(T) # expected number of descendants - div = zero(T) # ∑ (t_i - t_j) D_{ij} - for i in 2:n - w = (ψ * ω) / lambda_ts[i] # factor common to all j= max_iter && @warn "Maximum number of iterations reached without convergence." - - # 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{HawkesProcess}, h::History{H,M}; kwargs...) where {H<:Real,M} - T = promote_type(Float64, H) - return fit(HP{T}, h; kwargs...) -end - -function time_change(h::History{R,M}, hp::HawkesProcess) where {R<:Real,M} - T = float(R) - n = nb_events(h) - n == 0 && return History(T[], zero(T), T(hp.μ * duration(h)), event_marks(h)) - times = zeros(T, n) - - # In this step, `times` is the vector A in Ozaki (1979) - # A[1] = 0, A[i] = exp(-ω(tᵢ - tᵢ₋₁)) (1 + A[i - 1]) - for i in 2:n - times[i] = exp(-hp.ω * (h.times[i] - h.times[i - 1])) * (1 + times[i - 1]) - end - tmax = exp(-hp.ω * (h.tmax - h.times[end])) * (1 + times[end]) - - # Integral of the intensity corresponding to activation functions - # Λ(tₙ) - μtₙ = (α/ω) ((n-1) - A[n]) - for i in eachindex(times) - times[i] = (hp.α / hp.ω) * ((i - 1) - times[i]) # Add contribution of activation functions - end - tmax = (hp.α / hp.ω) * (n - tmax) # Add contribution of activation functions - - # Add integral of base intensity - times .+= T.(hp.μ .* (h.times .- h.tmin)) - tmax += T(hp.μ * duration(h)) - - return History(; times=times, marks=h.marks, tmin=zero(T), tmax=tmax, check=false) # A time re-scaled process starts at t=0 -end - -function ground_intensity(hp::HawkesProcess, h::History, t) - activation = sum(exp.(hp.ω .* (@view h.times[1:(searchsortedfirst(h.times, t) - 1)]))) - return hp.μ + (hp.α * activation / exp(hp.ω * t)) -end - -function integrated_ground_intensity(hp::HawkesProcess, h::History, tmin, tmax) - times = event_times(h, h.tmin, tmax) - integral = 0 - for ti in times - # Integral of activation function. 'max(tmin - ti, 0)' corrects for events that occurred - # inside or outside the interval [tmin, tmax]. - integral += (exp(-hp.ω * max(tmin - ti, 0)) - exp(-hp.ω * (tmax - ti))) - end - integral *= hp.α / hp.ω - integral += hp.μ * (tmax - tmin) # Integral of base rate - return integral -end +#= + process_mark(distribution, mark) -function DensityInterface.logdensityof(hp::HawkesProcess, h::History) - A = zeros(nb_events(h)) # Vector A in Ozaki (1979) - for i in 2:nb_events(h) - A[i] = exp(-hp.ω * (h.times[i] - h.times[i - 1])) * (1 + A[i - 1]) - end - return sum(log.(hp.μ .+ (hp.α .* A))) - # Value of intensity at each event - (hp.μ * duration(h)) - # Integral of base rate - ((hp.α / hp.ω) * sum(1 .- exp.(-hp.ω .* (duration(h) .- h.times)))) # Integral of each kernel -end +For all algorithms, an unmarked process and a process +whose marks are all equal to 1 are equivalent. +This unifies the calculation of intensities and avoids +code repetition. +=# +process_mark(::Any, ::Nothing) = 1.0 +process_mark(::Any, m::Real) = m +process_mark(::Type{Dirac{Nothing}}, ::Real) = 1.0 #= -Internal function for simulating Hawkes processes -The first generation, gen_0, is the `immigrants`, which is a set of event times. -For each t_g ∈ gen_n, simulate an inhomogeneous Poisson process over the interval [t_g, T] -with intensity λ(t) = α exp(-ω(t - t_g)) with the inverse method. -gen_{n+1} is the set of all events simulated from all events in gen_n. -The algorithm stops when the simulation from one generation results in no further events. + update_A(α, ω, t, ti_1, mi_1, Ai_1, distribution) + +Used for calculating the elements in the vector A in Ozaki (1979). There, the +definition of A is + A[0] = 0 + A[i] = exp(-ω (tᵢ - tᵢ₋₁)) * (1 + A[i - 1]) +and λ(tᵢ) = μ + α A[i]. +Here we make two modifications. +First, if α is not fixed for all events, that is, we have α₁, ..., αₙ +corresponding to each of the events t₁, ..., tₙ (in the marked case, +αᵢ = αmᵢ, mᵢ the i-th mark), then we define + A[0] = 0 + A[i] = exp(-ω (tᵢ - tᵢ₋₁)) * (αᵢ + A[i - 1]) +so λ(tᵢ) = μ + A[i]. +The second is that we allow A to be updated with an arbitrary time t, because we +can calculate λ(t) as λ(t) = μ + At, with + At = exp(-ω (t - tₖ)) * (αₖ + A[k]) +where tₖ is the last event time of the process before t. =# -function generate_descendants( - rng::AbstractRNG, immigrants::Vector{T}, tmax, α, ω -) where {T<:Real} - descendants = T[] - next_gen = immigrants - while !isempty(next_gen) - # OPTIMIZE: Can this be improved by avoiding allocations of `curr_gen` and `next_gen`? Or does the compiler take care of that? - curr_gen = copy(next_gen) # The current generation from which we simulate the next one - next_gen = eltype(immigrants)[] # Gathers all the descendants from the current generation - for parent in curr_gen # Generate the descendants for each individual event with the inverse method - activation_integral = (α / ω) * (one(T) - exp(ω * (parent - tmax))) - sim_transf = simulate_poisson_times(rng, one(T), zero(T), activation_integral) - @. sim_transf = parent - (inv(ω) * log(one(T) - ((ω / α) * sim_transf))) # Inverse of integral of the activation function - append!(next_gen, sim_transf) - end - append!(descendants, next_gen) - end - return descendants +function update_A(α::R, ω::R, t::Real, ti_1::Real, Ai_1)::float(R) where {R<:Real} + return exp(-ω * (t - ti_1)) * (α + Ai_1) end diff --git a/src/hawkes/multivariate/intensity.jl b/src/hawkes/multivariate/intensity.jl new file mode 100644 index 0000000..d0b0063 --- /dev/null +++ b/src/hawkes/multivariate/intensity.jl @@ -0,0 +1,128 @@ +#= +OPTMIZE: If instead of calculating for a single time `t` we want to return + the intensities for a ORDERED vector fo times, these methods calculating + be refactored to calculate all the intensities in a single pass. +=# + +function ground_intensity(hp::MultivariateHawkesProcess, h::History{<:Real,<:Int}, t::Real) + return mapreduce(m -> intensity(hp, m, t, h), +, support(hp.mark_dist)) +end + +#= +OPTMIZE: Calculating `ground_intensity` from `intensity` is cleaner, but + it is more efficient to calculate directly in one pass. +=# +function intensity( + hp::MultivariateHawkesProcess{R}, m::Int, t::Real, h::History{<:Real,<:Int} +) where {R<:Real} + T = float(R) + m in support(hp.mark_dist) || return zero(T) + + # If t is outside of history, intensity is 0 + (t < h.tmin || t > h.tmax) && return zero(T) + + μ = T(hp.μ * probs(hp.mark_dist)[m]) + α = T.(hp.α[:, m]) + ω = T(hp.ω[m]) + + # If t is before first event, intensity is just the base intensity + t <= h.times[1] && return μ + + # Calculate the contribution of the activation functions using Ozaki (1979) + λt = zero(T) + ind_h = 2 + while ind_h <= nb_events(h) && t > h.times[ind_h] + mi_1 = h.marks[ind_h - 1] + λt = update_A(α[mi_1], ω, h.times[ind_h], h.times[ind_h - 1], λt) + ind_h += 1 + end + mi_1 = h.marks[ind_h - 1] + λt = update_A(α[mi_1], ω, t, h.times[ind_h - 1], λt) + + # Add the base intensity + λt += μ + + return λt +end + +#= +OPTMIZE: Calculating `integrated`ground`intensity` from `integrated_intensity` is cleaner, but + it is more efficient to calculate directly in one pass. +=# +function integrated_intensity( + hp::MultivariateHawkesProcess{R}, m::Int, h::History{<:Real,<:Int}, t::Real +) where {R<:Real} + T = float(R) + m in support(hp.mark_dist) || return zero(T) + + # If t occurs before the history starts, then the integral is 0 + t < h.tmin && return zero(T) + + μ = T(hp.μ * probs(hp.mark_dist)[m]) + α = T.(hp.α[:, m]) + ω = T(hp.ω[m]) + + Λt = μ * T(t) + + # If t is before first event, intensity is just the base intensity + t <= h.times[1] && return Λt + + # Calculate the contribution of the activation functions using Ozaki (1979) + ind_h = 2 + sum_marks = zero(T) + A = zero(T) + while ind_h <= nb_events(h) && t > h.times[ind_h] + mi_1 = h.marks[ind_h - 1] + sum_marks += α[mi_1] + A = update_A(α[mi_1], ω, h.times[ind_h], h.times[ind_h - 1], A) + ind_h += 1 + end + mi_1 = h.marks[ind_h - 1] + sum_marks += α[mi_1] + A = update_A(α[mi_1], ω, t, h.times[ind_h - 1], A) + + # Add contribution of intensity functions + Λt += inv(ω) * (sum_marks - A) + + return Λt +end + +function integrated_intensity( + hp::MultivariateHawkesProcess, m::Int, h::History{<:Real,<:Int}, a, b +) + return integrated_intensity(hp, m, h, b) - integrated_intensity(hp, m, h, a) +end + +function integrated_ground_intensity(hp, h, t) + return mapreduce(m -> integrated_intensity(hp, m, h, t), +, support(hp.mark_dist)) +end + +function integrated_ground_intensity(hp, h, a, b) + return mapreduce(m -> integrated_intensity(hp, m, h, a, b), +, support(hp.mark_dist)) +end + +function DensityInterface.logdensityof( + hp::MultivariateHawkesProcess{R}, h::History{<:Real,<:Int} +) where {R<:Real} + T = float(R) + + isempty(h.times) && return T(hp.μ * duration(h)) + + μ = T.(hp.μ .* probs(hp.mark_dist)) + α = T.(hp.α) + ω = T.(hp.ω) + + A = zero(μ) + sum_marks = α[h.marks[1], :] + sum_logλ = log(μ[h.marks[1]]) + for i in 2:nb_events(h) + mi_1 = h.marks[i - 1] + mi = h.marks[i] + sum_marks .+= α[mi, :] + A .= update_A.(α[mi_1, :], ω, h.times[i], h.times[i - 1], A) + sum_logλ += log(μ[mi] + A[mi]) + end + A .= update_A.(α[h.marks[end], :], ω, h.tmax, h.times[end], A) + ΛT = (μ .* duration(h)) .+ inv.(ω) .* (sum_marks .- A) + return sum_logλ - sum(ΛT) +end diff --git a/src/hawkes/multivariate/multivariate_hawkes.jl b/src/hawkes/multivariate/multivariate_hawkes.jl new file mode 100644 index 0000000..74af781 --- /dev/null +++ b/src/hawkes/multivariate/multivariate_hawkes.jl @@ -0,0 +1,85 @@ +""" + 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. +μ and ω are vectors of length M with elements μₘ, ωₘ. α is a M×M +matrix with elements αₘₗ. + +The process is represented as a marked process, where each marginal +process m = 1, 2, ..., M is corresponds to events (tᵐᵢ, m), tᵐᵢ being +the i-th element with mark m. + +This is the 'Linear Multidimensional EHP' in Section 2 of +[Bonnet, Dion-Blanc and Perrin (2024)](https://arxiv.org/pdf/2410.05008v1) + +# Fields +- `μ::Vector{<:Real}`: baseline intensity (immigration rate) +- `α::Matrix{<:Real}`: Jump size. +- `ω::Vector{<:Real}`: Decay rate. +""" +struct MultivariateHawkesProcess{T} <: HawkesProcess{T,Categorical{Float64,Vector{Float64}}} + μ::T + α::Matrix{T} + ω::Vector{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 + +function HawkesProcess( + μ::Vector{<:Real}, α::Matrix{<:Real}, ω::Vector{<:Real}; check_args::Bool=true +) + check_args && check_args_Hawkes(μ, α, ω) + R = promote_type(eltype(μ), eltype(α), eltype(ω)) + return MultivariateHawkesProcess( + R(sum(μ)), R.(α), R.(ω), Categorical(Float64.(μ / sum(μ))) + ) +end + +# For M independent marginal processes +function HawkesProcess( + μ::Vector{<:Real}, α::Vector{<:Real}, ω::Vector{<:Real}; check_args::Bool=true +) + return HawkesProcess(μ, diagm(α), ω; check_args=check_args) +end + +function check_args_Hawkes(μ::Vector{<:Real}, α::Matrix{<:Real}, ω::Vector{<:Real}) + if (length(μ) != size(α)[2]) + throw(DimensionMismatch("α must have size $(length(μ))×$(length(μ))")) + end + if (length(μ) != length(ω)) + throw(DimensionMismatch("α and ω must have the same length")) + end + if any(α ./ ω' .>= 1) + @warn """HawkesProcess: There are parameters αᵢⱼ and ωⱼ with + αᵢⱼ / ωⱼ >= 1. This may cause problems, especially in simulations.""" + end + if any(ω .<= zero(ω)) + throw( + DomainError( + "ω = $ω", "HawkesProcess: All elements of ω must be strictly positive." + ), + ) + end + if any(μ .< zero(μ)) || any(α .< zero(α)) + throw( + DomainError( + "μ = $μ, α = $α", + "HawkesProcess: All elements of μ, α and ω must be non-negative.", + ), + ) + end + return nothing +end diff --git a/src/hawkes/multivariate/simulation.jl b/src/hawkes/multivariate/simulation.jl new file mode 100644 index 0000000..5bbe08f --- /dev/null +++ b/src/hawkes/multivariate/simulation.jl @@ -0,0 +1,10 @@ +function one_parent!( + rng::AbstractRNG, h::History{T}, mh::MultivariateHawkesProcess, parent_t, parent_m +) where {T<:Real} + for m in support(mh.mark_dist) + sim_transf = descendants(rng, parent_t, mh.α[parent_m, m], mh.ω[m], h.tmax) + append!(h.times, sim_transf) + append!(h.marks, fill(m, length(sim_transf))) + end + return nothing +end diff --git a/src/hawkes/multivariate/time_change.jl b/src/hawkes/multivariate/time_change.jl new file mode 100644 index 0000000..0485a73 --- /dev/null +++ b/src/hawkes/multivariate/time_change.jl @@ -0,0 +1,38 @@ +""" + time_change(hp::MultivariateHawkesProcess, h::History{R}) where {R<:Real} + +Time rescaling for multivariate Hawkes processes. + +Notice that the duration of the rescaled processes might be different. The output +is a single `History` on the interval from 0 to the length of the longest rescaled +process. + +Assume there are M processes, each of them with Nₘ events, m = 1, ..., M. Then the +m-th rescaled process will be defined on the interval [0, Nₘ). +""" +function time_change(hp::MultivariateHawkesProcess, h::History{R,<:Int}) where {R<:Real} + T = float(R) + μ = T.(hp.μ .* (probs(hp.mark_dist))) # Base intensity for each process + A = zero(μ) # For calculating the compensator for each process + sum_marks = zero(μ) + times_transformed = zero(h.times) # The rescaled event times + times_transformed[1] = (h.times[1] - h.tmin) * μ[h.marks[1]] # No activations before first event + for i in 2:nb_events(h) + mi_1 = h.marks[i - 1] + mi = h.marks[i] + α = hp.α[mi_1, :] + sum_marks .+= α + A .= update_A.(α, hp.ω, h.times[i], h.times[i - 1], A) + times_transformed[i] = + ((h.times[i] - h.tmin) * μ[mi]) + # Integral of base rate + inv(hp.ω[mi]) * (sum_marks[mi] - A[mi]) # Integral of activations + end + mn = h.marks[end] + α = hp.α[mn, :] + sum_marks .+= α + A .= update_A.(α, hp.ω, h.tmax, h.times[end], A) + ΛT = (μ .* duration(h)) .+ inv.(hp.ω) .* (sum_marks .- A) # Length of the intervals of rescaled processes + return History(; + times=T.(times_transformed), tmin=zero(T), tmax=T(maximum(ΛT)), marks=event_marks(h) + ) +end diff --git a/src/hawkes/simulation.jl b/src/hawkes/simulation.jl new file mode 100644 index 0000000..f62851f --- /dev/null +++ b/src/hawkes/simulation.jl @@ -0,0 +1,49 @@ +function simulate(rng::AbstractRNG, hp::HawkesProcess, tmin::Real, tmax::Real) + h = simulate(rng, PoissonProcess(hp.μ, hp.mark_dist), tmin, tmax) + generate_descendants!(rng, h, hp) # Recursively generates descendants from first events + perm = sortperm(h.times) + h.times .= h.times[perm] + h.marks .= h.marks[perm] + return h +end + +function simulate(hp::HawkesProcess, tmin::Real, tmax::Real) + return simulate(default_rng(), hp, tmin, tmax) +end + +#= +Internal function for simulating Hawkes processes +The history `h` contains the event times generated from the base rate of the process, the +"parents" of the process. This is the starting generation, gen_0. +For each t_g ∈ gen_n, simulate an inhomogeneous Poisson process over the interval [t_g, T] +with intensity λ(t) = α exp(-ω(t - t_g)) using the inverse method. +gen_{n+1} is the set of all events simulated from all events in gen_n. +The algorithm stops when the simulation from one generation results in no further events. +=# +function generate_descendants!( + rng::AbstractRNG, h::History{T,M}, hp::HawkesProcess +) where {T<:Real,M} + gen_start = 1 # Starting index of generation gen_0 + gen_end = nb_events(h) # Last index of generation gen_0 + while gen_start <= gen_end # As long as the current genertion is not empty, generate more descendants + for parent_idx in gen_start:gen_end + # generates the descendants of one single element of the current generation + one_parent!(rng, h, hp, h.times[parent_idx], h.marks[parent_idx]) # dispatch on the type of `hp` + end + # gen_i+1 are the descendants of gen_i + gen_start = gen_end + 1 + gen_end = nb_events(h) + end + return nothing +end + +# generates the descendants of one single parent using the inverse transform method +function descendants(rng::AbstractRNG, parent::R, α::Real, ω::Real, tmax::R) where {R<:Real} + T = float(R) + αT = T(α) + ωT = T(ω) + activation_integral = (αT / ωT) * (one(T) - exp(ωT * (parent - tmax))) + sim_transf = simulate_poisson_times(rng, one(T), zero(T), activation_integral) + @. sim_transf = parent - (inv(ωT) * log(one(T) - ((ωT / αT) * sim_transf))) # Inverse of integral of the activation function + return sim_transf +end diff --git a/src/hawkes/univariate/fit.jl b/src/hawkes/univariate/fit.jl new file mode 100644 index 0000000..71fa00b --- /dev/null +++ b/src/hawkes/univariate/fit.jl @@ -0,0 +1,123 @@ +""" + StatsAPI.fit(::Type{<:UnivariateHawkesProcess{R}}, h::History; step_tol::Float64=1e-6, max_iter::Int=1000, rng::AbstractRNG=default_rng()) where {R<:Real} + +Adaptation of the Expectation-Maximization algorithm from [E. Lewis, G. Mohler (2011)](https://arxiv.org/pdf/1801.08273)) +for marked processes. The relevant calculations are in page 4, equations 6-13. + +Let (t₁, m₁), (t₂, m₂), ..., (tₙ, mₙ) be a marked event history with (t₁ < ... < tₙ) 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 α mᵢ 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ᵢⱼ = ψ ω mⱼ 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 / ∑ mᵢ + ω = D / div +5. If convergence criterion is met, return updated parameters, otherwise, back to step 2. + +Notice that, for numerical stability, 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( + HP::Type{UnivariateHawkesProcess{R,D}}, + h::History; + step_tol::Float64=1e-6, + max_iter::Int=1000, + rng::AbstractRNG=default_rng(), +) where {R,D} + div_ψ = D == Dirac{Nothing} ? nb_events(h) : sum(h.marks) + params = fit_hawkes_params(HP, h, div_ψ; step_tol=step_tol, max_iter=max_iter, rng=rng) + d = D == Dirac{Nothing} ? Dirac(nothing) : fit(D, h.marks) + return HawkesProcess(params..., d) +end + +#= +Method for fitting the parameters of marked and unmarked Hawkes +processes. +`div_ψ` must be passed as an argument. For unmarked processes, +`div_ψ` is equal to the number of events in the evet history, for +marked processes, `div_ψ` is equal to the sum of all the marks. +=# +function fit_hawkes_params( + ::Type{UnivariateHawkesProcess{R,Dist}}, + h::History, + div_ψ::Real; + step_tol::Float64=1e-6, + max_iter::Int=1000, + rng::AbstractRNG=default_rng(), +) where {R<:Real,Dist} + T = float(R) + div_ψ = T(div_ψ) + n = nb_events(h) + N = 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= step_tol) && (n_iters < max_iter) + # compute A, S, and λ + for i in 2:n + Δ = norm_ts[i] - norm_ts[i - 1] + e = exp(-ω * Δ) + mi_1 = process_mark(Dist, h.marks[i - 1]) + A[i] = e * (mi_1 + A[i - 1]) # Different updates for real marks and no marks + S[i] = (Δ * A[i]) + (e * S[i - 1]) + lambda_ts[i] = μ + (ψ * ω) * A[i] + end + + # E-step aggregates + D = zero(T) # expected number of descendants + div = zero(T) # ∑ (t_i - t_j) D_{ij} + for i in 2:n + w = (ψ * ω) / lambda_ts[i] # factor common to all j= max_iter && @warn "Maximum number of iterations reached without convergence." + + # 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 (μ * (N / tmax), ψ * ω * (N / tmax), ω * (N / tmax)) +end diff --git a/src/hawkes/univariate/intensity.jl b/src/hawkes/univariate/intensity.jl new file mode 100644 index 0000000..73ec66d --- /dev/null +++ b/src/hawkes/univariate/intensity.jl @@ -0,0 +1,99 @@ +function ground_intensity( + hp::UnivariateHawkesProcess{R,D}, h::History, t::Real +) where {R<:Real,D} + T = float(R) + μ, α, ω = T.((hp.μ, hp.α, hp.ω)) + + # If t is outside of history, intensity is 0 + (t < h.tmin || t > h.tmax) && return zero(T) + + # If t is before first event, intensity is just the base intensity + t <= h.times[1] && return μ + + # Calculate the contribution of the activation functions using Ozaki (1979) + λt = zero(T) + ind_h = 2 + while ind_h <= nb_events(h) && t > h.times[ind_h] + mi_1 = T(process_mark(D, h.marks[ind_h - 1])) + λt = update_A(mi_1 * α, ω, h.times[ind_h], h.times[ind_h - 1], λt) + ind_h += 1 + end + mi_1 = T(process_mark(D, h.marks[ind_h - 1])) + λt = update_A(mi_1 * α, ω, t, h.times[ind_h - 1], λt) + + # Add the base intensity + λt += hp.μ + + return λt +end + +function intensity(hp::UnivariateHawkesProcess{R}, m, t::Real, h::History) where {R<:Real} + T = float(R) + m in support(hp.mark_dist) || return zero(T) + return ground_intensity(hp, h, t) * T(pdf(hp.mark_dist, m)) # `pdf` always returns a `Float64` +end + +function intensity(hp::UnmarkedUnivariateHawkesProcess, _, t::Real, h::History) + return ground_intensity(hp, h, t) +end + +function integrated_ground_intensity( + hp::UnivariateHawkesProcess{R,D}, h::History, t +) where {R<:Real,D} + T = float(R) + μ, α, ω = T.((hp.μ, hp.α, hp.ω)) + + # If t is outside of history, intensity is 0 + t <= h.tmin && return zero(T) + + # Add the base intensity + Λt = μ * T(t) + + # If t is before first event, intensity is just the base intensity + t <= h.times[1] && return Λt + + # Calculate the contribution of the activation functions using Ozaki (1979) + A = zero(T) + sum_marks = zero(T) + ind_h = 2 + while ind_h <= nb_events(h) && t > h.times[ind_h] + mi_1 = T(process_mark(D, h.marks[ind_h - 1])) + sum_marks += mi_1 + A = update_A(mi_1 * α, ω, h.times[ind_h], h.times[ind_h - 1], A) + ind_h += 1 + end + mi_1 = T(process_mark(D, h.marks[ind_h - 1])) + sum_marks += mi_1 + A = update_A(mi_1 * α, ω, t, h.times[ind_h - 1], A) + + # Add intgral of activation functions + Λt += inv(ω) * (sum_marks - A) + + return Λt +end + +function integrated_ground_intensity(hp::UnivariateHawkesProcess, h::History, a, b) + return integrated_ground_intensity(hp, h, b) - integrated_ground_intensity(hp, h, a) +end + +function DensityInterface.logdensityof( + hp::UnivariateHawkesProcess{R,D}, h::History +) where {R<:Real,D} + T = float(R) + μ, α, ω = T.((hp.μ, hp.α, hp.ω)) + + A = zero(T) + sum_marks = zero(T) + sum_logλ = zero(T) + for i in 2:nb_events(h) + mi_1 = T(process_mark(D, h.marks[i - 1])) + sum_marks += mi_1 + A = update_A(mi_1 * α, ω, h.times[i], h.times[i - 1], A) + sum_logλ += log(μ + α * A) + end + mn = T(process_mark(D, h.marks[end])) + sum_marks += mn + A = update_A(mn * α, ω, h.tmax, h.times[end], A) + Λ_T = (μ * duration(h)) + inv.(ω) * (sum_marks - A) + return sum_logλ - sum(Λ_T) +end diff --git a/src/hawkes/univariate/simulation.jl b/src/hawkes/univariate/simulation.jl new file mode 100644 index 0000000..da5c98d --- /dev/null +++ b/src/hawkes/univariate/simulation.jl @@ -0,0 +1,17 @@ +function one_parent!( + rng::AbstractRNG, h::History, hp::UnivariateHawkesProcess, parent_t, parent_m +) + sim_transf = descendants(rng, parent_t, (parent_m * hp.α), hp.ω, h.tmax) + append!(h.times, sim_transf) + append!(h.marks, [rand(rng, mark_distribution(hp)) for _ in 1:length(sim_transf)]) + return nothing +end + +function one_parent!( + rng::AbstractRNG, h::History, hp::UnmarkedUnivariateHawkesProcess, parent_t, _ +) + sim_transf = descendants(rng, parent_t, hp.α, hp.ω, h.tmax) + append!(h.times, sim_transf) + append!(h.marks, fill(nothing, length(sim_transf))) + return nothing +end diff --git a/src/hawkes/univariate/time_change.jl b/src/hawkes/univariate/time_change.jl new file mode 100644 index 0000000..33f48d6 --- /dev/null +++ b/src/hawkes/univariate/time_change.jl @@ -0,0 +1,38 @@ +function time_change(hp::UnivariateHawkesProcess{<:Real,D}, h::History{R}) where {R<:Real,D} + T = float(R) + n = nb_events(h) + A = zero(T) + μ, α, ω = T.((hp.μ, hp.α, hp.ω)) + sum_marks = zero(T) + times_transformed = zeros(T, n) + + # Loop to calculate i -> ∑_{j= ω + @warn """HawkesProcess: α / ω >= 1. This may cause problems, + especially in simulations.""" + end +end + +function check_args_Hawkes(μ::Real, α::Real, ω::Real, mark_dist) + mean_α = mark_dist isa Dirac{Nothing} ? α : mean(mark_dist) * α + check_args_Hawkes(μ, mean_α, ω) + 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 diff --git a/src/history.jl b/src/history.jl index b0463e4..7121d7f 100644 --- a/src/history.jl +++ b/src/history.jl @@ -16,14 +16,8 @@ struct History{T<:Real,M} tmax::T marks::Vector{M} - function History( - times::Vector{<:Real}, - tmin::Real, - tmax::Real, - marks=fill(nothing, length(times)); - check=true, - ) - if check + function History(times, tmin, tmax, marks=fill(nothing, length(times)); check_args=true) + if check_args if tmin >= tmax throw( DomainError( @@ -55,8 +49,8 @@ struct History{T<:Real,M} end end -function History(; times, tmin, tmax, marks=fill(nothing, length(times)), check=true) - History(times, tmin, tmax, marks; check) +function History(; times, tmin, tmax, marks=fill(nothing, length(times)), check_args=true) + History(times, tmin, tmax, marks; check_args) end function Base.show(io::IO, h::History{T,M}) where {T,M} @@ -65,6 +59,19 @@ function Base.show(io::IO, h::History{T,M}) where {T,M} ) end +""" + length(h) + +Alias for `nb_events(h)`. +""" +Base.length(h::History) = nb_events(h) + +Base.firstindex(h::History) = 1 +Base.lastindex(h::History) = length(h) +Base.getindex(h::History, i::Int) = (h.times[i], h.marks[i]) +Base.getindex(h::History, I::Vararg{Int}) = [(h.times[i], h.marks[i]) for i in I] +Base.getindex(h::History, I::AbstractVector{<:Int}) = getindex(h, I...) + """ event_times(h) @@ -118,13 +125,6 @@ Count events in `h`. """ nb_events(h::History) = length(h.marks) -""" - length(h) - -Alias for `nb_events(h)`. -""" -Base.length(h::History) = nb_events(h) - """ nb_events(h, tmin, tmax) @@ -162,8 +162,8 @@ duration(h::History) = max_time(h) - min_time(h) Add event `(t, m)` inside the interval `[h.tmin, h.tmax)` at the end of history `h`. """ -function Base.push!(h::History, t::Real, m; check=true) - if check +function Base.push!(h::History, t::Real, m; check_args=true) + if check_args @assert h.tmin <= t < h.tmax @assert (length(h) == 0) || (h.times[end] <= t) end @@ -177,8 +177,8 @@ end Append events `(ts, ms)` inside the interval `[h.tmin, h.tmax)` at the end of history `h`. """ -function Base.append!(h::History, ts::Vector{<:Real}, ms; check=true) - if check +function Base.append!(h::History, ts::Vector{<:Real}, ms; check_args=true) + if check_args perm = sortperm(ts) ts .= ts[perm] ms .= ms[perm] @@ -206,7 +206,7 @@ function Base.cat(h1::History, h2::History) ) times = [h1.times; h2.times] marks = [h1.marks; h2.marks] - return History(; times=times, tmin=h1.tmin, tmax=h2.tmax, marks=marks, check=false) + return History(; times=times, tmin=h1.tmin, tmax=h2.tmax, marks=marks, check_args=false) end """ @@ -219,7 +219,9 @@ function time_change(h::History, Λ) new_marks = copy(event_marks(h)) new_tmin = Λ(min_time(h)) new_tmax = Λ(max_time(h)) - return History(; times=new_times, marks=new_marks, tmin=new_tmin, tmax=new_tmax) + return History(; + times=new_times, marks=new_marks, tmin=new_tmin, tmax=new_tmax, check_args=false + ) end """ @@ -236,7 +238,7 @@ function split_into_chunks(h::History{T,M}, chunk_duration) where {T,M} for (a, b) in zip(limits[1:(end - 1)], limits[2:end]) times = [t for t in event_times(h) if a <= t < b] marks = [m for (t, m) in zip(event_times(h), event_marks(h)) if a <= t < b] - chunk = History(; times=times, marks=marks, tmin=a, tmax=b, check=false) + chunk = History(; times=times, marks=marks, tmin=a, tmax=b, check_args=false) push!(chunks, chunk) end return chunks diff --git a/src/poisson/fit.jl b/src/poisson/fit.jl index 4329c11..b7b7863 100644 --- a/src/poisson/fit.jl +++ b/src/poisson/fit.jl @@ -19,7 +19,7 @@ function StatsAPI.fit( return PoissonProcess(λ, mark_dist) end -function StatsAPI.fit(pptype::Type{PoissonProcess{R,D}}, args...; kwargs...) where {R,D} +function StatsAPI.fit(pptype::Type{<:PoissonProcess}, args...; kwargs...) ss = suffstats(pptype, args...) return fit(pptype, ss) end diff --git a/src/poisson/poisson_process.jl b/src/poisson/poisson_process.jl index 54cc8b8..6ead0ec 100644 --- a/src/poisson/poisson_process.jl +++ b/src/poisson/poisson_process.jl @@ -64,6 +64,7 @@ function Base.show(io::IO, pp::MultivariatePoissonProcess) end ## Constructors +### MarkedPoissonProcess function PoissonProcess(λ::Vector{R}; check_args::Bool=true) where {R<:Real} if check_args if any(λ .< zero(λ)) @@ -83,6 +84,7 @@ function PoissonProcess(λ::Vector{R}; check_args::Bool=true) where {R<:Real} return PoissonProcess(sum(λ), Categorical(λ / sum(λ)); check_args=check_args) end +### UnivariatePoissonProcess function PoissonProcess(λ::R; check_args::Bool=true) where {R<:Real} return PoissonProcess(λ, Dirac(nothing); check_args=check_args) end diff --git a/src/poisson/simulation.jl b/src/poisson/simulation.jl index aa22d8d..e4a66c1 100644 --- a/src/poisson/simulation.jl +++ b/src/poisson/simulation.jl @@ -1,5 +1,5 @@ function simulate(rng::AbstractRNG, pp::PoissonProcess, tmin::T, tmax::T) where {T<:Real} times = simulate_poisson_times(rng, ground_intensity(pp), tmin, tmax) marks = [rand(rng, mark_distribution(pp)) for _ in 1:length(times)] - return History(; times=times, marks=marks, tmin=tmin, tmax=tmax) + return History(; times=times, marks=marks, tmin=tmin, tmax=tmax, check_args=false) end diff --git a/src/simulation.jl b/src/simulation.jl index 2c39803..9844166 100644 --- a/src/simulation.jl +++ b/src/simulation.jl @@ -5,6 +5,7 @@ Simulate the event times of a homogeneous Poisson process with parameter λ on t Internal function to use in all other simulation algorithms. =# function simulate_poisson_times(rng::AbstractRNG, λ, tmin, tmax) + λ == 0 && return Float64[] N = rand(rng, Poisson(λ * (tmax - tmin))) times = [rand(rng, Uniform(tmin, tmax)) for _ in 1:N] # rand(rng, Uniform(tmin, tmax), N) always outputs a `Float64` sort!(times) @@ -19,9 +20,8 @@ Simulate a temporal point process `pp` on interval `[tmin, tmax)` using Ogata's # Technical Remark To infer the type of the marks, the implementation assumes that there is method of `mark_distribution` without the argument `h` such that it corresponds to the distribution of marks in case the history is empty. """ -function simulate_ogata( - rng::AbstractRNG, pp::AbstractPointProcess, tmin::T, tmax::T -) where {T<:Real} +function simulate_ogata(rng::AbstractRNG, pp::AbstractPointProcess, tmin::Real, tmax::Real) + T = promote_type(typeof(tmin), typeof(tmax)) M = typeof(rand(mark_distribution(pp, tmin))) h = History(; times=T[], marks=M[], tmin=tmin, tmax=tmax) t = tmin @@ -31,12 +31,12 @@ function simulate_ogata( if τ > L t = t + L elseif τ <= L - U_max = ground_intensity(pp, t + τ, h) / B + U_max = ground_intensity(pp, h, t + τ) / B U = rand(rng, typeof(U_max)) if U < U_max m = rand(rng, mark_distribution(pp, t + τ, h)) if t + τ < tmax - push!(h, t + τ, m; check=false) + push!(h, t + τ, m; check_args=false) end end t = t + τ diff --git a/test/bounded_point_process.jl b/test/bounded_point_process.jl index 2996869..4286673 100644 --- a/test/bounded_point_process.jl +++ b/test/bounded_point_process.jl @@ -6,11 +6,11 @@ h = simulate(rng, bpp) @test min_time(bpp) == 0.0 @test max_time(bpp) == 1000.0 -@test ground_intensity(bpp, 0, h) == sum(intensities) -@test mark_distribution(bpp, 100.0, h) == Categorical(intensities / sum(intensities)) +@test ground_intensity(bpp, h, 0) == sum(intensities) +@test mark_distribution(bpp, h, 100.0) == Categorical(intensities / sum(intensities)) @test mark_distribution(bpp, 0.0) == Categorical(intensities / sum(intensities)) -@test intensity(bpp, 1, 0, h) == intensities[1] -@test log_intensity(bpp, 2, 1.0, h) == log(intensities[2]) +@test intensity(bpp, 1, h, 0) == intensities[1] +@test log_intensity(bpp, 2, h, 1.0) == log(intensities[2]) @test ground_intensity_bound(bpp, 243, h) == (sum(intensities), typemax(Int)) @test ground_intensity_bound(bpp, 243.0, h) == (sum(intensities), Inf) diff --git a/test/hawkes.jl b/test/hawkes.jl deleted file mode 100644 index abc6fda..0000000 --- a/test/hawkes.jl +++ /dev/null @@ -1,57 +0,0 @@ -# Constructor -@test HawkesProcess(1, 1, 2) isa HawkesProcess{Int} -@test HawkesProcess(1, 1, 2.0) isa HawkesProcess{Float64} -@test_throws DomainError HawkesProcess(1, 1, 1) -@test_throws DomainError HawkesProcess(-1, 1, 2) - -hp = HawkesProcess(1, 1, 2) -h = History(; times=[1.0, 2.0, 4.0], marks=["a", "b", "c"], tmin=0.0, tmax=5.0) -h_big = History(; - times=BigFloat.([1, 2, 4]), marks=["a", "b", "c"], tmin=BigFloat(0), tmax=BigFloat(5) -) - -# Time change -h_transf = time_change(h, hp) -integral = - (hp.μ * duration(h)) + - (hp.α / hp.ω) * sum([1 - exp(-hp.ω * (h.tmax - ti)) for ti in h.times]) -@test isa(h_transf, typeof(h)) -@test h_transf.marks == h.marks -@test h_transf.tmin ≈ 0 -@test h_transf.tmax ≈ integral -@test h_transf.times ≈ [1, (2 + (1 - exp(-2)) / 2), 4 + (2 - exp(-2 * 3) - exp(-2 * 2)) / 2] -@test isa(time_change(h_big, hp), typeof(h_big)) - -# Ground intensity -@test ground_intensity(hp, h, 1) == 1 -@test ground_intensity(hp, h, 2) == 1 + hp.α * exp(-hp.ω * 1) - -# Integrated ground intensity -@test integrated_ground_intensity(hp, h, h.tmin, h.tmax) ≈ integral -@test integrated_ground_intensity(hp, h, 2, 3) ≈ - hp.μ + - ((hp.α / hp.ω) * (exp(-hp.ω) - exp(-hp.ω * 2))) + - ((hp.α / hp.ω) * (1 - exp(-hp.ω))) - -# Rand -h_sim = simulate(hp, 0.0, 10.0) -@test issorted(h_sim.times) -@test isa(h_sim, History{Float64,Nothing}) -@test isa(simulate(hp, BigFloat(0), BigFloat(10)), History{BigFloat,Nothing}) - -# Fit -Random.seed!(123) -params_true = (100.0, 100.0, 200.0) -model = HawkesProcess(params_true...) -h_sim = simulate(model, 0.0, 50.0) -model_est = fit(HawkesProcess, h_sim) -params_est = (model_est.μ, model_est.α, model_est.ω) -@test isa(model_est, HawkesProcess) -@test all((params_true .* 0.9) .<= params_est .<= (params_true .* 1.1)) -@test isa(fit(HawkesProcess, h_big), HawkesProcess{BigFloat}) -@test isa(fit(HawkesProcess{Float32}, h_big), HawkesProcess{Float32}) - -# logdensityof -@test logdensityof(hp, h) ≈ - sum(log.(hp.μ .+ (hp.α .* [0, exp(-hp.ω), exp(-hp.ω * 2) + exp(-hp.ω * 3)]))) - - integral diff --git a/test/hawkes_multivariate.jl b/test/hawkes_multivariate.jl new file mode 100644 index 0000000..b7505ef --- /dev/null +++ b/test/hawkes_multivariate.jl @@ -0,0 +1,69 @@ +# Constructor +@test HawkesProcess(rand(Float32, 3), rand(Float32, 3), rand(Float32, 3) .+ 1) isa + MultivariateHawkesProcess{Float32} +@test_throws DomainError HawkesProcess(rand(3) .- 1, rand(3, 3), rand(3) .+ 1) +@test_throws DimensionMismatch HawkesProcess(rand(3), rand(2, 2), rand(3) .+ 1) +@test_warn r"may cause problems" HawkesProcess(rand(3), rand(3, 3) .+ 1, rand(3)) + +hp = HawkesProcess([1, 2, 3], [1 0.1 0.1; 0.2 2 0.2; 0.3 0.3 3], [2, 4, 6]) +h = History(; times=[1.0, 2.0, 4.0], marks=[3, 2, 1], tmin=0.0, tmax=5.0) +h_big = History(; + times=BigFloat.([1, 2, 4]), marks=[3, 2, 1], tmin=BigFloat(0), tmax=BigFloat(5) +) + +# Ground intensity +@test ground_intensity(hp, h, 1) ≈ 6 +@test ground_intensity(hp, h, 2) ≈ 6 + sum(hp.α[3, :] .* exp.(-hp.ω .* 1)) +# @test ground_intensity(hp, h, [1, 2, 3, 4]) ≈ +# [ground_intensity(hp, h, t) for t in [1, 2, 3, 4]] + +# intensity +@test intensity(hp, 1, 1, h) ≈ 1 +@test intensity(hp, 1, 2, h) ≈ 1 + hp.α[3, 1] * exp(-hp.ω[1] * 1) +@test intensity(hp, 4, 2, h) == 0 +# @test intensity(hp, 1, [1, 2, 3, 4], h) ≈ [intensity(hp, 1, t, h) for t in [1, 2, 3, 4]] + +h1 = History([1], 0, 3, [2]) +h2 = History([1, 2], 0, 3, [2, 1]) + +# Integrated ground intensity +@test integrated_ground_intensity(hp, h1, 0, 1) ≈ 6 +@test integrated_ground_intensity(hp, h1, 0, 1000) ≈ 6000 + sum(hp.α[2, :] ./ hp.ω) +@test integrated_ground_intensity(hp, h2, 0, 1000) ≈ + 6000 + sum(hp.α[2, :] ./ hp.ω) + sum(hp.α[1, :] ./ hp.ω) +# @test integrated_ground_intensity(hp, h, [1, 2, 3, 4]) ≈ +# [integrated_ground_intensity(hp, h, t) for t in [1, 2, 3, 4]] + +# Log likelihood (logintensityof) +@test logdensityof(hp, h1) ≈ + log(2) - (hp.μ * duration(h1)) - sum((hp.α[2, :] ./ hp.ω) .* (1 .- exp.(-2 .* hp.ω))) +@test logdensityof(hp, h2) ≈ + log(2) + log(1 + hp.α[2, 1] * exp(-hp.ω[1])) - (hp.μ * duration(h2)) - + sum((hp.α[2, :] ./ hp.ω) .* (1 .- exp.(-2 .* hp.ω))) - + sum((hp.α[1, :] ./ hp.ω) .* (1 .- exp.(-hp.ω))) + +# time change +h_transf = time_change(hp, h1) +integral = maximum( + ((hp.μ .* probs(hp.mark_dist)) .* duration(h1)) .+ + (hp.α[2, :] ./ hp.ω) .* (1 .- exp.(-2 .* hp.ω)), +) + +@test isa(h_transf, History{Float64,Int}) +@test h_transf.marks == h1.marks +@test h_transf.tmin ≈ 0 +@test h_transf.tmax ≈ integral +@test h_transf.times ≈ [2.0] +@test isa(time_change(hp, h_big), typeof(h_big)) + +# simulate +h_sim = simulate(hp, 0.0, 10.0) +hp_mult = HawkesProcess([1, 1], [1, 1], [2, 2]) +hp_univ = HawkesProcess(1, 1, 2) +n_sims = 10000 +mean_mult = sum([nb_events(simulate(hp_mult, 0, 10)) for _ in 1:n_sims]) / (2 * n_sims) +mean_univ = sum([nb_events(simulate(hp_univ, 0, 10)) for _ in 1:n_sims]) / n_sims + +@test issorted(h_sim.times) +@test isa(h_sim, History{Float64,Int64}) +@test isapprox(mean_mult, mean_univ, atol=1) diff --git a/test/hawkes_univariate.jl b/test/hawkes_univariate.jl new file mode 100644 index 0000000..a9afbf2 --- /dev/null +++ b/test/hawkes_univariate.jl @@ -0,0 +1,81 @@ +# Constructor +@test HawkesProcess(1, 1, 2, Uniform()) isa UnivariateHawkesProcess{Int64,Uniform{Float64}} +@test_throws DomainError HawkesProcess(1, 1, 2, Uniform(-1, 1)) +@test_warn r"may cause problems" HawkesProcess(1, 1, 2, Uniform(10, 11)) + +hp = HawkesProcess(1, 1, 2, Uniform(0, 3)) +h = History(; times=[1.0, 2.0, 4.0], marks=[3, 2, 1], tmin=0.0, tmax=5.0) +h_big = History(; + times=BigFloat.([1, 2, 4]), marks=[3.0, 2.0, 1.0], tmin=BigFloat(0), tmax=BigFloat(5) +) + +# Ground intensity +@test ground_intensity(hp, h, 1) ≈ 1 +@test ground_intensity(hp, h, 2) ≈ 1 + 3 * hp.α * exp(-hp.ω * 1) +# @test ground_intensity(hp, h, [1, 2, 3, 4]) ≈ +# [ground_intensity(hp, h, t) for t in [1, 2, 3, 4]] + +# Intensity +@test intensity(hp, 1, 1, h) ≈ ground_intensity(hp, h, 1) * (1/3) +@test intensity(hp, 1, 2, h) ≈ ground_intensity(hp, h, 2) * (1/3) +@test intensity(hp, 4, 2, h) == 0 +# @test intensity(hp, 1, [1, 2, 3, 4], h) ≈ [intensity(hp, 1, t, h) for t in [1, 2, 3, 4]] + +h1 = History([1], 0, 3, [2]) +h2 = History([1, 2], 0, 3, [2, 1]) + +# Integrated ground intensity +@test integrated_ground_intensity(hp, h1, 0, 1) ≈ 1 +@test integrated_ground_intensity(hp, h1, 0, 1000) ≈ 1001 +@test integrated_ground_intensity(hp, h2, 0, 1000) ≈ 1001.5 +# @test integrated_ground_intensity(hp, h, [1, 2, 3, 4]) ≈ +# [integrated_ground_intensity(hp, h, t) for t in [1, 2, 3, 4]] + +# Log likelihood (logintensityof) +@test logdensityof(hp, h1) ≈ log(1) - (hp.μ * duration(h1)) - (1 - exp(-4)) +@test logdensityof(hp, h2) ≈ + log(1) + log(1 + 2 * exp(-2)) - (hp.μ * duration(h2)) - ((1 - exp(-4))) - + ((1 - exp(-2)) / 2) + +# time change +h_transf = time_change(hp, h) +integral = + (hp.μ * duration(h)) + + (hp.α / hp.ω) * + sum([h.marks[i] * (1 - exp(-hp.ω * (h.tmax - h.times[i]))) for i in 1:length(h.times)]) + +@test isa(h_transf, typeof(h)) +@test h_transf.marks == h.marks +@test h_transf.tmin ≈ 0 +@test h_transf.tmax ≈ integral +@test h_transf.times ≈ + [1, 2 + (3 / 2) * (1 - exp(-2)), 4 + (3 / 2) * (1 - exp(-6)) + (1 - exp(-4))] +@test isa(time_change(hp, h_big), typeof(h_big)) + +# Fit +Random.seed!(1) +model_true = HawkesProcess(params_true..., Uniform()) +h_sim = simulate(model_true, 0.0, 50.0) +model_est = fit(UnivariateHawkesProcess{Float32,Uniform}, h_sim) +params_est = (model_est.μ, model_est.α, model_est.ω) +a_est = model_est.mark_dist.a +b_est = model_est.mark_dist.b + +@test isa(model_est, UnivariateHawkesProcess) +@test params_true[1] * 0.9 <= params_est[1] <= params_true[1] * 1.1 +@test params_true[2] * 0.9 <= params_est[2] <= params_true[2] * 1.1 +@test params_true[3] * 0.9 <= params_est[3] <= params_true[3] * 1.1 +@test a_est < 0.01 +@test b_est > 0.99 +@test typeof(model_est.μ) == Float32 + +# simulate +hp_univ = HawkesProcess(1, 1, 2, Uniform(0, 2)) +hp_unmk = HawkesProcess(1, 1, 2) +n_sims = 10000 +mean_univ = sum([nb_events(simulate(hp_univ, 0, 10)) for _ in 1:n_sims]) / n_sims +mean_unmk = sum([nb_events(simulate(hp_unmk, 0, 10)) for _ in 1:n_sims]) / n_sims +h_sim = simulate(hp, 0.0, 10.0) + +@test issorted(h_sim.times) +@test isa(h_sim, History{Float64,Float64}) diff --git a/test/hawkes_unmarked.jl b/test/hawkes_unmarked.jl new file mode 100644 index 0000000..3c23f4b --- /dev/null +++ b/test/hawkes_unmarked.jl @@ -0,0 +1,71 @@ +# Constructor +@test HawkesProcess(1, 1, 2) isa UnmarkedUnivariateHawkesProcess{Int64} +@test HawkesProcess(1, 1, 2.0) isa UnmarkedUnivariateHawkesProcess{Float64} +@test_throws DomainError HawkesProcess(-1, 1, 2) +@test_warn r"may cause problems" HawkesProcess(1, 1, 1) + +# Ground intensity +hp = HawkesProcess(1, 1, 2) +h = History(; times=[1.0, 2.0, 4.0], marks=[3, 2, 1], tmin=0.0, tmax=5.0) +h_big = History(; + times=BigFloat.([1, 2, 4]), marks=[3.0, 2.0, 1.0], tmin=BigFloat(0), tmax=BigFloat(5) +) + +@test ground_intensity(hp, h, 1) ≈ 1 +@test ground_intensity(hp, h, 2) ≈ 1 + hp.α * exp(-hp.ω * 1) +# @test ground_intensity(hp, h, [1, 2, 3, 4]) ≈ +# [ground_intensity(hp, h, t) for t in [1, 2, 3, 4]] + +# intensity +@test intensity(hp, "a", 1, h) ≈ ground_intensity(hp, h, 1) +@test intensity(hp, 1, 2, h) ≈ ground_intensity(hp, h, 2) +# @test intensity(hp, 1, [1, 2, 3, 4], h) ≈ [intensity(hp, 1, t, h) for t in [1, 2, 3, 4]] + +# Integrated ground intensity +h1 = History([1], 0, 3, [2]) +h2 = History([1, 2], 0, 3, [2, 1]) + +@test integrated_ground_intensity(hp, h1, 0, 1) ≈ 1 +@test integrated_ground_intensity(hp, h1, 0, 1000) ≈ 1000.5 +@test integrated_ground_intensity(hp, h2, 0, 1000) ≈ 1001 +# @test integrated_ground_intensity(hp, h, [1, 2, 3, 4]) ≈ +# [integrated_ground_intensity(hp, h, t) for t in [1, 2, 3, 4]] + +# Log likelihood (logintensityof) +@test logdensityof(hp, h1) ≈ log(1) - (hp.μ * duration(h1)) - (1 - exp(-4)) / 2 +@test logdensityof(hp, h2) ≈ + log(1) + log(1 + exp(-2)) - (hp.μ * duration(h2)) - ((1 - exp(-4)) / 2) - + ((1 - exp(-2)) / 2) + +# time change +h_transf = time_change(hp, h) +integral = + (hp.μ * duration(h)) + + (hp.α / hp.ω) * sum([1 - exp(-hp.ω * (h.tmax - ti)) for ti in h.times]) + +@test isa(h_transf, typeof(h)) +@test h_transf.marks == h.marks +@test h_transf.tmin ≈ 0 +@test h_transf.tmax ≈ integral +@test h_transf.times ≈ [1, (2 + (1 - exp(-2)) / 2), 4 + (2 - exp(-2 * 3) - exp(-2 * 2)) / 2] +@test isa(time_change(hp, h_big), typeof(h_big)) + +# Fit +# Random.seed!(1) +params_true = (100.0, 100.0, 200.0) +model_true = HawkesProcess(params_true...) +h_sim = simulate(model_true, 0.0, 50.0) +model_est = fit(UnmarkedUnivariateHawkesProcess{Float64}, h_sim) +params_est = (model_est.μ, model_est.α, model_est.ω) + +@test isa(model_est, UnmarkedUnivariateHawkesProcess) +@test all((params_true .* 0.9) .<= params_est .<= (params_true .* 1.1)) +@test isa( + fit(UnmarkedUnivariateHawkesProcess{Float32}, h_sim), + UnmarkedUnivariateHawkesProcess{Float32}, +) + +# simulate +h_sim = simulate(hp, 0.0, 10.0) +@test issorted(h_sim.times) +@test isa(h_sim, History{Float64,Nothing}) diff --git a/test/history.jl b/test/history.jl index f59e269..7052fba 100644 --- a/test/history.jl +++ b/test/history.jl @@ -13,6 +13,11 @@ h = History([0.2, 0.8, 1.1], 0.0, 2.0, ["a", "b", "c"]); @test event_marks(h) == ["a", "b", "c"] @test event_marks(h, 0.2, 0.8) == ["a"] @test event_marks(h, 0.8, 0.2) == [] +@test h[begin] == (0.2, "a") +@test h[end] == (1.1, "c") +@test h[2] == (0.8, "b") +@test h[1, 2] == [(0.2, "a"), (0.8, "b")] +@test h[1:3] == [(0.2, "a"), (0.8, "b"), (1.1, "c")] push!(h, 1.7, "d") diff --git a/test/runtests.jl b/test/runtests.jl index 382405c..6144ed5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -40,7 +40,15 @@ DocMeta.setdocmeta!(PointProcesses, :DocTestSetup, :(using PointProcesses); recu end end @testset verbose = true "Hawkes" begin - include("hawkes.jl") + @testset verbose = true "Unmarked Univariate" begin + include("hawkes_unmarked.jl") + end + @testset verbose = true "Univariate" begin + include("hawkes_univariate.jl") + end + @testset verbose = true "Multivariate" begin + include("hawkes_multivariate.jl") + end end @testset verbose = true "PPTests" begin include("pp_tests.jl")