Title: | Renewal Hawkes Process |
---|---|
Description: | The renewal Hawkes (RHawkes) process (Wheatley, Filimonov, and Sornette, 2016 <doi:10.1016/j.csda.2015.08.007>) is an extension to the classical Hawkes self-exciting point process widely used in the modelling of clustered event sequence data. This package provides functions to simulate the RHawkes process with a given immigrant hazard rate function and offspring birth time density function, to compute the exact likelihood of a RHawkes process using the recursive algorithm proposed by Chen and Stindl (2018) <doi:10.1080/10618600.2017.1341324>, to compute the Rosenblatt residuals for goodness-of-fit assessment, and to predict future event times based on observed event times up to a given time. A function implementing the linear time RHawkes process likelihood approximation algorithm proposed in Stindl and Chen (2021) <doi:10.1007/s11222-021-10002-0> is also included. |
Authors: | Feng Chen [aut, cre] , Tom Stindl [ctb] |
Maintainer: | Feng Chen <[email protected]> |
License: | GPL (>= 2) |
Version: | 1.0 |
Built: | 2025-01-10 05:19:08 UTC |
Source: | https://github.com/cran/RHawkes |
The renewal Hawkes (RHawkes) process (Wheatley, Filimonov, and Sornette, 2016 <doi:10.1016/j.csda.2015.08.007>) is an extension to the classical Hawkes self-exciting point process widely used in the modelling of clustered event sequence data. This package provides functions to simulate the RHawkes process with a given immigrant hazard rate function and offspring birth time density function, to compute the exact likelihood of a RHawkes process using the recursive algorithm proposed by Chen and Stindl (2018) <doi:10.1080/10618600.2017.1341324>, to compute the Rosenblatt residuals for goodness-of-fit assessment, and to predict future event times based on observed event times up to a given time. A function implementing the linear time RHawkes process likelihood approximation algorithm proposed in Stindl and Chen (2021) <doi:10.1007/s11222-021-10002-0> is also included.
The DESCRIPTION file:
Package: | RHawkes |
Type: | Package |
Title: | Renewal Hawkes Process |
Version: | 1.0 |
Date: | 2022-5-4 |
Authors@R: | c(person(given="Feng", family="Chen", role = c("aut", "cre"), email = "[email protected]", comment = c(ORCID="0000-0002-9646-3338") ), person(given="Tom", family="Stindl", role = "ctb", email="[email protected]", comment=c(ORCID="0000-0001-8627-9337") ) ) |
Description: | The renewal Hawkes (RHawkes) process (Wheatley, Filimonov, and Sornette, 2016 <doi:10.1016/j.csda.2015.08.007>) is an extension to the classical Hawkes self-exciting point process widely used in the modelling of clustered event sequence data. This package provides functions to simulate the RHawkes process with a given immigrant hazard rate function and offspring birth time density function, to compute the exact likelihood of a RHawkes process using the recursive algorithm proposed by Chen and Stindl (2018) <doi:10.1080/10618600.2017.1341324>, to compute the Rosenblatt residuals for goodness-of-fit assessment, and to predict future event times based on observed event times up to a given time. A function implementing the linear time RHawkes process likelihood approximation algorithm proposed in Stindl and Chen (2021) <doi:10.1007/s11222-021-10002-0> is also included. |
License: | GPL (>= 2) |
Depends: | R (>= 2.10), IHSEP |
NeedsCompilation: | no |
Packaged: | 2022-05-04 13:17:53 UTC; z3243864 |
Author: | Feng Chen [aut, cre] (<https://orcid.org/0000-0002-9646-3338>), Tom Stindl [ctb] (<https://orcid.org/0000-0001-8627-9337>) |
Maintainer: | Feng Chen <[email protected]> |
Date/Publication: | 2022-05-05 14:40:09 UTC |
Repository: | https://fchenunswms.r-universe.dev |
RemoteUrl: | https://github.com/cran/RHawkes |
RemoteRef: | HEAD |
RemoteSha: | 9a033bc7a3c829df6aa494c2bf79501c7152df3d |
Index of help topics:
EM1partial Partial EM algorithm for the RHawkes process, version 1 EM2partial Partial EM algorithm for the RHawkes process, version 2 RHawkes-package Renewal Hawkes Process damllRH Dynamically approxomated minus loglikelihood of a RHawkes model mllRH Minus loglikelihood of a RHawkes model mllRH1 Minus loglikelihood of a RHawkes model with parent probabilities mllRH2 Minus loglikelihood of a RHawkes model with Rosenblatt residuals pred.den RHawkes predictive density function pred.haz RHawkes predictive hazard function quake An RHawkes earthquake data set sim.pred Simulate a fitted RHawkes process model sim.pred1 Simulate a fitted RHawkes process model for prediction purposes simRHawkes Simulate a renewal Hawkes (RHawkes) process simRHawkes1 Simulate a renewal Hawkes (RHawkes) process tms mid-price change times of the AUD/USD exchange rate
Feng Chen [aut, cre] (<https://orcid.org/0000-0002-9646-3338>), Tom Stindl [ctb] (<https://orcid.org/0000-0001-8627-9337>)
Maintainer: Feng Chen <[email protected]>
Calculates an apprximation to the minus loglikelihood of a
RHawkes model with given immigration hazard function ,
offspring birth time density function
and branching ratio
relative to event times
tms
on interval .
damllRH(tms, cens, par, q=0.999, qe=0.999, h.fn = function(x, p) dexp(x, rate = 1 / p), mu.fn = function(x, p) { exp(dweibull(x, shape = p[1], scale = p[2], log = TRUE) - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE)) }, H.fn = function(x, p) pexp(x, rate = 1 / p), Mu.fn = function(x, p) { -pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE) }, keepB=FALSE, H.inv=function(x,p)qexp(x,rate=1/p) )
damllRH(tms, cens, par, q=0.999, qe=0.999, h.fn = function(x, p) dexp(x, rate = 1 / p), mu.fn = function(x, p) { exp(dweibull(x, shape = p[1], scale = p[2], log = TRUE) - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE)) }, H.fn = function(x, p) pexp(x, rate = 1 / p), Mu.fn = function(x, p) { -pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE) }, keepB=FALSE, H.inv=function(x,p)qexp(x,rate=1/p) )
tms |
A numeric vector, with values sorted in ascending order. Event times to fit the RHawkes point process model. |
cens |
A numericl scalar. The censoring time. |
par |
A numeric vector containing the parameters of the model,
in order of the immigration parameters, in |
q |
A numeric scalar in (0,1] and close to 1, which controls how far we look back when truncating the distribution of the most recent immigrant. |
qe |
A numeric scalar in (0,1] and close to 1, which controls how to truncation is used in the offspring birth time distribution. |
h.fn |
A (vectorized) function. The offspring birth time density function. |
mu.fn |
A (vectorized) function. The immigrant waiting time hazard function. |
H.fn |
A (vectorized) function. Its value at |
Mu.fn |
A (vectorized) function. Its value at |
keepB |
A boolean scalar, indicating whether the looking back
values |
H.inv |
A (vectorized) function, giving the inverse function of the integral of the excitation. |
A scalar giving the value of the (approximate) negative
log-likelihood, when keepB
is FALSE (the default); A list with
components mll
, whhich contains the value of the negative
log-likelihood, Bs
, which gives the look-back order of the
truncation of the distribution of the last immigrant, and Bes
,
which gives the look-forward order in determining how far into the
future the excitation effect is allowed to last.
Feng Chen <[email protected]>
## Not run: ## earthquake times over 96 years data(quake); tms <- sort(quake$time); # add some random noise to the simultaneous occurring event times tms[213:214] <- tms[213:214] + sort(c(runif(1, -1, 1)/(24*60), runif(1, -1, 1)/(24*60))) ## calculate the minus loglikelihood of an RHawkes with some parameters ## the default hazard function and density functions are Weibull and ## exponential respectively mllRH(tms, cens = 96*365.25 , par = c(0.5, 20, 1000, 0.5)) damllRH(tms, cens = 96*365.25 , par = c(0.5, 20, 1000, 0.5),q=1,qe=1) ## calculate the MLE for the parameter assuming known parametric forms ## of the immigrant hazard function and offspring density functions. system.time(est <- optim(c(0.5, 20, 1000, 0.5), mllRH, tms = tms, cens = 96*365.25, mu.fn=function(x,p)p[1]/p[2]*(x/p[2])^(p[1]-1), Mu.fn=function(x,p)(x/p[2])^p[1], control = list(maxit = 5000, trace = TRUE), hessian = TRUE) ) system.time(est1 <- optim(c(0.5, 20, 1000, 0.5), function(p){ if(any(p<0)||p[4]<0||p[4]>=1) return(Inf); damllRH(tms = tms, cens = 96*365.25, mu.fn=function(x,p)p[1]/p[2]*(x/p[2])^(p[1]-1), Mu.fn=function(x,p)(x/p[2])^p[1], par=p,q=0.999999,qe=0.999999) }, control = list(maxit = 5000, trace = TRUE), hessian = TRUE) ) ## point estimate by MLE est$par est1$par ## standard error estimates: diag(solve(est$hessian))^0.5 diag(solve(est1$hessian))^0.5 ## End(Not run)
## Not run: ## earthquake times over 96 years data(quake); tms <- sort(quake$time); # add some random noise to the simultaneous occurring event times tms[213:214] <- tms[213:214] + sort(c(runif(1, -1, 1)/(24*60), runif(1, -1, 1)/(24*60))) ## calculate the minus loglikelihood of an RHawkes with some parameters ## the default hazard function and density functions are Weibull and ## exponential respectively mllRH(tms, cens = 96*365.25 , par = c(0.5, 20, 1000, 0.5)) damllRH(tms, cens = 96*365.25 , par = c(0.5, 20, 1000, 0.5),q=1,qe=1) ## calculate the MLE for the parameter assuming known parametric forms ## of the immigrant hazard function and offspring density functions. system.time(est <- optim(c(0.5, 20, 1000, 0.5), mllRH, tms = tms, cens = 96*365.25, mu.fn=function(x,p)p[1]/p[2]*(x/p[2])^(p[1]-1), Mu.fn=function(x,p)(x/p[2])^p[1], control = list(maxit = 5000, trace = TRUE), hessian = TRUE) ) system.time(est1 <- optim(c(0.5, 20, 1000, 0.5), function(p){ if(any(p<0)||p[4]<0||p[4]>=1) return(Inf); damllRH(tms = tms, cens = 96*365.25, mu.fn=function(x,p)p[1]/p[2]*(x/p[2])^(p[1]-1), Mu.fn=function(x,p)(x/p[2])^p[1], par=p,q=0.999999,qe=0.999999) }, control = list(maxit = 5000, trace = TRUE), hessian = TRUE) ) ## point estimate by MLE est$par est1$par ## standard error estimates: diag(solve(est$hessian))^0.5 diag(solve(est1$hessian))^0.5 ## End(Not run)
Calculates the RHawkes model parameters via a partial Expectation-Maximization (EM1) algorithm of Wheatley, Filimonov and Sornette (2016).
EM1partial(tms, cens, pars, maxiter = 1000, tol = 1e-8, h.fn = function(x, p) dexp(x, rate = 1 / p), mu.fn = function(x, p){ exp(dweibull(x, shape = p[1], scale = p[2], log = TRUE) - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE)) }, H.fn = function(x, p) pexp(x, rate = 1 / p), logg.fn = function(x, p){ dweibull(x, shape = p[1], scale = p[2], log = TRUE) - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE) - (x / p[2])^p[1]}, Mu.fn = function(x, p){ - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE) })
EM1partial(tms, cens, pars, maxiter = 1000, tol = 1e-8, h.fn = function(x, p) dexp(x, rate = 1 / p), mu.fn = function(x, p){ exp(dweibull(x, shape = p[1], scale = p[2], log = TRUE) - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE)) }, H.fn = function(x, p) pexp(x, rate = 1 / p), logg.fn = function(x, p){ dweibull(x, shape = p[1], scale = p[2], log = TRUE) - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE) - (x / p[2])^p[1]}, Mu.fn = function(x, p){ - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE) })
tms |
A numeric vector, with values sorted in ascending order. Event times to fit the RHawkes point process model. |
cens |
A scalar. The censoring time. |
pars |
A numeric vector containing the parameters of the model, in order of the
immigration parameters |
maxiter |
The maximum number of iterations to perform. |
tol |
The algorithm stops when the difference between the previous iteration and
current iteration parameters sum is less than |
h.fn |
A (vectorized) function. The offspring density function. |
mu.fn |
A (vectorized) function. The immigration hazard function. |
H.fn |
A (vectorized) function. Its value at |
Mu.fn |
A (vectorized) function. Its value at |
logg.fn |
A (vectorized) function. The log of the immigrant distribution function. |
iterations |
The number of iterations until convergence |
diff |
The absolute sum of the difference between the final two parameter estimates |
pars |
The parameter estimates from the EM algorithm |
Feng Chen <[email protected]> Tom Stindl <[email protected]>
## Not run: ## simulated data tms <- sort(runif(100,0,100)) ## the slower version of the EM algorithms on simulated data with default ## immigrant hazard function ## and offspring density system.time( est1 <- EM1partial(tms, 101, c(2,1,0.5,1)) ) ## End(Not run)
## Not run: ## simulated data tms <- sort(runif(100,0,100)) ## the slower version of the EM algorithms on simulated data with default ## immigrant hazard function ## and offspring density system.time( est1 <- EM1partial(tms, 101, c(2,1,0.5,1)) ) ## End(Not run)
Calculates the RHawkes model parameters via a partial Expectation-Maximization (EM2) algorithm of Wheatley, Filimonov and Sornette (2016).
EM2partial(tms, cens, pars, maxiter = 1000, tol = 1e-8, h.fn = function(x, p) dexp(x, rate = 1 / p), mu.fn = function(x, p){ exp(dweibull(x, shape = p[1], scale = p[2], log = TRUE) - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE)) }, H.fn = function(x, p) pexp(x, rate = 1 / p), logg.fn = function(x, p){ dweibull(x, shape = p[1], scale = p[2], log = TRUE) - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE) - (x / p[2])^p[1]}, Mu.fn = function(x, p){ - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE) })
EM2partial(tms, cens, pars, maxiter = 1000, tol = 1e-8, h.fn = function(x, p) dexp(x, rate = 1 / p), mu.fn = function(x, p){ exp(dweibull(x, shape = p[1], scale = p[2], log = TRUE) - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE)) }, H.fn = function(x, p) pexp(x, rate = 1 / p), logg.fn = function(x, p){ dweibull(x, shape = p[1], scale = p[2], log = TRUE) - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE) - (x / p[2])^p[1]}, Mu.fn = function(x, p){ - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE) })
tms |
A numeric vector, with values sorted in ascending order. Event times to fit the RHawkes point process model. |
cens |
A scalar. The censoring time. |
pars |
A numeric vector containing the parameters of the model, in order of the
immigration parameters |
maxiter |
The maximum number of iterations to perform. |
tol |
The algorithm stops when the difference between the previous iteration and
current iteration parameters sum is less than |
h.fn |
A (vectorized) function. The offspring density function. |
mu.fn |
A (vectorized) function. The immigration hazard function. |
H.fn |
A (vectorized) function. Its value at |
Mu.fn |
A (vectorized) function. Its value at |
logg.fn |
A (vectorized) function. The log of the immigrant distribution function. |
iterations |
The number of iterations until convergence |
diff |
The absolute sum of the difference between the final two parameter estimates |
pars |
The parameter estimates from the EM algorithm |
Feng Chen <[email protected]> Tom Stindl <[email protected]>
## Not run: ## simulated data tms <- sort(runif(100,0,100)) ## the quicker version on simulated data with default immigrant hazard function ## and offspring density system.time( est2 <- EM2partial(tms, 101, c(2,1,0.5,1)) ) ## End(Not run)
## Not run: ## simulated data tms <- sort(runif(100,0,100)) ## the quicker version on simulated data with default immigrant hazard function ## and offspring density system.time( est2 <- EM2partial(tms, 101, c(2,1,0.5,1)) ) ## End(Not run)
Calculates the minus loglikelihood of a RHawkes model with given
immigration hazard function , offspring density function
and branching ratio
for event times
tms
on interval .
mllRH(tms, cens, par, h.fn = function(x, p) dexp(x, rate = 1 / p), mu.fn = function(x, p) { exp(dweibull(x, shape = p[1], scale = p[2], log = TRUE) - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE)) }, H.fn = function(x, p) pexp(x, rate = 1 / p), Mu.fn = function(x, p) { -pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE) })
mllRH(tms, cens, par, h.fn = function(x, p) dexp(x, rate = 1 / p), mu.fn = function(x, p) { exp(dweibull(x, shape = p[1], scale = p[2], log = TRUE) - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE)) }, H.fn = function(x, p) pexp(x, rate = 1 / p), Mu.fn = function(x, p) { -pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE) })
tms |
A numeric vector, with values sorted in ascending order. Event times to fit the RHawkes point process model. |
cens |
A scalar. The censoring time. |
par |
A numeric vector containing the parameters of the model, in order of the
immigration parameters |
h.fn |
A (vectorized) function. The offspring density function. |
mu.fn |
A (vectorized) function. The immigration hazard function. |
H.fn |
A (vectorized) function. Its value at |
Mu.fn |
A (vectorized) function. Its value at |
The value of the negative log-likelihood.
Feng Chen <[email protected]> Tom Stindl <[email protected]>
## Not run: ## earthquake times over 96 years data(quake); tms <- sort(quake$time); # add some random noise to the simultaneous occurring event times tms[213:214] <- tms[213:214] + sort(c(runif(1, -1, 1)/(24*60), runif(1, -1, 1)/(24*60))) ## calculate the minus loglikelihood of an RHawkes with some parameters ## the default hazard function and density functions are Weibull and ## exponential respectively mllRH(tms, cens = 96*365.25 , par = c(0.5, 20, 1000, 0.5)) ## calculate the MLE for the parameter assuming known parametric forms ## of the immigrant hazard function and offspring density functions. system.time(est <- optim(c(0.5, 20, 1000, 0.5), mllRH, tms = tms, cens = 96*365.25, control = list(maxit = 5000, trace = TRUE), hessian = TRUE) ) ## point estimate by MLE est$par ## standard error estimates: diag(solve(est$hessian))^0.5 ## End(Not run)
## Not run: ## earthquake times over 96 years data(quake); tms <- sort(quake$time); # add some random noise to the simultaneous occurring event times tms[213:214] <- tms[213:214] + sort(c(runif(1, -1, 1)/(24*60), runif(1, -1, 1)/(24*60))) ## calculate the minus loglikelihood of an RHawkes with some parameters ## the default hazard function and density functions are Weibull and ## exponential respectively mllRH(tms, cens = 96*365.25 , par = c(0.5, 20, 1000, 0.5)) ## calculate the MLE for the parameter assuming known parametric forms ## of the immigrant hazard function and offspring density functions. system.time(est <- optim(c(0.5, 20, 1000, 0.5), mllRH, tms = tms, cens = 96*365.25, control = list(maxit = 5000, trace = TRUE), hessian = TRUE) ) ## point estimate by MLE est$par ## standard error estimates: diag(solve(est$hessian))^0.5 ## End(Not run)
Calculates the minus loglikelihood of a RHawkes model with given
immigration hazard function , offspring density function
and branching ratio
for event times
tms
on interval . The same as
mllRH
although this
version also returns the parent probabilities.
mllRH1(tms, cens, par, h.fn = function(x, p) dexp(x, rate = 1/p), mu.fn = function(x, p) { exp(dweibull(x, shape = p[1], scale = p[2], log = TRUE) - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE)) }, H.fn = function(x, p) pexp(x, rate = 1/p), Mu.fn = function(x, p) { -pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE) })
mllRH1(tms, cens, par, h.fn = function(x, p) dexp(x, rate = 1/p), mu.fn = function(x, p) { exp(dweibull(x, shape = p[1], scale = p[2], log = TRUE) - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE)) }, H.fn = function(x, p) pexp(x, rate = 1/p), Mu.fn = function(x, p) { -pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE) })
tms |
A numeric vector, with values sorted in ascending order. Event times to fit the RHawkes point process model. |
cens |
A scalar. The censoring time. |
par |
A numeric vector containing the parameters of the model, in order of the
immigration parameters |
h.fn |
A (vectorized) function. The offspring density function. |
mu.fn |
A (vectorized) function. The immigration hazard function. |
H.fn |
A (vectorized) function. Its value at |
Mu.fn |
A (vectorized) function. Its value at |
mll |
minus log-likelihood |
log.p |
parent probabilities |
n |
number of events |
Feng Chen <[email protected]> Tom Stindl <[email protected]>
mllRH
tmp <- mllRH1(sort(runif(1000,0,1000)), 1001, c(2,1,0.5,1)) for(i in 1:tmp$n) cat(exp(tmp$log.p[i*(i - 1)/2 + 1:i]), "\n")
tmp <- mllRH1(sort(runif(1000,0,1000)), 1001, c(2,1,0.5,1)) for(i in 1:tmp$n) cat(exp(tmp$log.p[i*(i - 1)/2 + 1:i]), "\n")
Calculates the minus loglikelihood of a RHawkes model with given
immigration hazard function , offspring density function
and branching ratio
for event times
tms
on interval . The same as
mllRH
although this
version also returns the Rosenblatt residuals.
mllRH2(tms, cens, par, h.fn = function(x, p) dexp(x, rate = 1/p), mu.fn = function(x, p) { exp(dweibull(x, shape = p[1], scale = p[2], log = TRUE) - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE))}, H.fn = function(x, p) pexp(x, rate = 1/p), Mu.fn = function(x, p) { -pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE) })
mllRH2(tms, cens, par, h.fn = function(x, p) dexp(x, rate = 1/p), mu.fn = function(x, p) { exp(dweibull(x, shape = p[1], scale = p[2], log = TRUE) - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE))}, H.fn = function(x, p) pexp(x, rate = 1/p), Mu.fn = function(x, p) { -pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE) })
tms |
A numeric vector, with values sorted in ascending order. Event times to fit the RHawkes point process model. |
cens |
A scalar. The censoring time. |
par |
A numeric vector containing the parameters of the model, in order of the
immigration parameters |
h.fn |
A (vectorized) function. The offspring density function. |
mu.fn |
A (vectorized) function. The immigration hazard function. |
H.fn |
A (vectorized) function. Its value at |
Mu.fn |
A (vectorized) function. Its value at |
Calculate the RHawkes point process Rosenblatt residuals
mll |
minus log-likelihood |
U |
Rosenblatt residual of observed event time |
n |
number of events |
Feng Chen <[email protected]> Tom Stindl <[email protected]>
mllRH
## Not run: tmp <- mllRH2(sort(runif(1000,0,1000)),1001,c(2,1,0.5,1)) par(mfrow=c(1,2)) qqunif<-function(dat,...){ dat<-sort(as.numeric(dat)); n<-length(dat); pvec<-ppoints(n); plot(pvec,dat,xlab="Theoretical Quantiles", ylab="Sample Quantiles",main="Uniform Q-Q Plot",...) } qqunif(tmp$U) acf(tmp$U) ks.test(tmp$U,"punif") ## End(Not run)
## Not run: tmp <- mllRH2(sort(runif(1000,0,1000)),1001,c(2,1,0.5,1)) par(mfrow=c(1,2)) qqunif<-function(dat,...){ dat<-sort(as.numeric(dat)); n<-length(dat); pvec<-ppoints(n); plot(pvec,dat,xlab="Theoretical Quantiles", ylab="Sample Quantiles",main="Uniform Q-Q Plot",...) } qqunif(tmp$U) acf(tmp$U) ks.test(tmp$U,"punif") ## End(Not run)
Calculates the predictive density of the next observed event time after the
censoring time cens
based on observations over the interval
[0,cens
].
pred.den(x, tms, cens, par, h.fn = function(x, p) dexp(x, rate = 1 / p), mu.fn = function(x, p) { exp(dweibull(x, shape = p[1], scale = p[2], log = TRUE) - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE))}, H.fn = function(x, p) pexp(x, rate = 1 / p), Mu.fn = function(x, p) { -pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE) })
pred.den(x, tms, cens, par, h.fn = function(x, p) dexp(x, rate = 1 / p), mu.fn = function(x, p) { exp(dweibull(x, shape = p[1], scale = p[2], log = TRUE) - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE))}, H.fn = function(x, p) pexp(x, rate = 1 / p), Mu.fn = function(x, p) { -pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE) })
x |
A scalar. The amount of time after the censoring time |
tms |
A numeric vector, with values sorted in ascending order. The event times to fit the RHawkes point process model. |
cens |
A scalar. The censoring time. |
par |
A numeric vector. Contains the parameters of the model, in order of the
immigration parameters |
h.fn |
A (vectorized) function. The offspring density function. |
mu.fn |
A (vectorized) function. The immigration hazard function. |
H.fn |
A (vectorized) function. Its value at |
Mu.fn |
A (vectorized) function. Its value at |
The predictive density of the next event evaluated at x
.
Feng Chen <[email protected]> Tom Stindl <[email protected]>
data(quake); tms <- sort(quake$time); # add some random noise to the identical event times tms[213:214] <- tms[213:214] + sort(c(runif(1, -1, 1)/(24*60), runif(1, -1, 1)/(24*60))) curve(pred.den(x, tms = tms, cens = 35064, par= c(0.314, 22.2, 1266, 0.512)) ,0 ,2000, col = 2, lty = 2)
data(quake); tms <- sort(quake$time); # add some random noise to the identical event times tms[213:214] <- tms[213:214] + sort(c(runif(1, -1, 1)/(24*60), runif(1, -1, 1)/(24*60))) curve(pred.den(x, tms = tms, cens = 35064, par= c(0.314, 22.2, 1266, 0.512)) ,0 ,2000, col = 2, lty = 2)
Calculates the predictive hazard function of the next observed event time after
the censoring time cens
based on observations over the interval
[0,cens
].
pred.haz(x, tms, cens, par, h.fn = function(x, p) dexp(x, rate = 1 / p), mu.fn = function(x, p) { exp(dweibull(x, shape = p[1], scale = p[2], log = TRUE) - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE)) }, H.fn = function(x, p) pexp(x, rate = 1/p), Mu.fn = function(x, p) { - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE) })
pred.haz(x, tms, cens, par, h.fn = function(x, p) dexp(x, rate = 1 / p), mu.fn = function(x, p) { exp(dweibull(x, shape = p[1], scale = p[2], log = TRUE) - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE)) }, H.fn = function(x, p) pexp(x, rate = 1/p), Mu.fn = function(x, p) { - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE) })
x |
A scalar. The amount of time after the censoring time |
tms |
A numeric vector, with values sorted in ascending order. The event times to fit the RHawkes point process model. |
cens |
A scalar. The censoring time. |
par |
A numeric vector. Contains the parameters of the model, in order of the
immigration parameters |
h.fn |
A (vectorized) function. The offspring density function. |
mu.fn |
A (vectorized) function. The immigration hazard function. |
H.fn |
A (vectorized) function. Its value at |
Mu.fn |
A (vectorized) function. Its value at |
The predictive hazard rate of the next event evaluated at x
.
Feng Chen <[email protected]> Tom Stindl <[email protected]>
data(quake); tms <- sort(quake$time); # add some random noise to the identical event times tms[213:214] <- tms[213:214] + sort(c(runif(1, -1, 1)/(24*60), runif(1, -1, 1)/(24*60))) curve(pred.haz(x, tms = tms, cens = 35064, par= c(0.314, 22.2, 1266, 0.512)) ,0 ,2000, col = 2, lty = 2)
data(quake); tms <- sort(quake$time); # add some random noise to the identical event times tms[213:214] <- tms[213:214] + sort(c(runif(1, -1, 1)/(24*60), runif(1, -1, 1)/(24*60))) curve(pred.haz(x, tms = tms, cens = 35064, par= c(0.314, 22.2, 1266, 0.512)) ,0 ,2000, col = 2, lty = 2)
An earthquake data set containing the earthquake occurrence times near the Japan region previously examined by Ogata (1998).
data(quake)
data(quake)
The format is a vector of the arrival/birth times of earthquakes.
Times of arrivals of earthquake occurrences in a vector in ascending order.
Simulated by a call to the function simHawkes1
.
## Not run: data(quake) ## number of earthquake occurrences nrow(quake) ## End(Not run)
## Not run: data(quake) ## number of earthquake occurrences nrow(quake) ## End(Not run)
Simulate a fitted RHawkes process model after the censoring time cens
to a future time point cens.tilde
.
sim.pred(tms, re.dist = rweibull, par, par.redist = list(shape = par[1], scale = par[2]), h.fn = function(x, p) dexp(x, rate = 1 / p), p.ofd = par[3], branching.ratio = par[4], cens, cens.tilde = cens * 1.5, mu.fn = function(x, p) { exp(dweibull(x, shape = p[1], scale = p[2], log = TRUE) - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE)) })
sim.pred(tms, re.dist = rweibull, par, par.redist = list(shape = par[1], scale = par[2]), h.fn = function(x, p) dexp(x, rate = 1 / p), p.ofd = par[3], branching.ratio = par[4], cens, cens.tilde = cens * 1.5, mu.fn = function(x, p) { exp(dweibull(x, shape = p[1], scale = p[2], log = TRUE) - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE)) })
tms |
A numeric vector, with values sorted in ascending order. Event times to fit the RHawkes point process model. |
re.dist |
A (vectorized) function. The immigrant renewal distribution function. |
par |
A numeric vector, giving the parameters of the model with the
immigration parameters |
par.redist |
A numeric vector. The parameters of the immigrant renewal distribution. |
h.fn |
A (vectorized) function. The offspring density function. |
p.ofd |
A (named) list. The parameters of the offspring density. |
branching.ratio |
A scalar. The branching ratio parameter. |
cens |
A scalar. The censoring time. |
cens.tilde |
A scalar. The future time that the simulation run until. |
mu.fn |
A (vectorized) function. The immigration hazard function. |
A numeric vector that contains the simulated event times from censoring time
cens
up until cens.tilde
Feng Chen <[email protected]> Tom Stindl <[email protected]>
N <- 5; i <- 0; data(quake); tms <- sort(quake$time); # add some random noise the simultaneous occurring event times tms[213:214] <- tms[213:214] + sort(c(runif(1, -1, 1)/(24*60), runif(1, -1, 1)/(24*60))) # simulate future event time based on MLE fitted RHawkes model times <- replicate(N, {cat(i<<-i+1,'\n'); sim.pred(tms = tms, par = c(0.314, 22.2, 1266, 0.512), cens=35063) }) plot(NA,NA,xlim=c(0,35063*1.5),ylim=c(0,max(lengths(times))+nrow(quake)), xlab="time",ylab="Sample path") lines(c(0,quake$time),0:nrow(quake),type="s") for(i in 1:N) lines(c(tail(quake$time,1),times[[i]]),nrow(quake)+0:length(times[[i]]), type="s",lty=2)
N <- 5; i <- 0; data(quake); tms <- sort(quake$time); # add some random noise the simultaneous occurring event times tms[213:214] <- tms[213:214] + sort(c(runif(1, -1, 1)/(24*60), runif(1, -1, 1)/(24*60))) # simulate future event time based on MLE fitted RHawkes model times <- replicate(N, {cat(i<<-i+1,'\n'); sim.pred(tms = tms, par = c(0.314, 22.2, 1266, 0.512), cens=35063) }) plot(NA,NA,xlim=c(0,35063*1.5),ylim=c(0,max(lengths(times))+nrow(quake)), xlab="time",ylab="Sample path") lines(c(0,quake$time),0:nrow(quake),type="s") for(i in 1:N) lines(c(tail(quake$time,1),times[[i]]),nrow(quake)+0:length(times[[i]]), type="s",lty=2)
Simulate a fitted RHawkes process model from the censoring time cens
to a future time point cens.tilde
, conditional on the observed
event times until the censoring time.
sim.pred1(tms, par, re.dist = rweibull, par.redist = list(shape = par[1], scale = par[2]), of.dis="exp", par.ofdis = list(rate=par[3]), branching.ratio = par[4], cens=tail(tms,1)+mean(diff(tms))/2, cens.tilde = cens * 1.5, mu.fn = function(x, p) { exp(dweibull(x, shape = p[1], scale = p[2], log = TRUE) - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE)) })
sim.pred1(tms, par, re.dist = rweibull, par.redist = list(shape = par[1], scale = par[2]), of.dis="exp", par.ofdis = list(rate=par[3]), branching.ratio = par[4], cens=tail(tms,1)+mean(diff(tms))/2, cens.tilde = cens * 1.5, mu.fn = function(x, p) { exp(dweibull(x, shape = p[1], scale = p[2], log = TRUE) - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE)) })
tms |
A numeric vector, with values sorted in ascending order. Event times to fit the RHawkes point process model. |
par |
A numeric vector, giving the parameters of the model with the
immigration parameters |
re.dist |
A (vectorized) function. The function to simulate from the immigrant waiting times distribution. |
par.redist |
A (named) list, giving the parameters of the immigrant waiting time distribution. |
of.dis |
A character string, for the name of the offspring birth time distribution. |
par.ofdis |
A (named) list, giving the parameters of the offspring birth time distribution. |
branching.ratio |
A scalar in [0,1), the branching ratio parameter. |
cens |
A scalar. The censoring time. |
cens.tilde |
A scalar. The future time to run the simulation to. |
mu.fn |
A (vectorized) function. The immigration hazard function. |
A numeric vector that contains the simulated event times from censoring time
cens
up until cens.tilde
Feng Chen <[email protected]> Tom Stindl <[email protected]>
N <- 5; i <- 0; data(quake); tms <- sort(quake$time); # add some random noise the simultaneous occurring event times tms[213:214] <- tms[213:214] + sort(c(runif(1, -1, 1)/(24*60), runif(1, -1, 1)/(24*60))) # simulate future event time based on MLE fitted RHawkes model times <- replicate(N, {cat(i<<-i+1,'\n'); sim.pred1(tms = tms, par = c(0.314, 22.2, 1266, 0.512), cens=35063) }) plot(NA,NA,xlim=c(0,35063*1.5),ylim=c(0,max(lengths(times))+nrow(quake)), xlab="time",ylab="Sample path") lines(c(0,quake$time),0:nrow(quake),type="s") for(i in 1:N) lines(c(tail(quake$time,1),times[[i]]),nrow(quake)+0:length(times[[i]]), type="s",lty=2)
N <- 5; i <- 0; data(quake); tms <- sort(quake$time); # add some random noise the simultaneous occurring event times tms[213:214] <- tms[213:214] + sort(c(runif(1, -1, 1)/(24*60), runif(1, -1, 1)/(24*60))) # simulate future event time based on MLE fitted RHawkes model times <- replicate(N, {cat(i<<-i+1,'\n'); sim.pred1(tms = tms, par = c(0.314, 22.2, 1266, 0.512), cens=35063) }) plot(NA,NA,xlim=c(0,35063*1.5),ylim=c(0,max(lengths(times))+nrow(quake)), xlab="time",ylab="Sample path") lines(c(0,quake$time),0:nrow(quake),type="s") for(i in 1:N) lines(c(tail(quake$time,1),times[[i]]),nrow(quake)+0:length(times[[i]]), type="s",lty=2)
Simulate a renewal Hawkes (RHawkes) process with given renewal immigration distribution function, offspring density function and branching ratio.
simRHawkes(re.dist = rweibull, par.redist = list(shape = 1, scale = 1), ofspr.den = function(x, p.ofd) 1 / p.ofd * exp(-x / p.ofd), p.ofd = 1, branching.ratio = 0.5, cens = 1, B = 10, B0 = 50, max.ofspr.den = max(optimize(ofspr.den, c(0, cens), maximum = TRUE, p = p.ofd)$obj, ofspr.den(0, p.ofd), ofspr.den(cens, p.ofd)) * 1.1)
simRHawkes(re.dist = rweibull, par.redist = list(shape = 1, scale = 1), ofspr.den = function(x, p.ofd) 1 / p.ofd * exp(-x / p.ofd), p.ofd = 1, branching.ratio = 0.5, cens = 1, B = 10, B0 = 50, max.ofspr.den = max(optimize(ofspr.den, c(0, cens), maximum = TRUE, p = p.ofd)$obj, ofspr.den(0, p.ofd), ofspr.den(cens, p.ofd)) * 1.1)
re.dist |
A (vectorized) function. The immigrant renewal distribution function. |
par.redist |
A numeric vector. The parameters of the immigrant renewal distribution. |
ofspr.den |
A (vectorized) function. The offspring density function. |
p.ofd |
A numeric vector. The parameters of the offspring density. |
branching.ratio |
A scalar. The branching ratio parameter. |
cens |
A scalar. The censoring time. |
B |
A scalar. Tuning parameter for simulation of further immigrants. |
B0 |
A scalar. Tuning parameter for simulation of initial immigrants. |
max.ofspr.den |
A scalar. The maximum value of the offspring density function from 0
to |
The function works by simulating the arrival times of immigrants according to the renewal immigration distribution. The birth times of offspring from each immigrant are then simulated according to an inhomogeneous Poisson processes with appropriate intensity functions.
A numeric vector of all pooled events (immigration/birth) times of all generations 0, 1, ...
Feng Chen <[email protected]> Tom Stindl <[email protected]>
B <- 10; i <- 0; tms <- replicate(B, {cat(i<<-i+1,'\n'); simRHawkes(par.redist = list(shape = 3, scale = 1), p.ofd = 0.5, branching.ratio = 0.5, cens = 100) })
B <- 10; i <- 0; tms <- replicate(B, {cat(i<<-i+1,'\n'); simRHawkes(par.redist = list(shape = 3, scale = 1), p.ofd = 0.5, branching.ratio = 0.5, cens = 100) })
Simulate a renewal Hawkes (RHawkes) process with given renewal immigration distribution function, offspring density function and branching ratio.
simRHawkes1(re.dist = rweibull, par.redist = list(shape = 1, scale = 1), of.dis = "exp", par.ofdis = list(rate=1), branching.ratio = 0.5, cens = 1, B = 10, B0 = 50, flatten=TRUE)
simRHawkes1(re.dist = rweibull, par.redist = list(shape = 1, scale = 1), of.dis = "exp", par.ofdis = list(rate=1), branching.ratio = 0.5, cens = 1, B = 10, B0 = 50, flatten=TRUE)
re.dist |
A (vectorized) function. The immigrant renewal distribution function. |
par.redist |
A numeric vector. The parameters of the immigrant renewal distribution. |
of.dis |
A character string indicating the distribution for the offspring
birth times, which has to be the
'distname' part of the |
par.ofdis |
A list with named elements, giving the list of parameters of the offspring distribution, such as list(rate=1), list(shape=1,scale=1), etc. |
branching.ratio |
A scalar between 0 and 1, the branching ratio parameter. |
cens |
A scalar. The censoring time. |
B |
A scalar. Tuning parameter for simulation of further immigrants. |
B0 |
A scalar. Tuning parameter for simulation of initial immigrants. |
flatten |
A boolean scalar, which indicates whether the output events times should be flattened into an increasing sequence of times, or not (in which case the output is the immigrant arrival times, and the offspring birth times for different immigrants). |
The function works by simulating the arrival times of immigrants according to the renewal immigration distribution. The birth times of offspring from each immigrant are then simulated according to an inhomogeneous Poisson processes with appropriate intensity functions.
A numeric vector of pooled event (immigration/offspring birth) times of
all generations 0, 1, ..., if flatten=TRUE
; A list with two
components: immitimes
for immigrant arrival times, and
offspringtimes
for birth times of offspring due to different
immigrants.
Feng Chen <[email protected]> Tom Stindl <[email protected]>
tms <- simRHawkes1(par.redist = list(shape = 3, scale = 1), par.ofdis = list(rate=0.5), branching.ratio = 0.5, cens = 50) plot(stepfun(tms,0:length(tms)),do.points=FALSE,vertical=FALSE,xlim=c(0,50)) tms.clust <- simRHawkes1(par.redist = list(shape = 3, scale = 1), par.ofdis = list(rate=0.5), branching.ratio = 0.5, cens = 50,flatten=FALSE) plot(c(0,50),c(0, 1+(nt <-length(it <- tms.clust$immitimes))), type="n",xlab="time",ylab="cluster") segments(x0 = it,y0=-0.2,y1=0.2) for(i in 1:nt) segments(x0 = c(it[i],it[i]+tms.clust$offspringtimes[[i]]), y0=i-0.2,y1=i+0.2) abline(h=0:(nt+1),col="light gray",lty=2) segments(x0=unlist(lapply(1:nt,function(i)c(it[i],it[i]+tms.clust$offspringtimes[[i]]))), y0=nt+1-0.2,y1=nt+1+0.2)
tms <- simRHawkes1(par.redist = list(shape = 3, scale = 1), par.ofdis = list(rate=0.5), branching.ratio = 0.5, cens = 50) plot(stepfun(tms,0:length(tms)),do.points=FALSE,vertical=FALSE,xlim=c(0,50)) tms.clust <- simRHawkes1(par.redist = list(shape = 3, scale = 1), par.ofdis = list(rate=0.5), branching.ratio = 0.5, cens = 50,flatten=FALSE) plot(c(0,50),c(0, 1+(nt <-length(it <- tms.clust$immitimes))), type="n",xlab="time",ylab="cluster") segments(x0 = it,y0=-0.2,y1=0.2) for(i in 1:nt) segments(x0 = c(it[i],it[i]+tms.clust$offspringtimes[[i]]), y0=i-0.2,y1=i+0.2) abline(h=0:(nt+1),col="light gray",lty=2) segments(x0=unlist(lapply(1:nt,function(i)c(it[i],it[i]+tms.clust$offspringtimes[[i]]))), y0=nt+1-0.2,y1=nt+1+0.2)
A financial data set containing the mid-price changes of the AUD/USD foreign exchange rate during the trading week from 20:00:00 GMT on Sunday 19/07/2015 to 21:00:00 GMT Friday 24/07/2015.
data(tms)
data(tms)
The format is a list of the arrival times of mid price changes that occur every hour in 121 non-overlapping windows.
Times of arrivals of mid-price changes is listed together in ascending order.
Simulated by a call to the function simHawkes1
.
data(tms) ## number of non over-lapping hourly windows length(tms)
data(tms) ## number of non over-lapping hourly windows length(tms)