problem formulation
Given a discrete probability density $\{p(i)\}_{i=0}^\infty$, let $G(z)$ be the following power series (called probability generating function) \begin{equation}G(z)\triangleq \sum_{i=0}^\infty p(i)\cdot z^i\end{equation} and abbreviate as $G^{(n)}(z)$ the $n$-th derivative of $G(z)$ with respect its argument $z\in\mathbb{R}$, \begin{equation}G^{(n)}(z)\triangleq \frac{\text{d}^{n}\,G(z)}{\text{d}z^n}\end{equation}
I have to resolve the following problem: compute the $n$-th derivative, where $n$ is an arbitrary integer, of the ratio $G^{(1)}(z)/G(z)$, i.e. I have to find the explicit expression of \begin{equation}\zeta_n(z)\triangleq \left(\frac{G^{(1)}(z)}{G(z)}\right)^{(n)}\end{equation} in function of the given density $p(i)$.
partial solution
a quick search (e.g. see here) tells me that a possible approach consist into consider $G^{(1)}/G$ as $G^{(1)}\cdot 1/G$ and then to apply the Leibnitz rule, so \begin{equation}\begin{aligned} \zeta_n(z)&=\sum_{k=0}^n\binom{n}{k}\cdot\left(G^{(1)}(z)\right)^{(k)}\cdot\left(\frac{1}{G(z)}\right)^{(n-k)}\\ &=\sum_{k=0}^n\binom{n}{k}\cdot G^{(1+k)}(z)\cdot\left(\frac{1}{G(z)}\right)^{(n-k)}\\ \end{aligned}\end{equation} here the first derivative $G^{(1+k)}(z)$ is straightforward because $G(z)$ is a power series, in fact for any integer $k$ \begin{equation}\begin{aligned} G^{(k)}(z)&\triangleq\frac{\text{d}^k}{\text{d}z^k}\left[\sum_{i=0}^\infty p(i)\cdot z^i\right]= \sum_{i=0}^\infty p(i)\cdot\frac{\text{d}^k}{\text{d}z^k}\left[z^i\right]= \sum_{i=0}^\infty p(i)\cdot \left(P_k^i\cdot z^{i-k}\right)\\ \end{aligned}\end{equation} where I have introduced the permutation coefficient \begin{equation}P_k^i\triangleq \begin{cases}\frac{i!}{(i-k)!} & \text{if } i \geq k \\ 0 & \text{otherwise}\end{cases}\end{equation} the derivative $G^{(k)}(z)$ is still a power series, indeed by exploiting the fact that $P_k^i=0$ when $i<k$, \begin{equation}\begin{aligned} G^{(k)}(z)&=\sum_{i=k}^\infty p(i)\cdot \left(P_k^i\cdot z^{i-k}\right)=\sum_{i=0}^\infty p(i+k)\cdot \left(P_k^{i+k}\cdot z^{i}\right)\\ &=\sum_{i=0}^\infty p(i+k)\cdot \left((i+k)!\cdot z^{i}\right)=\sum_{i=0}^\infty p_k(i)\cdot z^{i} \end{aligned}\end{equation} where I have introduced the non-negative sequence (I have a strong suspect, but I'm not sure, that in general it is not a discrete density anymore - however this is an off topic problem) \begin{equation}p_k(i)\triangleq (i+k)!\cdot p(i+k)\end{equation}
question
the real problem is to compute for any $k$ the derivative \begin{equation}\left(\frac{1}{G(z)}\right)^{(k)}\tag{1}\end{equation} this derivative can be computed by using the Faa' di Bruno's formula, which if I'm not wrong is a chain rule of arbitrary order, because $1/G(z)$ can be written as the composition \begin{equation}\frac{1}{x}\bigg|_{x=G(z)}\end{equation} however, before to dig into the rabbit hole of complex computations, I would like to have some external suggestion about how to compute $(1)$. My question is the following:
it is really necessary to use Faa' di Bruno or, due to the special form of $G(z)$ (which I know it is not a generic function but a nice power series), it is possible to use some clever method?
