Given a random process (for instance) given by a stochastic differential equation $$\ddot X + \mu \dot X + g X = \xi(t)$$ (where $\xi$ is white noise) or equivalently $$d\mathbf{X} =\mathbf{V} dt$$ $$d\mathbf{V} = (-\mu \mathbf{V} -g \mathbf{X}) dt + d\mathbf{W}$$
($d\mathbf{W} =\xi dt=$ increment of Wiener process) suppose we find a stationary probability distribution
$$P_{\infty}(X,V).$$
Is it correct to interpret $P_{\infty}(0,V)$ as the distribution of velocities of the zero-crossings of the random process? Is it correct to interpret $P_{\infty}(X,0)$ as the distribution of positions of all local extrema of the random process?
(I ask this perhaps obvious question because I find numerically that this almost certainly is not the case.)
First, let me clarify the notation. I will use $P_{\infty}(x,v)$ as the stationary density defined on $\mathbb{R}^2$.
Second, note that unless the density of $(X_0,V_0)$ is given by $P_{\infty}$, the process is not stationary. So you have to make sure that your numerical simulations are either initialized according to the stationary distribution, or are long enough for the transient effects to become irrelevant. The latter is only expected to work if the process is ergodic, i.e. $P_{\infty}$ is unique as a stationary probability measure, and it might take very long unless you can show that the convergence to the stationary probability measure is quick (e.g. exponential).
In order to have a single (long) sample trajectory of $(X,V)$ represent the stationary distribution, you need the process to be ergodic.
As for your question about the distribution of velocities of the zero-crossings of the random process, provided that the process is indeed ergodic, this is given by the conditional density
$$P_{\infty}(v|X=0)=\frac{P_{\infty}(0,v)}{\int_{-\infty}^{\infty}P_{\infty}(0,v')dv'}.$$
Similarly, the distribution of positions of all local extrema of the random process can be expressed as
$$P_{\infty}(x|V=0)=\frac{P_{\infty}(x,0)}{\int_{-\infty}^{\infty}P_{\infty}(x',0)dx'}.$$