In general, Mathematica is not always the best tool for this type of simulation. That said, it is not necessary to calculate sampling paths to calculate the mean or the variance.

For a Markov chain with probability transition matrix $ P_ {i, j}: = mathbb {P} (X_1 = j | X_0 = i) $, you can calculate the $ n $ -step transition matrix with $ mathbb {P} (X_n = j | X_0 = i) = P ^ n $. Therefore, given the initial distribution $ x $, the probability of being in a given state is given by the corresponding input of $ xP ^ n $.

Given a discrete random variable $ Y $, we can calculate the expectation using the formula,

$$ mathbb {E}[Y] = sum_ {k = 0} ^ { infty} k mathbb {P} (Y = k) $$

In your case, we will calculate $ mathbb {E}[xP^n]$ for each $ n $. Since you have only three states, there will only be three entries in this sum (for $ k = $ 0, $ k = $ 3 and $ k = $ 1).

In your case, with your weighting of the state space, you can do it like this:

```
ListPlot @ Table[Dot[{1, 0, 0}.MatrixPower[P, n], {0, 3, 1}], {n, 0, 500}]
```

Similarly, the variance is computed by $ mathbb {E}[(xP^n)^2]- mathbb {E}[xP^n]^ 2 $. You can calculate this by:

```
ListPlot @ Table[
Dot[{1, 0, 0}.MatrixPower[P, n], {0, 3, 1} ^ 2 -
Point[{1, 0, 0}.MatrixPower[P, n], {0, 3, 1}]^ 2], {n, 0, 500}]
```

These two solutions could be made more efficient by using the previous iteration of the matrix power to calculate the next, but since $ P $ is small, these commands are executed instantly and the efficiency is therefore not not a problem.

Regarding your code:

It looks like you're calculating an average of 500 traces and then drawing 100. Do you want to calculate these 100 quantities or did you want one? When I run your code my parcel of `bad`

only goes up to a little over 2.5. I do not know why your plot has so many.

You also use the replacement after calculating the string. You can switch the inputs of the probability transition matrix in advance to avoid this and save a little calculation.