Adjustment models for multivariate geodetic time series with vector-autoregressive errors

1 Introduction

Geodetic observation models of surveyed phenomena often include unknown quantities in terms of parameters to be estimated by adjustment computations. The most general structure of observation model is described by non-linear condition equations in which multiple observations and unknown parameters may be linked to each other [1]. This general case of adjustment calculus was called Gauss-Helmert model (GHM) by Wolf [2]. There exist numerous special cases of this model that have been studied extensively, e. g., the errors-in-variables (EIV) model adjusted by the method of total least squares [3] or by its various generalizations [4], [5], [6], [7]. Another special case of the GHM is the Gauss-Markov model (GMM), in which the observations are still linked to functional model parameters but occur in separate equations [8], [9], [10]. Such models may also contain variance-covariance information about the measurements; this information defines the so-called stochastic model. Since adjustment calculus has been elaborated on the basis of matrix calculus, variance-covariance information regarding the observables is usually arranged as a matrix, called the “variance-covariance matrix” (VCM).

Since the dimension of the VCM equals the number of observations, the storage and processing of this matrix in the course of the adjustment computations can be challenging, especially because modern sensors such as accelerometers, laser scanners, laser tracker or satellite-based sensors tend to provide huge data sets consisting of millions or even billions of observations even within relatively short periods of time. It is usually not an option to omit the variance-covariance information within the adjustment computations since the standard parameter estimators lose desirable quality characteristics when this information is neglected. For instance, the least-squares-estimator in the GMM ceases to be the best linear unbiased estimator (BLUE) when the VCM is dropped in the estimation equation [11]. Furthermore, it is often desirable in geodetic data analysis to obtain a realistic description of the uncertainties or variance-covariance information of the estimated parameters, which information is unavailable when an adequate stochastic model of the observables is missing. The “curse of dimensionality” concerning the VCM is aggravated by the increasingly practiced combination of multiple observation series stemming, e. g., from 3D sensors or from sensor networks. Time series of such measurements, e. g., GPS and laser scanner measurements, frequently have been found to be correlated in such a way that the VCM is fully populated, e. g., [12], [13], [17], [16], [14], [15].

To alleviate the computational burden or even computational infeasibility in connection with the storage and processing of a huge VCM, an alternative approach to the stochastic modeling of variance-covariance information can be applied. The goal of that approach is to fully de-correlate the observables by transforming the observation model in such a way that the VCM of the transformed observables becomes a diagonal matrix, e. g., the identity matrix [18]. A diagonal VCM can easily be handled within the adjustment procedure [19], e. g., through vectorization of the occurring matrix products. In geodetic applications the correlations of the observables may frequently be described as colored noise [20], in which case the stochastic model can be expressed as an autoregressive (AR) or autoregressive moving average (ARMA) process and the transformation of the observation equations achieved by means of a computationally efficient digital de-correlation filter [21], [22]. Such processes have also been estimated by means of total least squares based on an EIV model [23]. To model auto- and cross-correlations of multivariate time series, the aforementioned AR and ARMA processes have been extended to vector-autoregressive (VAR) and vector-autoregressive moving average (VARMA) processes [24]. The model order of these recursive processes can be specified to define how far the correlations reach into the past. Thus, VAR(MA) processes can be used to model quite detailed patterns of auto- and cross-correlations. The innovation of the current contribution is to incorporate VAR processes into the GHM as colored noise models.

The white-noise components of (V)AR or (V)ARMA processes are usually assumed to be normally distributed. However, this assumption may be unrealistic or at least questionable in a practical situation, e. g., due to numerous outliers afflicting the measurements. It may then be safer to replace the normal distribution by a larger family of distributions defined by probability density functions having thicker tails than the “ideal” normal distributions. For instance, the usage of generalized t-distributions and of scaled t-distributions has been proposed in connection with ARMA processes [25]. The scaled, multivariate t-distributions have also been applied in connection with VAR processes within rather simple GMM [26], [27]. A key feature of using multivariate t-distributions is that the associated degree of freedom (df) can be estimated alongside the functional parameters, the VAR coefficients and the scale or cofactor matrix; such a method may be called a self-tuning estimator [28] since it does not involve a tuning constant for classifying the in- and outliers, in contrast to for instance Huber’s M-estimator [29]. The df provides evidence of deviations of the measurements’ actual probability distribution from a normal distribution: For large df the multivariate t-distribution is similar to a multivariate normal distribution, whereas the t-distribution has much thicker tails than the latter for a small df, which may indicate the presence of numerous outliers.

The estimation procedure is implemented as an iteratively reweighted least-squares (IRLS) method. The weights are treated as random variables, which are numerically determined by the expectation step within an expectation maximization (EM) algorithm [30]. The various types of model parameters are estimated group-by-group within the maximization step. To deal with the condition equations within the GHM, constrained EM in the spirit of [31] is applied. Due to the down-weighting of observations with extreme errors,^[1] this method can be expected to yield a robust estimator for the GHM. Koch introduced the GHM with t-distributed errors and subsequently extended the model to include variance components [32], [33]. This type of GHM is now extended to include a VAR model with multivariate-t-distributed errors. This yields an approach complementary to that of modeling a small number of outliers as additive, unknown offsets within the afflicted observations, as done in [34], [35]. Therefore, the proposed adjustment model is directed at applications in which rather large numbers of outliers can be expected. To model the auto- and cross-correlations, we do not employ the more general VARMA processes since we found the estimation of the moving average part by means of an EM algorithm to diverge.

The paper is structured as follows. In Section 2, we specify the multivariate observation model, the correlation model, and the stochastic model, which jointly define the GHM with VAR and multivariate t-distributed errors. In Section 3, we formulate the optimization principle employed to adjust this model, and we derive the normal equations to be solved within the EM algorithm. We show that this algorithm estimates the functional parameters, the VAR coefficients and the parameters defining the multivariate t-distribution essentially via IRLS. The Monte Carlo simulation described in Section 4 serves the purpose of validating the algorithm, by showing the biases to be expected in the practical situation of data approximation by a circle based on a GHM with VAR and multivariate t-distributed errors. In Section 5 we demonstrate the application of the algorithm for the adjustment of a real, multivariate accelerometer data set, whose offset and drift are modeled as part of a GMM (viewed as a special case of the GHM) with VAR and multivariate t-distributed errors. In a second applied scenario real laser tracker measurements with outliers are adjusted to estimate the parameters of a sphere employing the proposed GHM with VAR and multivariate t-distributed errors. Finally, some limitations and consequently potential extensions of the presented methodology to be explored in the future are outlined in the conclusions and outlook.

2 Adjustment models for multivariate time series

The purpose of this section is to develop adjustment models suitable for the fitting of multivariate geodetic time series. According to the standard categorization in geodetic adjustment theory these models have two main components: the deterministic model and the stochastic model. The following sub-sections demonstrate some possible structures of these models that have been found useful in geodetic time series analysis. These structures allow for certain generalizations of the standard GHM and GHM by allowing their random deviations to be auto- and cross-correlated through (V)AR processes while potential deviations from a normal distribution are modeled stochastically by a multivariate t-distribution.

3 Estimation

Maximum likelihood estimation of the model parameters θ GHM based on the log-likelihood function (81) requires values for the unobservable data w. Since a stochastic model has been specified for these quantities, such values can be determined as part of an EM algorithm [30]. A similar approach has previously been carried out for numerous models that can be identified as special cases of the current GHM. In all previous instances the EM algorithm was found to be suitable in view of its stable convergence and its computationally simple form as IRLS (see, e. g., [47], [48]), where the entries of w play the role of the weights. The iterations are indexed by s in the following. In the expectation step (“E-step”) of the EM or IRLS algorithm the weights are computed as conditional expectations of the random variables W based on their stochastic model, using the current solution θ ˆ GHM ( s ) . In the subsequent maximization step (“M-step”) a new solution θ ˆ GHM ( s + 1 ) is computed by maximizing the conditional expectation of the log-likelihood function (“Q-function”) using the current weights.

The Q-function for the GHM with VAR and multivariate t-distributed errors is given by (82) in the Appendix. Evaluation of the conditional expectation of the random weight variables W yields the computational formula

The tilde on the conditional expectation w ˜ t ( s ) is used to distinguish this “theoretical” type of quantity from estimates. To maximize the resulting Q-function (85) under the constraints the Lagrangian is defined as

which involves as additional parameters an ( r × 1 )-vector λ of unknown Lagrange multipliers. This approach leads to a “constrained EM algorithm” [31]. In the frequently encountered case that h ( β , μ ) cannot be expressed “linearly” in the form of

for some numerically fixed matrix X, matrix B and vector m, linearization

is applied as shown for instance in [32]. Here, we define the ( m × 1 )-vector of “incremental parameters”

the ( r × 1 )-vector of “misclosures”

the ( r × 1 )-vector of “pseudo-misclosures”

the ( r × m )-Jacobi matrix

and the ( r × N n )-Jacobi matrix

Whereas Δ β is unknown, the quantities (44)–(47) are computable with the given parameter estimates of iteration step s. For notational brevity, the dependency of the latter quantities on the iteration step will not be indicated in the following; thus, we write m, m ‾, X and B instead of m ( s ) , m ‾ ( s ) , X ( s ) and B ( s ) . Forming now the partial derivatives of the linearized Lagrangian (86) given in the Appendix with respect to the different groups of parameters and setting these equal to zero, one obtains

1. the normal equations for the incremental functional parameters Δ β:

   (48) X T λ = 0 [ m × 1 ] ;

2. the normal equations for the location parameters μ:

   (49) ∑ t = 1 n J ¯ t T w ˜ t ( s ) Ψ − 1 [ ℓ ¯ t − J ¯ t μ ) ] + B T λ = 0 [ N n × 1 ]

   denoting

   (50) J ¯ t = A ( L ) J t ,

   (51) ℓ ¯ t = A ( L ) ℓ t ;

3. the normal equations for the VAR coefficient matrices A 1 , … , A p :

   (52) E W ˜ ( s ) e 1 T ⋮ e n T − E W ˜ ( s ) E T A 1 T ⋮ A p T = 0 N ⋮ 0 N ,

   denoting

   (53) E = e 0 ⋯ e n − 1 ⋮ ⋮ e 1 − p ⋯ e n − p

   with e i = 0 3 × 1 for all i ≤ 0 and letting W ˜ ( s ) be the n-dimensional diagonal matrix with diagonal entries w ˜ 1 ( s ) , … , w ˜ n ( s ) ;

4. the normal equations for the inverse cofactor matrix Ψ − 1 :

   (54) Ψ − 1 n ∑ t = 1 n w ˜ t ( s ) [ A ( L ) ( ℓ t − μ t ) ] [ A ( L ) ( ℓ t − μ t ) ] T = 0 N ;

5. the normal equation for the df ν (cf. Sect. 5.8.2 [47]):

   (55) − dg ν 2 + log ν 2 + 1 + 1 n ∑ t = 1 n log w ˜ t ( s ) − w ˜ t ( s ) + dg ν ( s ) + N 2 − log ν ( k ) + N 2 = 0 ,

   involving the digamma function dg;

6. the normal equations for the Lagrange multipliers λ:

   (56) X Δ β + B ( μ − ℓ ) + m ‾ = 0 [ r × 1 ] .

(supposing that the inverse exists), these parameters can be eliminated from (56) via substitution, which results in

Together with (48), this gives the normal equation system

The parameter groups ( λ , Δ β ), ( A 1 , … , A p ), Ψ − 1 and ν are estimated in that order by successively solving (59), (52), (54) and (55). The parameter group μ can be determined through (57) after solving for ( λ , Δ β ). While equation systems (59), (52) and (54) are linear, the solution of (55) requires a zero search. This step-wise solution approach defines an ECM algorithm, which extends previously established ECM algorithms in simpler models.

Having completed s iterations, the first ECM-step of the new iteration s + 1 yields a new solution

from (59), where the unknown matrices Ψ and

are fixed by substituting (i. e., “conditioning on”) the previous solutions Ψ ( s ) and A 1 ( s ) , … , A p ( s ) . Substituting then in addition λ ( s + 1 ) for the variables λ in (57) gives the new solution

Moreover, the new solution

via substitution into the rearranged equation (43).

The solution ( β ( s + 1 ) , μ ( s + 1 ) ) might not fulfill the constraints h ( β , μ ) = 0 exactly due to linearization errors. Therefore, we apply the principle of constrained EM in the spirit of [31] and iterate the estimation of these parameters until either a maximum number of iteration steps is reached or the maximum absolute misclosure falls below a specified threshold (say, 10 − 14 ). To do this, the quantities m, m ‾, X and B in (44)–(47) are recomputed using β [ 0 ] : = β ( s + 1 ) and μ [ 0 ] : = μ ( s + 1 ) , and new solutions λ [ 1 ] , β [ 1 ] , μ [ 1 ] are computed via evaluation of (60)–(62) by substituting again the currently available estimates Ψ ( s ) and A 1 ( s ) , … , A p ( s ) . When one of the stop criteria is fulfilled, say after t additional iterations, one redefines λ ( s + 1 ) : = λ [ t ] , β ( s + 1 ) : = β [ t ] and μ ( s + 1 ) : = μ [ t ] .

Next, the new solution

of the VAR coefficient matrices is computed from (52), based on the new correlated residuals

computed from the model (7) and the corresponding residual matrix

with the usual boundary values e i ( s + 1 ) = 0 3 × 1 for all i ≤ 0. The current VAR model defines a new de-correlation filter A ( s + 1 ) ( L ), by which the correlated residuals e ( s + 1 ) are transformed into the white noise residuals

in view of (33) and (64). Conditional on these values, equation system (54) gives rise to the solution for the cofactor matrix

Finally, ν ( s + 1 ) may be found by a zero search to satisfy (55). It is generally known that the ECM algorithm can be sped up by solving instead the equation

where

Equation (68) is obtained by defining the log-likelihood function directly on the joint pdf (30) of the white-noise series based on the multivariate t-distribution, and setting its partial derivative with respect to ν equal to zero. This modification of the ECM algorithm has been called ECM either (ECME) [49], [50]. The zero search can be performed with mathematical standard software such as the INTLAB library ([51]) or using a grid search. The estimation of the df concludes the current M-step. The E- and M-step (in the form of an ECME algorithm) are repeated until certain stop criteria are fulfilled, e. g., a combination of a maximum number of iterations and thresholds

for the minimum absolute difference between the current and the previous solution. We defined the threshold d ν by a greater value than the other threshold since the df ν fluctuates quite strongly and is the parameter most challenging to estimate. In the following, the solution resulting from the final iteration step is indicated by using “hats”, that is, β ˆ, μ ˆ, etc. This EM algorithm was tested in the Monte Carlo simulation of Sect. 4.1. In the special case of a GMM with VAR and multivariate t-distributed errors, a similar unconstrained algorithm can be derived [26]. The methodology is applied in the real data study of Sect. 4.2 and Sect. 4.3.

4 Results

4.1 Gauss-Helmert model: Adjustment of simulated data of a 2D circle

In the first part of this section, the performance of the EM algorithm derived in previous Sect. 3 in producing correct estimates within a Monte Carlo (MC) simulation is studied. For this purpose, generated measurements of N = 2 time series are approximated by a 2D circle defined by the nonlinear equation

(70) h t ( β , μ ) = x t − c x 2 + y t − c y 2 − r 2 = 0 .

On the one hand, the vector β of functional parameters consists of the circle center coordinates c x , c y and radius r, whose simulated true values are specified by

(71) β = c x c y r = 1 1 2 .

On the other hand, the vector μ = [ μ 1 T , … , μ n T ] T of location parameters consists of n measured points

μ t = μ x , t μ y , t = x t y t ,

where n was set to 250, 500, 2500 and 5000 in different simulation runs. For greater clarity of notation, we write μ x , t , μ y , t instead of μ 1 , t , μ 2 , t in the sequel. The true point coordinates were defined by

x t = c x + r sin Φ t y t = c y + r cos Φ t

with

Φ t = t · 10 n ( t ∈ { 1 , … , n } ) .

Thus, the measurements are equidistantly distributed along the circle, and one half of the circle was measured twice. The numerical realizations e t = [ e x , t , e y , t ] T of the corresponding random errors are assumed to be determined by a VAR(1) process

(72) e x , t e y , t = 0.5653 − 0.0197 − 0.0431 0.7577 e x , t − 1 e y , t − 1 + u x , t u y , t ,

where random numbers u t = [ u x , t , u y , t ] T for the white-noise components were generated from the multivariate t-distribution

U x , t U y , t ∼ ind . t ( 0 , Ψ , ν ) = t 0 0 , σ 0 2 · Ψ x Ψ x , y Ψ x , y Ψ y , ν

with σ 0 2 = 0.001 2 , Ψ x = 1.0, Ψ y = 4.0 and Ψ x , y = 0.5. The given cofactor matrix reflects the assumption that the entries of U t are pairwise correlated by correlation coefficient ρ = 0.25. Simulations were carried out for the df ν = 4 (referred to as “stochastic model A”) and ν = 7 (“stochastic model B”).

In total 500 MC samples of the white-noise components for stochastic model A and B as well as for the different sample sizes n were randomly generated using MATLAB’s mvtrnd routine. The resulting simulated observations in models A and B were then adjusted by means of the EM algorithm. The maximum number of iterations was set to itermax = 50, and the convergence thresholds to ε = 10 − 8 and ε ν = 0.001. The unknown parameters β and Ψ of the simulated observations of Model A were additionally estimated using the classical least-squares approach to the GHM (“LS”), without iterative re-weighting and without imposing a VAR model (cf. [18]). The comparison between the results of LS and EM algorithm demonstrates the effect of assuming a simplified stochastic model.

4.1.1 Analysis of the simplified stochastic model on the estimation results

The results of this comparison are visualized in Fig. 1 and Fig. 2 by means of box plots. In each box plot, the middle line shows the median. The bottom and the top border of the boxes are the 25^th percentile (Q1) and 75^th percentile (Q3), respectively. The distance Q3-Q1 is the inter-quartile range (IQR). The additional lines above and below the box are the whiskers defined by Q 1 − 1.5 · IQR and Q 3 + 1.5 · IQR, respectively. The red crosses indicate the outliers in the sense that their value is more than 1.5 times the IQR away from Q1 and Q3. The results of Fig. 1 show that the precision of the estimated parameters β ˆ decreases with increasing number of observations. Overall, the mean values of the estimated parameters are similar for LS and EM. However, Fig. 1 shows for numbers of observations less than 5000 that LS produces numerous extreme MC radius estimates outside of the inter-quartile range far away from the true radius, whereas all radius values estimated by EM lie within that range. Thus, EM is significantly more accurate than LS as far as the estimation of the radius is concerned. Figure 2 shows the estimated scale parameters Ψ ˆ for 500 MC runs by means of EM and LS estimation. The number of outliers (red crosses) resulting from the LS estimate is significantly larger than the number of outliers from the EM algorithm. It can be clearly seen that the bias in the EM-estimated scale parameters is much smaller compared to the LS-estimated scale parameters. As the number of observations increases, the bias of the estimated scale parameters decreases and the dispersion becomes smaller. The bias of all estimated parameter can be calculated for each MC run by

(73) Δ Φ ˆ ( i ) = Φ ˆ ( i ) − Φ ( i = 1 , … , 500 ) .

The statistical moments of the calculated bias Δ β and Δ Ψ of the simulations with 5000 observations are presented in Table 1. For the functional parameters β the true value lies in the 95 % confidence interval for all estimates (LS, EM model A and EM model B). In contrast, the true values of Ψ lie outside the estimated 95 %-confidence interval of the estimated cofactor matrix Ψ ˆ when LS is employed. Applying EM, this statement is valid only for the 95 %-confidence interval for Ψ ˆ x under stochastic model A and for Ψ ˆ y under stochastic model B. Here, the biases of the estimator Ψ ˆ with respect to the true simulation parameter values applying EM is less than the corresponding bias when employing LS. The greater deviation concerning the VCM estimated by LS is the result of not taking the correlations and the heavy-tailedness of the noise into account.

4.1.2 Analysis of the degree of freedom on the estimation results

The influence of the df on the estimation results using the EM algorithm is displayed in Fig. 3 and Fig. 4. Figure 3 shows the estimated β ˆ x , estimated scale parameters Ψ ˆ and the differences of the estimated ν to the true values for 500 MC runs by means of EM model A and EM model B estimation. The values of the estimated functional parameters β ˆ under stochastic model A and stochastic model B show no significant differences. It can be clearly seen that the bias in the EM model B estimated scale parameters is much smaller compared to EM model A estimated scale parameters. The number of outliers (red crosses) resulting from the EM model A estimate is significantly larger than the number of outliers from the EM model B algorithm. For both models (A and B) applies that if the number of observations is increased, the bias and the dispersion becomes smaller. The bias of the estimated df ν ˆ under model A is approximately 2.25, which is significantly larger than the bias under model B (approximately 0.44). It can clearly be seen that increasing the number of observations has no significant influence on the bias of the estimated df ν ˆ for both models. The estimated VAR coefficients for both models (EM model A and EM model B) are shown in Fig. 4. The estimated MC solutions for both models show similar spreading behavior around the true value. The number of observations only influences the precision of the estimation but does not lead to a reduction of the bias. The statistical moments of the calculated bias Δ α ˆ and Δ ν ˆ of the simulations with 5000 observations are presented in Table 2. For the VAR parameters only α 1 ; 2.2 lies in the 95 % confidence interval.

Figure 1

MC simulation results of the estimated parameters for stochastic model A. Dash line: true value of simulation.

Figure 2

MC simulation results of the estimated VCM for stochastic model A.

Figure 3

Comparison of the MC simulation results between the stochastic model A and B. The maximum deviation of Δ ν ˆ for the stochastic model B with n = 250 is 113.

Figure 4

Comparison of VAR coefficients between the MC simulation results of stochastic model A and B.

Table 1

Statistics of MC simulation results for the bias of the estimated parameters β ˆ and Ψ ˆ. The results were calculated for the simulation with n = 5000 observations. The boundaries of the 95 % interval estimates are defined by the empirical 2.5 and 97.5 percentiles.

	LS Model A	EM Model A	EM Model B
Δ β ˆ c x
mean	3.7 · 10 − 6	3.0 · 10 − 6	1.2 · 10 − 6
std	1.1 · 10 − 4	9.8 · 10 − 5	9.0 · 10 − 5
minimum	− 2.4 · 10 − 4	− 2.4 · 10 − 4	− 2.2 · 10 − 4
maximum	3.5 · 10 − 4	3.1 · 10 − 4	2.5 · 10 − 4
95 % interval	[ − 2.0 · 10 − 4 , 2.1 · 10 − 4 ]	[ − 1.9 · 10 − 4 , 1.9 · 10 − 4 ]	[ − 1.7 · 10 − 4 , 1.8 · 10 − 4 ]
Δ β ˆ c y
mean	− 6.9 · 10 − 7	6.0 · 10 − 7	7.1 · 10 − 6
std	2.0 · 10 − 4	1.9 · 10 − 4	1.7 · 10 − 4
minimum	− 7.1 · 10 − 4	− 6.5 · 10 − 4	− 4.7 · 10 − 4
maximum	5.5 · 10 − 4	5.0 · 10 − 4	4.2 · 10 − 4
95 % interval	[ − 3.9 · 10 − 4 , 4.1 · 10 − 4 ]	[ − 3.8 · 10 − 4 , 3.9 · 10 − 4 ]	[ − 3.3 · 10 − 4 , 3.4 · 10 − 4 ]
Δ β ˆ r
mean	6.3 · 10 − 6	3.3 · 10 − 6	3.8 · 10 − 7
std	9.9 · 10 − 5	9.0 · 10 − 5	8.3 · 10 − 5
minimum	− 3.1 · 10 − 4	− 2.7 · 10 − 4	− 2.1 · 10 − 4
maximum	2.8 · 10 − 4	2.6 · 10 − 4	2.7 · 10 − 4
95 % interval	[ − 1.9 · 10 − 4 , 1.9 · 10 − 4 ]	[ − 1.9 · 10 − 4 , 1.7 · 10 − 4 ]	[ − 1.5 · 10 − 4 , 1.6 · 10 − 4 ]
Δ Ψ ˆ x
mean	1.2 · 10 − 6	3.5 · 10 − 7	− 6.1 · 10 − 8
std	2.8 · 10 − 7	1.6 · 10 − 7	4.8 · 10 − 8
minimum	1.2 · 10 − 7	− 2.5 · 10 − 7	− 1.9 · 10 − 7
maximum	2.9 · 10 − 6	1.6 · 10 − 6	1.1 · 10 − 7
95 % interval	[ 7.0 · 10 − 7 , 1.8 · 10 − 6 ]	[ 1.2 · 10 − 7 , 6.7 · 10 − 7 ]	[ − 1.5 · 10 − 7 , 3.2 · 10 − 8 ]
Δ Ψ ˆ x , y
mean	2.8 · 10 − 7	7.4 · 10 − 8	− 9.2 · 10 − 8
std	4.7 · 10 − 7	2.1 · 10 − 7	6.5 · 10 − 8
minimum	− 3.1 · 10 − 6	− 1.4 · 10 − 6	− 3.1 · 10 − 7
maximum	4.0 · 10 − 6	1.7 · 10 − 6	1.1 · 10 − 7
95 % interval	[ − 4.6 · 10 − 7 , 1.3 · 10 − 6 ]	[ − 3.0 · 10 − 7 , 4.6 · 10 − 7 ]	[ − 2.2 · 10 − 7 , 3.5 · 10 − 8 ]
Δ Ψ ˆ y
mean	4.8 · 10 − 6	− 2.6 · 10 − 7	− 1.4 · 10 − 6
std	9.4 · 10 − 7	3.6 · 10 − 7	1.3 · 10 − 7
minimum	2.8 · 10 − 6	− 9.7 · 10 − 7	− 1.8 · 10 − 6
maximum	1.2 · 10 − 5	2.6 · 10 − 6	− 7.1 · 10 − 7
95 % interval	[ 3.4 · 10 − 6 , 6.7 · 10 − 6 ]	[ − 7.5 · 10 − 7 , 5.1 · 10 − 7 ]	[ − 1.6 · 10 − 6 , − 1.2 · 10 − 6 ]

Table 2

Statistics of MC simulation results for the bias of the estimated parameters ν ˆ and α ˆ. The results were calculated for the simulation with n = 5000 observations. The boundaries of the 95 % interval estimates are defined by the empirical 2.5 and 97.5 percentiles.

	EM Model A	EM Model B
Δ ν ˆ
mean	2.25	0.44
std	0.21	0.21
minimum	1.70	−0.19
maximum	3.12	1.16
95 % interval	[ 1.87 , 2.70 ]	[ 0.03 , 0.89 ]
Δ α ˆ 1 ; 1 , 1
mean	6.6 · 10 − 2	6.3 · 10 − 2
std	1.9 · 10 − 2	2.1 · 10 − 2
minimum	0.8 · 10 − 2	− 2.6 · 10 − 2
maximum	11.7 · 10 − 2	12.4 · 10 − 2
95 % interval	[ 2.9 · 10 − 2 , 10.1 · 10 − 2 ]	[ 2.2 · 10 − 2 , 10.3 · 10 − 2 ]
Δ α ˆ 1 ; 1 , 2
mean	− 4.4 · 10 − 2	− 4.8 · 10 − 2
std	2.4 · 10 − 2	2.4 · 10 − 2
minimum	− 14.2 · 10 − 2	− 11.2 · 10 − 2
maximum	3.9 · 10 − 2	3.0 · 10 − 2
95 % interval	[ − 8.7 · 10 − 2 , − 0.2 · 10 − 2 ]	[ − 9.5 · 10 − 2 , − 0.3 · 10 − 2 ]
Δ α ˆ 1 ; 2 , 1
mean	3.3 · 10 − 2	3.4 · 10 − 2
std	1.0 · 10 − 2	1.0 · 10 − 2
minimum	4.8 · 10 − 3	6.5 · 10 − 3
maximum	6.9 · 10 − 2	6.8 · 10 − 2
95 % interval	[ 1.2 · 10 − 2 , 5.3 · 10 − 2 ]	[ 1.6 · 10 − 2 , 5.2 · 10 − 2 ]
Δ α ˆ 1 ; 2 , 2
mean	− 7.3 · 10 − 4	6.3 · 10 − 4
std	1.3 · 10 − 2	1.4 · 10 − 2
minimum	− 3.4 · 10 − 2	− 3.9 · 10 − 2
maximum	4.6 · 10 − 2	3.5 · 10 − 2
95 % interval	[ − 2.9 · 10 − 2 , 2.4 · 10 − 2 ]	[ − 2.7 · 10 − 2 , 2.8 · 10 − 2 ]