IHSEP (version 0.3.1)

mloglik1d: Minus loglikelihood of an IHSEP model

Description

Calculates the minus loglikelihood of an IHSEP model with given baseline intensity function \(\nu\) and excitation function \(g(x)=\sum a_i exp(-b_i x)\) for event times jtms on interval [0,TT].

Usage

mloglik1d(jtms, TT, nu, gcoef, Inu)

Value

The value of the negative log-liklihood.

Arguments

jtms

A numeric vector, with values sorted in ascending order. Jump times to fit the inhomogeneous self-exciting point process model on.

TT

A scalar. The censoring time, or the terminal time for observation. Should be (slightly) greater than the maximum of jtms.

nu

A (vectorized) function. The baseline intensity function.

gcoef

A numeric vector (of 2k elements), giving the parameters (a_1,...,a_k,b_1,...,b_k) of the exponential excitation function \(g(x)=\sum_{i=1}^k a_i*exp(-b_i*x)\).

Inu

A (vectorized) function. Its value at t gives the integral of the baseline intensity function \(\nu\) from 0 to t.

Author

Feng Chen <[e-mail]>

Details

This function calculates the minus loglikelihood of the inhomegeneous Hawkes model with background intensity function \(\nu(t)\) and excitation kernel function \(g(t)=\sum_{i=1}^{k} a_i e^{-b_i t}\) relative to continuous observation of the process from time 0 to time TT. Like mloglik1c, it takes advantage of the Markovian property of the intensity process and uses external C++ code to speed up the computation.

See Also

mloglik1c

Examples

## simulated data of an IHSEP on [0,1] with baseline intensity function
## nu(t)=200*(2+cos(2*pi*t)) and excitation function
## g(t)=8*exp(-16*t)
data(asep)

## get the birth times of all generations and sort in ascending order 
tms <- sort(unlist(asep))
## calculate the minus loglikelihood of an SEPP with the true parameters 
mloglik1d(tms,TT=1,nu=function(x)200*(2+cos(2*pi*x)),
          gcoef=8*1:2,
          Inu=function(y)integrate(function(x)200*(2+cos(2*pi*x)),0,y)$value)
## calculate the MLE for the parameter assuming known parametric forms
## of the baseline intensity and excitation functions  
if (FALSE) {
system.time(est <- optim(c(400,200,2*pi,8,16),
                         function(p){
                           mloglik1d(jtms=tms,TT=1,
                                     nu=function(x)p[1]+p[2]*cos(p[3]*x),
                                     gcoef=p[-(1:3)],
                                     Inu=function(y){
                                      integrate(function(x)p[1]+p[2]*cos(p[3]*x),
                                                0,y)$value
                                     })
                         },hessian=TRUE,control=list(maxit=5000,trace=TRUE),
                         method="BFGS")
            )
## point estimate by MLE
est$par
## standard error estimates:
diag(solve(est$hessian))^0.5
}