Imagine a two-state Markov chain which hops between the states $\pm 1$ with probability $p<1/2$, so that the autocorrelation function after $k$ steps is
$\rho_k = (2p-1)^k$
If I take an exponential moving average of this series with weighting parameter $\lambda$, what does the distribution of values of the new series look like?
Probably the answer is "gaussian, centered on 0" but what is the variance? Is there a known result that makes this computation trivial?
The answer is certainly not "Gaussian". What you describe is often called Bernoulli convolutions and, even in the independent case (in your setting, $p=1/2$), the limiting object is quite complicated and interesting since it involves some deep number theoretic properties of $\lambda$.
To begin with, let $(X_n)_{n\in \mathbb{Z}}$ denote the $\pm1$-valued Markov chain with probability $p$ of switching states. Let $(Y_n)_{n\in \mathbb{Z}}$ denote the exponential moving average of parameter $\lambda$ with $0<\lambda < 1$ you are interested in, that is, $$ Y_n=\sum_{k=0}^{+\infty}\lambda(1-\lambda)^kX_{n-k}. $$ The Markov chain is centered and has correlation $E(X_nX_{n+k})=(1-2p)^k$ for every integers $n$ and $k\ge0$. (Hence you should check your formula.) From there, one sees that $Y_n$ is centered and one can compute its variance. If I am not mistaken, one finds something like $$ E(Y_n^2)=\frac{1-2p(1-\lambda)/(2-\lambda)}{1+2p(1-\lambda)/\lambda}. $$ The stationary distribution of the moving average is a different story. It is often best described as a measure-valued fixed point problem, as follows. First, $Y_n=X_nY_+$ where $Y_+$ and $X_n$ are independent, and $Y_+$ is distributed like $Y_n$ conditioned on $[X_n=+1]$. Second, $Y_+$ is distributed like $\lambda+(1-\lambda)ZY_+$, where $Z=\pm1$, $P(Z=+1)=1-p$, $P(Z=-1)=p$, and $Z$ independent of $Y_+$.
This indirect description of the stationary distribution is often the most useful tool to get some information on it.
As regards your original "Gaussian" hint, note that conditioning on $(X_{n-k})_{0\le k\le N-1}$ for a given $N$ yields that $Y_n$ is in one of $2^N$ intervals of length $2(1-\lambda)^N$. If $2(1-\lambda)<1$, this simple remark shows that the distribution of $Y_n$ is concentrated on a Cantor set of Lebesgue measure zero (hence this probability distribution does not even have a density with respect to the Lebesgue measure, and it has no atom either). The argument uses only the fact that each $X_n=\pm1$ almost surely and not the structure of the process $(X_n)_n$.
Another easy case is when $\lambda=1/2$. Then, if $(X_n)_n$ is in fact independent ($p=1/2$), one recognises the usual binary expansion of a random number hence $Y_n$ is uniformly distributed on $[-1,1]$, but for every other value of $p$, the distribution of $Y_n$ is concentrated on a subset of $[-1,1]$ of Lebesgue measure zero.
For much more on the stationary distributions of moving averages like the ones which interests you, some starting points could be the paper Sixty years of Bernoulli convolutions by Peres, Schlag and Solomyak, and the book Some random series of functions by Kahane.
The measure $\mu_\lambda$ on an interval whose distribution is given by the random variable $$\sum_{n=1}^\infty \epsilon_n\lambda^n,$$ where the $\epsilon_n$ assume the values 0 and 1 (or $\pm1$) independently with probabilities $(p,1-p)$ is called a biased Bernoulli convolution.
If one assumes $\lambda\in(0,1/2)$, then $\mu_\lambda$ is supported by a Cantor set and is consequently singular. If $\lambda=1/2$, then it is a well known singular measure on $[0,1]$. (It is invariant and ergodic under the doubling map $\tau x=2x\ \bmod 1$ and so is the Lebesgue measure.)
The most interesting case is $\lambda\in(1/2,1)$. Here if $p\in[1/3,2/3]$, then for a.e. $\lambda$ the measure $\mu_\lambda$ is known to be equivalent to the Lebesgue measure. If $\lambda^{-1}$ is a Pisot number, then it is singular. (Which was essentially proved by Erdős in 1939.)This is almost all that is known about these measures.
For more detail see, e.g.,
B. Solomyak, Notes on Bernoulli convolutions
By exponential moving average you mean something like $b_k = C \sum_{i=k}^{\infty} {\frac{a_i}{\lambda^i}}$ ? If $\lambda$ is small enough then this would depend hugely on $a_k$ and so it wouldn't be gaussian right?
