<figure><img src="image_1.png" /></figure>

HAL
open science

# Spatial constraints improve filtering of measurement noise from animal tracks

Alexandre Delporte, Susanne Ditlevsen, Adeline Samson

► To cite this version:

Alexandre Delporte, Susanne Ditlevsen, Adeline Samson. Spatial constraints improve filtering of measurement noise from animal tracks. 2025. hal-05390302

HAL Id: hal-05390302

https://hal.science/hal-05390302v1

Preprint submitted on 1 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
Spatial constraints improve filtering of
measurement noise from animal tracks

Alexandre Delporte*¹, Susanne Ditlevsen², and Adeline Samson¹

¹Laboratoire Jean Kuntzmann, Université Grenoble-Alpes

²Department of Mathematical Sciences, University of Copenhagen

Abstract

Advances in tracking technologies for animal movement require new statistical tools to better exploit the increasing amount of data. Animal positions are usually calculated using the GPS or Argos satellite system and include potentially complex non-Gaussian and heavy-tailed measurement error patterns. Errors are usually handled through a Kalman filter algorithm, which can be sensitive to non-Gaussian error distributions. In this paper, we introduce a realistic latent movement model through an underdamped Langevin stochastic differential equation (SDE) that includes an additional drift term to ensure that the animal remains in a known spatial domain of interest. This can be applied to aquatic animals moving in water or terrestrial animals moving in a restricted zone delimited by fences or natural barriers. We demonstrate that the incorporation of these spatial constraints into the latent movement model improves the accuracy of filtering for noisy observations of the positions. We implement an Extended Kalman Filter as well as a particle filter adapted to non-Gaussian error distributions. Our filters are based on solving the SDE through splitting schemes to approximate the latent dynamic.

# 1 Introduction

Animal tracking allows us to uncover movement patterns of wild species in their natural environment. Having accurate trajectories is of particular interest to better infer space use and define conservation areas [Hays et al., 2021]. Position data typically comes with measurement error whose distribution depends on the type of device attached to the animal and the method used to estimate the position. One of the most common methods for tracking animals at sea is via the Argos telemetry system, with more than 40000 individuals tracked since 2007 [Jønsen et al., 2020]. In this system, the estimation of the animal position, which relies on Doppler shift calculations, can exhibit heavy-tailed and anisotropic errors of several kilometers [Hoenner et al., 2012]. Alternatively, new Fastloc-GPS loggers deployed in wild animals have provided more accurate positions with heavy-tailed errors ranging from ten to a few hundreds meters [Wensveen et al., 2015].

To estimate true positions from noisy observations, it is critical to accurately model the error distribution, as well as the latent dynamics of motion. In the literature, the estimation of

<footnote>*alexandre.delporte@univ-grenoble-alpes.fr</footnote>
the true position given the observations is typically performed using a Kalman filter algorithm
[Jonsen et al., 2020; Michelot et al., 2021]. This algorithm may be sensitive to non-Gaussian
errors and requires to pre-filter the data to ensure that it does not deviate too much from
Gaussianity [Patterson et al., 2010].

In addition to a precise error model, the latent dynamics must also be specified correctly for the algorithm to provide good filtering estimates of the trajectory. The dynamics of motion is often influenced by complex landscape boundaries that the animals cannot cross. Brost et al. [2015] showed that knowledge of environmental boundaries can be leveraged to improve inference in the context of resource selection functions. They modelled the latent dynamics using a point process, and the observations are Argos positions with heavy-tailed errors. Although stochastic differential equations (SDE) are becoming more important in movement ecology, no work has been done yet on the influence of these spatial constraints when the latent process is modelled by an SDE.

Recently, Michelot and Hanks [2025] applied an under-damped stochastic Langevin equation linking motion to environmental covariates in the context of animal movement and provided parameter inference tools, although without boundary constraints nor measurement noise. Cholaquidis et al. [2020] used a reflected Brownian motion with drift to estimate home ranges and core areas of the animals, but they did not consider measurement error in the observations. Inspired by this work, Etoré et al. [2025] defined a second-order reflected SDE that includes directional persistence and is therefore more suitable for modelling animal movement than a reflected Brownian motion. However, the paper focuses on the theoretical properties of the model, no simulation, filtering or inference procedures have been developed yet. Hanks et al. [2017] defined a reflected linear SDE and inferred the movement parameters based on Argos data with non-Gaussian location errors. In their model, the spatial constraints are included by projecting the positions that fall outside the domain onto its boundary. The SDE is solved numerically by replacing the velocity with the finite differences between consecutive positions inside an Euler-Maruyama step. Delporte et al. [2025] also modeled spatial constraints through a smooth rotation of the velocity as positions approach the domain's boundary. They only considered a linear dynamic with Gaussian measurement error and inferred movement parameters through a classical Kalman Filter algorithm.

In contrast, in this paper, we propose a new penalized version of the Langevin SDE that smoothly enforces spatial constraints through a simple drift correction term rather than a projection onto the boundary or a rotation. This approach ensures that the position remains close to the ecological domain of interest. Moreover, our model goes beyond linear dynamics: we specify a non-linear drift as a mixture of Gaussians and solve the resulting SDE using splitting schemes adapted to that nonlinearity. This allows for efficient simulation of the constrained dynamics. In addition, we implement filters for state estimation, accounting for domain constraint: first an extended Kalman filter under Gaussian errors, and second a particle filter to handle non-Gaussian measurement errors such as heavy-tailed Student's t-distributions.

Through a thorough simulation study, we compare the performances of naive Kalman
filtering, extended Kalman filtering, and particle filtering, both with and without spatial
penalization. We demonstrate that incorporating domain knowledge substantially improves
filtering accuracy, especially when the observation frequency is high or the measurement error
distribution deviates from Gaussianity.
The remainder of the paper is structured as follows. In Section 2, we introduce the penalized Langevin SDE and discuss the relationship to other existing models. In Section 3, we describe the measurement error models that we consider throughout the paper. Section 4 presents numerical schemes for simulating the dynamics. Section 5 details the filtering algorithms under different measurement error specifications. It includes Kalman filtering and Extended Kalman filtering for Gaussian measurement error, as well as particle filtering for Student errors. Section 6 is dedicated to a simulation study in different scenarios. We vary the measurement error distribution as well as the time step between consecutive observations and compare the performances of the filtering algorithms. We conclude with a discussion of potential applications and extensions.

## 2 Penalized Langevin SDE in a bounded domain

Let $\mathcal{D} \subset \mathbb{R}^2$ be an open bounded domain. We consider the following penalized Langevin SDE:

$$ \begin{cases} dX(t) = V(t)dt \\ dV(t) = -AV(t) - \nabla H(X(t))dt + \sigma dW(t) - \beta_{\lambda}(X(t))dt \end{cases} \qquad (2.1) $$

where $A$ is a positive definite matrix parametrised as

$$ A = \begin{pmatrix} c & -\omega \\ \omega & c \end{pmatrix} $$

with $c > 0$ a damping (or autocorrelation) parameter and $\omega$ an angular velocity parameter. $X = (X_1 \ X_2)^\top$ represents the horizontal position of the animal, $V = (V_1 \ V_2)^\top$ the horizontal velocity of the animal, and $W(t), t \ge 0$, is a two-dimensional Brownian motion.

The function $H$ defines the potential surface that controls areas of attraction (or repulsion) within the domain $\mathcal{D}$. Examples of potential surfaces $H$ can be found in Preisler et al. [2004], we use a mixture of Gaussian surfaces below.

The penalization term $\beta_{\lambda}(X(t))$ is given by

$$ \beta_{\lambda}(x) = \frac{x - \pi(x)}{\lambda} \quad \text{and} \quad \pi(x) = \underset{x' \in \bar{\mathcal{D}}}{\operatorname{argmin}} \|x - x'\|, \quad x \in \mathbb{R}^2, $$

where $\pi(x)$ is the projection of the position $x$ into $\bar{\mathcal{D}}$. The parameter $\lambda$ controls the strength of the penalization applied to keep the process inside $\mathcal{D}$. When the domain is convex, the penalty simplifies to [Liu, 1995]

$$ \beta_{\lambda}(x) = \frac{1}{2\lambda} \nabla_x \inf_{x' \in \bar{\mathcal{D}}} \|x - x'\|^2. $$

This model differs from the one proposed in Michelot and Hanks [2025] in two ways: the penalization that restricts the position $X$ to the domain $\mathcal{D}$, and the rotational component in the velocity via the parameter $\omega$. The model proposed by [Gurarie et al., 2017] corresponds to a flat surface $H$. Here, we extend their model with the penalty, corresponding to

$$ \begin{cases} dX(t) = V(t)dt \\ dV(t) = -AV(t)dt + \frac{2\nu}{\sqrt{\pi\tau}} dW_t - \beta_{\lambda}(X(t))dt \end{cases} \qquad (2.2) $$
where $\tau = 1/c$ is the movement persistence and $\nu = \sigma\sqrt{\pi\tau}/2$ controls the velocity norm.

In the sequel, we consider a potential surface $H$ in eq. (2.1) defined as a mixture of Gaussian potentials:

$$ \forall x \in \mathbb{R}^2: \quad H(x) = -\sum_{j=1}^{J} H_j(x) \quad \text{with} \quad H_j(x) = \alpha_j \exp(-(x - x_j^*)^\top B_j(x - x_j^*)) $$

with weights $\alpha_j$, precision matrices $B_j$ and centers $x_j^* \in \mathcal{D}$ for $j \in \{1, \dots, J\}$. We denote $e_j(x) = \exp(-(x - x_j^*)^\top B_j(x - x_j^*))$. The gradient and Hessian matrices are computed as follows, see details in Appendix A:

$$ \nabla H(x) = \sum_{j=1}^{J} 2\alpha_j B_j(x - x_j^*) e_j(x) \quad (2.3) $$

$$ D^2 H(x) = \sum_{j=1}^{J} 2\alpha_j e_j(x) B_j \left( I_2 - 2(x - x_j^*)(x - x_j^*)^\top B_j \right). \quad (2.4) $$

From now on, we also suppose that $\mathcal{D}$ is a polygon with vertices $(v_k)_{1 \le k \le p}$. This polygon can represent an area delimited by fences that the animal cannot cross [Brillinger, 2003], as well as water areas surrounded by coastlines that restrict the movement of marine mammals [Hanks et al., 2017]. When the position is within the domain $\mathcal{D}$, then we simply have $\pi(x) = x$, and the penalization is zero. Otherwise, there exists a vertex $v_k$ such that $\pi(x) = v_k + \gamma(x)(v_{k+1} - v_k)$ where $\gamma(x) = \max(0, \min(1, \frac{\langle x - v_k, v_{k+1} - v_k \rangle}{l^2}))$ and $l = \|v_{k+1} - v_k\|$. In theory, the probability that $\pi(x) \in \{v_k, v_{k+1}\}$ is zero, in which case we have $\gamma(x) = \frac{\langle x - v_k, v_{k+1} - v_k \rangle}{l^2}$. Hence, the projection operator $\pi$ is piecewise linear in $x$.

## 3 Measurement error models

In this section, we describe the statistical models used to represent measurement errors in the observed animal positions. Accurate modeling of these errors is essential for reliable state estimation, as their distribution strongly depends on the tracking technology and can substantially affect the performance of filtering algorithms.

In practice, the velocity process $V$ is not observed, and only the position process $X$ is observed with measurement error. We denote $(Y_k)_{1 \le k \le n}$ the discrete-time noisy observations of the positions. The most common assumption for animal tracks is Gaussian error distributions [Johnson et al., 2008]. The reason it is so popular is that under this assumption, filtering as well as parameters estimation are easily done with the Kalman filter [Michelot et al., 2021]. The Gaussian measurement error model is

$$ Y_k = LU_k + \varepsilon_k, \quad \varepsilon_k \sim \mathcal{N}(0, \sigma_{obs}^2 I_2), \quad (3.1) $$

where $L = (I_2 \quad 0_{2,2})$ and $(U_k)_{1 \le k \le n} = (X_k \quad V_k)_{1 \le k \le n}$ is the vector of true positions and velocities at observation times. Here $I_2$ denotes the 2-dimensional identity matrix and $0_{2,2}$ denotes a 2 × 2 null-matrix. The parameter $\sigma_{obs}$ quantifies the scale of the errors in the directions Northing and Easting (isotropic).
However, GPS tracks often exhibit heavy-tailed error distributions [Wensveen et al., 2015], which are best modelled with a location-scale Student's t-distribution, $t_d$, where $d$ is the degrees of freedom. The measurement error model is then

$$Y_k = LU_k + \sigma_{obs}\varepsilon_k, \quad \varepsilon_k \sim t_d, \qquad (3.2)$$

where $\sigma_{obs}$ quantifies the isotropic error scale in Northing and Easting. In practice, both $\sigma_{obs}$ and $d$ may vary with the number of satellites that capture the signal [Wensveen et al., 2015]. An increased number of satellites results in more accurate positions with a distribution closer to a Gaussian (i.e., larger $d$). For simplicity, we consider a constant value of $\sigma_{obs}$ and $d$ in this study, however, it can easily be generalized.

For tracks relying on the Argos system, Doppler shift is used to estimate animal positions and the errors have generally higher variance and a complex shape favouring errors in the North-West/South-East and North-East/South-West directions [Brost et al., 2015; Hanks et al., 2017]. This can be modeled as follows:

$$Y_k = LU_k + \varepsilon_k \qquad (3.3)$$

where

$$\varepsilon_k \sim \begin{cases} t_d(0, \Sigma) & \text{with probability } p, \\ t_d(0, \tilde{\Sigma}) & \text{with probability } 1-p, \end{cases}$$

with $\Sigma = \sigma_{obs}^2 \begin{pmatrix} 1 & \rho\sqrt{a} \\ \rho\sqrt{a} & 1 \end{pmatrix}$ and $\tilde{\Sigma} = \sigma_{obs}^2 \begin{pmatrix} 1 & -\rho\sqrt{a} \\ -\rho\sqrt{a} & 1 \end{pmatrix}$. Here, $t_d(0, \Sigma)$ denotes a bivariate Student's t-distribution with $d$ degrees of freedom, mean vector 0 and scale matrix $\Sigma$. The parameter $p$ assigns weights to the North-West/South-East and North-East/South-West directions, $\rho \in (-1, 1)$ is a correlation parameter and $a > 0$ is an anisotropy parameter. When $a = 0$, we obtain an isotropic Student distribution similar to (3.2). Obviously, only the product $\rho\sqrt{a}$ can be identified, not each parameter instead.

Figure 1 shows a trajectory of (2.1) with observations simulated with Gaussian, Student's t, and Argos error distributions. In the Gaussian case, the observed positions are distributed circularly around the true position, reflecting small and isotropic errors. For the case of Student's t, the majority of the observations remain close to the true path, but occasional outliers appear far from the animal's actual position, illustrating the heavy-tailed nature of the distribution. In the Argos case, the errors are both large and directional: observations spread preferentially along the North-West/South-East and North-East/South-West axes, producing stretched error patterns that differ strongly from the other two models.

# 4 Splitting schemes

To efficiently simulate the non-linear dynamics of the penalized Langevin SDE in eq. (2.1), we need numerical approximation methods. With no penalization and a flat potential surface, the SDE is linear with an explicit Gaussian solution. However, both the penalization term and the potential function $H$ introduce non-linearity into the equation and the solution is not explicit. We need a numerical scheme to approximate the solution. Classical schemes such as Euler-Maruyama may not be precise enough to capture the non-linearity. We instead employ a splitting method to approximate the solution of the penalized SDE [Pilipovic et al.,
<figure><img src="image_2.png" /><figcaption>Figure 1: Part of a trajectory simulated from model (2.1) with observations from three error distributions. The black line is the true trajectory, with time step of 5 minutes between consecutive positions. The grey thick line is one side of the polygon defining the domain $\mathcal{D}$. Orange points correspond to observations with Gaussian errors ($\sigma_{obs} = 0.2$) with line segments connecting true and observed positions. Blue points illustrate observations with heavy-tailed Student's t errors with scale parameter $\sigma_{obs} = 0.2$ and degree of freedom $d = 3$. Red points show observations obtained from a mixture of multivariate Student's t-distributions with $\sigma_{obs} = 0.2$, $\rho = 0.7$, $a = 0.9$, $p = 0.5$, and degree of freedom $d = 3$, representing typical Argos data errors.</figcaption></figure>

2024], which is particularly suitable for non-linear SDEs. The idea is to split the process into a linear SDE, for which we know the distribution, and an ordinary differential equation (ODE), which can be (explicitly) solved.

In practice, it is recommended to split the drift in the SDE around one of its fixed points [Pilipovic et al., 2025]. In the interior of the domain $\mathcal{D}$, the penalization term is zero and we can choose a fixed point $u_l^* = (x_l^* \ 0_{2,1})^\top$, $l \in \{1, \dots, J\}$, where $x_l^*$ is one of the centers of the Gaussians in the potential $H$. Hence, we rewrite the SDE as

$$d \begin{pmatrix} X(t) \\ V(t) \end{pmatrix} = \begin{pmatrix} 0_{2,2} & I_2 \\ -2\alpha_l B_l & -A \end{pmatrix} \left[ \begin{pmatrix} X(t) \\ V(t) \end{pmatrix} - \begin{pmatrix} x_l^* \\ 0_{2,1} \end{pmatrix} \right] dt + \begin{pmatrix} 0_{2,2} \\ \sigma I_2 \end{pmatrix} dW(t) + \begin{pmatrix} 0_{2,1} \\ -2\alpha_l(e_l(X(t)) - 1)B_l(X(t) - x_l^*) - \nabla H_{-l}(X(t)) - \beta_\lambda(X(t)) \end{pmatrix} dt \quad (4.1)$$

where $H_{-l}(x) = - \sum_{j=1, j \neq l}^J H_j(x)$. We then define the splitting through the two processes:
$$
\left\{
\begin{array}{l}
d \begin{pmatrix} X^{(1)}(t) \\ V^{(1)}(t) \end{pmatrix} = \begin{pmatrix} 0_{2,2} & I_2 \\ -2\alpha_l B_l & -A \end{pmatrix} \left[ \begin{pmatrix} X^{(1)}(t) \\ V^{(1)}(t) \end{pmatrix} - \begin{pmatrix} x_l^* \\ 0_{2,1} \end{pmatrix} \right] dt + \begin{pmatrix} 0_{2,2} \\ \sigma I_2 \end{pmatrix} dW(t) \\
d \begin{pmatrix} X^{(2)}(t) \\ V^{(2)}(t) \end{pmatrix} = \begin{pmatrix} 0_{2,1} \\ -2\alpha_l(e_l(X^{(2)}(t)) - 1)B_l(X^{(2)}(t) - x_l^*) - \nabla H_{-l}(X^{(2)}(t)) - \beta_\lambda(X^{(2)}(t)) \end{pmatrix} dt
\end{array}
\right.
\quad (4.2)
$$

The splitting approximation combines the solutions of the two equations. There are several ways to combine these solutions that each give different approximations. Denote the solutions of the SDE and ODE with initial condition $u = (x \ v)^\top \in \mathbb{R}^4$ by $\phi_h^{(1)}(u)$ and $\phi_h^{(2)}(u)$, respectively, over a time-interval of length $h$. The Lie-Trotter splitting provides the solution $U^{LT}$ iteratively defined by:

$$U^{LT}(t+h) = \phi_h^{(1)} \circ \phi_h^{(2)}(U^{LT}(t)).$$

The Strang splitting provides the solution $U^S$ defined iteratively by:

$$U^S(t+h) = \phi_{\frac{h}{2}}^{(2)} \circ \phi_h^{(1)} \circ \phi_{\frac{h}{2}}^{(2)}(U^S(t)).$$

In the following, we detail the solutions $\phi_h^{(1)}$ and $\phi_h^{(2)}$ of the two systems, see details in Appendix B.

The first process in (4.2) is a multidimensional Ornstein-Uhlenbeck process with explicit solution

$$\begin{pmatrix} X^{(1)}(t+h) \\ V^{(1)}(t+h) \end{pmatrix} = e^{\tilde{A}h} \left( \begin{pmatrix} X^{(1)}(t) \\ V^{(1)}(t) \end{pmatrix} - \begin{pmatrix} x_l^* \\ 0_{2,1} \end{pmatrix} \right) + \begin{pmatrix} x_l^* \\ 0_{2,1} \end{pmatrix} + \eta(h)$$

with $\tilde{A} = \begin{pmatrix} 0_{2,2} & I_2 \\ -2\alpha_l B_l & -A \end{pmatrix}$, $\eta(h) \sim \mathcal{N}(0, \tilde{Q}(h))$ and

$$\tilde{Q}(h) = \int_{0}^{h} e^{-\tilde{A}(u-h)} \Gamma e^{-\tilde{A}^{\top}(u-h)} du; \quad \Gamma = \begin{pmatrix} 0_{2,2} & 0_{2,2} \\ 0_{2,2} & \sigma^2 I_2 \end{pmatrix}.$$

The covariance matrix $\tilde{Q}(h)$ is defined as an integral which can be computed it explicitly. Writing $C$ for the matrix such that $\mathrm{vec}(C) = (\tilde{A} \oplus \tilde{A})^{-1}\mathrm{vec}(\Gamma)$, we have, see details in Appendix C,

$$\tilde{Q}(h) = e^{\tilde{A}h} C e^{\tilde{A}^{\top}h} - C.$$

The second process in (4.2) is an ODE with solution $\phi_h^{(2)}(u) = u - hg_\lambda(u)$ where

$$g_{\lambda}(u) = g_{\lambda}(x, v) = \begin{pmatrix} 0_{2,1} \\ 2\alpha_l(e_l(x) - 1)B_l(x - x_l^*) + \nabla H_{-l}(x) + \beta_{\lambda}(x) \end{pmatrix}.$$

Note that the function $g_\lambda$ depends only on $x$ and not on $v$.

The Lie-Trotter approximation at time $t+h$, given the current value $U^{LT}(t)$ at time $t$ is

$$U^{LT}(t+h) = e^{\tilde{A}h} (U^{LT}(t) - u_l^* - hg_\lambda(U^{LT}(t))) + u_l^* + \eta(h).$$
The Strang approximation is defined by

$$ U^S(t+h) = Z - \frac{h}{2}g_{\lambda}(Z), \text{ with } Z = e^{\tilde{\lambda}h}(\hat{U}_{\frac{h}{2}} - u_l^*) + u_l^* + \eta(h) \text{ and } \hat{U}_{\frac{h}{2}} = U^S(t) - \frac{h}{2}g_{\lambda}(U^S(t)). \quad (4.3) $$

Conditionally on $U^S(t)$, $Z := (X_Z \ V_Z)^\top$ is Gaussian with mean $\tilde{m}(h) = e^{\tilde{\lambda}h}(\hat{U}_{\frac{h}{2}} - u_l^*) + u_l^*$ and covariance matrix $\tilde{Q}(h)$ [Pilipovic et al., 2024]. Since the first two coordinates of $g_\lambda$ are zero, we have that $X^S(t+h) = X_Z$, implying that the marginal distribution of $X^S(t+h)$ conditionally on $U^S(t)$ is Gaussian with $\tilde{m}_x$ and covariance $\tilde{Q}_{xx}$, where we denote $\tilde{m}(h) = (\tilde{m}_x(h) \ \tilde{m}_v(h))$ and $\tilde{Q}(h) = \begin{pmatrix} \tilde{Q}_{xx} & \tilde{Q}_{xv} \\ \tilde{Q}_{vx} & \tilde{Q}_{vv} \end{pmatrix}$. Moreover,

$$ V^S(t+h) = V_Z - \frac{h}{2}2\alpha_l(e_l(X_Z) - 1)B_l(X_Z - x_l^*) + \nabla H_{-l}(X_Z) + \beta_\lambda(X_Z). $$

Therefore, conditionally on $U^S(t)$ and $X^S(t+h) = x$, $V^S(t+h)$ is also Gaussian with mean $\tilde{m}_{v|x}$ and covariance matrix $\tilde{Q}_{v|x}$ given by

$$ \tilde{m}_{v|x} = \tilde{m}_v + \tilde{Q}_{vx}\tilde{Q}_{xx}^{-1}(x - \tilde{m}_x) - \frac{h}{2}(2\alpha_l(e_l(x) - 1)B_l(x - x_l^*) + \beta_\lambda(x) + \nabla_{-l}H(x)), $$

$$ \tilde{Q}_{v|x} = \tilde{Q}_{vv} - \tilde{Q}_{vx}\tilde{Q}_{xx}^{-1}\tilde{Q}_{xv}. $$

In practice, we choose the fixed point $x_l^*$ depending on the current location $X(t)$. We choose the center as $l = \underset{j \in \{1, \dots, J\}}{\operatorname{argmax}} \log(||\nabla H_j(X(t))||^2)$ where

$$ \log(||\nabla H_j(X(t))||^2) = 2\log(2) + 2\log(\alpha_j) - 2(x - x_j^*)^\top B_j(x - x_j^*) + \log((x - x_j^*)^\top B_j^\top B_j(x - x_j^*)). $$

Figure 2 shows a trajectory of (2.1) with a mixture of two Gaussian potentials around $x_1^* = (25, 5)$ and $x_2^* = (35, 15)$ obtained with the Lie-Trotter approximation with time step $h = 1/3600$ hour (1 second between consecutive observations) and penalization $\lambda = h^{0.8}$ over the time interval $[0, 72]$. The trajectory has been subsampled to 1 minute between observation times. The parameters of the SDE are $\tau = 1$, $\nu = 5$ and $\omega = 0.1$. The weights in the potential are $\alpha_1 = 70$ and $\alpha_2 = 50$. The precision matrices $B_1$ and $B_2$ are chosen as

$$ B_1 = \begin{pmatrix} \frac{1}{9} & \frac{1}{40} \\ \frac{1}{40} & \frac{1}{4} \end{pmatrix}; \quad B_2 = \begin{pmatrix} \frac{1}{36} & -\frac{1}{100} \\ -\frac{1}{100} & \frac{1}{100} \end{pmatrix}. \quad (4.4) $$

On average, the animal spends more time near the potential centers $x_1^*$ and $x_2^*$, which represent the preferred habitat. The movement is persistent in time due to the autocorrelation parameter $\tau > 0$, and areas outside the domain $\mathcal{D}$ are very rarely visited as the penalization reorients the velocity inward each time the position falls outside the boundaries of $\mathcal{D}$. The model thus captures key features of wildlife movement: directional persistence, attraction toward biologically meaningful sites such as feeding or mating areas, and confinement to the appropriate habitat. This combination produces a realistic and interpretable movement model consistent with many ecological applications.
<figure><img src="image_3.png" /><figcaption>Figure 2: Trajectory from (2.1) with a mixture of Gaussian potentials using the Lie-Trotter approximation scheme. Parameters are $h = 1/3600$, $\lambda = h^{0.8}$, $\tau = 1$, $\nu = 5$, $\omega = 0.1$.</figcaption></figure>

# 5 Filtering algorithms

In this section, we present algorithms for estimating the latent states of the movement process from noisy observations. We first consider the classical Extended Kalman Filter for Gaussian measurement errors and then generalize to particle filtering approaches capable of handling heavy-tailed or complex error distributions such as Student's t or Argos X-shaped errors.

## 5.1 Extended Kalman filter for Gaussian measurement errors

We assume that the true positions $X_k = X(kh)$, $k \in \{1, \dots, n\}$ are observed with Gaussian measurement noise. We use the Lie-Trotter approximation of the state process to obtain a Gaussian state-space model for subsequent Kalman filtering. The non-linear Gaussian state-space is then

$$ \left\{ \begin{aligned} U_{k+1}^{LT} &= e^{\tilde{\lambda}h} U_k^{LT} - e^{\tilde{\lambda}h} u_l^* - he^{\tilde{\lambda}h} g_{\lambda}(U_k^{LT})) + u_l^* + \eta_k; \quad \eta_k \sim \mathcal{N}(0, \tilde{Q}) \\ Y_k &= LU_k^{LT} + \varepsilon_k; \quad \varepsilon_k \sim \mathcal{N}(0, \sigma_{\text{obs}}^2 I_2) \end{aligned} \right. \qquad (5.1) $$

where $L = (I_2, I_2)$ and the other matrices are as before. The subscript $h$ is suppressed for convenience. The non-linearity comes from the potential $H$ as well as the penalization $\beta_{\lambda}$ encoded in the function $g_{\lambda}$. We now detail the Extended Kalman Filter algorithm that produces state estimates $\hat{U}_k^{LT}$, $k \in \{1, \dots, n\}$ approximating both the position and the velocity given the observations.
First, initialize $\hat{U}_0$ and $R_0$. Then iterate as follows:

$$ \hat{U}_{k-}^{LT} = e^{\bar{A}h} \hat{U}_{k-1}^{LT} - e^{\bar{A}h} u_l^* - he^{\bar{A}h} g_{\lambda}(\hat{U}_{k-1}^{LT}) + u_l^*, \quad (5.2) $$

$$ R_{k-}^{LT} = e^{\bar{A}h} R_{k-1}^{LT} e^{\bar{A}^{\top}h} + h^2 e^{\bar{A}h} G_{k-1} R_{k-1}^{LT} G_{k-1}^{\top} e^{\bar{A}^{\top}h} + \tilde{Q}, \quad (5.3) $$

where $G_{k-1} = Dg_{\lambda}(\hat{U}_{k-1}^{LT}) = \left( \frac{g_{\lambda}^{(i)}}{\partial u_j} (\hat{U}_{k-1}^{LT}) \right)_{1 \le i, j \le 4}$. Let $g_{v,\lambda}$ denote the velocity components for the function $g_{v,\lambda}$, that is

$$ g_{v,\lambda}(x, v) = 2\alpha_l(e_l(x) - 1)B_l(x - x_l^*) + \nabla H_{-l}(x) + \beta_{\lambda}(x). $$

Here,

$$ G_{k-1} = \begin{pmatrix} 0_2 & 0_2 \\ D_x g_{v,\lambda}(\hat{U}_{k-1}^{LT}) & 0_2 \end{pmatrix} $$

where $D_x$ is the Jacobian matrix with respect to the $x$ components.

The Hessian of the potential surface (2.4) is given by:

$$ D_x g_{v,\lambda}(x) = \frac{I_2 - D\pi(\hat{X}_{k-1}^{LT})}{\lambda} + D^2 H_{-l}(x) + 2\alpha_l \left( (e_l(x) - 1)B_l - 2e_l(x)B_l(x - x_l^*)(x - x_l^*)^\top B_l \right), $$

where $D\pi$ is the Jacobian matrix of the projection operator. Recall that the projection $\pi$ is piecewise linear in $x$. The Jacobian is thus

$$ D\pi(x) = \frac{(v_{k+1} - v_k)(v_{k+1} - v_k)^\top}{l^2}. $$

Then, the correction step of the extended Kalman filter is

$$ K_k^{LT} = R_{k-}^{LT} L^\top (L R_{k-}^{LT} L^\top + \sigma_{obs}^2 I_2)^{-1}, \quad (5.4) $$

$$ \hat{U}_k^{LT} = \hat{U}_{k-}^{LT} + K_k^{LT} (Y_k - L \hat{U}_{k-}^{LT}), \quad (5.5) $$

$$ R_k^{LT} = (I_4 - K_k^{LT} L) R_{k-}^{LT}. \quad (5.6) $$

## 5.2 Particle filter for Student's t-distribution measurement errors

Observed animal positions often suffer from non-Gaussian measurement errors, for example for Fastloc-GPS observations [Wensveen et al., 2015], which are better described by a Student's t-distribution. Then Kalman filtering techniques may be unreliable, especially if the error distribution differs a lot from a Gaussian distribution. To address this, we can resort to particle filtering techniques, also called Sequential Monte Carlo (SMC) [Doucet and Johansen, 2011]. We develop two particle filtering algorithms depending on the splitting scheme (Lie-Trotter or Strang), which is used to approximate the solution to the SDE.

### 5.2.1 Particle filter based on Lie-Trotter scheme

The state-space updates derived from the Lie-Trotter approximation of the SDE are:

$$ \begin{cases} U_{k+1}^{LT} = e^{\bar{A}h} U_k^{LT} - e^{\bar{A}h} u_l^* - he^{\bar{A}h} g_{\lambda}(U_k^{LT}) + u_l^* + \eta_k; & \eta_k \sim \mathcal{N}(0, \tilde{Q}), \\ Y_k = L U_k^{LT} + \sigma_{obs} \varepsilon_k; & \varepsilon_k \sim t_d. \end{cases} $$
The Student's t-distribution can be approximated by a Gaussian distribution with standard deviation $\sigma'_{obs} = \frac{d-2}{d}\sigma_{obs}$ to transform the equations into a Gaussian state-space model and then resort to Kalman filtering algorithm.

However, without this simplification, the filtering distributions $p(U_k|Y_{1:k})$ are not explicit and a particle filter is needed to approximate it. Particle filtering is described in Algorithm 1.

**Algorithm 1 SMC algorithm**

**Require:** Data y

Sample $U_0^{(k)} \sim \pi_0()$ for $k \in \{1, \dots, K\}$.

Compute and normalize the weights

$$ \omega_0(U_0^{(k)}) = p(y_0, U_0^{(k)}); \quad W_0(U_0^{(k)}) = \frac{w_0(U_0^{(k)})}{\sum_{k=1}^{K} \omega_0(U_0^{(k)})} $$

**for** $j = 1, \dots, J$ **do**

**for** $k = 1, \dots, K$ **do**

Sample $U_{0:j-1}^{(k)}$ from the distribution $\Psi_{j-1}^K = \sum_{k=1}^K W_{j-1}^{(k)} \delta_{U_{0:j-1}^{(k)}}$ for $k \in \{1, \dots, K\}$.

Sample $U_j^{(k)} \sim q(\cdot|y_j, U_{0:j-1}^{(k)})$.

Set $U_{0:j}^{(k)} = (U_{0:j-1}^{(k)} \quad U_j^{(k)})$

Compute and normalize the weights

$$ w_j^{(k)} = \frac{p(y_j|U_j^{(k)})p(U_j^{(k)}|U_{j-1}^{(k)})}{q(U_j^{(k)}|y_j, U_{0:j-1}^{(k)})}; \quad W_j^{(k)} = \frac{w_j^{(k)}}{\sum_{k=1}^{K} w_j^{(k)}} $$

**end for**

**end for**

After running the algorithm, the conditional mean $\hat{E}(X_j|Y_{1:j-1})$ is approximated by

$$ \hat{E}(X_j|Y_{1:j}) = \sum_{k=1}^{K} w_j^k X_j^{(k)}. $$

A proposal distribution $q$ has to be chosen in the particle filter. The conditional distribution $p(U_j|U_{j-1}, y_j)$ is known to be the optimal proposal, that is, it minimizes the variance of the particle weights [Doucet and Johansen, 2011]. However, this density cannot be computed explicitly in our model. We propose to approximate it as follows. Layes's formula gives

$$ p(U_j|U_{j-1}, y_j) \propto p(U_j|U_{j-1})p(y_j|U_j). $$

We replace $p(U_j|U_{j-1})$ in this equation by an approximation $p^{LT}(U_j|U_{j-1})$ obtained from the Lie-Trotter splitting scheme. We also propose to approximate $p(y_j|U_j)$ in the proposal by a Gaussian density $\tilde{p}(y_j|U_j)$ with mean $LU_j$ and covariance matrix $\sigma_{obs}^2 \frac{d}{d-2}I_2$. Hence, the new proposal $q(U_j|U_{j-1}, y_j) \propto p^{LT}(U_j|U_{j-1})\tilde{p}(y_j|U_j)$ is a product of Gaussian densities. It is therefore a Gaussian proposal with covariance matrix and mean given by

$$ \Gamma^{LT} = \left( \tilde{Q}^{-1} + \frac{1}{\sigma_{obs}^2} \frac{d-2}{d} L^\top L \right)^{-1}, $$
$$m_j^{LT} = \Gamma^{LT} \left( \tilde{Q}^{-1} (e^{\tilde{A}h}(U_{j-1} - u_l^* - hg_\lambda(U_{j-1})) + u_l^*) + \frac{1}{\sigma_{obs}^2} \frac{d-2}{d} L^\top y_j \right).$$

It is then easy to propagate the particles using the Gaussian proposal. The weights are also easily computed as follows:

$$w_j^{(k)} = \frac{p(y_j|U_j^{(k)})p^{LT}(U_j^{(k)}|U_{j-1}^{(k)})}{\mathcal{N}(U_j^{(k)}; m_j^{LT}, \Gamma^{LT})}.$$

In this algorithm, the penalization plays a role in the propagation of the particles as well as in the weighting of each particle through the Gaussian kernel $q(U_j|U_{j-1}, y_j)$.

*Remark.* An alternative to modifying the underlying dynamics by penalizing the SDE is to include constraints by setting the weights to zero when the propagated particle does not belong to the domain. This modify the weight equations as follows,

$$w_j^{(k)} = \mathbb{1}_D(X_j^{(k)}) \frac{p(y_j|U_j^{(k)})p(U_j^{(k)}|U_{j-1}^{(k)})}{q(U_j^{(k)}|y_j, U_{0:j-1}^{(k)})}.$$

However, this may discard too many particles due to the hard decision threshold and create weight degeneracy when a corner of the domain is reached. Moreover, it would require precise knowledge of the domain boundary (for instance, shoreline maps) not to drop particles that actually are in the domain. With our method, the weights are smooth as a function of the position of the propagated particle (up to resampling), and the penalization allows for unprecise knowledge of the boundary through the parameter $\lambda$.

### 5.2.2 Particle filter based on the Strang approximation

We now introduce a particle filter based on the Strang splitting scheme,

$$U^S(t+h) = (X_Z, V_Z - \frac{h}{2}(\beta_\lambda(X_Z) + \nabla H(X_Z)))$$

where $Z = (X_Z, V_Z)$ is given in (4.3).

We apply Algorithm 1 and we need to specify the proposal $q$. The optimal proposal in the particle filter can be rewritten as

$$p(U_j|U_{j-1}, y_j) \propto p(X_j|U_{j-1})p(V_j|X_j, U_{j-1})p(y_j|U_j).$$

We define the proposal distribution as

$$q(U_j|U_{j-1}, y_j) \propto p^S(V_j|X_j, U_{j-1})p^S(X_j|U_{j-1})\tilde{p}(y_j|U_j) \quad (5.7)$$

where $\tilde{p}(y_j|U_j)$ is a Gaussian density with mean $LU_j$ and covariance matrix $\sigma_{obs}^2 \frac{d}{d-2} I_2$. Hence, the product $p^S(X_j|U_{j-1})\tilde{p}(y_j|U_j)$ is Gaussian with covariance matrix and mean

$$\Gamma^S = (\tilde{Q}_{xx}^{-1} + \frac{1}{\sigma_{obs}^2} \frac{d-2}{d} I_2)^{-1}; \quad \hat{x}_j^S = \Gamma^S (\tilde{Q}_{xx}^{-1} \tilde{m}_x + \frac{1}{\sigma_{obs}^2} \frac{d-2}{d} y_j).$$
Note that $\tilde{p}(y_j|U_j) = \tilde{p}(y_j|X_j)$ so that the normalisation constant in (5.7) is given by
$\int p(X_j|U_{j-1})\tilde{p}(y_j|X_j)dX_j$ and we obtain the proposal

$$q(U_j|U_{j-1}, y_j) = p^S(V_j|X_j, U_{j-1})\mathcal{N}(X_j; \hat{x}_j^S, \Gamma^S).$$

The propagation of the particles is thus done in two steps. First, we propagate the position based on the current position and velocity with a Gaussian kernel. Then we propagate the velocity given the new position. As mentioned in section 4, $p^S(V_j|X_j, U_{j-1})$ is explicit and Gaussian. This produces a simple formula for the updates of weights

$$w_j^{(k)} = \frac{p(y_j|U_j^{(k)})p^S(X_j^{(k)}|U_{j-1}^{(k)})}{\mathcal{N}(X_j^{(k)}; \hat{x}_j^S, \Gamma^S)}.$$

## 5.3 Particle filter for X-shaped measurement errors

For Argos telemetry data, the error distribution has a complex X-shape that can be modelled with a mixture of multivariate Student's t-distributions [Brost et al., 2015] detailed in section 3. We present a particle filter algorithm adapted to such data, either based on the Lie-Trotter or the Strang approximation.

### 5.3.1 Particle filter based on the Lie-Trotter approximation

For the proposal in the SMC algorithm, we also approximate the Bayes formula defining the posterior distribution. We consider a mixture $\tilde{p}(y_j|U_j)$ of Gaussian densities with covariance matrices $\Sigma$ and $\tilde{\Sigma}$ to approximate the term $p(y_j|U_j)$. For Lie-Trotter, $p^{LT}(U_j|U_{j-1})\tilde{p}(y_j|U_j)$ is a product of Gaussian densities. Hence, the proposal $q$ is a mixture of two multivariate Gaussians, which is easy to simulate. The covariance matrices are

$$\Gamma^{LT} = (\tilde{Q}^{-1} + \frac{d-2}{d} L^{\top} \Sigma L)^{-1}; \quad \tilde{\Gamma}^{LT} = (\tilde{Q}^{-1} + \frac{d-2}{d} L^{\top} \tilde{\Sigma} L)^{-1}$$

and the means are

$$m_j^{LT} = \Gamma^{LT}(\tilde{Q}^{-1}(e^{\tilde{A}h}(U_{j-1} - u_l^* - hg_\lambda(U_{j-1})) + u_l^*) + \frac{d-2}{d}L^{\top}\Sigma^{-1}y), \quad (5.8)$$

$$\tilde{m}_j^{LT} = \tilde{\Gamma}^{LT}(\tilde{Q}^{-1}(e^{\tilde{A}h}(U_{j-1} - u_l^* - hg_\lambda(U_{j-1})) + u_l^*) + \frac{d-2}{d}L^{\top}\tilde{\Sigma}^{-1}y). \quad (5.9)$$

The update of the weights is then:

$$w_j^{(k)} = \frac{p(y_j|U_j^{(k)})p^{LT}(U_j^{(k)}|U_{j-1}^{(k)})}{p\mathcal{N}(U_j^{(k)}; m_j^{LT}T, \Gamma^{LT}) + (1-p)\mathcal{N}(U_j^{(k)}; \tilde{m}_j^{LT}, \tilde{\Gamma}^{LT})}.$$

### 5.3.2 Particle filter based on the Strang approximation

Similar to the previous method with Student measurement errors, we approximate $p(y_j|U_j)$ in the proposal by a mixture of Gaussian densities with covariance matrices $\Sigma, \tilde{\Sigma}$. The product $p^S(X_j|U_{j-1})\tilde{p}(y_j|U_j)$ is thus a mixture of Gaussian densities with covariance matrices

$$\Gamma^S = (Q_{xx}^{-1} + \frac{d-2}{d}\Sigma^{-1})^{-1}; \quad \tilde{\Gamma}^S = (Q_{xx}^{-1} + \frac{d-2}{d}\tilde{\Sigma}^{-1})^{-1}$$
and means

$$ \hat{x}_j^S = \Gamma^S(\tilde{Q}_{xx}^{-1}\tilde{m}_x + \frac{d-2}{d}\Sigma^{-1}y_j) \quad \tilde{x}_j^S = \tilde{\Gamma}^S(\tilde{Q}_{xx}^{-1}\tilde{m}_x + \frac{d-2}{d}\tilde{\Sigma}^{-1}y_j). $$

The update of the weights is

$$ w_j^{(k)} = \frac{p(y_j|U_j^{(k)})p^S(X_j^{(k)}|U_{j-1}^{(k)})}{p\mathcal{N}(X_j^{(k)}; \hat{x}_j^S, \Gamma^S) + (1-p)\mathcal{N}(X_j^{(k)}; \tilde{x}_j^S, \tilde{\Gamma}^S)}. $$

# 6 Simulation study

We start by evaluating the filtering methods in a standard scenario where measurement errors are Gaussian. This setting allows us to isolate the effects of domain penalization and sampling frequency on the filtering accuracy. In particular, we study how the different filters perform when trajectories are simulated within a complex polygonal domain and corrupted by Gaussian noise. We then assess the performance of particle filters for the non-Gaussian, heavy-tailed error models introduced in Section 3, and compare their accuracy with that of the Kalman and Extended Kalman filters, in which the non-Gaussian errors are approximated by Gaussian errors with matching variances.

## 6.1 Gaussian errors

The domain $\mathcal{D}$ is a polygon with 19 vertices, see Figure 2. The values of the movement parameters and the potential parameters are given at the end of Section 4.

We simulate 50 trajectories over a period $[0, T]$ with $T = 72$ hours and constant time step $h = 1$ second based on the Lie-Trotter approximation. The penalization is fixed to $\lambda = h^{0.8}$. We add Gaussian measurement errors with standard deviation $\sigma_{obs} = 0.2$ to create noisy observations of the trajectories.

The data is then subsampled to time steps 1 min, 3 min, 5 min and 20 min. We test the filtering algorithms for the four sampling intervals.

We compare the results of four methods: 1) Naive Kalman filtering without domain penalization ($\lambda = +\infty$) where the non-linear part in the state space is treated as a constant term (KF); 2) Extended Kalman filtering without domain penalization ($\lambda = +\infty$) (EKF); 3) Naive Kalman filtering with domain penalization (Penalized KF) and 4) Extended Kalman filtering with domain penalization (Penalized EKF).

The results are shown in Figure 3. Figure 3a shows boxplots of the root mean squared error (RMSE) for each value of the time step $h$ using KF, EKF, Penalized KF and Penalized EKF. The theoretical RMSE before filtering is equal to $\sqrt{\frac{\pi}{2}}\sigma_{obs} \approx 0.25$. The incorporation of spatial constraints through the domain penalization reduces the RMSE for small time steps. For a time step of $h = 20$ minutes, the Kalman filters with or without domain penalization have slightly poorer performance than the extended Kalman filters. The reason is that the approximation of constant potential over each time step that is used in the naive Kalman filter is very poor for large time steps. The extended Kalman filter then provides a slightly better approximation of the potential surface and hence slightly better accuracy.

Figure 3b shows boxplots of the maximum error between the true positions and the filtered positions on the 50 trajectories for each value of the time step $h$ with the four filters. The
incorporation of domain knowledge improves the performance of the filter for high frequency observations ($h = 1, 3$ or $5$ min). Again, the Kalman filter has slightly poorer performance than the Extended Kalman filter for large time steps ($h = 20$ min). Mean values of the maximum error and RMSEs over the 50 trajectories are shown in Table 1. The simulation study was performed on a personal laptop and took only 35 minutes to run.

Overall, the penalized filters avoid large errors in the reconstructed trajectories and ensure that they remain within or close enough to the spatial domain of interest, which is critical in applications to animal movements. The benefits of domain penalization are most pronounced for small time steps (high-frequency observations), where the filter has enough temporal resolution to adjust to the local shape of the domain.

<table><caption>Table 1: Max error and root mean squared error (RMSE) averaged over the 50 trajectories with Gaussian measurement error for different values of the time step $h$. The smallest value in each row is marked in bold.</caption><thead><tr><th>$h$</th><th>KF</th><th>EKF</th><th>Penalized KF</th><th>Penalized EKF</th><th>Before filter</th></tr></thead><tbody><tr><td colspan="6" style="text-align: center;">Average max error (km)</td></tr><tr><td>1 min</td><td>0.540</td><td>0.540</td><td><strong>0.437</strong></td><td><strong>0.437</strong></td><td>0.851</td></tr><tr><td>3 min</td><td>0.643</td><td>0.643</td><td><strong>0.561</strong></td><td>0.562</td><td>0.792</td></tr><tr><td>5 min</td><td>0.656</td><td>0.656</td><td>0.657</td><td><strong>0.626</strong></td><td>0.761</td></tr><tr><td>20 min</td><td>0.694</td><td>0.682</td><td>1.087</td><td><strong>0.680</strong></td><td>0.686</td></tr><tr><td colspan="6" style="text-align: center;">Average RMSE (km)</td></tr><tr><td>1 min</td><td>0.136</td><td>0.136</td><td><strong>0.133</strong></td><td><strong>0.133</strong></td><td>0.251</td></tr><tr><td>3 min</td><td>0.185</td><td>0.185</td><td>0.182</td><td><strong>0.181</strong></td><td>0.250</td></tr><tr><td>5 min</td><td>0.207</td><td>0.207</td><td>0.205</td><td><strong>0.204</strong></td><td>0.251</td></tr><tr><td>20 min</td><td>0.249</td><td><strong>0.247</strong></td><td>0.256</td><td><strong>0.247</strong></td><td>0.250</td></tr></tbody></table>

## 6.2 Student's t-distributed errors

We consider the same model and domain $\mathcal{D}$. We keep the same values for parameters $\tau$, $\nu$, and $\omega$. We simulate 20 trajectories over a period $[0, T]$ with $T = 48$ based on the Strang approximation. The penalization parameter is again set to $\lambda = h^{0.8}$, where $h = \frac{1}{3600}$ denotes the simulation time step. The data is subsampled as previously. We add $t_d$ measurement errors with scale parameter $\sigma_{obs} = 0.2$ and degrees of freedom $d = 3$ to create noisy observations of the trajectories. We run an EKF and a Penalized EKF where the $t_d$ is approximated by a Gaussian distribution with variance $\frac{d}{d-2}\sigma_{obs}^2$. Then, we run a particle filter without domain penalization based on the Lie-Trotter scheme (Lie-Trotter PF), a particle filter with domain penalization based on the Lie-Trotter scheme (Penalized Lie-Trotter PF), a particle filter without domain penalization based on the Strang scheme (Strang PF) and a particle filter with domain penalization based on the Strang scheme (Penalized Strang PF). We choose $K = 1000$ particles in all cases. The calculations are performed in parallel over 48 CPU cores. The particle filters are notably longer to run than the Kalman filters. The whole simulation took about 15 hours.

Figure 4a shows boxplots of the RMSE for each value of the time step $h$. The incorporation of domain knowledge associated with a particle filter to handle non-Gaussian error increases
<figure><img src="image_4.png" /><figcaption>Figure 3: Results of filtering for 50 simulated trajectories with additive Gaussian noise. (a) Boxplot of RMSE for KF, EKF, Penalized KF and Penalized EKF. (b) Boxplot of max errors.</figcaption></figure>

the performance of the filter for time steps $h = 1, 3$ and $5$ min. Notably, the RMSE for the particle filters without domain penalization is higher than the RMSE obtained with EKF. The misspecification of the underlying dynamics in the particle filter likely produces trajectories that progressively drift away from the true track each time we linger outside the boundary. Figure 4b shows the same boxplots for the maximum error. The incorporation of domain
knowledge also increases the performance of the filter for time steps $h = 1, 3$ and $5$ min. For $h = 1$ min, the maximum error is reduced by about 40% for the penalized particle filters compared to the EKF and penalized EKF. For large time steps ($h = 20$ min), the penalized filters give worse results than the non-penalized ones. The explanation is as follows. When the time step is short, the term $\beta_{\lambda}$ in the SDE is computed frequently enough to adapt to the boundary and decreases as we approach the boundary. When the time step becomes too large, the term $\beta_{\lambda}$ can become dominant in the drift and push the process too strongly inward on the interval $[t, t+h]$, which results in inconsistent proposed particles in the algorithm. This phenomenon was already discussed in [Pettersson, 1997] for the approximation of reflected SDE solutions. To address this, the penalty parameter $\lambda$ might be adapted to the frequency of the data.

From a practical perspective this demonstrates that including even simple geometric information about the movement domain can improve the reliability of track reconstruction from noisy telemetry data, particularly when measurement errors are heavy-tailed.

Mean values of the maximum error and RMSE over the 20 trajectories are shown in Table 2.

<table><caption>Table 2: Max error and root mean squared error (RMSE) averaged over the 20 trajectories with Student's t-distributed measurement errors for different values of the time step $h$. The smallest value in each row is marked in bold. The last column shows the average error in the observations (before filtering) compared to the true positions.</caption><thead><tr><th>$h$</th><th>EKF</th><th>Penalized EKF</th><th>Lie-Trotter PF</th><th>Penalized Lie-Trotter PF</th><th>Strang PF</th><th>Penalized Strang PF</th><th>Before filter</th></tr></thead><tbody><tr><td colspan="8" style="text-align:center;">Average max error (km)</td></tr><tr><td>1 min</td><td>1.287</td><td>1.278</td><td>2.304</td><td><strong>0.754</strong></td><td>1.793</td><td>0.757</td><td>5.657</td></tr><tr><td>3 min</td><td>2.010</td><td>2.076</td><td>2.708</td><td>1.334</td><td>2.592</td><td><strong>1.245</strong></td><td>4.694</td></tr><tr><td>5 min</td><td>1.980</td><td>2.037</td><td>2.335</td><td>1.879</td><td>2.334</td><td><strong>1.711</strong></td><td>3.587</td></tr><tr><td>20 min</td><td>2.193</td><td><strong>2.188</strong></td><td>2.387</td><td>2.767</td><td>2.303</td><td>2.522</td><td>2.582</td></tr><tr><td colspan="8" style="text-align:center;">Average RMSE (km)</td></tr><tr><td>1 min</td><td>0.202</td><td>0.194</td><td>0.229</td><td><strong>0.187</strong></td><td>0.217</td><td><strong>0.187</strong></td><td>0.361</td></tr><tr><td>3 min</td><td>0.272</td><td><strong>0.264</strong></td><td>0.328</td><td>0.272</td><td>0.323</td><td>0.271</td><td>0.362</td></tr><tr><td>5 min</td><td>0.305</td><td><strong>0.297</strong></td><td>0.352</td><td>0.323</td><td>0.344</td><td>0.318</td><td>0.357</td></tr><tr><td>20 min</td><td>0.365</td><td>0.365</td><td>0.363</td><td>0.386</td><td><strong>0.356</strong></td><td>0.380</td><td>0.356</td></tr></tbody></table>

## 6.3 Argos errors

We consider the same setup. We add X-shaped Argos measurement errors with scale parameters $\sigma_{obs} = 0.2$, $\rho = 0.7$, $a = 0.4$, $p = 0.5$, and degrees of freedom $d = 3$ to create noisy observations of the trajectories. We run an EKF, a Penalized EKF for which the error distribution is approximated by an isotropic Gaussian distribution with variance $\frac{d}{d-2}\sigma_{obs}^2$. Then, we run a Lie-Trotter PF, a Penalized Lie-Trotter PF, a Strang PF and a Penalized Strang PF with $K = 1000$ particles. As previously, the calculations were performed in parallel over
<figure><img src="image_5.png" /><figcaption>Figure 4: Results of filtering for 20 simulated trajectories with additive Student noise. (a) Boxplots of root mean squared error (RMSE) for EKF, Penalized EKF, Strang PF, Penalized Strang PF, Lie-Trotter PF, Penalized Lie-Trotter PF. (b) Boxplots of max error.</figcaption></figure>
48 CPU cores and took about 18 hours. Figure 5a shows boxplots of the RMSEs. The incorporation of domain knowledge significantly improves the performance of the filter for $h = 1$, 3 and 5 minutes. Figure 4b shows the same boxplots for the maximum error. Overall, the maximum error is reduced by more than 50% with the penalized particle filters compared to the EKF and penalized EKF for $h = 1$ min.

Mean values of the maximum error and RMSE over the 20 trajectories are shown in Table 3.

<table><caption>Table 3: Max error and root mean squared error (RMSE) averaged over the 20 trajectories with Argos X-shaped measurement errors for each value of the time step $h$. The smallest value in each row is marked in bold. The last column shows the average error in the observations (before filtering) compared to the true positions.</caption><thead><tr><th><i>h</i></th><th>EKF</th><th>Penalized<br>EKF</th><th>Lie-<br>Trotter<br>PF</th><th>Penalized<br>Lie-<br>Trotter<br>PF</th><th>Strang PF</th><th>Penalized<br>Strang PF</th><th>Before<br>filter</th></tr></thead><tbody><tr><td colspan="8" style="text-align:center; font-weight:bold;">Average max error (km)</td></tr><tr><td>1 min</td><td>1.493</td><td>1.499</td><td>1.892</td><td><b>0.623</b></td><td>2.03</td><td>0.630</td><td>6.538</td></tr><tr><td>3 min</td><td>2.4</td><td>2.43</td><td>2.25</td><td>0.961</td><td>2.131</td><td><b>0.935</b></td><td>5.575</td></tr><tr><td>5 min</td><td>2.260</td><td>2.310</td><td>2.55</td><td>1.262</td><td>2.386</td><td><b>1.225</b></td><td>3.998</td></tr><tr><td>20 min</td><td>2.265</td><td>2.263</td><td>2.169</td><td>2.128</td><td><b>2.125</b></td><td>5.665</td><td>2.350</td></tr><tr><td colspan="8" style="text-align:center; font-weight:bold;">Average RMSE (km)</td></tr><tr><td>1 min</td><td>0.197</td><td>0.189</td><td>0.193</td><td><b>0.156</b></td><td>0.189</td><td><b>0.156</b></td><td>0.340</td></tr><tr><td>3 min</td><td>0.267</td><td>0.257</td><td>0.259</td><td><b>0.218</b></td><td>0.264</td><td><b>0.218</b></td><td>0.341</td></tr><tr><td>5 min</td><td>0.295</td><td>0.287</td><td>0.294</td><td><b>0.257</b></td><td>0.295</td><td>0.258</td><td>0.340</td></tr><tr><td>20 min</td><td>0.348</td><td>0.347</td><td>0.408</td><td>0.353</td><td>0.396</td><td><b>0.344</b></td><td>0.331</td></tr></tbody></table>

These results illustrate particularly well that a precise modelling of the measurement error, together with an accurate consideration of the spatial constraints on the movement, can help to make the most of available data and reconstruct reliable trajectories. Particle filters with domain penalization reduce both RMSE and maximum error compared to Kalman and Extended Kalman filters. However, this gain comes at a computational cost due to the Monte Carlo nature of the particle filter approximation. Overall, the results suggest that EKF-based methods remain effective baseline approaches for Argos-type telemetry, while penalized particle filters offer a way to improve accuracy when computational resources allow and when reliably handling highly non-Gaussian, anisotropic errors is essential.

# 7 Discussion and conclusion

We introduced a penalized Langevin SDE suitable for modeling the dynamics of animals moving within a constrained domain. This includes potential applications to seals [Hanks et al., 2017], killer whales [Lin et al., 2025], narwhals [Delporte et al., 2025] or elks moving in fenced areas [Brillinger, 2003]. The effect of landscape boundaries is incorporated through an additional drift term that penalizes trajectories stepping too far outside the domain of interest. This approach was originally proposed by Lions and Sznitman [1984] and Liu [1995]
<figure><img src="image_6.png" /><figcaption>Figure 5: Results of filtering for 20 simulated trajectories with additive Argos X-shaped noise. (a) Boxplots of RMSE (b) Boxplots of max error.</figcaption></figure>
to approximate first-order reflected SDEs. Compared with Lin et al. [2025], we include the
domain constraint directly in the latent dynamics and therefore do not require post-processing
of the observed trajectories to ensure land avoidance for aquatic animals.

Simulation of our model relies on splitting the SDE into a non-linear ODE and a linear
SDE, both of which can be solved explicitly [Pilipovic et al., 2024, 2025]. This splitting
scheme is straightforward and easy to implement. Based on this approximate solution to
the non-linear SDE, we introduced filtering methods for noisy observations of the tracks. We
considered the classical Gaussian distribution for the measurement error, as well as more chal-
lenging alternatives such as Student's t or mixtures of multivariate Student's t-distributions,
which are suitable to model Argos telemetry data errors.

For Gaussian errors, we show that the classical Kalman filter and the Extended Kalman
filter perform similarly for high-frequency data. We found that incorporating spatial con-
straints into the latent movement model within the Kalman filters greatly improves filtering
results. Additionally, we proposed a particle filter algorithm to recover true trajectories
from observations with heavy-tailed error distributions. We designed Gaussian proposals
that approximate the optimal proposal in each case, depending on the splitting scheme used
to approximate the SDE. Although they greatly improve filtering accuracy in the case of
non-Gaussian measurement error, we may emphasize that particle filters come at a consid-
erable computational cost compared with Kalman filtering and may therefore be unsuitable
for real-time (online) applications.

The results presented here are application-oriented, and we did not provide theoretical
guarantees for the convergence of the Monte Carlo estimation of filtering distributions. Such
guarantees can be found in Del Moral et al. [2001] for the specific case of first-order SDEs
solved via the Euler scheme. One of the main requirements for convergence to hold is that
the density of the approximation of the solution to the latent SDE converges to the density of
the true process. Establishing such a result is particularly difficult, especially when the latent
process solves a hypoelliptic SDE (an SDE with degenerate noise), which is constrained in a
domain.

A critical parameter in the model is the penalty λ. A small value of λ implies a hard
constraint, as the process is immediately pushed back into the domain when leaving it,
whereas a large value of λ allows the process to leave the boundary for some time and slowly
return to the domain of interest. The choice of λ should primarily depend on the frequency
of the available data. From a theoretical point of view, in first-order reflected SDEs, λ = √h
is often chosen to minimize the RMSE between the true reflected process and its penalized
approximation [Pettersson, 1997]. The choice of λ may also depend on the level of confidence
in the landscape boundaries. If we know with certainty that the animal cannot cross these
boundaries, it is sensible to choose a low value of λ to ensure that the process remains very
close to the ecological domain of interest. Conversely, if there is some uncertainty, a higher
value of λ can allow the process to leave the domain temporarily.

In this work, we assumed that all parameters were known, including the movement pa-
rameters τ, ν, and ω, as well as the potential surface parameters x*, α, and B, and the
parameters of the measurement error model. In practice, the movement parameters and
the potential surface need to be estimated from the data, whereas the measurement error
parameters are generally known, depending on the type of tracking device attached to the
animal. Future work will focus on estimating these parameters from data using Monte Carlo
estimates of the filtering distributions obtained from the particle filters.

Finally, in the latent dynamic, the potential surface *H* may depend on environmental covariates that link areas of attraction with feeding habitats. For instance, the presence of whales in a region is strongly influenced by the availability of zooplankton, which is closely associated with patterns of chlorophyll-a productivity [Panigada et al., 2024]. Such covariate might be included in the latent movement model through resource selection functions [Michelot and Hanks, 2025].
References

Christoffer Moesgaard Albertsen. Generalizing the first-difference correlated random walk for marine animal movement data. *Scientific Reports*, 9(1):4017, March 2019.

David R. Brillinger. Simulating constrained animal motion using stochastic differential equations. In *Institute of Mathematical Statistics Lecture Notes - Monograph Series*, pages 35–48. Institute of Mathematical Statistics, Beachwood, OH, 2003.

Brian M. Brost, Mevin B. Hooten, Ephraim M. Hanks, and Robert J. Small. Animal movement constraints improve resource selection inference in the presence of telemetry error. *Ecology*, 96(10):2590–2597, October 2015.

Alejandro Cholaquidis, Ricardo Fraiman, Ernesto Mordecki, and Cecilia Papalardo. Level sets and drift estimation for reflected Brownian motion with drift. *Statistica Sinica*, 2020. Publisher: Statistica Sinica (Institute of Statistical Science).

Pierre Del Moral, Jean Jacod, and Philip Protter. The Monte-Carlo method for filtering with discrete-time observations. *Probability Theory and Related Fields*, 120(3):346–368, July 2001.

Alexandre Delporte, Susanne Ditlevsen, and Adeline Samson. Varying coefficients correlated velocity models in complex landscapes with boundaries applied to narwhal responses to noise exposure. *The Annals of Applied Statistics*, 2025. To appear.

Arnaud Doucet and Adam Johansen. A tutorial on particle filtering and smoothing: fifteen years later. In *The Oxford Handbook of nonlinear filtering*, pages 656–705. Oxford University Press, 2011.

Pierre Etoré, José Rafael León, and Clémentine Prieur. About hypoelliptic SDE equations with boundary conditions. working paper or preprint, 2025. URL https://hal.science/ hal-04984931.

Eliezer Gurarie, Christen H. Fleming, William F. Fagan, Kristin L. Laidre, Jesús Hernández-Pliego, and Otso Ovaskainen. Correlated velocity models as a fundamental unit of animal movement: synthesis and applications. *Movement Ecology*, 5(1):13, December 2017.

Ephraim M. Hanks, Devin S. Johnson, and Mevin B. Hooten. Reflected Stochastic Differential Equation Models for Constrained Animal Movement. *Journal of Agricultural, Biological and Environmental Statistics*, 22(3):353–372, September 2017.

Graeme C. Hays, Jeanne A. Mortimer, Alex Rattray, Takahiro Shimada, and Nicole Esteban. High accuracy tracking reveals how small conservation areas can protect marine megafauna. *Ecological Applications*, 31(7):e02418, October 2021.

Xavier Hoenner, Scott D. Whiting, Mark A. Hindell, and Clive R. McMahon. Enhancing the Use of Argos Satellite Data for Home Range and Long Distance Migration Studies of Marine Animals. *PLoS ONE*, 7(7):e40713, July 2012.

Devin S. Johnson, Joshua M. London, Mary-Anne Lea, and John W. Durban. Continuous-time random walk model for animal telemetry data. *Ecology*, 89(5):1208–1215, 2008.
Ian D. Jonsen, Toby A. Patterson, Daniel P. Costa, Philip D. Doherty, Brendan J. Godley, W. James Grecian, Christophe Guinet, Xavier Hoenner, Sarah S. Kienle, Patrick W. Robinson, Stephen C. Votier, Scott Whiting, Matthew J. Witt, Mark A. Hindell, Robert G. Harcourt, and Clive R. McMahon. A continuous-time state-space model for rapid quality control of argos locations from animal-borne tags. *Movement Ecology*, 8(1):31, December 2020.

Teng-Wei Lin, Michael Dowd, and Ruth Joy. Forecasting trajectories of Southern Resident killer whales with stochastic movement models incorporating direction modification. *Ecological Modelling*, 509:111254, October 2025.

P. L. Lions and A. S. Sznitman. Stochastic differential equations with reflecting boundary conditions. *Communications on Pure and Applied Mathematics*, 37(4):511–537, July 1984.

Yingjie Liu. Discretization of a class of reflected diffusion processes. *Mathematics and Computers in Simulation*, 38(1-3):103–108, May 1995.

Théo Michelot and Ephraim M. Hanks. Multiscale modelling of animal movement with persistent dynamics, 2025. URL https://arxiv.org/abs/2406.15195.

Théo Michelot, Richard Glennie, Catriona Harris, and Len Thomas. Varying-Coefficient Stochastic Differential Equations with Applications in Ecology. *Journal of Agricultural, Biological and Environmental Statistics*, 26(3):446–463, September 2021.

Viola Panigada, Thomas W. Bodey, Ari Friedlaender, Jean-Noël Druon, Luis A. Huckstädt, Nino Pierantonio, Eduard Degollada, Beatriu Tort, and Simone Panigada. Targeting fin whale conservation in the North-Western Mediterranean Sea: insights on movements and behaviour from biologging and habitat modelling. *Royal Society Open Science*, 11(3): 231783, March 2024.

Toby A. Patterson, Bernie J. McConnell, Mike A. Fedak, Mark V. Bravington, and Mark A. Hindell. Using GPS data to evaluate the accuracy of state-space methods for correction of Argos satellite telemetry error. *Ecology*, 91(1):273–285, January 2010.

Roger Pettersson. Penalization Schemes for Reflecting Stochastic Differential Equations. *Bernoulli*, 3(4):403, December 1997.

Predrag Pilipovic, Adeline Samson, and Susanne Ditlevsen. Parameter Estimation in Nonlinear Multivariate Stochastic Differential Equations Based on Splitting Schemes. *The Annals of Statistics*, 52(2), April 2024.

Predrag Pilipovic, Adeline Samson, and Susanne Ditlevsen. Strang Splitting for Parametric Inference in Second-order Stochastic Differential Equations. *Stochastic Processes and their Applications*, 187:104650, September 2025.

Haiganoush K. Preisler, Alan A. Ager, Bruce K. Johnson, and John G. Kie. Modeling animal movements using stochastic differential equations. *Environmetrics*, 15(7):643–657, November 2004.
Paul J. Wensveen, Len Thomas, and Patrick J. O. Miller. A path reconstruction method integrating dead-reckoning and position fixes applied to humpback whales. *Movement Ecology*, 3(1):31, 2015.
## A Derivation of the gradient and Hessian of the potential

We derive the gradient and Hessian matrix of the potential

$$H(x) = -\sum_{j=1}^{J} H_j(x) \text{ with } H_j(x) = \alpha_j \exp(-(x - x_j^*)^\top B_j(x - x_j^*)).$$

We use the following formulas for matrix and vector differentiation. For any $\phi: \mathbb{R} \to \mathbb{R}$, $a: \mathbb{R}^d \to \mathbb{R}$, $v: \mathbb{R}^d \to \mathbb{R}^d$, $B \in M_{d,d}(\mathbb{R})$,

$$\nabla \phi(a(x)) = \phi'(a(x)) \nabla a(x) \qquad (A.1)$$

$$\nabla(x^\top B x) = (B + B^\top)x \qquad (A.2)$$

$$D(a(x)v(x)) = v(x)\nabla a(x)^\top + a(x)Dv(x) \qquad (A.3)$$

where D denotes the Jacobian matrix.

By (A.1) and (A.2) we get

$$\begin{aligned} \nabla H_j(x) &= -\alpha_j(B_j + B_j^\top)(x - x_j^*) \exp(-(x - x_j^*)^\top B_j(x - x_j^*)) \\ &= -2\alpha_j B_j(x - x_j^*) \exp(-(x - x_j^*)^\top B_j(x - x_j^*)) \end{aligned}$$

where we used that $B_j$ is symmetric. Then, by (A.3),

$$\begin{aligned} D^2 H_j(x) &= D \nabla H_j(x) = -2\alpha_j \left( B_j(x - x_j^*) \times \left( -2B_j(x - x_j^*) \exp(-(x - x_j^*)^\top B_j(x - x_j^*)) \right)^\top \right. \\ & \qquad \left. + \exp(-(x - x_j^*)^\top B_j(x - x_j^*)) \times B_j \right) \\ &= -2\alpha_j \left( -2B_j(x - x_j^*)(x - x_j^*)^\top B_j^\top \exp(-(x - x_j^*)^\top B_j(x - x_j^*)) \right. \\ & \qquad \left. + \exp(-(x - x_j^*)^\top B_j(x - x_j^*)) \times B_j \right) \\ &= -2\alpha_j \exp(-(x - x_j^*)^\top B_j(x - x_j^*)) (B_j - 2B_j(x - x_j^*)(x - x_j^*)^\top B_j) \\ &= -2\alpha_j \exp(-(x - x_j^*)^\top B_j(x - x_j^*)) B_j (I_2 - 2(x - x_j^*)(x - x_j^*)^\top B_j). \end{aligned}$$

## B Full calculation for SDE splitting

Let $l \in \{1, \dots, J\}$. We want to split the following non-linear SDE:

$$\begin{cases} dX(t) = V(t)dt \\ dV(t) = -AV(t) + \nabla H_l(X(t))dt - \nabla H_{-l}(X(t))dt + \sigma dW(t) - \beta_{\lambda}(X(t))dt \end{cases} \qquad (B.1)$$

where $H_{-l}(x) = -\sum_{j=1, j \neq l}^{J} H_j(x)$. We define

$$\begin{aligned} F(x, v) &= -Av + \nabla H_l(x)dt - \nabla H_{-l}(x) - \beta_{\lambda}(x) \\ &= -Av - 2\alpha_l B_l(x - x_l^*)e_l(x) - \beta_{\lambda}(x) \end{aligned}$$
where $e_l(x) = \exp(-(x - x_l^*)^\top B_l(x - x_l^*))$.

(B.1) can be rewritten as

$$dU(t) = \tilde{F}(U(t))dt - \begin{pmatrix} 0_2 \\ \nabla H_{-l}(X(t)) \end{pmatrix} dt + \tilde{\Sigma}dW(t) \quad (B.2)$$

with $\tilde{F}(x, v) = (v \ F(x, v))^\top$

Hence, we can decompose the equation as

$$dU(t) = \tilde{A}(U(t) - u^*)dt + \tilde{N}(x, v) + \tilde{\Sigma}dW(t)$$

with $\tilde{A} = \begin{pmatrix} 0_{2,2} & I_2 \\ A_x & A_v \end{pmatrix}$, $\tilde{N}(x, v) = \begin{pmatrix} 0_2 \\ N(x, v) - \nabla H_{-l}(X(t)) \end{pmatrix}$, $u^* = \begin{pmatrix} x^* \\ 0_2 \end{pmatrix}$ and $A_x, A_v, x^*$, $N(x, v)$ such that

$$F(x, v) = A_x(x - x^*) + A_v v + N(x, v)$$

It is clear that $u_l^* = (x_l^*, 0_{1,2})$ is a zero of $\tilde{F}$. We split the drift as follows: $\tilde{F}(x, v) = \tilde{A}(u - u_l^*) + \tilde{N}(x, v)$ with

$$\tilde{A} = D_u \tilde{F}(u_l^*) = \begin{pmatrix} 0_{2,2} & I_2 \\ -2\alpha_l B_l & -A \end{pmatrix}, \quad \tilde{N}(x, v) = \begin{pmatrix} 0_{2,1} \\ N(x, v) \end{pmatrix}$$

and

$$\begin{aligned} N(x, v) &= F(x, v) - A_x(x - x_l^*) - A_v v \\ &= -Av - 2\alpha_l B_l(x - x_l^*)e_l(x) - \beta_\lambda(x) + 2\alpha_l B_l(x - x_l^*) + Av \\ &= -2\alpha_l B_l(x - x_l^*)(e_l(x) - 1) - \beta_\lambda(x) \end{aligned}$$

Hence we obtain equation (4.1).

## C Covariance matrix for the OU process

The covariance matrix of the OU process in the splitting is:

$$\tilde{Q}(h) = \int_0^h e^{\tilde{A}u} \Gamma e^{-\tilde{A}^\top u} du.$$

Following Albertsen [2019], we vectorize the covariance matrix,

$$\begin{aligned} \mathrm{vec}(\tilde{Q}(h)) &= \int_0^h e^{\tilde{A}u} \otimes e^{\tilde{A}u} \mathrm{vec}(\Gamma) du \\ &= \int_0^h e^{(\tilde{A} \oplus \tilde{A})u} du \mathrm{vec}(\Gamma) \\ &= (\tilde{A} \oplus \tilde{A})^{-1}(e^{(\tilde{A} \oplus \tilde{A})h} - I_{16}) \mathrm{vec}(\Gamma) \\ &= (\tilde{A} \oplus \tilde{A})^{-1} e^{(\tilde{A} \oplus \tilde{A})h} \mathrm{vec}(\Gamma) - (\tilde{A} \oplus \tilde{A})^{-1} \mathrm{vec}(\Gamma) \\ &= e^{(\tilde{A} \oplus \tilde{A})h} (\tilde{A} \oplus \tilde{A})^{-1} \mathrm{vec}(\Gamma) - (\tilde{A} \oplus \tilde{A})^{-1} \mathrm{vec}(\Gamma) \end{aligned}$$
where we used the following properties of the vectorization operator and Kronecker sum and product: for any matrices A,B and C,

$$ \mathrm{vec}(ABC) = (C^{\top} \otimes A)\mathrm{vec}(B) $$

$$ e^{M} \otimes e^{N} = e^{M \oplus N}. $$

Let C be the matrix such that $\mathrm{vec}(C) = (\tilde{A} \oplus \tilde{A})^{-1}\mathrm{vec}(\Gamma)$. Then

$$ \mathrm{vec}(\tilde{Q}(h)) = e^{(\tilde{A} \oplus \tilde{A})h} \mathrm{vec}(C) - \mathrm{vec}(C) = e^{\tilde{A}h} \otimes e^{\tilde{A}h} \mathrm{vec}(C) - \mathrm{vec}(C) = \mathrm{vec}(e^{\tilde{A}h}Ce^{\tilde{A}^{\top}h} - C) $$

$$ \text{Hence } \tilde{Q}(h) = e^{\tilde{A}h}Ce^{\tilde{A}^{\top}h} - C. $$