I need to simulate a compound Poisson Process in R, however I am not clear with the algorithm to generate it. I have conceptual gaps.

I know by definition that: A compound Poisson process is the process: $$Z_t=\sum_{i=0}^{N_t}X_i$$ where $X_i$ are i.i.d random variables with a specific distribution. In my case I consider Weibull random variables with shape parameter $\tau$ and scale parameter $k$.

According to my self-study I know that $N_t$ should be a Poisson Process. Also, I know according to the theory that, the interarrival times of a Poisson process follows a exponential distribution with parameter $\lambda$. With these facts, my algorithm proposed is the next one.

i)$t=0$ and generate a random number that follows a Poisson distribution with parameter $\lambda t$

ii) Generate $n$ Weibull random variables $X_i$, $i=0,...,n$, with $X_0=0$

iii)Generate $n$ interarrival times ($E_i$) exponential random variables with rate $lambda$.

iv) Generate the vector $T=(0,E_1,E_1+E_2,....,E_1+E_2+...+E_n)$, in this way I would get the arrival times.

v) Now, for every $t\in T$ I can generate the Variable $Z_t=\displaystyle\sum_{i=0}^{N_t} X_i$

vi) $t=t+1$ and start again.

but, I am confused. First, If I generate at the beginning the Poisson random number, according this algorithm for every step I would get different $X_i$ that is not that I want, because I need the same $X_i$ (in the context that I am study it represents random losses)

Also, If I follow this algorithm, I think that for every step I would get different interarrival times and also this is not desired.

I have seen some algorithms on internet deffining a tmax and then generate the poisson random number, but I am very confused.

The lambda of the interarrival times is the same of the Poisson random variable? I am not clear with it too.

arrivaltimes <- c (0); # Array of arrival times
cumtime <- rexp (1, lambda_arr ); # Time of next arrival
level <- c (0); # Level of the compound Poisson process
while ( cumtime < maxtime )
 arrivaltimes <- c( arrivaltimes , cumtime );

 # Draw a new jump from the Weibull dsitribution distribution
 # and add it to the level .
 oldLevel <- level [ length ( level )];
 newLevel <- oldLevel + rweibull(1 ,scale = scale_loss,shape=shape_loss);
 level <- c(level , newLevel );

 # Draw a new interarrival time and add it
 # to the cumulative arrival time
 cumtime <- cumtime + rexp (1, lambda_arr );
 Z<-level# level is the vector that is the sum of random loss(Weibull)

1 Answer 1


You've sort of mixed up two distinct approaches to this problem, hence your understandable confusion.

One approach is to ignore the time dimension and simply generate a sample of $Z_i, i \in \{1, \dots, I\}$ random variables that have the appropriate distribution. This can be done quite easily by:

  1. Generate $I$ Poisson variates $N_i$ from the Poisson$(\lambda)$ distribution.
  2. For $i = 1, \dots, I$, generate $Z_i$ as the sum of $N_i$ random numbers, possibly zero, from the appropriately-parameterized Weibull distribution.

If you care about the time dimension, you can follow the second approach:

  1. Set $t = 0$.
  2. Generate an Exponential$(1/\lambda)$ interarrival time $\Delta t$.
  3. Generate $X_{t+\Delta t}$ from the appropriately-parameterized Weibull distribution
  4. Set $t := t + \Delta t$
  5. If not done, go to step 2, otherwise exit.

This will give you a sample from a Poisson process with jumps at the various (random) $t$ of (random) size determined by the Weibull distribution. If you pick an arbitrary time $\tau$, the distribution of the sum of $X_{t \leq \tau}$ is a compound Poisson with Poisson parameter $\lambda \tau$.

EDIT in response to questions in comments:

Now the question arises, how large a sample should you take? That depends upon the details of the problem you face. If, for example, you have a fixed time $T$ that you are supposed to generate a collection $\{X_t\}_{t \leq T}$ for, the most straightforward way to do this with the second approach is to stop after step 2 whenever $t + \Delta t > T$ (because your next $X_t$ would be for a $t > T$). With the first approach, if you have $T$ time intervals of length $1$ and you are interested in generating $Z_t$ for each time interval, you'd just set $I = T$ and proceed as above. If you do the latter, you can skip generating the $X_t$ altogether.

Second EDIT:

Given that the Weibull distribution is continuous rather than discrete, there is no counting process in the usual sense. Of course, the number of "arrivals" / Poisson events is still a counting process. For more on Poisson processes, you can check out the Wikipedia page https://en.wikipedia.org/wiki/Poisson_point_process.

  • $\begingroup$ but using the second approach, What value of the counting process I should take? It is ok to have as input $n:$ the number of interarrival times?, or better take a maximum value of time and find the number $n$ that follows a Poisson distribution, and then generate that number of interarrival times?, What is exactly $n$ this should be an input or should be generated as a Poisson random variable? $\endgroup$
    – Boris
    Nov 29, 2017 at 20:46
  • $\begingroup$ For instance I found in this post stackoverflow.com/questions/45084537/compound-poisson-process that at some point use : n = qpois(1-eps, lambda = lambda * t) but I don´t have clear why?. $\endgroup$
    – Boris
    Nov 29, 2017 at 20:46
  • $\begingroup$ but why is called Compound Poisson process? According to the algorithms that I have read and the algorithm that you suggest me, Only is important to take in account that the interarrival times follows exponential distribution,but where is the Poisson distribution or something related? $\endgroup$
    – Boris
    Nov 29, 2017 at 23:02
  • 1
    $\begingroup$ The Exponential distribution is the interarrival time distribution of a Poisson process. The Poisson distribution is that of the number of arrivals that occur in a unit of time given that the interarrival time distribution is Exponential. It's a compound Poisson process because the size of the "jumps" at each arrival is not always equal to 1 but has a distribution in its own right, so the sum of all of the "jumps" over a given time period isn't distributed Poisson even though the number of jumps is. $\endgroup$
    – jbowman
    Nov 29, 2017 at 23:24
  • $\begingroup$ But in the second approach that you suggested me, you never generates a Poisson random variable, only the exponentials? do you mean that in the Poisson process this is the condition? I mean that to generate the Poisson process is sufficient to generate the interarrival times that follows a Exponential distribution? $\endgroup$
    – Boris
    Nov 29, 2017 at 23:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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