Self-exciting point process modeling of crime
G. O. Mohler∗

M. B. Short†

P. J. Brantingham‡
G. E. Tita¶

F. P. Schoenberg§

Abstract
Highly clustered event sequences are observed in certain types of crime data, such as
burglary and gang violence, due to crime specific patterns of criminal behavior. Similarclustering patterns are observed by seismologists, as earthquakes are well known toincrease the risk of subsequent earthquakes, or aftershocks, near the location of an
initial event. Space-time clustering is modeled in seismology by self-exciting point
processes and the focus of this paper is to show that these methods are well suited for
criminological applications. We first review self-exciting point processes in the context
of seismology. Next, using residential burglary data provided by the Los Angeles Police
Department, we illustrate the implementation of self-exciting point process models in
the context of urban crime. For this purpose we use a fully non-parametric estimation
methodology to gain insight into the form of the space-time triggering function andtemporal trends in the background rate of burglary.

1

Introduction

Criminological research has shown that crime can spread through local environments via acontagion-like process (Johnson 2008). For example, burglars will repeatedly attack clusters∗

Department of Mathematics and Computer Science, Santa Clara University
Department of Mathematics, University of California, Los Angeles‡
Department of Anthropology, University of California, Los Angeles§
Department of Statistics, University of California, Los Angeles¶
Department of Criminology, University of California, Irvine†

1

of nearby targets because local vulnerabilities are well known to the offenders (Bernasco andNieuwbeerta, 2005). A gang shooting may incite waves of retaliatory violence in the local set
space (territory) of the rival gang (Tita and Ridgeway, 2007; Cohen and Tita, 1999). The
local, contagious spread of crime leads to the formation of crime clusters in space and time.
Similarly, the occurrence of an earthquake is well known to increase the likelihood of
another earthquake nearby in space and time. For example, we plot in Figure 1 a histogramof the times between “nearby earthquakes”, pairs of earthquake events separated in space by
110 kilometers or less, for all recorded earthquakes of magnitude 3.0 or greater in Southern
California during 2004-2005. The histogram shows a spike at short times, indicating an
increased likelihood of another event in the days following each earthquake. For a stationary
Poisson process the distribution of times between pairs of events would be approximately
uniform when the length of the time window is much larger than the longest time bin of thehistogram.
In the case of residential burglary, evidence indicates that an elevated risk exists for botha house that has been recently burgled and its neighboring houses (Farrell and Pease, 2001;
Short et al., 2009; Johnson et al., 2007). To illustrate this point further, we plot in Figure
1 a histogram of the times between “nearby burglaries”, residential burglaries separated inspace by 200 meters or less, for all recorded residential burglaries within an 18km by 18km
region of the San Fernando Valley in Los Angeles during 2004-2005. Again we observe a
spike at short times, indicating an increased likelihood of victimization within a few hundred
meters and several days of each burglary.
Self-excitation is also found in gang violence data, as an event involving rival gangs canlead to retaliatory acts of violence. In Figure 2, we plot the times of all recorded violent
crimes between the gang known as “Locke Street” and the rival gang known as “Lowell
Street” occurring between 2000 and 2002 in the Los Angeles police district of Hollenbeck.
Here we observe clear clustering patterns suggestive of self-excitation in the rate at which
2Figure 1: On the left, histogram of times (less than 300 days) between Southern California earthquake events of magnitude 3.0 or greater separated by 110 kilometers or less. On the right, histogram of times (less than 50 days) between burglary events separated by 200 meters or less.

the two rival gangs attack each other.

Figure 2: Times of violent crimes between two rivalry gangs in Los Angeles.
We propose that self-exciting point processes can be adapted for the purpose of crimemodeling and are well suited to capture the spatial-temporal clustering patterns observedin crime data. More specifically, spatial heterogeneity in crime rates can be treated using
background intensity estimation and the self-exciting effects detected in crime data can bemodeled with a variety of kernels developed for seismological applications or using non3

parametric methods. In Section 2, we review self-exciting point processes in the context ofseismological modeling. In Section 3, we present a model for residential burglary based upon
non-parametric methods for Epidemic Type Aftershock-Sequences models of earthquakes.
Our methodology combines the idea of stochastic declustering with Kernel Density Estimation in a novel way. In Section 4, we compare the predictive accuracy of our methodology
with prospective crime hotspot maps. The results illustrate how crime hotspot maps canbe improved using the self-exciting point process framework. We validate the methodology
with a simulated point process in the Appendix.

2

Self-exciting point process models in seismology

A space-time point process N (t, x, y) is typically characterized via its conditional intensity
λ(t, x, y), which may be defined as the limiting expected rate of the accumulation of points
around a particular spatial-temporal location, given the history Ht of all points up to timet (Daley and Vere-Jones, 2003):

E N {(t, t + ∆t) × (x, x + ∆x) × (y, y + ∆y)}|Ht.
λ(t, x, y) =lim∆t,∆x,∆y↓0∆t∆x∆y(1)

In Seismology a mark Mk , the magnitude of the earthquake, is associated with each event(tk , xk , yk ) and the conditional intensity often takes the form,

λ(t, x, y, M ) = j(M )λ(t, x, y),λ(t, x, y) = µ(x, y) +Xg(t − tk , x − xk , y − yk ; Mk ).

(2)

(3)

{k: tk <t}Models of this type, referred to as Epidemic Type Aftershock-Sequences (ETAS) models,
work by dividing earthquakes into two categories, background events and aftershock events.

4

Background events occur independently according to a stationary Poisson process µ(x, y),with magnitudes distributed independently of µ according to j(M ). Each of these earthquakes then elevates the risk of aftershocks and the elevated risk spreads in space and time
according to the kernel g(t, x, y, M ).
Many forms for g have been proposed in the literature, though in general the kernel ischosen such that the elevated risk increases with earthquake magnitude and decreases inspace and time away from each event. For example, the isotropic kernel,g(t, x, y; M ) =eα(M −M0 )K0·,(t + c)p (x2 + y 2 + d)q(4)

is one of a variety of kernels reviewed in (Ogata, 1998). Here K0 , M0 , and α are parametersthat control the number of aftershocks, c and d are parameters that control the behavior of
the kernel at the origin, and p and q are parameters that give the (power law) rate of decayof g.
Standard models for the background intensity µ(x, y) include spline, kernel smoothing,and Voronoi estimation (Ogata and Katsura, 1988; Silverman, 1986; Okabe et al., 2000). In
the case of fixed bandwidth kernel smoothing, the background intensity is estimated by,µ(x, y) = µ ·Xu(x − xk , y − yk ; σ),(5)kwhere µ is a parameter controlling the overall background rate. The events (tk , xk , yk , Mk )are assumed to be background events and in practice can be obtained through a declustering
algorithm (Zhuang, Ogata, and Vere-Jones, 2002).
The appropriate selection of parameter values is as critical to the modeling process as
specifying accurate forms for µ, g, and j. The distance in space and time over which therisk spreads, the percentage of background events vs. aftershocks, the dependence of the5increased risk on magnitude size, etc. all can have a great impact on the predictive power ofa point process model. Parameter selection for ETAS models is most commonly accomplishedthrough maximum likelihood estimation, where the log-likelihood function (Daley and VereJones, 2003),

l(θ) =X Z T ZZ
log λ(tk , xk , yk ; θ) −λ(t, x, y; θ)dydxdt,

k

0



(6)

S

is maximized over all parameter sets θ. Here S × [0, T ] is the space-time window of observation.
More recently, non-parametric methods have been introduced for self-exciting point process estimation (Zhuang, 2006; Marsan and Lenglin, 2008). Consider space-time point data{(tk , xk , yk )}N
k=1 and a general self-exciting point process model of the form,Xλ(t, x, y) = µ(t, x, y) +g(t − tk , x − xk , y − yk ).

(7){k: tk <t}Assuming model correctness, the probability that event i is a background event, pii , is given
by,
pii =

µ(ti , xi , yi ),λ(ti , xi , yi )(8)and the probability that event j triggered event i, pji , is given by,pji =g(ti − tj , xi − xj , yi − yj ),λ(ti , xi , yi )(9)(Zhuang, Ogata, and Vere-Jones, 2002). Let P denote the matrix with entries pji (note thatthe columns sum to one). Then stochastic declustering can be used in the following way.
Given an initial guess P0 of the matrix P , we then have N (N + 1)/2 probabilistic data points{(tk , xk , yk , pkk )}N
k=1 and {(ti − tj , xi − xj , yi − yj , pji )}i>j . Given this data, a non-parametric6density estimation procedure can be used to estimate µ from {(tk , xk , yk , pkk )}N
k=1 and gfrom {(ti − tj , xi − xj , yi − yj , pji )}i>j , providing estimates µ0 and g0 . We can then proceed
iteratively as follows until convergence is achieved,

Step 1) Estimate µn and gn from Pn−1 .
Step 2) Update Pn from µn and gn using (8) and (9).

For example, a simple histogram estimator is used in (Marsan and Lenglin, 2008) in step 1.

3

A self-exciting point process model of burglary

For the purpose of modeling burglary we consider an unmarked self-exciting model for theconditional intensity of the form,

λ(t, x, y) = ν(t)µ(x, y) +Xg(t − tk , x − xk , y − yk ).

(10){k: tk <t}Here we neglect spatially-localized temporal fluctuations in the background rate and assume
that the fluctuations occur globally (for example due to weather, seasonality, time of day,etc.) In the case of seismology, research over a number of decades was needed to refine the(parametric) form of the triggering function g. For this reason, non-parametric methods areappealing in the context of crime in order to quickly gain insight into the forms of ν, µ and g.
For this purpose we use the iterative procedure outlined in the previous section to estimate
the model, with several modifications.
Because the probabilistic data {(tk , xk , yk , pkk )}N
k=1 and {(ti − tj , xi − xj , yi − yj , pji )}i>j isboth 3-dimensional and the number of data points is O(N 2 ) (where N is typically O(1000)
for earthquake and crime data sets), the estimation step for µ and g is computationally7expensive. The dimensionality prevents straight forward use of binning methods such as the
Average Shifted Histogram (Marsan and Lenglin use a logarithmically scaled histogram ona coarse grid), as many bins may have extremely small, but non-zero, values (since the data
is probabilistic, the count in each bin can be less than 1). Alternatively, the large size of thedata set prevents efficient use of off-grid methods such as Kernel Density Estimation. To getaround these issues we use the following Monte-Carlo based iterative procedure,b
Step 1) Sample background events {(tbi , xbi , yib )}N
i=1 and offspring/parent inter-point disotances {(toi , xoi , yio )}N
i=1 from Pn−1 .

Step 2) Estimate νn , µn and gn from the sampled data using variable bandwidth Kernel
Density Estimation.
Step 3) Update Pn from νn , µn and gn using (8) and (9).

Because Nb + No = N , the size of the sampled data at each iteration allows for the useof Kernel Density Estimation. Another issue is that the number of background and offspring events, Nb and No , is changing at each iteration. Thus a fixed bandwidth for any
density estimation technique (kernel smoothing, histogram, etc.) will over smooth at someiterations and under smooth at others. Therefore we employ variable bandwidth KDE (alternatively Cross Validation could be used). We give further details of our approach andprovide validation using a simulated point process in the Appendix.

Results

We fit the model given by Equation (10) to a dataset collected by the Los Angeles Police
Department of 5376 reported residential burglaries in an 18km by 18km region of the San
Fernando Valley in Los Angeles occurring during the years 2004 and 2005. Each burglary is
8associated with a reported time window over which it could have occurred, often a few hourspan (for instance, the time span over which a victim was at work), and we define the timeof burglary to be the midpoint of each burglary window.
o
In Figure 3, we plot the sampled inter-point distances {(toi , xoi , yio )}N
i=1 on the 75th iter-

ation of the stochastic declustering algorithm (see Appendix for convergence verification).
The number of sampled (offspring) events is 706 (13.1% of all events) and of these events
approximately 63% are exact-repeats (occurring at the same house). On the left, the spatial
inter-point distances are plotted showing that elevated crime risk travels around 50m−100mfrom the house of an initial burglary to the location of direct offspring events. As discussedin (Marsan and Lenglin, 2008), the overall distance near-repeat risk travels is several times
further due to the cascading property of self-exciting point processes. Note also that therisk travels vertically and horizontally (along streets), more so than it does in other directions. On the right, we plot the spatial (x-coordinate) inter-point distances against the timeinter-point distances. Here exact-repeat burglaries, those occurring at the same house, areapparent along the x-axis.

Figure 3:

Spatial (left) and space-time (right) offspring/parent inter-point distanceso{(toi , xoi , yio )}N
i=1 sampled from P75 .

9

In Figure 4, we plot (on a logarithmic scale) the estimated marginals
Z Z
g75 (t) =g75 (t, x, y)dxdy
x

y

and
Z Z
g75 (t, x, y)dydt

g75 (x) =tycomputed from the KDE estimate of g at the 75th iteration. Here the presence of exactrepeat events can again be seen, as g75 (x) appears to approximate a delta distribution at theorigin. The spike around 1-2 days in the plot of g75 (t) is due to the presence of fast “crime
sprees,” where most likely the same burglar visited several neighboring houses within a timespan of a few minutes to several days. There are also several bumps in the elevated risk ofburglary, for example around 7 days. Here one possible explanation is that the routine of theburglar and/or the victim is such that a particular day of the week is a preferable time tocommit the burglary. After 7-10 days, the elevated risk of repeat/near-repeat victimization
drops to an intermediate level and stays relatively flat for a time span on the order of several
hundred days before decaying back to baseline rates. These results are consistent withprevious quantitative studies of exact-repeat burglaries (Short et al, 2009).

Figure 4: Marginal g75 (t) (left) and marginal g75 (x) (right) estimated using KDE based upon
offspring/parent inter-point distances sampled from P75 .

10In (Johnson et al, 2008), the authors discuss the need to model risk heterogeneity andin general it is a difficult task to separate clustering due to background heterogeneity andclustering due to self-excitation. One benefit of using the non-parametric approach outlined
above is that temporal and spatial changes in the rate of crime are automatically separated
into those stemming from exogenous effects and those due to self-excitation. In Figure 5, weoplot the estimated marginals ν75 (t) and µ75 (x, y) estimated using KDE from {(toi , xoi , yio )}N
i=1

at the 75th iteration. Here the estimated background rate exhibits temporal fluctuations on
a time scale of months/days, separate from the fluctuations due to self-excitation. Thesefluctuations are likely caused by a number of factors such as seasonal, economic, and demographic changes, as well as temporal variations in burglar routine activities (Felson, 1998).
For example, residential burglary tends to have a higher weekday rate (when victims are atwork) compared to weekends.
Similarly, the background rate is also spatially variable, which is consistent with fixed environmental heterogeneity in crime opportunities, as well as variability in population density
through space (Bernasco et al, 2005). In seismology, declustered catalogs are of great interestas they can be used in estimating the background rate of major earthquakes. Declusteredcrime catalogs could potentially be used by police to distinguish between areas of a city withintrinsically high crime rates and areas with temporarily high crime rates (due to near-repeat
effects). As the former arises due to structural properties of a given neighborhood and thelatter from behavioral characteristics of individual burglars, police and community responses
would likely need to be different in each case.

4

Crime forecasting: point processes vs. hotspot maps

Crime hotspot maps are a well established tool for visualization of space-time crime patterns
and can be used as a method for prediction of near-repeat crimes. Given space-time crime11Figure 5: Background rate time marginal ν75 (t) (left) and space marginal µ75 (x, y) (right) estimated using KDE from the background events sampled from P75 .

observations (tk , xk , yk ), crime hotspot maps are generated for a time interval [t − T, t] byoverlaying a density plot of the function,

λ(t, x, y) =Xg(t − tk , x − xk , y − yk ),(11)t−T <tk <tonto a city map, where g(t, x, y) is a space-time kernel. By flagging the areas of the citywhere λ takes on its highest values, crime hotspot maps can be used to indicate which areasin the city are likely to contain future crimes (Bowers et al., 2004; Chainey et al., 2008).
For example, in (Bowers et al., 2004) a coarse grained kernel is used that decays inversely
proportional to spatial and temporal distance. In particular, with spatial distance d in unitsof 1/2 cell widths and time t in units of weeks, the kernel in (11) is specified asg(t, d) =1(1 + t)(1 + d)

(12)

on the domain (t, d) ∈ [0, 2 months] × [0, 400 meters] and 0 otherwise. Such a crime hotspot
map is referred to as “prospective”, as it uses past crimes, coupled with the contagious spreadof crime (modeled by g), to estimate future relative crime risk across the city. It should be12noted that the risk is relative because (11) is not a point process intensity.
Here we compare the predictive accuracy of the self-exciting point process model developed in Section 3 to the prospective crime hotspot map given by (11)-(12). Because crime
is localized in small regions of the city (a commercial zone with no residential burglary may
be located 100 meters from a neighborhood), we find that for predictive purposes variablebandwidth KDE is less accurate than fixed bandwidth KDE. We therefore estimate µ(x, y)
in Equation (10) using fixed bandwidth Gaussian KDE, with 20-fold cross validation used
to select the bandwidth (Silverman, 1986).
For every day k of 2005, each model assesses the risk of burglary within each of M 2 cells
partitioning an 18km by 18km region of the San Fernando Valley in Los Angeles. Based on
the data from the beginning of 2004 up through day k, the N cells with the highest risk
(value of λ) are flagged yielding a prediction for day k + 1. The percentage of crimes fallingwithin the flagged cells on day k + 1 is then recorded and used to measure the accuracy ofeach model.
In Figure 6, on the left we plot the percentage of crimes predicted averaged over theforecasting year against the percentage of flagged cells for the self-exciting point process
and the prospective hotspot strategy. For example, with 10% of the city flagged the point
process and prospective hotspot correctly predict 660 and 547 crimes (respectively) out of2627. The difference in accuracy between the two methodologies can be attributed to thecrime hotspot maps failure to account for the background rate of crime. While prospective
crime hotspot maps used for crime prediction attempt to quantify the contagious spread ofcrime following past events, they fail to assess the likelihood of future “background” events,
the initial events that trigger crime clusters.
In order to disentangle the dependence of model accuracy on parameter selection, in
Figure 6 on the right we repeat the same prediction exercise but with parameters of eachmodel selected to yield the highest number of crimes predicted (L1 norm over 1 through 15
13

% of cells flagged). The optimal cutoff parameters for the prospective hotspot map are 200
meters and 39 weeks. With these parameter values, in particular the slow decay of g in time,
Equation (11) is closer to Poisson estimation. For the point process model we only optimize
the bandwidth used for µ(x, y) as the computational cost of the stochastic declustering algorithm is relatively high. Whereas the bandwidth is estimated to be approximately 300 meters
using cross validation, a smaller bandwidth, 130 meters, provides a higher level of predictive
accuracy. This can be attributed to the spatially localized features of neighborhoods, andhence burglary.
For all percentages of cells flagged the prospective hotspot map underperforms the point
process, though for certain percentages the relative underperformance is less. On the left
in Figure 6, the prospective hotspot map performs better (relative to the point process) forsmaller percentages of cells flagged, as the parameters are selected to account for near-repeat
effects. On the right, the prospective hotspot map performs better for larger percentages offlagged cells, since for these parameter values the model is more accurately estimating fixedenvironmental heterogeneity. For crime types such as robbery and auto theft, where nearrepeat effects play less of a role, prospective hotspot maps tailored for near-repeat effects
are likely to be outperformed by simple Poisson estimation. The advantage of models of
the form (10) is that the balance between exogenous and endogenous contributions to crime
rates is inferred from the data as opposed to being imposed apriori.

5

Discussion

We showed how self-exciting point processes from seismology can be used for the purpose ofcrime modeling. In the future it may be desirable to tailor point process models specificallyfor crime, taking into account the crime type and the local geography of the city. Based upon
the insights provided by nonparametric estimates, parametric models can be constructed that14

Figure 6: Forecasting strategy comparison. Average daily percentage of crimes predicted plottedagainst percentage of cells flagged for 2005 burglary using 200m by 200m cells. Error bars correspond to the standard error. Prospective hotspot cutoff parameters are 400 meters and 8 weeks(left) and optimal parameters (right) are 200 meters and 39 weeks. Spatial background intensityµ(x, y) smoothing bandwidth for the point process is 300 meters (left) selected by cross validationand 130 meters (right) selected to optimize the number of crimes predicted.

have advantages with respect to model fitting and simulation. Background rates can also be
improved by incorporating other data types (in Johnson et al, 2008 housing density is used
to improve models of repeat victimization). In the case of gang violence, a hybrid networkpoint process approach may be useful for capturing the self-exciting effects stemming fromgang retaliations. Here increased risk may not diffuse in geographic space, but instead may
travel through the network space of gang rivalry relations.
The methodology used in this study can be implemented for other applications as well,for example refining point process models of earthquakes. It could potentially be adapted,
more generally, to other second order models of point processes. The stochastic declustering
algorithm opens up the door to a plethora of density estimation techniques (Silverman, 1986;
Scott, 1992; Eggermont and LaRiccia, 2001) that could be used to explore point processes
in a way parametric methods do not allow.
In (Marsan and Lenglin, 2010) it is shown that the method is an Expectation-Maximization(EM) type algorithm. At the maximization step the complete data log-likelihood function

15

decouples in terms of the background and triggering functions, which is why at each iteration
the problem reduces to several decoupled density estimation problems. Several issues could
potentially arise here, one being that the method could converge to a local (but not global)
minimum of the observed data log-likelihood function. Another, as pointed out in (Sornette
and Utkin, 2009), is that the sample size and domain size (relative to the support of thetriggering kernel) play a key role in the accuracy of stochastic declustering. In numericaltests we have found that at least O(1000) data points are needed in 3 dimensions for theiterates to converge to the right densities and the domain needs to be several times larger
than the support of the triggering kernel. Similar to analytical results for standard density
estimation, it would be useful to have convergence results relating sample size, branchingratio, domain size, and the bandwidth of the density estimators to the solution of the fixed
point iteration.

6

Appendix

Given point data (tk , xk , yk )N
k=1 and a self-exciting point process model of the form,

λ(t, x, y) = ν(t)µ(x, y) +Xg(t − tk , x − xk , y − yk ),(13){k: tk <t}

we iterate the following until convergence:b
Step 1) Sample background events {(tbi , xbi , yib )}N
i=1 and offspring/parent inter-point disotances {(toi , xoi , yio )}N
i=1 from Pn−1 .

Step 2) Estimate νn , µn and gn from the sampled data.
Step 3) Update Pn from νn , µn and gn using (8) and (9).

In order to estimate νn , µn and gn from the sampled data, we use variable bandwidth16oKernel Density Estimation. To estimate gn , we first scale the data {(toi , xoi , yio )}N
i=1 to haveunit variance in each coordinate and based upon the rescaled data compute Di , the kthnearest neighbor distance (3-dimensional Euclidean distance) to data point i. We thentransform the data back to its original scale and, letting σx , σy , and σt be the sample
standard deviation of each coordinate, estimate the triggering function as,
N

o
1 X
1
exp
gn (t, x, y) =N i=1 σx σy σt (2π)(3/2) Di3(x − xoi )2 (y − yio )2 (t − toi )2−−−.
2σx2 Di2
2σy2 Di2
2σt2 Di2

The background rate is estimated similarly, where 1-dimensional and 2-dimenional Gaussian
kernels are used to estimate νn and µn respectively. In (Zhuang, Ogata, and Vere-Jones,
2002), the authors recommending using the 10th-100th nearest neighbor distance for Di .
Throughout we compute Di corresponding to ν using the 100th nearest neighbor distanceand in higher dimensions we use the 15th nearest neighbor distance for Di corresponding toµ and g.
We validate the method by simulating (13) withµ
ν(t)µ(x, y) =exp(2π)(4.5)2

x2−
2(4.5)2




exp

y2
−
2(4.5)2and
g(t, x, y) = θω exp






x2
y2
− ωt exp − 2 exp − 22σx2σyand comparing the estimates supplied by the method with the known distribution. The
simulation was carried out by first simulating all background events according to the Poissonprocess νµ. The rest of the simulation was carried out iteratively, where each point of eachgeneration generates its own offspring according to the Poisson process g centered at theparent point. The process terminates at the nth generation when all events of the nthgeneration lie outside of the time window under consideration. In order to have a realization17

of the point process at steady state, the first and last 2000 points were disregarded in eachsimulation.
In Figure (7), we plot the L2 error kPn − Pn−1 k2 at the nth iteration and the numberof sampled background events Nb at the nth iteration against the true number of background events for one realization of the known point process. Here we observe that theerror converges quickly for the first 10 iterations and then stabilizes as the error introducedby estimating the point process through sampling P cannot be reduced further (unless adeterministic iterative procedure is employed). We also verify that the method applied to
the 5376 burglary events in Section 3 reached convergence in Figure (7). Here we observe asimilar rate of convergence for the crime data as with the simulated point process.
In Table 1, we list the exact parameter values used for the simulated point process andthe estimates averaged over the final 10 iterations of the stochastic declustering algorithm
for each of 5 simulations of the point process. The parameter values were selected to yield
point patterns with scales similar to those observed in crime data. The parameter estimateso
are computed using the sample variances of the coordinates of {(toi , xoi , yio )}N
i=1 and the valuesof Nb , No . As some error is due to sample variation, we plot in the last two columns theestimated number of background events vs. the actual number of background events ineach of the 5 simulations to assess the ability of the method to reconstruct the realized
branching structure. In Figure (8), we plot the estimated marginals of g(t, x, y) against
the actual distributions on the 75th iteration of the stochastic declustering algorithm. The
estimated time marginal density deviates from the true density at the origin due to the jump
discontinuity of the exponential distribution. However, the estimate of the parameter ω isstill close to the true value (see Table 1).

18

Figure 7: L2 error kPn − Pn−1 k2 (top left) and Nb , the number of sampled background events,(top right) at the nth iteration for known point process model. L2 error kPn − Pn−1 k2 (bottom
left) and Nb , the number of sampled background events, (bottom right) at the nth iteration for the
method applied to the 5376 burglary events in Section 3.

References Cited[1] Bernasco, W. and Nieuwbeerta, P. (2005). How do residential burglars select target
areas? A new approach to the analysis of criminal location choice. Brit. J. Criminology
45, 296–315.
[2] Bowers, K. J., Johnson, S. D., and Pease, K. (2004). Prospective Hot-Spotting: The
Future of Crime Mapping? Brit. J. Criminol., 44, 641–658.
[3] Chainey, S., Tompson, L., and Uhlig S. (2008). The Utility of Hotspot Mapping for
Predicting Spatial Patterns of Crime, Security Journal, 21, 4–28.
19

Table 1: Parameter Value Estimates
True Values
Run 1 Est.
Run 2 Est.
Run 3 Est.
Run 4 Est.
Run 5 Est.

w−1
10.00
11.08
12.20
11.76
13.30
11.27

σx
0.01000.01760.01560.01500.01350.0147σy0.10000.14330.12960.12950.14070.1317θ
0.20000.20010.19670.19970.20490.2102µNb Est.
5.71005.6921 3999.7
5.7768 4016.5
5.6711 4001.5
5.6185 3975.3
5.7652 3948.9

Nb True
4041
4026
4017
4015
3977

[4] Cohen, J. and Tita, G. (1999). Spatial Diffusion in Homicide: Exploring a General
Method of Detecting Spatial Diffusion Processes. Journal of Quantitative Criminology
15 (4), 451–493.
[5] Daley, D. and Vere-Jones, D. (2003). An Introduction to the Theory of Point Processes,
2nd edition. New York: Springer.
[6] Eggermont, P. P. B. and LaRiccia, V. N. (2001). Maximum Penalized Likelihood Estimation: Volume I: Density Estimation, New York: Springer.
[7] Farrell, G. and Pease, K. (Ed.) (2001). Repeat Victimization, Criminal Justice Press,
New York.
[8] Felson, M. (1998). Crime and Everyday Life, Pine Forge Press, Thousand Oaks.
[9] Johnson, S. D. (2008). Repeat burglary victimisation: a tale of two theories. Journal of
Experimental Criminology 4, 215–240.
[10] Johnson, S. D., Bernasco, W., Bowers, K. J., Elffers, H., Ratcliffe, J., Rengert, G.,
and Townsley, M. (2007). Space-time patterns of risk: A cross national assessment ofresidential burglary victimization, J. Quant. Crim., 23, 201-219.

20

Figure 8: Estimated (circles) and actual (solid line) marginals of g(t, x, y) on the 75th iteration.
Upper left is the marginal g(x), upper right is the marginal g(y), and lower figure is marginal g(t).

[11] Marsan D. and Lenglin, O. (2008). Extending Earthquakes Reach Through Cascading,
Science, 319: 1076–1079.
[12] Marsan D. and Lenglin, O. (2010). A new estimation of the decay of aftershock density with distance to the mainshock. Preprint (http://www.lgit.univsavoie.fr/dmars/Publications.html.
[13] Ogata, Y. (1998). Space-time point process models for earthquake occurrences, Ann.
Inst. Statist. Math., 50 (2), 379–402.
[14] Ogata, Y. and Katsura, K. (1988). Likelihood analysis of spatial inhomogeneity formarked point patterns, Ann. Inst. Statist. Math., 40, 20–39.
21

[15] Okabe, A., Boots, B., Sugihara, K. and Chiu, S. (2000). Spatial Tessellations, 2nd ed.
Wiley, Chichester.
[16] Scott, D. W. (1992). Multivariate Density Estimation, New York: Wiley.
[17] Short, M.B. et al.. (2009). Measuring and Modeling Repeat and Near-Repeat Burglary
Effects, J. Quant. Criminol. 25 (3), 325.
[18] Silverman, B. W. (1986). Density Estimation for Statistics and Data Analysis, London:
Chapman and Hall.
[19] Sornette, D. and Utkin, S. (2009). Limits of declustering methods for disentangling
exogenous from endogenous events in time series with foreshocks, main shocks andaftershocks, Phys. Rev. E, 79, 061110.
[20] Tita, G. and Ridgeway, G. (2007). The Impact of Gang Formation on Local Patternsof Crime. Journal of Research on Crime and Delinquency 44 (2), 208–237.
[21] Zhuang, J., Ogata, Y., and Vere-Jones, D. (2002). Stochastic Declustering of SpaceTime Earthquake Occurences, JASA, 97 (458), 369–380.
[22] Zhuang, J. (2006). Second-order residual analysis of spatiotemporal point processes andapplications to model evaluation. Journal of the Royal Statistical Society: Series B, 68(4), 635–653.

22