An equivalent problem is to show the positivity of the connection factors $c^1_{i,j}$ in the expansions
$$p_i(t)p_j(t) = \sum^{i+j}_{n=1}\; c^n_{i,j}p_n(t),$$
where $p_n(t)$ are cycle index partition polynomials of the symmetric groups (A036039) with the indeterminates $x_n = (-1)^{n-1}h_{n-1}t$ and $h_n$ are the complete homogeneous symmetric polynomials with all of their indeterminates positive. The $c^1$ are essentially the coefficients of the FGL expansion Strickland displays.
Jair, in your Sage computations, if the coefficients of f(x) are expressed as its Taylor series coefficients, i.e., they are normalized by the factorials, it is easier to recognize them as A145271, the refined Eulerian numbers. Same for the polynomials p, and if the coefficients of p are grouped together by powers of t and expressed as the elementary symmetric polynomials/functions $e_n=n!t_n=n!\phi_n \;$, it is easy to see they are signed A036039 with the appropriate determinates given above, e.g., $3! p_3 = 2(e_1^2-e_2)t - 3e_1t^2+t^3 = 2h_2t-3h_1t^2+t^3$.