Finite sum of beta prime iid random variables

Question

The beta prime distribution is infinitely divisible, as proved in Steutel and van Harn, 2003 (Appendix B). Sadly, in this book, there is no expression of the parameters of the distribution of $n$ variables iid following a beta prime distribution, as it is done for some other distributions (normal, gamma, etc).

I cannot find how to derive it by myself. As an example, I am trying to derive the probability of the sum of two variables, and I am stuck in solving the basic convolution integral:

$$ f_{X+Y}(z)=\int_{-\infty}^{+\infty}f_X(x)f_Y(z-x)dx=\int_0^{z}\frac{1}{B(a,b)^2}\frac{(x(z-x))^{a-1}}{((1+x)(1+z-x))^{a+b}}dx $$

This distribution must be beta prime and I want to compute the parameters of this joint distribution. Is there a mean to get it from the convolution of PDF? Is there another mean to get directly the parameters of the beta prime pdf of the sum of n random variables iid following a beta prime distribution?

Answer 1

Let $\{X_i\}_{i=1}^n$ with $X_i\overset{\text{i.i.d}}{\sim}\beta^\prime(\alpha,\beta)$ and $Z=X_1+\dots+X_n$. It follows from linearity of the expected value that $$ \mathsf EZ=\sum_{i=1}^n\mathsf EX_i=n\mathsf EX=\frac{n\alpha}{\beta-1},\quad\beta>1. $$ Furthermore, by mutual independence of the $X_i$'s we have $$ \mathsf{Var}Z=\sum_{i=1}^n\mathsf{Var}X_i=n\mathsf{Var}X=\frac{n\alpha(\alpha+\beta-1)}{(\beta-2)(\beta-1)^2},\quad\beta>2. $$ Since we know $Z\sim\beta^\prime(\gamma,\delta)$ we may write the system of equations $$ \begin{aligned} \frac{n\alpha}{\beta-1} &=\frac{\gamma}{\delta-1}\\ \frac{n\alpha(\alpha+\beta-1)}{(\beta-2)(\beta-1)^2} &=\frac{\gamma(\gamma+\delta-1)}{(\delta-2)(\delta-1)^2}. \end{aligned} $$ Solving this system for $\gamma$ and $\delta$ subsequently yields the following result: $$ \begin{aligned} \gamma &=\frac{\alpha n \left(\alpha +\beta ^2-2 \beta +\alpha \beta n-2 \alpha n+1\right)}{(\beta -1) (\alpha +\beta -1)}\\ \delta &=\frac{2 \alpha +\beta ^2-\beta +\alpha \beta n-2 \alpha n}{\alpha +\beta -1}. \end{aligned} $$ With some algebra you may be able to simplify these expressions.

Update:

Based on the discussion surrounding the exactness/correctness of the results I decided to perform an experiment in MATLAB. Here is the code used which performs the simulation for $\alpha=\beta=15$ and $n=5$:

a = 15; %alpha
b = 15; %beta
n = 5;
c = (a*n*(a-2*b-2*a*n+b^2+a*b*n+1))/((b-1)*(a+b-1)); %gamma
d = (2*a-b-2*a*n+b^2+a*b*n)/(a+b-1); %delta

Xdata = 1./betarnd(a,b,1e6,n)-1;
Xn = sum(Xdata,2); %Xn = X_1+X_2+...+X_n
ax = linspace(0,max(Xn),256);
f_Xn = @(x) x.^(c-1).*(1+x).^(-c-d)/beta(c,d);

figure
hold on
histogram(Xn,64,'normalization','pdf')
plot(ax,f_Xn(ax),'Color',[0,0,0],'LineWidth',1.5)
xlabel(['X_' num2str(n)])
ylabel('density')
box on
hold off

Here we see that the histogram does indeed agree with the theoretical beta prime distribution with parameters $\gamma$ and $\delta$ as derived above. I tried other values of $\alpha$, $\beta$ and $n$ with similar results.