Simulate a composite Poisson process in r

Question

I am trying to simulate a complex Poisson process in r. The process is defined by $ \ sum_ {j = 1} ^ {N_t} Y_j $, where $ Y_n $ is iid sequence-independent values ​​$ N (0,1) $, and $ N_t $ is a Poisson process with parameter $ 1 $. I am trying to mimic this in r with no luck. I have an algorithm to calculate this: Simutale 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 simulate this in r so that I can build a process? I tried, but can't do it.

Answer 1

Use cumsumfor cumulative amounts defining the times N_t as well as X_t. This illustrative code indicates the number of times to simulate n,, simulates the time in n.tand the values ​​in x, and (to show what it did) determines the path.

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, works quickly: the six-year-old system I tested on will generate more than three million (time, value) pairs per second.

---

, , T. , . , T. . , () , T .

.

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)))