Gillespie Stochastic Simulation Algorithm
Algorithm 1
For $\mathbf{x}_{t} = i$, $$ P \left( \mathbf{X}_{t + \tau} = j | \mathbf{x}_{t} = i \right) = a_{j} \left( \mathbf{x}_{t} \right) \qquad , j = 1, \cdots , n $$ the method for simulating the realization of a continuous Markov chain $\left\{ \mathbf{X}_{t} \right\}$ is called the Gillespie Stochastic Simulation Algorithm, or simply SSA.
Input
Receives initial values $t=0$ and $\mathbf{X}_{0} = \mathbf{x}$.
Step 1.
Calculate $\sum_{j} a_{j} \left( \mathbf{x}_{t} \right)$.
Step 2.
Sample $\tau$ from the exponential distribution $\exp \left( \sum_{j} a_{j} \left( \mathbf{x}_{t} \right) \right)$. The larger $\sum_{j} a_{j} \left( \mathbf{x}_{t} \right)$ is, the smaller the expected value of $\tau$ becomes, and events occur more frequently.
Step 3.
Select $j$ with the probability of $\displaystyle {{a_{j}} \over {\sum_{j} a_{j} \left( \mathbf{x}_{t} \right)}}$. The larger $a_{j}$ is, the higher the probability that event $j$ will occur.
Step 4.
Record $\left( t + \tau, \mathbf{x}_{t + \tau} \right)$ and update as $$ \begin{align*} t \gets & t + \tau \\ \mathbf{x}_{t} \gets & \mathbf{x}_{t + \tau} \end{align*} $$. Since $\left\{ \mathbf{X}_{t} \right\}$ is a continuous-time Markov chain, one can return to Step 1 and repeat. Stop the algorithm when enough data has been obtained.
Output
Obtains the realization of given $\left\{ \mathbf{X}_{t} \right\}$.
Explanation
Example
The algorithm might look complicated with all the general expressions, but it’s actually pretty simple.
Consider a simple Birth-Death Process. This process starts with $N$ individuals and shows the increase or decrease in population according to the birth rate $b > 0$ and the death rate $d > 0$. Whether $b$ is greater or smaller than $d$ will determine the pattern, but the frequency of events itself depends on their absolute values. It is plausible that events, whether births or deaths, happen more frequently in proportion to the current number of individuals, and the interval at which events occur when the population is $i$ follows $$ \tau \sim \exp \left( (b+d) i \right) $$. Once the timing of an event is decided, we choose which event, birth or death, will occur. At $t$, if the population is $i$, to breed at $t + \tau$ means $j = i+1$, and to die means $j = i-1$. These are proportional to the rate $b, d$ and the population at $t$, $i$, respectively, so $$ a_{i+1} (x) = b i \\ a_{i-1} (x) = d i $$. According to Step 3, decide which event will occur with the probability calculated by normalizing these.
Code
The following is the Julia code that directly translates the above example. $$ b = 2 \\ d = 1 $$ With the given parameters, repeat the simulation ten thousand times starting with the initial population $N=2$ and compute the average at $t = 2$ to compare with the theoretical prediction. It is written to be easy to read rather than efficient, so let’s check if the algorithm is implemented as described in the example.
using Distributions
using Plots
b = 2.0
d = 1.0
θ = b/(b+d)
result = plot(xlims = (0,2), ylims = (0, Inf))
Y = Int64[]
@time for _ in 1:10000
y = [2]
t = [0.0]
while t[end] < 2 && y[end] > 0
push!(t, t[end] + rand(Exponential(1/((b+d) * y[end]))))
if rand() < θ
push!(y, y[end] + 1)
else
push!(y, y[end] - 1)
end
end
if t[end] > 2
push!(Y, y[end-1])
else
push!(Y, 0)
end
plot!(result, t,y, label = :none, alpha = 0.01, color = :black)
end
plot!(result, 0:0.1:2,2exp.((b-d) .* 0:0.1:2), label = :none, color = :red, linewidth = 4, linestyle = :dot)
result
abs(mean(Y) - 2exp.(2(b-d)))
Visualizing the result looks like the above. The difference from the theoretical prediction $2\exp(2(b-d)) = 14.78$ is about 0.1, confirming that it has been calculated quite accurately.