Medium- and Long-Term Prediction of Polar Motion Using Weighted Least Squares Extrapolation and Vector Autoregressive Modeling

The x_p and y_p components contain a few well-known periodicities such as CW and AW. A priori model hence consists of a purely linear part (offset and drift parameters a, b, c, d) and, two periodical elliptical motions representing CW (parameters f₁, A_1,x, B_1,x, A_1,y, B_1,y) and AW (parameters f₂, A_2,x, B_2,x, A_2,y, B_2,y), as written in the following. 1 xp(t)=a+bt+∑i=12Ai,xcos(2πfit)+Bi,xsin(2πfit) 2 yp(t)=c+dt+∑i=12Ai,ycos(2πfit)+Bi,ysin(2πfit)

The semimajor and semiminor axes of an elliptical motion correspond to the amplitude of circular motion. For completeness, bias and drift of the linear part and parameters of the CW and AW in this a priori model are estimated simultaneously from x_p and y_p observations by a WLS fit as 3 β=(CTWC)−1CTWL where β represents the vector of estimated parameter (Eq. (4)), L represents the matrix of x_p and y_p obervations (Eq. (5)), in which n is the length of observed time-series, C represents the coefficient matrix (Eq. (6)), and W represents a diagonal matrix of the weights, W = diag[w₁, w₂, ⋯ , w_n], in which the weights w_i are greater than zero, as will be described in the next paragraph. 4 β=[abA1,xB1,xA2,xB2,xcdA1,yB1,yA2,yB2,y]T 5 L=[xp(t1)xp(t2)⋯xp(tn)yp(t1)yp(t2)⋯xp(tn)]T 6 C=[1t1cos(2πf1t1)sin(2πf1t1)cos(2πf2t1)sin(2πf2t1)1t2cos(2πf1t2)sin(2πf1t2)cos(2πf2t2)sin(2πf2t2)⋮1tncos(2πf1tn)sin(2πf1tn)cos(2πf2tn)sin(2πf2tn)]

The deterministic models as shown in Eq. (1) and (2) aim to extrapolate the present trend, CW and AW within x_p and y_p time-series. This can be carried out by the WLS extrapolation. This deterministic models are subsequently utilized for two purposes: (1) to extrapolate the deterministic components in time-series and (2) to get stochastic residuals for WLS misfit (the difference between the deterministic models and data themselves).

A key issue in WLS fit is the appropriate choice of a weighting matrix W, the determination of a proper weight w_i assigned to the observation at time t_i. Schuh et al. (2001) calculated the weights by estimating empirical variance components. Sun and Xu (2012) developed a power weighting function with the optimal power. Wu et al. (2018) determined the weights with the a priori variance of observables. When observations are closer to the epochs of predictions, the greater the effect on predictions are. Considering the effect of intervals between the observation points and prediction epochs on results of LS extrapolation, the weights w_i are selected according to the following two weighting functions in current work. (1)

Equal weighting function: w_i = 1, i = 1, 2, ⋯, n. This simple weighting function is widely used for LS analysis, but does not take into consideration the impact of the geometry.

The VAR model developed by Love and Zicchino (2006) allows one to account for unobserved individual heterogeneity for the entire time-series via the introduction of fixed effects that enhance the consistency and coherence of measurements. Moreover, this model has some practical advantages that make it a very suitable method to model x_p and y_p together. First, unlike multivariate autoregressive (MAR) process, the VAR model makes no distinction between exogenous and endogenous variables, in accordance with the interdependence reality. Instead, all variables are treated as endogenous mutually. Each variable that entered a VAR model not only depends on its historical realization but also correlates with other variables, conforming with real simultaneities and a priori correlation between x_p and y_p. Second, VAR is very uncomplicated for efficient and coherent estimations for a small number of variables, especially for two variables in case of x_p and y_p. Let us denote the bivariate x_p and y_p residual data as a vector R(t) = (r_x(t), r_y(t))′, where r_x(t) and r_y(t) are x_p and y_p residual time-series, respectively. Modeling the residual bivariate time-series R(t) may be performed utilizing the VAR technique.

Given a vector of residuals R(t) = (r_x(t), r_y(t))′ at time t, the basic VAR model of a p-order without deterministic components has the form 8 R(t)=α1R(t−1)+α2R(t−2)+⋯+αpR(t−p)+U(t) where α_i is the 2×2 matrix of VAR coefficients at lags i (i = 1, 2, ⋯ , p), and U(t) = (u_x(t), u_y(t))′ is an unobservable zero mean white noise vector with time-invariant positive definite covariance matrix, in which u_x(t) and u_y(t) are independent or serially uncorrelated. The model is referenced to briefly as a VAR(p) process as the number of lags is p.

There are a lot of model selection criteria to select an order p. They proceed by fitting VAR(m) models with different orders m = 1, 2, ⋯ , p_max and then select an estimator of the p-order that minimizes some criterion. In this work, we apply the popular Akaike information criterion (AIC) to choose an appropriate order. This criterion for order selection is based on the following statistics (Akaike, 1971). 9 AIC(p)=In det(∑u˜(m))+2mk2n where det(·) represents the determinant, ∑_u(m) denotes the estimator of the residual covariance matrix for a VAR(m) model, n is the sample data size, and k denotes the dimension of the vector for modeling, specially k = 2 in the work. Following Johansen (1995), this work estimates values of the coefficient matrix α_i resulting from fitting the VAR(p) model to the two-dimensional residuals using maximum likelihood.

The fitted VAR model is applied in the process of predicting x_p and y_p together. The one-step prediction is computed as 10 R^(t+1)=α^1R(t)+α^2R(t−1)+⋯+α^pR(t−p+1) where R^(t+1) is the predicted vector at time t + 1, α^i , i = 1, 2, ⋯ ,p, are the estimated coefficient matrices. Multistep predictions of residual vector for t + 2, t + 3, t + T are calculated recursively using Eq. (10) by taking place at time t by t + 1.

The IERS regularly releases daily values of the EOP C04 series with a 30-day delay after the combination of solutions of space geodetic technologies including Very Long Baseline Interferometry (VLBI), Satellite Laser Ranging (SLR), Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS), and Global Navigation Satellite System (GNSS). The combined EOP C04 series is routinely published at one-day intervals at the IRES official website (https://www.iers.org/IERS/EN/DataProducts/data.html). The series consists of observations of the x_p and y_p pole coordinates from 1962 until 30 days from now with a high accuracy. The mean and standard deviation for pole coordinates of the differences between combined solutions derived by the RS/PC and C04 over 2018 are less than 0.06 ms (Dick and Thaller, 2020), and therefore have a negligible impact on PM predictions. The IERS directing board adopted the EOP 14 C04 as the IERS reference series in February 1, 2017 (Bizouard et al., 2019). In this work, the EOP 14 C04 series is taken as data source of the x_p and y_p pole coordinates.

The predicted values of x_p and y_p are calculated as the sum of the WLS extrapolations and univariate or vector predictions of x_p and y_p. The mean absolute error (MAE), used as the official 1st EOP PCC statistics, is selected for evaluating the prediction accuracy. The MAE is expressed for the ith day into the future as follows (Kalarus et al., 2010). 11 MAEi=1N∑j=1N|Oi,j−Pi,j| where N is the number of predictions, O_i,j and P_i,j are the observed pole coordinates and their ith point of jth prediction, respectively.

Predictions of x_p and y_p are computed by means of two prediction techniques: LS+AR and WLS+VAR. The range of all predictions is equal to 365 days. Totally, 500 predictions are generated every 7 days at different starting prediction days from January 1, 2012 to July 25, 2021 to calculate the representative statistics, the number N of prediction is 500. The parameters of the linear part, CW and AW in LS and WLS fit are solved by a sliding window analysis for every 365-day prediction. The window size for LS and WLS fit is set to 12 years, approximately twice the beating period of the CW and AW. The above two weighting functions have been applied in WLS fit to extrapolate the linear trend, CW and AW and obtain the residuals of misfit. The comparison of the predictions is carried out by means of the analysis of: (1) the absolute differences between the x_p and y_p data and their predictions (Figure 1 and Table 1), and (2) MAE statistics (Figure 2 and Table 2).

Figure 1.

Figure 2.

The absolute values of the difference between the observed x_p and y_p pole coordinates and their one-year predictions starting at different days do not exceed 85 mas (Figure 1). The larger absolute differences are for the chosen predictions up to prediction horizons greater than 8 months. Those inaccurate predictions with the maximum absolute difference as given in Table 1 are generated for starting prediction epochs in the early days of 2014 (Figure 1), which might correspond to the El Niño event 2014/2015. An El Niño phenomena could cause extreme stochastic perturbations in the seasonal harmonic oscillations and would hence decrease the PM prediction accuracy in time to the El Niño event. In addition, the predictions beginning in the middle of 2016 are inaccurate. During these prediction starting days, the maximum absolute differences are 78.07 mas, 77.84, mas and 66.13 mas respectively calculated by LS+AR, WLS+VAR with equal weights and WLS+VAR with peicewise weights. The interpretation can be correlated with the La Niña 2016/2017. There are also less accurate predictions beginning in the last half of 2020. This can be linked to the La Niña 2020. For the span of tested data, the ”uncorrelated” LS+AR solution generates the predictions of x_p and y_p with the maximum value of difference between the data themselves and predictions of more than 84 mas, while the absolute difference for the ”correlated” WLS+VAR techniques is less than 77 mas. These differences and extreme perturbations are substantially smaller than those calculated for the conventional solution, although the worst predictions obtained by the ”correlated” techniques also correspond to EI Niño or La Niña events. The cause of the poor accuracy of PM predictions can be the period or phase variation of the annual oscillation in the occurrence of such events. (Kosek et al., 2001). The prediction accuracy may be enhanced in such a way that the period or phase of the annual oscillation is treated as a variable rather than a constant and modeled from short-term PM time-series in the LS adjustment (Guo et al., 2013; Zhang et al., 2012). In the case of WLS+VAR, one can assign larger weights to recent observations to fit the variable annual oscillation.

Table 1 gives the maximum absolute differences for all considered predictions of x_p and y_p for 10-day, 30-day, 150-day, 270-day, and 365-day into the future. It can be seen that the maximum errors of the medium- and long-term predictions can be noticeably decreased by the proposed WLS+VAR approach, indicting that the error increasing of both x_p and y_p predictions will be alleviated if two components are modeled and predicted together.

The essential results regarding the accuracy of different techniques come from an analysis of the MAE for the x_p and y_p predictions (Figure 2 and Table 2). The results show that the accuracy of three techniques is comparable for the predictions up to 30 days into the future, while the accuracy of the predictions out 30 days computed by WLS+VAR is better than that by the conventional LS+AR approach. In particular, the enhancement in accuracies is mostly seen in the case of the medium- and long-term predictions for horizons of 60 to 365 days based on WLS+VAR. One key reason is that VAR works better than AR in modeling and prediction of the residuals for WLS misfit. Furthermore, it can be seen that the choice of appropriate weights is crucial. The equal weights w_i = 1 ignore the effect of intervals between the observation points and prediction epochs on the predicted values. The maximum values of absolute difference between the observations and predictions and MAE values obtained by the equal weights are 76.66 mas and 23.40 mas for horizons of 1 to 365 days respectively, larger than the maximum absolute difference of less than 71 mas and MAE of no more than 22 mas by the piecewise weights as shown in Tables 1 and 2. This is owing to the assignment of different weights to observations that can improve the LS fit for the time-variant CW and AW. That is, if appropriate weights are assigned to PM measurement, WLS fit will work better. The figure and statistics of the MAE not only demonstrates that PM predictions would be biased if the correlation between x_p and y_p was ignored, but also confirms that WLS+VAR with piecewise weights is a more accurate method for PM prediction than two others. The accuracies of the x_p predictions are 13.97 mas, 18.47 mas, and 20.52 mas respectively for 150, 270, and 365 days in terms of the MAE statistics, which are 36%, 24.8%, and 33.5% higher than LS+AR respectively and 8.2%, 9.2% and 9.6% higher than WLS+VAR with equal weights respectively for 150, 270, and 365 days. For the y_p predictions, the 150-, 270- and 365-day accuracies are 15.41 mas, 21.17 mas, and 21.82 mas respectively, which are 27.4%, 11.9%, and 21.8% higher than LS+AR respectively and 3.1%, 7.5%, and 6.8% higher than WLS+VAR with equal weights respectively for the corresponding prediction days.

In October 2005, the 1st EOP PCC was launched for making an examination of the property of existing prediction algorithms. There were 13 participants who contributed to the campaign and submitted EOP predictions. At the end of February 2008, this campaign collected about 6,500 submissions of EOP predictions calculated by 20 prediction methods and techniques denoted by ID number. The 1st EOP PCC is a variety of contest where the validation strategy and prediction period were specified clearly in advance. The campaign is refereed to Kalarus et al. (2010) for detail. This work uses WLS+VAR with the piecewise weights to generate the x_p and y_p predictions up to 365 days into the future between January 1, 2005 and February 28, 2008, a same prediction period as that in the 1st EOP PCC, so as to compare the predictions with those submitted to this campaign in a common and objective way. Using the EOP 05 C04 series as a reference, Figure 3 and Table 2 compare the accuracy of the short-term prediction for 30 days ahead and medium- and long-term predictions of x_p and y_p pole coordinates up 365 days into the future based on WLS+AR (piecewise weights) with the results submitted to the campaign by different methods and techniques. What can be said with the information from Figure 3 is that LS extrapolation+AR prediction by ID Kalarus 061 and ID Kosek 051 is the most accurate technique for medium- and long-term x_p, y_p predictions as demonstrated in the 1st EOP PCC, while pure AR prediction (Zotov 091) and neural networks (Zotov 093) are more accurate methods to predict short-term x_p and y_p in addition to LS extrapolation+AR prediction.

Figure 3.

Figure 3 and Table 2 indicate that the accuracy of the short-term x_p and y_p predictions up to 30 days obtained by WLS+VAR is comparable with that by LS extrapolation+AR prediction (Kalarus 061 and ID Kosek 051) as well as by pure AR prediction (Zotov 091), and better than the prediction accuracy by neural networks (Zotov 093) participating in the campaign. Moreover, the accuracy of the medium- and long-term predictions up to 365 days is inferior to the best presently available methods in the 1st EOP PCC, LS extrapolation+AR prediction by Kalarus 061 and ID Kosek 051. However, the predictions of both x_p and y_p are more accurate than those submitted to the 1st EOP PCC computed by the other methods and technologies. The MAE values of the x_p and y_p predictions up to 365 days by WLS+VAR are less than 29 mas and 33 mas respectively, while the MAE values of the 365-day horizon are more than 37 mas and 40 mas respectively for the x_p and y_p predictions by the participants except Kosek 051 and Karalus 061.