knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
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)
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)
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 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)
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)))
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)
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.
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)
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)
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)
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.