Nothing
Example : emhakwes package "
In emhawkes: Exponential Multivariate Hawkes Model
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
Univariate Hawkes process
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)
Bivariate Hawkes model
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)
Parameter setting
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))
More complicated model
Linear impact model
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)
Hawkes flocking model
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)
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, 2021, 9:06 a.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 = "#>" )
Univariate Hawkes process
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)
Bivariate Hawkes model
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)
Parameter setting
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))
More complicated model
Linear impact model
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)
Hawkes flocking model
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)
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.