Manually simulating Poisson Process in R

The following problem tells us to generate a Poisson process step by step from ρ (inter-arrival time), and τ (arrival time).

One of the theoretical results presented in the lectures gives the following direct method for simulating Poisson process:

• Let τ₀ = 0.
• Generate i.i.d. exponential random variables ρ1, ρ2, . . ..
• Let τ_n = ρ₁ + . . . + ρ_n for n = 1, 2, . . . .
• For each k = 0, 1, . . ., let N_t = k for τ_k ≤ t < τ_k+1.

1. Using this method, generate a realization of a Poisson process (N_t)_t with λ = 0.5 on the interval [0, 20].
2. Generate 10000 realizations of a Poisson process (N_t)_t with λ = 0.5 and use your results to estimate E(N_t) and Var(N_t). Compare the estimates with the theoretical values.

My attempted solution:

First, I have generated the values of ρ using rexp() function in R.

rhos <-function(lambda, max1)
{
    vec <- vector()

    for (i in 1:max1) 
    {
        vec[i] <- rexp(0.5)
    }

    return (vec)
}

then, I created τs by progressive summing of ρs.

taos <- function(lambda, max)
{
    rho_vec <- rhos(lambda, max)
    #print(rho_vec)

    vec <- vector()
    vec[1] <- 0
    sum <- 0
    for(i in 2:max)
    {
        sum <- sum + rho_vec[i]
        vec[i] <- sum
    }

    return (vec)
}

The following function is for finding the value of N_t=k when the value of k is given. Say, it is 7, etc.

Ntk <- function(lambda, max, k)
{
    tao_vec <- taos(lambda, max)
    val <- max(tao_vec[tao_vec < k])
}

y <- taos(0.5, 20)
x <- seq(0, 20-1, by=1)

plot(x,y, type="s")

Output:

As you can see, the plot of the Poisson process is blank rather than a staircase.

If I change rexp to exp, I get the following output:

.. which is a staircase function but all steps are equal.

Why is my source code not producing the expected output?

It looks like you're using max1 to indicate how many times to sample the exponential distribution in your rhos function. I would recommend something like this:

rhosGen <- function(lambda, maxTime){
  rhos <- NULL
  i <- 1
  while(sum(rhos) < maxTime){
    samp <- rexp(n = 1, rate = lambda)
    rhos[i] <- samp
    i <- i+1
  }
  return(head(rhos, -1))
}

This will continue to sample from the exponential until the sum of these holding times is larger than the length of the given interval. head the removes the last sample so that all of the events that we keep track of definitely occur in our time interval of interest. From here you have to generate the taos by summing the previous holding times (rhos):

taosGen <- function(lambda, maxTime){
  rhos <- rhosGen(lambda, maxTime)
  taos <- NULL
  cumSum <- 0
  for(i in 1:length(rhos)){
    taos[i] <- sum(rhos[1:i])
  }
  return(taos)
}

Now that you have the taos we know at what time each event in the time interval (0,maxTime) occurs. This leads us to generating the associated Poisson Process by finding the value of the Nt for each t in the time interval:

ppGen <- function(lambda, maxTime){
  taos <- taosGen(lambda, maxTime)
  pp <- NULL
  for(i in 1:maxTime){
    pp[i] <- sum(taos <= i)
  }
  return(pp)
}

This generates the value of the Poisson Process at each integer time in the interval. I suspect that part of your issue was trying to put the tao values on the y-axis instead of the count of events that had occurred already. The following code worked for me to produce a random looking stair case, similar to your example.

y <- ppGen(0.5, 20)
x <- seq(0, 20-1, by=1)

plot(x,y, type="s")