I came across the following probability problem:
There are four states, named: S, A, B, C.
Transition probability is given by:
$$
\begin{matrix}
& S & A & B & C \\
S & 0 & \frac{1}{3} & \frac{1}{3} & \frac{1}{3} \\
A & 0 & 0 & \frac{1}{2} & \frac{1}{2} \\
B & 0 & \frac{1}{2} & 0 & \frac{1}{2} \\
C & 0 & \frac{1}{2} & \frac{1}{2} & 0 \\
\end{matrix}
$$
The system starts in state $S$ and is run for a total of $t$ steps. What is the probability distribution of spending a certain amount of time $u$ in state $A$ (By symmetry the same probability would be true for the other states, $B$ and $C$)?
Originally, I was only interested in small numbers for the value of $t$. E.g. for the first four:
$$
\begin{array}{|c|c|c|}
\hline
t/p & 1 & 2 & 3 & 4 \\
\hline
0 & \frac{2}{3} & \frac{1}{3} & \frac{1}{6} & \frac{1}{12} \\
1 & \frac{1}{3} & \frac{2}{3} & \frac{2}{3} & \frac{1}{2} \\
2 & 0 & 0 & \frac{1}{6} & \frac{5}{12} \\
\hline
\end{array}
$$
A useful check for this table is that:
$$
\sum_{j=0}^{\lceil\frac{1}{2}n\rceil} ju_j(t) = \frac{t}{3}
$$
as generally the chance for A to be the state at any random time past the starting state should be $\frac{1}{3}$, as the matrix is both symmetric and has steady state of each of A, B, C being equal in probability
The first two case, $t=1$ is trivial. For $t=2$, observe that the only time the second value can be $A$ is if the first is not, so you can just add the two probabilities, $\frac{2}{3} \cdot \frac{1}{2} + \frac{1}{3} = \frac{2}{3}$. For $t=3$ and $t=4$, I wrote out tables of all possible orders, which are equal in their probability of occurrence
$$P = \frac{1}{3} \cdot (\frac{1}{2})^t = \frac{1}{2} u_0(t)$$
This gives me also the result for $u_0(t)$ as there are only two patterns that can lead to this case, $BCBCBC\dots$ and $CBCBCB\dots$.
I wonder if there's a more general way of computing the table further out with more than exponential efficiency though, as with each increase in $t$ by one, the amount of possible orders doubles, and brute force won't go much further than 20-30 using a computer.
Best Answer
A dynamic programming approach can compute this in polynomial time, even for arbitrary Markov chains. The key to avoid the exponential blow-up you encountered with your approach is to use the Markov property, i.e., that you do not need to look at the entire initial run fragment of the Markov Chain to have full information about the coninuation of the run, but only the last state.
For a set $X$ let $Dist(X)$ denote the set of probaility distributions over $X$. Let $M = (S,i,T)$ be a Markov chain with state space $S$, initial distribution $i \in Dist(S)$ and transition matrix $T$. We define a function $f \colon S \times S \times \mathbb{N} \rightarrow Dist(\mathbb{N})$ where $f(s_0,s_g,t)$ is the probability distribution over the number of times visiting $s_g$ in a run of time $t$ starting in $s_0$. The probability of visiting $s_g$ exactly $k$ times can then be expressed as $f(s_0,s_g,t)(k)$
The property you're interesed in, i.e., the probability distribution over the times visiting a state $s_g$ in a finite run of the Markov chain of length $t$, which I'll denote as $d(M,s_g,t)$ can then be defined for every $k\in \mathbb{N}$ as follows:
$$ d(M,s_g,t)(k) = \sum_{s\in S} i(s) f(s,s_g,t)(k) $$
Claim: $f$ (and thus $d$) is computable in polynomial time in $\lvert S \rvert$ and $t$.
Let $eq \colon S \times S \rightarrow \{0,1\}$ be defined such that $eq(s_1,s_2)=1$ iff $s_1 = s_2$.
We define $f$ inducdively on $t$.
For the base case, define $f(s,s_g,0)(1) = eq(s,s_g)$ and $f(s,s_g,0)(0) = 1-eq(s,s_g)$ (and consequently $0$ for all other inputs), i.e., the probability distribution puts all probability on one either $0$ or $1$, depending on whether the goal state $s_g$ is equal to the starting state.
For the recursive case define
$$f(s_0,s_g,t+1)(k) = \sum_{s_1\in S} T(s_0,s_1) \cdot f(s_1,s_g,t)(k-eq(s_0,s_g))$$
To see that this is correct, notice that we here split up any possible run $s_0s_1\dots s_{t+1}$ starting in $s$ and performing $t+1$ transitions into the segments $s_0$ and $s_1\dots s_{t+1}$. Notice that for the latter the distribution of visiting $s_g$ can be expressed as $f(s_1,s_g,t)$. We then weight this function by the probability of transitioning from $s$ to $s'$ and "shift" the arguments of the distribution by $1$ in case $s_0=s_g$ so that the count remains accurate for the latter run fragment.
To see that this is polynomially time-computable, notice that for any $k>t$ we clearly have $f(s_0,s_g,t)(k) = 0$. This means that we only have to compute $f(s_0,s_g,t)(k)$ for inputs in the set $S\times \{ s_g \}\times \{ 0,\dots,t \}\times \{ 0,\dots,t \} $ which has size in $O(\lvert S \rvert t^2)$.
Naively, one might think that the recursion is exponential since you have up to $\lvert S \rvert$ branching cases, but they just reduce to the same cases. E.g. in your example if we choose $t=4$ we need to compute $f(S,A,4)$ which in turn needs to compute $f(A,A,3)$, $f(B,A,3)$ and $f(C,A,3)$, i.e., has 3 branches. Now to compute $f(A,A,3)$ we need to compute $f(B,A,2)$ and $f(C,A,2)$ - two more branches. But note that $f(B,A,2)$ is also required to compute $f(C,A,3)$. Similarly, $f(C,A,2)$ is also required to compute $f(B,A,3)$. So while there are 6 cases from the $3*2$ branches corresponding to all 6 possible 2-step runs (and generally exponentially many), they reduce to the same subproblems.