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.
