Let $X=(X_{t},t \in T)$ be a non-homogeneous, continuous time Markov process with a finite state space S={1,...,K}.

Let $\alpha_{i,j}(t)$ be the hazard rates of some $\varGamma$-distributed random variables.

$\textbf{Question}:$ How can I simulate random paths of the Markov process with a transition intensity matrix which is built with the above hazard rates $\alpha_{i,j}(t)$?

For example: S={1,2,3,4} with given $\alpha_{i,j}(t), i,j \in S$

Any idea is greatly appreciated, with best regards, Wolfgang