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")
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)
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")
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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.