What are the properties of this polynomial sequence?
Here is an algebraic formalism that seems to finish off the problem. Let ${\mathbb Z}[x]$ be the ring of integer polynomials in $x$, and let ${\mathbb Z}[x]^{\mathbb N}$ denote the ring of functions $f: {\mathbb N} \to {\mathbb Z}[x]$ that take natural numbers $n \in {\mathbb N} = \{0,1,2,\dots\}$ to polynomials. On this ring we have a shift homomorphism $T: {\mathbb Z}[x]^{\mathbb N} \to {\mathbb Z}[x]^{\mathbb N}$ defined by $Tf(n) := f(n-1)$ (with the convention that $f(n)=0$ for negative $n$). Given a sign pattern $\epsilon \in \{-1,+1\}^{\mathbb N} \subset {\mathbb Z}[x]^{\mathbb N}$, solving the recurrence $$ a(0) = 1$$ $$ a(n) = x a(n-1) + \epsilon(n) a(n-2) \hbox{ for } n \geq 1$$ (with the convention $a(n)=0$ for negative $n$) is equivalent to locating an element $a \in {\mathbb Z}[x]^{\mathbb N}$ obeying the equation $$ (1 - x T - \epsilon T^2) a = \delta$$ where $\delta \in {\mathbb Z}[x]^{\mathbb N}$ is the Kronecker delta, defined by setting $\delta(0)=1$ and $\delta(n) = 0$ for $n > 0$. This has a unique solution that we can write as $a = (1 - xT - \epsilon T^2)^{-1} \delta$, where the inverse is expanded by formal Neumann series.
We introduce the ${\mathbb Z}[x]$-linear operators $A, B_\epsilon: {\mathbb Z}[x]^{\mathbb N} \to {\mathbb Z}[x]^{\mathbb N}$ by $A := xT$ and $B_\epsilon := \epsilon T^2$, thus we are trying to find sign patterns $\epsilon$ such that $(1 - A - B_\epsilon)^{-1} \delta$ takes values in polynomials with coefficients in $\{-1,0,+1\}$.
Call a sign pattern $\epsilon \in \{-1,+1\}^N$ good if it obeys the following properties:
- $\epsilon(2n+2)=-\epsilon(2n+3)$ for all $n$.
- $\epsilon(2^m(2n+3)) = - \epsilon(2^m(2n+1))$ for all $n$ and all $m \geq 1$.
One can make a good sign pattern by the formula $\epsilon(2^m(2n+1)) = \sigma_m (-1)^n$ for all $m \geq 1$ and all $n$ with arbitrary signs $\sigma_m \in \{-1,+1\}$, and then setting $\epsilon(2n+1) = -\epsilon(2n)$ for all $n$.
The point of a good sign pattern is the following: if $\epsilon$ is good, then one has the anti-commutativity property $(A^{2^m} B_\epsilon^{2^m} + B_\epsilon^{2^m} A^{2^m}) f = 0$ whenever $m \geq 0$ and $f \in {\mathbb Z}[x]^{2^{m+1} {\mathbb N}}$ is supported on the multiples $2^{m+1} {\mathbb N}$ of ${\mathbb N}$. Indeed, this is equivalent to the identity $$ \prod_{j=1}^{2^m} \epsilon(2^{m+1} n + 2j) + \prod_{j=1}^{2^m} \epsilon(2^{m+1} n + 2^m + 2j) = 0$$ for any $n$. For $m=0$ this is exactly property 1 of a good sequence. For $m>0$ we cancel off a common factor of $\prod_{j=2^{m-1}+1}^{2^m} \epsilon(2^{m+1}+2j)$ to write the identity as $$ \prod_{j=1}^{2^{m-1}} \epsilon(2^{m+1} n + 2j) + \prod_{j=1}^{2^{m-1}} \epsilon(2^{m+1} (n+1) + 2j) = 0$$ and then observe from property 2 that $\epsilon(2^{m+1}(n+1) + 2j) = \epsilon(2^{m+1} n + 2j)$ when $1 \leq j < 2^{m-1}$ and $\epsilon(2^{m+1}(n+1)+ 2j) = -\epsilon(2^{m+1} n + 2j)$ for $j = 2^{m-1}$.
Using this anticommutativity and induction we obtain the Frobenius type identity $$ (A + B_\varepsilon)^{2^m} f = (A^{2^m} + B_\varepsilon^{2^m}) f$$ whenever $m \geq 0$ and $f \in {\mathbb Z}[x]^{2^m {\mathbb N}}$ (note that $A^{2^m}$ and $B_\varepsilon^{2^m}$ map ${\mathbb Z}[x]^{2^{m+1} {\mathbb N}}$ to ${\mathbb Z}[x]^{2^{m} {\mathbb N}}$). A similar induction then leads to the identity $$ (1-A-B_\epsilon) \prod_{i=0}^{m-1} (1 + A^{2^i} + B_\epsilon^{2^i}) f = (1 - A^{2^m} - B_\epsilon^{2^m}) f$$ whenever $m \geq 0$ and $f \in \mathbb{Z}[x]^{2^m {\mathbb N}}$, where the product is ordered from left to right, thus $$ \prod_{i=0}^{m-1} (1 + A^{2^i} + B_\epsilon^{2^i}) = (1 + A + B_\epsilon) (1 + A^2 + B_\epsilon^2) \dots (1 + A^{2^{m-1}} + B_\epsilon^{2^{m-1}}).$$ Specialising to $f=\delta$, applying $(1-A-B_\varepsilon)^{-1}$ and then sending $m$ to infinity we obtain the formula $$ (1-A-B_\epsilon)^{-1} \delta = \prod_{i=0}^{\infty} (1 + A^{2^i} + B_\epsilon^{2^i}) \delta$$ (where the product converges pointwise). One can check that every term in this pointwise product gives a different monomial located at a different point with coefficient $\pm 1$, so this indeed gives a sequence $a(n) \in {\mathbb Z}[x]$ for $n \in {\mathbb N}$ with the stated properties, with the explicit form $$ a = (1 + A + B_\epsilon) (1 + A^2 + B_\epsilon^2) (1 + A^4 + B_\epsilon^4) \dots \delta.$$
To give some sense of this formula, if we write $a$ as a sequence $a(0),a(1),a(2),\dots$, then $$ \delta = 1, 0, 0, \dots$$ $$ (1+A+B_\epsilon) \delta = 1, x, \epsilon(2), 0, \dots$$ $$ (1+A+B_\epsilon) (1+A^2+B_\epsilon^2) \delta = 1, x, \epsilon(2) + x^2, x^3, \epsilon(2)\epsilon(4) + \epsilon(2) x^2, \epsilon(2)\epsilon(4) x, \epsilon(2)\epsilon(4)\epsilon(6),0, \dots$$ and one converges to $$ a(0) = 1$$ $$ a(1) = x$$ $$ a(2) = \epsilon(2) + x^2$$ $$ a(3) = x^3$$ $$ a(4) = \epsilon(2) \epsilon(4) + \epsilon(2) x^2 + x^4 $$ $$ a(5) = \epsilon(2) \epsilon(4) x + x^5 $$ $$ a(6) = \epsilon(2)\epsilon(4) \epsilon(6) + \epsilon(6)x^4 + x^6 $$ $$ \dots$$ which is a reparameterisation of the previous solutions.
One can show that these are in fact the only sequences of polynomials $a$ that maintain their coefficients in $-1,+1$, but the proof of this uniqueness is a somewhat tedious induction and this answer is already quite long, so I'll leave it as an exercise.
No proofs, just experimental facts (hence cw):
Let $a_{n+1}=xa_n+\varepsilon_{n+1}a_{n-1}$, with $\varepsilon_k=\pm1$; then $\varepsilon_{2^k}$ might be arbitrary, while all others are uniquely determined by them. For the first few cases we have $$ \begin{aligned} \varepsilon_{3}&=-\varepsilon_{2}\\ \varepsilon_{5}&=-\varepsilon_{4}\\ \varepsilon_{6}&=-\varepsilon_{2}\\ \varepsilon_{7}&=\varepsilon_{2}\\ \varepsilon_{9}&=-\varepsilon_{8}\\ \varepsilon_{10}&=\varepsilon_{2}\\ \varepsilon_{11}&=-\varepsilon_{2}\\ \varepsilon_{12}&=-\varepsilon_{4}\\ \varepsilon_{13}&=\varepsilon_{4}\\ \varepsilon_{14}&=-\varepsilon_{2}\\ \varepsilon_{15}&=\varepsilon_{2}\\ \varepsilon_{17}&=-\varepsilon_{16}\\ \varepsilon_{18}&=\varepsilon_{2}\\ \varepsilon_{19}&=-\varepsilon_{2}\\ \varepsilon_{20}&=\varepsilon_{4}\\ \varepsilon_{21}&=-\varepsilon_{4}, \end{aligned} $$ the polynomials being, respectively, $$ \begin{aligned} a_1&=x \\ a_2&=x^2+\varepsilon_{2} \\ a_3&=x^3 \\ a_4&=x^4+\varepsilon_{4}x^2 +\varepsilon_{2} \varepsilon_{4} \\ a_5&=x^5+\varepsilon_{2} \varepsilon_{4}x \\ a_6&=x^6-\varepsilon_{2}x^4 -\varepsilon_{4} \\ a_7&=x^7 \\ a_8&=x^8+\varepsilon_{8}x^6 -\varepsilon_{2} \varepsilon_{8}x^4 -\varepsilon_{4} \varepsilon_{8} \\ a_9&=x^9-\varepsilon_{2} \varepsilon_{8}x^5 -\varepsilon_{4} \varepsilon_{8}x \\ a_{10}&=x^{10}+\varepsilon_{2}x^8 -\varepsilon_{8}x^4 - \varepsilon_{4} \varepsilon_{8}x^2-\varepsilon_{2}\varepsilon_{4} \varepsilon_{8} \\ a_{11}&=x^{11}-\varepsilon_{4} \varepsilon_{8}x^3 \\ a_{12}&=x^{12}-\varepsilon_{4}x^{10} -\varepsilon_{2} \varepsilon_{4}x^8 + \varepsilon_{8}x^2+\varepsilon_{2}\varepsilon_{8} \\ a_{13}&=x^{13}-\varepsilon_{2} \varepsilon_{4}x^9 +\varepsilon_{2} \varepsilon_{8}x \\ a_{14}&=x^{14}-\varepsilon_{2}x^{12} +\varepsilon_{4}x^8 - \varepsilon_{8} \\ a_{15}&=x^{15} \\ a_{16}&=x^{16}+\varepsilon_{16}x^{14} -\varepsilon_{2} \varepsilon_{16}x^{12} +\varepsilon_{4} \varepsilon_{16}x^8 -\varepsilon_{8} \varepsilon_{16} \\ a_{17}&=x^{17}-\varepsilon_{2} \varepsilon_{16}x^{13} + \varepsilon_{4} \varepsilon_{16}x^9 - \varepsilon_{8} \varepsilon_{16}x \\ a_{18}&=x^{18}+\varepsilon_{2}x^{16} -\varepsilon_{16}x^{12} +\varepsilon _{4} \varepsilon_{16}x^{10} +\varepsilon_{2}\varepsilon_{4} \varepsilon_{16}x^8 - \varepsilon_{8} \varepsilon_{16}x^2 -\varepsilon_{2} \varepsilon_{8} \varepsilon_{16} \\ a_{19}&=x^{19}+\varepsilon_{4} \varepsilon_{16}x^{11} - \varepsilon_{8} \varepsilon_{16}x^3 \\ a_{20}&=x^{20}+\varepsilon_{4}x^{18} +\varepsilon_{2} \varepsilon_{4}x^{16} +\varepsilon_{16}x^{10} +\varepsilon_{2}\varepsilon_{16}x^8 -\varepsilon_{8} \varepsilon_{16}x^4 -\varepsilon_{4}\varepsilon_{8} \varepsilon _{16}x^2 -\varepsilon_{2}\varepsilon_{4}\varepsilon_{8} \varepsilon_{16} \\ a_{21}&=x^{21}+\varepsilon_{2} \varepsilon_{4}x^{17} +\varepsilon_{2} \varepsilon_{16}x^9 - \varepsilon_{8} \varepsilon_{16}x^5- \varepsilon_{2} \varepsilon_{4} \varepsilon_{8} \varepsilon_{16}x \end{aligned} $$
More experimental facts. I think it should not be too hard prove them by induction, taking into account that we know exactly when the coefficients of $a_n$ are non-zero (see my comment to მამუკა ჯიბლაძე's answer).
Given a sequence $\epsilon_{2^n}\in\{-1,+1\}$ for $n\ge1$, let's define inductively the sequence $(\epsilon_j)_{j\ge1}$ for all positive integers according to the recurrence relations:
$$\epsilon_{2^n+1}=-\epsilon_{2^n}\qquad\text{for }n\ge1 $$
$$\epsilon_{2^n+j}=\epsilon_{2^n-j+1}\qquad\text{for }2\le j< 2^n\ ,$$
and also define
$$ \eta_j:=\begin{cases} (-1)^{ j\over 2},& \text {for even }\ j \\ \\ (-1)^{ j-1\over 2}\epsilon_{j+2}, &\text {for odd } j \ .\\ \end{cases} $$ and $$ \delta_j:=\begin{cases} -1,& \text {for }\ j=0 \\ \\ 1, &\text {for } j\neq0 \ .\\ \end{cases} $$
Then (experimentally) the polynomials $a_k= a_k(x,\epsilon_2,\epsilon_4,\dots,\epsilon_{2^n})$ are determined inductively by $a_{-1}:=0, a_0:=1$, and by the recurrence relation, for any $0\le n$ and any $0\le j< 2^n\ $
$$a_{2^n+j}=x^{2^n}b_j+ {\delta_{2^{n-1}-1-j}}\ \eta_j\ \epsilon_{2^n}\ a_{2^n-2-j},$$
where
$$b_j:=a_j(x,\epsilon_2,\epsilon_4,\dots,\epsilon_{2^{n-2}},-\epsilon_{2^{n-1}})\ .$$
For, instance the above gives $$a_{51}={x}^{51}-\epsilon_{{4}}\epsilon_{{16}}{x}^{43}+ \epsilon_{{8}} \epsilon_{{16}}{x}^{35}+ \epsilon_{{4}}\epsilon_{{32}}{x}^{11}-\epsilon_ {{8}}\epsilon_{{32}}x^3$$