Package 'MRHawkes' reference manual

Title:	Multivariate Renewal Hawkes Process
Description:	Simulate a (bivariate) multivariate renewal Hawkes (MRHawkes) self-exciting process, with given immigrant hazard rate functions and offspring density function. Calculate the likelihood of a MRHawkes process with given hazard rate functions and offspring density function for an (increasing) sequence of event times. Calculate the Rosenblatt residuals of the event times. Predict future event times based on observed event times up to a given time. For details see Stindl and Chen (2018) <doi:10.1016/j.csda.2018.01.021>.
Authors:	Tom Stindl and Feng Chen
Maintainer:	Tom Stindl <[email protected]>
License:	GPL (>= 2)
Version:	1.0
Built:	2025-02-12 04:40:46 UTC
Source:	https://github.com/cran/MRHawkes

Multivariate Renewal Hawkes Process

Description

Simulate the (bivariate) multivariate renewal Hawkes (MRHawkes) process with a given distribution for the two waiting times between immigrants, given offspring density functions, and also the branching ratios. Calculation of the likelihood of a MRHawkes process model with a given sequence of (distinct) event times and labels up to a censoring time. Calculate the Rosenblatt residuals of fitting an MRHawkes process model to a sequence of event times and labels. Calculate the (filtering) distribution of the index of the most recent immigrant. Predict the time of the next event since the censoring time. Predict event times from the censoring time to a future time point.

Details

The DESCRIPTION file:

Package:	MRHawkes
Type:	Package
Title:	Multivariate Renewal Hawkes Process
Version:	1.0
Date:	2018-08-15
Author:	Tom Stindl and Feng Chen
Maintainer:	Tom Stindl <[email protected]>
Description:	Simulate a (bivariate) multivariate renewal Hawkes (MRHawkes) self-exciting process, with given immigrant hazard rate functions and offspring density function. Calculate the likelihood of a MRHawkes process with given hazard rate functions and offspring density function for an (increasing) sequence of event times. Calculate the Rosenblatt residuals of the event times. Predict future event times based on observed event times up to a given time. For details see Stindl and Chen (2018) <doi:10.1016/j.csda.2018.01.021>.
License:	GPL (>= 2)
Depends:	IHSEP, stats
NeedsCompilation:	no
Packaged:	2018-08-20 01:02:54 UTC; tomstindl
Date/Publication:	2018-08-20 11:20:02 UTC
Repository:	https://tomstindl.r-universe.dev
RemoteUrl:	https://github.com/cran/MRHawkes
RemoteRef:	HEAD
RemoteSha:	e9b8c8344331ba53ad1df61b30a5f29d3ac95643

Index of help topics:

MRHawkes                Multivariate Renewal Hawkes Process
TmllMRH                 Minus loglikelihood of an (bivariate) MRHawkes
                        model with truncated most recent immigrant
                        probabilities
fivaqks                 Fiji and Vanuatu Earthquake Data
mllMRH                  Minus loglikelihood of an (bivariate) MRHawkes
                        model
mllMRH1                 Minus loglikelihood of an (bivariate) MRHawkes
                        model with most recent immigrant probabilities
mllMRH2                 Minus loglikelihood of an (bivariate) MRHawkes
                        model with Rosenblatt residuals
predDen                 MRHawkes (bivariate) predictive density
                        function
simMRHawkes             Simulate an (bivariate) renewal Hawkes
                        (MRHawkes) process
simNSMHP                Simulate a (bivariate) non-stationary
                        multivariate Hawkes process (NSMHP)
simpred                 Simulate a fitted (bivariate) MRHawkes process
                        model
typeRes                 Minus loglikelihood of an (bivariate) MRHawkes
                        model with Universal residuals

Author(s)

Tom Stindl <[email protected]> Feng Chen <[email protected]>

References

Chen, F. and Stindl, T. (2017). Direct Likelihood Evaluation for the Renewal Hawkes Process. Journal of Computational and Graphical Statistics.

Stindl, T. & Chen, F. (2018). Likelihood based inference for the multivariate renewal Hawkes process. Computational Statistics and Data Analysis.

Wheatley, S., Filimonov, V., and Sornette, D. (2016) The Hawkes process with renewal immigration & its estimation with an EM algorithm. Computational Statistics and Data Analysis. 94: 120-135.

Fiji and Vanuatu Earthquake Data

Description

Times and magnitudes (Richter scale) of earthquakes in the regions of Fiji and Vanuatu, for the period 1990-2015.

Usage

data(fivaqks)data(fivaqks)

Format

A data set containing 22 variables for earthquakes around Fiji and Vanuatu from 1991 to 2015.

time: Time of quake
latitude: Latitude of the quake
longitude: Longitude of the quake
mag: Magnitude of the quake

Source

United States Geological Survey (USGS)

Examples

  data(fivaqks)
data(fivaqks)

Minus loglikelihood of an (bivariate) MRHawkes model

Description

Calculates the minus loglikelihood of an (bivariate) MRHawkes model with given immigration hazard functions $\mu$ , offspring density functions $h$ and bracnhing ratios $\eta$ for the event times and types data on the interval $[0,cens]$ .

Usage

mllMRH(data, cens, par,
       h1.fn = function(x, p) 1 / p * exp( - x / p),
       h2.fn = function(x, p) 1 / p * exp( - x / p),
       mu1.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))
       },
       mu2.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))
       },
       H1.fn = function(x, p) pexp(x, rate = 1 / p),
       H2.fn = function(x, p) pexp(x, rate = 1 / p),
       Mu1.fn = function(x, p){
         - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, 
                   log.p = TRUE)
       },
       Mu2.fn = function(x, p){
         - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, 
                    log.p = TRUE)
       })
mllMRH(data, cens, par,
       h1.fn = function(x, p) 1 / p * exp( - x / p),
       h2.fn = function(x, p) 1 / p * exp( - x / p),
       mu1.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))
       },
       mu2.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))
       },
       H1.fn = function(x, p) pexp(x, rate = 1 / p),
       H2.fn = function(x, p) pexp(x, rate = 1 / p),
       Mu1.fn = function(x, p){
         - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, 
                   log.p = TRUE)
       },
       Mu2.fn = function(x, p){
         - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, 
                    log.p = TRUE)
       })

Arguments

`data`	A two column matrix. The first column contains the event times sorted in ascending order. The second column contains the corresponding event type with the label one or two.
`cens`	A scalar. The censoring time.
`par`	A numeric vector. Contains the ten parameters of the model, in order of the immigration parameters $\mu(.)$ for the two renewal distributions, the two offspring parameters $h(.)$ and lastly the four branching ratios $\eta$ .
`h1.fn`	A (vectorized) function. The offspring density function for type one events.
`h2.fn`	A (vectorized) function. The offspring density function for type two events.
`mu1.fn`	A (vectorized) function. The immigration hazard function for events of type one.
`mu2.fn`	A (vectorized) function. The immigration hazard function for events of type two.
`H1.fn`	A (vectorized) function. Its value at `t` gives the integral of the offspring density function from 0 to `t` for type one events.
`H2.fn`	A (vectorized) function. Its value at `t` gives the integral of the offspring density function from 0 to `t` for type two events.
`Mu1.fn`	A (vectorized) function. Its value at `t` gives the integral of the immigrant hazard function from 0 to `t` for type one events.
`Mu2.fn`	A (vectorized) function. Its value at `t` gives the integral of the immigrant hazard function from 0 to `t` for type two events.

Value

Value of the negative log-liklihood.

Author(s)

Tom Stindl <[email protected]> Feng Chen <[email protected]>

Examples

  ## Magnitude 5.5 or greater earthquakes over the 25 year period from 
  ## 01/01/1991 to 31/12/2015.
  data(fivaqks); 
  near.fiji <- grep("Fiji", fivaqks$place)
  near.vanuatu <- grep("Vanuatu", fivaqks$place)
  t.beg <- strptime("1991-01-01 00:00:00", "%Y-%m-%d %H:%M:%S", tz = "UTC")
  t.end <- strptime("2015-12-31 23:59:59", "%Y-%m-%d %H:%M:%S", tz = "UTC")
  t0 <- 0
  t1 <- as.numeric(t.end - t.beg)
  tms <- strptime(fivaqks$time, "%Y-%m-%dT%H:%M:%OSZ", tz = "UTC")
  ts <- as.numeric(tms[-1] - t.beg)
  ts <- c(as.numeric(tms[1] - t.beg)/24, ts)
  ts.fi <- ts[near.fiji]; ts.fi <- ts.fi[ts.fi >= 0 & ts.fi <= t1]
  ts.va <- ts[near.vanuatu]; ts.va <- ts.va[ts.va >=0 & ts.va <= t1]
  ts.c <- c(ts.fi, ts.va)
  z.c <- c(rep(1, times = length(ts.fi)), rep(2, times = length(ts.va)))
  o <- order(ts.c)
  data <- cbind(ts.c[o], z.c[o])
  
  ## calculate the minus loglikelihood of an (bivariate) MRHawkes with some 
  ## parameters the default hazard functions and density functions are Weibull 
  ## and exponential respectivley
  mllMRH(data, cens = t1, par = c(0.488, 20.10, 0.347, 9.53, 461, 720,
                                  0.472, 0.293, 0.399, -0.0774))
                                  
  ## 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.488, 20.10, 0.347, 9.53, 461, 720,
                                  0.472, 0.293, 0.399, -0.0774), 
                             mllMRH, data = data, cens = t1,
                             control = list(maxit = 5000, trace = TRUE),
                             hessian = TRUE)
    )
    ## point estimate by MLE
    est$par
    ## standard error estimates:
    diag(solve(est$hessian))^0.5
  
## Magnitude 5.5 or greater earthquakes over the 25 year period from 
  ## 01/01/1991 to 31/12/2015.
  data(fivaqks); 
  near.fiji <- grep("Fiji", fivaqks$place)
  near.vanuatu <- grep("Vanuatu", fivaqks$place)
  t.beg <- strptime("1991-01-01 00:00:00", "%Y-%m-%d %H:%M:%S", tz = "UTC")
  t.end <- strptime("2015-12-31 23:59:59", "%Y-%m-%d %H:%M:%S", tz = "UTC")
  t0 <- 0
  t1 <- as.numeric(t.end - t.beg)
  tms <- strptime(fivaqks$time, "%Y-%m-%dT%H:%M:%OSZ", tz = "UTC")
  ts <- as.numeric(tms[-1] - t.beg)
  ts <- c(as.numeric(tms[1] - t.beg)/24, ts)
  ts.fi <- ts[near.fiji]; ts.fi <- ts.fi[ts.fi >= 0 & ts.fi <= t1]
  ts.va <- ts[near.vanuatu]; ts.va <- ts.va[ts.va >=0 & ts.va <= t1]
  ts.c <- c(ts.fi, ts.va)
  z.c <- c(rep(1, times = length(ts.fi)), rep(2, times = length(ts.va)))
  o <- order(ts.c)
  data <- cbind(ts.c[o], z.c[o])
  
  ## calculate the minus loglikelihood of an (bivariate) MRHawkes with some 
  ## parameters the default hazard functions and density functions are Weibull 
  ## and exponential respectivley
  mllMRH(data, cens = t1, par = c(0.488, 20.10, 0.347, 9.53, 461, 720,
                                  0.472, 0.293, 0.399, -0.0774))
                                  
  ## 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.488, 20.10, 0.347, 9.53, 461, 720,
                                  0.472, 0.293, 0.399, -0.0774), 
                             mllMRH, data = data, cens = t1,
                             control = list(maxit = 5000, trace = TRUE),
                             hessian = TRUE)
    )
    ## point estimate by MLE
    est$par
    ## standard error estimates:
    diag(solve(est$hessian))^0.5

Minus loglikelihood of an (bivariate) MRHawkes model with most recent immigrant probabilities

Description

Calculates the minus loglikelihood of an (bivariate) RHawkes model with given immigration hazard functions $\mu$ , common offspring density functions $h$ and bracnhing ratios $\eta$ for event times and event types data on interval $[0,cens]$ . The same as mllMRH although this version also returns the most recent immigrant probabilities at the censoring.

Usage

mllMRH1(data, cens, par,
       h1.fn = function(x, p) 1 / p * exp( - x / p),
       h2.fn = function(x, p) 1 / p * exp( - x / p),
       mu1.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))
       },
       mu2.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))
       },
       H1.fn = function(x, p) pexp(x, rate = 1 / p),
       H2.fn = function(x, p) pexp(x, rate = 1 / p),
       Mu1.fn = function(x, p){
         - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, 
                   log.p = TRUE)
       },
       Mu2.fn = function(x, p){
         - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, 
                    log.p = TRUE)
       })
mllMRH1(data, cens, par,
       h1.fn = function(x, p) 1 / p * exp( - x / p),
       h2.fn = function(x, p) 1 / p * exp( - x / p),
       mu1.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))
       },
       mu2.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))
       },
       H1.fn = function(x, p) pexp(x, rate = 1 / p),
       H2.fn = function(x, p) pexp(x, rate = 1 / p),
       Mu1.fn = function(x, p){
         - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, 
                   log.p = TRUE)
       },
       Mu2.fn = function(x, p){
         - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, 
                    log.p = TRUE)
       })

Arguments

`data`	A two column matrix. The first column contains the event times sorted in ascending order. The second column contains the corresponding event type with the label one or two.
`cens`	A scalar. The censoring time.
`par`	A numeric vector. Contains the ten parameters of the model, in order of the immigration parameters $\mu(.)$ for the two renewal distributions, the two offspring parameters $h(.)$ and lastly the four branching ratios $\eta$ .
`h1.fn`	A (vectorized) function. The offspring density function for type one events.
`h2.fn`	A (vectorized) function. The offspring density function for type two events.
`mu1.fn`	A (vectorized) function. The immigration hazard function for events of type one.
`mu2.fn`	A (vectorized) function. The immigration hazard function for events of type two.
`H1.fn`	A (vectorized) function. Its value at `t` gives the integral of the offspring density function from 0 to `t` for type one events.
`H2.fn`	A (vectorized) function. Its value at `t` gives the integral of the offspring density function from 0 to `t` for type two events.
`Mu1.fn`	A (vectorized) function. Its value at `t` gives the integral of the immigrant hazard function from 0 to `t` for type one events.
`Mu2.fn`	A (vectorized) function. Its value at `t` gives the integral of the immigrant hazard function from 0 to `t` for type two events.

Value

`mll`	minus log-likelihood
`p`	most recent immigrant probabilities at the censoring time
`n`	number of events

Author(s)

Tom Stindl <[email protected]> Feng Chen <[email protected]>

Examples


  data <- cbind(sort(runif(1000,0,1000)), 
                sample(1:2, size = 1000, replace = TRUE))
  tmp <- mllMRH1(data = data, cens = 1001, 
                 par = c(3,1.2,1/3,0.2,1,1,0.5,0.2,0.2,0.3))
  ## last immigrant probabilities
  lip <- tmp$p
  ## sample last immigrant at censoring time for component one and 
  ## component two respectively
  c(sample(0:1000, 1, replace = TRUE, prob = rowSums(lip)), 
  sample(0:1000, 1, replace = TRUE, prob = colSums(lip)))
  
data <- cbind(sort(runif(1000,0,1000)), 
                sample(1:2, size = 1000, replace = TRUE))
  tmp <- mllMRH1(data = data, cens = 1001, 
                 par = c(3,1.2,1/3,0.2,1,1,0.5,0.2,0.2,0.3))
  ## last immigrant probabilities
  lip <- tmp$p
  ## sample last immigrant at censoring time for component one and 
  ## component two respectively
  c(sample(0:1000, 1, replace = TRUE, prob = rowSums(lip)), 
  sample(0:1000, 1, replace = TRUE, prob = colSums(lip)))

Minus loglikelihood of an (bivariate) MRHawkes model with Rosenblatt residuals

Description

Calculates the minus loglikelihood of an (bivariate) RHawkes model with given immigration hazard functions $\mu$ , common offspring density functions $h$ and bracnhing ratios $\eta$ for event times and event types data on interval $[0,cens]$ . The same as mllMRH although this version also returns the Rosenblatt residuals for goodness-of-fit assessment of the event times.

Usage

mllMRH2(data, cens, par,
       h1.fn = function(x, p) 1 / p * exp( - x / p),
       h2.fn = function(x, p) 1 / p * exp( - x / p),
       mu1.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))
       },
       mu2.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))
       },
       H1.fn = function(x, p) pexp(x, rate = 1 / p),
       H2.fn = function(x, p) pexp(x, rate = 1 / p),
       Mu1.fn = function(x, p){
         - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, 
                   log.p = TRUE)
       },
       Mu2.fn = function(x, p){
         - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, 
                    log.p = TRUE)
       })
mllMRH2(data, cens, par,
       h1.fn = function(x, p) 1 / p * exp( - x / p),
       h2.fn = function(x, p) 1 / p * exp( - x / p),
       mu1.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))
       },
       mu2.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))
       },
       H1.fn = function(x, p) pexp(x, rate = 1 / p),
       H2.fn = function(x, p) pexp(x, rate = 1 / p),
       Mu1.fn = function(x, p){
         - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, 
                   log.p = TRUE)
       },
       Mu2.fn = function(x, p){
         - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, 
                    log.p = TRUE)
       })

Arguments

`data`	A two column matrix. The first column contains the event times sorted in ascending order. The second column contains the corresponding event type with the label one or two.
`cens`	A scalar. The censoring time.
`par`	A numeric vector. Contains the ten parameters of the model, in order of the immigration parameters $\mu(.)$ for the two renewal distributions, the two offspring parameters $h(.)$ and lastly the four branching ratios $\eta$ .
`h1.fn`	A (vectorized) function. The offspring density function for type one events.
`h2.fn`	A (vectorized) function. The offspring density function for type two events.
`mu1.fn`	A (vectorized) function. The immigration hazard function for events of type one.
`mu2.fn`	A (vectorized) function. The immigration hazard function for events of type two.
`H1.fn`	A (vectorized) function. Its value at `t` gives the integral of the offspring density function from 0 to `t` for type one events.
`H2.fn`	A (vectorized) function. Its value at `t` gives the integral of the offspring density function from 0 to `t` for type two events.
`Mu1.fn`	A (vectorized) function. Its value at `t` gives the integral of the immigrant hazard function from 0 to `t` for type one events.
`Mu2.fn`	A (vectorized) function. Its value at `t` gives the integral of the immigrant hazard function from 0 to `t` for type two events.

Details

Calculate the MRHawkes point process Rosenblatt residuals

Value

`mll`	minus log-likelihood
`W`	Rosenblatt residuals of observed event times

Author(s)

Tom Stindl <[email protected]> Feng Chen <[email protected]>

Examples


  n <- 1000
  data <- cbind(sort(runif(n,0,1000)), 
                sample(1:2, size = n, replace = TRUE))
  tmp <- mllMRH2(data = data, cens = 1001, 
                 par = c(1,1,1,1,1,1,0.5,0.2,0.2,0.3))              
  pp <- ppoints(n)
  par(mfrow=c(1,2))
  plot(quantile(tmp$W,prob=pp),pp,type="l",
       main="Uniform QQ plot",
       xlab="Sample quantiles",ylab="Theoretical quantiles")
  abline(a = 0, b = 1, col = 2)
  a <- acf(tmp$W, main = "ACF Plot")
  ks.test(tmp$W,"punif")
  Box.test(tmp$W,lag=tail(a$lag,1))
  
n <- 1000
  data <- cbind(sort(runif(n,0,1000)), 
                sample(1:2, size = n, replace = TRUE))
  tmp <- mllMRH2(data = data, cens = 1001, 
                 par = c(1,1,1,1,1,1,0.5,0.2,0.2,0.3))              
  pp <- ppoints(n)
  par(mfrow=c(1,2))
  plot(quantile(tmp$W,prob=pp),pp,type="l",
       main="Uniform QQ plot",
       xlab="Sample quantiles",ylab="Theoretical quantiles")
  abline(a = 0, b = 1, col = 2)
  a <- acf(tmp$W, main = "ACF Plot")
  ks.test(tmp$W,"punif")
  Box.test(tmp$W,lag=tail(a$lag,1))

MRHawkes (bivariate) predictive density function

Description

Calculates the predictive density of the next event time after the censoring time cens based on the observations over the interval [0,cens].

Usage

predDen(x, data, cens, par, 
        h1.fn = function(x, p) 1 / p * exp( - x / p),
        h2.fn = function(x, p) 1 / p * exp( - x / p),
        mu1.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))
        },
        mu2.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))
        },
        H1.fn = function(x, p) pexp(x, rate = 1 / p),
        H2.fn = function(x, p) pexp(x, rate = 1 / p),
        Mu1.fn = function(x, p){
         - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, 
                    log.p = TRUE)
        },
        Mu2.fn = function(x, p){
         - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, 
                    log.p = TRUE)
        })
predDen(x, data, cens, par, 
        h1.fn = function(x, p) 1 / p * exp( - x / p),
        h2.fn = function(x, p) 1 / p * exp( - x / p),
        mu1.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))
        },
        mu2.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))
        },
        H1.fn = function(x, p) pexp(x, rate = 1 / p),
        H2.fn = function(x, p) pexp(x, rate = 1 / p),
        Mu1.fn = function(x, p){
         - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, 
                    log.p = TRUE)
        },
        Mu2.fn = function(x, p){
         - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, 
                    log.p = TRUE)
        })

Arguments

`x`	A scalar. The amount of time after the censoring tine `cens`.
`data`	A two column matrix. The first column contains the event times sorted in ascending order. The second column contains the corresponding event type with the label one or two.
`cens`	A scalar. The censoring time.
`par`	A numeric vector containing the twelve parameters of the model, in order of the immigration parameters $\mu(.)$ for the two component distributions, the four offspring parameters $h(.)$ and lastly the four branching ratios $\eta$ .
`h1.fn`	A (vectorized) function. The offspring density function for type one events.
`h2.fn`	A (vectorized) function. The offspring density function for type two events.
`mu1.fn`	A (vectorized) function. The immigration hazard function for events of type one.
`mu2.fn`	A (vectorized) function. The immigration hazard function for events of type two.
`H1.fn`	A (vectorized) function. Its value at `t` gives the integral of the offspring density function from 0 to `t` for type one events.
`H2.fn`	A (vectorized) function. Its value at `t` gives the integral of the offspring density function from 0 to `t` for type two events.
`Mu1.fn`	A (vectorized) function. Its value at `t` gives the integral of the immigrant hazard function from 0 to `t` for type one events.
`Mu2.fn`	A (vectorized) function. Its value at `t` gives the integral of the immigrant hazard function from 0 to `t` for type two events.

Value

The predictive density of the next event time evaluated at x.

Author(s)

Tom Stindl <[email protected]> Feng Chen <[email protected]>

Examples


  ## Magnitude 5.5 or greater earthquakes over the 25 year period from 
  ## 01/01/1991 to 31/12/2015.  
  data(fivaqks); 
  near.fiji <- grep("Fiji", fivaqks$place)
  near.vanuatu <- grep("Vanuatu", fivaqks$place)
  t.beg <- strptime("1991-01-01 00:00:00", "%Y-%m-%d %H:%M:%S", tz = "UTC")
  t.end <- strptime("2015-12-31 23:59:59", "%Y-%m-%d %H:%M:%S", tz = "UTC")
  t0 <- 0
  t1 <- as.numeric(t.end - t.beg)
  tms <- strptime(fivaqks$time, "%Y-%m-%dT%H:%M:%OSZ", tz = "UTC")
  ts <- as.numeric(tms[-1] - t.beg)
  ts <- c(as.numeric(tms[1] - t.beg)/24, ts)
  ts.fi <- ts[near.fiji]; ts.fi <- ts.fi[ts.fi >= 0 & ts.fi <= t1]
  ts.va <- ts[near.vanuatu]; ts.va <- ts.va[ts.va >=0 & ts.va <= t1]
  ts.c <- c(ts.fi, ts.va)
  z.c <- c(rep(1, times = length(ts.fi)), rep(2, times = length(ts.va)))
  o <- order(ts.c)
  data <- cbind(ts.c[o], z.c[o])
  curve(predDen(x, data = data, cens = t1, 
                 par = c(0.488, 20.10, 0.347, 9.53, 461, 720, 
                         0.472, 0.293, 0.399, -0.0774)) 
        ,0 ,200, col = "red", lwd = 2, ylab = "Density")
  
## Magnitude 5.5 or greater earthquakes over the 25 year period from 
  ## 01/01/1991 to 31/12/2015.  
  data(fivaqks); 
  near.fiji <- grep("Fiji", fivaqks$place)
  near.vanuatu <- grep("Vanuatu", fivaqks$place)
  t.beg <- strptime("1991-01-01 00:00:00", "%Y-%m-%d %H:%M:%S", tz = "UTC")
  t.end <- strptime("2015-12-31 23:59:59", "%Y-%m-%d %H:%M:%S", tz = "UTC")
  t0 <- 0
  t1 <- as.numeric(t.end - t.beg)
  tms <- strptime(fivaqks$time, "%Y-%m-%dT%H:%M:%OSZ", tz = "UTC")
  ts <- as.numeric(tms[-1] - t.beg)
  ts <- c(as.numeric(tms[1] - t.beg)/24, ts)
  ts.fi <- ts[near.fiji]; ts.fi <- ts.fi[ts.fi >= 0 & ts.fi <= t1]
  ts.va <- ts[near.vanuatu]; ts.va <- ts.va[ts.va >=0 & ts.va <= t1]
  ts.c <- c(ts.fi, ts.va)
  z.c <- c(rep(1, times = length(ts.fi)), rep(2, times = length(ts.va)))
  o <- order(ts.c)
  data <- cbind(ts.c[o], z.c[o])
  curve(predDen(x, data = data, cens = t1, 
                 par = c(0.488, 20.10, 0.347, 9.53, 461, 720, 
                         0.472, 0.293, 0.399, -0.0774)) 
        ,0 ,200, col = "red", lwd = 2, ylab = "Density")

Simulate an (bivariate) renewal Hawkes (MRHawkes) process

Description

Simulate an (bivairate) renewal Hawkes (MRHawkes) process with given renewal immigration distribution functions $\mu$ , offspring density functions $h$ and branching ratios $\eta$ using the cascading structure of the process.

Usage

simMRHawkes(re.dist1 = rweibull, par.redist1 = list(shape = 3, scale = 1.2), 
            re.dist2 = rweibull, par.redist2 = list(shape = 1/3, scale = 0.2), 
            h1.fn = function(x, p.h1) 1/p.h1 * exp(-x/p.h1), 
            h2.fn = function(x, p.h2) 1/p.h2 * exp(-x/p.h2), 
            p.h1 = 1, p.h2 = 1, 
            eta11 = 0.3, eta12 = 0.1, eta21 = 0.1, eta22 = 0.3, cens = 100, 
            B = 10, B0 = 50, 
            max.h1 = max(optimize(h1.fn, c(0, cens), maximum = TRUE, p = p.h1)$obj, 
                         h1.fn(0, p.h1), h1.fn(cens, p.h1)) * 1.1, 
            max.h2 = max(optimize(h2.fn, c(0, cens), maximum = TRUE, p = p.h2)$obj, 
                         h2.fn(0, p.h2), h2.fn(cens, p.h2)) * 1.1)
simMRHawkes(re.dist1 = rweibull, par.redist1 = list(shape = 3, scale = 1.2), 
            re.dist2 = rweibull, par.redist2 = list(shape = 1/3, scale = 0.2), 
            h1.fn = function(x, p.h1) 1/p.h1 * exp(-x/p.h1), 
            h2.fn = function(x, p.h2) 1/p.h2 * exp(-x/p.h2), 
            p.h1 = 1, p.h2 = 1, 
            eta11 = 0.3, eta12 = 0.1, eta21 = 0.1, eta22 = 0.3, cens = 100, 
            B = 10, B0 = 50, 
            max.h1 = max(optimize(h1.fn, c(0, cens), maximum = TRUE, p = p.h1)$obj, 
                         h1.fn(0, p.h1), h1.fn(cens, p.h1)) * 1.1, 
            max.h2 = max(optimize(h2.fn, c(0, cens), maximum = TRUE, p = p.h2)$obj, 
                         h2.fn(0, p.h2), h2.fn(cens, p.h2)) * 1.1)

Arguments

`re.dist1`	The renewal distribution for type one events.
`re.dist2`	The renewal distribution for type two events.
`par.redist1`	A numeric list. The parameters of the renewal distribution for type one events.
`par.redist2`	A numeric list. The parameters of the renewal distribution for type two events.
`h1.fn`	A (vectorized) function. The offspring density function for type one events.
`h2.fn`	A (vectorized) function. The offspring density function for type two events.
`p.h1`	A numeric vector. The paramters of the offspring density for type one events.
`p.h2`	A numeric vector. The paramters of the offspring density for type two events.
`eta11`	A numeric scalar. The self-exciting branching ratio for type one events.
`eta12`	A numeric scalar. The cross-exciting branching ratio for type one events due to the effects of a type two event.
`eta21`	A numeric scalar. The cross-exciting branching ratio for type two events due to the effects of a type one event.
`eta22`	A numeric scalar. The self-exciting branching ratio for type two events.
`cens`	A scalar. The censoring time.
`B`	A numeric scalar. Tuning parameter
`B0`	A numeric scalar. Tuning parameter
`max.h1`	A numeric scalar. The maximum value of the offspring density for type one events.
`max.h2`	A numeric scalar. The maximum value of the offspring density for type two events.

Details

The function works by simulating the arrival times of immigrants accoridng to the respective renewal immigration distribution for each event type. The birth times ofoffspring from each immigrant are then simulated according to an non-stationary multivariate Hawkes Process (NSMHP) with appropriate baseline and excitation functions.

Value

A numeric matrix. The row coloumn contains the event times in ascending order while the second coloumn contains the corresponding event type.

Author(s)

Tom Stindl <[email protected]> Feng Chen <[email protected]>

Examples

  B <- 10; i <- 0;
  data <- replicate(B, 
                    {cat(i<<-i+1,'\n'); 
                    simMRHawkes(re.dist1 = rweibull, 
                                par.redist1 = list(shape = 3, scale = 1.2),
                                re.dist2 = rweibull, 
                                par.redist2 = list(shape = 1 / 3, scale = 0.2),
                                p.h1 = 1, p.h2 = 1,
                                eta11 = 0.3, eta12 = 0.1, 
                                eta21 = 0.1, eta22 = 0.3,
                                cens = 100)
                    })
B <- 10; i <- 0;
  data <- replicate(B, 
                    {cat(i<<-i+1,'\n'); 
                    simMRHawkes(re.dist1 = rweibull, 
                                par.redist1 = list(shape = 3, scale = 1.2),
                                re.dist2 = rweibull, 
                                par.redist2 = list(shape = 1 / 3, scale = 0.2),
                                p.h1 = 1, p.h2 = 1,
                                eta11 = 0.3, eta12 = 0.1, 
                                eta21 = 0.1, eta22 = 0.3,
                                cens = 100)
                    })

Simulate a (bivariate) non-stationary multivariate Hawkes process (NSMHP)

Description

Simulate a bivariate non-stationary multivariate Hawkes process (NSMHP) with given given baseline intensity functions and self-excitation functions using the cascading structure of the process.

Usage

simNSMHP(TT = 100,
        nu1 = function(t) 0.6*exp(-t),
        nu2 = function(t) 0.2*exp(-t),
        g11 = function(t) 0.6*exp(-t),
        g12 = function(t) 0.2*exp(-t),
        g21 = function(t) 0.1*exp(-t),
        g22 = function(t) 0.5*exp(-t))
simNSMHP(TT = 100,
        nu1 = function(t) 0.6*exp(-t),
        nu2 = function(t) 0.2*exp(-t),
        g11 = function(t) 0.6*exp(-t),
        g12 = function(t) 0.2*exp(-t),
        g21 = function(t) 0.1*exp(-t),
        g22 = function(t) 0.5*exp(-t))

Arguments

`TT`	A scalar. The censoring time.
`nu1`	Basline intensity function for type one events.
`nu2`	Basline intensity function for type two events.
`g11`	Self-exciting function for type one events given the parent is a type two event.
`g12`	Cross-exciting function for type one events given the parent is a type two event.
`g21`	Cross-exciting function for type two events given the parent is a type one event.
`g22`	Self-exciting function for type two events given the parent is a type two event.

Details

The function works by simulating generation 0 events according to independent Poisson processes with the baseline intensity functions; then keep simulating future generation events as long as the number of the previous generation events of any type is non-zero. For each event type, we simulate these events according to M independent Poisson processes with the appropriate excitation intensity. When this recursive process stops, return events of all generations with their respective type labels as the events of the NSMHP on the interval (0,T].

Value

`offspr1`	All offspring events of type one
`offspr2`	All offspring events of type two

Author(s)

Tom Stindl <[email protected]> Feng Chen <[email protected]>

Examples

  B <- 10; i <- 0;
  data <- replicate(B, 
                    {cat(i<<-i+1,'\n'); 
                    simNSMHP(TT = 100)
                    })
B <- 10; i <- 0;
  data <- replicate(B, 
                    {cat(i<<-i+1,'\n'); 
                    simNSMHP(TT = 100)
                    })

Simulate a fitted (bivariate) MRHawkes process model

Description

Simulate a fitted bivariate MRHawkes process model after the censoring time cens to a future time point cens.tilde using the cascading structure of the process.

Usage

  simpred(data, par, cens, cens.tilde, 
          re.dist1 = rweibull, 
          par.redist1 = list(shape = par[1], scale = par[2]),
          re.dist2 = rweibull, 
          par.redist2 = list(shape = par[3], scale = par[4]),
          h1.fn = function(x, p.h1) 1 / p.h1 * exp( - x / p.h1),
          h2.fn = function(x, p.h2) 1 / p.h2 * exp( - x / p.h2),
          p.h1 = par[5], p.h2 = par[6],
          eta11 = par[7], eta12 = par[8], eta21 = par[9], eta22 = par[10],
          B = 10, B0 = 50, pnp1 = NULL,
          max.h1 = max(optimize(h1.fn, c(0, cens.tilde - cens), maximum = TRUE,
                               p = p.h1)$obj, h1.fn(0, p.h1), 
                      h1.fn(cens.tilde - cens, p.h1)) * 1.1,
          max.h2 = max(optimize(h2.fn, c(0, cens.tilde - cens), maximum = TRUE,
                               p = p.h2)$obj, h2.fn(0, p.h2), 
                      h2.fn(cens.tilde - cens, p.h2)) * 1.1)
simpred(data, par, cens, cens.tilde, 
          re.dist1 = rweibull, 
          par.redist1 = list(shape = par[1], scale = par[2]),
          re.dist2 = rweibull, 
          par.redist2 = list(shape = par[3], scale = par[4]),
          h1.fn = function(x, p.h1) 1 / p.h1 * exp( - x / p.h1),
          h2.fn = function(x, p.h2) 1 / p.h2 * exp( - x / p.h2),
          p.h1 = par[5], p.h2 = par[6],
          eta11 = par[7], eta12 = par[8], eta21 = par[9], eta22 = par[10],
          B = 10, B0 = 50, pnp1 = NULL,
          max.h1 = max(optimize(h1.fn, c(0, cens.tilde - cens), maximum = TRUE,
                               p = p.h1)$obj, h1.fn(0, p.h1), 
                      h1.fn(cens.tilde - cens, p.h1)) * 1.1,
          max.h2 = max(optimize(h2.fn, c(0, cens.tilde - cens), maximum = TRUE,
                               p = p.h2)$obj, h2.fn(0, p.h2), 
                      h2.fn(cens.tilde - cens, p.h2)) * 1.1)

Arguments

`data`	A two column matrix. The first column contains the event times sorted in ascending order. The second column contains the corresponding event type with the label one or two.
`par`	A numeric vector. Contains the ten parameters of the model, in order of the immigration parameters $\mu(.)$ for the two renewal distributions, the two offspring parameters $h(.)$ and lastly the four branching ratios $\eta$ .
`cens`	A scalar. The censoring time.
`cens.tilde`	A scalar. The time that the simulation run uptil.
`re.dist1`	The renewal distribution for type one events.
`re.dist2`	The renewal distribution for type two events.
`par.redist1`	A numeric list. The parameters of the renewal distribution for type one events.
`par.redist2`	A numeric list. The parameters of the renewal distribution for type two events.
`h1.fn`	A (vectorized) function. The offspring density function for type one events.
`h2.fn`	A (vectorized) function. The offspring density function for type two events.
`p.h1`	A numeric vector. The paramters of the offspring density for type one events.
`p.h2`	A numeric vector. The paramters of the offspring density for type two events.
`eta11`	A numeric scalar. The self-exciting branching ratio for type one events.
`eta12`	A numeric scalar. The cross-exciting branching ratio for type one events due to the effects of a type two event.
`eta21`	A numeric scalar. The cross-exciting branching ratio for type two events due to the effects of a type one event.
`eta22`	A numeric scalar. The self-exciting branching ratio for type two events.
`B`	A numeric scalar. Tuning parameter
`B0`	A numeric scalar. Tuning parameter
`pnp1`	A numeric square matrix. The joint last immigrant probabilities.
`max.h1`	A numeric scalar. The maximum value of the offspring density for type one events.
`max.h2`	A numeric scalar. The maximum value of the offspring density for for type two events.

Value

A numeric matrix that contains the simulated event times from censoring time cens up until cens.tilde and the corresponding event types.

Author(s)

Tom Stindl <[email protected]> Feng Chen <[email protected]>

Examples

  
    ## Magnitude 5.5 or greater earthquakes over the 25 year period from 
    ## 01/01/1991 to 31/12/2015.
    data(fivaqks); 
    near.fiji <- grep("Fiji", fivaqks$place)
    near.vanuatu <- grep("Vanuatu", fivaqks$place)
    t.beg <- strptime("1991-01-01 00:00:00", "%Y-%m-%d %H:%M:%S", tz = "UTC")
    t.end <- strptime("2015-12-31 23:59:59", "%Y-%m-%d %H:%M:%S", tz = "UTC")
    t0 <- 0
    t1 <- as.numeric(t.end - t.beg)
    tms <- strptime(fivaqks$time, "%Y-%m-%dT%H:%M:%OSZ", tz = "UTC")
    ts <- as.numeric(tms[-1] - t.beg)
    ts <- c(as.numeric(tms[1] - t.beg)/24, ts)
    ts.fi <- ts[near.fiji]; ts.fi <- ts.fi[ts.fi >= 0 & ts.fi <= t1]
    ts.va <- ts[near.vanuatu]; ts.va <- ts.va[ts.va >=0 & ts.va <= t1]
    ts.c <- c(ts.fi, ts.va)
    z.c <- c(rep(1, times = length(ts.fi)), rep(2, times = length(ts.va)))
    o <- order(ts.c)
    data <- cbind(ts.c[o], z.c[o])
    # simulate future event time based on MLE fitted Rhawkes model
    N <- 100; i <- 0;
    data.pred <- replicate(N, 
                          {cat(i<<-i+1,'\n'); 
                          simpred(data = data,
                                    par = c(0.488, 20.10, 0.347, 9.53, 
                                            461, 720,
                                            0.472, 0.293, 0.399, -0.0774), 
                                            cens = t1, cens.tilde = t1 + 1000)
                                            })
  
## Magnitude 5.5 or greater earthquakes over the 25 year period from 
    ## 01/01/1991 to 31/12/2015.
    data(fivaqks); 
    near.fiji <- grep("Fiji", fivaqks$place)
    near.vanuatu <- grep("Vanuatu", fivaqks$place)
    t.beg <- strptime("1991-01-01 00:00:00", "%Y-%m-%d %H:%M:%S", tz = "UTC")
    t.end <- strptime("2015-12-31 23:59:59", "%Y-%m-%d %H:%M:%S", tz = "UTC")
    t0 <- 0
    t1 <- as.numeric(t.end - t.beg)
    tms <- strptime(fivaqks$time, "%Y-%m-%dT%H:%M:%OSZ", tz = "UTC")
    ts <- as.numeric(tms[-1] - t.beg)
    ts <- c(as.numeric(tms[1] - t.beg)/24, ts)
    ts.fi <- ts[near.fiji]; ts.fi <- ts.fi[ts.fi >= 0 & ts.fi <= t1]
    ts.va <- ts[near.vanuatu]; ts.va <- ts.va[ts.va >=0 & ts.va <= t1]
    ts.c <- c(ts.fi, ts.va)
    z.c <- c(rep(1, times = length(ts.fi)), rep(2, times = length(ts.va)))
    o <- order(ts.c)
    data <- cbind(ts.c[o], z.c[o])
    # simulate future event time based on MLE fitted Rhawkes model
    N <- 100; i <- 0;
    data.pred <- replicate(N, 
                          {cat(i<<-i+1,'\n'); 
                          simpred(data = data,
                                    par = c(0.488, 20.10, 0.347, 9.53, 
                                            461, 720,
                                            0.472, 0.293, 0.399, -0.0774), 
                                            cens = t1, cens.tilde = t1 + 1000)
                                            })

Minus loglikelihood of an (bivariate) MRHawkes model with truncated most recent immigrant probabilities

Description

Calculates the minus loglikelihood of an (bivariate) MRHawkes model with given immigration hazard functions $\mu$ , offspring density functions $h$ and bracnhing ratios $\eta$ for event times and event types data on interval $[0,cens]$ with truncated most recent immigrant probabilities looking back at only the previous B event times.

Usage

TmllMRH(data, cens, par, B = 25,
        h1.fn = function(x, p) 1 / p * exp( - x / p),
        h2.fn = function(x, p) 1 / p * exp( - x / p),
        mu1.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))
        },
        mu2.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))
        },
        H1.fn = function(x, p) pexp(x, rate = 1 / p),
        H2.fn = function(x, p) pexp(x, rate = 1 / p),
        Mu1.fn = function(x, p){
          - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, 
                    log.p = TRUE)
        },
        Mu2.fn = function(x, p){
          - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, 
                     log.p = TRUE)
        })
TmllMRH(data, cens, par, B = 25,
        h1.fn = function(x, p) 1 / p * exp( - x / p),
        h2.fn = function(x, p) 1 / p * exp( - x / p),
        mu1.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))
        },
        mu2.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))
        },
        H1.fn = function(x, p) pexp(x, rate = 1 / p),
        H2.fn = function(x, p) pexp(x, rate = 1 / p),
        Mu1.fn = function(x, p){
          - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, 
                    log.p = TRUE)
        },
        Mu2.fn = function(x, p){
          - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, 
                     log.p = TRUE)
        })

Arguments

`data`	A two column matrix. The first column contains the event times sorted in ascending order. The second column contains the corresponding event type with the label one or two.
`cens`	A scalar. The censoring time.
`par`	A numeric vector. Contains the ten parameters of the model, in order of the immigration parameters $\mu(.)$ for the two renewal distributions, the two offspring parameters $h(.)$ and lastly the four branching ratios $\eta$ .
`B`	A scalar. The number of previous events that are considered to be possible last immigrant arrivals.
`h1.fn`	A (vectorized) function. The offspring density function for type one events.
`h2.fn`	A (vectorized) function. The offspring density function for type two events.
`mu1.fn`	A (vectorized) function. The immigration hazard function for events of type one.
`mu2.fn`	A (vectorized) function. The immigration hazard function for events of type two.
`H1.fn`	A (vectorized) function. Its value at `t` gives the integral of the offspring density function from 0 to `t` for type one events.
`H2.fn`	A (vectorized) function. Its value at `t` gives the integral of the offspring density function from 0 to `t` for type two events.
`Mu1.fn`	A (vectorized) function. Its value at `t` gives the integral of the immigrant hazard function from 0 to `t` for type one events.
`Mu2.fn`	A (vectorized) function. Its value at `t` gives the integral of the immigrant hazard function from 0 to `t` for type two events.

Value

Value of the negative log-liklihood.

Author(s)

Tom Stindl <[email protected]> Feng Chen <[email protected]>

Examples

  ## Magnitude 5.5 or greater earthquakes over the 25 year period from 
  ## 01/01/1991 to 31/12/2015.
  data(fivaqks); 
  near.fiji <- grep("Fiji", fivaqks$place)
  near.vanuatu <- grep("Vanuatu", fivaqks$place)
  t.beg <- strptime("1991-01-01 00:00:00", "%Y-%m-%d %H:%M:%S", tz = "UTC")
  t.end <- strptime("2015-12-31 23:59:59", "%Y-%m-%d %H:%M:%S", tz = "UTC")
  t0 <- 0
  t1 <- as.numeric(t.end - t.beg)
  tms <- strptime(fivaqks$time, "%Y-%m-%dT%H:%M:%OSZ", tz = "UTC")
  ts <- as.numeric(tms[-1] - t.beg)
  ts <- c(as.numeric(tms[1] - t.beg)/24, ts)
  ts.fi <- ts[near.fiji]; ts.fi <- ts.fi[ts.fi >= 0 & ts.fi <= t1]
  ts.va <- ts[near.vanuatu]; ts.va <- ts.va[ts.va >=0 & ts.va <= t1]
  ts.c <- c(ts.fi, ts.va)
  z.c <- c(rep(1, times = length(ts.fi)), rep(2, times = length(ts.va)))
  o <- order(ts.c)
  data <- cbind(ts.c[o], z.c[o])
  
  ## calculate the minus loglikelihood of an (bivariate) MRHawkes with some 
  ## parameters the default hazard functions and density functions are Weibull 
  ## and exponential respectivley
  TmllMRH(data, cens = t1, par = c(0.488, 20.10, 0.347, 9.53, 461, 720,
                                  0.472, 0.293, 0.399, -0.0774), B = 200)
                                  
  ## 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.488, 20.10, 0.347, 9.53, 461, 720,
                                  0.472, 0.293, 0.399, -0.0774), 
                             TmllMRH, data = data, cens = t1, B = 200,
                             control = list(maxit = 5000, trace = TRUE),
                             hessian = TRUE)
    )
    ## point estimate by MLE
    est$par
    ## standard error estimates:
    diag(solve(est$hessian))^0.5
  
## Magnitude 5.5 or greater earthquakes over the 25 year period from 
  ## 01/01/1991 to 31/12/2015.
  data(fivaqks); 
  near.fiji <- grep("Fiji", fivaqks$place)
  near.vanuatu <- grep("Vanuatu", fivaqks$place)
  t.beg <- strptime("1991-01-01 00:00:00", "%Y-%m-%d %H:%M:%S", tz = "UTC")
  t.end <- strptime("2015-12-31 23:59:59", "%Y-%m-%d %H:%M:%S", tz = "UTC")
  t0 <- 0
  t1 <- as.numeric(t.end - t.beg)
  tms <- strptime(fivaqks$time, "%Y-%m-%dT%H:%M:%OSZ", tz = "UTC")
  ts <- as.numeric(tms[-1] - t.beg)
  ts <- c(as.numeric(tms[1] - t.beg)/24, ts)
  ts.fi <- ts[near.fiji]; ts.fi <- ts.fi[ts.fi >= 0 & ts.fi <= t1]
  ts.va <- ts[near.vanuatu]; ts.va <- ts.va[ts.va >=0 & ts.va <= t1]
  ts.c <- c(ts.fi, ts.va)
  z.c <- c(rep(1, times = length(ts.fi)), rep(2, times = length(ts.va)))
  o <- order(ts.c)
  data <- cbind(ts.c[o], z.c[o])
  
  ## calculate the minus loglikelihood of an (bivariate) MRHawkes with some 
  ## parameters the default hazard functions and density functions are Weibull 
  ## and exponential respectivley
  TmllMRH(data, cens = t1, par = c(0.488, 20.10, 0.347, 9.53, 461, 720,
                                  0.472, 0.293, 0.399, -0.0774), B = 200)
                                  
  ## 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.488, 20.10, 0.347, 9.53, 461, 720,
                                  0.472, 0.293, 0.399, -0.0774), 
                             TmllMRH, data = data, cens = t1, B = 200,
                             control = list(maxit = 5000, trace = TRUE),
                             hessian = TRUE)
    )
    ## point estimate by MLE
    est$par
    ## standard error estimates:
    diag(solve(est$hessian))^0.5

Minus loglikelihood of an (bivariate) MRHawkes model with Universal residuals

Description

Calculates the minus loglikelihood of an (bivariate) RHawkes model with given immigration hazard functions $\mu$ , common offspring density functions $h$ and bracnhing ratios $\eta$ for event times and event types data on interval $[0,cens]$ . The same as mllMRH although this version also returns the Universal residuals for goodness-of-fit assessment of the event types.

Usage

  typeRes(data, cens, par, U = runif(length(data[,1])),
          h1.fn = function(x, p) 1 / p * exp( - x / p),
          h2.fn = function(x, p) 1 / p * exp( - x / p),
          mu1.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))
          },
          mu2.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))
          },
          H1.fn = function(x, p) pexp(x, rate = 1 / p),
          H2.fn = function(x, p) pexp(x, rate = 1 / p),
          Mu1.fn = function(x, p){
            - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, 
                       log.p = TRUE)
          },
          Mu2.fn = function(x, p){
            - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, 
                       log.p = TRUE)
          })
typeRes(data, cens, par, U = runif(length(data[,1])),
          h1.fn = function(x, p) 1 / p * exp( - x / p),
          h2.fn = function(x, p) 1 / p * exp( - x / p),
          mu1.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))
          },
          mu2.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))
          },
          H1.fn = function(x, p) pexp(x, rate = 1 / p),
          H2.fn = function(x, p) pexp(x, rate = 1 / p),
          Mu1.fn = function(x, p){
            - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, 
                       log.p = TRUE)
          },
          Mu2.fn = function(x, p){
            - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, 
                       log.p = TRUE)
          })

Arguments

`data`	A two column matrix. The first column contains the event times sorted in ascending order. The second column contains the corresponding event type with the label one or two.
`cens`	A scalar. The censoring time.
`par`	A numeric vector. Contains the ten parameters of the model, in order of the immigration parameters $\mu(.)$ for the two renewal distributions, the two offspring parameters $h(.)$ and lastly the four branching ratios $\eta$ .
`U`	A numeric vector. Contains auxillary uniform random varables on the unit interval.
`h1.fn`	A (vectorized) function. The offspring density function for type one events.
`h2.fn`	A (vectorized) function. The offspring density function for type two events.
`mu1.fn`	A (vectorized) function. The immigration hazard function for events of type one.
`mu2.fn`	A (vectorized) function. The immigration hazard function for events of type two.
`H1.fn`	A (vectorized) function. Its value at `t` gives the integral of the offspring density function from 0 to `t` for type one events.
`H2.fn`	A (vectorized) function. Its value at `t` gives the integral of the offspring density function from 0 to `t` for type two events.
`Mu1.fn`	A (vectorized) function. Its value at `t` gives the integral of the immigrant hazard function from 0 to `t` for type one events.
`Mu2.fn`	A (vectorized) function. Its value at `t` gives the integral of the immigrant hazard function from 0 to `t` for type two events.

Details

Calculate the MRHawkes point process Universal residuals

Value

`mll`	minus log-likelihood
`V`	Universal residuals of observed event types

Author(s)

Tom Stindl <[email protected]> Feng Chen <[email protected]>

Examples


  n <- 1000
  data <- cbind(sort(runif(n,0,1000)), 
                sample(1:2, size = n, replace = TRUE))
  tmp <- typeRes(data = data, cens = 1001, 
                 par = c(1,1,1,1,1,1,0.5,0.2,0.2,0.3))              
  pp <- ppoints(n)
  par(mfrow=c(1,2))
  plot(quantile(tmp$V,prob=pp),pp,type="l",
       main="Uniform QQ plot",
       xlab="Sample quantiles",ylab="Theoretical quantiles")
  abline(a = 0, b = 1, col = 2)
  a <- acf(tmp$V, main = "ACF Plot")
  ks.test(tmp$V,"punif")
  Box.test(tmp$V,lag=tail(a$lag,1))
  
n <- 1000
  data <- cbind(sort(runif(n,0,1000)), 
                sample(1:2, size = n, replace = TRUE))
  tmp <- typeRes(data = data, cens = 1001, 
                 par = c(1,1,1,1,1,1,0.5,0.2,0.2,0.3))              
  pp <- ppoints(n)
  par(mfrow=c(1,2))
  plot(quantile(tmp$V,prob=pp),pp,type="l",
       main="Uniform QQ plot",
       xlab="Sample quantiles",ylab="Theoretical quantiles")
  abline(a = 0, b = 1, col = 2)
  a <- acf(tmp$V, main = "ACF Plot")
  ks.test(tmp$V,"punif")
  Box.test(tmp$V,lag=tail(a$lag,1))

Package 'MRHawkes'

Help Index

Multivariate Renewal Hawkes Process

Description

Details

Author(s)

References

Fiji and Vanuatu Earthquake Data

Description

Usage

Format

Source

Examples

Minus loglikelihood of an (bivariate) MRHawkes model

Description

Usage

Arguments

Value

Author(s)

Examples

Minus loglikelihood of an (bivariate) MRHawkes model with most recent immigrant probabilities

Description

Usage

Arguments

Value

Author(s)

See Also

Examples

Minus loglikelihood of an (bivariate) MRHawkes model with Rosenblatt residuals

Description

Usage

Arguments

Details

Value

Author(s)

See Also

Examples

MRHawkes (bivariate) predictive density function

Description

Usage

Arguments

Value

Author(s)

Examples

Simulate an (bivariate) renewal Hawkes (MRHawkes) process

Description

Usage

Arguments

Details

Value

Author(s)

Examples

Simulate a (bivariate) non-stationary multivariate Hawkes process (NSMHP)

Description

Usage

Arguments

Details

Value

Author(s)

Examples

Simulate a fitted (bivariate) MRHawkes process model

Description

Usage

Arguments

Value

Author(s)

Examples

Minus loglikelihood of an (bivariate) MRHawkes model with truncated most recent immigrant probabilities

Description

Usage

Arguments

Value

Author(s)

Examples

Minus loglikelihood of an (bivariate) MRHawkes model with Universal residuals

Description

Usage

Arguments

Details

Value