simulation of inhomogeneous poisson process

NHPP.R

# Simulation of NHPP (CTDE style)
# based on: https://freakonometrics.hypotheses.org/724
# analytical functions from Mathematica

rm(list=ls());gc()

# intensity function
lambda <- function(x){
  100*(sin(x*pi)+1)
}

# cumulative (integrated) intensity
Lambda <- function(t){
  integrate(f=lambda,lower=0,upper=t)$value
}
Lambda <- Vectorize(Lambda,"t")

# analytic integrated intensity
Lambda1 <- function(t){
  (100 * (1 + (pi*t) - cos(pi*t)))/ pi
}

# CDF of time to next event
Ft <- function(t,x){
  1-exp(-Lambda(t+x)+Lambda(t))
}

# CDF of time to next event (Faster)
Ft1 <- function(t,x){
  1-exp(-Lambda1(t+x) + Lambda1(t))
}

# CDF of time to next event (fastest; all analytic)
Ft2 <- function(t,x){
  1 - exp(-(100 * ((pi*x) + cos(pi*t) - cos(pi*(t+x)) )) / pi)
}

# inverse CDF to sample the time to next event
Ftinv1 <- function(u,t,tmax){
  fn <- function(uu){
    Ft2(t,uu) - u
  }
  uniroot(f = fn,lower = 0,upper = tmax)$root
}

# simulation
Tmax <- 10
t <- 0
X <- t
while(tail(X,1) < Tmax){
  delta_t <- Ftinv1(u = runif(1),t = t,tmax = Tmax)
  t <- t + delta_t
  X <- c(X,t)
}

hist(X,breaks=seq(0,max(X)+1,by=.1),col="steelblue",
     main="NHPP Simulation",xlab="Time",ylab="Jump Density")
u=seq(0,max(X),by=.02)
lines(u,lambda(u)/10,lwd=2,col="firebrick3")