Getting Parameter of Translated Gamma Distribution from Monte Carlo

Spin-off from here.

(Edit) Main question: What do I do about a parameter whose suggested values range quite vastly?

(Edit) Backstory: I am given data of loss values and the dates that correspond to when each loss was incurred. I am to fit a distribution for the aggregate losses: I must first simulate using poisson or negative binomial the frequency of loss (some positive integer, usually less than 15) and then simulate losses given the frequency (e.g. simulate 15 loss values) that follow some distribution e.g. loglogistic, mixture, lognormal. I have to them sum up those losses and that counts as the first aggregate loss. I have to do this 2000 times and then fit those 2000 values into a distribution.

Richard referred to me an article that tells me how to get parameters of a translated gamma distribution to which I should consider fitting simulated aggregated loss values.

The parameters depend on moments of S (or. in Richard's terms, L):

When I simulate the parameters, as we would expect I get different values each time. The values for $\alpha$ and $\beta$ don't vary much and seem to be close to zero, but $x_0$ seems to vary each time. Iirc, I got values ranging from -100,000 to -600,000. How should I know what $x_0$ to use? Do I get an average of 1000 x_0's?

This would seems impractical even if it were 100 x_0's since each x_0 is obtained from 2000 simulations (the requirement of the project).

Btw, I am assuming the $E(S^n)$'s can be approximated with mean($S^n$)'s. Is that right?

If you estimate $E[S^n]$ by averaging the $S_k,k=1,\ldots,K$ of the sample that you have generated by MC simply means that you apply the empirical distribution function in order to calculate this quantity. I.e. $$E[S^n] \approx \sum_{k=1}^K p_k S^n_k = \frac1K \sum_{k=1}^K S^n_k$$ and $p_k = 1/K$ for all $k$.