API reference

PointProcesses

A package for temporal point process modeling, simulation and inference.

History{T<:Real, M}

Linear event histories with temporal locations of type T and marks of type M.

Fields

• times::Vector{T}: sorted vector of event times
• tmin::T: start time
• tmax::T: end time
• marks::Vector{M}: associated vector of event marks

event_marks(h)

Return the vector of event marks for h, sorted according to their event times.

event_marks(h, tmin, tmax)

Return the sorted vector of marks of events between tmin and tmax in h.

event_times(h)

Return the sorted vector of event times for h.

event_times(h, tmin, tmax)

Return the sorted vector of event times between tmin and tmax in h.

min_time(h)

Return the starting time of h (not the same as the first event time).

max_time(h)

Return the end time of h (not the same as the last event time).

nb_events(h)

Count events in h.

nb_events(h, tmin, tmax)

Count events in h during the interval [tmin, tmax).

has_events(h)

Check the presence of events in h.

has_events(h, tmin, tmax)

Check the presence of events in h during the interval [tmin, tmax).

length(h)

Alias for nb_events(h).

length(pp::MultivariatePoissonProcess)

Return the number of marks (dimensions) in a multivariate Poisson process.

duration(h)

Compute the difference h.tmax - h.tmin.

min_mark(h; [init])

Return the smallest event mark if it is smaller than init, and init otherwise.

max_mark(h; [init])

Return the largest event mark if it is larger than init, and init otherwise.

push!(h, t, m)

Add event (t, m) inside the interval [h.tmin, h.tmax) at the end of history h.

append!(h, ts, ms)

Append events (ts, ms) inside the interval [h.tmin, h.tmax) at the end of history h.

cat(h1, h2)

If h1 and h2 are consecutive event histories, i.e., the end of h1 coincides with the beginning of h2, then create a new event history by concatenating h1 and h2.

time_change(h, Λ)

Apply the time rescaling t -> Λ(t) to history h.

split_into_chunks(h, chunk_duration)

Split h into a vector of consecutive histories with individual duration chunk_duration.

AbstractPointProcess

Common interface for all temporal point processes.

BoundedPointProcess{P,T} <: AbstractPointProcess{}

Temporal point process P with pre-defined start and end times.

Implements some fallbacks for the AbstractPointProcess interface which accept fewer arguments.

Fields

• pp::P: underlying point process
• tmin::T: start time
• tmax::T: end time

intensity(pp, m, t, h)

Compute the conditional intensity for a temporal point process pp applied to history h and event (t, m).

The conditional intensity function λ(t,m|h) quantifies the instantaneous risk of an event with mark m occurring at time t after history h.

ground_intensity(pp, h, t)

Compute the ground intensity for a temporal point process pp applied to history h at time t.

The ground intensity quantifies the instantaneous risk of an event with any mark occurring at time t after history h:

λg(t|h) = Σₘ λ(t,m|h)

log_intensity(pp, m, t, h)

Compute the logarithm of the conditional intensity for a temporal point process pp applied to history h and event (t, m).

intensity_vector(pp<:MultivariatePoissonProcess)

Compute the vector of the marginal intensities λ for a multivariate Poisson process.

mark_distribution(pp, t, h)

Compute the distribution of marks for a temporal point process pp knowing that an event takes place at time t after history h.

simulate_ogata(rng, pp, tmin, tmax)

Simulate a temporal point process pp on interval [tmin, tmax) using Ogata's algorithm.

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.

simulate([rng,] pp, tmin, tmax)

Alias for simulate_ogata.

simulate([rng,] bpp::BoundedPointProcess)

Simulate a point process on a predefined time interval.

logdensityof(pp, h)

Compute the log probability density function for a temporal point process pp applied to history h:

ℓ(h) = Σₖ log λ(tₖ|hₖ) - Λ(h)

The default method uses a loop over events combined with integrated_ground_intensity, but it should be reimplemented for specific processes if faster computation is possible.

integrated_ground_intensity(pp, h, a, b)

Compute the integrated ground intensity (or compensator) Λ(t|h) for a temporal point process pp applied to history h on interval [a, b):

Λ(h) = ∫ λg(t|h) dt

ground_intensity_bound(pp, t, h)

Compute a local upper bound on the ground intensity for a temporal point process pp applied to history h at time t.

Return a tuple of the form (B, L) satisfying λg(t|h) ≤ B for all u ∈ [t, t+L).

fit(::Type{PP}, h)
fit(::Type{PP}, histories)

Fit a point process of type PP to one or several histories.

Not implemented by default.

fit_map(::Type{PP}, h, prior)
fit_map(::Type{PP}, histories, prior)

Fit a point process of type PP to one or several histories using maximum a posteriori with a prior.

Not implemented by default.

PoissonProcess{R,D}

Homogeneous temporal Poisson process with arbitrary mark distribution.

Fields

• λ::R: ground intensity.
• mark_dist::D: mark distribution.

Constructor

PoissonProcess(λ, mark_dist)

UnivariatePoissonProcess{R}

Homogeneous univariate temporal Poisson process with scalar intensity λ::R.

UnivariatePoissonProcess{R} is simply a type alias for PoissonProcess{R,Dirac{Nothing}}.

MultivariatePoissonProcess{R}

Homogeneous multivariate temporal Poisson process with marginal intensities of type R.

MultivariatePoissonProcess{R} is simply a type alias for PoissonProcess{R,Categorical{Float64,Vector{Float64}}}.

MultivariatePoissonProcessPrior{R1,R2}

Gamma prior on all the event rates of a MultivariatePoissonProcess.

Fields

• α::Vector{R1}
• β::R2

InhomogeneousPoissonProcess{F,M,C}

Inhomogeneous temporal Poisson process with time-varying intensity.

Fields

• intensity_function::F: callable intensity function λ(t).
• mark_dist::M: mark distribution.
• integration_config::C: configuration for numerical integration.

Constructor

InhomogeneousPoissonProcess(intensity_function, mark_dist; integration_config=IntegrationConfig())

Examples

# Linear intensity
pp = InhomogeneousPoissonProcess(PolynomialIntensity([1.0, 0.5]), Normal())

# Sinusoidal intensity
pp = InhomogeneousPoissonProcess(SinusoidalIntensity(5.0, 2.0, 2π), Categorical([0.3, 0.7]))

# Custom intensity function
pp = InhomogeneousPoissonProcess(t -> 1.0 + 0.5*sin(t), Uniform())

# Custom integration settings
pp = InhomogeneousPoissonProcess(
    my_intensity,
    Normal(),
    integration_config=IntegrationConfig(abstol=1e-10)
)

ParametricIntensity

Abstract trait for intensity functions that can be parameterized for MLE fitting.

Any intensity type that implements this interface must provide:

• from_params(::Type{F}, params): Construct intensity function from parameters

The parameter space should be unconstrained (e.g., use log-transforms for positive parameters).

PolynomialIntensity{R<:Real,L} <: ParametricIntensity

Polynomial intensity function with optional link function.

Fields

• coefficients::Vector{R}: polynomial coefficients [a₀, a₁, ..., aₙ].
• link::L: link function applied to the polynomial (:identity or :log).

Constructor

PolynomialIntensity(coefficients; link=:identity)

When link=:identity: λ(t) = a₀ + a₁t + a₂t² + ... + aₙtⁿ When link=:log: λ(t) = exp(a₀ + a₁t + a₂t² + ... + aₙtⁿ)

The log link ensures positivity of the intensity function.

Examples

# Linear identity: λ(t) = 2 + 3*t (may be negative!)
PolynomialIntensity([2.0, 3.0])

# Linear log: λ(t) = exp(2 + 3*t) (always positive)
PolynomialIntensity([2.0, 3.0]; link=:log)

# Quadratic: λ(t) = 1 + 2*t + 0.5*t²
PolynomialIntensity([1.0, 2.0, 0.5])

ExponentialIntensity{R<:Real} <: ParametricIntensity

Exponential intensity function: λ(t) = aexp(bt).

Fields

• a::R: scaling factor (must be positive).
• b::R: exponential rate.

Constructor

ExponentialIntensity(a, b)

Examples

# Increasing intensity
ExponentialIntensity(2.0, 0.1)

# Decreasing intensity
ExponentialIntensity(5.0, -0.05)

SinusoidalIntensity{R<:Real} <: ParametricIntensity

Sinusoidal intensity function: λ(t) = a + bsin(ωt + φ).

To ensure positivity, we require a >= |b| so that λ(t) >= 0 for all t.

Fields

• a::R: baseline intensity (must satisfy a >= |b|).
• b::R: amplitude.
• ω::R: angular frequency.
• φ::R: phase shift.

Constructor

SinusoidalIntensity(a, b, ω, φ=0.0)

Examples

# Valid: a=5, b=2, so a >= |b|
SinusoidalIntensity(5.0, 2.0, 2π)

# Valid: a=5, b=-3, so a >= |-3| = 3
SinusoidalIntensity(5.0, -3.0, 2π)

# Invalid: a=2, b=3, so a < |b| (will error)

PiecewiseConstantIntensity{R<:Real}

Piecewise constant intensity function.

Fields

• breakpoints::Vector{R}: sorted vector of breakpoints (including tmin and tmax).
• rates::Vector{R}: intensity values for each interval.

Constructor

PiecewiseConstantIntensity(breakpoints, rates)

The intensity is rates[i] for t ∈ [breakpoints[i], breakpoints[i+1]).

LinearCovariateIntensity{R<:Real,F}

Linear combination of covariate functions: λ(t) = β₀ + β₁x₁(t) + β₂x₂(t) + ... + βₙ*xₙ(t).

Fields

• intercept::R: intercept term β₀.
• coefficients::Vector{R}: coefficients [β₁, β₂, ..., βₙ].
• covariates::Vector{F}: covariate functions [x₁, x₂, ..., xₙ], each callable with signature xᵢ(t).

Constructor

LinearCovariateIntensity(intercept, coefficients, covariates)

Examples

# With time and sin(time) as covariates
LinearCovariateIntensity(1.0, [0.5, 2.0], [t -> t, t -> sin(t)])

# With custom covariate functions
temp_func = t -> 20 + 5*sin(2π*t/365)  # seasonal temperature
wind_func = t -> 10 + 2*rand()          # wind speed
LinearCovariateIntensity(0.1, [0.05, 0.02], [temp_func, wind_func])

from_params

Method used in optimization, where parameters are returned as vectors. Can be used to perform a transformation in the parameter space. Example:

struct Constant{R} <: ParametricIntensity
    a::R
end where {R<:Real}
(f::Constant, t::Real) = f.a

# Optimization procedure uses this method to calculate the objective.
# This is not constrained to positive parameter values anymore
function from_params(::Constant, params)
    return Constant(exp(params[1]))
end

from_params(::Type{PolynomialIntensity{R}}, params; link=:identity)

Construct PolynomialIntensity from parameters.

from_params(::Type{ExponentialIntensity{R}}, params)

Construct ExponentialIntensity from unconstrained parameters: params = [log(a), b].

from_params(::Type{SinusoidalIntensity{R}}, params; ω=2π)

Construct SinusoidalIntensity from unconstrained parameters.

Parameters: [p₁, p₂, p₃] where a = exp(p₁), b = tanh(p₂)*a, φ = p₃. The ω (angular frequency) must be specified separately.

IntegrationConfig

Configuration for numerical integration of intensity functions.

Fields

• solver: Integration solver from Integrals.jl (e.g., QuadGKJL(), HCubatureJL())
• abstol::Float64: Absolute tolerance for integration
• reltol::Float64: Relative tolerance for integration
• maxiters::Int: Maximum number of iterations

Constructor

IntegrationConfig(; solver=QuadGKJL(), abstol=1e-8, reltol=1e-8, maxiters=1000)

Examples

# Default configuration
config = IntegrationConfig()

# Higher precision
config = IntegrationConfig(abstol=1e-12, reltol=1e-12)

# Different solver
using Integrals
config = IntegrationConfig(solver=HCubatureJL())

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 (Lewis and Mohler, 2011).

Statistic

Interface for all test statistics.

statistic(::Statistic, pp::AbstractPointProcess, h::History)

Compute the value of the test statistic with respect to pp and h.

Arguments

• ::Statistic: The type of test statistic to be computed
• pp::AbstractPointProcess: null-hypothesis model for the event history h
• h::History: the observed event history

Returns

• Float64: the resulting test statistic

Example

#=
Calculate the Kolmogorov-Smirnov distance between the distribution of the
time-changed event times of `h` and a standard exponential.
=#
ks_stat = statistic(KSDistance{Exponential}, hawkes_process, history)

PointProcessTest

Interface for all goodness-of-fit tests

pvalue(test::PointProcessTest)

Calculate the p-value of a goodness-of-fit test on a process.

Arguments

• ::PointProcessTest: the test result object

Returns

• Float64: p-value in [0, 1], where small values provide evidence against the null hypothesis

KSDistance{T<:UnivariateDistribution}

A Kolmogorov-Smirnov distance statistic for testing goodness-of-fit of point processes against a specified distribution D after appropriate time rescaling. This test statistic is suitable only for non-marked processes, as it ignores the marks.

Type parameter

• D<:UnivariateDistribution: the target distribution to test against (e.g., Exponential, Uniform)

Available test statistics

• KSDistance{Exponential} Kolmogorov-Smirnov distance between the time-changed interevent times and a standard exponential

• KSDistance{Uniform} Kolmogorov-Smirnov distance between the time-changed event times and a uniform distribution

Example

BootstrapTest(KSDistance{Exponential}, HawkesProcess, history)

BootstrapTest <: PointProcessTest

An object containing the results of a bootstrap-based goodness-of-fit test. The p-value of the test is calculated as p = (count(simstats ≥ stat) + 1) / (nsims + 1).

Fields

• n_sims::Int: number of bootstrap simulations performed
• stat::Float64: observed test statistic value
• sim_stats::Vector{Float64}: test statistics from bootstrap simulations

BootstrapTest(S::Type{<:Statistic}, pp::AbstractPointProcess, h::History; n_sims::Int=1000, rng::AbstractRNG=default_rng())

Perform a goodness-of-fit test using simulation with bootstrap resampling, comparing the test statistic computed on the observed data against the distribution of the same statistic computed on data simulated from the fitted model.

If λ₀(t) is the true intensity function of the process that generated the observed history, and λ(t; θ) is a a parametrization of the intensity, then the null hypothesis is

H₀: There exists parameters θₒ such that λ₀(t) = λ(t; θ₀)

This procedure is specifically aimed for testing hypotheses where parameters need to be estimated. Details are provided in (Kling and Vetter, 2025).

Arguments

• S::Type{<:Statistic}: the type of test statistic to use
• pp::Type{<:AbstractPointProcess}: the null hypothesis model family
• h::History: the observed event history
• n_sims::Int=1000: number of bootstrap simulations to perform
• rng::AbstractRNG=default_rng(): Random number generator

Returns

• BootstrapTest: test result object containing the observed statistic, bootstrap statistics, and test metadata

Example

# Bootstrap test for Hawkes process model adequacy
test = BootstrapTest(KSDistance(Exponential), HawkesProcess, history; n_sims=1000)
p = pvalue(test)

MonteCarloTest <: PointProcessTest

An object containing the results of a non-bootstrap based goodness-of-fit test. The p-value of the test is calculated as p = (count(simstats ≥ stat) + 1) / (nsims + 1).

Fields

• n_sims::Int: number of simulations performed
• stat::Float64: observed test statistic value
• sim_stats::Vector{Float64}: test statistics from simulated data

MonteCarloTest(S::Type{<:Statistic}, pp::AbstractPointProcess, h::History; n_sims::Int=1000, rng::AbstractRNG=default_rng())

Perform a goodness-of-fit test using simulation without bootstrap resampling, comparing the test statistic computed on the observed data against the distribution of the same statistic computed on data simulated from the fitted model.

If λ₀(t) is the true intensity function of the process that generated the observed history, and λ(t; θ) is a a parametrization of the intensity, then there are two forms for the null hypothesis:

1. H₀: λ₀(t) = λ(t; θ₀)
2. H₀: There exists parameters θₒ such that λ₀(t) = λ(t; θ₀)

If pp is an instance of an AbstractPointProcess, the null hypothesis 1 is considered, if a pp is a Type{<:AbstractPointProcess}, the method uses null hypothesis 2.

Notice that this test is better suited when the parameter θ₀ is known (form 1), since this procedure does not account for parameter estimation error. For more details on this, see (Babu and Rao, 2004), (Reynaud-Bouret et al., 2014), (Kling and Vetter, 2025).

Arguments

• S::Type{<:Statistic}: the type of test statistic to use
• pp::Union{AbstractPointProcess, Type{<:AbstractPointProcess}}: the null hypothesis model family
• h::History: the observed event history
• n_sims::Int=1000: number of simulations to perform for the test
• rng::AbstractRNG=default_rng(): Random number generator

Returns

• MonteCarloTest: test result object containing the observed statistic, simulated statistics, and test metadata

Example

# Test null hypothesis of form 1. Known θ₀
test = MonteCarloTest(KSDistance(Exponential), HawkesProcess(1, 1, 2), history; n_sims=1000)
p = pvalue(test)
# Test null hypothesis of form 2. Unknown θ₀
test = MonteCarloTest(KSDistance(Exponential), HawkesProcess, history; n_sims=1000)
p = pvalue(test)