PtProcess (version 3.3-10)

etas_spatial: Template Function for Spatial ETAS

Description

This is a template function for a spatial ETAS model. The spatial component is based on the shape of the bivariate normal density function. This function is currently in development and may change, see Warnings below.

Usage

etas_normal0(data, evalpts, params, fixedparams, TT=NA,
             tplus=FALSE)

Arguments

data

a data frame containing the event history, where each row represents one event. There must be columns named "time", usually the number of days from some origin; "magnitude" which is the event magnitude less the magnitude threshol

evalpts

a matrix or data.frame containing columns named "time", "magnitude", "latitude", and "longitude".

params

vector of parameter values in the following order: $(\mu, A, \alpha, c, p, d_x, d_y, \beta)$.

fixedparams

vector of values defining the spatial boundaries (degrees) of the region, assumed to be rectangular: $(X_1, X_2, Y_1, Y_2)$.

TT

vector of length 2, being the time interval over which the integral of the ground intensity function is to be evaluated.

tplus

logical, $\lambda(t,x,y|{\cal H}_t)$ is evaluated as $\lambda(t^+,x,y|{\cal H}_t)$ if TRUE, else $\lambda(t^-,x,y|{\cal H}_t)$.

Value

• The function has two forms of usage as described in gif. The first usage returns a vector containing the values of $\lambda(t,x,y)$ evaluated at the specified points. In the second usage, it returns the value of the space-time integral.

Warnings

This is function is still in development and may often change. The code may well be inefficient, though the current purpose is to keep it transparent and hence easier to change.

When this function is used to build a model object (mpp), the model object will work correctly with the generic functions summary, residuals, neglogLik and logLik, but will not yet work with simulate and plot.

Details

Let $$f(t) = (1+t/c)^{-p}$$ and $$g(x,y) = \exp\left(-\frac{x^2}{2\sigma_x^2} - \frac{y^2}{2\sigma_y^2} \right) \,,$$ where $\sigma_x^2 = d_x e^{\alpha(M_i-M_0)}$ and $\sigma_y^2 = d_y e^{\alpha(M_i-M_0)}$. Then the conditional intensity function is $$\lambda(t,x,y|{\cal H}_t) = \mu + A \sum_{i: t_i < t} e^{\beta(M_i-M_0)} g(x-x_i, y-y_i) f(t-t_i) \,.$$ Note that the parameters $\mu$ and $A$ are not the same as those in etas_gif. That function is the spatial integral of the above, and so, for example, $\widetilde{\mu} = (X_2-X_1)(Y_2-Y_1)\mu$ where $\widetilde{\mu}$ is the $\mu$ in etas_gif. Ogata & Zhuang (2006) have discussed this and other formulations of the spatial ETAS models.

This can be thought of as a template function, i.e. a function that can be used to define the required spatial intensity function. It is probably most efficient to embed the boundaries of the region under analysis into the function. Say we want the scaling for latitude and longitude to be the same (i.e. $d=d_x=d_y$), $\beta=0$, and the analysed region is contained within $89^\circ$E--$105^\circ$E, and $5^\circ$S--$16^\circ$N. This is then the same as case (5) in Ogata & Zhuang (2006) with their $S_i$ being the identity matrix. The conditional intensity function would then be defined as follows. etas_tmp <- function(data, evalpts, params, TT=NA, tplus=FALSE) etas_normal0(data, evalpts, params=c(params, params[6], 0), fixedparams=c(89, 105, -5, 16), TT=TT, tplus=tplus) The vector params of the new function etas_tmp is $(\mu, A, \alpha, c, p, d)$.

References

Cited references are listed on the PtProcess manual page.

See Also

General details about the structure of ground intensity functions are given in the topic gif.

Examples

data(Phuket)
#    magnitudes are rounded to 1 dp
Phuket$magnitude <- Phuket$magnitude - 4.95

TT <- c(0, 1827)
params <- c(0.0001, 73, 1.3, 0.01, 1.0, 0.01)

etas_tmp <- function(data, evalpts, params, TT=NA, tplus=FALSE){
    # etas_normal0 params = c(mu, A, alpha, c, p, dx, dy, beta)
    # etas_tmp params = c(mu, A, alpha, c, p, d)
    # want: d=dx=dy and alpha=beta
    etas_normal0(data, evalpts, params=c(params, params[6], 0),
                 fixedparams=c(89, 105, -5, 16),
                 TT=TT, tplus=tplus)
}

x <- mpp(data=Phuket,
         gif=etas_tmp,
         marks=list(NULL, NULL),
         params=params,
         gmap=expression(params),
         mmap=NULL,
         TT=TT)

allmap <- function(y, p){
    y$params <- exp(p)
    return(y)
}

#    Note: the "Not run" blocks below are not run during package checks
#    as the makeSOCKcluster definition is specific to my network,
#    modify accordingly if you want parallel processing.

cl <- NULL
if (require(snow)){
    cl <- makeSOCKcluster(c("horoeka.localdomain", "horoeka.localdomain",
                            "kowhai.localdomain", "kowhai.localdomain",
                            "localhost", "localhost"))
    clusterExport(cl, c("etas_normal0"))
}

initial <- log(params)
#  set iterlim larger to converge satisfactorily
#  needs about 47 iterations
z <- nlm(neglogLik, initial, object=x, pmap=allmap,
         print.level=2, iterlim=1, SNOWcluster=cl)

if (!is.null(cl)){
    stopCluster(cl)
    rm(cl)
}

x <- allmap(x, z$estimate)
print(z$code)

print(x$params)
#  MLE = (1.139316e-04 7.258019e+01 1.304519e+00 8.839539e-03
#         9.999858e-01 9.299687e-03)

print(logLik(x))
# -2954.28

#   The integral term should be equal to the number
#   of events in TT (1248)
print(sum((x$data$time > TT[1]) & (x$data$time < TT[2])))
print(etas_tmp(x$data, NULL, x$params, TT=TT))