Wavelet-Variance-Based Estimation for Composite Stochastic Processes.

This article presents a new estimation method for the parameters of a time series model. We consider here composite Gaussian processes that are the sum of independent Gaussian processes which, in turn, explain an important aspect of the time series, as is the case in engineering and natural sciences. The proposed estimation method offers an alternative to classical estimation based on the likelihood, that is straightforward to implement and often the only feasible estimation method with complex models. The estimator furnishes results as the optimization of a criterion based on a standardized distance between the sample wavelet variances (WV) estimates and the model-based WV. Indeed, the WV provides a decomposition of the variance process through different scales, so that they contain the information about different features of the stochastic model. We derive the asymptotic properties of the proposed estimator for inference and perform a simulation study to compare our estimator to the MLE and the LSE with different models. We also set sufficient conditions on composite models for our estimator to be consistent, that are easy to verify. We use the new estimator to estimate the stochastic error's parameters of the sum of three first order Gauss-Markov processes by means of a sample of over 800,000 issued from gyroscopes that compose inertial navigation systems. Supplementary materials for this article are available online.

INTRODUCTION

Let F, ∈ Θ ⊆ ℝ be the model associated to the univariate Gaussian time series {Y; t ∈ ℤ} that is stationary or nonstationary but with stationary backward differences of order d, and let {y, t = 1, …, T} be the corresponding observed outcome. In this article, we propose a new estimation method for the parameter's vector that is based on the matching between the empirical and model-based wavelet variances (WV). We also set the conditions on F under which our estimator is consistent. In particular, we consider for the empirical WV the Maximal Overlap Discrete Wavelet Transform (MODWT) WV estimator (Greenhall 1991; Percival and Guttorp 1994) based on Haar wavelet filters for which d = 1. In this case, we show that for a model made of the sum of independent Gaussian white-noise (WN), drift, quantization noise (QN), random walk (RW), and a finite number of autoregressive models of order 1 (AR(1)), our estimator is consistent (see Corollaries 2 and 3 in the supplementary materials, Section F). This model encompasses most of the stochastic models used in engineering and natural sciences applications.

The processes can be represented without loss of generality, by a state-space model of the form

with measurements

where x is a q × 1 system state vector at time t, Φ is a q × q coefficient or state-transition matrix from t to t + 1, w is a q × 1 multivariate WN vector, that is, w (0, Q), u is a q × 1 deterministic input vector, y is the one-dimensional output variable, h is a 1 × q design vector which maps the state vector into the output y, and v is a one-dimensional noise such that v (0, σ2WN).

The books by Harvey (1989) and Durbin and Koopman (2001) contain extensive accounts of state-space models and their applications. For linear and/or Gaussian state-space models, the Maximum Likelihood Estimator (MLE) is a natural choice for the estimation of the model's parameters. The Kalman filter, MCMC, or other simulation smoother, together with the specification of a prior to start the latent process, are typically used for computing predictors of the state-variables and one-step-ahead predictors of the observations which are then used in an Expectation-Maximization (EM) algorithm of Dempster, Laird, and Rubing (1977) to compute the MLE (see also, e.g., Shumway and Stoffer 1982). The MLE is based on the density f(y) = ∂ F(y)/∂ y, with the p × 1 vector of parameters containing the unknown elements of Φ, Q, h and possibly also u together with σ2WN. The maximization step can be very complex and finding the MLE is not always a simple task. Moreover, even if new smoothers are regularly proposed (see, e.g., Durbin and Koopman 2002; McCausland, Miller, and Pelletier 2011) that improve the computational efficiency, the task becomes even more challenging when the observed process size is large and the model is, for example, a mixture of three first-order Gauss–Markov (GM) random processes (a reparameterization of AR(1) processes) with different parameters, a model that is used in the case study Section 5. We also found in a simulation exercise in Section 4 that our estimator, unlike the classical MLE or the least squares estimator (LSE), is importantly less biased when the roots of an AR(1) process lie near the unit circle. However, a formal and more general treatment of the unit root problem is left for future research.

An important useful tool for understanding underlying features of a process is the Power Spectral Density (PSD), which for stationary processes, is given by for |f| ≤ f = 1/(2Δt) (f is the Nyquist frequency) with C(τ) = cov [Y, Y] being the Auto-Covariance Function (ACF) of {Y} such that . The WV are related to the PSD since they decompose the variance process var as thoroughly described in Percival and Walden (2000) (see also Section 2). The WV can be estimated from the sample {y, t = 1, …, T} using, for example, the maximal-overlap discrete wavelet transform estimator defined in Serroukh, Walden, and Percival (2000) (see also Percival 1995 and Equation (6) next), and for the matching of these empirical estimates to the WV implied by the model, we propose an optimization criterion based on a standardized distance between both WV estimates, as is done, for example, with the Generalized Method of Moments (GMM) (Hansen 1982). We call the new estimator the Generalized Method of Wavelet Moments (GMWM) estimator. We provide, among others, sufficient conditions on the model F for the GMWM estimator to be consistent and asymptotically normal. Moreover, it can be computed using simulations as is done with indirect inference (Gourieroux, Monfort, and Renault 1993; Smith 1993; Gallant and Tauchen 1996), so that it is very straightforward to implement in practice.

While in this article we use wavelets (through their variances) for parametric estimation, wavelets have been used in many different statistical problems mainly for nonparametric estimation. To cite only a few, wavelets have been used in nonparametric density estimation (see, e.g., Doukhan and Leon 1990; Kerkyacharian and Picard 1992; Walter 1992; Donoho et al. 1996), inverse problems (see, e.g., Donoho 1995), nonparametric regression or curve fitting (see, e.g., Donoho and Johnstone 1994 and Donoho et al. 1995) also with correlated errors (see, e.g., Johnstone and Silverman 1997), or with locally stationary noise (see, e.g., Von Sachs and Macgibbon 2000).

This article is organized as follows. In Section 2, we introduce the WV and explain how it is used mainly in the engineering literature. In Section 3, we present the estimator based on the WV and derive its statistical properties. A simulation study is then presented in Section 4 that compares the finite sample performance of the GMWM estimator to the MLE and the LSE. In Section 5, we apply the new methodology to estimate parameters for the stochastic error model in inertial sensors with a real dataset of size 833,685.

THE WAVELET VARIANCE

Basically, as pointed out by Percival and Guttorp (1994), the WV can be interpreted as the variance of a process after it has been subject to an approximate bandpass filter. WV can be built using wavelet coefficients issued from a modified Discrete Wavelet Transform (DWT) (Mallat 1999; Percival and Walden 2000) called the Maximal Overlap DWT (MODWT); see Greenhall (1991), Percival and Guttorp (1994). The wavelet coefficients are built using wavelet filters {, J which for j = 1 and for the MODWT satisfy

where for l < 0 and l ≥ L1, L1 is the length of , and m is a nonzero integer. Let also be the transfer function of . To obtain the jth-level wavelet filters {j, t} of length L = (2 − 1)(L1 − 1) + 1, one computes the inverse discrete Fourier Transform of

The MODWT filter is actually a rescaled version of the DWT filter h, that is, . Filtering an infinite sequence {Y; t ∈ ℤ} using the wavelet filters {} yields the MODWT wavelet coefficients

We define the WV at dyadic scales τ = 2, as the variances of {}, that is

Notice that the WV are assumed not to depend on time. The condition for this property to hold is that the integration order d for the series {Y} has to be stationary such that d ≤ L1/2 and {} is based on a Daubechies (Daubechies 1992) wavelet filter; see Percival and Walden (2000), chapter 8. This is because Daubechies wavelet filters of width L1 contain an embedded backward difference filter of order L1/2. In such a case, the series of wavelet coefficients {} is stationary with PSD S(f) = |(f)|2S(f), | · | denoting the modulus. Hence (see Serroukh, Walden, and Percival 2000)

There is therefore an implicit link between the WV and the parameters of the data generating model F. We exploit this connection when defining an estimator for , namely by matching a sample estimate of the WV ν2(τ) together with the model-based expression of the WV given by the left side of Equation (5). For WV based on Haar wavelet filters (see Equation (9) next) and for the WN, RW, Drift, QN, AR(1), ARMA(1,1), and ARIMA(0,1,1) models, the integral in Equation (5) is solved and given in the supplementary materials (see Section A), based on the results of Zhang (2008). WV for other models can be computed using the same methodology.

For a finite (observed) process {y; t = 1, …, T}, the MODWT is given by

with and M = T – L + 1, is a consistent estimator for ν2(τ); see Serroukh, Walden, and Percival (2000) who also show that under suitable conditions (see also Theorem 3 in the supplementary materials), is asymptotically normal with mean 0 and variance

Equation (7) can be estimated by means of

and the asymptotic properties of Equation (8) are given in Percival and Walden (2000), p. 312.

A particular choice for the wavelet filter is the Haar wavelet filter whose first DWT filter (j = 1) is

with length L1 = 2. If the process is stationary with backward differences of order d > 1, then other wavelet filters such as Daubechies wavelet filters can be used (Daubechies 1992). When the WV is evaluated with Haar wavelet filters, it is actually equal to half the Allan variance (AV) (Allan 1966). In 1998, the IEEE Standard 1293 (1998) introduced the AV technique as a noise identification method which can be used for determining the characteristics of underlying stochastic processes affecting signals. AV and WV are often used in engineering disciplines and physical sciences as a graphical approach for model building purposes. In the supplementary materials (see Section B), we provide an example of the graphical use of WV. For a more detailed description, see, for example, Percival and Walden (2000).

When a graphical representation is suitable, with Haar WV, the parameters of stochastic processes are sometimes estimated by means of linear regression on preidentified linear regions in a log-log plot of scale τ versus WV ν2(τ). For example, this approach has been used for over 30 years as a standard routine measure of frequency stability in lasers (Fukuda, Tachikawa, and Kinoshita 2003) or atomic clocks (Allan 1987). More recently, the WV has also been used with optical sensors (Kebabian, Herndon, and Freedman 2005), various types of gas monitoring spectrometers (Werle, Mücke, and Slemr 1993; Bowling et al. 2003; Skrínskı et al. 2009), sonic anemometer-thermometers (Loescher et al. 2005), inertial sensors (El-Sheimy, Hou, and Niu 2008; Guerrier 2009), and radio-astronomical instrumentation (Schieder and Kramer 2001). The WV was also used in Percival and Guttorp (1994) to analyze vertical ocean shear measurements. In Fadel et al. (2004), it was employed to study the variability in heart-beat intervals. In Whitcher (2004), discrete wavelet packet transforms are used to estimate one of the parameters of a seasonal long memory process for the analysis of atmospheric and economic time series. In Gebber, Orer, and Barman (2006), WV is exploited to study nerve activities. WV has also been applied in Earth orientation metrology in Gambis (2002) and with other types of (geo)physical data. However, the linear regression on identified linear regions of the WV plots provides reasonably estimated parameters only for a limited number of processes and is often biased (Stebler et al. 2011). In our simulated example presented in Figure 4, the graphical estimation of first-order GM process parameters mixed with other processes like a WN is not suitable.

In the following section, we propose instead a criterion based on a standardized distance between sample and model-based WV that provides consistent estimators of the model's parameters for a wide range of models.

GMWM ESTIMATOR

We propose to estimate the model's parameters using an estimator which combines on the one hand the WV and on the other hand the Generalized Least Squares (GLS) principle, using the relationship given in Equation (5). More precisely, we propose to find such that the WV implied by the model, say (), matches the empirical WV, say , and solves the following GLS optimization problem

in which Ω, a positive-definite weighting matrix, is chosen in a suitable manner (see next). Equation (10) defines the GMWM estimator. (·) = [ϕ(·)] is a binding function between and ν2 = [ν2(τ)] such that () = 2, and and are two estimators.

In Theorem 3 given in the supplementary materials (see Section C), we use the results of Giraitis and Taqqu (1997) on limit theorems for bivariate Appell polynomials to set the conditions for to be asymptotically multivariate normal. In particular, S(f) ≤ C|f|−α, f ∈ [−1/2, 1/2], α < 1/2. This result generalizes the results of Serroukh, Walden, and Percival (2000) to the multivariate case, while Serroukh and Walden (2000a, 2000b) generalize the results of Serroukh, Walden, and Percival (2000) to wavelet covariances.

In Theorem 1 next, we state the conditions on , , and () (hence on F) under which is a consistent estimator of . For that, we follow the methodology used for GMM (see, e.g., Harris and Mátyás 1999) and use the results of Komunjer (2012). Let () = (1(), 2()) with dim(1()) = p (and dim(2()) = J – p), and let also ) (i.e., the same split as with ()). Consider the following conditions:

If Ω is estimated by , then , that is, the weighting matrix converges to a population weighting matrix.

is a consistent estimator of the WV ν2 = [ν2(τ)], …, J.

Let ∈ Θ with Θ being an open subset of ℝ. Then, there exists a bijective, twice continuous differentiable function : ℝ → Θ, such that (i.e., the determinant of the Jacobian matrix is nonnegative) and if p = 2, (∂/∂γ)k(γ) > 0, ∀ j (i.e., the Jacobian matrix is positive-definite definite).

1() is twice continuously differentiable.

|(∂/∂)1()| is nonnegative.

||1()|| →| ∞ whenever |||| → ∞.

For every c ∈ ℝ, the equation 1() = c has count-ably many (possibly zero) solutions in ℝ.

For p = 2, (∂/∂)1() is positive definite, and for p > 2, the set of points for which rank[(∂/∂)1()] < p − 1 is bounded.

For Ω satisfying (C.1), under (C.2), for () or –() satisfying (C.4)–(C8), and if there exists a function satisfying (C.3), then defined in Equation (10) is consistent.

The proof is given in the supplementary materials (see Section D). Condition (C.3) defines a reparameterization from ∈ Θ to ∈ ℝ. Hence, in practice, for a given model F one has first to check that a function = [k] exists that satisfies (C.3). For variance parameters 0 < σ2, the k's are the logarithm, and for other restricted parameters with restrictions of the form a < < b, one can use the logit function on (θ – a)/(b – a), that is, k(γ) = (b – a) exp(γ)((1 + exp(γ)) + a.

Condition (C.3) is convenient but actually too restrictive. The following corollary actually extends the conditions of Theorem 1.

If the bijective, twice continuous differentiable function does not satisfy the conditions on the Jacobian matrix in (C.3), then defined in (10) is consistent under conditions (C.1), (C.2), and (C.4)–(C8) in which (∂/∂)1() is replaced by (∂/∂)1()(∂/∂)().

The proof follows directly from the proof of Theorem 1.

In Theorem 2 below, we provide the consistency conditions of the GMWM when F is made of the sum of L independent processes. Let their respective WV be ϕ((() and partition as = [(], we need the following set of conditions:

The functions ϕ((() are linearly independent with all linear combinations of the other ϕ( ((), l ≠ l′, l, l′ = 1, …, L, and for all τ, that is, one cannot find a,b ∈ ℝ, m = 1, …, L such that .

For all j, ϕ( (() = ϕ( ((∗) if and only if ( = (∗ for all l.

Each element of , where ( = [(], m = 1, …, dim(() depends on τ for at least L − 1 of the l = 1, …, L.

Let F be the stochastic process associated to {Y ∈ ℝ} made of the sum of L independent processes with respective WV ϕ( ((), be the GMWM estimator of , then for Ω satisfying (C.1), satisfying (C.2), and ϕ( (() satisfying conditions (C.9)–(C.11), is a consistent estimator of .

The proof is given in the supplementary materials (see Section E). Condition (C10) can be verified for each subprocess WV using the conditions of Theorem 1. We do that in the supplementary materials (see Section F) for some of the processes presented in Table 3 (see Section A), and for the composite processes treated in this article, we show the consistency using the conditions in Theorem 2.

Obviously, the number of scales J should be J ≥ p but at the same time, as discussed in Section 4, a too large J can introduce variability in the estimator. The study of a possible optimal value for J is out of the scope of the present article. However, we may notice that in Guerrier et al. (2013), the authors investigated through simulation studies different choices for the number of scales J for mixtures of WN, QN, drift, RW, and AR(1) models, and found out that from a given number of scales, the variability of the GMWM does not really improve. This suggests that one could choose a number J∗ such that J∗ = min(J, a · p), with a a bounded (positive) integer.

Then, under the conditions of Theorem 1 (or possibly Corollary 1) or Theorem 2 as well as Theorem 3, we have that is asymptotically normal, with asymptotic covariance matrix given by

where

and where D = ∂()/∂, and is the asymptotic covariance matrix of with elements given in (A-3) (see Section C). When Ω = I, then . The most efficient estimator is obtained by choosing , leading then to . In practice, the matrix D is computed at .

The estimation of the WV covariance σ2 of given in (A-3) is in general not straightforward. In the supplementary materials (see Section G), we show that under the assumption of a Gaussian process for Y, a suitable estimator is given by

where M = min(M, M) (see Equation (6)). Alternatively, when the process is not Gaussian or when the sample size is very large as is the case with the dataset analyzed in Section 4 so that the computation of Equation (13) is infeasible, one can use a parametric bootstrap to estimate cov(. Q samples of size T are simulated from on which Q WV and , q = 1, …, Q, are computed and σ2 is estimated by their empirical covariance.

When J > p, that is, the number of WV is greater than the dimension of the parameter vector , one can assess the goodness of fit of the model F to the data by testing the hypotheses using the χ2-test statistic

which is asymptotically χ2 under H0 (see Hansen 1982). It will be used to assess the fit of the postulated model in the case study in Section 5. The investigation of the finite sample properties of Equation (14) is left for future research.

As a possible extension of the GMWM when analytical expressions for () in Equation (10) are too complicated to compute, one can resort to simulations to compute () and hence place the GMWM in the framework of indirect inference (Gourieroux, Monfort, and Renault 1993; Smith 1993; Gallant and Tauchen 1996). Basically, given a sample of observations {y, t = 1, …, T} and a hypothetical model F, is defined as the WV estimated from the sample using Equation (6). Let also with being the WV estimates computed on simulated series {y∗( (), t = 1, …, T}, r = 1, …, R. Then, and are used in Equation (10) to obtain an estimate of , whose properties are described in, for example, Gourieroux, Monfort, and Renault (1993). In particular, for R sufficiently large, . In that case, B can be computed numerically.

A step-by-step algorithm for computing the GMWM estimator (10) given a dataset of size T and a choice for F, is given by the following steps:

Compute the empirical WV for a suitable number of scales J ≥ dim() using Equation (6).

Construct the matrix Ω. We recommend using for Ω the inverse of a diagonal matrix whose elements are ŜW (0) given in Equation (8) for the considered scales. This choice ensures that condition (C.1) is satisfied.

Compute the model-based WV ().

Minimize the quadratic form defined in Equation (10) in to obtain the GMWM estimator .

Compute confidence intervals for by parametric bootstrap or using Equation (11) with replaced by .

Hence, the method is rather easy to implement in practice. Assuming that we have a time series, say y of length T and a model F made of the sum of processes satisfying Corollary 3, the estimation of starts with the computation (only once) of the empirical WV using Equation (6) and their estimated variances using Equation (8) (that are used in Ω). The model-based WV is then analytically built using the sum of the model-based WV (in Table 3) of the independent subprocesses in F. Given all these ingredients, the quadratic form defined in Equation (10) is then optimized. This optimization problem is very simple. Confidence intervals for are, however, generally computed by parametric bootstrap. This task is quite easy to implement, but of course implies simulating data and repeating several times Steps 1–4. This can take in practice about a minute with massive datasets such as the one considered in the case study section.

If a simulation-based approach (indirect inference) is preferred, then Step 3 is not necessary and Step 4 is replaced by the following steps:

Choose a sufficiently large value for R (R = 100 is usually a sensible choice)

Minimize the quadratic form defined in Equation (10) in which () is replaced by ∗(), each element being computed using Equation (6) in which the observations are replaced by a simulated sample from F (for a given current value of ) of size RT.

can be then used for inference, for sufficiently large R (see, e.g., Genton and Ronchetti 2003).

SIMULATION STUDY

This section is dedicated to the finite sample performance evaluation of the GMWM estimator compared to the MLE and the LSE. For the GMWM estimator, one has to make a choice about Ω in Equation (10). The choice is not straightforward in practice, as it is the case for the moment-based estimator GMM or Indirect-Inference-based estimators. Asymptotically, it is well known that the optimal choice for Ω is any matrix proportional to . This does not mean that this choice is better in finite samples, especially because the asymptotic covariance matrix needs to be estimated. A simple choice is the identity matrix, and in our case an intermediate solution is the inverse of the diagonal matrix made of the variance estimates of the MODWT using Equation (8). Indeed, in our simulations, we found that choosing may lead in some cases to numerical instability, probably due to the estimation of the covariances using Equation (13). Anyway, asymptotically these three choices for Ω lead to a consistent GMWM (provided the other conditions are also satisfied). When optimizing Equation (10), we use a quasi-Newton optimization method. For the MLE, we use the EM algorithm together with the Kalman smoother (EM-KF) as proposed by Shumway and Stoffer (1982) (see also Holmes 2010).

An important engineering application is sensor calibration in which the behavior of the errors affecting signals has to be well understood and modeled. In particular, gyroscopes and accelerometers which are part of inertial navigation systems used for space, aeronautical, ground, and underwater applications are subject to random errors that largely affect the quality of the positioning and navigation solution over time. A widely accepted model for the error behavior of an inertial sensor is given by the following state-space model

where Δt is the time interval between two consecutive measurements, u = ωΔt, w (0, q) with q = σ2GM(1 − e−2), and υ (0, σ2WN). The observed series is hence made up of the sum of a first-order GM process with correlation time β−1 and variance σ2GM, a drift with slope ω, and a WN with variance σ2WN. Even with such a relatively simple model, the estimation task is nontrivial.

For our simulation study, we consider simulated processes {y, t = 1, … T] from model (15) as well as an AR(1) with autocorrelation parameter near the boundary value of 1. From model (15), we actually generate three types of samples which correspond to three different (sub)models. In Model 1, we set σ2WN = ω = 0, hence we have only a first-order GM, in Model 2, we set ω = 0, hence we remove the drift only from the complete model (15) (Model 3). We generate 100 processes of size T = 6000 with Δt = 1 and = (σ2WN, σ2GM, β, ω) = (4, 16, 0.05, 0.005) at the complete model. For the submodels, the parameters are constrained accordingly and they are not estimated.

For both the GMWM (simulation based or not) and the EM-KF, the initial values for the optimizations were set to (start) = (1.0, 1.0, 1.0, 0.0), which is relatively far away from the true simulation values. We found that the choice for the starting values is not a serious issue for the computation of the GMWM. The root mean-squared errors (RMSE) as well as the relative RMSE (relative to the true parameter value) are presented in Table 1 for Models 1–3, and for the GMWM and the EM-KF. We also tried a simulation-based version of the GMWM with R = 100 and found out that the RMSE are of the same order as for the GMWM with analytical WV (results not presented here). The sample WV were computed for J = 12(< log(6000)/log(2) ≈ 12.55) scales for all models. The results show that for the smaller models 1 and 2, the RMSE is smaller for the EM-KF than for the GMWM estimator, while the RMSE of the EM-KF explodes for the complete model 3. This last feature is actually well known with models with a drift component. For example, Hinrichsen and Holmes (2009) consider a multivariate state-space model composed of a drift with unknown rate and a random walk process to model the growth of an ecological population with observations that are subject to measurement error. They actually use a two-step estimation procedure by which the drift is first estimated by linear regression and then removed from the state equation leading to a simpler model that is well estimated by means of the EM-KF (see also Stebler et al. 2011).

Table 1.

RMSE and relative RMSE (R-RMSE) of the GMWM and EM-KF estimators for 100 simulated processes of size T = 6000 from model (15) (Model 3) with Δt = 1 and = (σ2WN, σ2GM, β, ω) = (4, 16, 0.05, 0.005) and from submodels with σ2WN = ω = 0 (Model 1) and with ω = 0 (Model 2)

		GMWM		EM-KF
		RMSE	R-RMSE	RMSE	R-RMSE
Model 1	σ2GM	2.00	0.13	1.26	0.08
	β	8.14 · 10−3	0.16	4.09 · 10−3	0.08
Model 2	σ2GM	1.92	0.12	1.55	0.10
	β	1.05 · 10−2	0.21	4.94 · 10−3	0.10
	σ2GM	0.13	0.03	0.14	0.04
Model 3	σ2GM	0.96	0.06	74.57	4.66
	β	4.63 · 10−3	0.09	0.04	0.85
	σ2GM	0.11	0.03	0.16	0.04
	ω	2.79 · 10−4	0.06	0.12	23.58

When the EM-KF behaves well (models 1 and 2), it has a better performance in terms of RMSE than the GMWM estimator. However, one can further improve the efficiency of the latter by decreasing the number of scales J at which the WV are estimated. Indeed, in this example, J = 12 scales are used to estimate two or three parameters, and if more scales are added (supposing a larger T), this only introduces more variability in the GMWM estimator. An optimal choice (in terms of GMWM estimator's efficiency) of the scales and their number is beyond the scope of the present article. The number of scales is obviously a function of the number of parameters p, but their choice (among the possible ones) depends on the model that is considered. For example, with model (15) without the drift element (i.e., Model 2), according to the WV graph in Figure 4, if the last three to four scales are ignored, then the WV are still able to capture information about the other model components. Actually, removing the last four scales improves the efficiency of the GMWM estimator in the simulation study for Model 2 (results not shown here).

An important feature of the GMWM estimator is that it can be used for models for which the MLE (or other classical estimators) is seriously biased. For example, maximum likelihood estimation of ARMA processes, and/or least squares estimation (based on the Yule–Walker equations) of AR processes, when the roots lie near the unit circle, provide estimators that can be seriously biased (see, e.g., Andrews 1993). With the GMWM estimator with Haar coefficients for the WV, we do not have this boundary problem, since, for example, with an AR(1) process, with autocorrelation parameter ρ, when the latter tends to the value of 1, we get a random walk which is a nonstationary process. However, the wavelet coefficients W actually are backward differences of order at least 1, so that the problem of nonstationarity disappears. Note also that taking first differences of a stationary AR(1) process is known to lead to the problem of overdifferencing, since the resulting series is not invertible (it becomes an ARMA(1,1)) which has more parameters than the original one. Consequently, its parameters will be difficult to estimate, and it will tend to degrade the quality of forecasts, so differencing is not always a suitable method to overcome the estimation problem.

As an illustration, we simulated 1000 AR(1) processes of size T = 5000 with ρ ranging from 0.8 to 1 (and with residual variance σ2 = 4) and estimated the LSE based on the Yule-Walker equations together with the GMWM (J = 12). Figure 1 presents the relative RMSE of the LSE and GMWM estimator for both parameters. The GMWM estimator of ρ is less precise than the LSE, except when ρ is close to 1. For σ2, the RMSE of both estimators are quite similar for ρ smaller than 0.9 but for higher values the bias and variance of the LSE explodes showing that there is a boundary effect that affects the LSE. The MLE provides RMSE that are sensibly smaller than the LSE (not presented here) but for values of ρ larger than 0.95 the MLE is difficult if not impossible to compute.

Figure 1.

Relative RMSE based on 1000 simulations from AR(1) processes with ρ ranging from 0.8 to 1 and with residual variance σ2 = 4, of the LSE relative to the GMWM estimator of ρ (“□”) and of σ2 (“o”).

CASE STUDY: INERTIAL SENSORS

An inertial navigation system exploits observation from usually an orthogonal triad of gyroscopes and accelerometers measuring angular rates and specific forces, respectively. Once initialized, the navigation is a dead-reckoning process (i.e., the solution at one epoch is computed from current observations as well as from the solution of the previous epoch) in which gyroscope signals (after suppressing Earth rotation) are integrated once to yield orientation, and accelerometers (after accounting for gravity and apparent forces) are integrated twice to get the velocity, and finally the position in three-dimensional space. The advantage of inertial navigation is in the autonomy (i.e., no external infrastructure is needed) and the large sensor bandwidth (unlike satellite positioning system), while its weakness lies in the time-dependent error behavior (due to the integration process). Indeed, the sensor signals are corrupted by random errors, making the resulting positioning and attitude (or orientation) error increase rapidly with time. To bound the error growth, inertial sensors are combined with other sensors (e.g., GPS, odometer, and altimeter) through optimal (e.g., Kalman) filtering. Such combined navigation systems are not only the base for manned and unmanned vehicle (e.g., spacecrafts, air-crafts, cars, and robots) trajectory and orientation determination, but also for pedestrian and indoor localization devices. The gyroscope and accelerometer error behavior is modeled through state augmentation in the system where the augmented states account for sensor errors usually modeled as stochastic processes. This requires that the parameters of the stochastic processes in the state-space model are carefully set up prior to filtering. The quality of this setup (or estimation) strongly influences the quality of the predicted trajectory. Traditionally, the Allan variance with ad hoc graphical estimation (see, e.g., IEEE Standard 1293 1998), or the KF-Self-Tuning approach (see Niu et al. 2007; Waegli 2009; Stebler et al. 2012), or even the MLE on state-space models (with simple models) (Stebler et al. 2011) are used to both determine the error's model structure and to estimate their parameters. The GMWM estimator offers an alternative approach that is especially suitable when the stochastic error process is rather complex.

In this section, we consider the angular rate signal issued from a real Micro-Electro-Mechanical System (MEMS) gyroscope where the spectral structure of errors is often complex. Such sensors have a great potential in navigation due to ergonomic (smaller and lighter equipment) and economical considerations (El-Sheimy and Niu 2007). Therefore, MEMS sensors are increasingly used in a very wide range of applications such as 3D input devices, robotic, virtual reality, vehicle stability control and so forth, and research focuses on modeling and compensating for their large and variable measurement errors. To validate the use of the GMWM for modeling the MEMS’ errors, we first analyze the signal produced in static conditions at a frequency of 100 Hz (833, 685 measurements). We provide GMWM estimates on a rather complex model that cannot be estimated by an alternative standard method. To show the impact of the importance of model structure and estimation precision of its parameters in positioning, we then present an emulation study in which the trajectory of a helicopter performing airborne laser scanning is used and for which a nearly perfect (i.e., reference) trajectory can be computed. The emulated measurements along the real trajectory are corrupted with an error signal observed in static condition, and the three methods (AV, the KF-Self-Tuning Approach, and GMWM) are used to model the error signal. The predicted trajectories based on the corrupted measurements and estimated error model, together with GPS measurements, are then compared to determine the positioning error.

After removing the mean in the signal made of 833, 685 measurements produced in static conditions at a frequency of 100 Hz, the time-varying part of the sensor errors is directly available and is presented in the time domain in Figure 3 (upper panel), together with the Haar WV of this process with 95% confidence intervals (lower panel). The WV computed on the original signal give an indication of the underlying stochastic processes that are summed up to build up a composite process. A possible model is a mixture of three different first-order GM processes , with and Δt = 0.01 s, the sampling interval. In such a case, the parameters to be estimated are = {σ21, β1, σ22, β2, σ23, β3} and the GMWM is identifiable since a first order-Markov random process is a one-to-one reparameterization of an AR(1) process whose sum allows identifiable GMWM (see Corollaries 2 and 3). The GMWM estimates of the parameter set and its corresponding 95% confidence intervals using the WV covariances estimated by means of a parametric bootstrap (400 samples) are given in Table 2. The suitability of the estimated model for the data at hand can be judged graphically by a matching of the empirical WV and the parametric WV using the estimates in Table 2, as it is done in Figure 3 (lower panel). It can be seen that the estimated process matches the observed one very well across WV scales, yielding an indication of a suitable fit. Moreover, one can also compute the goodness of fit test statistic in Equation (14) for competing models and test to what extent the chosen one yields a suitable fit. For that, we also considered a single GM process and the sum of two GM processes. The test statistics and corresponding p-values are 384363.9 (on 17 df) and p-value < 10−5 for the single GM process, 84787.01 (on 15 df) and p-value < 10−5 for a sum of two GM processes, and 16.81637 (on 13 df) and p-value = 0.208 for a sum of three GM processes, respectively, clearly indicating that the last model is suitable. The corresponding estimated WV are drawn in Figure 3 (lower panel).

Figure 3.

Gyroscope-observed error process (top panel) and graphical comparison (log-log scale) between the Haar WV (line “o”) computed from the observed signal and the analytical signal using the estimated parameters of, respectively, the sum of three GM processes (line “▴”), the sum of two GM processes (line “•”), and one GM process (line “▪”).

Table 2.

Estimated parameters with associated 95% confidence intervals for the mixture of three first-order GM random processes with the Gyroscope signal data

	Estimates	IC(·,0.95)
β1	2.1720 · 102	(2.1465 · 102; 2.1906 · 102)
σ1	7.4521 · 10−3	(7.4477 · 10−3; 7.4688 · 10−3)
β2	6.0693 · 10−1	(2.7890 · 10−1; 7.8655 · 10−1)
σ2	2.9691 · 10−4	(2.9494 · 10−4; 2.9870 · 10−4)
β3	3.5563 · 10−3	(1.5333 · 10−3; 3.9086 · 10−3)
σ3	5.5127 · 10−4	(5.4038 · 10−4; 5.6223 · 10−4)

It should be noted that the sum of three first-order GM random processes can be reparameterized as an ARMA(3,2) process (see, e.g., Granger and Morris 1976; Terasvirta 1977), so that one could, in principle, estimate the latter instead of the former. However, one should then be able to invert the estimated ARMA since very often, and in particular in the example at hand, a sum of first-order GM random processes is the one that we believe explains the real underlying process, which for gyroscopes reflect the short- and long-term correlated noise due to quantization and mechanical thermal noise phenomena present in MEMS sensors (see Kittel 1958; Drexler 1992). To recover the sum of GM processes' parameters from an estimated ARMA process, and also their standard errors, several conditions need to be satisfied. First, the roots of the processes must lie outside the unit circle. Second, the Jacobian matrix of the transformation between the two parameterizations must be invertible to apply the delta method (Rao 1973; Benichou and Gail 1989). With the inertial sensor data, the estimated processes have roots that are near the unit circle, and moreover, the Jacobian matrix of the transformation evaluated at the estimated parameters is not invertible. In that case at least, estimating an ARMA process and converting the estimated model and inference to the sum of first-order GM processes is infeasible.

To illustrate the impact of the importance of model structure and estimation precision in device positioning, we also performed an emulation study that was presented in Stebler et al. (2012). Using extremely reliable equipments, the trajectory of a helicopter performing airborne laser scanning can be precisely determined. From this trajectory, we computed the theoretical perfect inertial measurements (i.e., accelerations and angular rotations) that should be observed along this trajectory. These perfect measurements were then corrupted with a real static error signal acquired such as the one analyzed above. Then, these “pseudo” measurements were used together with the real GPS measurements to compute a trajectory. At one point, we introduced relatively short artificial gaps in GPS observations of 1 min duration (note that, in practice, such gaps are very common) and compared the deviations of the trajectories obtained with different inertial error models estimated from the data. The two first estimated models are based on the AV and on the KF-Self-Tuning approach. A detailed discussion on these methods can be found for example in Stebler et al. (2012). The GMWM was also used as an alternative estimator and model building approach. The trajectories are depicted in Figure 2 (adapted from Stebler et al. 2012) along with the “true” trajectory. It can be seen that the GMWM-based model (black-dashed line) limits significantly the error growth during the GPS-signal outage compared with the other two benchmark methods which diverged from the “true” trajectory by several thousand of meters! The poor performances of the standard methods explains the recent explosion of the research conducted to determine the stochastic modeling of MEMS-type inertial sensors.

Figure 2.

Comparison between a reference trajectory (black dotted line) issued from a mapping flight in which a GPS outage was introduced, with estimated error model based on the AV (light-gray line), the KF-Self-Tunning (dark gray line), and the GMWM (black dashed line).

CONCLUDING REMARKS

In this article, we present a new estimator for the parameters of composite stochastic processes which is shown to be consistent for the class made of the sum of independent white-noise, drift, quantization noise, random walk, and k < ∞ AR(1). As demonstrated in Stebler et al. (2012), the GMWM has many advantages over existing alternative methods for applications in engineering or natural sciences. The enlargement of the class of models for which the GMWM is consistent involves verifying the conditions provided in Theorem 2 and more generally in Theorem 1 (with Corollary 1) and is left for future research. Also, one could use in principle another consistent estimator for the Haar WV, such as the one proposed in Nason, von Sachs, and Kroisandt (2000).

As an anonymous referee pointed out, one could wonder why to use wavelets instead of making inference directly on the spectral density for which a similar estimation procedure could be used. We believe that our approach is more suitable for the following reasons. First, inference on the PSD would make the optimization of a least-squares-type measure (between the empirical and model-based PSD) more difficult to solve when the PSD has large variability over very narrow frequency bands. As shown in Percival and Walden (2000), the wavelet coefficients at scale τ are associated with frequencies in the interval [1/2, 1/2] and relation (5) can be approximated by

This means that the WV summarizes the information in the PSD using just one value per octave frequency band. Thus, it is particularly useful when the PSD is relatively featureless within each octave band. In the case of the widely used pure power law processes (S(f) ∝ |f|α), for example from Equation (16) one gets ν2(τ) ∝ τ−α−1, meaning that no information is lost in using the PSD summary given by the WV. Second, the computation of empirical WV is more straightforward than nonparametric PSD; for example, the periodogram is an inconsistent estimator of SY(f) and can be badly biased even for large samples sizes (because of frequency leakage effects) so that more sophisticated PSD estimators and/or smoothing techniques such as prewhitening or tapering should be used. Third, the PSD of two important processes in sensor error models, namely the drift and the random walk, cannot be distinguished (both have slope of − 2 in a log–log representation of the PSD). Finally, the MODWT on which the WV computation is based requires a number of multiplications of order T log2 T, which is the same order as the widely used fast Fourier transform algorithm, so that the use of WV does not increase the computational burden.

SUPPLEMENTARY MATERIALS

The proofs and some additional results are given in the supplementary materials. This document is organized as follow. In Section A we provide the analytical version of the Haar WV for the processes considered in this article. In Section B we explain de graphical use of the WV for model building purposes. Theorem 3 which demonstrates the multivariate normality of the WV is given in Section C. The proofs of Theorems 1 and 2 are given in Sections D and E. Corollaries 2 and 3 which provides some additional results regarding the consistency of composite stochastic processes are given in Section F Finally, in Section G we explain the construction of the covariance estimator of the WV presented given (13).