Statistical Analysis of Point Processes

Description

R functions for statistical analysis of point processes

Details

This package provides functions for statistical analysis of series of events and seismicity.

For overview of point process models, 'Statistical Analysis of Point Processes with R' is available in the package vignette using the vignette function (e.g., vignette("SAPP")).

References

Ogata, Y., Katsura, K. and Zhuang, J. (2006) Computer Science Monographs, No.32, TIMSAC84: STATISTICAL ANALYSIS OF SERIES OF EVENTS (TIMSAC84-SASE) VERSION 2. The Institute of Statistical Mathematics. https://www.ism.ac.jp/editsec/csm/

Ogata, Y. (2006) Computer Science Monographs, No.33, Statistical Analysis of Seismicity - updated version (SASeies2006). The Institute of Statistical Mathematics. https://www.ism.ac.jp/editsec/csm/

The Occurrence Times Data

Description

This data consists of the occurrence times of 627 blastings at a certain stoneyard with a very small portion of microearthquakes during a past 4600 days.

Usage

data(Brastings)

Format

A numeric vector of length 627.

Source

Ogata, Y., Katsura, K. and Zhuang, J. (2006) Computer Science Monographs, No.32, TIMSAC84: Statistical Analysis of Series of events (TIMSAC84-SASE) Version 2. The Institute of Statistical Mathematics.

The Point Process Data

Description

The point process test data for linsim and linlin.

Usage

data(PProcess)

Format

A numeric vector of length 72.

Source

Ogata, Y., Katsura, K. and Zhuang, J. (2006) Computer Science Monographs, No.32, TIMSAC84: STATISTICAL ANALYSIS OF SERIES OF EVENTS (TIMSAC84-SASE) VERSION 2. The Institute of Statistical Mathematics.

Poisson Data

Description

Poisson test data for ptspec.

Usage

data(PoissonData)

Format

A numeric vector of length 2553.

Source

Ogata, Y., Katsura, K. and Zhuang, J. (2006) Computer Science Monographs, No.32, TIMSAC84: STATISTICAL ANALYSIS OF SERIES OF EVENTS (TIMSAC84-SASE) VERSION 2. The Institute of Statistical Mathematics.

Self-Exciting Point Process Data

Description

Self-exciting point process test data for linlin.

Usage

data(SelfExcit)

Format

A numeric vector of length 99.

Source

Ogata, Y., Katsura, K. and Zhuang, J. (2006) Computer Science Monographs, No.32, TIMSAC84: STATISTICAL ANALYSIS OF SERIES OF EVENTS (TIMSAC84-SASE) VERSION 2. The Institute of Statistical Mathematics.

Maximum Likelihood Estimates of Intensity Rates

Description

Compute the maximum likelihood estimates of intensity rates of either exponential polynomial or exponential Fourier series of non-stationary Poisson process models.

Usage

eptren(data, mag = NULL, threshold = 0.0, nparam, nsub, cycle = 0,
       tmpfile = NULL, nlmax = 1000, plot = TRUE)

Arguments

Details

This function computes the maximum likelihood estimates (MLEs) of the coefficients A_1, A_2,\ldots A_n is an exponential polynomial

f(t) = exp(A_1 + A_2t + A_3t^2 + ... )

or A_1, A_2, B_2, ..., A_n, B_n in a Poisson process model with an intensity taking the form of an exponential Fourier series

f(t) = exp\{ A_1 + A_2cos(2\pi t/p) + B_2sin(2\pi t/p) + A_3cos(4\pi t/p) + B_3sin(4\pi t/p) +... \}

which represents the time varying rate of occurrence (intensity function) of earthquakes in a region.

These two models belong to the family of non-stationary Poisson process. The optimal order n can be determined by minimize the value of the Akaike Information Criterion (AIC).

Value

References

Ogata, Y., Katsura, K. and Zhuang, J. (2006) Computer Science Monographs, No.32, TIMSAC84: STATISTICAL ANALYSIS OF SERIES OF EVENTS (TIMSAC84-SASE) VERSION 2. The Institute of Statistical Mathematics.

Ogata, Y. (2006) Computer Science Monographs, No.33, Statistical Analysis of Seismicity - updated version (SASeies2006). The Institute of Statistical Mathematics.

Examples

## The Occurrence Times Data of 627 Blastings
data(Brastings)

# exponential polynomial trend fitting
eptren(Brastings, nparam = 10, nsub = 1000)

# exponential Fourier series fitting
eptren(Brastings, nparam = 10, nsub = 1000, cycle = 1)

## Poisson Process data
data(PoissonData)

# exponential polynomial trend fitting
eptren(PoissonData, nparam = 10, nsub = 1000)

# exponential Fourier series fitting
eptren(PoissonData, nparam = 10, nsub = 1000, cycle = 1)

## The aftershock data of 26th July 2003 earthquake of M6.2
data(main2003JUL26)
x <- main2003JUL26

# exponential polynomial trend fitting
eptren(x$time, mag = x$magnitude, nparam = 10, nsub = 1000)

# exponential Fourier series fitting
eptren(x$time, mag = x$magnitude, nparam = 10, nsub = 1000, cycle = 1)   

Residual Point Process of the ETAS Model

Description

Compute the residual data using the ETAS model with MLEs.

Usage

etarpp(time, mag, threshold = 0.0, reference = 0.0, parami, zts = 0.0, tstart,
       zte, ztend = NULL, plot = TRUE)

etarpp2(etas, threshold = 0.0, reference = 0.0, parami, zts = 0.0, tstart, zte,
        ztend = NULL, plot = TRUE)

Arguments

Details

The cumulative number of earthquakes at time t since t_0 is given by the integration of \lambda(t) ( see etasap ) with respect to the time t,

\Lambda(t) = \mu(t-t_0) + K \Sigma_i \exp[\alpha(M_i-M_z)]\{c^{(1-p)}-(t-t_i+c)^{(1-p)}\}/(p-1),

where the summation of i is taken for all data event. The output of etarpp2 is given in a res-format dataset which includes the column of \{\Lambda(t_i), i = 1, 2, ..., N\}.

Value

References

Ogata, Y. (2006) Computer Science Monographs, No.33, Statistical Analysis of Seismicity - updated version (SASeies2006). The Institute of Statistical Mathematics.

Examples

data(main2003JUL26)  # The aftershock data of 26th July 2003 earthquake of M6.2

## output transformed times and cumulative numbers
x <- main2003JUL26
etarpp(time = x$time, mag = x$magnitude, threshold = 2.5, reference = 6.2,
       parami = c(0, 0.63348e+02, 0.38209e-01, 0.26423e+01, 0.10169e+01),
       tstart = 0.01, zte = 7, ztend = 18.68)

## output a res-format dataset
etarpp2(main2003JUL26, threshold = 2.5,  reference = 6.2,
        parami = c(0, 0.63348e+02, 0.38209e-01, 0.26423e+01, 0.10169e+01),
        tstart = 0.01, zte = 7, ztend = 18.68)

Maximum Likelihood Estimates of the ETAS Model

Description

Compute the maximum likelihood estimates of five parameters of ETAS model. This function consists of two (exact and approximated) versions of the calculation algorithm for the maximization of likelihood.

Usage

etasap(time, mag, threshold = 0.0, reference = 0.0, parami, zts = 0.0, 
       tstart, zte, approx = 2, tmpfile = NULL, nlmax = 1000, plot = TRUE)

Arguments

Details

The ETAS model is a point-process model representing the activity of earthquakes of magnitude M_z and larger occurring in a certain region during a certain interval of time. The total number of such earthquakes is denoted by N. The seismic activity includes primary activity of constant occurrence rate \mu in time (Poisson process). Each earthquake ( including aftershock of another earthquake) is followed by its aftershock activity, though only aftershocks of magnitude M_z and larger are included in the data. The aftershock activity is represented by the Omori-Utsu formula in the time domain. The rate of aftershock occurrence at time t following the ith earthquake (time: t_i, magnitude: M_i) is given by

n_i(t) = K exp[\alpha(M_i-M_z)]/(t-t_i+c)^p,

for t>t_i where K, \alpha, c, and p are constants, which are common to all aftershock sequences in the region. The rate of occurrence of the whole earthquake series at time t becomes

\lambda(t) = \mu + \Sigma_i n_i(t).

The summation is done for all i satisfying t_i < t. Five parameters \mu, K, c, \alpha and p represent characteristics of seismic activity of the region.

Value

References

Ogata, Y. (2006) Computer Science Monographs, No.33, Statistical Analysis of Seismicity - updated version (SASeies2006). The Institute of Statistical Mathematics.

Examples

data(main2003JUL26)  # The aftershock data of 26th July 2003 earthquake of M6.2 
x <- main2003JUL26
etasap(x$time, x$magnitude, threshold = 2.5, reference = 6.2, 
       parami = c(0, 0.63348e+02, 0.38209e-01, 0.26423e+01, 0.10169e+01),
       tstart = 0.01, zte = 18.68)

Simulation of Earthquake Dataset Based on the ETAS Model

Description

Produce simulated dataset for given sets of parameters in the point process model used in ETAS.

Usage

etasim1(bvalue, nd, threshold = 0.0, reference = 0.0, param)

etasim2(etas, tstart, threshold = 0.0, reference = 0.0, param)

Arguments

Details

There are two versions; either simulating magnitude by Gutenberg-Richter's Law etasim1 or using magnitudes from etas dataset etasim2. For etasim1, b-value of G-R law and number of events to be simulated are provided. stasim2 simulates the same number of events that are not less than threshold magnitude in the dataset etas, and simulation starts after a precursory period depending on the same history of events in etas in the period.

Value

etasim1 and etasim2 generate a etas-format dataset given values of 'no.', 'magnitude' and 'time'.

References

Ogata, Y. (2006) Computer Science Monographs, No.33, Statistical Analysis of Seismicity - updated version (SASeies2006). The Institute of Statistical Mathematics.

Examples

## by Gutenberg-Richter's Law 
etasim1(bvalue = 1.0, nd = 999, threshold = 3.5, reference = 3.5,
        param = c(0.2e-02, 0.4e-02, 0.3e-02, 0.24e+01, 0.13e+01))

## from a etas-format dataset
data(main2003JUL26)  # The aftershock data of 26th July 2003 earthquake of M6.2
etasim2(main2003JUL26, tstart = 0.01, threshold = 2.5, reference = 6.2,
        param = c(0, 0.63348e+02, 0.38209e-01, 0.26423e+01, 0.10169e+01))

Maximum Likelihood Estimates of Linear Intensity Models

Description

Perform the maximum likelihood estimates of linear intensity models of self-exciting point process with another point process input, cyclic and trend components.

Usage

linlin(external, self.excit, interval, c, d, ax = NULL, ay = NULL, ac = NULL,
       at = NULL, opt = 0, tmpfile = NULL, nlmax = 1000)

Arguments

Details

The cyclic part is given by the Fourier series, the trend is given by usual polynomial. The response functions of the self-exciting and the input are given by the Laguerre type polynomials (lgp), where the scaling parameters in the exponential function, say c and d, can be different. However, it is advised to estimate c first without the input component, and then to estimate d with the fixed c (this means that the gradient corresponding to the c is set to keep 0), which are good initial estimates for the c and d of the mixed self-exciting and input model.

Note that estimated intensity sometimes happen to be negative on some part of time interval outside the neighborhood of events. this take place more easily the larger the number of parameters. This causes some difficulty in getting the m.l.e., because the negativity of the intensity contributes to the seeming increase of the likelihood.

Note that for the initial estimates of ax(1), ay(1) and at(1), some positive value are necessary. Especially 0.0 is not suitable.

Value

References

Ogata, Y., Katsura, K. and Zhuang, J. (2006) Computer Science Monographs, No.32, TIMSAC84: STATISTICAL ANALYSIS OF SERIES OF EVENTS (TIMSAC84-SASE) VERSION 2. The Institute of Statistical Mathematics.

Ogata, Y. and Akaike, H. (1982) On linear intensity models for mixed doubly stochastic Poisson and self-exciting point processes. J. royal statist. soc. b, vol. 44, pp. 102-107.

Ogata, Y., Akaike, H. and Katsura, K. (1982) The application of linear intensity models to the investigation of causal relations between a point process and another stochastic process. Ann. inst. statist. math., vol. 34. pp. 373-387.

Examples

data(PProcess)  # point process data 
data(SelfExcit) # self-exciting point process data
linlin(PProcess[1:69], SelfExcit, interval = 20000, c = 0.13, d = 0.026,
       ax = c(0.035, -0.0048), ay = c(0.0, 0.00017), at = c(0.007, -.00000029))

Simulation of a Self-Exciting Point Process

Description

Perform simulation of a self-exciting point process whose intensity also includes a component triggered by another given point process data and a non-stationary Poisson trend.

Usage

linsim(data, interval, c, d, ax, ay, at, ptmax)

Arguments

Details

This function performs simulation of a self-exciting point process whose intensity also includes a component triggered by another given point process data and non-stationary Poisson trend. The trend is given by usual polynomial, and the response functions to the self-exciting and the external inputs are given the Laguerre-type polynomials (lgp), where the scaling parameters in the exponential functions, say c and d, can be different.

Value

References

Ogata, Y., Katsura, K. and Zhuang, J. (2006) Computer Science Monographs, No.32, TIMSAC84: STATISTICAL ANALYSIS OF SERIES OF EVENTS (TIMSAC84-SASE) VERSION 2. The Institute of Statistical Mathematics.

Ogata, Y. (1981) On Lewis' simulation method for point processes. IEEE information theory, vol. it-27, pp. 23-31.

Ogata, Y. and Akaike, H. (1982) On linear intensity models for mixed doubly stochastic Poisson and self-exciting point processes. J. royal statist. soc. b, vol. 44, pp. 102-107.

Ogata, Y., Akaike, H. and Katsura, K. (1982) The application of linear intensity models to the investigation of causal relations between a point process and another stochastic process. Ann. inst. statist math., vol. 34. pp. 373-387.

Examples

data(PProcess)   ## The point process data
linsim(PProcess, interval = 20000, c = 0.13, d = 0.026, ax = c(0.035, -0.0048), 
       ay = c(0.0, 0.00017), at = c(0.007, -0.00000029), ptmax = 0.007)

The Aftershock Data

Description

The aftershock data of 26th July 2003 earthquake of M6.2 at the northern Miyagi-Ken Japan.

Usage

data(main2003JUL26)

Format

main2003JUL26 is a data frame with 2305 observations and 9 variables named no., longitude, latitude, magnitude, time (from the main shock in days), depth, year, month, and day.

Source

Ogata, Y. (2006) Computer Science Monographs, No.33, Statistical Analysis of Seismicity - updated version (SASeies2006). The Institute of Statistical Mathematics.

Maximum Likelihood Estimates of Parameters in the Omori-Utsu (Modified Omori) Formula

Description

Compute the maximum likelihood estimates (MLEs) of parameters in the Omori-Utsu (modified Omori) formula representing for the decay of occurrence rate of aftershocks with time.

Usage

momori(data, mag = NULL, threshold = 0.0, tstart, tend, parami,
       tmpfile = NULL, nlmax = 1000)

Arguments

Details

The modified Omori formula represent the delay law of aftershock activity in time. In this equation, f(t) represents the rate of aftershock occurrence at time t, where t is the time measured from the origin time of the main shock. B, K, c and p are non-negative constants. B represents constant-rate background seismicity which may be included in the aftershock data.

f(t) = B + K/(t+c)^p

In this function the negative log-likelihood function is minimized by the Davidon-Fletcher-Powell algorithm. Starting from a given set of initial guess of the parameters parai, momori() repeats calculations of function values and its gradients at each step of parameter vector. At each cycle of iteration, the linearly searched step (lambda), negative log-likelihood value (-LL), and two estimates of square sum of gradients are shown (process=1).

The cumulative number of earthquakes at time t since t_0 is given by the integration of f(t) with respect to the time t,

F(t) = B(t-t_0) + K\{c^{1-p}-(t-t_i+c)^{1-p}\} / (p-1)

where the summation of i is taken for all data event.

Value

References

Ogata, Y. (2006) Computer Science Monographs, No.33, Statistical Analysis of Seismicity - updated version (SASeies2006). The Institute of Statistical Mathematics.

Examples

data(main2003JUL26)  # The aftershock data of 26th July 2003 earthquake of M6.2 
x <- main2003JUL26
momori(x$time, x$magnitude, threshold = 2.5, tstart = 0.01, tend = 18.68,
       parami = c(0,0.96021e+02, 0.58563e-01, 0.96611e+00))

Graphical Outputs for the Point Process Data Set

Description

Provide the several graphical outputs for the point process data set.

Usage

pgraph(data, mag, threshold = 0.0, h, npoint, days, delta = 0.0, dmax = 0.0,
       separate.graphics = FALSE)

Arguments

Value

References

Ogata, Y., Katsura, K. and Zhuang, J. (2006) Computer Science Monographs, No.32, TIMSAC84: STATISTICAL ANALYSIS OF SERIES OF EVENTS (TIMSAC84-SASE) VERSION 2. The Institute of Statistical Mathematics.

Ogata, Y. (2006) Computer Science Monographs, No.33, Statistical Analysis of Seismicity - updated version (SASeies2006). The Institute of Statistical Mathematics.

Ogata, Y. and Shimazaki, K. (1984) Transition from aftershock to normal activity: The 1965 Rat islands earthquake aftershock sequence. Bulletin of the seismological society of America, vol. 74, no. 5, pp. 1757-1765.

Examples

## The aftershock data of 26th July 2003 earthquake of M6.2
data(main2003JUL26)
x <- main2003JUL26
pgraph(x$time, x$magnitude, h = 6, npoint = 100, days = 10)

## The residual point process data of 26th July 2003 earthquake of M6.2
data(res2003JUL26)
y <- res2003JUL26
pgraph(y$trans.time, y$magnitude, h = 6, npoint = 100, days = 10)

The Periodogram of Point Process Data

Description

Provide the periodogram of point process data with the significant band (0.90, 0.95 and 0.99) of the maximum power in searching a cyclic component, for stationary Poisson Process.

Usage

ptspec(data, nfre, prdmin, prd, nsmooth = 1, pprd, interval, plot = TRUE)

Arguments

Value

References

Ogata, Y., Katsura, K. and Zhuang, J. (2006) Computer Science Monographs, No.32, TIMSAC84: STATISTICAL ANALYSIS OF SERIES OF EVENTS (TIMSAC84-SASE) VERSION 2. The Institute of Statistical Mathematics.

Examples

data(Brastings)  # The Occurrence Times Data of 627 Blastings
ptspec(Brastings, nfre = 1000, prdmin = 0.5, prd = 1.0, pprd = c(2.0, 1.0, 0.5),
       interval = 4600)

data(PoissonData)  # to see the contrasting difference
ptspec(PoissonData, nfre = 1000, prdmin = 0.5, prd =  1.0, pprd = c(2.0, 1.0, 0.5),
       interval = 5000)

The Residual Point Process Data

Description

The residual point process data of 26th July 2003 earthquake of M6.2 at the northern Miyagi-Ken Japan.

Usage

data(res2003JUL26)

Format

res2003JUL26 is a data frame with 553 observations and 7 variables named no., longitude, latitude, magnitude, time (from the main shock in days), depth, Ft (transformed time).

Source

Ogata, Y. (2006) Computer Science Monographs, No.33, Statistical Analysis of Seismicity - updated version (SASeies2006). The Institute of Statistical Mathematics.

The Residual Point Process of the ETAS Model

Description

Compute the residual of modified Omori Poisson process and display the cumulative curve and magnitude v.s. transformed time.

Usage

respoi(time, mag, param, zts, tstart, zte, threshold = 0.0, plot = TRUE)

respoi2(etas, param, zts, tstart, zte, threshold = 0.0, plot = TRUE)

Arguments

Details

The function respoi and respoi2 compute the following output for displaying the goodness-of-fit of Omori-Utsu model to the data. The cumulative number of earthquakes at time t since t_0 is given by the integration of f(t) with respect to the time t,

F(t) = B(t-t_0) + K\{ c^{(1-p)}-(t-t_i+c)^{(1-p)} \}/(p-1)

where the summation of i is taken for all data event.

respoi2 is equivalent to respoi except that input and output forms are different. When a etas-format dataset is given, respoi2 returns the dataset with the format as described below.

Value

References

Ogata, Y. (2006) Computer Science Monographs, No.33, Statistical Analysis of Seismicity - updated version (SASeies2006). The Institute of Statistical Mathematics.

Examples

data(main2003JUL26)  # The aftershock data of 26th July 2003 earthquake of M6.2 

# output transformed times and cumulative numbers
x <- main2003JUL26
respoi(x$time, x$magnitude, param = c(0,0.96021e+02, 0.58563e-01, 0.96611e+00),
       zts = 0.0, tstart = 0.01, zte = 18.68, threshold = 2.5)

# output a res-format dataset
respoi2(main2003JUL26, param = c(0,0.96021e+02, 0.58563e-01, 0.96611e+00),
        zts = 0.0, tstart = 0.01, zte = 18.68, threshold = 2.5)

Simulation of Bi-Variate Hawkes' Mutually Exciting Point Processes

Description

Perform the simulation of bi-variate Hawkes' mutually exciting point processes. The response functions are parameterized by the Laguerre-type polynomials.

Usage

simbvh(interval, axx = NULL, axy = NULL, axz = NULL, ayx = NULL, 
       ayy = NULL, ayz = NULL, c, d, c2, d2, ptxmax, ptymax)

Arguments

Value

References

Ogata, Y., Katsura, K. and Zhuang, J. (2006) Computer Science Monographs, No.32, TIMSAC84: STATISTICAL ANALYSIS OF SERIES OF EVENTS (TIMSAC84-SASE) VERSION 2. The Institute of Statistical Mathematics.

Ogata, Y. (1981) On Lewis' simulation method for point processes. IEEE Information Theory, IT-27, pp.23-31.

Examples

simbvh(interval = 20000,
       axx = 0.01623,
       axy = 0.007306,
       axz = c(0.006187, -0.00000023),
       ayz = c(0.0046786, -0.00000048, 0.2557e-10),
       c = 0.4032, d = 0.0219, c2 = 1.0, d2 = 1.0,
       ptxmax = 0.0062, ptymax = 0.08)