Wondering about the question I posted on Friday (on StackExchange, no satisfactory answer so far!), I looked further at the special case of the gamma distribution I suggested at the end. Starting from the moment conditions,
the solution is (hopefully) given by the system
The resolution of this system obviously imposes conditions on those moments, like
So I ran a small R experiment checking when there was no acceptable solution to the system. I started with five moments that satisfied the basic Stieltjes and determinant conditions
# basically anything mu=runif(2,0,10) # Jensen inequality sig=c(mu^2/runif(1),mu^2/runif(1)) # my R code returning the solution if any sol(mu,c(sig,runif(1,-sqrt(prod(sig)),sqrt(prod(sig)))))
and got a fair share (20%) of rejections, e.g.
> sol(mu,c(sig,runif(1,-sqrt(prod(sig)),sqrt(prod(sig))))) $solub  FALSE $alpha  0.8086944 0.1220291 -0.1491023 $beta  0.1086459 0.5320866
However, not being sure about the constraints on the five moments I am now left with another question: what are the necessary and sufficient conditions on the five moments of a pair of positive vectors?! Or, more generally, what are the necessary and sufficient conditions on the k-dimensional μ and Σ for them to be first and second moments of a positive k-dimensional vector?