Let $f$ be a formal (or convergent near the origin, as you wish) power series with variables $x_1$, ... , $x_n$ and complex coefficients. $$ f = \sum_I f_I \mathbf x^I $$
This power series is said to be differentially finite (or D-finite) if the $\mathbb C(\mathbf x)$-vector-space generated by all the derivatives $\partial^I f / \partial \mathbf x^I$ is finite dimensional.
Assume that $f_I$ is a product of integer power of factorial of integer valued linear forms of I (!) That is to say, there exists $\lambda_1$, ... , $\lambda_m$ linear forms over $\mathbb Q^n$ with integer coefficients, and $p_1$, ..., $p_m$, integers (possibly negative) such that $$ f_I = \prod_{k=1}^m \lambda_k(I)!^{p_k}.$$ In particular, if $I' = I + (0, \dotsc, 0,1,0,\dotsc )$, then $f_{I'}/f_I$ is a rational function of $I$.
With such an hypothesis, can we conclude that the power series $f$ is D-finite ?
Of course, we can translate the fact that $f_{I'}/f_I$ is an explicit rational function of $I$ into a first order recurrence with polynomial coefficients, and this recurrence can be translated into a differential equation for $f$. But I have no proofs that these differential equations are enough to get the D-finiteness of $f$. Experimentally, a Gröbner basis computation shows that they are enough.
It might be related to GKZ hypergeometric systems, as presented in the book Gröbner deformation of hypergeometric differential equations, but so far I haven't been able to apply this theory. Any reference, idea, or solution is welcome !