First, we repeat the arguments from this stackexchange answer. $P^{-1}$ is an $M$-matrix, and can thus be written as $s(I-A)$ for some positive $s$ and some $A$ with non-negative entries. As $P^{-1}$ is positive definite, the spectrum of $A$ lies to the left of $\{ z: \hbox{Re}(z) = 1 \}$, and hence by Perron-Frobenius the spectral radius of $A$ is less than $1$. Thus we have the absolutely convergent Neumann series $$ P = s^{-1} (I + A + A^2 + \dots )$$ and hence $$ PVPVP = s^{-3} \sum_{i=0}^\infty \sum_{j=0}^\infty \sum_{k=0}^\infty A^i V A^j V A^k.$$ It thus suffices to show that $$ \sum_{i+j+k=m} A^i V A^j V A^k \quad (1)$$ has non-negative coefficients for each $m \geq 0$ (where $i,j,k$ are understood to be non-negative integers). By change of variables, this is $$ \sum_{0 \leq q \leq r \leq m} A^q V A^{r-q} V A^{m-r}.$$ Writing $A = (a_{st})_{1 \leq s,t \leq n}$ and $V = \hbox{diag}(v_1,\dots,v_n)$, the $st$ coefficient of (1) can be expanded as $$ \sum_{s=s_0,s_1,\dots,s_m=t} a_{s_0 s_1} \dots a_{s_{m-1} s_m}\sum_{0 \leq q \leq r \leq m} v_{s_q} v_{s_r}.$$ But the quadratic form $$ \sum_{0 \leq q \leq r \leq m} x_q x_r = \frac{1}{2}(x_0+\dots+x_m)^2 + \frac{1}{2} x_0^2 + \dots + \frac{1}{2} x_m^2$$ is positive definite, and the $a_{st}$ are non-negative, and the claim follows.
[For the record, I found this argument while performing a perturbative analysis in the case where $P$ was close to $I$, or more precisely $P = (I-A)^{-1}$ for some $A$ with small non-negative entries.]
