My question is about finding the probability distribution of a variable x after running the following process:
x := 0
loop very many times
x -= k*x
x += sampleStandardNormal()
end loop
I have found empirically that the value of x is randomly distributed with a finite variance, provided that $0 < k < 2$. Using a program, I also found empirically that the variance of x is exactly
$$Var(x) = \frac{-1}{k(k-2)}$$
However, I do not know how to derive this result using real math. Furthermore, I do not know how to derive the actual probability density function for x.
(My program found the empirical solution by approximating a Taylor series and then finding the function with the best match.)
Looking at your code, we can see that $x$ is defined as a sequence :
$$X_n=-kX_{n-1}+U_n$$ with $n\geq1$, $X_0=0$, and $U_{k}$'s are standard normal and independent variables.
$$X_1=U_1$$ $$X_2=-kU_1+U_2$$ $$X_3=-kX_2+U_3=-k(-kU_1+U_2)+U_3=k^2U_1-kU_2+U_3$$
By induction :
$$X_n=\sum_{l=0}^{n-1}{(-1)^lk^lV_l}$$
with $V_l=U_{n-l}$ and $l\ge0$
$X_n$ is Gaussian as a linear combination of independent standard normal variables. It remains to derive its mean and variance in order to entirely define the density function .
We clearly have $E(X_n)=0$, and $$E(X_n^2)=\sum_{l=0}^{n-1}\sum_{m=0}^{n-1}{(-1)^{l+m}k^{l+m}E(V_lV_m)}=\sum_{l=0}^{n-1}{k^{2l}}=\frac{1-k^{2n}}{1-k^2}$$