Probability of a Markov chain in a finite set (2) This question is linked to the one here : Probability of markov chain in a finite set . My aim is to generalize it when we have $K$ states and compute it in the fastest way. Let's reintroduce the problem.
Let's $X$ be an homogeneous Markov chain with $K$ states $\{1,2,\dots,K\}$ et denote $(\pi_1,\,\pi_2,\,\dots,\,\pi_K)$ the initial probabilities and $P=(p_{ij})_{1\leq i,j\leq K}$ the transition matrix.
Let's consider a realization of $X$ in $T$ time step, meaning $X_1,\,X_2,\,\dots,\,X_T$. My aim is to compute for any $T$, the following probabilities  $P(N_1^T=n_1, N_2^T=n_2,\dots,N_{K-1}=n_{K-1})$ for any $\{(n_1,n_2,\dots,n_{K-1})\in(0,1,\dots,T)^{K-1} : n_1+n_2+\dots+n_{K-1}\leq T\}$, where $N_j^T$ ($j\in\{1,2,\,\dots,K-1\})$ is the number of time the Markov chain was at state $j$ between time $1$ to time $T$.
To simplify the notation, suppose $N_{1:j}^T$ is the vector $(N_1^T, N_2^T,\dots, N_{j})$ and $n_{1:j}$ is the vector $(n_1,\,n_2,\,\dots,\,n_{j})$. Here is what I proposed to resolve the problem. First, using Base formula, we have $$P(N_{1:(K-1)}^T=n_{1:(K-1)}) = P(N_1^T=n_1)P(N_2^T=n_2\mid N_1^T=n_1)\dots P(N_{K-1}^T=n_{K-1}\mid N_{1:(K-2)}=n_{1:(K-2)})\,.$$ My idea is to compute each term of this product using the Chapman-Kolmogorov equation : $$P(u) = 
 P^u\,,\quad \text{where}\quad P(u) = (p_{ij}(u))_{1\leq i,j\leq K} = (P[X_{t+u} = j \mid X_t=i])_{1\leq i,j\leq K}\,.$$ To do so I define the probability of the first time the state is $k$ and all the previous state between time $1$ to $T$ are greater than $j$ conditionnally on $X_0 = i$ by \begin{equation}
    f_{ik}^{(j)}(T) = P\left(X_T = k, \,X_{T-1}> j,\, \dots, \,X_{1}> j\;\middle\vert\;X_0=i\right)
\end{equation}and have shown that \begin{align}
    f_{ik}^{(j)}(1) &= P\left(X_{1} = k \mid X_0 = i\right) = p_{ik}\,,\\
    f_{ik}^{(j)}(T) &= p_{ik}(T) - \sum_{u=1}^{T-1}\sum_{l=1}^{j}f_{il}^{(j)}(u)\,p_{lk}(T-u)\,,\qquad  T=2,\,3,\,\dots\,.
\end{align}
I've also shown some recursive relationship for $P\left(N_{1:j}^{T} = \mathbf{0} \;\middle\vert\; X_0 = i\right)$, $P\left(N_1^{T} = n_1 \;\middle\vert\; X_0 = i \right)$, $P\left(N_j^{T} = n_1,\,N_{1:(j-1)}^{T} = \mathbf{0} \;\middle\vert\; X_0 = i \right)$ for $2\leq j\leq K-1$ and $1\leq i\leq K$. By exemple
$$P\left(N_j^{T} = n_1,\,N_{1:(j-1)}^{T} = \mathbf{0} \;\middle\vert\; X_0 = i \right) = \left\{\begin{array}{lcl}
\sum_{u=1}^{T-1} f_{ij}^{(j)}(u)P\left(N_{1:j}^{T-u} = \mathbf{0}\;\middle\vert\; X_t = j\right) + f_{ij}^{(j)}(T)\,, &\text{ if }& n_1 = 1 \\
\sum_{u=1}^{T-n_j+1} f_{ij}^{(j)}(u)P\left(N_j^{T-u} = n_j-1, \,N_{1:(j-1)}^{t,T-u} = \mathbf{0}\;\middle\vert\; X_t = j\right)\,, &\text{ if } & n_1\geq 2\\ \end{array}\right.$$
Where is the problem ?
My idea was to use the fact that $$P(N_{j}^T=n_{j}\mid N_{1:(j-1)}=n_{1:(j-1)}) = P(N_{j}^{T-\sum_{l=1}^{j-1}n_l}=n_{j}\mid N_{1:(j-1)}=\mathbf{0})\qquad\qquad (*)$$ and then by using the recursive relationship I established, I can compute the probability $P(N_{1:(K-1)}^T=n_{1:(K-1)})$. But, I finally realized that the relation $(*)$ is not true. If $(*)$ was true, the computation is quiet fast and take less time than computing $T^K$ operations.
My question
Can you help me finding a kind of relationship such as $(*)$ which is true and can help me achieve my goal in a similar logic as I did?
 A: I use a slightly different notation $N_t(j)$ for the number of visits of state $j$ during the time interval $[1,t]$.
A simpler approach is to note that the process $((X_t,N_t(1),\ldots,N_t(K)))_{t \ge 1}$ is a Markov chain on $\{1,\ldots,K\} \times \mathbb{Z}_+^K$ whose transitions are given by
$$q\big((i,n_1,\ldots,n_K),(j,n_1,\ldots,n_{j-1},n_j+1,n_{j+1},\ldots,n_K)\big) = p(i,j).$$
For all complex numbers $z_1,\ldots,z_n$, call $R_t(z_1,\ldots,z_K) \in \mathbb{C}^K$ the row matrix with components
$$f_{t,i}(z_1,\ldots,z_K) := \mathbb{E}[\mathbb{1}_{[X_t=i]} z_1^{N_t^1} \cdots z_K^{N_t^1}].$$
For every states $i$ and $j$,
$$f_{t+1,j}(z_1,\ldots,z_K) = \sum_{i=1}^K f_{t,i}(z_1,\ldots,z_K)p_{i,j}z_j.$$
Hence
$$R_{t+1}(z_1,\ldots,z_K) = R_{t}(z_1,\ldots,z_K) \times M \text{ with } M :=\big(p_{i,j}z_j \big)_{1 \le i,j \le K}.$$
By recursion,
$$R_T(z_1,\ldots,z_K) = (\pi_1z_1 \ldots \pi_Kz_K) \times M^{k-1}.$$
This a quite simple formula. Yet, in practice, computing the powers of $M$ is often difficult.
It is not very surprising. For example, for random walks on $\mathbb{Z}$ and for the Brownian motion in $\mathbb{R}$, a Ray-Knight theorem (which describes the joint law of the final position and local times at some random time $T$) are difficult to state for constant times.
