In ksublee/mhawkes: Simulate and estimate marked Hawkes process

Install pacakge

To install mHawkes package, first install devtools.

# install.packages("devtools")  #if devtools is not installed

Install mHawkes package from github.

#devtools::install_github("ksublee/mHawkes", build_vignettes=FALSE, force=TRUE)

Load mHawkes.

library("mHawkes")

one dimensional Hawkes process

First, create a mHSpec which defines the Hawkes model. Basically the model can be defined up to 9 dimension but currently fully supported for one and two diemsional model especillay for estimation. The parameter of the model is defined by a matrix but numeric vectors are also supported for one dimesional model.

The following is an example of one dimensional Hawkes model (without mark). Parameter inputs can be a numeric vector or 1 by 1 matrix. In the following case Jump and eta slots, which deteremine the mark, are default values.

MU1 <- 0.3
ALPHA1 <- 1.5
BETA1 <- 2
mHSpec1 <- new("mHSpec", MU=MU1, ALPHA=ALPHA1, BETA=BETA1)
show(mHSpec1)

To simulate, use function mHSim.

res1 <- mHSim(mHSpec1,  n=20)

The output res1 is a list of inter_arrival, arrival, jump_type, mark,N, Ng, lambda, and lambda_component. inter_arrival, arrival, jump_type, and mark are numeric vectors and N, Ng, lambda, and lambda_component are matrices.

res1

Note that the inter_arrival and arrival start at zero. Thus, inter_arrival[2] and arrival[2] are first arrival time.

Simle way to plot the realized processes:

# plot(res1$arrival, res1$N[,'N1'], 's', xlab="t", ylab="N")

Intensity process can be plotted by plotlambda function.

# plotlambda(res1$arrival, res1$N[,'N1'], BETA1)

The log-likelihood function is computed by logLik function. In this case, the inter-arrival times and mHSpec are inputs of the function.

inter_arrival1 <- res1$inter_arrival
logLik(mHSpec1, inter_arrival = inter_arrival1)

The likelihood estimation is performed using mHFit function. The specification of the initial values of the parameters, mHSpec0 is needed. In this example, mHSpec0 is set to be ``mHspec1 for simplicity.

mHSpec0 <- mHSpec1
mle <- mHFit(mHSpec0, inter_arrival = inter_arrival1)
summary(mle)

one dimensional Hawkes process with jump

Mark structure can be added with Jump slot. In addition, linear impact parameter ETA can be added.

ETA1 <- 0.15
JUMP1 <- function(n,...) rgeom(n, 0.65) + 1

mHSpec1 <- new("mHSpec", MU=MU1, ALPHA=ALPHA1, BETA=BETA1, ETA=ETA1, Jump =JUMP1)

res1 <- mHSim(mHSpec1,  n=10)

Plot the realized processes.

# plot(res1$arrival, res1$N[,'N1'], 's', xlab="t", ylab="N")

Two-dimensional Hawkes model

MU2 <- matrix(c(0.2), nrow = 2)
ALPHA2 <- matrix(c(0.75, 0.92, 0.92, 0.75), nrow = 2, byrow=TRUE)
BETA2 <- matrix(c(2.25, 2.25, 2.25, 2.25), nrow = 2, byrow=TRUE)
ETA2 <- matrix(c(0.19, 0.19, 0.19, 0.19), nrow = 2, byrow=TRUE)
JUMP2 <- function(n,...) rgeom(n, 0.65) + 1
LAMBDA0 <- matrix(c(0.1, 0.1, 0.1, 0.1), nrow = 2, byrow=TRUE)
mHSpec2 <- new("mHSpec", MU=MU2, ALPHA=ALPHA2, BETA=BETA2, ETA=ETA2, Jump =JUMP2)

mHSpec2

To simulate, use function mHSim.

res2 <- mHSim(mHSpec2,  n=100)
summary(res2)

Plot N.

# plot(res2$arrival[1:10], res2$N[1:10,1], 's')

# plot(res2)

Plot lambda.

# plotlambda(res2$arrival[1:10], res2$lambda[1:10,1], BETA2[1,1])

Frome the result, we get a vector of realized inter arrival times.

inter_arrival1 <- res1$inter_arrival

mle <- mHFit(mHSpec1, inter_arrival = inter_arrival1)

Frome the result, we get a vector of realized inter arrival times.

inter_arrival2 <- res2$inter_arrival
mark2 <- res2$mark
jump_type2 <- res2$jump_type

Log-likelihood is computed by a function logLik.

logLik(mHSpec2,  inter_arrival = inter_arrival2, jump_type = jump_type2, mark = mark2)

A log-likelihood estimation is performed using mHFit. mSpec0 is regarded as a starting point of the numerical optimization.

mHSpec0 <- mHSpec2
mle <- mHFit(mHSpec0, inter_arrival = inter_arrival2, jump_type = jump_type2, mark = mark2)
summary(mle)

The parameters to be estimated depends on mHSpec0. For example, if MU[1] and MU[2] are different, then both parameters are estimated. If MU[1] and MU[2] are equal, then two parameters are assumed to be the same.

MU0 <- matrix(c(0.2, 0.21), nrow = 2)
ALPHA0 <- matrix(c(0.75, 0.75, 0.75, 0.75), nrow = 2, byrow=TRUE)
BETA0 <- matrix(c(2.25, 2.251, 2.251, 2.25), nrow = 2, byrow=TRUE)
ETA0 <- matrix(c(0.19, 0.19, 0.19, 0.19), nrow = 2, byrow=TRUE)
JUMP0 <- function(n,...) rgeom(n, 0.65) + 1

mHSpec0 <- new("mHSpec", MU=MU0, ALPHA=ALPHA0, BETA=BETA0, ETA=ETA0, Jump =JUMP0)

mle <- mHFit(mHSpec0, inter_arrival = inter_arrival2, jump_type = jump_type2, mark = mark2)
summary(mle)

ksublee/mhawkes documentation built on May 3, 2019, 9:40 p.m.