Say we have a stochastic process described by a stochastic differential equation (in the Itô sense), and maybe we are able to find an explicit solution of it in terms of deterministic and Itô integrals (and maybe not).
In any of these cases, how can we compute the power spectral density of the process?
For instance, for the Ornstein-Uhlenbeck process, the steps would be
- Find autocorrelation function. This is easily done with the Itô isometry and the properties of the Wiener process.
- Compute Fourier transform
And that's it, because the Wiener-Khinchin theorem assures us that the process autocorrelation and the PSD are a Fourier transform pair as long as the process is wide sense stationary (and the Ornstein-Uhlenbeck process satisfies this). However, The Ornstein-Uhlenbeck process is too nice in reality.
What can be done to obtain the PSD in a more general, non stationary case? What can be then said about the autocorrelation function relation to the PSD? Could you provide any references with theory/worked examples?
For instance, I would really appreciate your help in the not-so-nice (but still nice) second order harmonic oscillator with additive noise term
$$
m\ddot{x} + \gamma\dot{x} + \delta x = \sigma \eta(t)
$$
with
$$
\eta(t) = \frac{\mathrm{d} W_t}{\mathrm{d}t}
$$
(forgive my abuse of notation).
Best Answer
I'll reply myself with a full derivation of the PSD of the harmonic oscillator. It's only half of the answer to my original question; the other half can be found here, after I asked the same question on mathoverflow
The full equation takes the expression $$ m \ddot{x} + \gamma \dot{x} + \delta x = \sigma \eta(t) $$ Performing a change of variables, $x = e^{\frac{-\gamma}{2m}t}x_1$, we get $$ m \ddot{x}_1 + \left(\delta - \frac{\gamma^2}{4m}\right)x_1 = \sigma e^{\frac{\gamma}{2m}t}\eta(t) $$ thus eliminating $\dot{x}$. Setting $a = \frac{\delta}{m} - \frac{\gamma^2}{4m^2}$, $b = \frac{\sigma}{m}$ and rewriting the equation as a first order linear system with $$ X = \begin{pmatrix} x_2\\ v_2 \end{pmatrix} $$ where $v_2 = \dot{x_2}$, we get in Ito's notation \begin{equation} X = \begin{pmatrix}\label{eq:complete} 0 & 1 \\ -a & 0 \end{pmatrix} \cdot X \mathrm{d}t + \begin{pmatrix} 0\\ b e^{\frac{\gamma}{2 m}t} \end{pmatrix} \cdot \mathrm{d} W_t \end{equation} The solution of a linear homogeneous SDE is $$ X_t = e^{\int_0^t A(t)\mathrm{d}t}\cdot X_0 + e^{\int_0^t A(t)\mathrm{d}t}\cdot \int_0^t e^{-\int A(s)\mathrm{d}s}\sigma(s)\mathrm{d}W_s $$ where $A(t)$ is the (generally vector) coefficient of $X$. For this SDE a fundamental matrix solution of the associated homogeneous noise-free system is $$ \Phi(t) = \begin{pmatrix} \cos \sqrt{a}t & \sin \sqrt{a}t/\sqrt{a} \\ -\sqrt{a} \sin\sqrt{a}t & \cos \sqrt{a}t \end{pmatrix} $$ The determinant of this matrix is 1, so its inverse matrix will be $$ \Phi^{-1}(t) = e^{-\int A(\tau)\mathrm{d}\tau} = \det \Phi(t)^{-1}\cdot \begin{pmatrix} \cos \sqrt{a}t & -\sin \sqrt{a}t/\sqrt{a} \\ \sqrt{a} \sin\sqrt{a}t & \cos \sqrt{a}t \end{pmatrix} = \begin{pmatrix} \cos \sqrt{a}t & -\sin \sqrt{a}t/\sqrt{a} \\ \sqrt{a} \sin\sqrt{a}t & \cos \sqrt{a}t \end{pmatrix} $$ and hence we can solve the complete system. We are interested in the first component of $X$, the position (we will only calculate the PSD of $x$, although this can be done considering the full matrix system)
\begin{equation*} \begin{split} x_1(t) = & \begin{pmatrix} \cos \sqrt{a}t & \sin \sqrt{a}t/\sqrt{a} \end{pmatrix} \cdot \begin{pmatrix} x_1(0)\\ v_1(0) \end{pmatrix} + \begin{pmatrix} \cos \sqrt{a}t & \sin \sqrt{a}t/\sqrt{a} \end{pmatrix} \cdot \int_0^{t} b e^{\frac{\gamma}{2m}r}\cdot \begin{pmatrix} -\sin \sqrt{a}r/\sqrt{a}\\ \cos \sqrt{a}r \end{pmatrix} \mathrm{d}W_r \end{split} \end{equation*} Finally, $x_1(t) = e^{\frac{\gamma t}{2m}}x(t)$, so \begin{equation} \begin{split} x(t) = & e^{-\frac{\gamma t}{2m}}\begin{pmatrix} \cos \sqrt{a}t & \sin \sqrt{a}t/\sqrt{a} \end{pmatrix} \cdot \begin{pmatrix} x(0)\\ v(0) + \frac{\gamma}{2m}x(0) \end{pmatrix} +\\ & e^{-\frac{\gamma t}{2m}} \begin{pmatrix} \cos \sqrt{a}t & \sin \sqrt{a}t/\sqrt{a} \end{pmatrix} \cdot \int_0^{t} b e^{\frac{\gamma}{2m}r}\cdot \begin{pmatrix} -\sin \sqrt{a}r/\sqrt{a}\\ \cos \sqrt{a}r \end{pmatrix} \mathrm{d} W_r \end{split} \end{equation}
We see that, after a transient time, only the term depending on $\mathrm{d} W_r$ remains, so the first moment of the process is zero. Now, applying Ito's isometry to calculate the covariance we get $$ \mathbb{E} \left[ \int_0^{T} f(u)\,\mathrm{d} W_u \int_0^S f(v) \,\mathrm{d} W_v \right] = $$ $$ b^2 e^{-\frac{\gamma (t + s)}{2m}} \begin{pmatrix} \cos \sqrt{a}t & \sin \sqrt{a}t/\sqrt{a} \end{pmatrix} \cdot \mathbb{E} \left[ \int_0^{\min(t,s)} e^{\frac{\gamma}{m}u} \begin{pmatrix} \frac{\sin^2 \sqrt{a}u}{a} & -\frac{\sin \sqrt{a}u \cos \sqrt{a}u}{\sqrt{a}}\\ -\frac{\sin \sqrt{a}u \cos \sqrt{a}u}{\sqrt{a}} & \cos^2 \sqrt{a}u \end{pmatrix} \, \mathrm{d} u \right] \cdot \begin{pmatrix} \cos \sqrt{a}s \\ \sin \sqrt{a}s/\sqrt{a} \end{pmatrix} $$
This is a quite uninteresting calculation. After simplification one gets $$ R(\tau) = \frac{b^2 m^2 e^{\frac{-\Gamma |\tau|}{2}}\left( 2 \sqrt{a}m \cos(\sqrt{a}|\tau|) + \gamma \sin(\sqrt{a}|\tau|) \right)}{\sqrt{a}(\gamma^3 + 4a\gamma m^2)} $$ where we have defined the normalized damping constant $\Gamma = \frac{\gamma}{m}$.
From the expression of the autocorrelation we see that $R(t,\tau) = R(\tau)$: therefore, the process is wide-sense stationary and the conditions to apply the Wiener-Khinchin theorem are satisfied. The Fourier transform of this autocorrelation function is the power spectral density $$ S(f) = \frac{16b^2}{\Gamma^4 + 8(a + 4\pi^2f^2)\Gamma^2 + 16(a-4\pi^2f^2)^2} $$ which, after replacing the variables and some rearranging takes the simpler expression $$ S(\omega) = \frac{\sigma^2/m^2}{(\omega_0^2 - \omega^2)^2 + \Gamma^2 \omega^2} $$ where we have replaced the unitary ordinary frequency Fourier transform (in terms of $f$) by the non-unitary angular frequency Fourier transform.