After seeing that some positivity problems get their solutions on MO, I am quite enthusiastic of posing my (and not only) problem of positive flavour.

In order to state it, I have to introduce the standard $q$-notation:
$$ (x ; q)_n = \prod _ {j = 1} ^n (1-q^{j-1} x) $$
which is seen to be a polynomial in both $x$ and $q$ for any *finite* $n=0,1,2,\dots$
(the empty product for $n=0$ has to be interpreted as $1$) but is also
meaningful for $n=\infty$ if $|q|<1$.

"My" polynomials are $$ P _n(q): =\sum _{r=1}^n \frac{ q^{r^2} (q;q) _{3n-r} (q^3;q^3) _{r-1} }{ (q;q) _r (q;q) _{2r-1} (q^3;q^3) _{n-r} }. $$ To convince you that they are indeed polynomials I write them as $$ P _n(q) =\sum _{r=1}^n(1+q^r)q^{r^2}\left[\begin{matrix} 3n-r \cr 2r\end{matrix}\right] _q \frac{(q^3;q^3) _{r-1}}{(q;q) _{r-1}}(q;q^3) _{n-r}(q^2;q^3) _{n-r} $$ where $$ \left[\begin{matrix} a \cr b\end{matrix}\right] _q =\frac{(q;q) _a}{(q;q) _b(q;q) _{a-b}} $$ are the $q$-binomial coefficients, or the Gaussian polynomials.

What I can show for the polynomials $P _n(q)$ is that their degree is $3n^2-1$, they are reciprocal (that is, $P _n(q)=q^{3n^2}P _n(1/q)$) and involve only powers of $q$ not divisible by 3.

What I cannot show is that the coefficients of these polynomials are nonnegative. Note that, for $r=1,\dots,n$, the polynomials $$ (1+q^r)q^{r^2}\left[\begin{matrix} 3n-r \cr 2r\end{matrix}\right] _q \frac{(q^3;q^3) _{r-1}}{(q;q) _{r-1}} $$ have nonnegative coefficients (as the Gaussian polynomials do) but the additional factor $$ (q;q^3) _{n-r}(q^2;q^3) _{n-r} =\prod _{j=1}^{n-r}(1-q^{3j-2})(1-q^{3j-1}) $$ changes the picture drastically.

An additional note in favour of the expected positivity is the limiting case $n\to\infty$: $$ P _\infty(q) =\frac{(q;q) _\infty}{(q^3;q^3) _\infty} \sum _{r=1}^\infty\frac{q^{r^2}(q^3;q^3) _{r-1}}{(q;q) _r(q;q) _{2r-1}}. $$ This series is the sum of two Virasoro characters and is thus a positive series. The latter can be also shown hypergeometrically and/or using the theory of modular forms (as $P _\infty(q)$ is a modular form).

So, why do the polynomials $P _n(q)$ have nonnegative coefficients? Any thoughts/links are greatly appreciated and acknowledged in advance.