Multivariate Signal Modelling with Applications to Inertial Sensor Calibration

TL;DR

This paper introduces a novel multivariate modeling approach for inertial sensor calibration that captures dependencies between sensors using Wavelet Cross-Covariance, enhancing calibration accuracy and enabling better navigation system performance.

Contribution

It develops a new method extending the Generalized Method of Wavelet Moments to model both individual signals and their dependencies, improving calibration and analysis of inertial sensors.

Findings

Method effectively captures sensor dependencies.

Improves calibration accuracy in simulations.

Enables testing of sensor dependence.

Abstract

The common approach to inertial sensor calibration for navigation purposes has been to model the stochastic error signals of individual sensors independently, whether as components of a single inertial measurement unit (IMU) in different directions or arrayed in the same direction for redundancy. These signals usually have an extremely complex spectral structure that is often described using latent (or composite) models composed by a sum of underlying models. A large amount of research in this domain has been focused on the latter aspect through the proposal of various methods that have been able to improve the estimation of these models both from a computational and a statistical point of view. However, the separate calibration of the individual sensors is still unable to take into account the dependence between each of them which can have an important impact on the precision of the navigation systems. In this paper, we develop a new approach to simultaneously model both the individual signals and the dependence between them by studying the quantity called Wavelet Cross-Covariance and using it to extend the application of the Generalized Method of Wavelet Moments. This new method can be used in other settings for time series modeling, especially in cases where the dependence among signals may be hard to detect. Moreover, in the field of inertial sensor calibration, this approach can deliver important contributions among which the possibility to test dependence between sensors, integrate their dependence within the navigation filter and construct an optimal virtual sensor that can be used to simplify and improve navigation accuracy. The advantages of this method and its usefulness for inertial sensor calibration are highlighted through a simulation study and an applied example with a small array of XSens MTi-G IMUs.

Equations

$${\mathbf{X}}_t\equiv \sum _{p=1}^P{\mathbf{A}}_p{\mathbf{X}}_{t-p}+\sum _{q=1}^Q{\mathbf{U}}_q{\mathbf{ϵ}}_{t-q}+{\mathbf{ϵ}}_t,$$

$${\mathbf{ϵ}}_t\stackrel{\text{i.i.d.}}{\sim }\mathcal{F}(0,\mathbf{\Sigma }),$$

$${\mathbf{X}}_t\equiv \sum _{k=1}^K{\mathbf{S}}_k{\mathbf{X}}_{k,t},$$

$$(\mathbf{S}_{k})_{i,j}=\left\{\begin{array}[]{ll}1&\mbox{if }X_{k,t}^{(j)}\mbox{ appears in }X_{t}^{(i)}\\ 0&\mbox{if }\ \text{otherwise}.\end{array}\right.$$

$${\mathbf{X}}_t=\sum _{k=1}^K{S}_k{\mathbf{X}}_{k,t},$$

$$S_{k}=\left\{\begin{array}[]{ll}1&\mbox{if }X_{k,t}\mbox{ appears in }X_{t}\\ 0&\mbox{if }\ \text{otherwise}.\end{array}\right.$$

$${\mathbf{X}}_{2,t}={\mathbf{X}}_{2,t-1}+{\mathbf{\iota }}_t,$$

$$X_{4,t}^{(i)}={\omega }_{i}t,$$

$${\mathbf{X}}_{k,t}=\mathbf{\Phi }{\mathbf{X}}_{k,t-1}+{\mathbf{\epsilon }}_t,$$

$$W_{j,t}^{(i)}\equiv \sum _{l=0}^{L_j-1}{h}_{j,l}{X}_{t-l}^{(i)},$$

$${\eta }_j^{(i)}({\mathbf{\theta }}^{(i)})\equiv \text{Var}(W_{j,t}^{(i)}),$$

$${\mathbf{\theta }}^{(i)}\equiv {\left[{\theta }_m^{(i)}\right]}_{m=1,\dots ,p^{(i)}}$$

$${\gamma }_{j,h}^{(i,i^{\prime })}({\mathbf{\theta }}^{(i,i^{\prime })})\equiv \text{Cov}\left(W_{j,t}^{(i)},W_{j,t+h}^{(i^{\prime })}\right),$$

$$\mathbf{\theta }\equiv {\left[{\left[{\mathbf{\theta }}^{(1)}\right]}^{⊺},\dots ,{\left[{\mathbf{\theta }}^{(I-1,I)}\right]}^{⊺}\right]}^{⊺},$$

$$p^{(i,i^{\prime })}=\sum _{m=1}^{p^{(i)}}\sum _{n=1}^{p^{(i^{\prime })}}\text{mathbbm}{1}_{\left\{{\theta }_m^{(i)}={\theta }_n^{(i^{\prime })}\right\}},$$

$${\hat{\gamma }}_{j,0}^{(i,i^{\prime })}\equiv \frac{1}{M_j}\sum _{t=L_j}^T{W}_{j,t}^{(i)}{W}_{j,t}^{(i^{\prime })},$$

$${\left\{\left(W_{j,t}^{(i)},W_{j,t}^{(i^{\prime })}\right)\right\}}_{t\in \mathbb{Z}}$$

$$\mathbf{\nu }({\mathbf{\theta }}_0)\equiv {\left[{\gamma }_j^{(i,i^{\prime })}\right]}_{\begin{array}{c}j=1,...,J\\ 1\le i\le i^{\prime }\le I\end{array}},$$

$$\hat{\mathbf{\nu }}\equiv {\left[{\hat{\gamma }}_j^{(i,i^{\prime })}\right]}_{\begin{array}{c}j=1,...,J\\ 1\le i\le i^{\prime }\le I\end{array}}.$$

$$W_{j,t}^{(i)}={\mathbf{c}}_j^{⊺}{\mathbf{\Delta }}_{j,t}^{(i)},$$

$$\mathbf{G}(\cdot )={\left(g^{(1)}(\cdot ),g^{(2)}(\cdot ),\dots ,g^{(I)}(\cdot )\right)}^{⊺},$$

$${\mathbf{\Delta }}_t\equiv {\left({\Delta }_t^{(1)},{\Delta }_t^{(2)},\dots ,{\Delta }_t^{(I)}\right)}^{⊺}=\mathbf{G}({\mathcal{F}}_t).$$

$${\mathbf{W}}_t\equiv {\left[W_{j,t}^{(i)}{W}_{j,t}^{(i^{\prime })}\right]}_{\begin{array}{c}j=1,...,J\\ 1\le i\le i^{\prime }\le I\end{array}}$$

$${\mathcal{P}}_t(\cdot )\equiv \mathbb{E}(\cdot |{\mathcal{F}}_t)-\mathbb{E}(\cdot |{\mathcal{F}}_{t-1}),$$

$$\sqrt{T}\left(\hat{\mathbf{\nu }}-\mathbf{\nu }({\mathbf{\theta }}_0)\right)\stackrel{\mathcal{D}}{\to }\mathcal{N}\left(0,\mathbf{V}\right),$$

$$\sqrt{T}(\tilde{\mathbf{\nu }}-\mathbf{\nu }({\mathbf{\theta }}_0))+\sqrt{T}(\hat{\mathbf{\nu }}-\tilde{\mathbf{\nu }}),$$

$${\tilde{\gamma }}_j^{(i,i^{\prime })}\equiv \frac{1}{T}\sum _{t=1}^T{W}_{j,t}^{(i)}{W}_{j,t}^{(i^{\prime })},$$

$${\hat{\mathbf{\theta }}}^{(i)}=\underset{\mathbf{\theta }\in {\mathbf{\Theta }}^{(i)}}{\mathrm{argmin}}{\left({\hat{\mathbf{\eta }}}^{(i)}-{\mathbf{\eta }}^{(i)}(\mathbf{\theta })\right)}^{⊺}{\mathbf{\Omega }}^{(i)}\left({\hat{\mathbf{\eta }}}^{(i)}-{\mathbf{\eta }}^{(i)}(\mathbf{\theta })\right),$$

$$\hat{\mathbf{\theta }}=\underset{\mathbf{\theta }\in \mathbf{\Theta }}{\mathrm{argmin}}{\left(\hat{\mathbf{\nu }}-\mathbf{\nu }(\mathbf{\theta })\right)}^{⊺}\mathbf{\Omega }\left(\hat{\mathbf{\nu }}-\mathbf{\nu }(\mathbf{\theta })\right).$$

$$\mathbf{\nu }({\mathbf{\theta }}_0)=\mathbf{\nu }({\mathbf{\theta }}_1)\text{if and only if}{\mathbf{\theta }}_0={\mathbf{\theta }}_1,$$

Full text

I Introduction

The modelling of multivariate time series is an important as well as challenging task in many applied domains of data analysis. Indeed, it is extremely common to measure phenomena over time which not only manifest autocorrelation with their past but also show a form of dependence between them. This setting can be observed, among others, in a variety of longitudinal studies [1, 2], in economic and financial research [3, 4, 5, 6] as well as in different fields of engineering [7, 8]. In the latter case, an increased attention has been given to the modelling of multivariate signals for the task of inertial sensor calibration where the error signals of the individual accelerometers and gyroscopes that compose the Inertial Measurement Units (IMU) are characterized by error signals that are dependent on each other [8]. To date, these error signals have been modeled independently from each other and this task alone has already been challenging to deal with. In fact, the individual error signals of these instruments are characterized by deterministic and stochastic components where the first can be dealt with through physical models (which are employed for calibration and compensation prior to use) while the latter are often described through complex stochastic models that are made by the sum of different underlying processes. The estimation of these stochastic models has always been complicated and many methods used for this purpose have either suffered from statistical limits and/or are computationally or numerically unstable. Indeed, in most latent variable modelling scenarios, the Maximum-Likelihood Estimator (MLE) cannot be directly applied since the latent stochastic components are unobservable and their marginal distribution is usually unknown. A feasible likelihood-based method is the Expectation-Maximization (EM) algorithm [9] whose implementation is nevertheless limited due to the complexity and, more importantly, the possible non-convexity of its associated optimization problems therefore severely limiting the use of its statistical properties in practice [10]. To avoid the problems of the likelihood-based methods, the linear regression based on the Allan Variance (AVLR) is widely used for estimation and prediction especially in the field of inertial sensor calibration [11]. Although, the AVLR is a computationally feasible method for estimating many complex models, its statistical properties are not optimal as discussed in [12] where the inconsistency of this method was shown for the majority of the latent models used for the purpose of sensor calibration.

Given the complexity of these models and the difficulties in estimating them, the signals have thus far been dealt with individually without taking into account the forms of dependence between them. As mentioned earlier, this approach is non-optimal since a multivariate approach to signal analysis would be more appropriate for different reasons. Firstly, it is reasonable to assume that there exists a dependence between individual sensors along the three axes of an IMU and, when using an array of sensors, it is even more reasonable to assume that there is dependence between error signals of sensors placed in the same direction since they are measuring the same phenomenon [13]. If the multivariate dependence is not considered, this can deliver different problems which include parameter misestimation and flawed testing procedures leading to incorrect modelling and degraded navigation performance. Considering this, the literature in the area of multivariate time series analysis has collected a variety of models to describe these settings as well as proposing methods that can estimate them. One of the first models to be proposed was the Vector Auto-Regression (VAR) model which, in some sense, generalizes the Auto-Regression (AR) model in order to allow lagged observations from some signals to influence current observations of others [14]. This class of models has been extended to the more general class of Vector Auto-Regression Moving Average (VARMA) models which also allow for correlation between innovation processes in each time series [14]. Although these models are quite flexible in considering different forms of autocorrelation and correlation between time series, their estimation can sometimes be numerically challenging and furthermore they may not be able to capture more complicated forms of multivariate dependence [15, 16]. Moreover, in practice, VARMA models are only used directly when the number of signals composing the multivariate time series is smaller than three since these models are often overparametrized thereby leading to identifiability issues [17]. In this case, a more appropriate approach could be represented by only considering a model (or some models) to describe the dependence between all or some subsets of the signals composing the multivariate time series. Indeed, the individual time series can often be characterized by a latent model structure in which different underlying processes are superposed and, within this setting, there can also be a dependence of these individual latent processes with those present in other signals [18, 19, 5]. We refer to these models as multivariate latent models which are particularly relevant for the task of inertial sensor calibration.

Considering the above setting, this paper proposes a new method that is able to estimate the mentioned multivariate latent models in a computationally efficient and numerically stable manner. This new method consists in an extension of the Generalized Method of Wavelet Moments (GMWM) which was initially proposed as a statistically consistent approach to estimating complex latent models in univariate settings, overcoming the theoretical and computational limitations of the existing methods mentioned earlier [18]. To do so the GMWM takes advantage of a quantity called the Wavelet Variance (WV) which is the variance of the wavelet coefficients resulting from a wavelet decomposition of a signal [20]. However, the WV does not take into account the variability in the signal that is explained by the dependence on other signals. For this reason, in this paper we will make use of the Wavelet Cross-Covariance (WCCV), which was proposed in [21] and whose asymptotic properties we develop further in this paper by reducing and weakening the relative conditions. Based on the latter properties, we extend the GMWM to include this quantity and deliver the Mulitvariate GMWM (MGMWM) whose statistical properties are studied in order to perform the required inference procedures. In order to achieve the latter properties, we will also discuss the identifiability of a wide class of latent multivariate time series models which is essential to obtain consistency of the MGMWM, thereby shedding light on what kind of latent multivariate models can be estimated and reducing the conditions necessary to deliver adequate inference procedures on them. Based on these properties, the proposed framework allows to model dependence between and within sensors as well as provide the tools to test the presence of multivariate dependence between stochastic error signals thereby facilitating sensor calibration and contributing to navigation accuracy.

This paper is organized as follows. In Section II we introduce the class of multivariate latent processes that we are going to propose and study, justifying their relevance for the task of inertial sensor calibration. In Section III we briefly introduce the wavelet decomposition to then define and study the WCCV thereby delivering the asymptotic properties of the WCCV vector. The latter properties are essential, among others, in order to obtain the asymptotic properties of the MGMWM which is introduced in Section IV where the identifiability of a large class multivariate latent models is studied. To highlight the good statistical properties of the new estimator as well as its usefulness, Section V presents a simulation study comparing the finite sample performance of the proposed estimator with a recently proposed alternative while Section VI presents an application in the field of inertial sensor calibration where this new approach can deliver an important contribution. Finally, Section VII concludes.

II Multivariate Latent Processes

Let us define a multivariate process as $(\mathbf{X}_{t})_{t\in\mathbb{Z}}\subset\mathbb{R}^{I}$, where $I\in\mathbb{N}$. As mentioned in the introduction, one of the most common approaches to modelling these multivariate process is through the use of VARMA models which have the following structure:

[TABLE]

where $P$ and $Q$ are the maximum lags of dependence for the AR and MA processes respectively among the individual time series composing the multivariate process, $\mathbf{A}_{p}$ and $\mathbf{U}_{q}$ are the coefficient matrices for the $p^{th}$ and $q^{th}$ lag and

[TABLE]

are independently distributed innovation vectors with multivariate distribution $\mathcal{F}$, expectation $\mathbf{0}$ and covariance matrix $\mathbf{\Sigma}$.

These multivariate time series models are very useful to describe a wide variety of phenomena and, at an individual level, can often be represented as a latent process made by the sum of underlying first-order autoregressive and white noise models [22]. However, these models hide the latent structure which can often be extremely useful for interpretation and, more importantly, are limited to a class of latent models in which many processes that are relevant to many domains cannot be included. For this reason, in this paper we consider a different class of multivariate latent models (which can possibly include reparametrizations of the aforementioned VARMA models) in order for practitioners to use these in a variety of settings where they are of great importance (e.g. inertial sensor calibration). To do so, let us define a multivariate time series composed by $K$ latent processes as

[TABLE]

where $\mathbf{X}_{t}\equiv[X_{t}^{(i)}]_{i=1,\ldots,I}$ is an $I$-dimensional vector and $\mathbf{X}_{k,t}\equiv[X_{k,t}^{(j)}]_{j=1,\ldots,D_{k}}$ is a $D_{k}$-dimensional vector with $k$ being the index for the type of univariate/multivariate model underlying some or all of the latent processes and $\mathbf{S}_{k}$ is an $I\times D_{k}$ matrix with either $1$ or [math] as its elements that we define as

[TABLE]

[TABLE]

where $\mathbf{X}_{t}\equiv[X_{t}^{(i)}]_{i=1,\ldots,I}$ is an $I$-dimensional vector and $\mathbf{X}_{k,t}\equiv[X_{k,t}^{(i)}]_{i=1,\ldots,I}$ is a $I$-dimensional vector, with $k$ being the index for the type of multivariate model underlying some or all of the latent processes and $S_{k}$ is a deterministic variable with either $1$ or [math] as its elements, which defined as

[TABLE]

Therefore, we have a set of $K$ multivariate time series models (which can characterize one, some or all of the univariate processes) that can be summed to compose a multivariate latent model. Given this, let us introduce the individual multivariate time series models, commonly used also for the task of inertial sensor calibration (see e.g. [23]), that will be considered in this paper to build a general multivariate latent model:

1. (T1)

White Noise (WN) which corresponds to the independently distributed innovation vectors $\bm{\epsilon}_{t}$ defined in (1). We denote the process characterized by this model as $\mathbf{X}_{1,t}$, hence $\bm{X}_{1,t}\overset{\text{i.i.d.}}{\sim}\mathcal{F}(\mathbf{0},\bm{\Sigma})$. 2. (T2)

Random Walk (RW) which we denote as $\mathbf{X}_{2,t}$ and is defined as

[TABLE]

where $\bm{\iota}_{t}\overset{\text{i.i.d.}}{\sim}\mathcal{F}(\mathbf{0},\bm{\Lambda})$. 3. (T3)

Quantization Noise (QN) (or rounding error, see e.g. [24]) which we denote as $\mathbf{X}_{3,t}$. For this process we do not consider a multivariate model and therefore each univariate process is characterized by the parameter $Q_{i}^{2}\in\{x\in\mathbb{R}|x>0\}$. 4. (T4)

Drift (DR) which we denote as $\mathbf{X}_{4,t}$ and, for each univariate process, is defined as

[TABLE]

where $\omega_{i}\in\{x\in\mathbb{R}|x>0\}$. 5. (T5)

First-Order Auto-Regressive Noise (AR1) which we denote as $\mathbf{X}_{k,t}$, $\forall\,\,k=5,\ldots,K$, where $K\in\mathbb{N}$, and defined as

[TABLE]

where $\bm{\Phi}$ is a diagonal matrix with diagonal elements $\phi_{k}^{(i)}$, $0<|\phi_{k}^{(i)}|<1$, for $i=1,\dots,D_{k}$ and $\bm{\varepsilon}_{t}\overset{\text{i.i.d.}}{\sim}\mathcal{F}(\mathbf{0},\bm{Z}_{k})$.

The covariance matrices $\bm{\Sigma}$, $\bm{\Lambda}$ and $\bm{Z}_{k}$ defined above are positive definite $D_{k}\times D_{k}$ matrices with $k=1,2,5,...,K$ respectively, with their respective elements being $[\sigma^{(i,i^{\prime})}]_{1\leq i\leq i^{\prime}\leq D_{1}}$, $[\lambda^{(i,i^{\prime})}]_{1\leq i\leq i^{\prime}\leq D_{2}}$ and $[z_{k}^{(i,i^{\prime})}]_{1\leq i\leq i^{\prime}\leq D_{k}}$. As can be seen from the above definitions, processes (T3) and (T4) do not have a multivariate structure (i.e. there is no dependence between univariate signals based on these models). This is a reasonable assumption in practice since process (T3) can be interpreted as a form of rounding error which is intuitively independent from the rounding error made on another signal while the (T4) process is a non-stochastic process whose behavior in time is deterministic in nature (hence the eventual dependence would be dealt with through deterministic models). Another aspect to underline regarding the above defined processes is that the dependence between univariate processes is only based on them sharing the same multivariate model (e.g. a signal cannot be dependent on another if they do not share at least one of these models). Although this assumption may restrict the choice of structures of multivariate dependence, this setting is nevertheless sufficiently general to at least well-approximate more complex forms of multivariate latent dependence while providing an important basis to extend this setting in future research. Based on this assumption, we state the following condition:

The covariance matrices $\bm{\Sigma}$, $\bm{\Lambda}$ and $\bm{Z}_{k}$ defined above are positive definite $D_{k}\times D_{k}$ matrices with $k=1,2,5,...,K$ respectively, with their respective elements being $[\sigma_{ij}]_{1\leq i\leq j\leq D_{1}}$, $[\lambda_{ij}]_{1\leq i\leq j\leq D_{2}}$ and $[z_{ij}]_{1\leq i\leq j\leq D_{k}}$ ( make the notations consistent with the ones in the appendix…). As can be seen from the above definitions, processes (T3) and (T4) do not have a multivariate structure (i.e. there is no dependence between univariate signals based on these models). This is a reasonable assumption in practice since the (T3) process is a rounding error which is intuitively independent from the rounding error made on another signal while the (T4) process is a non-stochastic process whose behavior in time is deterministic in nature (hence the eventual dependence would be dealt with through deterministic models). Another aspect to underline regarding the above defined processes is that the dependence between univariate processes is only based on them sharing the same multivariate model (e.g. a signal cannot be dependent on another if they don’t share at least one of these models). Although this assumption may restrict the choice of structures of multivariate dependence, this setting is nevertheless sufficiently general to at least well-approximate more complex forms of multivariate latent dependence while providing an important basis to extend this setting in future research. Based on this assumption, we state the following condition:

1. (C1)

$\bm{X}_{k,t}$ is independent of $\bm{X}_{k^{\prime},t}$, $\forall\,\,k\neq k^{\prime}$.

This condition implies that each model that composes a latent model (univariate or multivariate) is independent from the others (e.g. a (T1) process is independent from a (T5) process and two different (T5) processes are independent from each other). The reason for stating this condition is simply to provide proofs on the identifiability of these models for the method that is put forward in this paper (see Sec. IV-B). It is of course possible to estimate other forms of dependence in which the processes that compose a latent model are correlated, but this would imply deriving additional theoretical moments and, if not assumed, proving identifiability of these models on a case-by-case basis.

With this premise, in the next sections we present the quantities and relative asymptotic properties that are necessary to define the MGMWM. Using the modelling framework presented in this section, we will then be able to deliver the identifiability of these models through the proposed estimation method and consequently reduce the conditions for its asymptotic properties to some basic regularity conditions.

III The Wavelet Cross Covariance

In order to present the proposed MGMWM, we first describe and study the properties of the WCCV. For this purpose, let $J>0$ be a fixed integer representing the last level of wavelet decomposition of the univariate time series $(X_{t}^{(i)})_{t\in\mathbb{Z}}$. This allows us to define the wavelet coefficients for level $j=1,\dots,J$ as:

[TABLE]

where $[h_{j,l}]_{l=0,\dots,L_{j}-1}$ is the $j^{th}$ level wavelet filter with length $L_{j}$. In this paper we choose to use the Haar wavelet filter whose length is consequently $L_{j}=2^{j}$. The reason for using this particular filter resides mainly in the fact that, aside from being one of the most commonly used filters in practice, the theoretical form of the WV (see below) for many classes of time series models have been explicitly derived (see e.g. [zhang2008allan]) thereby facilitating the finding of theoretical results presented further on in the paper.

With the above discussion in mind, we can define the WV as:

[TABLE]

where

[TABLE]

is the parameter vector defining the latent model underlying the $i^{th}$ univariate process (of dimension $p^{(i)}$). Indeed, for each time series model there is a corresponding theoretical WV which can be used for estimation purposes in a method-of-moments fashion (as done by the GMWM). In a similar way, we can also define the WCCV [21] at the $j^{th}$ level as follows:

[TABLE]

where $h$ represents the lag in time between observations and $\bm{\theta}^{(i,i^{\prime})}$ represents the parameter vector defining the dependence between signals (of dimension $p^{(i,i^{\prime})}$). Given our definitions, the parameter vector of the overall model is therefore $\bm{\theta}\in\bm{\Theta}\subset\mathbb{R}^{p}$ which can be defined as

[TABLE]

where $p\equiv\sum_{i=1}^{I}\sum_{j=i}^{I}p^{(i,j)}$ and where $p^{(i,i)}\equiv p^{(i)}$.

Given our definitions, this implies that the dimension of the parameter vector $\bm{\theta}^{(i,i^{\prime})}$ is

[TABLE]

since two processes share the same model if and only if they share the same parameters.

I am not convinced by this definition, could you please double check it. We do not use it anymore, maybe we can remove it – HX

Being a covariance, we can see how this quantity is a generalization of the WV since, setting $h=0$ and $i=i^{\prime}$, Eq. (5) reduces to the WV $\eta_{j}^{(i)}(\bm{\theta}^{(i)})=\gamma_{j,0}^{(i,i)}(\bm{\theta}^{(i)})$.

Having defined these theoretical quantities, we can now define their empirical counterparts that are estimated on the observed signals $(X_{t}^{(i)})_{t=1,\dots,T}$ and $({X_{t}^{(i^{\prime})}})_{t=1,\dots,T}$, where $T\in\mathbb{N}$. While the unbiased estimator of WV is defined and studied in [20] and [25], here we define the unbiased estimator of WCCV [21] for $h=0$ as:

[TABLE]

where $M_{j}\equiv T-L_{j}+1$ is the number of wavelet coefficients generated at level $j$. In this paper we limit ourselves to studying the WCCV at lag $h=0$ since this lag is sufficient to identify and estimate the models considered for this work.

Having defined the theoretical and empirical WCCV, we will now study the asymptotic properties of $\hat{\gamma}_{j,0}^{(i,i^{\prime})}$ since they are essential to deliver inference when using the MGMWM which is presented in Sec. IV. In [26] the asymptotic normality of $\hat{\gamma}_{j,{0}}^{(i,i^{\prime})}$ was stated for any pair $(i,i^{\prime})$ and any level $j$. However, the conditions for the latter results to hold assume that the process

[TABLE]

is a bivariate Gaussian stationary process which may be a strong assumption to satisfy in general. In this paper we aim at reducing (or weakening) these conditions as well as clearly stating the multivariate asymptotic properties of the WCCV vector. To define the latter, we first simplify notation as follows: $\hat{\gamma}_{j}^{(i,i^{\prime})}\equiv\hat{\gamma}_{j,0}^{(i,i^{\prime})}$ and $\gamma_{j}^{(i,i^{\prime})}\equiv\gamma_{j,0}^{(i,i^{\prime})}(\bm{\theta}^{(i,i^{\prime})})$. Based on this notation, the WCCV vector is defined as:

[TABLE]

with the corresponding vector of estimated WCCV defined as:

[TABLE]

In addition, let us also define the first order difference of the $i^{th}$ signal as $(\Delta_{t}^{(i)})_{t\in\mathbb{Z}}$ (i.e. $\Delta_{t}^{(i)}\equiv X_{t}^{(i)}-X_{t-1}^{(i)}$). The reason for considering this process lies in the fact that the wavelet coefficients can be represented as particular linear combinations of the process $(\Delta_{t}^{(i)})$, i.e.

[TABLE]

where $\bm{\Delta}_{j,t}^{(i)}\equiv\left(\Delta_{t-L_{j}+2}^{(i)},\,\Delta_{t-L_{j}+3}^{(i)},\,\dots,\,\Delta_{t}^{(i)}\right)^{\intercal}$ and $\bm{c}_{j}$ represents a specific filter vector. In order to define the first condition, we also define

[TABLE]

as being an $\mathbb{R}^{I}$-valued measurable function as well as the filtration $\mathcal{F}_{t}=\sigma(\dots,\,\epsilon_{t-1},\,\epsilon_{t})$ where $\epsilon_{t}$ are i.i.d. random variables. We can now state the first condition needed for the asymptotic properties of the WCCV vector.

1. (C2)

We assume the multivariate process $(\bm{\Delta}_{t})_{t\in\mathbb{Z}}$ is well defined and can be represented as

[TABLE]

This condition basically implies that the multivariate process $(\bm{\Delta}_{t})$ is ergodic and strictly stationary. Hence, this condition is quite a strong one but is commonly assumed (and often implied by other conditions) to study asymptotic properties in the context of dependent processes. Moreover, this condition is respected by the processes defined in Sec. II and allows for a very general class of time series models e.g. [27].

In order to state the final conditions, we need to deliver additional definitions, starting from $\Delta_{t}^{\star\,(i)}=g^{(i)}(\mathcal{F}_{t}^{\star})$ where $\mathcal{F}_{t}^{\star}=\sigma(\dots,\,\epsilon_{0}^{\star},\,\dots,\,\epsilon_{t-1},\,\epsilon_{t})$, where also $\epsilon_{0}^{\star}$ is an i.i.d. random variable. Based on this definition, we can notice that for $t<0$ we have that $\bm{\Delta}_{t}^{\star}=\bm{\Delta}_{t}$. Finally, we define $\|Z\|_{p}\equiv\left(\mathbb{E}|Z|^{p}\right)^{\nicefrac{{1}}{{p}}}$, for $p>0$, to deliver the other required conditions.

1. (C3)

$\underset{i}{\max}\left\|\Delta_{t}^{(i)}\right\|_{4}<\infty$. 2. (C4)

$\underset{i}{\max}\sum_{t=0}^{\infty}\left\|\Delta_{t}^{(i)}-\Delta_{t}^{\star\,(i)}\right\|_{4}<\infty$.

Both conditions require the fourth moment of the difference process $(\Delta_{t}^{(i)})_{t\in\mathbb{Z}}$ (or some function of it) to be finite (for all $i$). More specifically, while Condition (C3) is straightforward to interpret in the latter sense, Condition (C4) ensures that the cumulative impact of $\epsilon_{0}^{\star}$ on future values of the process $(\Delta_{t}^{(i)})_{t\in\mathbb{Z}}$ is finite which allows us to interpret it as a short-range dependence condition [28]. We can now deliver the main result of this section and, for this purpose, we further define

[TABLE]

and the projection operator

[TABLE]

which measures how much the conditional expectation of a random variable changes when removing the information from the immediately previous time. Intuitively, this measure should not be large for a “stable” signal since these conditional expectations (that contribute to $\mathcal{P}_{t}(\cdot)$) should be arbitrarily close to the unconditional expectation. Using these last definitions, with “$\overset{\mathcal{D}}{\to}$” denoting convergence in distribution, we can finally state the asymptotic normality of the vector $\hat{\bm{\nu}}$.

Theorem 1:

Under Conditions (C2) to (C4), and assuming $I,J\in\mathbb{N}$, we have that

[TABLE]

where $\bm{V}=\mathbb{E}(\bm{D}_{0}\bm{D}_{0}^{\intercal})$ and $\bm{D}_{0}\equiv\sum_{t=0}^{\infty}\mathcal{P}_{0}\left(\bm{W}_{t}\right)$.

Proof: We can split the quantity $\sqrt{T}\left(\hat{\bm{\nu}}-\bm{\nu}(\bm{\theta}_{0})\right)$ as follows:

[TABLE]

where $\tilde{\bm{\nu}}$ is an alternative estimator of the WCCV whose elements are given by

[TABLE]

for all $1\leq i\leq i^{\prime}\leq I$ and $j=1,...,J$. Using the latter estimator and Conditions (C2) to (C4), we can directly apply the results of [28, Theorem 7] to show asymptotic normality of $\sqrt{T}(\tilde{\bm{\nu}}-\bm{\nu}(\bm{\theta}_{0}))$ with covariance matrix $\mathbb{E}(\bm{D}_{0}\bm{D}_{0}^{\intercal})$. Moreover, with $I$ and $J$ being fixed, we have that $\sqrt{T}(\hat{\bm{\nu}}-\tilde{\bm{\nu}})=\mathcal{O}_{p}(\nicefrac{{1}}{{T}})=o_{p}(1)$ and therefore, using Slutsky’s theorem, we have that asymptotically $\sqrt{T}\left(\hat{\bm{\nu}}-\bm{\nu}(\bm{\theta}_{0})\right)$ is equivalent in distribution to $\sqrt{\mathrm{T}}(\tilde{\bm{\nu}}-\bm{\nu}(\bm{\theta}_{0}))$ thus concluding the proof. A more detailed explanation of this proof can be found in App. LABEL:proof:thm:asy.nuhat. ∎

Having studied the asymptotic properties of the WCCV, the next section develops the MGMWM which makes use of this quantity and its properties for inference.

IV The Multivariate GMWM

Given the results presented in the previous section, we can now introduce the proposed approach to multivariate latent signal modelling which develops from the GMWM [18]. Indeed the latter approach considers a univariate process $(X_{t}^{(i)})_{t\in\mathbb{Z}}$ and, assuming that this process is generated from a parametric model $\bm{F}_{\bm{\theta}^{(i)}}$ known up to the parameter vector $\bm{\theta}^{(i)}$, the GMWM aims at estimating this parameter vector through the following minimization problem:

[TABLE]

where $\hat{\bm{\eta}}^{(i)}$ and $\bm{\eta}^{(i)}(\bm{\theta})$ are the vectors of estimated and theoretical WV (for the $i$-th process) respectively, $\bm{\Theta}^{(i)}\subset\mathbb{R}^{p^{(i)}}$ is the parameter space for $\bm{\theta}^{(i)}$ and $\bm{\Omega}^{(i)}\in\mathbb{R}^{J\times J}$ is a symmetric positive definite weighting matrix chosen in a suitable manner. For the latter matrix, one choice can be the inverse of the covariance matrix of $\hat{\bm{\eta}}^{(i)}$ [18]. The consistency and asymptotic normality of $\hat{\bm{\theta}}^{(i)}$ were shown in [18] for any suitable choice of $\bm{\Omega}^{(i)}$, thereby providing not only a statistically appropriate method for univariate latent time series modelling but also a computationally feasible and efficient approach.

The idea of the GMWM, being based on the Generalized Method of Moments (GMM) framework, is therefore to inverse the mapping from the WV vector back to the parameter vector $\bm{\theta}$. With this in mind, this idea can be extended by using the vector of WCCV which includes moments that contain information on the subset of parameters in $\bm{\theta}$ that contain information on the latent dependence between signals $i$ and $i^{\prime}$. For this reason, it is possible to build the MGMWM as an extension of the GMWM as follows:

[TABLE]

Hence, all that is needed to obtain the MGMWM is to postulate a multivariate latent model and compute the empirical WCCV for which we have studied the asymptotic properties in Sec. III. Moreover, the objective function defined in (9) can be used to test whether considering dependence between sensor errors better fits the observed signals by, for example, comparing it to the distribution of the objective function defined in (8) (which could be obtained via parametric bootstrap). However, substituting the WV in the GMWM with the WCCV in the MGMWM is not enough to obtain the necessary inference on the parameter vector $\bm{\theta}$. Indeed, we firstly would need to have an idea on which models it is possible to estimate using this approach and, secondly, use the latter information to then understand the asymptotic properties of this estimator.

Given the above, the next paragraphs will firstly study the so-called identifiability of a wide class of multivariate latent time series models and then define the asymptotic properties of the MGMWM.

IV-A Identifiability of Multivariate Latent Models

The condition of model identifiability is essential for any method that aims at estimating its corresponding parameters and is consequently necessary to then obtain its asymptotic properties for inference. In brief, model identifiability through the MGMWM can be defined as

[TABLE]

where $\bm{\theta}_{0},\bm{\theta}_{1}\in\bm{\Theta}\subset\mathbb{R}^{p}$. This basically means that there is a one-to-one correspondence between the WCCV and the parameter vector $\bm{\theta}$. This property, although essential, is often assumed in practice since it can be considerably challenging to prove.

In this paper however we present results on the identifiability of a class of multivariate latent processes defined in (3) which can also be used to assess model identifiability on a case-by-case basis. To do so, we need Condition (C1), which allows us to define the WCCV as the sum of the individual WCCV for each underlying model. In this manner, as mentioned earlier, it is possible to define a general class of multivariate latent models for which we can study the condition of identifiability. It must be noted that, in this paper, we use the word “class” to define an overall model which includes all possible subsets of model combinations that compose this overall model (or class). Given this definition, we now study the identifiability of the following general class of multivariate latent models:

1. (M1)

$\bm{X}_{t}\equiv\bm{X}_{1,t}+\bm{X}_{2,t}+\bm{X}_{3,t}+\bm{X}_{4,t}$. 2. (M2)

$\bm{X}_{t}\equiv\bm{X}_{1,t}+\bm{X}_{3,t}+\sum_{k=5}^{K}\bm{X}_{k,t}$.

In the univariate equivalent setting, these two classes of models are commonly used in various disciplines including finance, economics and engineering [11, 22, 29]. More specifically, the class of multivariate latent models (M1) consists of processes that are not second order stationary, although their first-order backward difference is indeed so, while the (M2) class is stationary based on the parameter definitions of the models given earlier. Considering these general classes of multivariate latent models (which therefore include all possible sub-model combinations), the following proposition states injectivity of $\bm{\nu}(\bm{\theta})$ for the first class of models.

Proposition 1:

Under Condition (C1), we have that $\bm{\nu}(\bm{\theta})$ is an injective function of $\bm{\theta}$ for the class of models (M1).

Proof: The proof is based on verifying the conditions for identifiability in [30, Theorem 2]. To do so, we start by defining $\bm{\nu}^{\star}(\bm{\theta})$ which is an appropriately chosen $p$-dimensional subvector of $\bm{\nu}(\bm{\theta})$ (in such a way that both $\bm{\nu}^{\star}(\bm{\theta})$ and $\bm{\theta}$ belong to $\mathbb{R}^{p}$). If the subvector $\bm{\nu}^{\star}(\bm{\theta})$ is identifiable, we have that also the entire vector $\bm{\nu}(\bm{\theta})$ is identifiable. Using log-transformations of the parameters (which is a one-to-one function), we can easily prove Assumptions A and C required by [30] which pertain to the behaviour of the function $\bm{\nu}^{\star}(\bm{\theta})$ (i.e. twice-differentiable and divergent when $\bm{\theta}$ diverges). Finally, we have that the determinant of the Jacobian of $\bm{\nu}^{\star}(\bm{\theta})$ is non-zero which verifies Assumptions B and D in [30] and proves the identifiability of $\bm{\nu}(\bm{\theta})$. A more detailed explanation of this proof can be found in App. LABEL:proof:prop:mod1.injective. ∎

This is an important result since it provides the fundamental basis for this class of models to be identifiable for the MGMWM. However, directly showing the injectivity of $\bm{\nu}(\bm{\theta})$ for the class of models (M2) is not obvious since explicitly verifying whether the determinant of the Jacobian of $\bm{\nu}(\bm{\theta})$ is non-zero is difficult. Nevertheless, we know that the WCCV has a direct relationship with the cross-covariance function through the cross-spectral density and we can therefore understand at least if the cross-covariance function is a one-to-one function thereby providing insight to the identifiability of $\bm{\nu}(\bm{\theta})$. For this purpose, let us denote the cross-covariance function as

[TABLE]

where $C_{l}^{(i,i^{\prime})}(\bm{\theta})\equiv\operatorname*{cov}\left(X_{t}^{(i)},\,X_{t+l}^{(i^{\prime})}\right)$ is the cross-covariance with respect to the directed lag $l$. Notice that $C_{l}^{(i,i^{\prime})}(\bm{\theta})\neq C_{-l}^{(i,i^{\prime})}(\bm{\theta})$ in general.

Considering these quantities, we firstly show the injectivity of $\bm{C}(\bm{\theta})$ for the class of models (M2) and, taking advantage of the one-to-one correspondence between $\bm{C}(\bm{\theta})$ and the cross-spectral density, denoted as $\bm{S}(\bm{\theta})$, study the one-to-one correspondence between $\bm{S}(\bm{\theta})$ and $\bm{\nu}(\bm{\theta})$. In order to do so, we first need to define the following condition where we let $\phi_{k}^{(i)}$ denote the parameter of the $k^{th}$ (T5) process in the $i^{th}$ signal.

1. (C5)

If the $i^{th}$ time series consists in the sum of $(K-4)$ (T5) processes, for $K\geq 5$, then we have that all $\phi_{k}^{(i)}$ are distinct and nonzero.

This condition ensures that the values of the autoregressive parameters of (T5) processes contributing to the $i^{th}$ signal are distinct since otherwise the theoretical cross-covariance function of this signal is not injective. Having given this condition, we can now deliver the following lemma which consists in the first step to proving injectivity of $\bm{\nu}(\bm{\theta})$.

Similar to LABEL:cond.nonflat, we also propose the following condition:

1. (C6)

Let $\bm{C}^{\star}(\bm{\theta})$ be the properly chosen $H(K-4)$-dimensional sub-vector of $\bm{C}^{o}(\bm{\theta})$, the CCV function of a sum of $(K-4)$ number of independent bivariate AR1 processes (see LABEL:proof:lem:ccv.injec in the Supplementary material for an example of choosing $\bm{C}^{\star}(\bm{\theta})$). For every $\bm{C}\in\mathbb{R}^{H(K-4)}$, the equation $\bm{C}^{\star}(\bm{\theta})=\bm{C}$ has countably many (possibly zero) solutions in $\mathbb{R}^{H(K-4)}$.

Lemma 1:

Under conditions (C1) and (C5), we have that the cross-covariance function $\bm{C}(\bm{\theta})$ of (M2) is injective.

Proof: The proof of this lemma is the same as that for Proposition 1 except that the steps are made for an appropriate sub-vector $C^{\star}(\bm{\theta})$ of the autocovariance sequence $C(\bm{\theta})$. For more details see App. LABEL:proofs.ident. ∎

This result is already an important finding since the cross-covariance function is used as a means to estimate these processes via estimators such as the MLE or the GMM. Hence, the identifiability of these models can already easily be proven based on Lemma 1. Nevertheless, as underlined earlier, proving the same property directly for $\bm{\nu}(\bm{\theta})$ is extremely complicated and we will therefore assume an additional condition for this reason. In this sense, let us express $\bm{\nu}(\bm{\theta})\equiv\bm{f}(\bm{S}(\bm{\theta}))$ where $\bm{f}(\cdot)$ is a known vector function e.g. [11].

1. (C6)

$\bm{f}(\cdot)$ is an injective function.

Regarding the above condition, existing results for univariate time series e.g. [31, 32] underline how this condition is not necessarily a strong one to assume in general. Indeed, [32] directly show the identifiability of many latent models including ARMA models and discuss general identifiability for other models supported by the findings of [31] who, in the continuous-time case, shows that this condition is almost always verified when using the Haar WV. With the addition of the latter condition, we can now deliver the following proposition.

Proposition 2:

Under conditions (C1), (C5) and (C6), $\bm{\nu}(\bm{\theta})$ is an injective function of $\bm{\theta}$ for the class of models (M2).

Proof: The proof of this proposition is a straightforward consequence of Lemma 1 since the cross-spectral density function is an injective function of the cross-covariance function $C(\bm{\theta})$ and, using Condition (C6), we prove Proposition 2. ∎

Based on Propositions 1 and Proposition 2, we have provided the results necessary to determine the asymptotic properties of the proposed MGMWM estimator to hold for two very general classes of models (based on the defined conditions) and identifiability can be reasonably assumed for other more specific classes of models if Condition (C6) holds.

IV-B Asymptotic Properties of the MGMWM

Having studied the properties of the WCCV in Sec. III and discussed the identifiability of different classes of multivariate latent models in Sec. IV-A, we can now study the asymptotic properties of the MGMWM. To do so, let us define

[TABLE]

as the objective function of the MGMWM in (9) where $\bm{\Omega}_{T}$ is an estimated positive definite matrix. In a similar way, we can define its theoretical counterpart as

[TABLE]

with $\bm{\theta}_{0}$ and $\bm{\Omega}$ being the true parameter vector and a fixed positive definite weighting matrix respectively.

Considering these definitions, let us further define $||\cdot||_{S}$ as being the matrix spectral norm and “$\overset{p}{\to}$” denoting convergence in probability which allow us to state a set of conditions that are necessary to obtain desirable asymptotic properties for the MGMWM.

1. (C7)

$\bm{\Theta}$ is compact. 2. (C8)

$||\bm{\Omega}_{T}-\bm{\Omega}||_{S}\overset{p}{\to}0$. 3. (C9)

$\bm{\nu}(\bm{\theta})=\bm{\nu}(\bm{\theta}_{0})$, if and only if $\bm{\theta}=\bm{\theta}_{0}$. 4. (C10)

$Q_{0}(\bm{\theta})$ is continuous $\forall\,\bm{\theta}\in\bm{\Theta}$.

Conditions (C7), (C8) and (C10) are standard regularity conditions that, among others, ensure the uniform convergence of $Q_{T}(\bm{\theta})$ while Condition (C9) is related to identifiability which was discussed for different classes of models in the previous section (hence this condition was proven for the general classes of Models (M1) and (M2)). Based on these conditions we can state the following theorem.

Theorem 2:

Under conditions (C2)-(C4) and (C7)-(C10), we have that

[TABLE]

Proof: By replacing the terms in the difference $Q_{T}(\bm{\theta})-Q_{0}(\bm{\theta})$ by their expressions and rearranging terms we can bound $\underset{\bm{\theta}}{\sup}|Q_{T}(\bm{\theta})-Q_{0}(\bm{\theta})|$ by a series of terms that depend on the parameter $\bm{\theta}$ and on the difference between weighting matrices $\bm{\Omega}_{T}$ and $\bm{\Omega}$ which, based on Theorem 1 as well as Conditions (C7) and (C8), represent bounded elements multiplied by a term that goes to zero. Hence, the uniform convergence of $Q_{T}(\bm{\theta})$ is verified and, using [33, Theorem 2.1], we combine this result with Conditions (C9) and (C10) thereby proving Theorem 2. A more detailed proof of this theorem can be found in App. LABEL:proof:thm:gmwm.consis ∎

Having proven consistency of the MGMWM, we now discuss its asymptotic normality. To do so, let us define

[TABLE]

and consider the following additional condition:

1. (C11)

$\bm{\theta}_{0}\mbox{ is an interior point of }\bm{\Theta}$, where $\bm{\Theta}\subset\mathbb{R}^{p}$ is convex. 2. (C12)

$\bm{A}(\bm{\theta}_{0})^{\intercal}\bm{\Omega}\bm{A}(\bm{\theta}_{0})$ is non-singular.

Condition (C11) allows us to make use of the multivariate mean value theorem while Condition (C12) is another regularity condition which allows an asymptotic covariance matrix of $\hat{\bm{\theta}}$ to be defined. Having these conditions, we can deliver the final theoretical result of this paper.

Theorem 3:

Under conditions (C2)-(C4) and (C7)-(C12), we have that

[TABLE]

*where $\bm{\Xi}$ is the asymptotic covariance matrix.

Proof: Having proven consistency in Theorem 2, we can use this result to prove asymptotic normality. Indeed, by definition we have that $\nicefrac{{\partial}}{{\partial\bm{\theta}}}\,Q_{T}(\bm{\theta})|_{\bm{\theta}=\hat{\bm{\theta}}}=\bm{0}$ and, making use of the multivariate mean value theorem, we can re-express and rearrange this last equality to obtain

[TABLE]

where $\bm{B}(\hat{\bm{\theta}})$ and $\bm{A}(\hat{\bm{\theta}},\bm{\theta}_{0})$ are both matrices that converge in probability to $\bm{A}(\bm{\theta}_{0})^{\intercal}$ and $\bm{A}(\bm{\theta}_{0})$ respectively based on Theorem 2. Using Slutsky’s theorem and Theorem 1 we obtain the final result of Theorem 3. A more detailed proof can be found in App. LABEL:proof:thm:gmwm.asy where the explicit form of the asymptotic covariance matrix $\bm{\Xi}$ is also given. ∎

With the above results we have now delivered an important framework to estimate and provide inference on a wide range of multivariate latent models that can considerably extend the modelling possibilities for inertial sensor calibration as well as for a variety of different applied settings. In order to make use of these results in the following sections, some computational approaches are needed to obtain the MGMWM estimator in practice. Various approaches can be used for this purpose, however we recommend using Algorithm 1.

Making use of this approach, the next sections highlight the good properties of the proposed MGMWM framework through some simulation studies and show how these results can be of great importance for considering dependence between sensors which, to date, has not been simple to model despite its relevance for enhanced navigation precision.

V Simulation study

In this section we present a simulation study which replicates the simulation setting considered in [8] who propose an approach to estimate a specific multivariate model composed of a (T1) and (T2) model where multivariate dependence is only present through the latter (a second simulation study is presented in App. LABEL:sim2 where a more complex setting is considered and for which the method in [8] cannot be applied). Within this setting we compare our approach with the estimators proposed in [8] where they consider their Generalized Least Square (GLS) approach using different scales of the Allan Variance (AV). In the latter perspective, for their simulations they consider the GLS discarding either the last two scales of AV (GLS(2)) or the last three scales (GLS(3)). As for the type of simulated process, [8] considered a multivariate process for an array of six gyroscopes constituted by, as mentioned above, the sum of (T1) and (T2) processes (i.e. multivariate WN and RW processes). For the sake of clarity and interpretation of results, in this section we only consider an array of three gyroscopes (corresponding to the first three gyroscopes in their simulation setting). For this reason, the covariance matrix $\bm{\Sigma}$ of the multivariate WN process is given by

[TABLE]

while the covariance matrix $\bm{\Lambda}$ of the multivariate RW process is given by

[TABLE]

The notation used for the parameters of these matrices has been slightly modified for simplicity (i.e. $\sigma_{i}\equiv\sigma^{(i,i)}$ and $\lambda_{ii^{\prime}}\equiv\lambda^{(i,i^{\prime})}$). Based on these definitions, the univariate WN processes are uncorrelated (hence three parameters to characterize their covariance matrix) while the RW processes are all correlated (hence a full symmetric matrix with six parameters to characterize it).

In order to study the finite sample performance of the estimators we make use of the empirical Mean Squared Error (MSE) which takes into account bias and variance of an estimator. This measure is based on 500 simulated replicates where the estimators are computed on each multivariate simulated process of length $T=10^{5}$. The empirical MSE of the estimators are shown in Fig. 1 where, for the GLS estimators, only the results for GLS(2) were represented since they were visually indistinguishable from the results of GLS(3).

Considering that the shaded areas represent the empirical confidence intervals of the MSE, it is evident from the plots that the MGMWM outperforms the GLS(2) (as well as the omitted GLS(3)) for the estimation of all the parameters considered in the simulation study (the corresponding boxplots of the results can be found in App. LABEL:sim2). As mentioned at the beginning of this section, another simulation study was performed on more complex composite processes for which other methods (including the GLS approach) are either unavailable or computationally unfeasible (these simulations can be found in App. LABEL:sim2).

VI Case study

In this section we apply the MGMWM to perform a multisensor calibration procedure on some real data. Sensor calibration is an extremely important procedure to ensure navigation performance, especially for low-cost IMUs whose employment is rapidly increasing but suffers from more complex error measurements that need to be accounted for. In general, this procedure is performed independently on each sensor composing the IMU which implicitly assumes that the error measurements of each sensor are independent of each other. The latter unrealistic assumption has nevertheless been necessary to deal with the complex error signals of these sensors which have been extremely difficult to characterize also at an individual level. Since the GMWM has allowed the problems to be greatly reduced in modelling the individual error measurements of each sensor, the MGMWM is an extension that allows complex multivariate dependencies between sensors to be taken into account and consequently make the modelling procedure more realistic and improve navigation precision. Moreover, this approach is starting to be taken increasingly into account as underlined in [8] where the construction of virtual gyroscopes based on sensor interdependence can considerably reduce errors in different navigation measurements.

For illustrative purposes, in this section we only consider two stochastic error measurements coming from the accelerometer X-axes of an array of two same-model IMUs (namely XSens MTi-G IMUs) sampled at 250Hz over roughly 1 hour with consequent sample size $T=873^{\prime}684$ (this procedure can of course take into account more sensors).

The diagonal plots of Fig. 2 (top-left and bottom-right) show the log-log plots of the empirical WV for each error signal of these sensors where the solid dotted lines represent the estimated WV (surrounded by a shaded area representing the confidence intervals). Simply looking at the lines in the diagonal plots we can see how they appear to have similar behaviour with different patterns at larger scales (as expected) where there is more variability. In order to select the appropriate multivariate model we first focus on the individual signals making use of visual matching of the WV plots and selection criteria based on the WV (see [guerrier2015automatic]). Based on this, we find that the sum of two (T5) models (AR(1)) and a (T2) model (RW) fits both signals well as seen by how close the dashed orange line follows the empirical WV. The contribution of the individual models is represented by the individual solid lines below the empirical WV. In practice, the first AR(1) process would usually be replaced by a (T1) model (WN) since the slope of the latter is close to the one observed at the first scales in the plots and the WN is more straightforward to integrate within a navigation filter. In addition, within the latter filter an AR(1) process with a very small correlation time is unlikely to be observable. However, for the purpose of this study we consider the models that better fit the observed error signals which obviously does not exclude the possibility of considering other models in practice. Having defined the individual models, we can now understand if any of these individual models contributes to the WCCV represented by the solid dotted line in the off-diagonal plots of Fig. 2 (top-right and bottom-left). Estimating all possible dependencies, we found that one of the (T5) models allowed the WCCV to be well described as shown by the orange dashed line representing the WCCV implied by the estimated multivariate model.

This case study highlights how the MGMWM allows to be increased the ease with which complex multivariate composite processes can be modelled, improving in terms of statistical performance and flexibility over existing approaches for multivariate sensor calibration. Indeed, its generality and computational feasibility can greatly contribute, for example, to facilitating and improving the creation of virtual sensors that can considerably enhance navigation performance. Indeed, this approach has already been used in [34] to construct an optimal virtual gyroscope with good results both in simulated and applied settings.

IMU4: RW(Q=1.756252e-11); 1AR(phi=1.300635e-01, Z=8.142854e-05); 2AR(phi=9.989909e-01, Z=1.612509e-10). IMU6: RW(Q=5.015933e-12); 1AR(phi=7.466659e-02, Z=1.255179e-04); 2AR(phi=9.999121e-01, Z=2.075570e-10). Cross: 1AR(phi1 = 0.1300635, phi2 = 0.07466659, Z.cross = -4.603401e-05)

VII Conclusions

This paper has presented a new approach to modelling a wide class of (composite) multivariate time series in order to take into account complex dependence structures between individual signals, thereby providing a contribution to improve multivariate time series modelling in general. Building upon the GMWM framework, we studied the WCCV and extended the asymptotic properties of its estimator to the more general multivariate case under weaker conditions compared to existing work. Based on these properties, we extended the GMWM framework by proposing the MGMWM which takes into account the WCCV to describe the dependence between individual time series. Benefitting from different identifiability results, we delivered the asymptotic properties of the MGMWM estimators under relatively weak conditions and highlighted its improved performance with respect to existing approaches. The case study underlined how the MGMWM can greatly support the calibration of the multivariate stochastic error signals coming from IMUs and consequently contribute to improving navigation performance by, for example, testing for dependence between sensors, modelling inter-axis dependence within an IMU or across sensors in an array of redundant sensors and, based on the latter, constructing optimal virtual sensors.