Calculating the expected number of visits to a state by a DTMC Suppose we have a DTMC $X$ : $\{X_n : n = 0, 1,2,\dots\}$, a transition probability matrix $P$, and state space $S = \{1,2,3\}$. Suppose I want to calculate the expected amount of times we visit state $3$ in $T$ periods given that we start our chain at $X_0=1$. Let's also assume we have some function $C$ that maps values of $X$ to some "cost". We will define $C$ as such: $C(1) = 0, C(2) = 0, C(3) = 1$.
Then to calculate expected amount of times we reach state $3$, using this cost funtion we have:
$$E[\sum\limits_{i=1}^{T}C(X_i)|X_0=1]$$
in which case we would get an answer of something like:
$$[P+P^2+\dots+P^T] \begin{bmatrix} C(1) \\ C(2) \\ C(3) \end{bmatrix}$$
If this is not the case, what flaws are in my logic and what would a proper setup look like?
 A: For a target state $s$ and any state $a,$ let $v_a(s,T)$ be the expected number of visits to state $s$ (not counting the current state) upon making $T\ge 0$ transitions from state $a.$  From basic properties of expectation and the Markov property, we know

*

*$v_a(s,0)=0$ for all $a.$


*$v_a(s,T) = \sum_{r\in S} P_{ar} v_r(s,T-1)+ \mathcal{I}(r=s) .$
This is conveniently expressed in matrix notation by letting $\mathbf{v}(s,T)$ be the vector $(v_a(s,0))_{a\in S}:$

*

*$\mathbf{v}(s,0) = \mathbf{0}.$


*$\mathbf{v}(s,T) = \mathbb{P} \mathbf{v}(s,T-1) + \mathbf{e}_s.$
$\mathbf{e}_s = (0,\ldots, 0,1,0,\ldots,0)^\prime$ has its unique $1$ at position $s.$
By inspection, the solution (easily verified) is
$$\mathbf{v}(s,T) = \mathbb{P}\left(\mathbb{P}\left(\cdots \mathbf{e}_s\right) + \cdots\right) + \mathbf{e}_s = \left(\mathbb{P}^T + \mathbb{P}^{T-1} + \cdots + \mathbb{P} + \mathbb{I}_3\right)\mathbf{e}_s.$$
In the question with $s=3$ and $a=1$ you would take the first entry in this vector, $v_1(3,T).$

To illustrate, here is brute-force R code to compute $\mathbf{v}(s,n):$
v <- e <- rep(0,dim(P)[1])
e[s] <- 1
for (i in seq_len(n)) v <- P %*% v + e

(For long sequences of transitions you would want to diagonalize $\mathbb{P}$ and sum the resulting geometric series appearing the diagonal--but that's standard Markov Chain machinery which we needn't discuss here.)
Here is R code to simulate 10,000 transitions starting at state a:
sim <- replicate(1e4, {
  x <- a
  count <- 0
  for (i in seq_len(n)) {
    if (x == s) count <- count + 1
    x <- sample.int(dim(P)[1], 1, prob=P[x,])
  }
  count
})
mean(sim)

Run these calculations with your favorite transition matrix to compare the results.  If you don't have a favorite, here is code to generate interesting ones (they have a bunch of zeros in them):
d <- 5 # Number of states
P <- cbind(rexp(d), matrix(ifelse(runif(d*(d-1)) <= 2/d, 0, rexp(d^2)), d))
P <- P / rowSums(P)
a <- 1 # Start state
s <- d # Target
n <- 8 # Transitions

For instance, after setting the seed with set.seed(17), the output was
> c(mean(sim), v[a])
[1] 2.132200 2.137539

showing good agreement between the observed frequency and computed value.
