Moments of generalized fractional polynomial processes

We derive a moment formula for generalized fractional polynomial processes, i.e., for polynomial-preserving Markov processes time-changed by an inverse L\'evy-subordinator. If the time change is inverse $\alpha$-stable, the time-derivative of the Kolmogorov backward equation is replaced by a Caputo fractional derivative of order $\alpha$, and we demonstrate that moments of such processes are computable, in a closed form, using matrix Mittag-Leffler functions. The same holds true for cross-moments in equilibrium, generalizing results of Leonenko, Meerschaert and Sikorskii from the one-dimensional diffusive case of second-order moments to the multivariate, jump-diffusive case of moments of arbitrary order. We show that also in this more general setting, fractional polynomial processes exhibit long-range dependence, with correlations decaying as a power law with exponent $\alpha$.