Equivalence of two state Markov chain and sampling via geometric distribution Let $\mathcal T = \{1,2,\ldots,T\}$ denote the set of points in time, $S = \{0,1\}$ the state space, $X = (X_t)_{t \in \mathcal T} \in S^\mathcal T$ a time series, $\alpha = \mathbb P(X_{t+1} = 0 \mid X_t = 0) \in [0, 1]$ the stationary probability that the state remains $0$, given it was $0$ in the previous period, and $\beta = \mathbb P(X_{t+1} = 1 \mid X_t = 1) \in [0, 1]$ the stationary probability that the state remains $1$, given it was $1$ in the previous period, respectively. The transition matrix of the Markov chain is then given by
\begin{align}
P = 
\begin{pmatrix}
\alpha & 1-\alpha \newline
1-\beta & \beta
\end{pmatrix}.
\end{align}
Let $\mathcal I = \{1,2,\ldots,n\}$ denote an index set. Assume we've have generated $n$ time series $(X^i)_{i \in \mathcal I}$ and let $\overline X(n) = \sum_{i \in \mathcal I}{X^i} / n$ denote the average time series.
I was thinking that we can generate $\overline X(n)$ via an alternative approach. Note the following: Given the state is $0$, the probability that it is going to be $1$ after $t$ periods is $f(t) = \alpha^{t-1}(1-\alpha)$. And given the state is $1$, the probability that it is going to be $0$ after $t$ periods is $g(t) = \beta^{t-1}(1-\beta)$. Now construct the time series $Y = (Y_t)_{t \in \mathcal T}$ as follows: Initialize $Y_0 \in S$. If $Y_0 = 0$, draw via $f(t)$ for how long it is going to stay $0$. Then draw via $g(t)$ for how long it its going to stay $1$ etc. until $t = T$. If $Y_0 = 1$ start with drawing from $g(t)$, then $f(t)$ etc. Generate $(Y^i)_{i \in \mathcal I}$ time series and let $\overline Y(n) = \sum_{i \in \mathcal I}{Y^i} / n$ denote the average time series.
Problem
Is $\lim_{n \to \infty} \overline X(n) = \lim_{n \to \infty} \overline Y(n)$ true?
Python
I implemented the idea in Python with $(\alpha, \beta, T, n) = (0.50, 0.75, 100, 10'000)$ and it seems that the two approaches are not equivalent. In the figure I plotted the average of $n = 10'000$ simulated time series. Is my logic flawed? Or is the implementation incorrect?

The code to produce the figure reads as follows:
import numpy as np
import matplotlib.pyplot as plt

# parameter
a = 0.5  # Pr(X(t+1) = 0 | X(t) = 0)
b = 0.75  # Pr(X(t+1) = 1 | X(t) = 1)
P = np.array([[a, 1-a], [1-b, b]])  # transition matrix
n = 10000  # number of time series
tf = 100  # length of time series
X = np.zeros((n, tf), dtype=int)  # Storage Matrix
Y = np.zeros((n, tf), dtype=int)  # Storage Matrix

# Simulation
# Markov-chain
X[:, 0] = np.random.choice(2, n, p=[0.5, 0.5])  # initialize
for i in range(n):  # loop over time series
    for t in range(tf-1):  # loop over time within time series
        if X[i, t] == 0:
            X[i, t+1] = np.random.choice(2, 1, p=P[0, :])
        else:
            X[i, t+1] = np.random.choice(2, 1, p=P[1, :])

# average
X_ave = X.sum(axis=0) / n

# geometric distribution
Y[:, 0] = np.random.choice(2, n, p=[0.5, 0.5])  # initialize
for i in range(n):  # loop over time series
    t = 0
    while t < tf:  # loop over time within time series
        if Y[i, t] == 0:
            z = np.random.geometric(1-a)
            if t+z < tf:
                Y[i, t+1: t+z-1] = 0
                Y[i, t+z] = 1
                t = t+z
            else:
                Y[i, t: tf] = 0
                t = tf
        else:
            z = np.random.geometric(1-b)
            if t+z < tf:
                Y[i, t+1: t+z-1] = 1
                Y[i, t+z] = 0
                t = t+z
            else:
                Y[i, t: tf] = 1
                t = tf

# average
Y_ave = Y.sum(axis=0) / n


# plot
fig, ax = plt.subplots()
ax.plot(range(tf), X_ave, label='Markov')
ax.plot(range(tf), Y_ave, label='geometric')
ax.set(xlabel='t', ylabel='X(t)/n, Y(t)/n')
plt.title(f'number of time series {n}')
plt.legend()
plt.show()

Edit I
There was mistake in the Python code due to (in my opinion) confusing accessing. It needs to read Y[i, t+1: t+z] and Y[i, t+1: tf+1] instead of Y[i, t+1: t+z-1] and Y[i, t+1: tf].
Edit II
There was another mistake. For z=1 we would not repeat the state. It therefore needs to read Y[i, t: t+z] and Y[i, t: tf+1] instead of Y[i, t+1: t+z] and Y[i, t+1: tf+1]. We should further note that Y[t]=Y[t: t+1].

 A: 
According to the documentation for numpy.random.geometric,

The geometric distribution models the number of trials that must be
run in order to achieve success. It is therefore supported on the
positive integers, k = 1, 2, ....

this means that z-1 should used as the number of repetitions of the current state, before moving to the other state. I thus wonder at the use of Y[i, t+1: t+z-1] when z=1 and at the replacement by Y[i, t+1: t+z] which seems to always include a repetition of the current state, while it should not when z=1.


Since the Markov chain is a sequence of 0 and 1, as eg
0100100010111010111001
updating the Markov chain one position at a time or updating the uninterrupted blocks of 0 and 1 all at once are equivalent. As noted in the question, when at a state 0 at time t, the number of subsequent 0 till the next 1 is a indeed Geometric random variable $\mathcal Geo(1-\alpha)$ [when defined as taking values in $\mathbb N$, i.e. zero is a possible value]. Hence generating using Geometric batch lengths rather than going one time step at a time is equivalent from a probability perspective (if not producing the same output when implemented on a computer).
I ran both ways of simulating the Markov chain in R and did not spot any discrepancy in the ergodic convergence of the chain average to the stationary mean $\frac{1-\alpha}{1-\alpha+1-\beta}$ (equal to $1/3$ with the chosen values):
T=1e5
a=.5;b=.75
x=1:T
for(t in 2:T)
  x[t]=(runif(1)>ifelse(x[t-1],1-b,a))
plot(cumsum(x)/(1:T),type="l",col="gold")
t=1
while(t<=T){
  y=rgeom(1,ifelse(x[t],1-b,1-a)) #warnin: number of failures
  x[t:(t+y)]=x[t];x[t<-t+y+1]=!x[t]}
lines(cumsum(x[1:T])/(1:T),lty=2,col="indianred")

