knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )

```
library(emhawkes)
```

This subsection explaines how to construct, simulate, and estimate a univariate Hawkes model.
First, create a `hspec`

which defines the Hawkes model.
S4 class `hspec`

contains slots of model parameters, `mu`

, `alpha`

, `beta`

, `dimens`

,
`rmark`

, and `impact`

.

The parameters of the model, `mu`

, `alpha`

, `beta`

, are defined by matrices.
The following is an example of a univariate model (without mark).

set.seed(1107) mu1 <- 0.3; alpha1 <- 1.2; beta1 <- 1.5 hspec1 <- new("hspec", mu=mu1, alpha=alpha1, beta=beta1) show(hspec1)

The function `hsim`

implements simulation where the input arguments are `hspec`

, `size`

and the initial values of intensity `lambda`

and Hawkes processes `N`

.
If the initial values of `lambda`

are omitted, the internally determined initial values are used.
The default initial value of `N`

is zero.

res1 <- hsim(hspec1, size=100) summary(res1)

The results of `hsim`

is an S3 class `hreal`

which consists of `hspec`

, `inter_arrival`

, `arrival`

, `type`

, `mark`

, `N`

, `Nc`

, `lambda`

, `lambda_component`

, `rambda`

, `rambda\_component`

.
The components `inter_arrival`

, `type`

, `mark`

, `N`

, `Nc`

, `lambda`

, `lambda_component`

, `rambda`

, `rambda_component`

can be excessed during simulation and estimation.
The column names of `N`

are `N1`

, `N2`

, `N3`

and so on.
Similarly, the column names of `lambda`

are `lambda1`

, `lambda2`

, `lambda3`

and so on.
The column names of `lambda_component`

are `lambda_component11`

, `lambda_component12`

, `lambda_component13`

and so on.
`inter_arrival`

, `type`

, `mark`

, `N`

, `Nc`

start at zero.
Using `summary`

function, one can print the first 20 elements of `arrival`

, `N`

and `lambda`

.

`hspec`

is the model specification, `inter_arrival`

is the inter-arrival time of every event, and `arrival`

is the cumulative sum of `inter_arrival`

.
`type`

is the type of events, i.e., $i$ for $N_i$, and `mark`

is a numeric vector which represents additional information for the event.
`lambda`

represents $\bm{\lambda}$ which is the left continuous and right limit version.
The right continuous version of intensity is `rambda`

.
`lambda_component`

represents $\lambda_{ij}$ and `rambda_component`

is the right continuous version.

The log-likelihood function is computed by `logLik`

method.
In this case, the inter-arrival times and `hspec`

are inputs of the function.

logLik(hspec1, inter_arrival = res1$inter_arrival)

The likelihood estimation is performed using `mhfit`

function.
The specification of the initial values of the parameters, `mhspec0`

is needed.
In the following example, `hspec0`

is set to be `hspec1`

, which is defined previously, for simplicity, but any candidates for the starting value of the numerical procedure can be used.

Note that only `inter_arrival`

is needed.
(Indeed, for more precise simulation, `lambda0`

, the inital value of lambda compoment, should be specified. If not, internally determined initial values are set.)

mu0 <- 0.5; alpha0 <- 1.0; beta0 <- 1.8 hspec0 <- new("hspec", mu=mu0, alpha=alpha0, beta=beta0) mle <- hfit(hspec0, inter_arrival = res1$inter_arrival) summary(mle)

In a bivariate model, the parameters, the slots of `hspec`

, are matrices.
`mu`

is 2-by-1, and `alpha`

and `beta`

are 2-by-2 matrices.
`rmark`

is a random number generating function.
`lambda0`

, 2-by-2 matrix, represents the initial values of `lambda_component`

, a set of `lambda11, lambda12, lambda21, lambda22`

.
The intensity processes are represented by

$$ \lambda_1(t) = \mu_1 + \lambda_{11}(t) + \lambda_{12}(t) $$

$$ \lambda_2(t) = \mu_2 + \lambda_{21}(t) + \lambda_{22}(t) $$

$\lambda_{ij}$ called lambda components and `lambda0`

is the time zero values of $lambda_{ij}$, i.e., $\lambda_{ij}(0)$.
`lambda0`

can be omitted and then internally determined initial values are used.

mu2 <- matrix(c(0.2), nrow = 2) alpha2 <- matrix(c(0.5, 0.9, 0.9, 0.5), nrow = 2, byrow = TRUE) beta2 <- matrix(c(2.25, 2.25, 2.25, 2.25), nrow = 2, byrow = TRUE) lambda0 <- matrix(c(0.1, 0.1, 0.1, 0.1), nrow = 2, byrow = TRUE) hspec2 <- new("hspec", mu=mu2, alpha=alpha2, beta=beta2) print(hspec2)

To simulate, use function `mhsim`

.

res2 <- hsim(hspec2, size=1000) summary(res2)

Frome the result, we get a vector of realized `inter_arrival`

and `type`

.
Bivariate model requires `inter_arrival`

and `type`

for estimation.

inter_arrival2 <- res2$inter_arrival type2 <- res2$type

Log-likelihood is computed by a function `logLik`

.

logLik(hspec2, inter_arrival = inter_arrival2, type = type2)

A log-likelihood estimation is performed using `hfit`

.
In the following, the values of parameter slots in `mhspec0`

, such as `mu, alpha, beta`

, are regarded as a starting point of the numerical optimization.
For simplicity, we use `hspec0 <- hspec2`

.
Since the true parameter values are not known in the actual problem, the initial value should be guessed.
The realized `inter_arrival`

and `type`

are used.

hspec0 <- hspec2 mle <- hfit(hspec0, inter_arrival = inter_arrival2, type = type2) summary(mle)

This subsection explains about the relation between parameter setting and estimation procedure in two-dimensional Hawkes model.
The number of parameters to be estimated in the model depends on how we set the parameter slots such as `alpha`

and `beta`

in `mhspec0`

, the sepcification for initial values..
Since the paremeter slot such as `alpha`

is a matrix, and the element in the matrix can be the same or different.
The number of parameters in the estimation varies depending on whether or not some of the elements in the initial setting are the same or different.

For example, if `alpha[1,1]`

and `alpha[1,2]`

in `mhspec0`

are different,
the numerical procedure tries to estimate both parameters of `alpha[1,1]`

and `alpha[1,2]`

.
If `alpha[1,1]`

and `alpha[1,2]`

are the same in the initial setting,
then the estimation procedure considered two parameters are the same in the model and hence only one of them is estimated.

Recall that th example in the previous section is of a symmetric Hawkes model
where `alpha`

is symmetric.

```
print(hspec2)
```

res2 <- hsim(hspec2, size=1000)

In the first example of estimation, `alpha0`

is a matrix where the all elements have the same value, 0.75.
In this setting, `mhfit`

considers that `alpha11 == alpha12 == alpha21 == alpha22`

in the model
(even though the actual parameters have different values).
Similarly for other parmater matrix `mu0`

and `beta0`

.

mu0 <- matrix(c(0.15, 0.15), nrow = 2) alpha0 <- matrix(c(0.75, 0.75, 0.75, 0.75), nrow = 2, byrow=TRUE) beta0 <- matrix(c(2.6, 2.6, 2.6, 2.6), nrow = 2, byrow=TRUE) hspec0 <- new("hspec", mu=mu0, alpha=alpha0, beta=beta0) summary(hfit(hspec0, inter_arrival = res2$inter_arrival, type = res2$type))

In the second example, `alpha0`

's elements are not same, but symmetric as in the original simulation.
We have `alpha11 == alpha22`

and `alpha11 == alpha22`

in `alpha0`

and hence
`alpha11`

and `alpha12`

will be estimated.

mu0 <- matrix(c(0.15, 0.15), nrow = 2) alpha0 <- matrix(c(0.75, 0.751, 0.751, 0.75), nrow = 2, byrow=TRUE) beta0 <- matrix(c(2.6, 2.6, 2.6, 2.6), nrow = 2, byrow=TRUE) hspec0 <- new("hspec", mu=mu0, alpha=alpha0, beta=beta0) summary(hfit(hspec0, inter_arrival = res2$inter_arrival, type = res2$type))

By simply setting `reduced = FALSE`

, all parameters are estimated (not recommended).

summary(hfit(hspec2, inter_arrival = res2$inter_arrival, type = res2$type, reduced=FALSE))

In the following, `impact`

represents the impact of mark.
It is a function, and the first argument is `param`

represents the parameter of the model.
`impact`

function can have additional arguments such as `alpha`

, `n`

, `mark`

, etc., which represents the dependence.
Do not miss `...`

as the ellipsis is omitted, an error occurs.

mu <- matrix(c(0.15, 0.15), nrow=2) alpha <- matrix(c(0.75, 0.6, 0.6, 0.75), nrow=2, byrow=T) beta <- matrix(c(2.6, 2.6, 2.6, 2.6), nrow=2) rmark <- function(param = c(p=0.65), ...){ rgeom(1, p=param[1]) + 1 } impact <- function(param = c(eta1=0.2), alpha, n, mark, ...){ ma <- matrix(rep(mark[n]-1, 4), nrow = 2) alpha * ma * matrix( rep(param["eta1"], 4), nrow=2) #alpha * ma * matrix( c(rep(param["eta1"], 2), rep(param["eta2"], 2)), nrow=2) } h1 <- new("hspec", mu=mu, alpha=alpha, beta=beta, rmark = rmark, impact=impact) res <- hsim(h1, size=1000, lambda0 = matrix(rep(0.1,4), nrow=2)) fit <- hfit(h1, inter_arrival = res$inter_arrival, type = res$type, mark = res$mark, lambda0 = matrix(rep(0.1,4), nrow=2)) summary(fit)

In the basic model, `alpha`

is a matrix, but it can be a function as in the following code.
The function `alpha`

simply return a $4\times4$ matrix but by doing so, we can fix some of elements as specific vales when estimating.
When estimating, the optimization only works for the specified parameters in `param`

.
In the case of simulation, there is no difference whether the parameter set is represented by a matrix or a function.

mu <- matrix(c(0.02, 0.02, 0.04, 0.04), nrow = 4) alpha <- function(param = c(alpha11 = 0, alpha12 = 0.3, alpha33 = 0.3, alpha34 = 0.4)){ matrix(c(param["alpha11"], param["alpha12"], 0, 0, param["alpha12"], param["alpha11"], 0, 0, 0, 0, param["alpha33"], param["alpha34"], 0, 0, param["alpha34"], param["alpha33"]), nrow = 4, byrow = T) } beta <- matrix(c(rep(0.7, 8), rep(1.1, 8)), nrow = 4, byrow = T) impact <- function(param = c(alpha1n=0, alpha1w=0.1, alpha2n=0.1, alpha2w=0.2), n=n, N=N, ...){ Psi <- matrix(c(0, 0, param['alpha1w'], param['alpha1n'], 0, 0, param['alpha1n'], param['alpha1w'], param['alpha2w'], param['alpha2n'], 0, 0, param['alpha2n'], param['alpha2w'], 0, 0), nrow=4, byrow=T) ind <- N[,"N1"][n] - N[,"N2"][n] > N[,"N3"][n] - N[,"N4"][n] + 0.5 km <- matrix(c(!ind, !ind, !ind, !ind, ind, ind, ind, ind, ind, ind, ind, ind, !ind, !ind, !ind, !ind), nrow = 4, byrow = T) km * Psi } h <- new("hspec", mu = mu, alpha = alpha, beta = beta, impact = impact) hr <- hsim(h, size=1000) fit <- hfit(h, hr$inter_arrival, hr$type) summary(fit)

**Any scripts or data that you put into this service are public.**

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.