Let $S_n = \sum_{i=1}^n X_i$, where the $X_i$ are IID with distribution $X$ with $\mathbb{P}(X = 1) = p $ and $ \mathbb{P}(X = 1) = q = 1-p $. I'm interested in computing the variance of $S_T$, where $T$ is the hitting time $$ T = \min\{ n : \lvert S_n \rvert = a \} . $$
I know (from the theory of stopping times, since $T$ is one), that we should have the identity $$ \operatorname{var}(S_T) = \mathbb{E}[T]\operatorname{var}(X) + \operatorname{var}(T) \mathbb{E}[X]^2 , \tag{*} $$ (see bottom of the answer for rigorous calculations in this case), but I can't seem to get this to work doing the calculations from first principles, using recurrence relations: I instead seem to find $$ \operatorname{var}(S_T) = \mathbb{E}[T]\operatorname{var}(X) -\operatorname{var}(T) \mathbb{E}[X]^2 . $$
Details of my computations follow.
First principles
We know, by solving the recurrence relation $$ P_x = p P_{x+1} + q P_{x-1} $$ for $P_x = \mathbb{P}(x+S_n \text{ reaches $a$ first})$, with boundary conditions $ P_a = 1 $ and $ P_{-a} = 0 $, that $$ \mathbb{P}( S_n \text{ reaches $a$ first} ) = \left. \frac{(q/p)^{a+x}-1}{(q/p)^{2a}-1} \right|_{x=0} = \frac{1}{(q/p)^a+1} , $$ or writing $ q/p = z $, $1/(1+z^a)$ (note that $z \in (0,\infty) $, the symmetric case being $z=1$). Similarly, $$ \mathbb{P}( S_n \text{ reaches $a$ first} ) = \frac{z^a}{z^a+1} $$
Since we know that $S_T= \pm a$, we can now use this to compute expectations of function of $S_T$ using the Law of Total Expectation: $$ \mathbb{E}[f(S_T)] = \mathbb{E}[f(S_T) \mid S_T = a] \mathbb{P}(S_T = a) + \mathbb{E}[f(S_T) \mid S_T = -a] \mathbb{P}(S_T = -a) = \frac{f(a) + f(-a) z^a}{1+z^a} . $$ In particular, $$ \mathbb{E}[S_T] = a\frac{1-z^a}{1+z^a} , \qquad \mathbb{E}[S_T^2] = a^2 , $$ and therefore $$ \operatorname{var}(S_T) = \frac{4a^2z^a}{(1+z^a)^2} . $$ So much for the left-hand side.
The easy things to compute on the right-hand side are the quantities relating to $X$. We have $$ \mathbb{E}[X] = \frac{1-z}{1+z} , \qquad \operatorname{var}(X) = \frac{4z}{(1+z)^2} $$ straightforwardly from the mass function. It remains to calculate $\mathbb{E}[T]$ and $\operatorname{var}(T)$.
We can use the recurrence relation $$ E_x = 1 + p E_{x+1} + q E_{x-1} $$ where $E_x = \mathbb{E}[\text{first time for }\lvert x+S_n \rvert = a ] $ and $E_{\pm a} = 0$ to derive the formula $$ \mathbb{E}[T] = a \frac{(1+z)(1-z^a)}{(1-z)(1+z^a)} $$ I'm pretty sure this is correct, because it agrees with the Wald identity for the stopping time $T$, $$ \mathbb{E}[S_T] = \mathbb{E}[X] \mathbb{E}[T] . $$
Now comes the problem. Suppose I try to find $\operatorname{var}(T)$ using $(*)$. Then $$ \operatorname{var}(T) = \frac{ \operatorname{var}(S_T) - \mathbb{E}[T]\operatorname{var}(X) }{ \mathbb{E}[X]^2 } \\ = \frac{(1+z)^2}{(1-z)^2} \left( \frac{4a^2z^a}{(1+z^a)^2} - a \frac{(1+z)(1-z^a)}{(1-z)(1+z^a)} \frac{4z}{(1+z)^2} \right) \tag{**} $$
But this is negative! To see this, one can check a few small values of $a$, or take the limit as $z \to 1$, which is the symmetric case, where we find $ -(2/3)a^2(a^2-1) < 0 $.
In fact, doing the calculation of $\operatorname{var}(T)$ from first principles, either using the generating function or this answer's idea, adapted to our nonsymmetric case, suggests that in fact the answer is actually the negative of $(**)$. I'm not sure there's much point in giving the details since $(**)$ already seems to indicate something's gone wrong.
Stopping time calculation
We can write $$ S_T = \sum_{n=1}^{\infty} X_n \mathbb{1}\{T \geq n\} . $$ Then using that $T$ is a stopping time, $$ \mathbb{E}[S_T] = \sum_{n=1}^{\infty} \mathbb{E}[ \mathbb{E}[ X_n \mathbb{1}\{T \geq n\} \mid X_1 , \dotsc , X_{n-1} ] ] \\ = \sum_{n=1}^{\infty} \mathbb{E}[ \mathbb{1}\{T \geq n\} \mathbb{E}[ X_n \mid X_1 , \dotsc , X_{n-1} ] ] \\ \sum_{n=1}^{\infty} \mathbb{E}[ \mathbb{1}\{T \geq n\} \mathbb{E}[ X_n ]] \\ = \sum_{n=1}^{\infty} \mathbb{P}[T \geq n] \mathbb{E}[ X ] = \mu \mathbb{E}[T] , $$ and following the idea in the comments on this question, $$ S_T^2 = \sum_{n=1}^{\infty} \sum_{m=1}^{\infty} X_n X_m \mathbb{1}\{T \geq \max\{n,m\}\} \\ = \sum_{n=1}^{\infty} X_n^2 \mathbb{1}\{T \geq n\} + 2\sum_{n=1}^{\infty} \sum_{m=1}^{n-1} X_n X_m \mathbb{1}\{T \geq n\} . $$ The first sum has expectation $ \mathbb{E}[X^2] \mathbb{E}[T] = (\sigma^2+\mu^2) \mathbb{E}[T]$. For the second sum, $$ 2\sum_{n=1}^{\infty} \sum_{m=1}^{n-1} \mathbb{E}[ \mathbb{E}[ X_n X_m \mathbb{1}\{T \geq m\} \mid X_1 , \dotsc , X_m ] ] = 2\sum_{n=1}^{\infty} \sum_{m=1}^{n-1} \mathbb{E}[ \mathbb{E}[ X_n X_m \mid X_1 , \dotsc , X_m ] \mathbb{1}\{T \geq n\} ] \\ = 2\sum_{n=1}^{\infty} \sum_{m=1}^{n-1} \mathbb{E}[ \mu^2 \mathbb{1}\{T \geq n\} ] \\ = 2\mu^2 \sum_{n=1}^{\infty} (n-1) \mathbb{P}( T \geq n ) \\ = 2\mu^2 \sum_{n=1}^{\infty} (n-1) \sum_{k=n}^{\infty} \mathbb{P}( T = k ) \\ = 2\mu^2 \sum_{k=1}^{\infty} \mathbb{P}( T = k ) \sum_{n=1}^{k} (n-1) \\ = \mu^2 \sum_{k=1}^{\infty} k(k-1) \mathbb{P}( T = k ) = \mu^2 (\mathbb{E}[T^2]-\mathbb{E}[T]) $$ and so $ \mathbb{E}(S_T^2) = \sigma^2 \mathbb{E}[T] + \mu^2 \mathbb{E}[T^2] $ and by the previous calculation $$ \operatorname{var}(S_T) = \mathbb{E}[T] \operatorname{var}(X) + \operatorname{var}(T) \mathbb{E}[X]^2 . $$
