Suppose I have a (left) stochastic matrix $P$, i.e. a non-negative matrix with column sums equal to 1. Its maximum eigenvalue will be equal to 1, and the corresponding eigenvector $\mathbf q$, if it's unique, can be interpreted the stationary distribution of a stochastic process with transition probabilities $P$.
Now suppose I have a family of stochastic matrices, $P(x)$, characterised by a real parameter $x$. I would like to know how the leading eigenvector $\mathbf{q}$ changes in response to a small change in $x$. In other words, given the values of $p_{ij}$ and $\frac{dp_{ij}}{dx}$, I would like to know $\frac{dq_i}{dx}$.
I'm aware that there is a general method for calculating the derivatives of eigenvectors. However, it's a little bit complicated, and I get the feeling it should simplify in this case, since we know that the leading eigenvalue must remain equal to 1 as the matrix changes.
Since $\mathbf{q}(x)$ is to be interpreted as a probability distribution, it's advantageous to assume it's normalised as $\sum_i q_i=1$, rather than $\mathbf{q}^T\mathbf{q}=1$.
I'm happy to accept answers that assume the leading eigenvector is unique (i.e. $P$ is primitive), but if there's an elegant way to handle the other cases, that would be great.
Assume that $P(x)$ has only one eigenvalue equal to $1$. Let $q(x)$ be the unique unit-norm leading eigenvector of $P(x)$ for all $x$. For convenience, I'll drop the "$(x)$" in the computations.
We require $Pq = q$, i.e. $(I-P)q = 0$ and $q^Tq = 0$ for all $x$.
Differentiation yields: $(I-P)\dfrac{dq}{dx} - \dfrac{dP}{dx}q = 0$ and $q^T\dfrac{dq}{dx} + \dfrac{dq}{dx}^Tq = 0$,
These can be rewritten as (1) $(I-P)\dfrac{dq}{dx} = \dfrac{dP}{dx}q$ and (2) $q^T\dfrac{dq}{dx} = 0$.
By assumption, $P$ has only one eigenvalue equal to $1$, so $I-P$ has only one $0$ eigenvalue. Thus, $I-P$ has a one dimensional nullspace, namely $\text{span}(q)$. So, if we can find one solution $\dfrac{dq}{dx} = v$ to (1), then all solutions to (1) will be of the form $\dfrac{dq}{dx} = v+tq$ for some $t \in \mathbb{R}$. Using the pseudoinverse, one solution is $\dfrac{dq}{dx} = (I-P)^+\dfrac{dP}{dx}q$. Therefore, the solutions to (1) are all in the form $\dfrac{dq}{dx} = (I-P)^+\dfrac{dP}{dx}q + tq$ for some $t \in \mathbb{R}$.
Condition (2) requires that $\dfrac{dq}{dx}$ is orthogonal to $q$. Since the solution $(I-P)^+\dfrac{dP}{dx}q$ is already orthogonal to the nullspace of $I-P$, which includes $q$, this is the only solution of (1) which also satisfies (2).
Therefore, the solution to (1) and (2) is: $\dfrac{dq}{dx} = (I-P)^+\dfrac{dP}{dx}q$.
If we instead wish to normalize $q$ such that $\sum_{i}q_i = 1$ instead of $q^Tq = 1$, then condition (2) becomes $\displaystyle\sum_{i}\dfrac{dq_i}{dx} = 0$. The solution will still be in the form $\dfrac{dq}{dx} = (I-P)^+\dfrac{dP}{dx}q + tq$ for some $t \in \mathbb{R}$. However, finding the value of $t$ to satisfy (2) might be harder.