Under the assumption that $f(x)$ is a monic polynomial, here's a sketch of how one could go about determining the coefficients of this expansion. It should be evident that the expression whose expansion is requested can be written as a polynomial of degree $D=p^2-p$. We write
$$g(x;\vec{\alpha})=P(x;\alpha_1)...P(x;\alpha_p)=\sum_{n=0}^D\Gamma_n(\vec{\alpha})x^n$$
where $P(x;r)=\sum_{n=0}^{p-1}r^{p-1-n}x^n$
It is easy to see by the definition of $g$ that it is invariant under any permutation $S=\{\sigma{(1)},...,\sigma{(p)}\}$ of the set $\{\alpha_1,..., \alpha_p\}$. We readily conclude that the coefficients have the same symmetry:
$$\Gamma_n(S\vec{\alpha})=\Gamma_n(\vec{\alpha})$$
Also the coefficients can be shown to obey the very important scaling relation
$$\Gamma_n(\mu \vec{\alpha})=\mu^{D-n}\Gamma_n(\vec{\alpha})$$
Now, also by the definition of $g$ and simple power counting one should be able to show that the coefficients are polynomials in multiple variables made up only of positive powers of the variables $\alpha_1,...,\alpha_p$. This along with the scaling relation and the symmetry requirement suffices to show that the coefficients must be symmetric polynomials of total degree $D-n$.
With this proven, one can now invoke the fundamental theorem of multivariable symmetric polynomials, which asserts that no matter the form of a symmetric polynomial in many variables, it always has an expression in terms of elementary symmetric polynomials, and since those can be computed to be integer valued by the assumption that the coefficients of the polynomial with roots $\alpha_1,...,\alpha_p$ are integers. Therefore, we conclude that $\Gamma_n \in \mathbb{Z}$.
It takes quite a bit more work to find explicit expressions for the coefficients defined above. Define a basis of symmetric polynomials of total degree $N$ as follows
$$s_{N}^{q_1,..., q_p}(\alpha_1, \alpha_2,... \alpha_p)=\sum_{\text{sym}}\alpha_1^{q_1}... \alpha_p^{q_p}~~,~~ \sum_{i=1}^p q_i=N$$
I conjecture that the coefficients can be expressed as an equal weight linear combination of the polynomials in the basis defined above as follows:
$$\Gamma_n(\vec{\alpha})=\sum_{r\in R} s^{r}_{D-n}(\vec{\alpha})$$
where $r$ is a $p$-index in the set defined by
$$R=\{(r_1,..., r_p):\sum_{i=1}^p{r_i}=D-n,0\leq r_i\leq p-1 \}$$
These sets of p-indices can be generated in Mathematica by evoking the function
r[p_, n_] := Select[IntegerPartitions[p^2 - n, {p}] - Table[Table[1, {i, 1, p}], {j, 1, Length[IntegerPartitions[p^2 - n, {p}]]}], Max[#] <= p - 1 &]
Interestingly enough, this symmetric polynomial basis is also good enough to solve the more general problem with
$$P(x;r)=\sum_{n=0}^{p-1}c_{p-1-n} x^n r^{p-1-n}$$
which admits a solution which is a simple generalization of the previous one
$$\Gamma_n(\vec{\alpha})=\sum_{r\in R} c_{r_1}...c_{r_p}s^{r_1,...,r_p}_{D-n}(\vec{\alpha})$$
Note that the coefficients will still be integral if the $c$'s are integral.