mloglik1e {IHSEP}

Minus loglikelihood of an IHSEP model

Description

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

Usage

mloglik1e(jtms, TT, nuvs, gcoef, InuT)

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.

nuvs

A numeric vector, giving the values of the baseline intensity function \nu at the jumptimes jtms.

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).

InuT

A numeric value (scalar) giving the integral of \nu on the interval [0,TT].

Details

This version of the mloglik function uses external C code to speedup the calculations. When given the analytical form of Inu or a quickly calculatable Inu, it should be (substantially) faster than mloglik1a when calculating the (minus log) likelihood when the excitation function is exponential. Otherwise it is the same as mloglik0, mloglik1a, mloglik1b.

Value

The value of the negative log-liklihood.

Author(s)

Feng Chen <[e-mail]>

See Also

mloglik0, mloglik1a and mloglik1b

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 
mloglik1e(tms,TT=1,nuvs=200*(2+cos(2*pi*tms)),
          gcoef=8*1:2,
          InuT=integrate(function(x)200*(2+cos(2*pi*x)),0,1)$value)
## 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){
                           mloglik1e(jtms=tms,TT=1,
                                     nuvs=p[1]+p[2]*cos(p[3]*tms),
                                     gcoef=p[-(1:3)],
                                     InuT=integrate(function(x)p[1]+p[2]*cos(p[3]*x),
                                                    0,1)$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

## End(Not run)