Strang Splitting for Parametric Inference in Second-order Stochastic Differential Equations

TL;DR

This work advances parameter inference for second-order SDEs by employing Strang splitting to construct a tractable pseudo-likelihood that remains effective under hypoellipticity and partial observations. The authors establish consistency and asymptotic normality for four estimators (full and rough objective variants, with complete and partial data), showing variance reduction for the diffusion estimator when the full objective is used with complete observations, and robust behavior under partial observations. They derive comprehensive corrections for finite-difference imputations of unobserved velocity to achieve asymptotic unbiasedness, and demonstrate practical performance via a Kramers oscillator simulation and a Greenland ice-core application. The approach combines theoretical guarantees with computational efficiency, and is complemented by open-source software to enable replication and application to real-world metastable-transition data.

Abstract

We address parameter estimation in second-order stochastic differential equations (SDEs), which are prevalent in physics, biology, and ecology. The second-order SDE is converted to a first-order system by introducing an auxiliary velocity variable, which raises two main challenges. First, the system is hypoelliptic since the noise affects only the velocity, making the Euler-Maruyama estimator ill-conditioned. We propose an estimator based on the Strang splitting scheme to overcome this. Second, since the velocity is rarely observed, we adapt the estimator to partial observations. We present four estimators for complete and partial observations, using the full pseudo-likelihood or only the velocity-based partial pseudo-likelihood. These estimators are intuitive, easy to implement, and computationally fast, and we prove their consistency and asymptotic normality. Our analysis demonstrates that using the full pseudo-likelihood with complete observations reduces the asymptotic variance of the diffusion estimator. With partial observations, the asymptotic variance increases as a result of information loss but remains unaffected by the likelihood choice. However, a numerical study on the Kramers oscillator reveals that using the partial pseudo-likelihood for partial observations yields less biased estimators. We apply our approach to paleoclimate data from the Greenland ice core by fitting the Kramers oscillator model, capturing transitions between metastable states reflecting observed climatic conditions during glacial eras.

Figures (4)

• Figure 1: Simulated trajectory from the Kramers oscillator. The trajectory is simulated using Euler-Maruyama scheme with parameters $\eta_0 = 62.5, a_0 = 297, b_0 = 219, \sigma_0^2 = 9125$ and with the step size $h = 0.02$. Left: Trajectories over time of state $X$ (top) and rate of change $V$ (bottom). Right: Empirical densities (black) and estimated invariant densities with confidence intervals (dark red).
• Figure 2: Simulated velocity vs. approximated velocity from the Kramers oscillator. The trajectory is simulated as in Figure \ref{['fig:sim_data']}. Left: Simulated trajectory over time of velocity $V$ (black) and approximated velocity using forward difference $\Delta_h X$ (green). Right: Empirical densities of $V$ (black), $\Delta_h X$ (green solid) and $\sqrt{3/2}\Delta_h X$ (green dashed).
• Figure 3: Parameter estimates in a simulation study for the Kramers oscillator, eq. \ref{['eq:KramersSDE']}. The color code is the same across all three figures. A) Normalized distributions of parameter estimation errors $(\hat{\bm{\theta}}_N - \bm{\theta}_0) \oslash \bm{\theta}_0$ (where $\oslash$ is the element-wise division) in the complete and partial observation cases, based on 500 simulated data sets with $h = 0.1$ and $N = 5000$. Each column corresponds to a different parameter, while the color indicates the type of estimator. Estimators are distinguished by superscripted objective functions ($\mathrm{F}$ for full and $\mathrm{R}$ for rough). B) Distribution of $\widehat{\sigma}_N^2$ estimators based on 1000 simulations with $h = 0.02$ and $N = 5000$ across different observation settings (complete or partial) and objective function choices (full or rough) using the Strang splitting scheme. The true value of $\sigma^2$ is set to $\sigma_0^2 = 100$. Theoretical normal densities are overlaid for comparison. Theoretical variances are calculated based on $C_{\sigma^2}(\bm{\theta}_0)$, eq. \ref{['eq:Csigma_Kramers']}. C) Median computing time in seconds for one estimation of various estimators based on 500 simulations with $h = 0.1$ and $N = 5000$. Shaded color patterns represent times in the partial observation case, while no color pattern indicates times in the complete observation case.
• Figure 4: Ice core data from Greenland.Left: Trajectories over time (in 1000 years) of the centered negative logarithm of the $\text{Ca}^{2+}$ measurements (top) and forward difference approximations of its rate of change (bottom). Forward differences are multiplied by $\sqrt{3/2}$ to compensate for underestimating the variance. The two vertical dark red lines represent the estimated stable equilibria of the double-well potential function. Green points denote up- and down-crossings of level $\pm 0.6$, conditioned on having crossed the other level. Green vertical lines indicate empirical occupancy estimates in either of the two metastable states. Right: Empirical densities (black) and estimated invariant densities with confidence intervals (dark red), prediction intervals (light red), and the empirical density of a simulated sample from the estimated model (blue). The empirical density of forward differences is also rescaled by $\sqrt{3/2}$ to match the variance of the theoretical density.

Theorems & Definitions (20)

• Definition 2.1
• Remark 1
• Remark 2
• Remark 3
• Theorem 3.1: Consistency of the estimators
• Remark 4
• Theorem 3.2
• Lemma 6.1
• Lemma 6.2
• Remark 5