# Simulate time series given temporal auto-correlation functions

Given a random process $$x[n] \in \mathbb{R}$$ (say of length $$N$$) and all correlation functions such as:

\begin{align} \langle x[i]\rangle\\ \langle x[i]x[j]\rangle\\ \langle x[i]x[j]x[k]\rangle\\ \vdots \end{align}

1. Is it possible to simulate a single realization of the random process and if so how? Thinking abstractly about the problem it seems like any trajectory could be mapped onto a vector in $$\mathbb{R}^n$$ and there is some probability distribution on $$\mathbb{R}^n$$. Sampling this distribution and then picking out the corresponding trajectory would then constitute of "simulation" of a single realization that I am looking for. It seems like given information about all of the auto-correlation functions it should be possible to determine the probability distribution. This seems related to the inverse moment problem.
2. What if we truncate and only have up to $$k$$-time correlation function with $$k? My guess is that it is not possible to uniquely generate the probability distribution but it should be possible to (non-uniquely) simulate a process with matching auto-correlation functions up to $$k$$.

Let's model your random process $$x[n]$$ as an auto-regressive (AR) process of order $$P$$, that is, $$x[n] + \sum_{p=1}^P a_px[n-p] = c + \epsilon[n]$$ where $$c$$ is some constant that I will set to 0 in what follows for the sake of simplicity, $$\{a_p\}_{p=1}^P$$ are the parameters of the AR model and $$\epsilon[n]$$ is some white noise of variance $$\sigma^2_\epsilon$$.
Multiplying each side of the above equation by $$x[n-l]$$ for some $$l \ge 0$$ leads us to $$\left(x[n] + \sum_{p=1}^P a_px[n-p]\right)x[n-l] = \epsilon[n]x[n-l].$$ Now, taking the expectation on each side, we obtain the Yule-Walker equations $$C_x[l] + \sum_{p=1}^P a_pC_x[l-p] = \begin{cases} 0 & \text{if} & l > 0, \\ \sigma^2_\epsilon & \text{if} & l = 0 \end{cases}$$ using the fact that $$\epsilon[n]$$ and $$x[n-l]$$ are uncorrelated for $$l > 0$$ (you can see that from the equation of the AR process above). Using this equation with $$l = 0, 1, 2, \dots, P$$, you have $$P+1$$ equations to solve for the $$P + 1$$ unknowns (i.e., $$\{a_p\}_{p=1}^P$$ and $$\sigma^2_\epsilon$$).