I've asked this question over in cross-validated over here:https://stats.stackexchange.com/q/610492/383970 But no answers. I was curious if somebody in this stack has an answer.
I've been learning Gillespie's algorithm to simulate continuous time Markov chains. I understand how the algorithm is derived from the reaction probability density function
$P(\tau, \mu)$ = probability at time $t$ that the event $\mu$ will occur in the time interval $(t + \tau, t + \tau + d\tau)$
What I don't get is how this guarantees us that we're sampling from the "master equation" or the Kolmogorov forward equations. Why is that sampling with this probability gives us trajectories that follow the distribution given by solving the forward equations?