Simulate Compound poisson process in r

Question

I'm trying to simulate a compound Poisson process in r. The process is defined by $ \sum_{j=1}^{N_t} Y_j $ where $Y_n$ is i.i.d sequence independent $N(0,1)$ values and $N_t$ is a Poisson process with parameter $1$. I'm trying to simulate this in r without luck. I have an algorithm to compute this as follows: Simutale the cPp from 0 to T:

Initiate: $ k = 0 $

Repeat while $\sum_{i=1}^k T_i < T$

Set $k = k+1$

Simulate $T_k \sim exp(\lambda)$ (in my case $\lambda = 1$)

Simulate $Y_k \sim N(0,1)$ (This is just a special case, I would like to be able to change this to any distribution)

The trajectory is given by $X_t = \sum_{j=1}^{N_t} Y_j $ where $N(t) = sup(k : \sum_{i=1}^k T_i \leq t )$

Can someone help me to simulate this in r so that I can plot the process? I have tried, but can't get it done.

Are you able to make use of `rnorm`, `rexp,` and `while`? It may be slow, but undifferent from any other programming language. What have you tried? — AdamO, Feb 09 '17 at 16:10
Hi, I have the same question here, how would you simulate Y from any distribution? — G-09, Jul 24 '21 at 08:22

whuber · Accepted Answer · 2017-02-09T16:46:39.770

Use cumsum for the cumulative sums that determine the times N_t as well as the X_t. This illustrative code specifies the number of times to simulate, n, simulates the times in n.t and the values in x, and (to display what it has done) plots the trajectory.

n <- 1e2
n.t <- cumsum(rexp(n))
x <- c(0,cumsum(rnorm(n)))
plot(stepfun(n.t, x), xlab="t", ylab="X")

This algorithm, since it relies on low-level optimized functions, is fast: the six-year-old system I tested it on will generate over three million (time, value) pairs per second.

That's usually good enough for simulation, but it doesn't quite satisfy the problem, which asks to generate a simulation out to time T. We can leverage the preceding code, but the solution is a little trickier. It computes a reasonable upper limit on how many times will occur in the Poisson process before time T. It generates the inter-arrival times. This is wrapped in a loop that will repeat the procedure in the (rare) event the time T is not actually reached.

The additional complexity doesn't change the asymptotic calculation time.

T <- 1e2            # Specify the end time
T.max <- 0          # Last time encountered
n.t <- numeric(0)   # Inter-arrival times
while (T.max < T) {
  #
  # Estimate how many random values to generate before exceeding T.
  #
  T.remaining <- T - T.max
  n <- ceiling(T.remaining + 3*sqrt(T.remaining))
  #
  # Continue the Poisson process.
  #
  n.new <- rexp(n)
  n.t <- c(n.t, n.new)
  T.max <- T.max + sum(n.new)
}
#
# Sum the inter-arrival times and cut them off after time T.
#
n.t <- cumsum(n.t)
n.t <- n.t[n.t <= T]
#
# Generate the iid random values and accumulate their sums.
#
x <- c(0,cumsum(rnorm(length(n.t))))
#
# Display the result.
#
plot(stepfun(n.t, x), xlab="t", ylab="X", sub=paste("n =", length(n.t)))

Thanks a lot, this helps me tremendously. Is it possible to get a nicer plot? i.e. one that does not display a circle around every point? — Slangers, Feb 09 '17 at 16:52
Yes: see the help page found at `?plot.stepfun`. It tells you to include the argument `do.points=FALSE` when you invoke `plot` on the last line and it describes other options for controlling the appearance of the plot. — whuber, Feb 09 '17 at 16:54

Simulate Compound poisson process in r

1 Answers1