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:

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 <t.stindl@unsw.edu.au> Feng Chen <feng.chen@unsw.edu.au>

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.

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)
        })

Arguments

Value

Value of the negative log-liklihood.

Author(s)

Tom Stindl <t.stindl@unsw.edu.au> Feng Chen <feng.chen@unsw.edu.au>

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
  

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)

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)

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)
       })

Arguments

Value

Value of the negative log-liklihood.

Author(s)

Tom Stindl <t.stindl@unsw.edu.au> Feng Chen <feng.chen@unsw.edu.au>

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
  

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)
       })

Arguments

Value

Author(s)

Tom Stindl <t.stindl@unsw.edu.au> Feng Chen <feng.chen@unsw.edu.au>

See Also

mllMRH

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)))
  

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)
       })

Arguments

Details

Calculate the MRHawkes point process Rosenblatt residuals

Value

Author(s)

Tom Stindl <t.stindl@unsw.edu.au> Feng Chen <feng.chen@unsw.edu.au>

See Also

mllMRH

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))
  

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)
        })

Arguments

Value

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

Author(s)

Tom Stindl <t.stindl@unsw.edu.au> Feng Chen <feng.chen@unsw.edu.au>

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")
  

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)

Arguments

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 <t.stindl@unsw.edu.au> Feng Chen <feng.chen@unsw.edu.au>

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)
                    })

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))

Arguments

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

Author(s)

Tom Stindl <t.stindl@unsw.edu.au> Feng Chen <feng.chen@unsw.edu.au>

Examples

  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)

Arguments

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 <t.stindl@unsw.edu.au> Feng Chen <feng.chen@unsw.edu.au>

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)
                                            })
  

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)
          })

Arguments

Details

Calculate the MRHawkes point process Universal residuals

Value

Author(s)

Tom Stindl <t.stindl@unsw.edu.au> Feng Chen <feng.chen@unsw.edu.au>

See Also

mllMRH

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))