Consider the following stochastic difference equation
\begin{equation} x(t+1) = x(t) + \nu(t+1) \end{equation} where, $x(t)\in\mathrm{R}$ be one dimensional and $\nu(t)$ be the disturbance with an unknown probability distribution. The first $2$ moments of the disturbance is known to us as a function of time and it is also known that its third central moment is non-zero, although its value is not known. Given the distribution of $x(0)$, my task is to estimate the distribution of $x(t)$ over time such that it is correct up to second order.
Here is what I did:
Let us call the mean and covariance of $\nu(t)$ as $m(t)$ and $Q(t)$. At time $t = 0$, I approximated the distribution $\nu(t)$ with a Gaussian distribution having mean and covariance as $m(0)$ and $Q(0)$. Since the difference is linear, the distribution of $x(t); t>0$ will be Gaussian and unscented transformation can be used to find the first two moments of $x(t)$.
In other words, irrespective of the actual distribution of the disturbance, if its first $n\ge2$ moments are known, we can obtain the correct covariance of $x(t)$ by assuming $\nu(t)$ to be normally distributed.
Is this approach correct?
Letting $\mu(t) = E[X(t)]$ we have $\mu(t) = \mu(t-1) + m(t) $
Hence $\mu(t) = m(t)+m(t-1)+m(t-2)+\cdots=\sum_{k=0}^\infty m(t-k)$
Letting $r(t,s)=E[ X(t)X(t-s)]$ , $r(t,0)=\sigma_x^2(t)+\mu(t)^2$ multiplying $x(t) = x(t-1) + \nu(t)$ by itself and assuming $\nu(t)$ is white noise (its samples are uncorrelated), we get
$$r(t,0) =r(t-1,0) + Q(t)+m(t)^2+2 \,m(t-1)\mu(t)$$
Can you go on from here?