R/rETASlp.R

In stopp: Spatio-Temporal Point Pattern Methods, Model Fitting, Diagnostics, Simulation, Local Tests

#' Simulation of a spatio-temporal ETAS (Epidemic Type Aftershock Sequence) model
#' on a linear network
#'
#'
#' @description
#' This function simulates a spatio-temporal ETAS
#' (Epidemic Type Aftershock Sequence) process on a linear network
#'  as a \code{stpm} object.
#'  
#' It is firstly introduced and employed for simulation studies in D'Angelo et al. (2021).
#'
#' It follows the generating scheme for simulating a pattern from an
#' Epidemic Type Aftershocks-Sequences (ETAS) process
#' (Ogata and Katsura 1988) with conditional intensity function (CIF) as in
#' Adelfio and Chiodi (2020), adapted for the space location of events
#' to be constrained on a linear network.
#'
#' The simulation on the network is guaranteed by the homogeneous spatial
#' Poisson processes being generated on the network.
#'
#'
#' @details
#' The CIF of an ETAS
#' process as in Adelfio and Chiodi (2020) can be written as \deqn{
#'   \lambda_{\theta}(t,\textbf{u}|\mathcal{H}_t)=\mu f(\textbf{u})+\sum_{t_j<t} \frac{\kappa_0 \exp(\eta_j)}{(t-t_j+c)^p} \{    (\textbf{u}-\textbf{u}_j)^2+d   \}^{-q}  ,
#'} where
#'
#' \eqn{\mathcal{H}_t} is the past history of the process up to time
#' \eqn{t}
#'
#' \eqn{\mu} is the large-scale general intensity
#'
#' \eqn{f(\textbf{u})} is
#' the spatial density
#'
#' \eqn{\eta_j=\boldsymbol{\beta}' \textbf{Z}_j} is a linear predictor
#'
#' \eqn{\textbf{Z}_j} the external known covariate vector, including the
#' magnitude
#'
#' \eqn{\boldsymbol{\theta}= (\mu, \kappa_0, c, p, d, q, \boldsymbol{\beta})}
#' are the parameters to be estimated
#'
#' \eqn{\kappa_0} is a
#' normalising constant
#'
#' \eqn{c} and \eqn{p} are characteristic parameters of the
#' seismic activity of the given region,
#'
#'  and \eqn{d} and \eqn{q} are two parameters
#' related to the spatial influence of the mainshock
#'
#' In the usual ETAS
#' model for seismic analyses, the only external covariate represents the magnitude,
#' \eqn{\boldsymbol{\beta}=\alpha}, as
#' \eqn{\eta_j = \boldsymbol{\beta}' \textbf{Z}_j = \alpha (m_j-m_0)}, where
#' \eqn{m_j} is the magnitude of the \eqn{j^{th}} event and \eqn{m_0} the threshold
#' magnitude, that is, the lower bound for which earthquakes with higher
#' values of magnitude are surely recorded in the catalogue.
#'
#'
#' @param pars A vector of parameters of the ETAS model to be simulated.
#'  See the 'Details' section.
#' @param betacov Numerical array. Parameters of the covariates ETAS model
#' @param m0 Parameter for the background general intensity of the ETAS model.
#' In the common seismic analyses it represents the threshold
#' magnitude.
#' @param b 1.0789
#' @param t.lag 200
#' @param tmin Minimum value of time.
#' @param covsim Default \code{FALSE}
#' @param L linear network
#' @param all.marks Logical value indicating whether to store
#' all the simulation information as marks in the \code{stlpm} object.
#' If \code{FALSE} (default option) only the magnitude is returned.
#' @return A \code{stlpm} object
#' @export
#'
#' @author Nicoletta D'Angelo and Marcello Chiodi
#'
#'
#' @examples
#'
#' set.seed(95)
#' X <- rETASlp(pars = c(0.1293688525, 0.003696, 0.013362, 1.2,0.424466,  1.164793),
#'      L = chicagonet)
#'
#' @references
#' Adelfio, G., and Chiodi, M. (2021). Including covariates in a space-time point process with application to seismicity. Statistical Methods & Applications, 30(3), 947-971.
#'
#' D’Angelo, N., Adelfio, G., and Mateu, J. (2021). Assessing local differences between the spatio-temporal second-order structure of two point patterns occurring on the same linear network. Spatial Statistics, 45, 100534.
#'
#' Ogata, Y., and Katsura, K. (1988). Likelihood analysis of spatial inhomogeneity for marked point patterns. Annals of the Institute of Statistical Mathematics, 40(1), 29-39.
#'

rETASlp <- function(pars = NULL
                        , betacov = 0.39
                        , m0 = 2.5
                        , b = 1.0789
                        , tmin = 0
                        , t.lag = 200
                        , covsim = FALSE
                        , L, all.marks = FALSE){
  if(is.null(pars)) stop("Please provide some parameters")
  
  cat = NULL

  xmin = L$window$xrange[1]
  xmax = L$window$xrange[2]
  ymin = L$window$yrange[1]
  ymax = L$window$yrange[2]

  mu	= pars[1]
  k0  	= pars[2]
  c   	= pars[3]
  p   	= pars[4]
  gamma   =0
  d       = pars[5]
  q       = pars[6]
  ncov    = length(betacov)
  beta    = log(10) * b
  mm      = (ymax - ymin) / (xmax - xmin)
  qq      = ymin - mm * xmin
  cat.pois = NULL
  cat.sim = NULL
  cat.new = NULL
  tmax    = tmin + t.lag
  ak      = k0 * c ^ (1-p) / (p-1)
  sk      = (pi * d ^ (1 - q)) / (q - 1)
  muback  = mu * (tmax - tmin)
  n0      = rpois(1, muback)
  nstart  = n0

  if (nstart > 0) {

    lgen       = array(0, nstart)
    ind     = array(0, nstart)
    father  = array(0, nstart)

    tback   = runif(n0, tmin, tmax)
    rand <- spatstat.random::datagen.runifpointOnLines(n0, spatstat.geom::as.psp(L))
    xback = rand$x
    yback = rand$y

    mback   = m0 + rexp(n0, rate = beta)
    zback   = array(0, n0)
    if (ncov == 1) {
      cov2back = array(0, n0)
    }
    else{
      if (covsim)
      {
        cov2back = rnorm(n0) ^ 2
      }
      else
      {
        cov2back = abs(yback - mm * xback - qq) / (sqrt(1 + mm * mm))
      }
    }


    cat.new = cbind(tback, xback, yback, mback, zback, cov2back)
    cat.pois = cat.new

    cat.new             = cbind(cat.new, lgen, ind, father)
    colnames(cat.new)   = c("time", "long", "lat", "magn1", "z", "cov2", "lgen", "ind", "father")
    cat.new             = as.data.frame(cat.new)
    cat.new = cat.new[order(cat.new$time), ]
    i = 0
    while(!prod(cat.new$ind)){
      i   = i + 1
      cat.new$ind[i]  = TRUE
      pred            = (cat.new$magn1[i] - m0) * betacov[1]
      if (ncov > 1) pred = pred + cat.new$cov2[i] * betacov[2]
      nexpected       = ak * sk * exp(pred)
      ni              = rpois(1, nexpected)
      if (ni > 0) {
        xy              = norm2_etas(ni, d, q, cat.new$long[i], cat.new$lat[i])
        t             = c * runif(ni) ^ ( - 1 / (p - 1)) - c + cat.new$time[i]
        ind.txy = (t > tmin)&(t < tmax)&(xy[,1] > xmin)&(xy[,1] < xmax)&(xy[,2] > ymin)&(xy[,2] < ymax)
        ntxy = sum(ind.txy)
        if (ntxy > 0) {
          m1      = m0 + rexp(ntxy, rate = beta)
          z1      = array(0, ntxy)
          cov2    = array(0, ntxy)
          lgen    = array(1, ntxy) + cat.new$lgen[i]
          ind     = array(0, ntxy)
          father  = array(i, ntxy)
          if (ncov > 1)
          {

            if (covsim)
            {
              cov2 = rnorm(ntxy) ^ 2
            }
            else
            {
              xo      = xy[ind.txy,1]
              yo      = xy[ind.txy,2]
              cov2    = abs(yo - mm * xo - qq) / (sqrt(1 + mm * mm))
            }


          }
          cat.son = cbind(t[ind.txy], xy[ind.txy, 1], xy[ind.txy, 2], m1, z1, cov2, lgen, ind, father)
          colnames(cat.son)   = c("time", "long", "lat", "magn1", "z", "cov2", "lgen", "ind", "father")
          cat.son = as.data.frame(cat.son)
          cat.new = rbind(cat.new, cat.son)
          cat.new = cat.new[order(cat.new$time), ]
          colnames(cat.new)   = colnames(cat.son) 

        }
      }
    }
  }
  nson = nrow(cat.new) - n0
  cat.new[, 1:3] <- cat.new[c(2, 3, 1)]
  if(all.marks){
    return(stpm(cat.new, names = colnames(cat.new)[-c(1:3)], L = L))
  } else {
    return(stpm(cat.new[, 1:4], names = "Magnitude", L = L))
  }
   
}