Why, empirically, does strict stationarity only hold *asymptotically* for an AR(1) process?

Question

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).

library(dplyr)
library(tidyverse)
library(simts)

# 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")

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.