Indeed, these formulas are standard. Their derivation in a slighly more general setting than yours is as follows. 

For every $t\ge0$, let 
$$
M_t=\displaystyle\exp\left(\int_0^tv(X_s)\mathrm{d}s\right),
$$ 
for a given function $v$ defined on the state space of the process $(X_t)$. (In your setting, every $X_t$ is real valued and $v(x)=\mathrm{i}x$ for every $x$ but these details are irrelevant.) For every states $x$ and $y$, let 
$$A_t(y)=[X_t=y],\qquad
Q_t(x,y)=E(M_t1_{A_t(y)}\vert X_0=x).
$$
Let $Q_t$ denote the associated square matrix (indexed by the state space, possibly infinite). For instance, $Q_0$ is the identity matrix. Note also that in the expression of $Q_t(x,y)$, $[X_0=x]$ appears as a <b>conditioning</b> while $A_t(y)=[X_t=y]$ is the event to which the expectation is <b>restricted</b> and that these are different operations hence your interpretation of Kubo's method should be rephrased.

The dynamics of $(Q_t)$ is driven by a linear differential equation $Q'_t=GQ_t$, where $G$ is a deformation of the infinitesimal generator of the process $(X_t)$. To identify $G$, one can compute $Q_{t+s}$ at the order $s$, for $s > 0$, when $s$ is small. 

To do so, call $r(x,y)$ the transition rate of $(X_t)$ from $x$ to $y\ne x$, and $c(x)$ the sum over $y\ne x$ of $r(x,y)$. (In Kubo's setting as reproduced in your post, $c(x)$ is your $c_x$ and $r(x,y)$ is your $c_xp_{xy}$. By the way, the sum over $y\ne x$ of your $p_{xy}$ should be $1$ instead of $0$ and you should make up your mind between the notations $p_{xy}$ and $P_{xy}$.)

Then, conditioning on $[X_0=x]$, one can decompose the expectation which defines $Q_{t+s}(x,y)$ along the values of $X_s$. This decomposition goes as follows. For every $z\ne x$, $X_s=z$ with probability $r(x,z)s+o(s)$, and $X_s=x$ with probability $1-c(x)s+o(s)$.
Furthermore, for every $z\ne x$, $M_{t+s}=(1+o(1))M_t$ on $[X_0=x,X_s=z]$.
And on $[X_0=X_s=x]$, the probability of a double transition in the time interval $[0,s]$ is $o(s)$, hence $M_{t+s}=(1+v(x)s+o(s))M_{t+s}/M_s$ where $M_{t+s}/M_s$ is distributed like $M_t$ conditional on $[X_0=x]$.

All this leads to
$$
Q_{t+s}(x,y)=Q_t(x,y)(1+v(x)s)(1-c(x)s)+\sum_zQ_t(z,y)r(x,z)s+o(s).
$$
When $s\to0$, one gets
$$
Q'_t(x,y)=(v(x)-c(x))Q_t(x,y)+\sum_zr(x,z)Q_t(z,y).
$$
In other words, $G(x,x)=v(x)-c(x)$ for every $x$
and $G(x,y)=r(x,y)$ for every $y\ne x$. These are the equations in Kubo's paper.