I am reading about Markov chain and I understand how to find stationary distribution of a Markov chain and the transition probability matrix at some time t. But what I fail to understand how can one use it to get a chain of observations.
So let's say I have 2 different states (s1,s2) some initial distribution and some transition probability matrix for these two states. And I want to output a string of states(s1, s2) of some desired length. Is there a formula for it? How can do that?
So, with the setup you have described, you can easily implement this using a simple matrix equation. For the following discussion, I'm going to assume that your transition matrix is time homogeneous, meaning that the transition matrix does not change over time. I'm also going to assume that we have a discrete-time process, as you describe this as a Markov chain.
Now, we can call your initial distribution $\pi_{0} = [\pi_{0}(s_{1}), \pi_{0}(s_{2})]$. The entries of the matrix of transition probabilities is usually described as a function of the current state, $x$, and the next state of interest, $y$:
$P(x, y) = \Bigg[\matrix{ P(s_{1}, s_{1}) & P(s_{1}, s_{2})\\ P(s_{2}, s_{1}) & P(s_{2}, s_{2})}\Bigg]$
Now, to perform the state transition, you simply perform matrix multiplication on your $\pi_{0}$ vector by the transition matrix:
$\pi_{1} = \pi_{0}\cdot P(x, y)$
where $\pi_{1}$ is the probability distribution of the two states at time $t=1$. This procedure can be iterated to produce the distribution over your two states at arbitrary finite time $t$: $\pi_{t} = \pi_{0} \cdot P^{t}(x, y)$, where $P^{t}(x, y)$ is the transition matrix raised to the $t$ power.
The description above should be easy enough to implement as a simple computer program in say Python. If you were interested in actual realizations, meaning your output string should be the specific state that is actually reached by the Markov chain at each time step, then you would just modify the above procedure to perform a random number draw at each time step. The random number draw determines which state is actually reached- say your current distribution at time $t$ is $[ 0.35, 0.65]$. You generate a random number $r$ in the interval $[0, 1]$, and if $r \leq 0.35$ then you set the next state to $1$. Then your distribution at the next time step $t + 1$ would be, for example if state $1$ were the actually currently occupied state, $\pi_{t + 1} = [1, 0]$. Perform the matrix multiplication, which here is just extracting the top row of your transition matrix, and generate your next random number $r$. Repeat as many times as required by your string's desired length. Now, if the chain is actually time inhomogeneous, meaning that the transition matrix can change at each iteration step, you would further have to modify the procedure so that the entries of the matrix $P(x, y)$ properly reflected this time dependence at each step.
