3
$\begingroup$

I am trying to simulate the compound Poisson process using the next algorithm that I found in a textbook on stochastic processes.

  1. Let $S_0 = 0$.
  2. Generate i.i.d. exponential random variables $X_1, X_2, \ldots$
  3. Let $S_n=X_1+\cdots + X_n,$ for $n = 1, 2, \ldots$
  4. For each $k = 0, 1, \ldots,$ let $N_t = k,$ for $S_k\le t\le S_{k+1}$
S <- vector(mode="integer", length=100)
S[1] = 0
## Generation of Exponential random variables with parameter lambda
X <- rexp(n=100, rate=0.1)

for(n in 1:100){
  S[n] = sum(X[1:n])
}    

But, I am not clear about how to write the step $4$, maybe I need to put the integer between two $S_k$ (the arrival times)? I am interested in the counting process $N_t$.

Furthermore, how do you plot it?

I have something like that, but I am not clear if the arrival times are ok, because, in the first plot on the left, I see that every arrival time is not necessarily in the vertical line of the process, and on the right all the arrival times have the same length, I think that it contradicts the independence of the arrival times.

nro<-10
S<-vector(mode="integer",length = nro)
S[1]=0

##Generation of Exponential random variables with parameter lambda
X<-rexp(n =nro,rate = 2)
for(n in 1:nro)
S[n]=sum(X[1:n])
S<-cumsum(X)
n_func <- function(t, S) sapply(t, function(t) sum(S < t))
t_series <- seq(0, max(S), by = max(S)/nro)

#Plot of the trajectory and add lines in the arrival times
par(mfrow=c(1,2))
plot(t_series, n_func(t_series, S),type = 
"s",ylab=expression(N[t]),xlab="t",las=1,cex.lab=0.8,main="Poisson 
   Process",cex.axis=0.8)
grid()
abline(v = S,col="red",lty=2)
plot(t_series, n_func(t_series, S),type = 
"s",ylab=expression(N[t]),xlab="t",las=1,cex.lab=0.8,main="Poisson 
 Process",cex.axis=0.8)
grid()
abline(v = t_series,col="blue",lty=4)

enter image description here

$\endgroup$
  • 1
    $\begingroup$ examine the functions in this carefully: n=5;plot(stepfun(cumsum(rexp(n,.1)),0:n)) ... Is there a statistical question here or is this just about the R implementation? If it's just how to do something in R it would likely be off topic for stats.stackexchange.com ... but if there is a stats issue you should edit to clarify that. $\endgroup$ – Glen_b Oct 19 '17 at 0:55
  • $\begingroup$ @Glen_b, although the question is clearly bound up w/ the code, I interpret this as a manifestation of not fully understanding the underlying statistical ideas. $\endgroup$ – gung Oct 27 '17 at 0:35
2
$\begingroup$

You want to get a function of $t$ that gives the count of events. So simply do

n_func <- function(t, S) sapply(t, function(t) sum(S <= t))
t_series <- seq(0, max(S), by = max(S)/100)
plot(t_series, n_func(t_series, S)

$S$ is basically the time stamps of each Poisson events in your sample. So you want to just count the number of events that have time stamped before $t$.

$\endgroup$
  • $\begingroup$ I am confused with something here. If S is the time stamps i.e the time arrivals, at the end of the code If I add vertical lines with abline(v=S), I think that I should get vertical lines in every arrival time, but I don´t get it. $\endgroup$ – Boris Oct 26 '17 at 23:26
2
$\begingroup$

The Poisson process is describing how events occur over time. The inter-arrival times (i.e., the times between each successive event) are distributed as an exponential. In the algorithm provided, $S_n$ is the total amount of time that has elapsed since you started recording until the $n$th event occurred. Since your $X_i$s are the times between each event, you sum them from the start up to the event of interest to get the total elapsed time until that event. The $N_t$s are the number of events that have occurred up to any given point in time. Thus, the plot would be a set of horizontal lines at successive integers, where the x-axis was time.

Here is a simple example, coded in R:

set.seed(3535)             # this makes the example exactly reproducible
                           # there's no need for an explicit s[0]
X = rexp(n=100, rate=0.1)  # step 2
S = cumsum(X)              # step 3
N = 0:100                  # step 4

windows()
  plot(stepfun(x=S, y=N), main="Counting Process",
       xlab="Time", ylab="Number of events")

enter image description here

$\endgroup$
  • $\begingroup$ I relate it in the insurance context, the Counting process $N_t$ gives me the number of claims occurred in $[0,t]$, Is possible to get something like that? If I call $N$ in the code above, I got, numbers from 0 to 100, in the insurance context, for instance, it means that the number of accidents every month was increasing by $1$. It seems no realistic. Could you explain how to interpret correctly? Maybe in the first code, I got in the counting process, numbers as $0,0,1,1,2,2,3..$ It means that sometimes there are no claims, it seems more realistic. I don´t understand how to explain it. $\endgroup$ – Boris Oct 27 '17 at 0:48
  • 1
    $\begingroup$ @Boris, N is just the number of accidents, incrementing by one accident at a time, from the time you start counting (0 accidents) until you stop counting. It can be nothing other than the integers up through the total when you stop (or the current moment). The number of accidents over time is a 2-dimensional value, the current count, $N_t$, paired with the current time $t$. What you are after isn't in either variable alone, but in their relation. $\endgroup$ – gung Oct 27 '17 at 0:57
  • 2
    $\begingroup$ @Boris the $t$ in that $N_t$ is an instant in time ($N_t$ being the number of events up to that time). For a process where events happen independently in continuous time the process can only increase by $1$ at any given instant (otherwise it remains unchanged). If you want how many events happen in a month that would be found by computing $N_t - N_{t-s}$ where $s$ corresponds to one month of time in whatever units $t$ is in (e.g. 30 if $t$ is in days, 1/12 if $t$ is in years). If you have a process where two events really do occur in exactly the same instant in time, it's not a Poisson process $\endgroup$ – Glen_b Oct 27 '17 at 1:09
  • 1
    $\begingroup$ @Boris, time is considered continuous, so if measured properly, there cannot be 2 claims at the same time. They may only differ by 1 second, or 1/1000000 of a second, but they are not simultaneous. If you are asking about discrete time (eg, only keeping track of days), you need to work with a modified process. $\endgroup$ – gung Oct 27 '17 at 1:11
  • 1
    $\begingroup$ @Boris Notice that for two events to happen together it would require that the exponential variate that represents the time between the two events be exactly 0. This is an event whose probability is zero. $\endgroup$ – Glen_b Oct 27 '17 at 1:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.