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