Package 'SAPP' reference manual

Title:	Statistical Analysis of Point Processes
Description:	Functions for statistical analysis of point processes.
Authors:	The Institute of Statistical Mathematics
Maintainer:	Masami Saga <[email protected]>
License:	GPL (>= 2)
Version:	1.0.9-1
Built:	2024-12-19 06:43:51 UTC
Source:	CRAN

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

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)
eptren(data, mag = NULL, threshold = 0.0, nparam, nsub, cycle = 0,
       tmpfile = NULL, nlmax = 1000, plot = TRUE)

Arguments

`data`	point process data.
`mag`	magnitude.
`threshold`	threshold magnitude.
`nparam`	maximum number of parameters.
`nsub`	number of subdivisions in either ( $0$ , $t$ ) or ( $0$ , `cycle`), where $t$ is the length of observed time interval of points.
`cycle`	periodicity to be investigated days in a Poisson process model. If zero (default) fit an exponential polynomial model.
`tmpfile`	a character string naming the file to write the process of minimizing by Davidon-Fletcher-Powell procedure. If "" print the process to the standard output and if `NULL` (default) no report.
`nlmax`	the maximum number of steps in the process of minimizing.
`plot`	logical. If `TRUE` (default) intensity rates are plotted.

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

`aic`	AIC.
`param`	parameters.
`aicmin`	minimum AIC.
`maice.order`	number of parameters of minimum AIC.
`time`	time ( `cycle` = 0 ) or superposed occurrence time ( `cycle` > 0 ).
`intensity`	intensity rates.

References

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

`time`	the time measured from the main shock(t = 0).
`mag`	magnitude.
`etas`	a etas-format dataset on 9 variables (no., longitude, latitude, magnitude, time, depth, year, month and days).
`threshold`	threshold magnitude.
`reference`	reference magnitude.
`parami`	initial estimates of five parameters $\mu$ , $K$ , $c$ , $\alpha$ and $p$ .
`zts`	the start of the precursory period.
`tstart`	the start of the target period.
`zte`	the end of the target period.
`ztend`	the end of the prediction period. If `NULL` (default) the last time of available data is set.
`plot`	logical. If `TRUE` (default) the graphs of cumulative number and magnitude against the ordinary time and transformed time are plotted.

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

trans.time

transformed time $\Lambda(t_i), i = 1, 2, ..., N$ .

no.tstart

data number of the start of the target period.

resData

a res-format dataset on 7 variables

(no., longitude, latitude, magnitude, time, depth and transformed time).

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)
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)
etasap(time, mag, threshold = 0.0, reference = 0.0, parami, zts = 0.0, 
       tstart, zte, approx = 2, tmpfile = NULL, nlmax = 1000, plot = TRUE)

Arguments

`time`	the time measured from the main shock(t=0).
`mag`	magnitude.
`threshold`	threshold magnitude.
`reference`	reference magnitude.
`parami`	initial estimates of five parameters $\mu$ , $K$ , $c$ , $\alpha$ and $p$ .
`zts`	the start of the precursory period.
`tstart`	the start of the target period.
`zte`	the end of the target period.
`approx`	> 0 : the level for approximation version, which is one of the five levels 1, 2, 4, 8 and 16. The higher level means faster processing but lower accuracy. = 0 : the exact version.
`tmpfile`	a character string naming the file to write the process of maximum likelihood procedure. If "" print the process to the standard output and if `NULL` (default) no report.
`nlmax`	the maximum number of steps in the process of minimizing.
`plot`	logical. If `TRUE` (default) the graph of cumulative number and magnitude of earthquakes against the ordinary time is plotted.

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 $i$ th 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

`ngmle`	negative max log-likelihood.
`param`	list of maximum likelihood estimates of five parameters $\mu$ , $K$ , $c$ , $\alpha$ and $p$ .
`aic2`	AIC/2.

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)
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)
etasim1(bvalue, nd, threshold = 0.0, reference = 0.0, param)

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

Arguments

`bvalue`	$b$ -value of G-R law if `etasim1`.
`nd`	the number of the simulated events if `etasim1`.
`etas`	a etas-format dataset on 9 variables (no., longitude, latitude, magnitude, time, depth, year, month and days).
`tstart`	the end of precursory period if `etasim2`.
`threshold`	threshold magnitude.
`reference`	reference magnitude.
`param`	five parameters $\mu$ , $K$ , $c$ , $\alpha$ and $p$ .

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))
## 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)
linlin(external, self.excit, interval, c, d, ax = NULL, ay = NULL, ac = NULL,
       at = NULL, opt = 0, tmpfile = NULL, nlmax = 1000)

Arguments

`external`	another point process data.
`self.excit`	self-exciting data.
`interval`	length of observed time interval of event.
`c`	exponential coefficient of lgp in self-exciting part.
`d`	exponential coefficient of lgp in input part.
`ax`	coefficients of self-exciting response function.
`ay`	coefficients of input response function.
`ac`	coefficients of cycle.
`at`	coefficients of trend.
`opt`	0 : minimize the likelihood with fixed exponential coefficient $c$ 1 : not fixed $d$ .
`tmpfile`	a character string naming the file to write the process of minimizing. If "" print the process to the standard output and if `NULL` (default) no report.
`nlmax`	the maximum number of steps in the process of minimizing.

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

`c1`	initial estimate of exponential coefficient of lgp in self-exciting part.
`d1`	initial estimate of exponential coefficient of lgp in input part.
`ax1`	initial estimates of lgp coefficients in self-exciting part.
`ay1`	initial estimates of lgp coefficients in the input part.
`ac1`	initial estimates of coefficients of Fourier series.
`at1`	initial estimates of coefficients of the polynomial trend.
`c2`	final estimate of exponential coefficient of lgp in self-exciting part.
`d2`	final estimate of exponential coefficient of lgp in input part.
`ax2`	final estimates of lgp coefficients in self-exciting part.
`ay2`	final estimates of lgp coefficients in the input part.
`ac2`	final estimates of coefficients of Fourier series.
`at2`	final estimates of coefficients of the polynomial trend.
`aic2`	AIC/2.
`ngmle`	negative max likelihood.
`rayleigh.prob`	Rayleigh probability.
`distance`	= $\sqrt(rwx^2+rwy^2)$ .
`phase`	phase.

References

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))
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)
linsim(data, interval, c, d, ax, ay, at, ptmax)

Arguments

`data`	point process data.
`interval`	length of time interval in which events take place.
`c`	exponential coefficient of lgp corresponding to simulated data.
`d`	exponential coefficient of lgp corresponding to input data.
`ax`	lgp coefficients in self-exciting part.
`ay`	lgp coefficients in the input part.
`at`	coefficients of the polynomial trend.
`ptmax`	an upper bound of trend polynomial.

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

`in.data`	input data for `sim.data`.
`sim.data`	self-exciting simulated data.

References

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.

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)
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)
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)
momori(data, mag = NULL, threshold = 0.0, tstart, tend, parami,
       tmpfile = NULL, nlmax = 1000)

Arguments

`data`	point process data.
`mag`	magnitude.
`threshold`	threshold magnitude.
`tstart`	the start of the target period.
`tend`	the end of the target period.
`parami`	the initial estimates of the four parameters $B$ , $K$ , $c$ and $p$ .
`tmpfile`	a character string naming the file to write the process of minimizing. If "" print the process to the standard output and if `NULL` (default) no report.
`nlmax`	the maximum number of steps in the process of minimizing.

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

`param`	the final estimates of the four parameters $B$ , $K$ , $c$ and $p$ .
`ngmle`	negative max likelihood.
`aic`	AIC = -2 $LL$ + 2*(number of variables), and the number = 4 in this case.
`plist`	list of parameters $t_i$ , $K$ , $c$ , $p$ and $cls$ .

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))
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)
pgraph(data, mag, threshold = 0.0, h, npoint, days, delta = 0.0, dmax = 0.0,
       separate.graphics = FALSE)

Arguments

`data`	point process data.
`mag`	magnitude.
`threshold`	threshold magnitude.
`h`	time length of the moving interval in which points are counted to show the graph.
`npoint`	number of subintervals in (0, days) to estimate a nonparametric intensity under the palm probability measure.
`days`	length of interval to display the intensity estimate under the palm probability.
`delta`	length of a subinterval unit in (0, dmax) to compute the variance time curve.
`dmax`	time length of an interval to display the variance time curve; this is less than (length of whole interval)/4. As the default setting of either delta = 0.0 or dmax = 0.0, set dmax = (length of whole interval)/4 and delta = dmax/100.
`separate.graphics`	logical. If `TRUE` a graphic device is opened for each graphics display.

Value

`cnum`	cumulative numbers of events time.
`lintv`	interval length.
`tau`	= time * (total number of events)/(time end).
`nevent`	number of events in [tau, tau+h].
`survivor`	log survivor curve with i*(standard error), i = 1,2,3.
`deviation`	deviation of survivor function from the Poisson.
`nomal.cnum`	normalized cumulative number.
`nomal.lintv`	U(i) = -exp(-(normalized interval length)).
`success.intv`	successive pair of intervals.
`occur`	occurrence rate.
`time`	time assuming the stationary Poisson process.
`variance`	Var(N(0,time)).
`error`	the 0.95 and 0.99 error lines assuming the stationary Poisson process.

References

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

Poisson Data

Description

Poisson test data for ptspec.

Usage

data(PoissonData)
data(PoissonData)

Format

A numeric vector of length 2553.

Source

The Point Process Data

Description

The point process test data for linsim and linlin.

Usage

data(PProcess)
data(PProcess)

Format

A numeric vector of length 72.

Source

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)
ptspec(data, nfre, prdmin, prd, nsmooth = 1, pprd, interval, plot = TRUE)

Arguments

`data`	data of events.
`nfre`	number of sampling frequencies of spectra.
`prdmin`	the minimum periodicity of the sampling.
`prd`	a periodicity for calculating the Rayleigh probability.
`nsmooth`	number for smoothing of periodogram.
`pprd`	particular periodicities to be investigated among others.
`interval`	length of observed time interval of events.
`plot`	logical. If `TRUE` (default) the periodogram is plotted.

Value

`f`	frequency.
`db`	D.B.
`power`	power.
`rayleigh.prob`	the probability of Rayleigh.
`distance`	= $\sqrt(rwx^2+rwy^2)$ .
`phase`	phase.

References

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)
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)
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)
respoi(time, mag, param, zts, tstart, zte, threshold = 0.0, plot = TRUE)

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

Arguments

`time`	the time measured from the main shock (t = 0).
`mag`	magnitude.
`etas`	an etas-format dataset on 9 variables (no., longitude, latitude, magnitude, time, depth, year, month and days).
`param`	the four parameters $B$ , $K$ , $c$ and $p$ .
`zts`	the start of the precursory period.
`tstart`	the start of the target period.
`zte`	the end of the target period.
`threshold`	threshold magnitude.
`plot`	logical. If `TRUE` (default) cumulative curve and magnitude v.s. transformed time $F(t_i)$ are plotted.

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

trans.time

transformed time $F(t_i), i = 1,2,...,N.$

cnum

cumulative number of events.

resData

a res-format dataset on 7 variables

(no., longitude, latitude, magnitude, time, depth and trans.time)

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

Self-Exciting Point Process Data

Description

Self-exciting point process test data for linlin.

Usage

data(SelfExcit)
data(SelfExcit)

Format

A numeric vector of length 99.

Source

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)
simbvh(interval, axx = NULL, axy = NULL, axz = NULL, ayx = NULL, 
       ayy = NULL, ayz = NULL, c, d, c2, d2, ptxmax, ptymax)

Arguments

`interval`	length of time interval in which events take place.
`axx`	coefficients of Laguerre polynomial (lgp) of the transfer function (= response function) from the data events x to x (trf; x –> x).
`axy`	coefficients of lgp (trf; y –> x).
`ayx`	coefficients of lgp (trf; x –> y).
`ayy`	coefficients of lgp (trf; y –> y).
`axz`	coefficients of polynomial for x data.
`ayz`	coefficients of polynomial for y data.
`c`	exponential coefficient of lgp corresponding to xx.
`d`	exponential coefficient of lgp corresponding to xy.
`c2`	exponential coefficient of lgp corresponding to yx.
`d2`	exponential coefficient of lgp corresponding to yy.
`ptxmax`	an upper bound of trend polynomial corresponding to xz.
`ptymax`	an upper bound of trend polynomial corresponding to yz.

Value

`x`	simulated data X.
`y`	simulated data Y.

References

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

Package 'SAPP'

Help Index

Statistical Analysis of Point Processes

Description

Details

References

The Occurrence Times Data

Description

Usage

Format

Source

Maximum Likelihood Estimates of Intensity Rates

Description

Usage

Arguments

Details

Value

References

Examples

Residual Point Process of the ETAS Model

Description

Usage

Arguments

Details

Value

References

Examples

Maximum Likelihood Estimates of the ETAS Model

Description

Usage

Arguments

Details

Value

References

Examples

Simulation of Earthquake Dataset Based on the ETAS Model

Description

Usage

Arguments

Details

Value

References

Examples

Maximum Likelihood Estimates of Linear Intensity Models

Description

Usage

Arguments

Details

Value

References

Examples

Simulation of a Self-Exciting Point Process

Description

Usage

Arguments

Details

Value

References

Examples

The Aftershock Data

Description

Usage

Format

Source

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

Description

Usage

Arguments

Details

Value

References

Examples

Graphical Outputs for the Point Process Data Set

Description

Usage

Arguments

Value

References

Examples

Poisson Data