generate a realization from a transition matrix

Question

Consider a markov chain of 4 states $\{S_1, S_2, S_3, S_4\}$ described by the transition matrix

$$ A = \begin{bmatrix} .25 & .20 & .25 & .30 \\ .20 & .30 & .25 & .30 \\ .25 & .20 & .40 & .10 \\ .30 & .30 & .10 & .30 \end{bmatrix} $$

I understand that given an initial state $x^{(0)}$ the probability that the system is in state $i$ after $n$ steps is

$$ x^{(n)} = A^n x^{(0)} $$

However, I want to know how to generate a realization of a sequence states. For example, from the matrix above one possible sequence of states over $n = 0, 1, 2, 3$ steps could be

$$ S_1 \rightarrow S_3 \rightarrow S_4 \rightarrow S_1 $$

Is there a way to programmatically generate a particular sequence or realization from the transition matrix, either based on a random initial state or given an initial state $x^{(0)}$?

Answer 1

Starting from state $i$, draw a sample according to the distribution defined by that column and update the state according to the output.

Here's an implementation in Python using $0$-index.

import numpy as np


A = np.array([[.25, .2, .25, .3], [.2, .3, .25, .3], [.25, .2, .4, .1], [.3, .3, .1, .3]])
xk = np.arange(len(A))

def generate_sample(cur_state):
    return np.random.choice(xk, 1,  p = A[:, cur_state])

initial_state = 0
sample_len = 100
output = [-1 for i in range(sample_len)]
output[0] = initial_state
for i in range(1, sample_len):
    output[i] = generate_sample(output[i-1])[0]
print(output)

Remark: In most literature, each row sums to $1$ but for your example, each column sums to $1$.