# Symmetric functions arising from continuous-time Markov chains

I'm interested in the following symmetric functions $s_k: \mathbb{R}_+^{k}\mapsto\mathbb{R_+}$: \begin{align*} s_k(x_1,x_2,\dots, x_k)&= \int\limits_{0=t_0<t_1<\dots<t_{k-1}<t_k=1} e^{-(t_1-t_0) x_1} e^{-(t_2-t_1)x_2}\dots e^{-(t_k-t_{k-1})x_k}dt_1 \dots dt_{k-1} \\ &=\sum_{i=1}^k \frac{e^{-x_i}} {\displaystyle\prod_{\substack{j=1,\dots,n\\ j\ne i}}(x_j-x_i)} \end{align*}

The formula on the second line holds whenever the $x_i$ are all distinct; however it becomes degenerate if $x_i=x_j$ for some $i\ne j$. Question: is there some "nice" pattern for the closed form expressions one obtains when some of the $x_i$ coincide? (Examples below.)

[ One motivation: consider a continuous-time Markov chain, with state space $\{1,2,3,\dots\}$, with rate $q_{ij}$ of jumping from state $i$ to state $j$ for $i\ne j$, and total rate $\lambda_i=\sum_{j\ne i} q_{ij}$ of leaving state $i$.

Suppose the chain starts in state $1$ at time $0$. Consider the event that during the time-interval $[0,1]$, the chain follows precisely the path $1\to2\to\dots\to (k-1)\to k$ and then stays in state $k$ until the end of the interval (making $k-1$ jumps in all).

This event has probability $q_{12}\dots q_{(k-1)k}s_k(\lambda_1, \dots, \lambda_k)$. ]

Anyway: everything is nice and smooth and closed-form expressions for cases with repeated entries can straightforwardly be derived (using, say, L'Hôpital or whatever) - however the details appear to get messy quite quickly. I'm wondering if there is some convenient or systematic way to write those expressions (other than just calculating each case).

For example $s_2(x, y)=(e^{-x}-e^{-y})/(y-x)$, and then $s_2(x, x)=e^{-x}$.

Taking e.g. $k=4$, we get (I think) \begin{align*} s_4(x,x,y,y)&=\frac{2(e^{-x}-e^{-y})}{(x-y)^3}+\frac{(e^{-x}+e^{-y})}{(x-y)^2} \\ s_4(x,x,x,y)&=\frac{e^{-y}-e^{-x}}{(x-y)^3}-\frac{e^{-y}}{(x-y)^2}+ \frac{e^{-y}}{2(x-y)} \\ s_4(x,y,z,z)&=\frac{e^{-x}}{(z-x)^2(y-x)} + \frac{e^{-y}}{(z-y)^2(x-y)} \\ &\,\,\,\,\,\,\,+e^{-z}\left\{\frac{(z-x)+(z-y)+(z-x)(z-y)+z^2}{(z-x)^2(z-y)^2}\right\}. \end{align*}

I'm aware that these things can be written in terms of multiple integrals with gamma densities; for example, if $k=k_1+\dots+k_m$ and there are $k_i$ arguments $x_i$ for each $i$, then I think \begin{equation*} s_k(\dots)=\int\limits_{0=t_0<t_1<\dots<t_{m-1}<t_m=1} \prod_{i=1}^m \left\{\frac{(t_i-t_{i-1})^{k_i-1}}{(k_i-1)!} e^{-(t_i-t_{i-1})x_i}\right\}dt_1\dots dt_{m-1}, \end{equation*} but I'm not sure that leads anywhere particularly nice.

Perhaps (beguiled by the beautiful form in the case of all $x_i$ distinct) I'm looking for something unreasonable, and these things just are what they are. Then again maybe someone has something nice to point out!

• Regarding your parenthetical remark: with one jump ($k=2$), the probability of the event you describe in square parenthesis is $q_{12}/\lambda_1$ times $P(t_1<1)$ times $P(t_2>1)$. Is that really equal to $q_{12} s_2(\lambda_1, \lambda_2)$? – Nawaf Bou-Rabee Oct 10 '16 at 19:57
• I don't think it's the product in the way you suggest, since the relevant events are not independent. But yes, following that idea you can write it as $\int_0^1 \frac{q_{12}}{\lambda_1}\lambda_1 e^{-\lambda_1 t} e^{-\lambda_2(1-t)}dt$ (here $t$ plays the role of the first jump time) which does seem to work out as claimed. – James Martin Oct 10 '16 at 20:43
• Yes, but what happened to the normalization factor $\lambda_2$? The first one cancels with the intensity of the first jump, but what about the second one? – Nawaf Bou-Rabee Oct 10 '16 at 20:50
• The probability that the holding time in state 2 is at least $(1-t)$ is $e^{-\lambda_2(1-t)}$; there is no factor of $\lambda_2$. – James Martin Oct 10 '16 at 20:53

Let the random variable $t_{k-1}$ be the sum of the holding times in states $1$ up to state $k-1$. Note that $t_{k-1}$ is a hypoexponential random variable. Let $f_{k-1}$ denote the PDF of this random variable. From the probabilistic interpretation given by the OP, the quantity of interest can be written as: \begin{align*} s_k &= \frac{1}{\lambda_1 \cdots \lambda_{k-1}} \int_0^1 f_{k-1}(z) \int_{1-z}^{\infty} \lambda_k e^{-\lambda_k s} ds dz \\ &=\frac{e^{-\lambda_k}}{\lambda_1 \cdots \lambda_{k-1}} \int_0^1 f_{k-1}(z) e^{z \lambda_k} dz \end{align*} To be sure, this basically comes from the probability that $t_{k-1}<1$ and that the holding time in state $k$ is at least $1-t_{k-1}$.

To avoid assuming that the jump rates are distinct, we will use the following matrix representation of the PDF $f_{k-1}$: $$f_{k-1}(z) = - \boldsymbol{\alpha}_{k-1} e^{z \Theta_{k-1}} \Theta_{k-1} \mathbf{1}_{k-1}$$ where I borrow the notation from here, except I have introduced subscripts to $\boldsymbol{\alpha}$, $\Theta$ and $\mathbf{1}$ in order to express the number of jump rates in the hypoexponential distribution. The matrix $\Theta_{k-1}$ is upper bidiagonal, and as long as all of the rates are positive, this matrix is invertible. Basically, we used the fact that a hypoexponential distribution is an example of a phase-type distribution.

Let $\mathbf{I}_{k-1}$ be the $(k-1) \times (k-1)$ identity matrix. By integrating by parts, it is straightforward to show that: $$\int_0^1 e^{z \Theta_{k-1}} e^{z \lambda_k} dz = (\Theta_{k-1} + \lambda_k \mathbf{I}_{k-1})^{-1} (e^{\Theta_{k-1}} e^{\lambda_k} - \mathbf{I}_{k-1} )$$ Hence, $$s_k = \frac{-e^{-\lambda_k}}{\lambda_1 \cdots \lambda_{k-1}} \boldsymbol{\alpha}_{k-1}(\Theta_{k-1} + \lambda_k \mathbf{I}_{k-1})^{-1} (e^{\Theta_{k-1}} e^{\lambda_k} - \mathbf{I}_{k-1} ) \Theta_{k-1} \mathbf{1}_{k-1}$$ Note that nowhere in this derivation did we assume that the jump rates are distinct. This formula is straightforward to implement numerically.

• Thanks for this answer. As far as numerical implementations go, one is rather spoilt for choice. For example, the matrix of time-$t$ transition probabilities is given by $P_t=\exp(tQ)$ where $Q$ is the generator of the Markov chain. As you point out it suffices to calculate this for a bi-diagonal matrix. However, what I'm specifically asking about (perhaps perversely) is the closed-form expression. – James Martin Oct 11 '16 at 10:31
• The formula given may be evaluated explicitly and does not require any time discretization. Granted though that it is not as simple as the formula you gave in the case of distinct rates. – Nawaf Bou-Rabee Oct 11 '16 at 12:18
• Also by closed form expression, I think you mean a combination of elementary functions. The formula given is one such combination, but written in terms of a matrices. If I had time, I would try to simplify that expression, but unfortunately I have other things to do right now. – Nawaf Bou-Rabee Oct 11 '16 at 12:45
• OK thanks! This is promising. Following your link to en.wikipedia.org/wiki/Hypoexponential_distribution#General_case seems to give something very relevant, which may be just what I need. – James Martin Oct 11 '16 at 13:03
• You're welcome. I am baffled why your question has not been upvoted much, and actually got downvoted (probably by accident). It's a great question. Thank you for asking it. – Nawaf Bou-Rabee Oct 11 '16 at 13:08

$\text{Dear James,}$ you can write your function as a(n infinite) sum of homogeneous symmetric functions \begin{equation*}%$h_\ell(x_1, \dots, x_k) := \sum_{1 \leqslant i_1 \leqslant \cdots \leqslant i_\ell \leqslant k} x_{i_1} \cdots x_{i_\ell} = [t^\ell]\prod_{j = 1}^k \frac{1}{1 - t x_j} \end{equation*} where$ [t^\ell]f(t) $is the$\ell$-th Fourier coefficient of$f$(for this last equality, see wikipedia or the book by Macdonald). Now, decomposing the rational fraction, we get \begin{equation*}%$ \prod_{j = 1}^k \frac{1}{1 - t x_j} = \sum_{j = 1}^k \frac{A_j(X)}{1 - t x_j}, \quad A_j(X) = \prod_{i \neq j} \frac{1}{1 - x_j^{-1} x_i} \end{equation*} thus, taking the coefficient of $[t^\ell]$, one gets \begin{equation*}%$h_\ell(x_1, \dots, x_k) = \sum_{j = 1}^k A_j(X) x_j^\ell = \sum_{j = 1}^k \frac{x_j^{\ell- k + 1} }{ \prod_{i \neq j} (x_j - x_i) } \end{equation*} Your function thus writes \begin{equation*}%$ s_k(x_1, \dots, x_k) = \sum_{\ell \geqslant 0} \frac{(-1)^\ell}{\ell !} h_{\ell + k - 1}(x_1, \dots, x_k) \end{equation*} and is hence projective due to the projectivity of the $h_i$'s. As a result, setting two variables equal amounts to get (for instance) $h_\ell(x_1, \dots, x_{k - 1}, x_{k - 1})$, which does not have such a simple expression. But we can continue the analysis, using the blog of Tao (https://terrytao.wordpress.com/2017/08/06/schur-convexity-and-positive-definiteness-of-the-even-degree-complete-homogeneous-symmetric-polynomials/). One can get a really nice representation of the complete homogeneous symmetric functions writing $\frac{1}{1 - t x} = \mathbb{E} \left( e^{ tx \mathbb{e} } \right)$ where $\mathbb{e} \sim \mathrm{Exp}(1)$. Considering a sequence of independent such exponential random variables $(\mathbb{e}_i)_{1 \leqslant i \leqslant k}$, one thus gets \begin{equation*}%$h_\ell(x_1, \dots, x_k) = [t^\ell]\prod_{j = 1}^k \frac{1}{1 - t x_j} = [t^\ell]\mathbb{E}\left( e^{ t \sum_{i = 1}^k x_i \mathbb{e}_i } \right) = \frac{1}{\ell !} \mathbb{E}\left( \left( \sum_{i = 1}^k x_i \mathbb{e}_i \right)^\ell \right) \end{equation*} Since$ x_i \in \mathbb{R}_+ $, your function admits a representation as the Bessel expectation \begin{equation*}%$ s_k(x_1, \dots, x_k) = \mathbb{E}\left( \sum_{\ell \geqslant 0} \frac{(-1)^\ell}{ \ell ! (\ell + k - 1)! } \left( \sum_{i = 1}^k x_i \mathbb{e}_i \right)^{\ell + k - 1} \right) = \mathbb{E}\left( J_{k + 1} \left(2 \sqrt{ \sum_{i = 1}^k x_i \mathbb{e}_i } \right)\right) \end{equation*}

You can transform again this form using the generating series of the Bessel functions and write $J_k(x)$ as the Fourier coefficient $[t^k] e^{ x (t - t^{-1})/2 }$ or as a Laplace transform. Last, you can even use a stable $1/2$ random variable to write the final result as $\mathbb{E}\left( \exp\left( Z \sum_{i = 1}^k x_i \mathbb{e}_i \right)\right)$ where $Z$ is a random variable independent of the $\mathbb{e}_i$'s.

Under this form, the operation of setting two variables equal amounts to change an exponential random variable by a sum of two exponential random variables.

Best, Y.