On Lewis' Simulation Method for Point Processes

and T. Cover, " Some equivalences between Shannon entropy and Kolmogorov complexity, " IEEE Trans. Abstract-A simple and efficient method of simulation is discwsed for point processes that are specified by their conditional intensities. Tbe metbud is based on tbe thinning algorithm which was introduced recently by Lewis and Shedler for the simulation of nonhomogeneous Poisson pmcesses. Algorithms are given for past dependent point processes con-tain@ multivariate processes. The simulations are performed for some parametric conditional intensity functions, and the accuracy of tbe sbnu-lated data is demonstrated by the liieliiood ratio test aud the minbuum Akaiie information criterion (AK) procedure.


I. INTRODUCTION
A NY point process (N,, 4, P) on a finite interval (0, T] is a submartingale and therefore by the Doob-Meyer decomposition may be written as N, =mt + A,, where m, is an (F,, P) martingale and A, is the natural increasing process.It is known that there is a predictable process (A,, F,), such that A, = jdh, ds, if and only if P is Similar results hold for multivariate or marked point processes [7], [8].
The main object of this paper is to discuss the applications of Lewis' thinning simulation algorithm to any point process which is absolutely continuous with respect to the standard Poisson process.Recently Ozaki [ll] generated simulation data for Hawkes' self-exciting processes by making use of a recursive structure.However his method is not fast enough unless the process has a simple structure, because given a past history of the process t,, t,; .* 7 t, and a uniform random number U,, i from the interval (0, l), we have to solve the equation U,, i = S(tn+*ltl,.-*,t,) by Newton's iterative method to get the next point tn+l, where S is the conditional survivor function qt/t,;.*,tn)=exp -1 J 'X(slt,; . ., t,)ds .

f" I
We do not need to solve this equation to get the next point.The idea of simulating these point processes by thinning is developed using algorithms due to Lewis and Shedler [9] for the simulation of nonhomogeneous Poisson processes.
In Section II we give the simulation method and a proof for past dependent point processes containing multivariate processes.Some typical algorithms also will be given.In Section III we give some examples of parametric intensity functions for the simulation and obtain their 0018-9448/81/0100-0023$00.75 0 1981 IEEE maximum-likelihood estimates from the simulated data.
Proof Define a random measure for the finite marked The accuracy of the simulated data will be discussed by, process (using the notation in [7]) by the likelihood ratio test or the minimum Akaike information criterion (AIC) procedure.to the past history 4. The main result, which is formally similar to the one given in [9], is as follows.
Let t:<t,*<;.-,<t&be the points in (0, T] of the process (N*; F, P).Delete the points t; with probability 1 -h,,/Xt forj= 1,2;. ., N;.Then the remaining points { ti} form a point process (N, F, P) with conditional intensity X= {A,} in the interval (0, T]. It is readily seen from the predictability of X* that the constructions of A*, t;, and ti should be performed sequentially in the following manner.1) Suppose that the last point before time t has just been obtained.Then construct AT which is 41-measurable, piecewise constant, and A: > A,, for t > ti.
2) Simulate homogeneous Poisson points t,T( > ti) according to the intensity X:.
3) For each of the points {t;}, the probability XtY/h; is given conditionally independent of ty under the past history J$.
4) ti+, is the first accepted point among t; (> ti).Details of the algorithm will be given later.By generalizing this result to a multivariate point process, we have the following proposition.i) Z(t, 4)=8,(q) with probability X:/AT for p, q= 0,1,2; -* ) m, where 6,(q) is a Dirac delta function.
ii) For fixed t,N*(dt) and {Z(t,p),p=O,l;--,m} Since the predictable random measure corresponds uniquely to the multivariate process (see [7], [S]), this completes the proof. Cl We now give some typical algorithms based on this proposition.

3) Simulate stationary
Poisson processes for each interval of constant intensity.Denote the points by t;, t;; * * ) t;,.Consider the case of a univariate self-exciting process.Since E;; = a{N,, 0 <s < t}, the intensity of the process is given by a function of t and the points ti before t, i.e.A, =A(tlt,; -* , t,).There are two types of intensity processes.
Consider first the case where the path function of the intensity process is decreasing if no more points occur.The predictability of h, implies left-continuity of the path functions.
We assume that the minimum value of the intensity function is p, the jump size at each point is not larger than cy, and the AT are values, of the piecewise constant function such that h(t It,, * . ., t,) < AT for t, < si < t<s,+, =G t,,,.If ua < T then put t, =~a.Otherwise stop.
Set i=j=k=O and n= 1.
Set j equal to j+ 1 and generate q.

Generate a homogeneous
Poisson process with intensity Xl on the interval (kr, (k+ l)r].
If the number of the points on the interval, say NT, is zero, go to step 11.
Denote the ordered points on the interval (ir,(i+

Hawkes' Self-Exciting Process
The intensity function is given by h(t)=A(tlt,,--.,t,)=p+'v(t-s)dN(s), s 0 where p>O, V(S) > 0, and /,"v(s)ds< 1 for asymptotic stationarity of the process.Hawkes and Oakes [5] first gave the author the idea that this process may be simulated through nonhomogeneous Poisson processes.Indeed they say that a Hawkes' self-exciting process is nothing but an immigrant-birth process which is composed of an homogeneous Poisson immigrant with rate p and nonhomogeneous Poisson descendants with rate Y(S).As a parametrization of Y(S), Hawkes [4] used an exponential Y(s)=(Y~ -ps.In this case we can apply Algorithm 2 for the simulation. Ozaki and Akaike [12]  The jump size for this case is Xip_oa,?(j//?)je -j.
Suppose t,, t,; . ., t, in (0, T] are the simulated data. Then the log-likelihood function is given by where Rj(i) and Sj(t) are given recursively in the follow-intensity function ing way.

Linear Wold Process
The intensity function is given by

X(t)=p+a,(t-t(,))+ i ak(t(k-l)-t(k)), k=2
where p and (Y are nonnegative parameters and t(,) is the kth last point before t.A point process with this intensity is always asymptotically stationary (see [3]).It is easily seen that this point process may be simulated by the simple relation

Stress-Release Process
Vere-Jones [ 141 has suggested models for a series of strong earthquakes.
One of these models is defined by the

Nonlinear Hawkes' Type Point Process
Consider an intensity function of the form  where R,(i), R2(i), and R3(i) are given recursively by R,(l)=O, R,(l)=l,

Bivariate Wold Process
The intensity functions are given by The coefficients must be nonnegative for asymptotic stationarity if y, and p2 are strictly positive.If p, = 0, or p2 =O, then there are explosive cases even if the other coefficients are nonnegative.For example, consider the simplest casep=l.
The 'gradients vector and Hessian matrix are given similarly by a recursive formula.Also the simulation algorithm should be performed recursively for the greatest efficiency.

IV. CONCLUDING REMARKS
In this section we discuss whether the simulation data are statistically accurate enough in each case.In [lo] a collection of regularity conditions is given to prove the following.
which is asymptotically x:-distributed.See [13] for an 3) w-v) -wm is asymptotically xi-distributed extensive discussion of the relation between the minimum as T+m, where k is the dimension of the parameter 8.
AIC procedure and the likelihood ratio test.The examples in the preceding section basically satisfy Using physically generated random numbers we perthese conditions, although the multivariate case is not formed simulation experiments five times for each examtreated there.So the adoption of the minimum AIC proce-ple.The maximum-likelihood estimates and the negatives dure [l] is justified.That is to say, we consider two of the log-likelihoods are listed in the tables in the Apcompeting models ZZ, and Hi, where ZZ,, is the model pendix.We used the Davidon-Flecher-Powell method supposed to have the true parameter t9,, and H, is any for the nonlinear optimization.AIC, = ( -2)(maximized value of log-likelihood) + 2k.

ACKNOWLEDGMENT
The AIC is an estimator of the expected negative entropy The idea of the thinning method was communicated by which is a natural measure of the discrimination between Professor P. A. W. Lewis, who kindly corrected the earlier the true distribution for the data and the estimated proba-versions of this paper.For this the author is very grateful.bility distribution.Therefore, if the simulated data are The author also would like to thank G. Kitagawa for correctly distributed according to H,,, we can expect providing a program in [2] and for valuable suggestions AIC, <AIC,.Another useful method is to adopt the such as the recursive calculation technique, and K. Katsura likelihood ratio test of Ho against H,.Under regular for his generous help in testing the algorithms.Finally the situations such as the nested sequence of models, the author would like to thank Professor Okamoto for his relationship between the AIC and the likelihood ratio many critical comments.
Consider a point process (N, F, P) = {N,, I;;, 0 < t < T, P} on a fixed interval with its F-predictable intensity process X= {A,}, where F= {e} is a family of right continuously increasing u-fields.Suppose we obtain a positive F-predictable piecewise constant process A* = {A:} which is constructed pathwise in such a way that X, < XT, almost surely (a.s.), O< t < T. Then X* can be an intensity process of a locally homogeneous Poisson process (N*, F, P)= { NF, F,, P} with piecewise constant intensity changing its rate according F, P),p=1,2;.., m, on an interval (0, T] with joint intensity (X, F) = {AT, F,}, p= 1; . ., m. Suppose we can find a one-dimensional F-predictable process AT which is defined pathwise satisfying p=l P-almost surely, and set A"=+ 5 x;.p=l Let t:, t2*; . ., t& ~(0, T] be the points of the process (N*, F, P) with intensity process A:.For each of the points, attach a mark p =O, 1,.* ., m with probability x$/A:.. Then the points with marks p= 1,2,. . ., m, provide a'multivariate point process which is the same as that given above.result of the conditions of the proposition, Z(t, p) has the following properties.
From the tables we can see other model with a fixed dimension k of the parameter 0. that the maximum-likelihood estimates get more accurate The values of the AIC for the two models are as the sample size (number of points) or the length of the AIC, = ( -2)(value of log-likelihood at 0,)) observed interval increases.Also for each sample size since the number of unknown parameters in Ho is zero, interval length, the AIC and log-likelihood ratio tests and work well and justify the accuracy of the simulations.