# Manually simulating Poisson Process in R

The following problem tells us to generate a Poisson process step by step from $$\rho$$ (inter-arrival time), and $$\tau$$ (arrival time).

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

• Let $$τ_0 = 0$$.
• Generate i.i.d. exponential random variables $$ρ_1, ρ_2, . . ..$$
• Let $$τ_n = ρ_1 + . . . + ρ_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 $$\rho$$ 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 $$\tau$$s by progressive summing of $$\rho$$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?

Also, I am not sure about my Ntk() function.

Your code already fails within the first function as rexp(0.5) does not generate a realization of an exponential random variable. In fact, the first argument of rexp() is the number of realizations you want to generate (n) followed by the parameter $$\lambda$$ (rate) as the second argument. So there is no need to define the function rhos() in the first place. In general, your code is quite cumbersome and redundant.

You may simply generate the arrival times as follows

cumsum(rexp(n, lambda)).

or alternatively using the inversion method

cumsum(-log(runif(n))/lambda)

for some specified n and lambda.

• You may simply generate the arrival times as follows- I must to go step by step from inter-arrival-time to arrival-time. This has to be implemented manually. So, I think, cumsum is not allowed here. – user366312 Apr 26 '19 at 12:51
• No one prevents you from doing it with a loop as you did. However, note in this case that your function taos() is incorrect as the first (inter-)arrival time $\tau_1$ is not contained in your vector. – Alex Blasberg Apr 26 '19 at 13:46