Let $X\sim Poi(5)$ and $N\sim Bin(60,1/2)$. We have portfolio given by: $$S_N=X_1+\dots+X_N$$ I already did some calculation in Maple and ended up with this density plot (on the graph the pink line is the approximation of $S_N$ by normal distribution):
Now i want to approximate it by distribution $Gamma(x_0,\alpha,\beta)$ Where: $$\alpha=4\frac{(VarS_i)^3}{(\mathbb{E}[(S_i-\mathbb{E}S_i)^3])^{2}},$$ $$\beta=2\frac{VarS_i}{\mathbb{E}[(S_i-\mathbb{E}S_i)^3]},$$ $$x_0=\mathbb{E}S_i-\frac{(VarS_i)^2}{\mathbb{E}[(S_i-\mathbb{E}S_i)^3]}$$
These parameters I also calculated and it shows that I need $Gamma(-66,356,0.82)$. Now here's my maple code:
A:=RandomVariable(GammaDistribution(365,.82))
A := _R1
CDFA := t -> CDF(A,t)
CDFA := proc (t) options operator, arrow; Statistics:-CDF(A, t) end proc
GammaXA := t -> CDF(A,t-66)
GammaXA := proc (t) options operator, arrow; Statistics:-CDF(A, t-66) end proc:
GammaXA(t)
int(-0.6935310692e-2*(-1+Heaviside(-_t))*exp(-(1/365)*_t)/_t^.18, _t = -infinity .. t-66)
GammaXdist := Statistics[Distribution](CDF =(t -> GammaXA(t)))
GammaXdist := _m194201600
Gamma11 := RandomVariable(GammaXdist)
Gamma11 := _R2
PDF(Gamma11,t)
-0.6935310692e-2*(-1+Heaviside(-t+66))*exp(-(1/365)*t+66/365)/(t-66)^.18
dense := DensityPlot(Gamma11):
display({plotN,dense})
The plot I get is this one:
What did I do wrong here? I must admit I'm a total newbie when it comes to Maple and I don't know if there's a smarter way to do this.
The hierarchical model is such that $$\operatorname{E}[S] = \operatorname{E}[\operatorname{E}[S \mid N]] = \operatorname{E}[5N] = 5(60)(1/2) = 150,$$ and $$\operatorname{Var}[S] = \operatorname{E}[\operatorname{Var}[S \mid N]] + \operatorname{Var}[\operatorname{E}[S \mid N]] = \operatorname{E}[N \operatorname{Var}[X]] + \operatorname{Var}[N \operatorname{E}[X]] \\ = \operatorname{E}[5N] + \operatorname{Var}[5N] = 5(60)(1/2) + 5^2 (60)(1/2)(1/2) = 525.$$
Then a gamma distribution with shape-rate parametrization $$f_Y(y) = \frac{\beta^\alpha y^{\alpha-1} e^{-\beta y}}{\Gamma(\alpha)}$$ has $$\operatorname{E}[Y] = \alpha/\beta, \quad \operatorname{Var}[Y] = \alpha/\beta^2,$$ hence matching these to the calculated mean and variance of the Poisson-binomial model, we get $$\alpha = \frac{300}{7}, \quad \beta = \frac{2}{7}.$$ If the parametrization is shape-scale, then $\theta = 1/\beta = 7/2$ where the density is $$f_Y(y) = \frac{y^{\alpha-1} e^{-y/\theta}}{\theta^\alpha \Gamma(\alpha)}.$$ In the case of Maple, the documentation would suggest that it is using a shape-scale parametrization. Mathematica also uses shape-scale. I have produced the following plot in the latter:
Here, the blue curve is the exact distribution of $S$, and the orange is the gamma approximation. As you can see, the fit is not very good, certainly worse than for the normal approximation. We can see it more clearly by computing the difference:
Here, the blue curve is the difference between the gamma distribution and $S$, and the orange curve is the difference between the normal distribution and $S$. It would also seem to suggest that a weighted average of the gamma and normal densities would perform much better than either of them alone; e.g., $$f_S(s) \approx \frac{1}{3}f_Y(s) + \frac{2}{3}f_Z(s)$$ where $Y$ is gamma and $Z$ is normal. But mixture densities are notoriously intractable for many purposes.