Perform maximum likelihood estimation

Description

Generic function hfit. A method for estimating the parameters of the exponential Hawkes model. The reason for being constructed as the S4 method is as follows. First, to represent the structure of the model as an hspec object. There are numerous variations on the multivariate marked Hawkes model. Second, to convey the starting point of numerical optimization. The parameter values assigned to the hspec slots become initial values. This function uses maxLik for the optimizer.

Usage

hfit(
  object,
  inter_arrival = NULL,
  type = NULL,
  mark = NULL,
  N = NULL,
  Nc = NULL,
  lambda_component0 = NULL,
  N0 = NULL,
  mylogLik = NULL,
  reduced = TRUE,
  grad = NULL,
  hess = NULL,
  constraint = NULL,
  method = "BFGS",
  verbose = FALSE,
  ...
)

## S4 method for signature 'hspec'
hfit(
  object,
  inter_arrival = NULL,
  type = NULL,
  mark = NULL,
  N = NULL,
  Nc = NULL,
  lambda_component0 = NULL,
  N0 = NULL,
  mylogLik = NULL,
  reduced = TRUE,
  grad = NULL,
  hess = NULL,
  constraint = NULL,
  method = "BFGS",
  verbose = FALSE,
  ...
)

Arguments

Value

maxLik object

See Also

hspec-class, hsim,hspec-method

Examples

# example 1
mu <- c(0.1, 0.1)
alpha <- matrix(c(0.2, 0.1, 0.1, 0.2), nrow=2, byrow=TRUE)
beta <- matrix(c(0.9, 0.9, 0.9, 0.9), nrow=2, byrow=TRUE)
h <- new("hspec", mu=mu, alpha=alpha, beta=beta)
res <- hsim(h, size=100)
summary(hfit(h, inter_arrival=res$inter_arrival, type=res$type))


# example 2

mu <- matrix(c(0.08, 0.08, 0.05, 0.05), nrow = 4)
alpha <- function(param = c(alpha11 = 0, alpha12 = 0.4, alpha33 = 0.5, alpha34 = 0.3)){
  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.6, 8), rep(1.2, 8)), nrow = 4, byrow = TRUE)

impact <- function(param = c(alpha1n=0, alpha1w=0.2, alpha2n=0.001, alpha2w=0.1),
                   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] + 0.5

  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
}
h <- new("hspec",
         mu = mu, alpha = alpha, beta = beta, impact = impact)
hr <- hsim(h, size=100)
plot(hr$arrival, hr$N[,'N1'] - hr$N[,'N2'], type='s')
lines(hr$N[,'N3'] - hr$N[,'N4'], type='s', col='red')
fit <- hfit(h, hr$inter_arrival, hr$type)
summary(fit)


# example 3

mu <- c(0.15, 0.15)
alpha <- matrix(c(0.75, 0.6, 0.6, 0.75), nrow=2, byrow=TRUE)
beta <- matrix(c(2.6, 2.6, 2.6, 2.6), nrow=2, byrow=TRUE)
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)
}
h1 <- new("hspec", mu=mu, alpha=alpha, beta=beta,
          rmark = rmark,
          impact=impact)
res <- hsim(h1, size=100, lambda_component0 = matrix(rep(0.1,4), nrow=2))

fit <- hfit(h1,
            inter_arrival = res$inter_arrival,
            type = res$type,
            mark = res$mark,
            lambda_component0 = matrix(rep(0.1,4), nrow=2))
summary(fit)

# For more information, please see vignettes.

Realization of Hawkes process

Description

hreal is the list of the following:

• hspec : S4 object hspec-class that specifies the parameter values.

• inter_arrival : the time between two consecutive events.

• arrival : cumulative sum of inter_arrival.

• type : integer, the type of event.

• mark : the size of mark, an additional information associated with event.

• N : counting process that counts the number of events.

• Nc : counting process that counts the number of events weighted by mark.

• lambda : intensity process, left-continuous version.

• lambda_component : the component of intensity process with mu not included.

• rambda : intensity process, right-continuous version.

• rambda_component : the right-continuous version of lambda_component.

Print functions for hreal are provided.

Usage

## S3 method for class 'hreal'
print(x, n = 20, ...)

## S3 method for class 'hreal'
summary(object, n = 20, ...)

Arguments

Simulate multivariate Hawkes process with exponential kernel.

Description

The method simulate multivariate Hawkes processes. The object hspec-class contains the parameter values such as mu, alpha, beta. The mark (jump) structure may or may not be included. It returns an object of class hreal which contains inter_arrival, arrival, type, mark, N, Nc, lambda, lambda_component, rambda, rambda_component.

Usage

hsim(object, size = 100, lambda_component0 = NULL, N0 = NULL, ...)

## S4 method for signature 'hspec'
hsim(object, size = 100, lambda_component0 = NULL, N0 = NULL, ...)

Arguments

Value

hreal S3-object, summary of the Hawkes process realization.

Examples

# example 1

mu <- 1; alpha <- 1; beta <- 2
h <- new("hspec", mu=mu, alpha=alpha, beta=beta)
hsim(h, size=100)


# example 2
mu <- matrix(c(0.1, 0.1), nrow=2)
alpha <- matrix(c(0.2, 0.1, 0.1, 0.2), nrow=2, byrow=TRUE)
beta <- matrix(c(0.9, 0.9, 0.9, 0.9), nrow=2, byrow=TRUE)
h <- new("hspec", mu=mu, alpha=alpha, beta=beta)
res <- hsim(h, size=100)
print(res)

An S4 class to represent an exponential marked Hawkes model

Description

This class represents a specification of a marked Hawkes model with exponential kernel. The intensity of the ground process is defined by:

\lambda(t) = \mu + \int_{(-\infty,t)\times E} ( \alpha + g(u, z) ) e^{-\beta (t-u)} M(du \times dz).

For more details, please see the vignettes.

Details

\mu is base intensity. This is generally a constant vector but can be extended to stochastic processes. Currently, piecewise constant mu is also possible. mu is left continuous.

\alpha is a constant matrix which represents impacts on intensities after events. It is represented by slot alpha.

g is for non-constant parts of the impact. It may depend on any information generated by N, \lambda, z and so on. It is represented by slot impact.

\beta is a constant matrix for exponential decay rates. It is represented by slot beta.

z is mark and represented by slot rmark.

mu, alpha and beta are required slots for every exponential Hawkes model. rmark and impact are additional slots.

Slots

mu

Numeric value or matrix or function, if numeric, automatically converted to matrix.

alpha

Numeric value or matrix or function, if numeric, automatically converted to matrix, exciting term.

beta

Numeric value or matrix or function, if numeric, automatically converted to matrix, exponential decay.

eta

Numeric value or matrix or function, if numeric, automatically converted to matrix, impact by additional mark.

dimens

Dimension of the model.

rmark

A function that generates mark for counting process, for simulation.

dmark

A density function for mark, for estimation.

impact

A function that describe the after impact of mark to lambda whose first argument is always param.

type_col_map

Mapping between type and column number of kernel used for multi-kernel model.

Examples

MU <- matrix(c(0.2), nrow = 2)
ALPHA <- matrix(c(0.75, 0.92, 0.92, 0.75), nrow = 2, byrow=TRUE)
BETA <- matrix(c(2.25, 2.25, 2.25, 2.25), nrow = 2, byrow=TRUE)
mhspec2 <- new("hspec", mu=MU, alpha=ALPHA, beta=BETA)
mhspec2

Compute Hawkes volatility

Description

This function computes Hawkes volatility. Only works for bi-variate Hawkes process.

Usage

hvol(
  object,
  horizon = 1,
  inter_arrival = NULL,
  type = NULL,
  mark = NULL,
  dependence = FALSE,
  lambda_component0 = NULL,
  ...
)

## S4 method for signature 'hspec'
hvol(
  object,
  horizon = 1,
  inter_arrival = NULL,
  type = NULL,
  mark = NULL,
  dependence = FALSE,
  lambda_component0 = NULL,
  ...
)

Arguments

Infer lambda process with given Hawkes model and realized path

Description

This method compute the inferred lambda process and returns it as hreal form. If we have realized path of Hawkes process and its parameter value, then we can compute the inferred lambda processes. Similarly with other method such as hfit, the input arguments are inter_arrival, type, mark, or equivalently, N and Nc.

Usage

infer_lambda(
  object,
  inter_arrival = NULL,
  type = NULL,
  mark = NULL,
  N = NULL,
  Nc = NULL,
  lambda_component0 = NULL,
  N0 = NULL,
  ...
)

## S4 method for signature 'hspec'
infer_lambda(
  object,
  inter_arrival = NULL,
  type = NULL,
  mark = NULL,
  N = NULL,
  Nc = NULL,
  lambda_component0 = NULL,
  N0 = NULL,
  ...
)

Arguments

Value

hreal S3-object, with inferred intensity.

Examples

mu <- c(0.1, 0.1)
alpha <- matrix(c(0.2, 0.1, 0.1, 0.2), nrow=2, byrow=TRUE)
beta <- matrix(c(0.9, 0.9, 0.9, 0.9), nrow=2, byrow=TRUE)
h <- new("hspec", mu=mu, alpha=alpha, beta=beta)
res <- hsim(h, size=100)
summary(res)
res2 <- infer_lambda(h, res$inter_arrival, res$type)
summary(res2)

Compute the log-likelihood function

Description

The log-likelihood of the ground process of the Hawkes model. (The log-likelihood for mark (jump) distribution is not provided.)

Usage

## S4 method for signature 'hspec'
logLik(
  object,
  inter_arrival,
  type = NULL,
  mark = NULL,
  N = NULL,
  Nc = NULL,
  N0 = NULL,
  lambda_component0 = NULL,
  ...
)

Arguments

See Also

hspec-class, hfit,hspec-method

Compute residual process

Description

Using random time change, this function compute the residual process, which is the inter-arrival time of a standard Poisson process. Therefore, the return values should follow the exponential distribution with rate 1, if model and rambda are correctly specified.

Usage

residual_process(
  component,
  inter_arrival,
  type,
  rambda_component,
  mu,
  beta,
  dimens = NULL,
  mark = NULL,
  N = NULL,
  Nc = NULL,
  lambda_component0 = NULL,
  N0 = NULL,
  ...
)

Arguments

Examples

mu <- c(0.1, 0.1)
alpha <- matrix(c(0.2, 0.1, 0.1, 0.2), nrow=2, byrow=TRUE)
beta <- matrix(c(0.9, 0.9, 0.9, 0.9), nrow=2, byrow=TRUE)
h <- new("hspec", mu=mu, alpha=alpha, beta=beta)
res <- hsim(h, size=1000)
rp <- residual_process(component = 1, res$inter_arrival, res$type, res$rambda_component, mu, beta)