Simulation of Gamma process (distribution of increments)

The gamma process is a Levy process $$X$$, where $$X_t$$ has gamma distribution with parameters $$at,b>0$$ and density $$f\left(x\right)=\frac{b^{at}}{\Gamma\left(at\right)}x^{at-1}e^{-bx}$$

I want to simulate gamma process by increments but what is the distribution of $$X_t - X_s$$? Of course gamma but with what parameters?

Suppose that we are given a probability space $$(\Omega, \mathcal{F}, P)$$. A Lévy process $$X = (X (t), t \geq 0)$$ taking values in $$\mathbb{R}^d$$ is essentially a stochastic process having stationary and independent increments; we always assume that $$X (0) = 0$$ with probability 1. So:
• each $$X (t) : \Omega \to \mathbb{R}^d$$;
• given any selection of distinct time-points $$0 \leq t_1 < t_2 < \ldots < t_n$$, the random vectors $$X(t_1), X(t_2) − X(t_1), X(t_3) − X(t_2), \ldots, X (t_n) − X(t_{n−1})$$ are all independent;
• given any two distinct times $$0 \leq s < t < \infty$$, the probability distribution of $$X(t) − X(s)$$ coincides with that of $$X(t − s)$$.
The Gamma distribution is scale invariant under summation, i.e. $$\sum_{i=1}^N X_i = \mathrm{Gamma}\left(\sum_{i=1}^N k_i, \theta\right)$$ so thanks to the third property $$X_t - X_s$$ has parameters $$a (t - s), b$$.