Let $a\in (0,1), \;\;\psi_a(x):=\prod_{j=0}^\infty (1-a^{2j+1}x).$
Question. Is it true that, for all $x\in [0,1]$ and all $k\in\mathbb{N},$ the following inequality holds: $$\frac{x^k}{(1-a)(1-a^3)\dots (1-a^{2k-1})}\leq\frac{1}{\psi_a(x)}\quad ?$$
For $a\leq 0.5,$ this is obvious. When $a$ approaches 1, it becomes difficult (for me).
The inequality arises in the study of the dependence of the limit $q$-Bernstein operator on parameter $q.$
It looks like the calculus part is doable (though I still don't like the way I approach it), so let me post the simple (analysis) part.
Consider the symmetric curve $e^{-s}+e^{-t}=1, s,t>0$ (the graph of $f(t)=-\log(1-e^{-t})$ The first observation is that for every $S,T>0$, we have $$ \int_S^\infty f(s)\,ds\ge \int_S^\infty\min(f(s),T)\,ds \\ =-\int_0^S\min(f(s),T)\,ds+\int_0^\infty\min(f(s),T)\,ds \\\ge -ST+\int_0^T f(t)\,dt $$ (we used the symmetry of the curve in the last line).
Now put $h=\log(1/a^2)$. Taking $T=kh$ and $S=\log(1/x)$, we see that what you want to prove (after taking the logarithm of both sides) is pretty much the same except the integrals are replaced by the midpoint Riemann sums with step $h$. We can also assume that the right hand side is made as large as it can be, i.e., $\frac{x}{1-a^{2k+1}}\le 1$ or, which is the same, $e^{-S}\le 1-e^{-T-\frac h2}$.
It will suffice to show that the deviation (down) of the Riemann sum for $\int_0^T$ from the true integral exceeds that for $\int_S^\infty$. The deviation on each step interval of length $h$ is given by the integral of $f''$ against some fixed symmetric positive kernel on $[-h/2,h/2]$ (something like $(1-|x|)^2$ rescaled). Thus, we need to investigate the second derivative of $f$, which is $$ f''(t)=\frac 1{e^t+e^{-t}-2}\,. $$ The first property we shall need is the inequality $$f''(t+s)\le e^{-s}f''(t)\,.$$ Indeed, it rewrites as $e^{t}-2e^{-s}+e^{-t-2s}\ge e^t-2+e^{-t}$. Now just note that the derivative of the LHS with respect to $s$ is $2e^{-s}-2e^{-t-2s}\ge 0$. Denote by $D(T)$ the deviation of the midpoint Riemann sum from $\int_T^\infty f(t)\,dt$. The inequality just proved implies immediately that $$ D(S)\le e^{-S}D(0)\,. $$ Thus, it will suffice to show that when $T$ is an integer multiple of $h$, we have $D(T)\le e^{-T-\frac h2}D(0)$. Since $D(T)\le e^{-T-h}D(h)$, it suffices to check that $D(h)\le e^{-\frac {3h}2}D(0)$. This can be done by pairwise interval comparison if we manage to show that the loss on $[0,h]$ is already at least $e^{\frac{3h}2}$ times as large as the losses on $[h,2h]$ and $[2h,3h]$ combined. After that we will compare only intervals with the shift $2h$ between them and gain $e^{-2h}$ every time.
Due to the symmetry of the kernel, we need the inequality $$ f''(\frac h2-s)+f''(\frac h2+s)\ge e^{\frac{3h}2}[f''(\frac {3h}2-s)+f''(\frac {3h}2+s)+f''(\frac {5h}2-s)+f''(\frac {5h}2+s)], \qquad 0\le s\le \frac h2\,. $$ Since $t\mapsto\frac{e^t-2+e^{-t}}{t^2}$ is increasing for $t>0$ (Taylor expansion), it suffices to prove the same inequality for $t\mapsto t^{-2}$ instead of $f''$. Also $w\mapsto \frac 1{(1-w)^2}+\frac 1{(1+w)^2}$ increases for $w>0$, so we can replace $s$ by $3s$ and $5s$ in the terms with the corresponding multiples of $h$. After that the desired inequality becomes $$ 1\ge e^{\frac{3h}2}[\frac 19+\frac 1{25}] $$, which gives the announced range of $a=e^{-\frac h2}$.
Edit The last part can be done better but we need to know that the kernel is the rescaled version of $(1-|x|)^2$. Again, we can replace $f''(t)$ by $1/t^2$ in the corresponding integral inequality, after which we will need to prove that $$ \int_0^1(1-x)^2\left[\frac 1{(1-x)^2}+\frac 1{(1+x)^2}\right]\,dx \\ \ge 8\int_0^1(1-x)^2\left[\frac 1{(3-x)^2}+\frac 1{(3+x)^2}+ \frac 1{(5-x)^2}+\frac 1{(5+x)^2}\right]\,dx\,. $$ However, the LHS equals $$ 1+\int_0^1\frac{(1-x)^2}{(1+x)^2}\,dx\ge 1+\int_0^1\frac{(1-x)^2}{4}=1+\frac 1{12} $$ while, replacing the expression in the brackets on the right hand side by its maximal value (attained at $x=1$), we can bound the RHS by $$ \frac 83\left[\frac 14+\frac 1{16}+\frac 1{16}+\frac 1{36}\right]=1+\frac{2}{27}\,. $$ This yields the desired result for all $a\ge \frac 12$.
Let me provide a suggestion based on Pietro Majer's expansion.
If it makes it any easier, it suffices to prove the following finite version of my claim: $$\prod_{j=1}^n\frac{x}{1-a^{2j-1}}\leq\sum_{k=0}^{2n}\prod_{j=1}^k\frac{ax}{1-a^{2j}},$$ and perhaps proceed by induction on $n\geq1$.
I'll give this a jumpt start for the base case $n=1$ to prove $$\frac{a^2x^2}{(1-a^2)(1-a^4)}+\frac{ax}{1-a^2}+1-\frac{x}{1-a} =\frac{a^2x^2-(1-a^4)x+(1-a^2)(1-a^4)}{(1-a^2)(1-a^4)}\geq0.$$ It suffices to verify $f(x;b)=bx^2-(1-b^2)x+(1-b)(1-b^2)\geq0$ where $0<b:=a^2<1$.
From the discriminant, $f(x;b)$ has two real roots provided $D=(1-b)(1-b^2)(1-3b)>0$; that means $0<b<\frac13$.
On the other hand, $f(x;b)$ has a minimum at $x_*=\frac{1-b^2}{2b}$. Since we anticipate $0<x_*<1$, it must be that $b^2+2b-1>0$ or $b>\sqrt{2}-1$.
The inequalities $b<\frac13$ and $b>\sqrt{2}-1$ are not compatible. hence, $f(x;b)\geq0$ for $0<a,x<1$.
This was meant to cover the "calculus part" in fedja's initial answer (now completed independently), $0\le a\le a_0:=(1/9+1/25)^{1/3}$, and may possibly be completed to an independent proof.
The functions $\phi(x):=\prod_{j=0}^\infty {1\over 1-a^{2j+1}x}$ and $x^{-k}$ are positive and logarithmically convex on $[0,1]$, therefore such is their product $x^{-k}\phi(x)$ too, which is therefore convex. Moreover $\big(x^{-k}\phi(x)\big)'_{x=1}=-k\phi(1)+\phi'(1)\le0$ iff $k\ge \phi'(1)/\phi(1)=\sum_{j=0}^\infty{a^{2j+1}\over 1-a^{2j+1}},$ in which case $x^{-k}\phi(x)$ has a minimum at $x=1$. We conclude that the improved inequality $$\prod_{j=0}^\infty {1\over 1-a^{2j+1}x}\ge \Big( \prod_{j=0}^\infty {1\over 1-a^{2j+1}}\Big)x^k,\qquad 0\le x\le 1$$ holds whenever $k\ge \sum_{j=0}^\infty{a^{2j+1}\over 1-a^{2j+1}};$ in particular, (by a simple numeric estimates of the sum) for all $k\ge 2$ and $0\le a\le a_0$. (For $k=1$ the inequality can be proven directly as shown in T. Amdeberhan's answer)
