Expected number of summations of positive iid discrete random variables Let $X_i$ be a sequence of discrete i.i.d. random variables taking values in $\{ 1, \dots, K \}$, $K \in \mathbb{N}$ with given distribution $p_j = P(X_1 = j)$. Set $S_n := \sum_{i=1}^n X_i$ and for $M \geq 0$ consider the number of summation steps $N := \inf \{ n \mid S_n \geq M \}$ until the sum reaches $[M, \infty)$. Is there some closed-form expression for $\mathbb{E}[N]$ in terms of the given data ($p_j$, $M$)?
I know that one can formalize this problem as a discrete phase-type distribution  comprising a discrete-time Markov chain on the finite state space $\{ 0, 1, \dots, M+K-1 \}$ with initial state $0$, absorbing states $M, \dots, M+K-1$ and transitions $k \to k+j$ with probability $p_j$ for $k = 0, \dots, M-1$. Then $N$ is the (discrete) time until absorption and there is an expression for $\mathbb{E}[N]$ in matrix notation as $\delta_0 (I - T)^{-1} \mathbf{1}$ where $T$ is the transition probability matrix for the transient part. My hope is that one can simplify this expression directly in terms of ($p_j$, $M$) by e.g. some matrix-geometric methods since $T$ has a simple upper triangular form:
$$ T = \begin{pmatrix} 0 & p_1 & \dots & p_K & 0 & 0 & \dots & 0 \\ & 0 & p_1 & \dots & p_K & 0 & \dots & 0 \\ & & \ddots & \ddots & & \ddots & \\ 0 & 0 & \dots & & & & \dots & 0 \end{pmatrix} \in \mathbb{R}^{(M-1) \times (M-1)} $$
(On each next row, the row vector $(p_1, \dots, p_K)$ is shifted to the right-hand side until the last row is completely $0$ (the only transitions from state $M-1$ are to the absorbing states $M, \dots, M+K-1$).)
 A: Let $E_M$ denote the expected number of terms it takes to reach $M$. We have the equation
$$
E_M=\begin{cases}
0 &  M\le 0
\\
1+p_1E_{M-1}+p_2E_{M_2}+\dots+p_KE_{M-k} &M >0
\end{cases}
$$
Let $e(x)=\sum_{n\ge 1} E_nx^n$ be the generation function for $(E_1,E_2,\dots,)$. If you multiply both side of the the above equation by $x^M$, and take the sum over $M\in \{1,2,\dots\}$, the result is
$$
e(x)=\frac x{1-x}+(p_1x+p_2x^2+\dots+p_Kx^K)e(x)
$$
This allows you to solve for $e(x)$, and you can conclude that
$$
E_M=[x^M]\frac x{(1-x)(1-p_1x-p_2x^2-\dots -p_Kx^K)}
$$
Here, $[x^n]f(x)$ means "the coefficient of $x^n$ in the power series $f(x)$.
To calculate this, you would need to factorize the denominator, split it up with partial fractions to get a bunch of terms like$$\frac{C_i}{(1-r_ix)^{p_i}},$$ and note that each term contributes $C_i \binom{p_i+M-1}{p_i-1}r_i^M$ to the coefficient $E_M$.

If you need an exact answer, above is pretty much the best you can do. In what follows, I give a way to get a good approximation with exponentially decreasing error. As long as $M$ is moderately large, this will be very close.
Using Wald's equation, we have
$$
E[X_N]=E[N]\cdot E[X_1]
$$
where $E[X_1]=p_1+2p_2+\dots+Kp_K$. To compute $E[N]$, we just need $X_N$. We can write
$$
X_N=M+\text{overshoot}
$$
so all that remains is to find the distribution of the overshoot in terms of $(p_1,\dots,p_K)$. It turns out this distribution approaches a limiting distribution as $M\to\infty$ (exponentially quickly), and I showed how to compute that distribution in this answer. Rather, that answer covers the case where $K=6$ and $p_1=\dots=p_6=1/6$, but the same method generalizes easily.
