Let $k$ be a positive integer and let $$h(n,k,q)=\frac{1-(1+q^{k})q^{2k(n-1)+1}+q^{2}}{1-q^{2n-1}}h(n-1,k,q)-\frac{(1-q^{k(2n-3)})(1-q^{2k(n-1)})q^2}{(1-q^{2n-1})(1-q^{2n-3})}h(n-2,k,q)$$ with initial values $h(0,k,q)=1$ and $h(1,k,q)=1+q^2{\frac{1-q^{k-1}}{1-q}}.$
Remark: This question is related to question 212926 since $h(n,k,q)$ can also be written as $\frac{\sum\limits_{j = 0}^{2n} {( - 1)}^{ j}q^{j} \binom{2n}{j}_{q^k}}{(q;q^2)_{n}}.$
The rational functions $h(n,k,q)$ converge for $q\to1$ to $h(n,k,1)=k^n$ and for $q\to-1$ to $h(n,2k,-1)=2^n$ and $h(n,2k+1,-1)=1.$
It seems that if $k=2^m$ is a power of $2$ then $h(n,2^m,q)$ is a polynomial with (non- negative) integer coefficients.
My question is: Which methods can be used to prove such facts?