This article was downloaded by: [Adelfio, Giada] On: 19 April 2011 Access details: Access Details: [subscription number 936584859] Publisher Taylor & Francis Informa Ltd Registered in England and Wales Registered Number: 1072954 Registered office: Mortimer House, 3741 Mortimer Street, London W1T 3JH, UK
Communications in Statistics - Simulation and Computation
Publication details, including instructions for authors and subscription information: http://www.informaworld.com/smpp/title~content=t713597237
Gamma Kernel Intensity Estimation in Temporal Point Processes
G. Marcona; G. Adelfioa; M. Chiodia a Dipartimento di Scienze Statistiche e Matematiche “S. Vianelli”, Università degli Studi di Palermo, Palermo, Italy Online publication date: 18 April 2011
To cite this Article Marcon, G. , Adelfio, G. and Chiodi, M.(2011) 'Gamma Kernel Intensity Estimation in Temporal Point
Processes', Communications in Statistics - Simulation and Computation, 40: 8, 1146 — 1162 To link to this Article: DOI: 10.1080/03610918.2011.563158 URL: http://dx.doi.org/10.1080/03610918.2011.563158
PLEASE SCROLL DOWN FOR ARTICLE Full terms and conditions of use: http://www.informaworld.com/terms-and-conditions-of-access.pdf This article may be used for research, teaching and private study purposes. Any substantial or systematic reproduction, re-distribution, re-selling, loan or sub-licensing, systematic supply or distribution in any form to anyone is expressly forbidden. The publisher does not give any warranty express or implied or make any representation that the contents will be complete or accurate or up to date. The accuracy of any instructions, formulae and drug doses should be independently verified with primary sources. The publisher shall not be liable for any loss, actions, claims, proceedings, demand or costs or damages whatsoever or howsoever caused arising directly or indirectly in connection with or arising out of the use of this material.
Communications in Statistics—Simulation and Computation® , 40: 1146–1162, 2011 Copyright © Taylor & Francis Group, LLC ISSN: 0361-0918 print/1532-4141 online DOI: 10.1080/03610918.2011.563158
Gamma Kernel Intensity Estimation in Temporal Point Processes
Downloaded By: [Adelfio, Giada] At: 09:44 19 April 2011
G. MARCON, G. ADELFIO, AND M. CHIODI Dipartimento di Scienze Statistiche e Matematiche “S. Vianelli”, Università degli Studi di Palermo, Palermo, Italy In this article, we propose a nonparametric approach for estimating the intensity function of temporal point processes based on kernel estimators. In particular, we use asymmetric kernel estimators characterized by the gamma distribution, in order to describe features of observed point patterns adequately. Some characteristics of these estimators are analyzed and discussed both through simulated results and applications to real data from different seismic catalogs. Keywords Bounded domain; Gamma kernel estimator; Intensity function; Temporal point processes. Mathematics Subject Classification 60G55; 62G07.
1. Introduction Nonparametric estimation based on kernel approach is often used to get a functional form of unknown functions in absence of any guidance or constraints from theory. In seismological context Normal (i.e., Gaussian) kernel estimators have been used mostly to describe the spatial component of the seismicity (Choi and Hall, 1998; Jackson and Kagan, 1999). Zhuang et al. (2002) used kernel approaches for the spatial distribution of the background component of ETAS model. A spacetime nonparametric approach based on Gaussian kernel estimators was introduced by Adelfio et al. (2006) to define the intensity function of a clustering point-process for seismic data; some extensions were provided in Adelfio et al. (2010) and Adelfio (2010). In this article, we provide estimates of the intensity function of temporal point processes based on kernel estimators. In particular, we focus on situations in which the use of symmetric kernel functions could not be appropriate to describe the features of the observed point patterns. For instance, to describe temporal distribution of the triggered seismicity, the use of symmetric kernels (such as the Received July 27, 2010; Accepted February 4, 2011 Address correspondence to G. Adelfio, Dipartimento di Scienze Statistiche e Matematiche “S. Vianelli”, Università degli Studi di Palermo, viale delle Scienze, ed. 13, 90128, Palermo, Italy; E-mail: adelfi
[email protected]
1146
Downloaded By: [Adelfio, Giada] At: 09:44 19 April 2011
Gamma Kernel Intensity In Temporal Point Processes
1147
Normal one) might be not appropriate since strong earthquakes are usually followed by sequences of further events, close both in time and in space, reflecting a gradual release of the seismic stress in the source fault area and producing an increasing of the seismic rate. In general, when the natural domain of definition of a density to be estimated is not the whole real line but a bounded interval, a Gaussian kernel estimator might be biased, although consistency of kernel estimators is well documented when the support of the underlying function is unbounded (Rosenblatt, 1956; Silverman, 1986). In this article, we propose the use of the Gamma kernel estimators introduced by Chen (2000) to describe the intensity function of temporal point processes whose supports are bounded from one end only. We focus on the estimation of the intensity function of an inhomogeneous time point process by using Gamma kernel estimator and comparing results with Normal kernel estimator. In Sec. 2, the definition of the intensity function for temporal point processes is reviewed, introducing the Normal kernel intensity estimator ( ) and the Gamma kernel intensity estimators ( ) for intensity functions and bandwidth selection procedures for both estimators. In Sec. 3, simulation results concerning finite sample properties of the estimators for various distributions and parameter values are reported. Empirical studies of the shape of highly skewed seismic distributions based on real datasets are provided in Sec. 4. Some conclusions are finally provided in Sec. 5.
2. Temporal Intensity Function and Kernel Estimators Point processes can be specified mathematically in several ways, for instance by considering joint distributions of the counts of points in arbitrary sets or defining a complete intensity function, that is a function of the points history that generalizes the rate function of a Poisson process. More formally, a time point process is a random collection of points, each one representing the time of a single event. Any analytic temporal point process is uniquely characterized by its associated conditional intensity function t t (Daley and Vere-Jones, 2003) defined as the frequency with which events are expected to occur around a particular location in time, conditional on the prior history t of the point process up to time t: t t = lim
t→0
ENt t + t t t
(1)
where t is time increment of time t and ENt t + t t is the history-dependent expected value of occurrence in t t + t . When · does not depend on the past history but just on the current time t, the expression in (1) identifies an inhomogeneous Poisson process (as analyzed in this article), while a constant intensity characterizes homogeneous Poisson processes. Kernel estimators provide a simple and effective tool of finding structure in datasets without the imposition of a parametric model, with the advantage of being very intuitive and relatively simple to analyze mathematically. Let t1 t2 tN be the N ordered observations on the fixed interval 0 T of an inhomogeneous Poisson process with intensity function as in (1), for which the conditional intensity is independent on the past history, such that t t = t, and t can be regarded as a mean intensity.
1148
Marcon et al.
The natural kernel estimate of t is simply an extension of kernel density estimator, except that there is no division by N , since t is based on an expected number rather than an expected proportion. Therefore, the standard kernel estimator for intensity function is: ˆ = t
N
t ti h
(2)
Downloaded By: [Adelfio, Giada] At: 09:44 19 April 2011
i=1
where is a kernel function which determines the shape of the bumps such that tdt = 1 and h is a positive number usually called the smoothing parameter or bandwidth. The kernel function is usually symmetric and its choice is often considered a problem of less importance than the smoothing parameters. ˆ as an estimator · the mean integrated To measure the global accuracy of · squared error (MISE) and the Kullback–Leibler distance (KLD) are widely used. They are defined by: ˆ = E t ˆ − t 2 dt MISE (3) ˆ = KLD
t log
t dt ˆ t
(4)
As a particular case of (2), the Normal intensity kernel estimator is defined by: N t − t 1 i ˆ t = h i=1 h
(5)
where is a standard Normal density function. The use of Normal kernel functions in some application fields is actually not the best solution, since in case of a bounded support and asymmetric functions boundary bias could be observed. This problem is due to the use of a fixed kernel which assigns weight outside the support when smoothing is carried out near the boundary (Bouezmarni and Scaillet, 2005). To solve this issue, Gamma kernel estimators for intensity function based on asymmetric kernels with flexible form and location on the non-negative real line are introduced. 2.1. The Gamma Kernel Intensity Estimators The Gamma kernel estimators for density functions have been introduced by Chen (2000); they are based on an asymmetric kernel function defined on + with varying shape and location parameters, since depending on points location. These estimators are not subject to boundary bias and achieve the optimal rate of convergence for the MISE as for functions of the class of positive kernel estimators. The first Gamma kernel intensity estimator is defined as in (5) by using the gamma kernel defined by: tti /h exp−t/h (6) hti /h+1 t/h + 1 that is a Gamma density function with thi + 1 and h the shape and scale parameters respectively, such that the ith kernel function has its mode corresponding on the ∗ t ti h =
Gamma Kernel Intensity In Temporal Point Processes
1149
observation ti . This leads to the Gamma kernel intensity estimator ˆ ∗ t, defined by: ˆ ∗ t =
N
∗ t ti h
(7)
Downloaded By: [Adelfio, Giada] At: 09:44 19 April 2011
i=1
It is similar to the normal estimator defined in (5), replacing the symmetric kernel function with the gamma kernel ∗ . To reduce bias effects of the estimator ˆ ∗ t, due to the fact that each observation ti corresponds to the mode of the kernel function rather than its mean, Chen (2000) suggested the use of the kernel t ti h, that is a gamma kernel function characterized by a shape parameter h ti , such that: if ti ≥ 2h t /h h ti = 1i 2 t /h + 1 if 0 ≤ ti < 2h 4 i This provides the following Gamma kernel intensity estimator: ˆ t =
N
t ti h
(8)
i=1
with t ti h =
h ti exp−t/h hh ti h ti
Extending the results for densities by Chen (2000), the estimator (8) has usually better global performance (in terms of MISE) and uses a smaller bandwidth than (7) when h is fixed. Therefore, throughout this article, the Gamma kernel estimator refers to the one defined in (8) unless otherwise specified. Denoting with f = /N the univariate distribution probability density function which we are trying to estimate, Chen (2000) derived the asymptotic MISE for optimal choice of h for the corresponding Gamma kernel density estimator: 5 MISEopt fˆ ∼ 4/5 4
4/5
2 1/5 −4/5 1 −1/2 xf x dx x fxdx N √ 2 0 0
(9)
In this article, we use an approximated version of the asymptotic formulation of MISE for densities, since we deal with intensities of a inhomogeneous Poisson Process and N is a Poisson random variable with expected value T (see Sec. 3). Therefore, from the fundamental relationships and using the standard firstorder Taylor approximation for the variance of a product: VarX1 X2 ≈ E 2 X1 VarX2 + E 2 X2 VarX1 then: 2 ˆ ˆ ˆ MSEx = Varx + Biasx
≈ E 2 NVarfˆ x + E 2 fˆ xVarN + T2 Bias2 fˆ x
1150
Marcon et al. = T2 Varfˆ x + E 2 fˆ xT + T2 Bias2 fˆ x = T2 MSEfˆ x + E 2 fˆ xT
Finally, replacing Efˆ x by fx in the second term and integrating on x, we obtain the following approximation: ˆ ≈ 2 MISEfˆ + T MISE T
f 2 xdx
(10)
0
Downloaded By: [Adelfio, Giada] At: 09:44 19 April 2011
2.2. Bandwidth Selection Kernel density estimator requires the specification of a bandwidth h. This choice is very important, as it is shown in Figs. 1 and 2, where the individual bumps are shown as well as the estimates of the intensity functions constructed by adding them up, both for the Normal estimator ( ) and the Gamma one ( ). It should be noted the effect of varying the window width: (a) if h is too small then spurious fine structure becomes visible and ˆ has a higher variance; (b) if h is too large then the real nature of the distribution is obscured and ˆ has a lower variance. Although in this article we focus on the choice of the proper kernel function instead of the selection of smoothing parameters, we briefly summarize the properties of the cross-validation based approach (Hall, 1982), that is the one used throughout the article. A generalization of cross-validation is provided in Adelfio and Chiodi (2011), by the definition of the Forward Likelihood-based Predictive (FLP) approach aimed to prediction and accounting for the temporal ordering of observations. The bandwidth selector based on cross-validation approach is among the earliest fully automatic and consistent ones: cross-validation involves partitioning a
Figure 1. Normal Kernel estimates showing individual kernels. Window widths: (a) h = 0 1, (b) h = 0 3.
Downloaded By: [Adelfio, Giada] At: 09:44 19 April 2011
Gamma Kernel Intensity In Temporal Point Processes
1151
Figure 2. Gamma Kernel estimates showing individual kernels. Window widths: (a) h = 0 1, (b) h = 0 3.
sample of data into complementary subsets, performing the analysis on one subset (called the training set), and validating the analysis on the other subset (called the validation set or testing set). Let ˆ −i ti be the intensity function computed for observation ti and estimated by using the remaining data; then logLh =
N i=1
log ˆ −i ti −
T
ˆ tdt
(11)
0
The likelihood cross-validated selection of h is the value of h which maximize the function (11). A considerable amount has been written about selection of bandwidth for density estimation, but often the proposed methods depend on the unknown density (e.g., Silverman, 1978), implying restrictive limits to the accuracy of any possible data-driven bandwidth selector. The cross-validation based methods are the most interesting automatic methods of choosing a bandwidth (Hall, 1983; Stone, 1984); in particular, while crossvalidation based on quadratic Kullback–Leibler loss was considered by Bowman et al. (1984) and (Rudemo, 1982), cross-validation for density estimation has been shown to be asymptotically optimal (Hall, 1983) rather then simply consistent (Chow et al., 1983). However, the method is often computationally expensive because of the large number of the training process iterations, and sometimes the theoretical performance of this bandwidth selector is somewhat disappointing (Hall and Marron, 1987); nevertheless, when the loss function is chosen correctly, crossvalidation can be relied on to perform well for large samples (Bowman et al., 1984).
3. Simulation Results A simulation study was designed and carried out to investigate the performance of the Gamma kernel estimators (both (7) and (8)) and the Normal one (5) to
1152
Marcon et al.
Downloaded By: [Adelfio, Giada] At: 09:44 19 April 2011
estimate the intensity functions defined in different scenarios. We simulated an inhomogeneous Poisson process with a smoothed true intensity function, named true t and sample size n,obtained generating a Poisson random variable N with T assigned parameter T = 0 true tdt. Therefore, the intensity function is defined by true t = ftT , where ft is one of the six distributions described below, with either bounded or unbounded support and either one or two modes: (a) Gamma mixture obtained by mixing the distributions 1 1 and 3 2 with mixing proportion 21 ; (b) Gamma density with shape = 1 and scale = 2; (c) Gamma density with shape = 3 and scale = 1; (d) Normal density with mean = 5 and standard deviation = 1; (e) Pareto density with shape = 5 and scale = 1; (f) Normal mixture obtained by mixing the distributions 5 1 and 8 2 with mixing proportion 21 . For each of the above cases for the theoretical intensity function true t we used three values of T (100, 200, 500) and for each combination of true t and T we generated 100 simulated samples. For each simulation, the bandwidth is chosen as the one that maximizes logL(h) in (11). We do not report results relative to the first Gamma kernel estimator because simulations confirm the better performance of the estimator (8) than (7), as mentioned in Sec. 2.1; therefore, we focus on the comparison between the second Gamma kernel estimator and the Normal one. Intensities have been simulated for different values of T 100 200 500, showing an improving performance for increasing sample size. To make figures less crowded, we provide the simulated intensities together with their average just for the case T = 500 (see the plots in Fig. 3). These plots give an idea of the finite sample bias of the two estimators: it seems that bias is always smaller for the Gamma kernel estimates (solid lines) in comparison with the Normal kernel estimates (dotted lines) when the intensity functions are highly skewed (in Fig. 3, see plots (a), (b), and (e)). At the same time the Normal kernel function performs well when the real intensities are symmetric and defined on boundless support (in Fig. 3, see plots (c), (d), and (f)). Let us see Figs. 4–5 in order to analyze bias and variance of each estimator in detail: in Fig. 4 the averages of the intensity functions estimated over the 100 simulated samples both using the Gamma and the Normal kernel estimators are showed; Fig. 5 reports the corresponding simulated Mean Squared Errors (MSE’s). Both groups of figures confirm previous results for almost all the sample sizes considered across the entire range. As shown before, behaviors of the Normal kernel estimator and the true intensity are strongly different in cases (a), (b), and (e); on the other hand, the Gamma kernel estimator exhibits the same behavior both at the boundary and elsewhere as the true intensity function. The main difference between the Gamma and the Normal kernel estimators is showed in Fig. 5 plots (a), (b), and (e), where the simulated MSE’s of the Gamma kernel are much lower near t = 0; whereas plots (c), (d), and (f) in Fig. 5 show similar patterns over the entire support. The mean of the cross-validation log-likelihood for the Gamma and the Normal estimators computed over 100 simulated samples with different sizes and corresponding t-tests for paired samples are reported in Table 1; the Gamma kernel
1153
Downloaded By: [Adelfio, Giada] At: 09:44 19 April 2011
Gamma Kernel Intensity In Temporal Point Processes
Figure 3. Simulated intensity functions over 100 samples with T = 500 by using the Gamma (black solid line) and the Normal kernel estimator (dotted line). The true intensity (grey solid line) is: (a) Gamma mixture; (b) Gamma(1, 2); (c) Gamma(3, 1); (d) Normal(5, 1); (e) Pareto(5, 1); and (f) Normal mixture.
Marcon et al.
Downloaded By: [Adelfio, Giada] At: 09:44 19 April 2011
1154
Figure 4. Simulated mean of the 100 intensity functions with T = 500 obtained by using the Gamma kernel estimator (black solid line) and the Normal kernel estimator (dotted line). The true intensity (grey solid line) is: (a) Gamma mixture; (b) Gamma(1, 2); (c) Gamma(3, 1); (d) Normal(5, 1); (e) Pareto(5, 1); and (f) Normal mixture.
1155
Downloaded By: [Adelfio, Giada] At: 09:44 19 April 2011
Gamma Kernel Intensity In Temporal Point Processes
Figure 5. Simulated MSE to compare Gamma (solid line) and Normal kernel estimator (dotted line) over 100 samples with T = 500. The true intensity is: (a) Gamma mixture; (b) Gamma(1, 2); (c) Gamma(3, 1); (d) Normal(5, 1); (e) Pareto(5, 1); and (f) Normal mixture.
1156
Marcon et al.
Downloaded By: [Adelfio, Giada] At: 09:44 19 April 2011
Table 1 Mean cross-validation log-likelihood for the Gamma kernel estimator (Mean logL ) and the Normal one (Mean logL ) and corresponding t-tests computed over 100 simulated samples with T = 100 200 500 for the 6 theoretical intensity functions (a)–(f) described in Sec. 3 Mean logL
Mean logL
t-Test
p-Value
T = 100
(a) (b) (c) (d) (e) (f)
120 86 193 36 180 31 225 51 406 50 123 30
109 24 186 24 180 37 226 93 386 95 124 03
16 12 15 07 −0 34 −13 08 9 70 −14 54
0.00 0.00 0.74 0.00 0.00 0.00
T = 200
(a) (b) (c) (d) (e) (f)
383 58 517 69 492 30 588 37 935 06 384 09
363 87 503 47 492 32 589 91 888 55 384 68
16 61 13 84 −0 13 −21 11 9 73 −6 98
0.00 0.00 0.90 0.00 0.00 0.00
T = 500
(a) (b) (c) (d) (e) (f)
1397 43 1769 70 1676 61 1921 69 2825 39 1399 43
1362 07 1738 55 1676 15 1923 95 2726 82 1400 11
11 71 11 88 1 72 −14 81 8 39 −7 39
0.00 0.00 0.09 0.00 0.00 0.00
estimator fits significantly better than the Normal one for cases (a), (b), and (e), while for the cases (c), (d), and (f) there are not significant differences between the two estimators. Table 2 summarizes the MISE’s defined in (10) and the Kullback–Leibler distances, defined in (4), computed for both estimators over the different simulated intensity functions: we can observe that the performance of the Gamma and the Normal smoothers are close when the true intensity is symmetric whereas their difference increases for highly skewed functions. In other words, the Gamma kernel adequately estimates asymmetric intensity functions and its performance improves for greater sample sizes. The similarity between the Gamma empirical and theoretical MISE (obtained replacing (9) in (10)), is a further result validating our estimation procedures, though the approximation of the results. Asymptotic MISE evaluation for the Normal kernel estimator (Silverman, 1986) is not reported since dealing with densities with bounded domains the assumptions required for its applicability (Rosenblatt, 1956) concerning the continuity of derivatives and symmetry around zero are not satisfied by the simulated functions. Finally, a simulation study was carried out to investigate the performance of the Gamma kernel estimator (8) and the Normal one (5) in a very extreme situation where the true intensity function of an inhomogeneous Poisson process is
Gamma Kernel Intensity In Temporal Point Processes
1157
Table 2 Empirical Mean Integrated Squared Error (MISE) and Kullback–Leibler distance (KLD) for the Gamma and the Normal kernel estimators computed over 100 simulated samples with T = 100 200 500 for the 6 theoretical intensity functions (a)–(f) described in Sec. 3. The corresponding asymptotic theoretical MISE for the Gamma kernel estimator is also reported MISE
Downloaded By: [Adelfio, Giada] At: 09:44 19 April 2011
T = 100
T = 200
T = 500
(a) (b) (c) (d) (e) (f) (a) (b) (c) (d) (e) (f) (a) (b) (c) (d) (e) (f)
MISE
KLD
Empirical
Theoretical
Empirical
61 58 61 20 97 60 121 22 623 43 72 43 145 82 183 55 194 78 205 01 1451 62 151 78 499 80 528 34 509 64 573 54 3819 51 460 25
119 29 115 17 93 65 141 04 1117 04 90 77 267 85 257 15 209 58 315 63 2498 69 204 30 783 31 747 04 610 47 919 41 7274 54 599 19
310 56 303 61 105 19 115 46 4977 54 79 88 1005 06 1135 82 224 19 193 26 24070 73 167 17 4638 90 6242 72 751 46 524 04 125742 54 515 81
79 09 68 81 672 86 977 05 452 31 325 76 81 39 91 63 455 64 267 97 492 55 295 32 97 58 121 70 326 72 305 09 488 64 158 85
181 86 242 81 836 47 425 21 2862 19 4062 83 247 81 348 05 1023 21 276 51 6537 75 1553 64 379 25 620 39 1204 03 278 68 11220 38 202 29
Figure 6. Simulated mean (on the left) and MSE (on the right) over 100 samples with T = 500 for the Gamma kernel estimator (solid line) and the Normal kernel estimator (dotted line), assuming a step intensity function as in Eq. (12).
1158
Marcon et al.
proportional to a step function such that: 1 0 < t < 10 true t = 5 10 ≤ t ≤ 20
(12)
Simulation results over 100 samples for different values of T 100 200 500 show that both the Normal and Gamma kernel estimators provide similarly results, reflecting the hardness of this extreme case, mainly related to the presence of a bounded support and discontinuity points, in correspondence of which very high values of the simulated MSE are obtained (see Fig. 6 for T = 500).
Downloaded By: [Adelfio, Giada] At: 09:44 19 April 2011
4. Application to Real Data The empirical illustration concerns the analysis of four different seismic datasets. The estimation of the time intensity function is carried out for four seismic catalogs of different geographical areas. Each seismic catalog is formed by the space-time coordinates of seismic events. In this article, we focus just on temporal values. (a) Colfiorito (Fig. 7). The Colfiorito earthquake sequence (n = 2201) occurred on September 26, 1997, at 2:33 am in the area between Umbria and Marche, near Colfiorito (source: INGV – Istituto Nazionale di Geofisica e Vulcanologia). The seismic sequence was characterized both by the two mainshocks (2:33 am and 11:40 am, latitude 43 0 N and longitude 12 9 E) and by several aftershocks until November 3, 1997. (b) Palermo (Fig. 8). On September 6, 2002, at 3:21 am, a M = 5 6 earthquake, occurring some tens of kilometers offshore from the Northern Sicilian coast (latitude 13 7 N and longitude 38 4 E), slightly damaged the city of Palermo and surroundings. The seismic sequence (the mainshock and the several aftershocks, occurred for more than 80 days) is constituted by 810 events occurred in the
Figure 7. sequence.
Gamma (solid line) and Normal (dotted line) kernel estimates for (a) Colfiorito
Downloaded By: [Adelfio, Giada] At: 09:44 19 April 2011
Gamma Kernel Intensity In Temporal Point Processes
Figure 8. sequence.
1159
Gamma (solid line) and Normal (dotted line) kernel estimates for (b) Palermo
area 13.18–14.45 E and 37.95–38.95 N; the cluster has been drawn from the Southern Tyrrhenian Sea catalog (source INGV – Istituto Nazionale di Geofisica e Vulcanologia). (c) Landers (Fig. 9). The Landers sequence (n = 751) is drawn from the Southern California catalogs (source: U.S.G.S. – U.S. Geological Survey), that represents earthquakes occurred between January 1984 and June 2004 in a rectangular area around Los Angeles, between longitude −112 and −114 E and latitude 32–37 N; the dataset consists of 6796 earthquakes, with magnitude not smaller than 3.0. Landers is the earthquake identified by latitude 34 2 N, longitude
Figure 9. sequence.
Gamma (solid line) and Normal (dotted line) kernel estimates for (c) Landers
Downloaded By: [Adelfio, Giada] At: 09:44 19 April 2011
1160
Figure 10. sequence.
Marcon et al.
Gamma (solid line) and Normal (dotted line) kernel estimates for (d) Japan
−116 4 E, magnitude M = 7 3, and time June 28 1992, 4:57:31 am (Adelfio et al., 2006; Schoenberg, 2003). (d) Japan (Fig. 10). We use the hypocenter data recorded by the Japan Meteorological Agency, and consider one of several datasets discussed in Ogata (1998) and Ogata and Zhuang (2006). These data cover the area of the east coast of Tohoku District, one of the most active plate boundary zones in and around Japan. The data (n = 4333) consist of earthquakes of magnitude 4.5 and larger that are chosen from the region 36–42 N and 141–145 E for depths down to 100 km and for the time span 1926–1995. During 1938, the Shioya-Oki great swarm, including several large earthquakes, occurred near the coast of Fukusima Prefecture. We examine the intensity rate changes in time at the epicenter of the largest earthquake of M = 7 5 that occurred on November 5, 1938, with longitude 142 18 E and latitude 37 33 N (Adelfio and Ogata, 2010).
Table 3 Values of the h and of the log-likelihood estimated for (a) Colfiorito, (b) Palermo, (c) Landers, and (d) Japan data by cross-validation Gamma Kernel
(a) Colfiorito (b) Palermo (c) Landers (d) Japan
Normal Kernel
h
logL
h
logL
0 024 0 019 0 193 0 174
7202 02 1759 55 974 15 −136 43
0 113 0 678 10 059 3 044
7332 12 1718 82 619 44 −113 25
Gamma Kernel Intensity In Temporal Point Processes
1161
Downloaded By: [Adelfio, Giada] At: 09:44 19 April 2011
The Gamma and the Normal kernel estimators are applied to each seismic catalog. The cross-validation approach is used to select the bandwidth parameters. Figures 8 and 9 show that for the Palermo and Landers sequences, the Gamma kernel estimator catches the highest intensity near the origin, whereas the Normal kernel estimator does not. On the other hand, Figs. 7 and 10 show similar fitting of the Gamma and the Normal kernel estimators for Colfiorito and Japanese data. From these figures, we can observe that, for highly skewed distributions, the Gamma kernel estimator performs better than the Normal one mostly near the origin, suggesting that its choice could be crucial in those cases where bias effects are relevant. These results are summarized in Table 3, reporting the logL and h values estimated for each sequence.
5. Summary and Conclusion We compared the Gamma and the Normal kernel intensity estimators both when the intensity of the generating process is asymmetric, and defined on a limited support 0 +, and when it is a symmetric function with unbounded domain. For simulated cases, results show that the Gamma kernel estimator is less biased than the Normal one near the origin for limited support 0 +, while the choice between the two competitive estimators is irrelevant when the true intensity is symmetric. This point was also illustrated through a nonparametric estimation of the seismic activity of four real different catalogs, observing again that the Gamma kernel estimator performs better when the time behavior of the sequences is described by a skewed distribution, while both the Gamma and the Normal kernel provide good results in other cases. Therefore, although it is usual to consider the choice of the kernel function less important than the bandwidth selection, there are some situations in which the choice of an appropriate kernel function is a crucial issue. As a direction for future work, we are considering methods for simultaneous choice of bandwidth parameters and kernel functions in space-time point patterns.
Acknowledgment The authors thank the reviewers for their useful comments and suggestions which improved the quality of the article.
References Adelfio, G. (2010). Kernel estimation and display of a five-dimensional conditional intensity function. Nonlinear Processes in Geophysics 17:237–244. Adelfio, G., Chiodi, M. (2009). Second-order diagnostics for space-time point processes with application to seismic events. Environmetrics 20:895–911. Adelfio, G., Chiodi, M. (2011). Kernel intensity for space-time point processes with application to seismological problems. In: Classification and Multivariate Analysis for Complex Data Structures. Series: Studies in Classification, Data Analysis, and Knowledge Organization. Heidelberg: Springer-Verlag, pp. 401–408. Adelfio, G., Ogata, Y. (2010). Hybrid kernel estimates of space-time earthquake occurrance rates using the etas model. Annals of the Institute of Statistical Mathematics 62:127–143.
Downloaded By: [Adelfio, Giada] At: 09:44 19 April 2011
1162
Marcon et al.
Adelfio, G., Chiodi, M., De Luca, L., Luzio, D. (2006). Nonparametric clustering of seismic events. In: Data Analysis, Classification and the Forward Search. Series: Studies in Classification, Data Analysis, and Knowledge Organization. Heidelberg: Springer-Verlag, pp. 397–404. Adelfio, G., Chiodi, M., De Luca, L., Luzio, D. (2010). An algorithm for earthquake clustering based on maximum likelihood. In: Data Analysis and Classification. Series: Studies in Classification, Data Analysis, and Knowledge Organization. New York: Springer, pp. 25–32. Bouezmarni, T., Scaillet, O. (2005). Consistency of asymmetric kernel density estimators and smoothed histograms with application to income data. Technical Report 0306, Interuniversity Attraction Pole. Economic Theory 21:390–412. Bowman, A. W., Hall, P., Titterington, D. M. (1984). Cross-validation in nonparametric estimation of probabilities and probability density. Biometrika Trust 71:341–351. Chen, S. X. (2000). Probability density function estimating using gamma kernels. Annals of the Institute of Statistical Mathematics 52:471–480. Choi, E., Hall, P. (1998). On bias reduction in local linear smoothing. Biometrika 85(2): 333–345. Chow, Y. S., German, S., Wu, L. D. (1983). Consistent cross-validated density estimation. The Annals of Statistics 11:25–38. Daley, D. J., Vere-Jones, D. (2003). An Introduction to the Theory of Point Processes. 2nd ed. New York: Springer-Verlag. Hall, P. (1982). Cross-validation in density estimation. Biometrika 69(2):383–390. Hall, P. (1983). Large samples optimality of least squares cross-validation in density estimation. The Annals of Statistics 11:1156–74. Hall, P., Marron, J. S. (1987). Extent to which least-squares cross-validation minimizes integrated squared error in nonparametric density estimation. Probability Theory Related Fields 74:567–581. Jackson, D. D., Kagan, Y. Y. (1999). Testable earthquake forecasts for 1999. Seismlogical Research Letters 70:393–403. Ogata, Y. (1998). Space-time point-process models for earthquake occurrences. Annals of the Institute of Statistical Mathematics 50(2):379–402. Ogata, Y., Zhuang, J. (2006). Space-time etas model and an improved extension. Tectonophysics 413:13–23. Rosenblatt, M. (21956). Remarks on some nonparametric estimates of a density function. Annals of Mathematical Statistics 27:832–837. Rudemo, M. (1982). Empirical choice of histograms and kernel density estimators. Scandinavian Journal of Statistics 9:65–78. Schoenberg, F. P. (2003). Multi-dimensional residual analysis of point process models for earthquake occurrences. Journal American Statistical Association 98(464):789–795. Silverman, B. W. (1978). Choosing the window width when estimating a density. Biometrika Trust 65:1–11. Silverman, B. W. (1986). Density Estimation for Statistics and Data Analysis. London: Chapman and Hall. Stone, C. J. (1984). An asymptotically optimal window selection rule for kernel density estimates. The Annals of Statistics 12:1285–1297. Zhuang, J., Ogata, Y., Vere-Jones, D. (2002). Stochastic declustering of space-time earthquake occurrences. Journal of the American Statistical Association 97(458):369–379.
