The issue with this approach is that the symmetric direct product doesn't commute, but it is a bit more subtle than that. Summarizing in somewhat less mathematical terms the answer given to my question on MathSE, the issue stems from the fact that once you move past doing a single direct product, there is no longer just a symmetric and an antisymmetric part, but rather $n+1$ parts, where $n$ is the number of direct products being taken. So, trying to drop any antisymmetric terms that arise at each step will still leave terms that aren't in the symmetric part of the representation. (To say this in mathematical language for those interested, the representation for $n$ direct products of a particular irreducible representation of a group can be decomposed into $n+1$ parts by projecting onto the symmetric group $S_{n+1}$).
Physically, the reason for needing only the symmetric part rather than the whole direct product is that we are filling quanta into one particular mode. For example, if we wanted to find a combination band composed of two different $\Pi_u$ modes, we could just take the direct product to find the representation. But, filling two into a single $\Pi_u$ mode limits the number of possible distributions of the quanta to only those given by the symmetric direct product.
As to why the recursive formula gives the symmetric direct product, I give a summary of the argument from Wilson, Decius, and Cross's Molecular Vibrations for the case of a doubly degenerate mode.
Let $Q_a$ and $Q_b$ be the normal coordinates of a doubly degenerate mode. We choose these coordinates such that for an operation $R$ $$Q_a \xrightarrow{R}R_aQ_a$$
$$Q_b \xrightarrow{R}R_bQ_b$$
While the same coordinates will not diagonalize each $R$ in the group, we can find a set of coordinates for each $R$ that will. The character of $R$ for the $(v-1)^{\mathrm{th}}$ overtone, $\chi_v(R)$ can then be determined by summing the contribution of $R$ acting on each possible distribution of quanta between the modes $$\chi_v(R)=R_a^v+R_a^{v-1}R_b+...+R_aR_b^{v-1}+R_b^v$$
We can see that recursion formula holds by looking at the form of the terms on the right hand side
$$\chi(R)\chi_{v-1}(R)=(R_a+R_b)(R_a^{v-1}+R_a^{v-2}R_b+...+R_aR_b^{v-2}+R_b^{v-1})=R_a^v+2R_a^{v-1}R_b+...+2R_aR_b^{v-1}+R_b^v$$
and
$$\chi(R^v)=R_a^v+R_b^v$$
Averaging these two terms, we obtain the recursion formula from the question.
$$\chi_v(\hat{R})=\frac{1}{2}[\chi({R})\chi_{v-1}({R})+\chi({R}^v)]$$
By a similar argument, we can obtain a recursive expression for triply degenerate states
$$\chi_v({R})=\frac{1}{3}\bigg[2\chi({R})\chi_{v-1}({R})+\chi({R}^v)+\frac{1}{2}\big(\chi({R^2})+[\chi({R}]^2\big)\chi_{v-2}({R})+\chi({R}^v)\bigg]$$
