Solved – Non homogenous Poisson process with simple rates

poisson processrsimulation

I am trying to stimulate number of claims in the next 12 months using a non-homogeneous poisson process. The rates are:

11.02 per day, during March, April and May
11.68 per day, during June, July and August
26.41 per day, during September, October, November
20.83 per day, during December, Jan and Feb

I came across the inversion method, which seems the simplest for simple rates like in my case (compared to the thinning method).

I tried to follow the algorithm:

s=0; v=seq(0,Tmax,length=1000)
   X=numeric(0)
   while(X[length(X)]<=Tmax){
     u=runif(1)
     s=s-log(u)
     t=min(v[which(Vectorize(Lambda)(v)>=s)])
     X=c(X,t)
   }

but how do I set my lambda to change depending on the month??

Best Answer

Your rate function,$\,\lambda(t)$, will look something like the image below if I understand what you're asking. To do the procedure I think you're referencing, you need the rate function $\lambda(t)$ and the cumulative rate function, $\Lambda(t) = \int_0^t \lambda(s)ds$. RateFunction

Generate a Homogeneous Poisson Process (NHPP) via Nonlinear Time Transformation
Task: generate $\{N(t),t\ge 0\}$ which is an NHPP with rate $\lambda(t)$
Given: target rate function $\lambda(t)$.

Let $N_0(t)$ be a rate-1 Poisson process (i.e. $\lambda = 1$). Let $\text{E}[N(t)] = m(t)$. Then $m(t) = \int_0^t \lambda(s)ds$.

Let $N(t) = N_0(m(t))$. Then the counting process $\{N(t),t\ge 0\}$ is a NHPP with rate $\lambda(t)$.

Procedure
1. Generate arrival times $S_1^{(0)},S_2^{(0)},S_3^{(0)},\ldots$ for a rate-1 Poisson Process with interarrival times $X_1,X_2,X_3,\ldots \sim \text{Exponential}(1)$.
2. Use the nonlinear time transformation to obtain event times: $0\le \tilde S_1 \le \tilde S_2 \le \tilde S_3 \le \cdots$ where $\tilde S_n = m^{-1}\left(S_n^{(0)}\right)$.
3. Construct the counting process $N(t) \equiv \text{max}\{n\ge 0:\tilde S_n \le t\}$.

If I can find time, I'll post some working code. Hope this helps.

I bootlegged the rate function like so (as an example only, the total number of days is off). All my implementations of this have been in MATLAB. 

Rate = [0;ones(90,1)*11.02;
        ones(90,1)*11.68;
        ones(90,1)*26.41;
        ones(90,1)*20.83];

Update: Added Thinning Method w/ short comparison here.

Related Question