Consider a Poisson process , with non-homogeneous intensity
. Here, we consider a deterministic function, not a stochastic intensity. Define the cumulated intensity
in the sense that the number of events that occurred between time and
is a random variable that is Poisson distributed with parameter
.
For example, consider here a cyclical Poisson process, with intensity
lambda=function(x) 100*(sin(x*pi)+1)
To compute the cumulated intensity, consider a very general function
Lambda=function(t) integrate(f=lambda,lower=0,upper=t)$value
The idea is to generate a Poisson process on a finite interval .
The first code is based on a proposition from Çinlar (1975),
- start with
- generate
- set
- set
denote
- deliver
- go to step 2.
In order to get the infinimum of , consider a code as
v=seq(0,Tmax,length=1000) t=min(v[which(Vectorize(Lambda)(v)>=s)])
(it might not be very efficient…. but it should work). Here, the code to generate that Poisson process is
s=0; v=seq(0,Tmax,length=1000) X=numeric(0) while(X[length(X)]<=Tmax){ u=runif(1) s=s-log(u) t=min(v[which(Vectorize(Lambda)(v)>=s)]) X=c(X,t) }
Here, we get the following histogram,
hist(X,breaks=seq(0,max(X)+1,by=.1),col="yellow") u=seq(0,max(X),by=.02) lines(u,lambda(u)/10,lwd=2,col="red")

Consider now another strategy. The idea is to use the conditional distribution before the next event, given that one occurred at time ,
- start with
- generate
- set
- deliver
- go to step 2.
Here the algorithm is simple. For the computational side, at each step, we have to compute and then
. To do so, since
is increasing with values in
, we can use a dichotomic algorithm,
Ft=function(x) 1-exp(-Lambda(t+x)+Lambda(t)) Ftinv=function(u){ a=0 b=Tmax for(j in 1:20){ if(Ft((a+b)/2)<=u){binf=(a+b)/2;bsup=b} if(Ft((a+b)/2)>=u){bsup=(a+b)/2;binf=a} a=binf b=bsup } return((a+b)/2) }
Here the code is the following
t=0; X=t while(X[length(X)]<=Tmax){ Ft=function(x) 1-exp(-Lambda(t+x)+Lambda(t)) Ftinv=function(u){ a=0 b=Tmax for(j in 1:20){ if(Ft((a+b)/2)<=u){binf=(a+b)/2;bsup=b} if(Ft((a+b)/2)>=u){bsup=(a+b)/2;binf=a} a=binf b=bsup } return((a+b)/2) } x=Ftinv(runif(1)) t=t+x X=c(X,t) }

The third code is based on a classical algorithm to generate an homogeneous Poisson process on a finite interval: first, we generate the number of events, then, we draw uniform variates, and we sort them. Here, the strategy is closed, except that is won’t be uniform any longer.
- generate the number of events on the time interval
- generate independently
where
- set
i.e. the ordered values
- deliver
‘s
This algorithm is extremely simple, and also very fast. This is one function to inverse, and it is not in the loop,
n=rpois(1,Lambda(Tmax)) Ft=function(x) Lambda(x)/Lambda(Tmax) Ftinv=function(u){ a=0 b=Tmax for(j in 1:20){ if(Ft((a+b)/2)<=u){binf=(a+b)/2;bsup=b} if(Ft((a+b)/2)>=u){bsup=(a+b)/2;binf=a} a=binf b=bsup } return((a+b)/2) } X0=rep(NA,n) for(i in 1:n){ X0[i]=Ftinv(runif(1)) } X=sort(X0)
Here is the associated histogram,

An alternative is based on a rejection technique. Actually, it was the algorithm mentioned a few years ago on this blog (well, the previous one). Here, we need an upper bound for the intensity, so that computations might be much faster. Here, consider
- start with
- generate
- set
- generate
(independent of
)
- if
then deliver
- go to step 2.
Here, consider a constant upper bound,
lambdau=function(t) 200 Lambdau=function(t) lambdau(t)*t
The code to generate a Poisson process is
t=0 X=numeric(0) while(X[length(X)]<=Tmax){ u=runif(1) t=t-log(u)/lambdau if(runif(1)<=lambda(t)/lambdau) X=c(X,t) }
The histogram is here

Finally, the last one is also based on a rejection technique, mixed with the second one. I.e. define
The good thing is that this function can easily be inverted
- start (as usual) with
- generate
- set
- generate
- if
then deliver
- goto step 2.
Here, the algorithm is simply
t=0 while(X[length(X)]<=Tmax){ Ftinvu=function(u) -log(1-x)/lambdau x=Ftinvu(runif(1)) t=t+x if(runif(1)<=lambda(t+x)/lambdau(t+x)) X=c(X,t) }

Obviously those five codes work, the first one being much slower than the other three. But it might be because my strategy to seek the infimum is not great. And the latter worked well since there were not much rejection, I guess it can be worst…
All those algorithms were mentioned in a nice survey written by Raghu Pasupathy and can be downloaded from http://web.ics.purdue.edu/~pasupath/…. In the paper, non-homogeneous spatial Poisson processes are also mentioned…






soit stationnaire, de moyenne
. Alors en prenant l’espérance de part et d’autre dans l’expression précédente,
, i.e.
, ou
. Plus généralement, avec un ARMA,
, en prenant la encore l’espérance des deux cotés, on en déduit
.



(afin de se débarrasser de cette tendance). Alors


ou encore
suit alors un processus ARIMA(1,1,0) sans constante cette fois.
avec



dans le modèle de
(qui vont correspondre aux durées entre sauts), et éventuellement aussi des amplitudes de sauts (e.g. des pertes qui suivent des lois exponentielles). On a ici
. On commence par générer
en notant que
sont Gaussiens (centrés et de variance
) et indépendants les uns des autres. Quant aux durées entre sauts, les
. Voilà le code qui permet de générer les trois composantes,

est la date où survient le
ième saut,
, alors conditionnellement au fait que
, les variables
correspondent aux statistiques d’ordres de
, i.e.
,
)



Via 




Yesterday, I wrote a post (
But as @
Now both stocks are expressed in the same currency. To compare returns volatility, a first idea can be to use GARCH models,






est une suite de variables i.i.d. et que l’on définie les durées entre records,
, avec
. Alfred Rényi a montré que


, avec les dates des catastrophes, et
leur coût.




.










Sur les
, correspondant aux dates d’arrivée des sinistres, ou bien par
.
sera un processus de Poisson homogène d’intensité
si les durées entre arrivées,
i.i.d. de loi exponentielle de paramètre
est alors
. La première idée pour simuler un tel processus est de supposer qu’il existe
.
qui vaut
