simpred {MRHawkes}

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

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 <[e-mail]> Feng Chen <[e-mail]>

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