Experimenting with a CAS suggests an induction. In order to handle the induction, we need to consider the forms of the numbers involved. $\frac{4^m-1}{3} = 1 + 2^2 + 2^4 + \cdots + 2^{2m-2}$ alternates $1$ and $0$ bits. The map $2n + 1 \to n - 2^{f(n)}$ drops the rightmost bit (which is $1$) and clears the next least significant bit. The map $2n \to n$ drops the rightmost bit (which is $0$). And the map $2n \to 2n - 2^{f(n)}$ moves the least significant bit right one place. So the numbers which occur in the evaluation tree of $\frac{4^m-1}{3}$ are of the form $2^i + 2^j + 2^{j+2} + \cdots + 2^{j+2k}$ where $i \le j - 2$ and $k \ge 0$; and (when we get down to one bit) the form $2^i$.
Starting with the simplest case, when $i > 0$, $a(2^i) = pa(2^{i-1}) + qa(2^{i-1})$ so by induction $a(2^i) = (p+q)^i$.
For the more general case, let $A(i,j,k) = a(2^i + 2^j + 2^{j+2} + \cdots + 2^{j+2k})$ subject to the aforementioned constraints on $i,j,k$.
For $i > 0$, we have an even argument and use the second case of the recursion:
\begin{eqnarray*}
A(i,j,k) &=& a(2^i + 2^j + 2^{j+2} + \cdots + 2^{j+2k}) \\
&=& pa(2^{i-1} + 2^{j-1} + 2^{j+1} + \cdots + 2^{j+2k-1}) + qa(2^{i-1} + 2^j + 2^{j+2} + \cdots + 2^{j+2k}) \\
&=& pA(i-1, j-1, k) + qA(i-1, j, k)
\end{eqnarray*}
Then it's a trivial proof by induction that $$A(i,j,k) = \sum_{u=0}^i \binom{i}{u} p^u q^{i-u} A(0, j-u, k)$$
For $i = 0$, we have an odd argument and use the third case of the recursion:
\begin{eqnarray*}
A(0,j,k) &=& a(1 + 2^j + 2^{j+2} + \cdots + 2^{j+2k}) \\
&=& a(2^{j+1} + \cdots + 2^{j+1+2(k-1)}) \\
&=& \begin{cases}
a(0) & \textrm{if } k=0 \\
a(2^{j+1}) & \textrm{if } k=1 \\
A(j+1, j+3, k-2) & \textrm{otherwise}
\end{cases} \\
&=& \begin{cases}
1 & \textrm{if } k=0 \\
(p+q)^{j+1} & \textrm{if } k=1 \\
\sum_{u=0}^{j+1} \binom{j+1}{u} p^u q^{j+1-u} A(0, j+3-u, k-2) & \textrm{otherwise}
\end{cases}
\end{eqnarray*}
Further CAS experimentation suggests that the theorem we need to prove is $$A(0, j, 2v-1) = A(0, j, 2v) = \left(p\frac{q^v-1}{q-1} + q^v\right)^{j+1} A(0, 2, 2v-2)$$
The first of those equalities is easy: $$A(0,j,2) = \sum_{u=0}^{j+1} \binom{j+1}{u} p^u q^{j+1-u} A(0, j+3-u, 0) = \sum_{u=0}^{j+1} \binom{j+1}{u} p^u q^{j+1-u} = (p+q)^{j+1}$$ so $A(0,j,1) = A(0,j,2)$ and since the only occurrence of $k$ in the third case is in the parameter $k-2$ the rest follows by induction on $v$.
The second equality is the interesting one. Again, by induction on $v$:
\begin{eqnarray*}
A(0,j,2v) &=& \sum_{u=0}^{j+1} \binom{j+1}{u} p^u q^{j+1-u} A(0, j+3-u, 2v-2) \\
&=& \sum_{u=0}^{j+1} \binom{j+1}{u} p^u q^{j+1-u} \left(p\frac{q^{v-1}-1}{q-1} + q^{v-1}\right)^{j+4-u} A(0, 2, 2v-4) \\
&=& \left[ \sum_{u=0}^{j+1} \binom{j+1}{u} p^u \left(pq\frac{q^{v-1}-1}{q-1} + q^v\right)^{j+1-u} \right] \color{Blue}{\left(p\frac{q^{v-1}-1}{q-1} + q^{v-1}\right)^3 A(0, 2, 2v-4)} \\
&=& \left( p + pq\frac{q^{v-1}-1}{q-1} + q^v \right)^{j+1} \color{Blue}{A(0, 2, 2v-2)} \\
&=& \left(p\frac{q^v-1}{q-1} + q^v\right)^{j+1} A(0, 2, 2v-2)
\end{eqnarray*}
as desired.
The answer to the original question now drops out: $a\left(\tfrac{4^m-1}{3}\right) = A(0, 2, m-2)$, but $A(0, 2, k)$ has the form of a cube times $A(0, 2, k-2)$ so by induction it's always a cube.
As proved in this answer, if we represent $n$ as $n=2^{t_1}(1+2^{t_2+1}(1+\dots(1+2^{t_m+1}))\dots)$ with $t_j\geq 0$, then
\begin{split}
a(n) = P(\ell)^{t_1}\prod_{j=1}^{\ell-1} P(\ell-j)^{t_{2j}+t_{2j+1}+1},
\end{split}
where $\ell:=\left\lfloor\frac{m+1}2\right\rfloor$ and
$$P(k) := q^k+p\frac{q^k-1}{q-1}.$$
Notice that this does not depend on $t_m$ when $m$ is even.
Grouping the summands in $s(n)$ by the number of unit bits (very much like in this other answer), for $n\geq 1$ we have
\begin{split}
s(n) &= \sum_{m=1}^n\ \sum_{t_1+t_2+\dots+t_{m}=n-m}\ P(\lfloor(m+1)/2\rfloor)^{t_1} \prod_{j=1}^{\lfloor(m-1)/2\rfloor} P(\lfloor(m+1)/2\rfloor-j)^{t_{2j}+t_{2j+1}+1}\\
&= \sum_{m=1}^n\ [x^{n-m}]\ \frac1{1-P(\lfloor(m+1)/2\rfloor)x}\cdot\prod_{j=1}^{\lfloor(m-1)/2\rfloor} \frac{P(j)}{(1-P(j)x)^2}\cdot\frac{1+\frac{(-1)^m-1}2x}{1-x}\\
&= [x^n]\ \sum_{m=1}^\infty\ \frac{x^m}{1-P(\lfloor(m+1)/2\rfloor)x}\cdot\prod_{j=1}^{\lfloor(m-1)/2\rfloor} \frac{P(j)}{(1-P(j)x)^2}\cdot\frac{1+\frac{(-1)^m-1}2x}{1-x}.
\end{split}
That is, the generating function for $s(n)$ is
$$\sum_{n\geq 0} s(n)x^n = 1+\frac1{1-x}\sum_{m=1}^\infty\ \frac{x^m(1+\frac{(-1)^m-1}2x)}{1-P(\lfloor(m+1)/2\rfloor)x}\cdot\prod_{j=1}^{\lfloor(m-1)/2\rfloor} \frac{P(j)}{(1-P(j)x)^2}.$$
Combining terms for $m=2k-1$ and $m=2k$, we can rewrite it as
$$\sum_{n\geq 0} s(n)x^n = 1+\frac{1}{1-x}\sum_{k=1}^\infty\ \frac{x^{2k-1}}{1-P(k)x}\cdot\prod_{j=1}^{k-1} \frac{P(j)}{(1-P(j)x)^2},$$
which is quite close to the conjectured formula.
Best Answer
If $n\in \mathbb N$ has binary expansion $n=\sum_{k\in S}2^k$ for a finite subset $S\subset \mathbb N$, then $$ \sum_{j=0}^n\Big[\big( {n \atop j}\big)\text{mod}\, 2 \Big]x^jy^{n-j}=\prod_{k\in S} (x^{2^k}+y^{2^k})$$
More generally, but not needed here (see below): for sums of $r$ indeterminates $x_i$, the binary coefficient homogeneous polynomial obtained as expansion of $\Big( \sum_{i=1}^r x_i\Big)^n$ with coefficients reduced $\text{mod}\,2$ has the factorisation (in $\mathbb Z[x_1,\dots,x_r]$)
$$ \sum_{\alpha\in\mathbb N^r,\, |\alpha|=n}\Big[\big( {n \atop \alpha }\big)\text{mod}\,2 \Big]x^\alpha=\prod_{k\in S} \sum_{i=1}^r x_i^{2^k}. $$
In particular (for $x=1$ and $y=m$) by definition of $f$ and $g$, this implies immediately the relation $a(n)=(1+m^{2^f})a(g(n))$, just because $f\in S$ and $\{2^k\}_{ k\in S\setminus\{f\}}$ are the powers of the binary expansion of $n-2^f$.
Details on the above factorization: immediately by induction on $k$ we have $$\Big( \sum_{i=1}^r x_i\Big)^{2^k}= \sum_{i=1}^r x_i^{2^k}\text{mod}\,2.$$ Then, from $n=\sum_{k\in S}2^k$ one has $$\Big( \sum_{i=1}^r x_i\Big)^n=\Big( \sum_{i=1}^r x_i\Big)^{\sum_{k\in S}2^k}= \prod_{k\in S} \Big( \sum_{i=1}^r x_i\Big)^{2^k} = \prod_{k\in S} \Big( \sum_{i=1}^r x_i^{2^k}\Big) \,\text{mod}\,2.$$
On the other hand, again by induction (on $|S|$) if we do the expansion $$\prod_{k\in S}\sum_{i\in [r]}c(k,i)=\sum_{\iota\in [r]^S}\prod_{k\in S}c(k,\iota(k))$$ for $\prod_{k\in S} \Big( \sum_{i=1}^r x_i^{2^k}\Big)$, we see that all the $r^{|S|}$ terms of the expansion are different from each other, so that it is a $\{0,1\}$-coefficient polynomial, namely $\displaystyle \sum_{\alpha\in\mathbb N^r,\, |\alpha|=n}\Big[\big( {n \atop \alpha }\big)\text{mod}\,2 \Big]x^\alpha$.
Notation: As usual, for $x=(x_1,\dots,x_r)\in \mathbb R^r$ and $\alpha=(\alpha_1,\dots,\alpha_r)\in \mathbb N^r$ we denote $x^\alpha:=x_1^{\alpha_1}\cdots x_r^{\alpha_r}$, $|\alpha|:=\alpha_1+\dots+\alpha_r$, and $\big({|\alpha| \atop \alpha }\big):=\frac{|\alpha|!}{\alpha_1!\cdots\alpha_r!}$.