A Gaussian AR(1) process with autocorrelation $|\phi|<1$ is strictly stationary, meaning that:

$$F_{X}(x_{t_1+\tau} ,\ldots, x_{t_n+\tau}) = F_{X}(x_{t_1},\ldots, x_{t_n}) \quad \text{for all } \tau,t_1, \ldots, t_n \in \mathbb{R} \text{ and for all } n \in \mathbb{N}$$

As we can see from the "for all $\tau \in \mathcal{R}$", strict stationarity by definition must hold even in finite samples.

However, as a previous question demonstrated, simulated Gaussian AR(1) processes clearly don't achieve stationarity, even weak stationarity, until after some number of draws has elapsed. (My own MWE appears below.) In the comment thread, Glen_b helpfully explains:

...when you get close to 𝜙=1...even though your process is stationary, the sample behaviour in samples can sometimes mimic a nonstationary process for a fair while (and increasingly so as you get closer); you need bigger and bigger samples to see it. Similarly, effects of initial values propagate further and further; with larger 𝜙 it takes a longer warmup before your series behaves like a stationary AR.

Given that the above definition of strict stationarity is clearly a finite-sample rather than an asymptotic property, why does stationarity only seem to hold asymptotically in simulation?

MWE showing finite-sample non-stationarity

Here's a quick example showing that the SD of the draws in increasing in $t$ until approximately the 20$^{th}$ draw. Here, I simulate 1,000 individual time series from the same underlying AR(1) process with $\phi=0.9$ and errors $\epsilon_{t} \sim N(0, 0.5)$. The plot shows the standard deviation across the 1,000 time series of random variable $X_t$, conditional on $t$ (labeled draw.index).


# number of time series to simulate
k = 1000

# number of draws in each series
draws = 100

# simulate Gaussian AR(1)'s with autocorrelation = 0.9
for ( i in 1:k ) {
  .d = data.frame( yi = as.numeric( gen_gts( draws, AR1(phi = 0.9, sigma2 = 0.5) ) ) )
  .d$iterate = i
  .d$draw.index = 1:nrow(.d)
  if ( i == 1 ) d = .d else d = bind_rows(d, .d)

agg = d %>% group_by(draw.index) %>%
  summarise( Mean = mean(yi),
             SD = sd(yi) )

plot(agg$draw.index, agg$SD, type="l")

enter image description here

  • 4
    $\begingroup$ It is helpful to distinguish among (1) a time series process, (2) a time series, and (3) finite-length samples of either. See stats.stackexchange.com/a/126830/919 for some discussion and a reference. $\endgroup$
    – whuber
    Commented Jan 31, 2022 at 17:20
  • $\begingroup$ Thanks, @whuber. These distinctions make sense, but given that what I'm plotting is an average over 1,000 realizations (i.e., time series) of this stochastic process, not a single realization, I'm not seeing how this resolves the question. $\endgroup$
    – half-pass
    Commented Feb 1, 2022 at 13:38
  • 1
    $\begingroup$ Because reading your code requires familiarity with R along with extra packages, It is difficult to determine exactly what you are plotting: a verbal description would be helpful. $\endgroup$
    – whuber
    Commented Feb 1, 2022 at 14:19
  • $\begingroup$ Right, thanks. I've edited to clarify. $\endgroup$
    – half-pass
    Commented Feb 1, 2022 at 15:01

1 Answer 1


strict stationarity by definition must hold even in finite samples

It must hold for a finite number of random variables $\{X_1, X_2, X_3, \dots, X_t\}$ that are part of the stochastic process. $F_X(x_1, x_2, \dots, x_t)$ is the joint CDF of all these $t$ random variables, and the $x_t$s are the possible values a particular random variable $X_t$ can take. So, the definition of stationarity is not about samples $x_t$, but about random variables $X_t$.

However, you don't know the CDF of any $X_t$ in one simulation since one realization of a stochastic process contains one sample $x_t$ per each random variable $X_t$. Thus, in order to say anything about the process, you'll need multiple realizations of it.

If your process is ergodic, you can learn more about it by observing more and more individual realizations $x_t$ from this one ongoing process, so you don't need multiple full realizations of the process. That's why you can "see" stationarity on the plot.

  • $\begingroup$ But this simulation does have multiple realizations (i.e., 1000) of $\{ X_1, \cdots, X_{100} \}$. $\endgroup$
    – half-pass
    Commented Feb 1, 2022 at 13:37
  • $\begingroup$ @half-pass, it has one realization of $X_1$, one realization of $X_2$, one realization of $X_3$ and so on. If you sample the process once more, you'll get the second realization of $X_1$, the second realization of $X_2$ and so on. Each $x_t$ in your sample $\{X_1, X_2, \dots\}$ comes from its own random variable $\{X_1, X_2, \dots\}$. Ah, okay, i get it. This code indeed samples the process multiple times $\endgroup$
    – ForceBru
    Commented Feb 1, 2022 at 13:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.