mloglik1a: Minus loglikelihood of an IHSEP model

In IHSEP: Inhomogeneous Self-Exciting Process

Minus loglikelihood of an IHSEP model

Description

Calculates the minus loglikelihood of an IHSEP model with given baseline inensity function ν and excitation function g for event times jtms on interval [0,TT].

Usage

mloglik1a(jtms, TT = max(jtms), nu, g,
          Ig = function(x) sapply(x, function(y) integrate(g, 0, y)$value),
          tol.abs = 1e-12, tol.rel = 1e-12, limit = 1000
          )

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.
g	A (vectorized) function. The excitation function.
Ig	A (vectorized) function. Its value at t gives the integral of the excitation function from 0 to t.
tol.abs	A small positive number. The tolerance of the absolute error in the numerical integral of ν.
tol.rel	A small positive number. The tolerance of the relative error in the numerical integral of ν.
limit	An (large) positive integer. The maximum number of subintervals allowed in the adaptive quadrature method to find the numerical integral of ν.

Details

This version of the mloglik function uses external C code to speedup the calculations. Otherwise it is the same as the mloglik0 function.

Value

The value of the negative log-liklihood.

Author(s)

Feng Chen <[e-mail]>

See Also

mloglik0

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 
mloglik1a(tms,TT=1,nu=function(x)200*(2+cos(2*pi*x)),
          g=function(x)8*exp(-16*x),Ig=function(x)8/16*(1-exp(-16*x)))
## calculate the MLE for the parameter assuming known parametric forms
## of the baseline intensity and excitation functions  
## Not run: 
system.time(est <- optim(c(400,200,2*pi,8,16),
                         function(p){
                           mloglik1a(jtms=tms,TT=1,
                                     nu=function(x)p[1]+p[2]*cos(p[3]*x),
                                     g=function(x)p[4]*exp(-p[5]*x),
                                     Ig=function(x)p[4]/p[5]*(1-exp(-p[5]*x)))
                         },
                         hessian=TRUE,control=list(maxit=5000,trace=TRUE))
            )
## point estimate by MLE
est$par
## standard error estimates:
diag(solve(est$hessian))^0.5

## End(Not run)

IHSEP documentation built on Sept. 17, 2022, 1:05 a.m.