Nothing
Example : emhawkes package "
In emhawkes: Exponential Multivariate Hawkes Model
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
Basic Hawkes model
Univariate Hawkes process
library(emhawkes)
This subsection outlines the steps for constructing, running simulations, and estimating a univariate Hawkes model.
To begin, create an hspec object which defines the Hawkes model.
S4 class hspec contains slots for the model parameters: mu, alpha, beta, dimens, rmark, and impact.
In a univariate model, the basic parameters of the model, mu, alpha, beta, can be given as numeric.
If numeric values are given, they will be converted to matrices.
Below is an example of a univariate Hawkes model without a 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 component process, lambda_component0, and the initial values of Hawkes processes, N0.
More precisely, the intensity process the basic univariate Hawkes model is represented by $$ \lambda(t) = \mu + \int_{-\infty}^t \alpha e^{-\beta (t-s)} d N(s) = \mu + \lambda_c(0) e^{-\beta t} + \int_0^t \alpha e^{-\beta (t-s)} d N(s) $$ where the lambda_component0 denotes $$ \lambda_c(0) = \int_{-\infty}^0 \alpha e^{\beta s} d N(s).$$ If lambda_component0 is not provided, the internally determined initial values for intensity process are used.
If size is sufficiently large, exact value of lambda_component0 may not be important.
The default initial value of counting process, N0, is zero.
res1 <- hsim(hspec1, size = 1000)
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.
-
hspec is the model specification
-
inter_arrival is the inter-arrival time of every event
-
arrival is the cumulative sum of inter_arrival
-
type is the type of events, i.e., $i$ for $N_i$, and used for multivariate model
-
mark is a numeric vector which represents additional information for the event
-
lambda represents $\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.
inter_arrival, type, mark, N, and Nc start at zero.
Using summary() function, one can print the first 20 elements of arrival, N and lambda.
print() function also can be used.
By the definition, we have lambda == mu + lambda_compoent:
# first and third columns are the same
cbind(res1$lambda[1:5], res1$lambda_component[1:5], mu1 + res1$lambda_component[1:5])
Except the first row, rambda == lambda + alpha .
# second and third columns are the same
cbind(res1$lambda[1:5], res1$rambda[1:5], res1$lambda[1:5] + alpha1)
Also check the exponential decaying:
# By definition, the following two are equal:
res1$lambda[2:6]
mu1 + (res1$rambda[1:5] - mu1) * exp(-beta1 * res1$inter_arrival[2:6])
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, hspec0 is needed.
Note that only inter_arrival is needed in this univariate model.
(Indeed, for more precise simulation, lambda0, the initial value of lambda component, should be specified. If not, internally determined initial values are used.) By default, it uses the BFGS method for numerical optimization.
# initial value for numerical optimization
mu0 <- 0.5; alpha0 <- 1.0; beta0 <- 1.8
hspec0 <- new("hspec", mu = mu0, alpha = alpha0, beta = beta0)
# the intial values are provided through hspec
mle <- hfit(hspec0, inter_arrival = res1$inter_arrival)
summary(mle)
Bivariate Hawkes model
The intensity process of basic bivariate Hawkes model is defined by
$$
\lambda_1(t) = \mu_1 + \int_{-\infty}^t \alpha_{11} e^{-\beta_{11}(t-s)} d N_1(s) + \int_{-\infty}^t \alpha_{12} e^{-\beta_{12}(t-s)} d N_2(s),
$$
$$
\lambda_2(t) = \mu_2 + \int_{-\infty}^t \alpha_{21} e^{-\beta_{21}(t-s)} d N_1(s) + \int_{-\infty}^t \alpha_{22} e^{-\beta_{22}(t-s)} d N_2(s).
$$
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:
$$
\mu = \begin{bmatrix} \mu_1 \ \mu_2 \end{bmatrix}, \quad
\alpha = \begin{bmatrix} \alpha_{11} & \alpha_{12} \ \alpha_{21} & \alpha_{22} \end{bmatrix}, \quad
\beta =
\begin{bmatrix} \beta_{11} & \beta_{12} \ \beta_{21} & \beta_{22} \end{bmatrix}
$$
rmark is a random number generating function for mark and is not used for non-mark model.
lambda_component0, 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 $\lambda_{ij}(0)$.
lambda_component0 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)
hspec2 <- new("hspec", mu=mu2, alpha=alpha2, beta=beta2)
print(hspec2)
To simulate, use function hsim.
res2 <- hsim(hspec2, size=1000)
summary(res2)
type is crucial in multi-variate models, which represents the type of event.
# Under bi-variate model, there are two types, 1 or 2.
res2$type[1:10]
The column names of N are N1, N2, N3 and so on, for multivariate models.
res2$N[1:3, ]
Similarly, the column names of lambda are lambda1, lambda2, lambda3 and so on.
res2$lambda[1:3, ]
The column names of lambda_component are lambda_component11, lambda_component12, lambda_component13 and so on.
res2$lambda_component[1:3, ]
By definition, the following two are the same:
mu2[1] + rowSums(res2$lambda_component[1:5, c("lambda11", "lambda12")])
res2$lambda[1:5, "lambda1"]
From the result, we get vectors 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 maximum log-likelihood estimation is performed using hfit.
In the following, the values of parameter slots in hspec0, such as mu, alpha, beta, serve as starting points of the numerical optimization.
For the purpose of illustration, we use hspec0 <- hspec2.
Since the true parameter values are not known in practical applications, the initial guess is used.
The realized inter_arrival and type are used for estimation.
hspec0 <- hspec2
mle <- hfit(hspec0, inter_arrival = inter_arrival2, type = type2)
summary(mle)
coef(mle)
miscTools::stdEr(mle)
Parameter setting
This subsection covers about the relation between parameter setting and estimation procedure in multi-variate Hawkes model.
The number of parameters to be estimated in the model depends on how we set the parameter slots such as mu, alpha and beta in hspec0, the specification for initial values.
Since the parameter 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 hspec0 are different in initial starting, the numerical procedure tries to estimate both parameters of alpha[1,1] and alpha[1,2] differently.
If alpha[1,1] and alpha[1,2] are the same in the initial setting, then the estimation procedure considered two parameters are identical in the model and hence only one value is estimated.
Recall that the example in the previous section is of a symmetric Hawkes model where the matrix alpha is symmetric.
In addition, the elements of beta are all the same.
print(hspec2)
res2 <- hsim(hspec2, size = 1000)
In the first example of estimation, the initial value of alpha0 is a matrix where the all elements have the same value of 0.75.
In this configuration, hfit assumes that alpha11 == alpha12 == alpha21 == alpha22 in the model (even if the actual parameters have different values).
Similarly, the other parameter matrices mu0 and beta0 are also treated in the same way.
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))
Note that in the above result, alpha1.1 is somewhere between original alpha1.1 = 0.5 and alpha1.2 = 0.9.
In the following second example, alpha0's elements are not same, but symmetric as in the original values in the simulation.
We have alpha11 == alpha22 and alpha11 == alpha22 in alpha0 and hence alpha11 and alpha12 will be estimated differently.
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))
In the third example, the estimation is performed under the assumption that mu1 and mu2 may also be different (even though they are the same in the original model).
mu0 <- matrix(c(0.15, 0.14), 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))
The same logic is applied to all the higher dimensional model.
Residual process
Residual process can be extracted by residual_process().
inter_arrival, type, rambda_component, mu,beta should be provided.
The component denotes the type of the process to be extracted for multivariate model.
For example, for a bi-variate model, we have $N_1$ and \$N_2\$.
component=1 is for the residual of $N_1$ and component=2 is for the residual of $N_2$.
hrp <- new("hspec", mu = 0.3, alpha = 1.2, beta = 1.5)
res_rp <- hsim(hrp, size = 1000)
rp <- residual_process(component = 1,
inter_arrival = res_rp$inter_arrival, type = res_rp$type,
rambda_component = res_rp$rambda_component,
mu = 0.3, beta = 1.5)
p <- ppoints(100)
q <- quantile(rp,p=p)
plot(qexp(p), q, xlab="Theoretical Quantiles",ylab="Sample Quantiles")
qqline(q, distribution=qexp,col="blue", lty=2)
In case that rambda_component is unknown, it can be inferred by infer_lambda().
For infer_lambda(), the model hspec, inter_arrival and type are required.
The above example is then:
# estimation
mle_rp <- hfit(new("hspec", mu = 0.2, alpha = 1, beta = 2),
res_rp$inter_arrival)
# construct hspec from estimation result
he <- new("hspec", mu = coef(mle_rp)["mu1"],
alpha = coef(mle_rp)["alpha1"], beta = coef(mle_rp)["beta1"])
# infer intensity
infered_res <- infer_lambda(he, res_rp$inter_arrival, res_rp$type)
# compute residuals where we use
rp2 <- residual_process(component = 1,
inter_arrival = res_rp$inter_arrival, type = res_rp$type,
rambda_component = infered_res$rambda_component,
mu = coef(mle_rp)["mu1"], beta = coef(mle_rp)["beta1"])
p <- ppoints(100)
q <- quantile(rp2, p=p)
plot(qexp(p), q, xlab="Theoretical Quantiles",ylab="Sample Quantiles")
qqline(q, distribution=qexp,col="blue", lty=2)
More complicated model
Multi-kernel model
In a multi-kernel Hawkes model, type_col_map is required for hspec.
type_col_map is a list that represents the mapping between type and column number.
For example, consider a bi-variate multi-kernel model: $$
\lambda_t = \mu + \int_{-\infty}^{t} h(t-u) d N(u)
$$ where $$ h = \sum_{k=1}^{K} h_k, \quad
h_k (t) = \alpha_k \circ \begin{bmatrix}
e^{-\beta_{k11} t} & e^{-\beta_{k12} t} \
e^{-\beta_{k21} t} & e^{-\beta_{k22} t}
\end{bmatrix}
$$
with matrix $\alpha_k$ and $k$ denoting kernel number.
For example, in a bi-variate Hawkes model with two kernels, the intensity processes are
$$
\begin{bmatrix} \lambda_1(t) \ \lambda_2(t) \end{bmatrix} =
\begin{bmatrix} \mu_1 \ \mu_2 \end{bmatrix} + \int_{-\infty}^{t} \begin{bmatrix}
\alpha_{111} e^{-\beta_{111} t} & \alpha_{112} e^{-\beta_{112} t} \
\alpha_{121}e^{-\beta_{121} t} & \alpha_{122}e^{-\beta_{122} t}
\end{bmatrix} \begin{bmatrix} d N_1(s) \ dN_2(s) \end{bmatrix}
+ \int_{-\infty}^{t} \begin{bmatrix}
\alpha_{211} e^{-\beta_{211} t} & \alpha_{212} e^{-\beta_{212} t} \
\alpha_{221}e^{-\beta_{221} t} & \alpha_{222}e^{-\beta_{222} t}
\end{bmatrix} \begin{bmatrix} d N_1(s) \ dN_2(s) \end{bmatrix}.
$$
The parameter matrix is defined by
$$
\alpha = \begin{bmatrix} \alpha_{111} & \alpha_{112} & \alpha_{211} & \alpha_{212} \ \alpha_{121} & \alpha_{122} & \alpha_{221} & \alpha_{222}
\end{bmatrix}, \quad
\beta = \begin{bmatrix} \beta_{111} & \beta_{112} & \beta_{211} & \beta_{212} \ \beta_{121} & \beta_{122} & \beta_{221} & \beta_{222}
\end{bmatrix} \quad
$$
and we should specify which columns of matrix are associated with which $N_i$.
mu <- matrix(c(0.02, 0.02), nrow=2)
beta_1 <- matrix(rep(10, 4), nrow=2)
beta_2 <- matrix(rep(1, 4), nrow=2)
beta <- cbind(beta_1, beta_2)
alpha_1 <- matrix(c(3, 2,
2, 3), nrow=2, byrow=TRUE)
alpha_2 <- matrix(c(0.3, 0.2,
0.2, 0.3), nrow=2, byrow=TRUE)
alpha <- cbind(alpha_1, alpha_2)
print(alpha)
Note that $d N_1(s)$ is multiplied by first and third columns of $\alpha$ and $dN_2(s)$ is multiplied by second and fourth columns of $\alpha$ and hence type_col_map is
type_col_map <- list(c(1,3), # columns for dN1
c(2,4)) # columns for dN2
type_col_map
where type i is associated with columns of type_col_map[[i]].
Thus,
cat("Part of alpha associated with N1: \n")
alpha[, type_col_map[[1]]] # associated with N1
cat("Part of alpha associated with N2: \n")
alpha[, type_col_map[[2]]] # associated with N2
cat("Part of beta associated with N1: \n")
beta[, type_col_map[[1]]] # associated with N1
cat("Part of beta associated with N2: \n")
beta[, type_col_map[[2]]] # associated with N2
h <- new("hspec", mu = mu, alpha = alpha, beta=beta, type_col_map = type_col_map)
h
In addition, lambda_component0 should be provided for simulation and estimation.
res_mk <- hsim(h, size = 2000,
# for an illustration purpose
lambda_component0 = matrix(seq(1, 1.7, 0.1), nrow = 2))
res_mk
summary(hfit(h, res_mk$inter_arrival, res_mk$type,
lambda_component0 = matrix(seq(1, 1.7, 0.1), nrow = 2)))
Synchronized intensity model
This model is basically two-kernel model and defined by little bit complicated reparameterization.
$$
\mu = \begin{bmatrix} \theta/(1 - \kappa)/2 + \tilde\theta/(1 + \kappa)/2 \
\theta/(1 - \kappa)/2 - \tilde\theta/(1 + \kappa)/2 \end{bmatrix}, \quad \theta = (\theta^- + \theta^+)/2,\quad \tilde\theta=(\theta^- -\theta^+)/2
$$
$$
\alpha = \begin{bmatrix} \zeta & \tilde\zeta & \zeta & -\tilde\zeta \
\zeta & -\tilde\zeta & \zeta & \tilde\zeta
\end{bmatrix}, \quad \zeta = (\eta + \nu) / 2, \quad \tilde \zeta = (\eta - \nu)/ 2
$$
$$
\beta = \begin{bmatrix} \beta_1 & \beta_2 & \beta_1 & \beta_2 \
\beta_1 & \beta_2 & \beta_1 & \beta_2 \end{bmatrix}, \quad \beta_1 = (\eta + \nu) / 2, \quad \beta_2 = (\eta - \nu)/2
$$
In order to handle complex re-parametrization, each slot is expressed as a function rather than a matrix.
The first argument param is a set of parameters.
mu <- function(param = c(theta_p = 0.15, theta_n = 0.21, kappa = 0.12)){
theta <- (param["theta_n"] + param["theta_p"])/2
theta_tl <- (param["theta_n"] - param["theta_p"])/2
matrix(c(theta/2/(1 - param["kappa"]) + theta_tl/2/(1 + param["kappa"]),
theta/2/(1 - param["kappa"]) - theta_tl/2/(1 + param["kappa"])), nrow=2)
}
alpha <- function(param = c(eta = 5, nu = 3)){
zeta <- (param["eta"] + param["nu"])/2
zeta_tl <- (param["eta"] - param["nu"])/2
matrix(c(zeta, zeta_tl, zeta, -zeta_tl,
zeta, -zeta_tl, zeta, zeta_tl), nrow=2, byrow=TRUE)
}
beta <- function(param = c(beta = 12, kappa = 0.12)){
beta1 <- param["beta"] * (1 - param["kappa"])
beta2 <- param["beta"] * (1 + param["kappa"])
matrix(c(beta1, beta2, beta1, beta2,
beta1, beta2, beta1, beta2), nrow = 2, byrow = TRUE)
}
type_col_map <- list(c(1,2), c(3,4))
h_sy <- new("hspec", mu = mu, alpha = alpha, beta = beta, type_col_map = type_col_map)
h_sy
# run simulation
res_sy <- hsim(h_sy, size = 2000, lambda_component0 = matrix(rep(1, 2 * 4), nrow=2))
summary(res_sy)
The estimation is based on function arguments param.
In addition, the initial values of the numerical optimization is the default values specified in param.
Note that the same name arguments are treated as the same parameter.
kappa is in both of mu and beta, but only one kappa appears in the estimation result.
fit_sy <- hfit(h_sy, inter_arrival=res_sy$inter_arrival,
type=res_sy$type,
lambda_component0 = matrix(rep(1, 2 * 4), nrow=2))
summary(fit_sy)
Extended model
The following family of extended multi-variate marked Hawkes models are implemented:
$$
\lambda(t) = \mu + \int_{(-\infty,t)\times E} h(t, u, z)M(du \times dz)
$$
where the kernel $h$ is represented by
$$
h(t, u, z) = (\alpha + g(t, z))\Gamma(t),
$$
and
-
$\alpha$ is a constant matrix,
-
$g(t, z)$ is additional impacts on intensities, which may depend on mark, or any information generated by underlying processes,
-
$\Gamma(t)$ is exponential decaying matrix such that $\Gamma_{ij}(t) = e^{-\beta_{ij}(t)}$,
-
$M$ denotes the random measures defined on the product of time and mark spaces.
Linear impact model
In the linear impact model,
$$
g(t, z) = \eta (z-1).
$$
impact represents $\Psi(z)$, the impact of mark on future intensity.
It is a function, and the first argument is param represents the parameter of the model.
impact() function can have additional arguments related to the model specification or generated path, such as n, mark, etc.
Do not miss ... as the ellipsis is omitted, an error occurs.
rmark() is a function that generate marks for simulation.
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)
ma * matrix( rep(param["eta1"], 4), nrow=2)
}
hi <- new("hspec", mu=mu, alpha=alpha, beta=beta,
rmark = rmark,
impact=impact)
hi
res_impact <- hsim(hi, size=1000, lambda_component0 = matrix(rep(0.1,4), nrow=2))
summary(res_impact)
fit <- hfit(hi,
inter_arrival = res_impact$inter_arrival,
type = res_impact$type,
mark = res_impact$mark,
lambda_component0 = matrix(rep(0.1,4), nrow=2))
summary(fit)
For a special case of linear impact function, the following implementation is recommended.
In a marked Hawkes model, the additional linear impact can be represented by slot eta.
In this model, the intensity process is
$$
\lambda(t) = \mu + \int_{(-\infty, t)\times E} (\alpha + \eta (z-1)) e^{-\beta(t-u)} M(dt \times dz).
$$
rmark <- function(param = c(p=0.65), ...){
rgeom(1, p=param[1]) + 1
}
h <- new("hspec", mu=0.15, alpha=0.7, beta=1.6, eta=0.3,
rmark = rmark)
h
res <- hsim(h, size = 1000)
summary(res)
fit <- hfit(h,
inter_arrival = res$inter_arrival,
type = res$type,
mark = res$mark)
summary(fit)
If you want to estimate the mark distribution also, then dmark slot that describes the density function of mark is required.
h_md <- h
h_md@dmark <- function(param = c(p = 0.1), n=n, mark=mark, ...){
dgeom(mark[n] - 1, prob = param["p"])
}
mle_md <- hfit(h_md,
inter_arrival = res$inter_arrival, type = res$type, mark = res$mark)
summary(mle_md)
Hawkes flocking model
The function $g$ is not necessarily depend on mark.
In the Hawkes flocking model, the kernel component is represented by$$
\alpha = \begin{bmatrix}\alpha_{11} & \alpha_{12} & 0 & 0 \\alpha_{12}& \alpha_{11} & 0 & 0 \0 & 0 & \alpha_{33} & \alpha_{34} \0 & 0 & \alpha_{34} & \alpha_{33} \end{bmatrix},
$$ $$
g = \begin{bmatrix} 0 & 0 & \alpha_{1w} 1_{{C_1(t) < C_2(t)}} & \alpha_{1n} 1_{{C_1(t) < C_2(t)}} \
0 & 0 & \alpha_{1n} 1_{{C_1(t) > C_2(t)}} & \alpha_{1w}1_{{C_1(t) > C_2(t)}} \
\alpha_{2w} 1_{{C_2(t) < C_1(t)}} & \alpha_{2n}1_{{C_2(t) < C_1(t)}} & 0 & 0 \
\alpha_{2n} 1_{{C_2(t) > C_1(t)}} & \alpha_{2w}1_{{C_2(t) > C_1(t)}} & 0 & 0 \end{bmatrix},
$$
where$$
C_1(t) = N_1(t) - N_2(t), \quad C_2(t) = N_3(t) - N_4(t).
$$
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 is only applied for the specified parameters in the argument 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.2, 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 = TRUE)
}
beta <- matrix(c(rep(0.7, 8), rep(1.1, 8)), nrow = 4, byrow = TRUE)
impact() function is little bit complicated, but it is nothing more than expressing the definition of the model to an R function.
Note that we specify N=N, n=n in the argument.
N is for counting process $N$ and n denotes the time step.
Both are needed to implement the function body and it is required to specify in the argument.
… also should not be omitted.
impact <- function(param = c(alpha1n=0.25, 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=TRUE)
ind <- N[,"N1"][n] - N[,"N2"][n] > N[,"N3"][n] - N[,"N4"][n]
km <- matrix(c(!ind, !ind, !ind, !ind,
ind, ind, ind, ind,
ind, ind, ind, ind,
!ind, !ind, !ind, !ind), nrow = 4, byrow = TRUE)
km * Psi
}
hspec_fl <- new("hspec",
mu = mu, alpha = alpha, beta = beta, impact = impact)
hspec_fl
hr_fl <- hsim(hspec_fl, size=1000)
summary(hr_fl)
fit_fl <- hfit(hspec_fl, hr_fl$inter_arrival, hr_fl$type)
summary(fit_fl)
Bid-ask price model
In this model, we use a system of counting processes with the corresponding conditional intensities to describe the bid-ask price processes:
$$ N_t = \begin{bmatrix}
N_1(t) \ N_2(t) \ N_3(t) \ N_4(t)
\end{bmatrix},
\quad
\lambda_t =
\begin{bmatrix} \lambda_1(t) \ \lambda_2(t) \ \lambda_3(t) \ \lambda_4(t)
\end{bmatrix}
$$
The ask price process $N_1(t) - N_2(t)$ and the bid price process is $N_3(t) - N_4(t)$.
The mid price process is $p(t) = N_1(t) + N_3(t) - N_2(t) - N_4(t)$ plus initial mid price level.
The base intensity process is
$$\mu = \begin{bmatrix} \mu_1 \ \zeta \ell(t-) \ \zeta \ell(t-) \ \mu_1 \end{bmatrix}, \quad \ell(t) = \frac{L(t)}{p(t)} $$
where $L(t) \in { 0, 1, 2, \cdots }$ is the absolute level of the bid-ask spread with $L(t)=0$ implying the minimum level.
Note that in the following code of the definition of mu, n is needed to represent time $t$ and Nc is needed to calculate the level and mid price.
# presumed initial bid and ask prices
initial_ask_price <- 1250 #cents
initial_bid_price <- 1150 #cents
initial_level <- round((initial_ask_price - initial_bid_price) - 1)
initial_mid_price <- (initial_bid_price + initial_ask_price) / 2
mu <- function(param = c(mu1 = 0.08, zeta1 = 0.10), n=n, Nc=Nc, ...){
if(n == 1){
level <- initial_level
mid <- initial_mid_price
} else {
level <- Nc[n-1,1] - Nc[n-1,2] - (Nc[n-1,3] - Nc[n-1,4]) + initial_level
ask <- initial_ask_price + (Nc[n-1,1] - Nc[n-1,2])
bid <- initial_bid_price + (Nc[n-1,3] - Nc[n-1,4])
mid <- (ask + bid) / 2
}
if(level <= 0){
matrix(c(param["mu1"], 0,
0, param["mu1"]), nrow = 4)
} else {
matrix(c(param["mu1"], param["zeta1"] * level / mid,
param["zeta1"]*level / mid, param["mu1"]), nrow = 4)
}
}
In addition, the kernel is represented by
$$h(t, u) =
\begin{bmatrix} \alpha_{s1} & \alpha_{m} & \alpha_{s2} & 0 \ \alpha_{w1} & \alpha_{n1}(u) & \alpha_{n1}(u) & \alpha_{w2} \ \alpha_{w2} & \alpha_{n2}(u) & \alpha_{n2}(u) & \alpha_{w1} \ 0 & \alpha_{s2} & \alpha_{m} & \alpha_{s1} \ \end{bmatrix},
$$
where
$$
\alpha_{n1}(u) = - \sum_{j=1}^4 \lambda_{2j}(u) + \xi \ell(u), \quad \alpha_{n2}(u) = - \sum_{j=1}^4 \lambda_{3j}(u) + \xi \ell(u),
$$
for constant $\xi \geq 0$ and $\lambda_{ij}$ is a component of $\lambda_i$ such that
$$\lambda_{ij}(t) = \int_{-\infty}^t h_{ij}(t, u) d N_j(u).$$
In the following code, we separate the constant part of $h$ as alpha and stochastic part as impact.
To represent $\lambda_{ij}$, we need lambda_component_n.
Note that
alpha <- function(param = c(alpha_s1=4, alpha_m=26, alpha_s2=5,
alpha_w1=11, alpha_w2=7)){
matrix(c(param["alpha_s1"], param["alpha_m"], param["alpha_s2"], 0,
param["alpha_w1"], 0, 0, param["alpha_w2"],
param["alpha_w2"], 0, 0, param["alpha_w1"],
0, param["alpha_s2"], param["alpha_m"], param["alpha_s1"]), nrow = 4, byrow = TRUE)
}
impact <- function(param = c(xi = 2.7), n=n, Nc=Nc, lambda_component = lambda_component, ... ){
if(n == 1){
level <- initial_level
# mid <- initial_mid_price
} else {
level <- Nc[n,1] - Nc[n,2] - (Nc[n,3] - Nc[n,4]) + initial_level
ask <- initial_ask_price + (Nc[n,1] - Nc[n,2])
bid <- initial_bid_price + (Nc[n,3] - Nc[n,4])
mid <- (ask + bid)/2
}
lambda_component_matrix <- matrix(lambda_component[n, ], nrow=4, byrow=TRUE)
l2 <- sum(lambda_component_matrix[2,]) # sum of second row
l3 <- sum(lambda_component_matrix[3,]) # sum of third row
im <- matrix(c(0, 0, 0, 0,
0, -l2 + param["xi"]*level/mid, -l2 + param["xi"]*level/mid, 0,
0, -l3 + param["xi"]*level/mid, -l3 + param["xi"]*level/mid, 0,
0, 0, 0, 0), nrow = 4, byrow = TRUE)
}
beta <- matrix(rep(50, 16), nrow = 4, byrow=TRUE)
rmark <- function(n=n, Nc=Nc, type, ...){
if(n == 1){
level <- initial_level
} else {
level <- Nc[n-1,1] - Nc[n-1,2] - (Nc[n-1,3] - Nc[n-1,4]) + initial_level
}
if (type[n] == 2 | type[n] == 3){
min(level, rgeom(1, p=0.65) + 1)
} else {
rgeom(1, p=0.65) + 1
}
}
h_ba <- new("hspec", mu = mu, alpha = alpha, beta = beta, impact=impact, rmark = rmark)
h_ba
hr_ba <- hsim(h_ba, size=1000, lambda_component0 = matrix(rep(1, 16), 4))
summary(hr_ba)
As a separate log-likelihood estimation performed, the parameter for mark distribution is not estimated.
mle_ba <- hfit(h_ba, inter_arrival = hr_ba$inter_arrival, type = hr_ba$type,
lambda_component0 = matrix(rep(1, 16), 4))
summary(mle_ba)
Try the emhawkes package in your browser
Any scripts or data that you put into this service are public.
emhawkes documentation built on Feb. 16, 2023, 9:02 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
Basic Hawkes model
Univariate Hawkes process
library(emhawkes)
This subsection outlines the steps for constructing, running simulations, and estimating a univariate Hawkes model.
To begin, create an hspec object which defines the Hawkes model.
S4 class hspec contains slots for the model parameters: mu, alpha, beta, dimens, rmark, and impact.
In a univariate model, the basic parameters of the model, mu, alpha, beta, can be given as numeric.
If numeric values are given, they will be converted to matrices.
Below is an example of a univariate Hawkes model without a 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 component process, lambda_component0, and the initial values of Hawkes processes, N0.
More precisely, the intensity process the basic univariate Hawkes model is represented by $$ \lambda(t) = \mu + \int_{-\infty}^t \alpha e^{-\beta (t-s)} d N(s) = \mu + \lambda_c(0) e^{-\beta t} + \int_0^t \alpha e^{-\beta (t-s)} d N(s) $$ where the lambda_component0 denotes $$ \lambda_c(0) = \int_{-\infty}^0 \alpha e^{\beta s} d N(s).$$ If lambda_component0 is not provided, the internally determined initial values for intensity process are used.
If size is sufficiently large, exact value of lambda_component0 may not be important.
The default initial value of counting process, N0, is zero.
res1 <- hsim(hspec1, size = 1000) 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.
-
hspecis the model specification -
inter_arrivalis the inter-arrival time of every event -
arrivalis the cumulative sum ofinter_arrival -
typeis the type of events, i.e., $i$ for $N_i$, and used for multivariate model -
markis a numeric vector which represents additional information for the event -
lambdarepresents $\lambda$ which is the left continuous and right limit version -
The right continuous version of intensity is
rambda -
lambda_componentrepresents $\lambda_{ij}$ andrambda_componentis the right continuous version.
inter_arrival, type, mark, N, and Nc start at zero.
Using summary() function, one can print the first 20 elements of arrival, N and lambda.
print() function also can be used.
By the definition, we have lambda == mu + lambda_compoent:
# first and third columns are the same cbind(res1$lambda[1:5], res1$lambda_component[1:5], mu1 + res1$lambda_component[1:5])
Except the first row, rambda == lambda + alpha .
# second and third columns are the same cbind(res1$lambda[1:5], res1$rambda[1:5], res1$lambda[1:5] + alpha1)
Also check the exponential decaying:
# By definition, the following two are equal: res1$lambda[2:6] mu1 + (res1$rambda[1:5] - mu1) * exp(-beta1 * res1$inter_arrival[2:6])
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, hspec0 is needed.
Note that only inter_arrival is needed in this univariate model.
(Indeed, for more precise simulation, lambda0, the initial value of lambda component, should be specified. If not, internally determined initial values are used.) By default, it uses the BFGS method for numerical optimization.
# initial value for numerical optimization mu0 <- 0.5; alpha0 <- 1.0; beta0 <- 1.8 hspec0 <- new("hspec", mu = mu0, alpha = alpha0, beta = beta0) # the intial values are provided through hspec mle <- hfit(hspec0, inter_arrival = res1$inter_arrival) summary(mle)
Bivariate Hawkes model
The intensity process of basic bivariate Hawkes model is defined by
$$ \lambda_1(t) = \mu_1 + \int_{-\infty}^t \alpha_{11} e^{-\beta_{11}(t-s)} d N_1(s) + \int_{-\infty}^t \alpha_{12} e^{-\beta_{12}(t-s)} d N_2(s), $$
$$ \lambda_2(t) = \mu_2 + \int_{-\infty}^t \alpha_{21} e^{-\beta_{21}(t-s)} d N_1(s) + \int_{-\infty}^t \alpha_{22} e^{-\beta_{22}(t-s)} d N_2(s). $$
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:
$$ \mu = \begin{bmatrix} \mu_1 \ \mu_2 \end{bmatrix}, \quad \alpha = \begin{bmatrix} \alpha_{11} & \alpha_{12} \ \alpha_{21} & \alpha_{22} \end{bmatrix}, \quad \beta = \begin{bmatrix} \beta_{11} & \beta_{12} \ \beta_{21} & \beta_{22} \end{bmatrix} $$
rmark is a random number generating function for mark and is not used for non-mark model.
lambda_component0, 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 $\lambda_{ij}(0)$.
lambda_component0 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) hspec2 <- new("hspec", mu=mu2, alpha=alpha2, beta=beta2) print(hspec2)
To simulate, use function hsim.
res2 <- hsim(hspec2, size=1000) summary(res2)
type is crucial in multi-variate models, which represents the type of event.
# Under bi-variate model, there are two types, 1 or 2. res2$type[1:10]
The column names of N are N1, N2, N3 and so on, for multivariate models.
res2$N[1:3, ]
Similarly, the column names of lambda are lambda1, lambda2, lambda3 and so on.
res2$lambda[1:3, ]
The column names of lambda_component are lambda_component11, lambda_component12, lambda_component13 and so on.
res2$lambda_component[1:3, ]
By definition, the following two are the same:
mu2[1] + rowSums(res2$lambda_component[1:5, c("lambda11", "lambda12")]) res2$lambda[1:5, "lambda1"]
From the result, we get vectors 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 maximum log-likelihood estimation is performed using hfit.
In the following, the values of parameter slots in hspec0, such as mu, alpha, beta, serve as starting points of the numerical optimization.
For the purpose of illustration, we use hspec0 <- hspec2.
Since the true parameter values are not known in practical applications, the initial guess is used.
The realized inter_arrival and type are used for estimation.
hspec0 <- hspec2 mle <- hfit(hspec0, inter_arrival = inter_arrival2, type = type2) summary(mle)
coef(mle)
miscTools::stdEr(mle)
Parameter setting
This subsection covers about the relation between parameter setting and estimation procedure in multi-variate Hawkes model.
The number of parameters to be estimated in the model depends on how we set the parameter slots such as mu, alpha and beta in hspec0, the specification for initial values.
Since the parameter 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 hspec0 are different in initial starting, the numerical procedure tries to estimate both parameters of alpha[1,1] and alpha[1,2] differently.
If alpha[1,1] and alpha[1,2] are the same in the initial setting, then the estimation procedure considered two parameters are identical in the model and hence only one value is estimated.
Recall that the example in the previous section is of a symmetric Hawkes model where the matrix alpha is symmetric.
In addition, the elements of beta are all the same.
print(hspec2)
res2 <- hsim(hspec2, size = 1000)
In the first example of estimation, the initial value of alpha0 is a matrix where the all elements have the same value of 0.75.
In this configuration, hfit assumes that alpha11 == alpha12 == alpha21 == alpha22 in the model (even if the actual parameters have different values).
Similarly, the other parameter matrices mu0 and beta0 are also treated in the same way.
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))
Note that in the above result, alpha1.1 is somewhere between original alpha1.1 = 0.5 and alpha1.2 = 0.9.
In the following second example, alpha0's elements are not same, but symmetric as in the original values in the simulation.
We have alpha11 == alpha22 and alpha11 == alpha22 in alpha0 and hence alpha11 and alpha12 will be estimated differently.
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))
In the third example, the estimation is performed under the assumption that mu1 and mu2 may also be different (even though they are the same in the original model).
mu0 <- matrix(c(0.15, 0.14), 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))
The same logic is applied to all the higher dimensional model.
Residual process
Residual process can be extracted by residual_process().
inter_arrival, type, rambda_component, mu,beta should be provided.
The component denotes the type of the process to be extracted for multivariate model.
For example, for a bi-variate model, we have $N_1$ and \$N_2\$.
component=1 is for the residual of $N_1$ and component=2 is for the residual of $N_2$.
hrp <- new("hspec", mu = 0.3, alpha = 1.2, beta = 1.5) res_rp <- hsim(hrp, size = 1000) rp <- residual_process(component = 1, inter_arrival = res_rp$inter_arrival, type = res_rp$type, rambda_component = res_rp$rambda_component, mu = 0.3, beta = 1.5) p <- ppoints(100) q <- quantile(rp,p=p) plot(qexp(p), q, xlab="Theoretical Quantiles",ylab="Sample Quantiles") qqline(q, distribution=qexp,col="blue", lty=2)
In case that rambda_component is unknown, it can be inferred by infer_lambda().
For infer_lambda(), the model hspec, inter_arrival and type are required.
The above example is then:
# estimation mle_rp <- hfit(new("hspec", mu = 0.2, alpha = 1, beta = 2), res_rp$inter_arrival) # construct hspec from estimation result he <- new("hspec", mu = coef(mle_rp)["mu1"], alpha = coef(mle_rp)["alpha1"], beta = coef(mle_rp)["beta1"]) # infer intensity infered_res <- infer_lambda(he, res_rp$inter_arrival, res_rp$type) # compute residuals where we use rp2 <- residual_process(component = 1, inter_arrival = res_rp$inter_arrival, type = res_rp$type, rambda_component = infered_res$rambda_component, mu = coef(mle_rp)["mu1"], beta = coef(mle_rp)["beta1"]) p <- ppoints(100) q <- quantile(rp2, p=p) plot(qexp(p), q, xlab="Theoretical Quantiles",ylab="Sample Quantiles") qqline(q, distribution=qexp,col="blue", lty=2)
More complicated model
Multi-kernel model
In a multi-kernel Hawkes model, type_col_map is required for hspec.
type_col_map is a list that represents the mapping between type and column number.
For example, consider a bi-variate multi-kernel model: $$
\lambda_t = \mu + \int_{-\infty}^{t} h(t-u) d N(u)
$$ where $$ h = \sum_{k=1}^{K} h_k, \quad
h_k (t) = \alpha_k \circ \begin{bmatrix}
e^{-\beta_{k11} t} & e^{-\beta_{k12} t} \
e^{-\beta_{k21} t} & e^{-\beta_{k22} t}
\end{bmatrix}
$$
with matrix $\alpha_k$ and $k$ denoting kernel number.
For example, in a bi-variate Hawkes model with two kernels, the intensity processes are
$$ \begin{bmatrix} \lambda_1(t) \ \lambda_2(t) \end{bmatrix} = \begin{bmatrix} \mu_1 \ \mu_2 \end{bmatrix} + \int_{-\infty}^{t} \begin{bmatrix} \alpha_{111} e^{-\beta_{111} t} & \alpha_{112} e^{-\beta_{112} t} \ \alpha_{121}e^{-\beta_{121} t} & \alpha_{122}e^{-\beta_{122} t} \end{bmatrix} \begin{bmatrix} d N_1(s) \ dN_2(s) \end{bmatrix} + \int_{-\infty}^{t} \begin{bmatrix} \alpha_{211} e^{-\beta_{211} t} & \alpha_{212} e^{-\beta_{212} t} \ \alpha_{221}e^{-\beta_{221} t} & \alpha_{222}e^{-\beta_{222} t} \end{bmatrix} \begin{bmatrix} d N_1(s) \ dN_2(s) \end{bmatrix}. $$
The parameter matrix is defined by
$$ \alpha = \begin{bmatrix} \alpha_{111} & \alpha_{112} & \alpha_{211} & \alpha_{212} \ \alpha_{121} & \alpha_{122} & \alpha_{221} & \alpha_{222} \end{bmatrix}, \quad \beta = \begin{bmatrix} \beta_{111} & \beta_{112} & \beta_{211} & \beta_{212} \ \beta_{121} & \beta_{122} & \beta_{221} & \beta_{222} \end{bmatrix} \quad $$
and we should specify which columns of matrix are associated with which $N_i$.
mu <- matrix(c(0.02, 0.02), nrow=2) beta_1 <- matrix(rep(10, 4), nrow=2) beta_2 <- matrix(rep(1, 4), nrow=2) beta <- cbind(beta_1, beta_2) alpha_1 <- matrix(c(3, 2, 2, 3), nrow=2, byrow=TRUE) alpha_2 <- matrix(c(0.3, 0.2, 0.2, 0.3), nrow=2, byrow=TRUE) alpha <- cbind(alpha_1, alpha_2) print(alpha)
Note that $d N_1(s)$ is multiplied by first and third columns of $\alpha$ and $dN_2(s)$ is multiplied by second and fourth columns of $\alpha$ and hence type_col_map is
type_col_map <- list(c(1,3), # columns for dN1 c(2,4)) # columns for dN2 type_col_map
where type i is associated with columns of type_col_map[[i]].
Thus,
cat("Part of alpha associated with N1: \n") alpha[, type_col_map[[1]]] # associated with N1 cat("Part of alpha associated with N2: \n") alpha[, type_col_map[[2]]] # associated with N2 cat("Part of beta associated with N1: \n") beta[, type_col_map[[1]]] # associated with N1 cat("Part of beta associated with N2: \n") beta[, type_col_map[[2]]] # associated with N2
h <- new("hspec", mu = mu, alpha = alpha, beta=beta, type_col_map = type_col_map) h
In addition, lambda_component0 should be provided for simulation and estimation.
res_mk <- hsim(h, size = 2000, # for an illustration purpose lambda_component0 = matrix(seq(1, 1.7, 0.1), nrow = 2)) res_mk
summary(hfit(h, res_mk$inter_arrival, res_mk$type, lambda_component0 = matrix(seq(1, 1.7, 0.1), nrow = 2)))
Synchronized intensity model
This model is basically two-kernel model and defined by little bit complicated reparameterization.
$$ \mu = \begin{bmatrix} \theta/(1 - \kappa)/2 + \tilde\theta/(1 + \kappa)/2 \ \theta/(1 - \kappa)/2 - \tilde\theta/(1 + \kappa)/2 \end{bmatrix}, \quad \theta = (\theta^- + \theta^+)/2,\quad \tilde\theta=(\theta^- -\theta^+)/2 $$
$$ \alpha = \begin{bmatrix} \zeta & \tilde\zeta & \zeta & -\tilde\zeta \ \zeta & -\tilde\zeta & \zeta & \tilde\zeta \end{bmatrix}, \quad \zeta = (\eta + \nu) / 2, \quad \tilde \zeta = (\eta - \nu)/ 2 $$
$$ \beta = \begin{bmatrix} \beta_1 & \beta_2 & \beta_1 & \beta_2 \ \beta_1 & \beta_2 & \beta_1 & \beta_2 \end{bmatrix}, \quad \beta_1 = (\eta + \nu) / 2, \quad \beta_2 = (\eta - \nu)/2 $$
In order to handle complex re-parametrization, each slot is expressed as a function rather than a matrix.
The first argument param is a set of parameters.
mu <- function(param = c(theta_p = 0.15, theta_n = 0.21, kappa = 0.12)){ theta <- (param["theta_n"] + param["theta_p"])/2 theta_tl <- (param["theta_n"] - param["theta_p"])/2 matrix(c(theta/2/(1 - param["kappa"]) + theta_tl/2/(1 + param["kappa"]), theta/2/(1 - param["kappa"]) - theta_tl/2/(1 + param["kappa"])), nrow=2) } alpha <- function(param = c(eta = 5, nu = 3)){ zeta <- (param["eta"] + param["nu"])/2 zeta_tl <- (param["eta"] - param["nu"])/2 matrix(c(zeta, zeta_tl, zeta, -zeta_tl, zeta, -zeta_tl, zeta, zeta_tl), nrow=2, byrow=TRUE) } beta <- function(param = c(beta = 12, kappa = 0.12)){ beta1 <- param["beta"] * (1 - param["kappa"]) beta2 <- param["beta"] * (1 + param["kappa"]) matrix(c(beta1, beta2, beta1, beta2, beta1, beta2, beta1, beta2), nrow = 2, byrow = TRUE) } type_col_map <- list(c(1,2), c(3,4)) h_sy <- new("hspec", mu = mu, alpha = alpha, beta = beta, type_col_map = type_col_map) h_sy
# run simulation res_sy <- hsim(h_sy, size = 2000, lambda_component0 = matrix(rep(1, 2 * 4), nrow=2)) summary(res_sy)
The estimation is based on function arguments param.
In addition, the initial values of the numerical optimization is the default values specified in param.
Note that the same name arguments are treated as the same parameter.
kappa is in both of mu and beta, but only one kappa appears in the estimation result.
fit_sy <- hfit(h_sy, inter_arrival=res_sy$inter_arrival, type=res_sy$type, lambda_component0 = matrix(rep(1, 2 * 4), nrow=2)) summary(fit_sy)
Extended model
The following family of extended multi-variate marked Hawkes models are implemented:
$$ \lambda(t) = \mu + \int_{(-\infty,t)\times E} h(t, u, z)M(du \times dz) $$
where the kernel $h$ is represented by
$$ h(t, u, z) = (\alpha + g(t, z))\Gamma(t), $$
and
-
$\alpha$ is a constant matrix,
-
$g(t, z)$ is additional impacts on intensities, which may depend on mark, or any information generated by underlying processes,
-
$\Gamma(t)$ is exponential decaying matrix such that $\Gamma_{ij}(t) = e^{-\beta_{ij}(t)}$,
-
$M$ denotes the random measures defined on the product of time and mark spaces.
Linear impact model
In the linear impact model,
$$ g(t, z) = \eta (z-1). $$
impact represents $\Psi(z)$, the impact of mark on future intensity.
It is a function, and the first argument is param represents the parameter of the model.
impact() function can have additional arguments related to the model specification or generated path, such as n, mark, etc.
Do not miss ... as the ellipsis is omitted, an error occurs.
rmark() is a function that generate marks for simulation.
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) ma * matrix( rep(param["eta1"], 4), nrow=2) } hi <- new("hspec", mu=mu, alpha=alpha, beta=beta, rmark = rmark, impact=impact) hi
res_impact <- hsim(hi, size=1000, lambda_component0 = matrix(rep(0.1,4), nrow=2)) summary(res_impact)
fit <- hfit(hi, inter_arrival = res_impact$inter_arrival, type = res_impact$type, mark = res_impact$mark, lambda_component0 = matrix(rep(0.1,4), nrow=2)) summary(fit)
For a special case of linear impact function, the following implementation is recommended.
In a marked Hawkes model, the additional linear impact can be represented by slot eta.
In this model, the intensity process is
$$ \lambda(t) = \mu + \int_{(-\infty, t)\times E} (\alpha + \eta (z-1)) e^{-\beta(t-u)} M(dt \times dz). $$
rmark <- function(param = c(p=0.65), ...){ rgeom(1, p=param[1]) + 1 } h <- new("hspec", mu=0.15, alpha=0.7, beta=1.6, eta=0.3, rmark = rmark) h
res <- hsim(h, size = 1000) summary(res)
fit <- hfit(h, inter_arrival = res$inter_arrival, type = res$type, mark = res$mark) summary(fit)
If you want to estimate the mark distribution also, then dmark slot that describes the density function of mark is required.
h_md <- h h_md@dmark <- function(param = c(p = 0.1), n=n, mark=mark, ...){ dgeom(mark[n] - 1, prob = param["p"]) } mle_md <- hfit(h_md, inter_arrival = res$inter_arrival, type = res$type, mark = res$mark) summary(mle_md)
Hawkes flocking model
The function $g$ is not necessarily depend on mark. In the Hawkes flocking model, the kernel component is represented by$$ \alpha = \begin{bmatrix}\alpha_{11} & \alpha_{12} & 0 & 0 \\alpha_{12}& \alpha_{11} & 0 & 0 \0 & 0 & \alpha_{33} & \alpha_{34} \0 & 0 & \alpha_{34} & \alpha_{33} \end{bmatrix}, $$ $$ g = \begin{bmatrix} 0 & 0 & \alpha_{1w} 1_{{C_1(t) < C_2(t)}} & \alpha_{1n} 1_{{C_1(t) < C_2(t)}} \ 0 & 0 & \alpha_{1n} 1_{{C_1(t) > C_2(t)}} & \alpha_{1w}1_{{C_1(t) > C_2(t)}} \ \alpha_{2w} 1_{{C_2(t) < C_1(t)}} & \alpha_{2n}1_{{C_2(t) < C_1(t)}} & 0 & 0 \ \alpha_{2n} 1_{{C_2(t) > C_1(t)}} & \alpha_{2w}1_{{C_2(t) > C_1(t)}} & 0 & 0 \end{bmatrix}, $$
where$$ C_1(t) = N_1(t) - N_2(t), \quad C_2(t) = N_3(t) - N_4(t). $$
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 is only applied for the specified parameters in the argument 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.2, 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 = TRUE) } beta <- matrix(c(rep(0.7, 8), rep(1.1, 8)), nrow = 4, byrow = TRUE)
impact() function is little bit complicated, but it is nothing more than expressing the definition of the model to an R function.
Note that we specify N=N, n=n in the argument.
N is for counting process $N$ and n denotes the time step.
Both are needed to implement the function body and it is required to specify in the argument.
… also should not be omitted.
impact <- function(param = c(alpha1n=0.25, 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=TRUE) ind <- N[,"N1"][n] - N[,"N2"][n] > N[,"N3"][n] - N[,"N4"][n] km <- matrix(c(!ind, !ind, !ind, !ind, ind, ind, ind, ind, ind, ind, ind, ind, !ind, !ind, !ind, !ind), nrow = 4, byrow = TRUE) km * Psi } hspec_fl <- new("hspec", mu = mu, alpha = alpha, beta = beta, impact = impact) hspec_fl
hr_fl <- hsim(hspec_fl, size=1000) summary(hr_fl)
fit_fl <- hfit(hspec_fl, hr_fl$inter_arrival, hr_fl$type) summary(fit_fl)
Bid-ask price model
In this model, we use a system of counting processes with the corresponding conditional intensities to describe the bid-ask price processes:
$$ N_t = \begin{bmatrix} N_1(t) \ N_2(t) \ N_3(t) \ N_4(t) \end{bmatrix}, \quad \lambda_t = \begin{bmatrix} \lambda_1(t) \ \lambda_2(t) \ \lambda_3(t) \ \lambda_4(t) \end{bmatrix} $$
The ask price process $N_1(t) - N_2(t)$ and the bid price process is $N_3(t) - N_4(t)$. The mid price process is $p(t) = N_1(t) + N_3(t) - N_2(t) - N_4(t)$ plus initial mid price level.
The base intensity process is
$$\mu = \begin{bmatrix} \mu_1 \ \zeta \ell(t-) \ \zeta \ell(t-) \ \mu_1 \end{bmatrix}, \quad \ell(t) = \frac{L(t)}{p(t)} $$
where $L(t) \in { 0, 1, 2, \cdots }$ is the absolute level of the bid-ask spread with $L(t)=0$ implying the minimum level.
Note that in the following code of the definition of mu, n is needed to represent time $t$ and Nc is needed to calculate the level and mid price.
# presumed initial bid and ask prices initial_ask_price <- 1250 #cents initial_bid_price <- 1150 #cents initial_level <- round((initial_ask_price - initial_bid_price) - 1) initial_mid_price <- (initial_bid_price + initial_ask_price) / 2 mu <- function(param = c(mu1 = 0.08, zeta1 = 0.10), n=n, Nc=Nc, ...){ if(n == 1){ level <- initial_level mid <- initial_mid_price } else { level <- Nc[n-1,1] - Nc[n-1,2] - (Nc[n-1,3] - Nc[n-1,4]) + initial_level ask <- initial_ask_price + (Nc[n-1,1] - Nc[n-1,2]) bid <- initial_bid_price + (Nc[n-1,3] - Nc[n-1,4]) mid <- (ask + bid) / 2 } if(level <= 0){ matrix(c(param["mu1"], 0, 0, param["mu1"]), nrow = 4) } else { matrix(c(param["mu1"], param["zeta1"] * level / mid, param["zeta1"]*level / mid, param["mu1"]), nrow = 4) } }
In addition, the kernel is represented by
$$h(t, u) = \begin{bmatrix} \alpha_{s1} & \alpha_{m} & \alpha_{s2} & 0 \ \alpha_{w1} & \alpha_{n1}(u) & \alpha_{n1}(u) & \alpha_{w2} \ \alpha_{w2} & \alpha_{n2}(u) & \alpha_{n2}(u) & \alpha_{w1} \ 0 & \alpha_{s2} & \alpha_{m} & \alpha_{s1} \ \end{bmatrix}, $$
where
$$ \alpha_{n1}(u) = - \sum_{j=1}^4 \lambda_{2j}(u) + \xi \ell(u), \quad \alpha_{n2}(u) = - \sum_{j=1}^4 \lambda_{3j}(u) + \xi \ell(u), $$
for constant $\xi \geq 0$ and $\lambda_{ij}$ is a component of $\lambda_i$ such that
$$\lambda_{ij}(t) = \int_{-\infty}^t h_{ij}(t, u) d N_j(u).$$
In the following code, we separate the constant part of $h$ as alpha and stochastic part as impact.
To represent $\lambda_{ij}$, we need lambda_component_n.
Note that
alpha <- function(param = c(alpha_s1=4, alpha_m=26, alpha_s2=5, alpha_w1=11, alpha_w2=7)){ matrix(c(param["alpha_s1"], param["alpha_m"], param["alpha_s2"], 0, param["alpha_w1"], 0, 0, param["alpha_w2"], param["alpha_w2"], 0, 0, param["alpha_w1"], 0, param["alpha_s2"], param["alpha_m"], param["alpha_s1"]), nrow = 4, byrow = TRUE) } impact <- function(param = c(xi = 2.7), n=n, Nc=Nc, lambda_component = lambda_component, ... ){ if(n == 1){ level <- initial_level # mid <- initial_mid_price } else { level <- Nc[n,1] - Nc[n,2] - (Nc[n,3] - Nc[n,4]) + initial_level ask <- initial_ask_price + (Nc[n,1] - Nc[n,2]) bid <- initial_bid_price + (Nc[n,3] - Nc[n,4]) mid <- (ask + bid)/2 } lambda_component_matrix <- matrix(lambda_component[n, ], nrow=4, byrow=TRUE) l2 <- sum(lambda_component_matrix[2,]) # sum of second row l3 <- sum(lambda_component_matrix[3,]) # sum of third row im <- matrix(c(0, 0, 0, 0, 0, -l2 + param["xi"]*level/mid, -l2 + param["xi"]*level/mid, 0, 0, -l3 + param["xi"]*level/mid, -l3 + param["xi"]*level/mid, 0, 0, 0, 0, 0), nrow = 4, byrow = TRUE) } beta <- matrix(rep(50, 16), nrow = 4, byrow=TRUE) rmark <- function(n=n, Nc=Nc, type, ...){ if(n == 1){ level <- initial_level } else { level <- Nc[n-1,1] - Nc[n-1,2] - (Nc[n-1,3] - Nc[n-1,4]) + initial_level } if (type[n] == 2 | type[n] == 3){ min(level, rgeom(1, p=0.65) + 1) } else { rgeom(1, p=0.65) + 1 } } h_ba <- new("hspec", mu = mu, alpha = alpha, beta = beta, impact=impact, rmark = rmark) h_ba
hr_ba <- hsim(h_ba, size=1000, lambda_component0 = matrix(rep(1, 16), 4)) summary(hr_ba)
As a separate log-likelihood estimation performed, the parameter for mark distribution is not estimated.
mle_ba <- hfit(h_ba, inter_arrival = hr_ba$inter_arrival, type = hr_ba$type, lambda_component0 = matrix(rep(1, 16), 4)) summary(mle_ba)
Try the emhawkes package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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.