Multivariate Signal Modelling
with Applications to Inertial Sensor Calibration
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 which are challenging to estimate. 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, albeit challenging, 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 as well as the dependence between them by studying the quantity called Wavelet CrossCovariance and using it to extend the application of the Generalized Method of Wavelet Moments to this setting. This new method can be used in many other settings for multivariate time series modelling, especially in cases where the dependence among signals may be hard to detect since it can be based on a shared underlying model that has a marginal contribution to the processes’ overall variance. 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 MTiG IMUs.
pinegreenrgb0.0, 0.47, 0.44
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 [see 24, 5], in economic and financial research [see 3, 26, 17, 21] as well as in different fields of engineering [see 7, 29]. 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 [see for example 29]. 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 MaximumLikelihood Estimator (MLE) cannot be directly applied since the latent stochastic components are unobservable and their marginal distribution is usually unknown. A feasible likelihoodbased method is the ExpectationMaximization (EM) algorithm [see 6] whose implementation is nevertheless limited due to the complexity and, more importantly, the possible nonconvexity of its associated optimization problems therefore severely limiting the use of its statistical properties in practice [see 1]. To avoid the problems of the likelihoodbased 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 [see 8]. Although, the AVLR is a computationally feasible method for estimating many complex models, its statistical properties are not optimal as discussed in [13] 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 nonoptimal 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 [see e.g. 2]. 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 AutoRegression (VAR) model which, in some sense, generalizes the AutoRegression (AR) model in order to allow lagged observations from some signals to influence current observations of others [see e.g. 14]. This class of models has been extended to the more general class of Vector AutoRegression Moving Average (VARMA) models which also allow for correlation between innovation processes in each time series [see 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 [see 18, 19]. 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 [see e.g. 4]. 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 [see 12, 15, 17]. We refer to these models as multivariate latent models which are particularly pertinent 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 [see 12]. To do so the GMWM takes advantage of a quantity called the Wavelet Variance (WV) which is the variance of the wavelet coefficients issued from a wavelet decomposition of a signal [see 23]. 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 CrossCovariance (WCCV), which was proposed by Whitcher et al. [31] 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 , where . 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:
where and are the maximum lags of dependence for the AR and MA processes respectively among the individual time series composing the mutlivariate process, and are the coefficient matrices for the and lag and
(1) 
are independently distributed innovation vectors with multivariate distribution , expectation and covariance matrix .
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 firstorder autoregressive and white noise models [see e.g. 9]. 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 make use of them 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 latent process as
(2) 
where is an dimensional vector and is a dimensional vector with being the index for the type of univariate/multivariate model underlying some or all of the latent processes and is an matrix with either or as its elements that we define as
Therefore, we have a set of 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 now introduce the processes that are each characterized by one of the multivariate models that we consider in the context of this paper and that constitute the basic models used for inertial sensor calibration [see e.g. 28]:

White Noise (WN) which corresponds to the independently distributed innovation vectors defined in (1). We denote the process characterized by this model as , hence .

Random Walk (RW) which we denote as and is defined as
where .

Quantization Noise (QN) (or rounding error, see e.g. Papoulis and Pillai [22]) which we denote as . For this process we do not consider a multivariate model and therefore each univariate process is characterized by the parameter .

Drift (DR) which we denote as and, for each univariate process, is defined as
where .

FirstOrder AutoRegressive Noise (AR1) which we denote as , , where , and defined as
where is a diagonal matrix with diagonal elements , , for and .
The covariance matrices , and defined above are positive definite matrices with repsectively, with their respective elements being , and . As can be seen from the above definitions, processes 3 and 4 don’t 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 3 can be interpreted as a form of rounding error which is intuitively independent from the rounding error made on another signal while the 4 process is a nonstochastic 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 wellapproximate 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:

is independent of , .
This condition implies that each model that composes a latent model (univariate or multivariate) is independent from the others (e.g. a 1 process is independent from a 5 process and two different 5 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. IVB). 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 casebycase 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 be a fixed integer representing the last level of wavelet decomposition of the univariate time series . This allows us to define the wavelet coefficients for level as:
where is the level wavelet filter with length . In this paper we choose to use the Haar wavelet filter whose length is consequently .
Based on the wavelet coefficients, we can define the WV as:
(3) 
where
is the parameter vector defining the latent model underlying the univariate process (of dimension ). Indeed, for each time series model there is a corresponding theoretical WV which can be used for estimation purposes in a methodofmoments fashion (as done by the GMWM). In a similar way, we can also define the WCCV [see 31] at the level as follows:
(4) 
where represents the lag in time between observations and represents the parameter vector defining the dependence between signals (of dimension ).
Having defined these theoretical quantities, we can now define their empirical counterparts that are estimated on the observed signals and , where . While the unbiased estimator of WV is defined and studied in [23] and [25], here we define the unbiased estimator of WCCV [see 31] for as:
(5) 
where is the number of wavelet coefficients generated at level . In this paper we limit ourselves to studying the WCCV at lag 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 since they are essential to deliver inference when using the MGMWM which is presented in Sec. IV. In [30] the asymptotic normality of was stated for any pair and any level . However, the conditions for the latter results to hold assume that the process
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: and . Based on this notation, the WCCV vector is defined as:
with the corresponding vector of estimated WCCV defined as:
In addition, let us also define the first order difference of the signal as (i.e. ). The reason for considering this process lies in the fact that the wavelet coefficients can be represented as particular linear combinations of the process , i.e.
(6) 
where and represents a specific filter vector. In order to define the first condition, we also define
as being an valued measurable function as well as the filtration where are i.i.d. random variables. We can now state the first condition needed for the asymptotic properties of the WCCV vector.

We assume the multivariate process is well defined and can be represented as
This condition basically implies that the multivariate process 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 [see e.g. 33].
In order to state the final conditions, we need to deliver additional definitions, starting from where , where also is an i.i.d. random variable. Based on this definition, we can notice that for we have that . Finally, we define , for , to deliver the other required conditions.

.

.
Both conditions require the fourth moment of the difference process (or some function of it) to be finite (for all ). More specifically, while Condition 1 is straightforward to interpret in the latter sense, Condition 2 ensures that the cumulative impact of on future values of the process is finite which allows us to interpret it as a shortrange dependence condition [see 32]. We can now deliver the main result of this section and, for this purpose, we further define
and the projection operator
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 ) should be arbitrarily close to the unconditional expectation. Using these last definitions, we can finally state the asymptotic normality of the vector .
Proof.
We can split the quantity as follows:
where is an alternative estimator of the WCCV whose elements are given by
for all and . Using the latter estimator and Conditions 1 to 2, we can directly apply the results of Theorem 7 in [32] to show asymptotic normality of with covariance matrix . Moreover, with and being fixed, we have that and therefore, using Slutsky’s theorem, we have that asymptotically is equivalent in distribution to 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 [see 12]. Indeed the latter approach considers a univariate process and, assuming that this process is generated from a parametric model known up to the parameter vector , the GMWM aims at estimating this parameter vector through the following minimization problem:
(7) 
where and are the vectors of estimated and theoretical WV respectively, is the parameter space for and 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 [see 12]. The consistency and asymptotic normality of were shown in Guerrier et al. [12] for any suitable choice of , 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 . 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 that contain information on the latent dependence between signals and . For this reason, it is possible to build the MGMWM as an extension of the GMWM as follows:
(8) 
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 (8) 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 (7) (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 . 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 socalled identifiability of a wide class of multivariate latent time series models and then define the asymptotic properties of the MGMWM.
Iva 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
where . This basically means that there is a onetoone correspondence between the WCCV and the parameter vector . 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 (2) which can also be used to assess model identifiability on a casebycase basis. To do so, we need Condition 1, 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:

.

.
In the univariate equivalent setting, these two classes of models are commonly used in various disciplines including finance, economics and engineering [see 8, 9, 27]. More specifically, the class of multivariate latent models 1 consists of processes that are not second order stationary, although their firstorder backward difference is indeed so, while the 2 class is stationary based on the parameter definitions of the models given earlier. Considering these general classes of mutlivariate latent models (which therefore include all possible submodel combinations), the following proposition states injectivity of for the first class of models.
Proposition 1:
Proof.
The proof is based on verifying the conditions for identifiability in Theorem 2 of [16]. To do so, we start by defining which is an appropriately chosen dimensional subvector of (in such a way that both and belong to ). If the subvector is identifiable, we have that also the entire vector is identifiable. Using logtransformations of the parameters (which is a onetoone function), we can easily prove Assumptions A and C required by [16] which pertain to the behaviour of the function (i.e. twicedifferentiable and divergent when diverges). Finally, we have that the determinant of the Jacobian of is nonzero which verifies Assumptions B and D in [16] and proves the identifiability of . 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 for the class of models 2 is not obvious since explicitly verifying whether the determinant of the Jacobian of is nonzero is difficult. Nevertheless, we know that the WCCV has a direct relationship with the crosscovariance function through the crossspectral density and we can therefore understand at least if the crosscovariance function is a onetoone function thereby providing insight to the identifiability of . For this purpose, let us denote the crosscovariance function as
where is the crosscovariance with respect to the directed lag . Notice that in general.
Considering these quantities, we firstly show the injectivity of for the class of models 2 and, taking advantage of the onetoone correspondence between and , study the onetoone correspondence between and . In order to do so, we first need to define the following condition where we let denote the parameter of the 5 process in the signal.

If the time series consists in the sum of 5 processes, for , then we have that all are distinct and nonzero.
This condition ensures that the values of the autoregressive parameters of 5 processes contributing to the signal are distinct since otherwise the theoretical crosscovariance 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 .
Proof.
The proof of this lemma is the same as that for Proposition 1 except that the steps are made for the an appropriate subvector of the autocovariance sequence . For more details see App. LABEL:proofs.ident. ∎
This result is already an important finding since the crosscovariance 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 is extremely complicated and we will therefore assume an additional condition for this reason. In this sense, let us express where is a known vector function [see e.g. 8].

is an injective function.
Regarding the above condition, existing results for univariate time series [see e.g. 10, 11] underline how this condition is not necessarily a strong one to assume in general. With the addition of the latter condition, we can now deliver the following proposition.
Proof.
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 6 holds.
IvB 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. IVA, we can now study the asymptotic properties of the MGMWM. To do so, let us define
as the objective function of the MGMWM in (8) where is an estimated positive definite matrix. In a similar way, we can define its theoretical counterpart as
with and being the true parameter vector and a fixed positive definite weighting matrix respectively.
Considering these definitions, let us define as being the matrix spectral norm which allows us to state a set of conditions that are necessary to obtain desirable asymptotic properties for the MGMWM.

is compact.

.

, if and only if .

is continuous .
Conditions 7, 8 and 10 are standard regularity conditions that, among others, ensure the uniform convergence of while Condition 9 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 1 and 2). Based on these conditions we can state the following theorem.
Proof.
By replacing the terms in the difference by their expressions and rearranging terms we can bound by a series of terms that depend on the parameter and on the difference between weighting matrices and which, based on Theorem 1 as well as Conditions 7 and 8, represent bounded elements multiplied by a term that goes to zero. Hence, the uniform convergence of is verified and, using Theorem 2.1 in [20], we combine this result with Conditions 9 and 10 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
and consider the following additional condition:

, where is convex.

is nonsingular.
Condition 1 allows us to make use of the multivariate mean value theorem to perform a MacLaurin expansion while Condition 2 is another regularity condition which allows to define an asymptotic covariance matrix of . Having these conditions, we can deliver the final theoretical result of this paper.
Proof.
Having proven consistency in Theorem 2, we can use this result to prove asymptotic normality. Indeed, by definition we have that and, making use of the multivariate mean value theorem, we can reexpress and rearrange this last equality to obtain
where and are both matrices that converge in probability to and 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 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 in a variety of different applied settings. 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 Vaccaro and Zaki [29] who propose an approach to estimate a specific multivariate model composed of a 1 and 2 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 Vaccaro and Zaki [29] cannot be applied). Within this setting we compare our approach with the estimators proposed in Vaccaro and Zaki [29] 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, Vaccaro and Zaki [29] considered a multivariate process for an array of six gyroscopes constituted by, as mentioned above, the sum of 1 and 2 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 of the multivariate WN process is given by
while the covariance matrix of the multivariate RW process is given by
The notation used for the parameters of these matrices has been slightly modified for simplicity (i.e. and ). Based on these definitions, the univariate WN processes are uncorrelated (hence three parameters to characterise their covariance matrix) while the RW processes are all correlated (hence a full symmetric matrix with six parameters to characterise 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 . The empricial 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 lowcost 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 to greatly reduce the problems in modelling the individual error measurements of each sensor, the MGMWM is an extension that allows to take into account complex multivariate dependencies between sensors 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 Vaccaro and Zaki [29] 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 Xaxes of an array of two samemodel IMUs (namely XSens MTiG IMUs) sampled at 250Hz over roughly 1 hour with consequent sample size (this procedure can of course take into account more sensors).
The diagonal plots of Fig. 2 (topleft and bottomright) show the loglog 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. After testing different models we find that the sum of two 5 models (AR(1)) and a 2 model (RW) fits both 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 1 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 offdiagonal plots of Fig. 2 (topright and bottomleft). Estimating all possible dependencies, we found that one of the 5 models allowed to well describe the WCCV 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 increase 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 Zhang et al. [34] to construct an optimal virtual gyroscope with good results both in simulated and applied settings.
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 interaxis dependence within an IMU or across sensors in an array of redundant sensors and, based on the latter, constructing optimal virtual sensors.
References
 Anandkumar et al. [2012] Animashree Anandkumar, Daniel J Hsu, and Sham M Kakade. A method of moments for mixture models and hidden markov models. In Conference on Learning Theory, volume 1, page 4, 2012.
 Bayard and Ploen [2002] DS Bayard and SR Ploen. Foundations of virtual gyroscope synthesis. Jet Propulsion Laboratory Internal Document, JPLD21656, 2002.
 Campbell et al. [1997] John Y Campbell, Andrew WenChuan Lo, and Archie Craig MacKinlay. The econometrics of financial markets. Princeton University press, 1997.
 Chang et al. [2018] Jinyuan Chang, Bin Guo, Qiwei Yao, et al. Principal component analysis for secondorder stationary vector time series. The Annals of Statistics, 46(5):2094–2124, 2018.
 Cho [2016] Hyunkeun Cho. The analysis of multivariate longitudinal data using multivariate marginal models. Journal of Multivariate Analysis, 143:481–491, 2016.
 Dempster et al. [1977] Arthur P Dempster, Nan M Laird, and Donald B Rubin. Maximum likelihood from incomplete data via the em algorithm. Journal of the Royal Statistical Society: Series B (Methodological), 39(1):1–22, 1977.
 Donoho [2006] David L Donoho. Compressed sensing. IEEE Transactions on information theory, 52(4):1289–1306, 2006.
 ElSheimy et al. [2008] Naser ElSheimy, Haiying Hou, and Xiaoji Niu. Analysis and modeling of inertial sensors using allan variance. IEEE Transactions on instrumentation and measurement, 57(1):140–149, 2008.
 Granger and Morris [1976] Clive WJ Granger and Michael J Morris. Time series modelling and interpretation. Journal of the Royal Statistical Society. Series A (General), pages 246–257, 1976.
 Greenhall [1998] Charles A Greenhall. Spectral ambiguity of allan variance. IEEE Transactions on Instrumentation and Measurement, 47(3):623–627, 1998.
 Guerrier and Molinari [2016] Stéphane Guerrier and Roberto Molinari. On the identifiability of latent models for dependent data. arXiv preprint arXiv:1607.05884, 2016.
 Guerrier et al. [2013] Stéphane Guerrier, Jan Skaloud, Yannick Stebler, and MariaPia VictoriaFeser. Waveletvariancebased estimation for composite stochastic processes. Journal of the American Statistical Association, 108(503):1021–1030, 2013.
 Guerrier et al. [2016] Stéphane Guerrier, Roberto Molinari, and Yannick Stebler. Theoretical limitations of allan variancebased regression for time series model estimation. IEEE Signal Processing Letters, 23(5):597–601, 2016.
 Hamilton [1994] James Douglas Hamilton. Time series analysis, volume 2. Princeton university press Princeton, 1994.
 Harvey [1990] Andrew C Harvey. Forecasting, structural time series models and the Kalman filter. Cambridge university press, 1990.
 Komunjer [2012] Ivana Komunjer. Global identification in nonlinear models with moment restrictions. Econometric Theory, 28(4):719–729, 2012.
 Koopman and Durbin [2000] Siem J Koopman and James Durbin. Fast filtering and smoothing for multivariate state space models. Journal of Time Series Analysis, 21(3):281–296, 2000.
 Lütkepohl [2006] Helmut Lütkepohl. Forecasting with varma models. Handbook of economic forecasting, 1:287–325, 2006.
 Metaxoglou and Smith [2007] Konstantinos Metaxoglou and Aaron Smith. Maximum likelihood estimation of varma models using a statespace em algorithm. Journal of Time Series Analysis, 28(5):666–685, 2007.
 Newey and McFadden [1994] Whitney K Newey and Daniel McFadden. Large sample estimation and hypothesis testing. Handbook of econometrics, 4:2111–2245, 1994.
 Pan and Yao [2008] Jiazhu Pan and Qiwei Yao. Modelling multiple time series via common factors. Biometrika, 95(2):365–379, 2008.
 Papoulis and Pillai [2002] Athanasios Papoulis and S Unnikrishna Pillai. Probability, random variables, and stochastic processes. Tata McGrawHill Education, 2002.
 Percival [1995] Donald B Percival. On estimation of the wavelet variance. Biometrika, 82(3):619–631, 1995.
 Roy and Lin [2000] Jason Roy and Xihong Lin. Latent variable models for longitudinal data with multiple continuous outcomes. Biometrics, 56(4):1047–1054, 2000.
 Serroukh et al. [2000] Abdeslam Serroukh, Andrew T Walden, and Donald B Percival. Statistical properties and uses of the wavelet variance estimator for the scale analysis of time series. Journal of the American Statistical Association, 95(449):184–196, 2000.
 Sims [1980] Christopher A Sims. Macroeconomics and reality. Econometrica: Journal of the Econometric Society, pages 1–48, 1980.
 Skrondal and RabeHesketh [2007] Anders Skrondal and Sophia RabeHesketh. Latent variable modelling: a survey. Scandinavian Journal of Statistics, 34(4):712–745, 2007.
 Titterton et al. [2004] David Titterton, John L Weston, and John Weston. Strapdown inertial navigation technology, volume 17. IET, 2004.
 Vaccaro and Zaki [2017] Richard J Vaccaro and Ahmed S Zaki. Reduceddrift virtual gyro from an array of lowcost gyros. Sensors, 17(2):352, 2017.
 Whitcher et al. [1999] Brandon Whitcher, Peter Guttorp, and Donald B Percival. Mathematical background for wavelet estimators of crosscovariance and crosscorrelation. 1999.
 Whitcher et al. [2000] Brandon Whitcher, Peter Guttorp, and Donald B Percival. Wavelet analysis of covariance with application to atmospheric time series. Journal of Geophysical Research, 105(D11):941–962, 2000.
 Wu [2011] Wei Biao Wu. Asymptotic theory for stationary processes. Statistics and its Interface, 4(2):207–226, 2011.
 Zhang et al. [2017] Danna Zhang, Wei Biao Wu, et al. Gaussian approximation for high dimensional time series. The Annals of Statistics, 45(5):1895–1919, 2017.
 Zhang et al. [2018] Yuming Zhang, Haotian Xu, Ahmed Radi, Roberto Molinari, Stéphane Guerrier, Mucyo Karemera, and Naser ElSheimy. An optimal virtual inertial sensor framework using wavelet cross covariance. In 2018 IEEE/ION Position, Location and Navigation Symposium (PLANS), pages 1342–1350. IEEE, 2018.
See pages  of suppmat.pdf