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_+$.
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 $\lambda (1-\lambda)^N$. 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. 1 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$\lambda (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.