No. 32

COMPUTER SCIENCE MONOGRAPHS
A Publicationof
The Institute of Statistical Mathematics

TIMSAC84: STATISTICAL ANALYSIS OF SERIES OF EVENTS(TIMSAC84-SASE) VERSION 2
byY. Ogata, K. Katsura and J. Zhuang

April 2006

Computer Science Monographs

Editors
Junji NANAKO(in chief )Satoshi ITO
Nobuhisa KASHIWAGI
Yumi TAKIZAWA

The series of Computer Science Monographs publishes the statisticalsoftware and the articles on computer applications originating at The
Institute of Statistical Mathematics.

All communications relating to thispublication should be addressed to the Editorial Office, The Institute of
Statistical Mathematics, 4-6-7 Minami-Azabu, Minato-ku, Tokyo 106-8569,
Japan.

The electronic version of this issue is downloadable at
http://www.ism.ac.jp/

1

TIMSAC84: Statistical Analysis of Series ofEvents (TIMSAC84-SASE)
Version 2
Ogata, Y., Katsura, K. and Zhuang, J.
The Institute of Statistical Mathematics
TIMSAC84: Statistical Analysis of Series of Events (TIMSAC84-SASE) ............................................1
Overview of the program [eptren]...................................................................................................3
Overview of the program [pgraph] ..................................................................................................5
Overview of the program [ptspec] .................................................................................................10
Overview of the program [linlin]...................................................................................................11
Overview of the program [simbvh]................................................................................................13
Overview of the program [linsim] .................................................................................................14This program package, TIMSAC84-SASE version2, is an updated version of TIMSAC84-SASE by
Ogata and Katsura, which is included in Part 2 of TIMSAC84 (Akaike et al., 1985). The programsin the original TIMSAC package for point-process analysis were developed for the PCDOS (Utsuand Ogata, 1997). The main features of the updation is that the calculation outputs are simplified byremoving the GPSL language for plotting figures, so that users can use their own devices or
languages for plotting the figures. The source files of the present versions are written in
FORTRAN77. For the convenience of users the execution files here were compiled by GNU
FORTRAN “g77” on the cygwin, and should be executable at least on Windows XP. Then, we
provide some source modules of the free graphical software R and postscript files in this library. All
the names are corresponding to the names of the FORTRAN sources and the output files.
All programs are downloadable in zip format from the "Computer Science Monographs" page in
WWW pages of the Institute of Statistical Mathematics. They can be also retrieved fromhttp://bemlar.ism.ac.jp/www2/SASeisUpCollection/TIMSAC84-SASE-v2/

or

http://www.ism.ac.jp/~ogata/Statsei4/download.html
The latter is the same as the one distributed at the 4th Statistical Seismology Workshop held at
Hayama Campus of the Graduate University of Advanced Studies, Japan, during 9-13 January 2006.

1

2

Some technical notes are as follows:

1. If the programs crash with an error message "cannot find cygwin1.dll", copy cygwin1.dll toc:¥windows¥system32.

2. To install updated version of R for windows, MacOS or some other UNIX based systems, please
go to "http://www.r-project.org/".

3. To install most recent version of cygwin, please go to "http://www.cygwin.com/".

About the Data

In order to demonstrate the outputs from each program, we use the following two datasets as input
examples.

The example data [Brasting.data] consists of the occurrence times of 627 blastings at

a certain stoneyard with a very small portion of microearthquakes during past 4600days. We can see
the long-term trend and the daily cyclic feature from fitting non-stationary Poisson process, theelementary features of point processes and the periodogram (spectral estimation), by using the
programs [eptrend], [pgraph] and [ptspec], respectively. Also, [Poisson.data] is prepared to see
the contrasting difference. These datasets are copied as [work.data] in the directory of[WorkDomain].

References
Akaike, H. Ozaki, T. Ishiguro, M., Ogata, Y., Kitagawa, G., Tamura, Y., Arahata, E., Katsura, K. andTamura, R. (1985) Time Series and Control Program Package, TIMSAC-84, Computer Science
Monograph, No. 22/23, The Institute of Statistical Mathematics, Tokyo, Japan.
Utsu, T. and Ogata, Y. (1997) Computer program package: Statistical Analysis of point processes
for Seismicity, SASeis, IASPEI Software Library for personal computers, the International
Association of Seismology and Physics of Earth's Interior in collaboration with the American
Seismological Society, Vol. 6, pp. 13-94.

2

3Overview of the program [eptren]
This program computes the maximum likelihood estimates (MLEs) of the coefficients A1, A2, …, An
in an exponential polynomial
f(t) = exp{A1 + A2 t + A3 t2 + … }(1)
or A1, A2, B2, … , An, Bn in a point process model with a conditional intensity taking the form of an
exponential Fourier seriesf(t) = exp{A1+A2cos(2πt/P)+B2sin(2πt/P)+A3cos(4πt/P)+B3sin(4πt/P)+… }(2)
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 processes. The optimal order ncan be determined by minimize the value of the Akaike Information Criterion (AIC) (cf., Akaike,1974).
B

B

B

B

The maximum number of coefficients is limited to 21 due to the dimension set in the FORTRAN
source. The coefficients A1, A2, . . . , A21 can be determined for equation (1) and A1, A2, B2, … , A10,B10 can be determined for equation (2) at the most.
B

B

Since a non-stationary Poisson process is assumed in the data, we need a careful interpretation of
the results when the model is applied to data containing earthquake clusters such as aftershocks,
foreshocks, and earthquake swarms.
R-modules are also attached in order to display the outputs of numerical data by plotting graphs ofthe estimated functions.
Structure of the program is the following.
[eptren]|---[inputs]|---[reduc1]|---[reduc2]|---[davidn]------------[funct]||---[hesian]||---[linear]----[funct]||---[fincal]|---[output]---[printr]---[trenfn]|-----[cyclfn]
1. The name of the input data file should be always [work.data] with required format. Before
analyzing any other data, one must record them in required format and name the data file as[work.data]. Two test data files [Brastings.data] and [Poisson.data] for this program are3

4

attached in the directory of [DataExamples]. When dataset includes remarkable clusters, theinterpretation of the trend and cycles has to be carefully made.
2. The input control file [eptren.open]selects either trend or cycle fitting and specify necessary
variables.
3. The output files of [eptern] are [out.eptren1] and [out.eptren2]. The R-module [r. eptrend] or[r. epcycle] is used to plot figures in both console and postscripts [eptrend.ps] and[epcycle.ps].

4

5

4. Calculated records of each program are stored by the name [*.print] in the directory of[Calculations].
This program was originally designed (January 1985) and revised (December 2005) by Yosihiko
Ogata, and programmed and also reprogrammed by Koichi Katsura. The R-module was designed
and programmed by Jiancang Zhuang and Yosihiko Ogata (December 2005).
References
Akaike, H. (1974). A new look at the statistical model identification, IEEE Trans. Automat. Control,
AC-19, 716-723.
Lewis, P.A.W. and G.S. Shedler (1976). Statistical analysis of non-stationary series in a databasesystem, IBM. J. Res. Develop., 20, 465-481.
Maclean, C.J. (1974). Estimation and testing of an exponential polynomial rate function within the
non-stationary Poisson process, Biometrika, 61, 81-86.

Overview of the program [pgraph]
This program makes computations to produce several graphical outputs for the point process data
set.

5

6

1. Subroutine structures of the FORTRAN program [pgraph.f].
[pgraph]|-------[input]|-------[count1]---[shimiz]|-------[surviv]----[errbr2]----[plsinv]||--------[unifrm]---[unitsq]||-------[palmpr]|-------[vtc]|-------[vtcprt]
2. Cumulative numbers of events versus time, and positions of spikes (subroutine cumlat) are givenin the output file [out.pgCumMT]. The R-module [r.pgCumMT] can be used to plot figures inboth the console and postscript [pgCumMT.ps].

3. The output file [out.pgPTnum] gives the series of counted events in the moving interval of a
fixed length. The R-module [r.pgPTnum] produces figures based on the output file[out.pgPTnum] in both the console and postscript [pgPTnum.ps]. In the figure dotted lines
indicate i x (standard error), i = 1, 2, 3, assuming the stationary Poisson process ( subroutine
count1 or count2 ).

6

7

4. The R-module [r.pgSurviv] produces figures based on the output files [out.pgSurviv] and[out.pgSurDev] in both the console and postscript [pgSurviv.ps]. The figures include the log
survivor curve with i x (standard error), i = 1, 2, 3, based on the assumptions of the stationary
Poisson process, and a similar graph in which (x, y) plots are rotated and shifted in such a waythat the standard error lines and expectation lines are parallel (subroutine surviv).

78

5. The R-module [r.pgInterP] plots two figures based on the output files [out.pgInter1] and[out.pgInter2] in both the console and postscript [pgInterP.ps]. The first diagram is on theempirical distribution of u(i) = exp{-μx(i)} where x(i) is the i-th interval length between
consecutive points, and μ = (number of events) / (length of whole interval). The dashed linesof .95 and .99 significance bands of the two-sided Kolmogorov-Smirnov test assuming theuniform distribution is given. The related graph of ( u(i), u(i+1) ) plots are also carried out(subroutine uniform).

6. The R-module [r.pgPalm] produces, in both the console and postscript [pgPalm.ps], a figure
estimated intensity p(t) under the palm probability based on the output file [out.pgPalm]. The
intensity p(t) is related to the covariance density c(t) by the relation of c(t) = μ{p(t) – μ}, where
μ is the mean intensity of the point process. The 0.95 and 0.99 error bands are provided
assuming the stationary Poisson process (cf., subroutine palmxy, palmpr).

89

7. The R-module [r.pgVTC] produces, in both the console and postscript [pgVTC.ps], a figure of
the estimated variance-time curve with the 0.95 and 0.99 error lines assuming the stationary
Poisson process (subroutine vtcxyp, vtcprt), based on the output file [out.pgVTC].

8. The input data file for this program should always be named [work.data]. Before analyzing anyother data, recorded them in the same format and name the data file as [work.data]. Two
example input data files [Brastings.data] and [Poisson.data] are attached in the directory of9

10[DataExamples]. When the dataset includes remarkable clusters, one must be careful when
making interpretation to the trend and cycles.
9. The input control file [pgraph.init] gives and explains necessary variables.
10. Calculation records are stored in the file [*.print] in the directory of [Calculations].
This program was originally designed (January 1985) and revised (December 2005) by Yosihiko
Ogata, and programmed and reprogrammed by Koichi Katsura. The R-module was designed andprogrammed by Jiancang Zhuang and Yosihiko Ogata (December 2005).
References
Cox, D.R. and Lewis, P.A.W. (1966) The Statistical Analysis of Series of Events, Methuen & Co.
Ltd, London; also published from Chapman & Hall, London.

Overview of the program [ptspec]
This program provides the periodogram (the spectrum estimate) of point process data with thesignificant bands (0.90, 0.95 and 0.99) of the maximum power in searching a cyclic component, forstationary Poisson processes. Powers of interested frequencies or periods are plotted by + sign (inred color).
Structure of the program:[ptspec]|---[input]|---[cycle]|---[period]|---[printr]|---[smooth]This is originally designed by Yosihiko Ogata and Tohru Ozaki, and programmed by Koichi Katsura
and T. Ozaki (January 1985), and modified by Katsura and Ogata (December 2005).
R-modules are also attached in order display the outputs of numerical data for plotting graphs of theestimated functions.
1. The data file for this program should be always named [work.data]. Before analyzing any otherdata, one should record them in required format and name the data file as [work.data]. Two test
data files [Brastings.data] and [Poisson.data] are attached in the directory of [DataExamples].
1011

2. The input control file is [ptspec.init], which provides necessary variables to read by the program.
3. The output file of [ptspec] is [out.ptspec]. The R-module [r. ptspec] displays the outputs byplotting graphs of the periodogram in both the console and postscript [ptspec.ps].

4. Calculation records of [ptspec] are stored by the name [ptspec.print] in the directory of[Calculations].
The R-module was designed and programmed by Jiancang Zhuang and Yosihiko Ogata (December2005).
Reference
Vere-Jones and Ozaki, T. (1982). Some examples of statistical estimation applied to earthquake data,
1. Cyclic Poisson and self-exciting models, Ann. Inst. Statist. Math., 34, pp. 189-207.

Overview of the program [linlin]
This program performs the maximum likelihood estimates of linear intensity models of self-exciting
point process with another point process input, cyclic and trend components. The cyclic part isgiven by the Fourier series, the trend is given by usual polynomials. The response functions to theself-exciting and the external inputs are given by the Laguerre-type polynomials, where the scaling

11

12

parameters in the exponential function, say c and d, can be different. However it is advised to fit amodel without the input component firstly to estimate an initial c, and then to fit a model with boththe self-exciting and the external input component to estimate an initial d with c being fixed (thismeans that the gradient corresponding to c is fixed at 0), which are good initial estimates for the cand d of the mixed self-exciting and input model. It is recommended to use Akaike Information
Criterion (AIC; Akaike, 1974) for selecting optimal orders of respective components. Finally, it
should be noted that estimated intensity sometimes happens to be negative on some parts of thetime interval outside the neighbourhood of events. This take place more easily the larger the numberof parameters. This causes some difficulty in getting the MLE, because the negativity of theintensity contributes to the seeming increase of the likelihood. See the references below.
Structure of the program[linlin]|---[input]|---[comfac]|---[cycle]|---[dav]-------[davidn]----------[funct1 or 2]|---[hesian]|---[linear]---[funct1 or 2]
Remarks:
1. The input data file for this program is always named [linlin.data]. Before analyzing any otherdata, record them in the same format and name the file [linlin.data].
2. The input control file is [linlin.init], in which you can select a model with a possible combinationof the trend, the self-exciting or/and the external response components, and in which you canspecify necessary variables and initial parameter values.
3. Calculation results of [linlin] are stored in the file [linlin.print]. In the directory of[Calculations], there is an example file [linlin.print] for your initial check of the program.
This program was originally designed (January 1985) and revised (December 2005) by Yosihiko
Ogata, and programmed and also reprogrammed by Koichi Katsura.
References
Akaike, H. (1974). A new look at the statistical model identification, IEEE Trans. Automat. Control,
AC-19, 716-723.
Ogata. Y. and Akaike, H. (1982). On linear intensity models for mixed doubly stochastic Poisson

12

13

and self-exciting point processes, J. Roy. Statist. Soc. Ser. B, 44, pp. 102-107.
Ogata, Y., Akaike, H. and Katsura, K. (1982). The application of linear intensity models to theinvestigation of causal relations between a point process and another stochastic process, Ann.
Inst. Statist. Math., 34. pp. 373-387.
Ogata, Y. and Katsura, K. (1986). Point process model with linearly parametrized intensity for theapplication to earthquake data, Essays in Time Series and Allied Processes (Festscrift for Prof.
E.J. Hannan), J. Gani and M.B. Priestley eds., J. Appl. Probab. 23A, pp.291-310.

Overview of the program [simbvh]
This program performs the simulation of bi-variate Hawkes' mutually exciting point processes.
Non-stationary Poisson trends can also be included in the both intensity functions. The trend isgiven by a usual polynomial and the response functions to the self-exciting and the external inputs
are given by the Laguerre-type polynomials, where the scaling parameters in the exponentialfunctions, say c, d, c2 and d2, can be different. More detailed explanation of parameters andvariables are given in the control file [simbvh.init].
The simulated data can be displayed in the console and be saved in the file [simbvh.print]. In order
to view the graphical features of the data, save each of the two components of the simulated data as
the name of [work.data] with the same format as the original [work.data], and use program[pgraph].
Finally, it should be noted that the parameters for simulation must be specified such that theconditional intensity function be positive throughout the simulation time interval.
Structure of the program[simbvh]|---[input]|---[initl]|---[simda]----[comfac]||------[duf]||------[unifor]||------[fx]-------ptrend]||---[output]

13

14

Remarks:
1. No input data file for this program is necessary.
2. The input control file [simbvh.init] selects a simulation of either self-exciting point processes orself and mutually exciting processes and also specifies necessary variables.
3. Calculatation records of [simbvh] are stored in [simbvh.print] in the directory of[Calculations].
This program was originally designed (January 1985) and revised (December 2005) by Yosihiko
Ogata, and programmed and also reprogrammed by Koichi Katsura.
References
Hawkes, A.G. (1971). Spectra of some self-exciting and mutually exciting point processes,
Biometrika, 58, pp. 83-90.
Ogata, Y. (1981). On Lewis' simulation method for point processes, IEEE Information Theory,
IT-27, pp. 23-31.
Ogata. Y. and Akaike, H. (1982). On linear intensity models for mixed doubly stochastic Poissonand self-exciting point processes, J. Roy. Statist. Soc. Ser. B, 44, pp. 102-107.
Ogata, Y., Akaike, H. and Katsura, K. (1982). The application of linear intensity models to theinvestigation of causal relations between a point process and another stochastic process, Ann.
Inst. Statist. Math., 34. pp. 373-387.
Ogata, Y. and Katsura, K. (1986). Point process model with linearly parametrized intensity for theapplication to earthquake data, Essays in Time Series and Allied Processes (Festscrift for Prof.
E.J. Hannan), J. Gani and M.B. Priestley eds., J. Appl. Probab. 23A, pp. 291-310.

Overview of the program [linsim]
This program performs simulation of a self-exciting point process whose intensity also includes acomponent triggered by another given point process data and a 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 by the Laguerre-type polynomials, where the scaling parameters in the exponentialfunctions, say c and d, can be different. More detailed explanation of parameters and variables are
given in the control file [linsim.init].
The simulated data are displayed on the console and also save in the file [linsim.print] with the

14

15

same format as the original [work.data]. In order to view the graphical features of the data by using
the program [pgraph], change the name of the file for the simulated data to [work.data]. Similarly,
the data file of external input [linsim.data] can be also viewed by change its name to [work.data].
Finally, it should be noted that the parameters for simulation must be specified such that theconditional intensity function be positive throughout the simulation time interval.
Structure of the program[linsim]|---[input]|---[initl]|---[simxx]----[comfac]||------[duf]||------[unifor]||------[fx]--------[ptrend]||---[output]
Remarks:
1. The input data file for this program should be always named [linsim.data]. Before using anyother data as the external input point process for the simulation, one must save the data in therequired format and name the data file as [linsim.data].
2. The input control file chooses the combination of the trend, the self-exiting and the external input
components and specifies necessary variables.
3. Calculation outputs of [linsim] are stored by the name [linsim.print] in the directory of[Calculations].
This program was originally designed (January 1985) and revised (December 2005) by Yosihiko
Ogata, and programmed and also reprogrammed by Koichi Katsura.
References
Ogata, Y. (1981). On Lewis' simulation method for point processes, IEEE Information Theory,
IT-27, pp. 23-31.
Ogata. Y. and Akaike, H. (1982). On linear intensity models for mixed doubly stochastic Poissonand self-exciting point processes, J. Roy. Statist. Soc. Ser. B, 44, pp. 102-107.
Ogata, Y., Akaike, H. and Katsura, K. (1982). The application of linear intensity models to the

15

16

investigation of causal relations between a point process and another stochastic process, Ann.
Inst. Statist. Math., 34. pp. 373-387.
Ogata, Y. and Katsura, K. (1986). Point process model with linearly parametrized intensity for theapplication to earthquake data, Essays in Time Series and Allied Processes (Festscrift for Prof.
E.J. Hannan), J. Gani and M.B. Priestley eds., J. Appl. Probab. 23A, pp. 291-310.

16

The Institute of Statistical Mathematics
4-6-7 Minami-Azabu, Minato-ku
Tokyo 106-8569, Japan