Polar motion prediction using the combination of SSA and Copula-based analysis

The real-time estimation of polar motion (PM) is needed for the navigation of Earth satellite and interplanetary spacecraft. However, it is impossible to have real-time information due to the complexity of the measurement model and data processing. Various prediction methods have been developed. However, the accuracy of PM prediction is still not satisfactory even for a few days in the future. Therefore, new techniques or a combination of the existing methods need to be investigated for improving the accuracy of the predicted PM. There is a well-introduced method called Copula, and we want to combine it with singular spectrum analysis (SSA) method for PM prediction. In this study, first, we model the predominant trend of PM time series using SSA. Then, the difference between PM time series and its SSA estimation is modeled using Copula-based analysis. Multiple sets of PM predictions which range between 1 and 365 days have been performed based on an IERS 08 C04 time series to assess the capability of our hybrid model. Our results illustrate that the proposed method can efficiently predict PM. The improvement in PM prediction accuracy up to 365 days in the future is found to be around 40% on average and up to 65 and 46% in terms of success rate for the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hbox{PM}}_{x}$$\end{document}PMx and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hbox{PM}}_{y}$$\end{document}PMy, respectively.


Introduction
Polar motion (PM) describes the movement of the Earth's rotation axis w.r.t the Earth surface. The study of PM provides valuable information for studying many geophysical and meteorological phenomena (Barnes et al. 1983;Wahr 1982Wahr , 1983Mathews et al. 1991;Gross et al. 2003;Chen and Wilson 2005;Gross 2015;Seitz and Schuh 2010;Schuh and Böhm 2011).
Since the 1960s, highly accurate PM coordinates can be obtained by different space geodesy techniques. These techniques include: Satellite Laser Ranging (SLR) (Coulot et al. 2010), Lunar Laser Ranging (LLR) (Dickey et al. 1985), Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS) (Angermann et al. 2010), Global Navigation Satellite Systems (GNSS) (Dow et al. 2009;Byram and Hackman 2012), and very-long-baseline interferometry (VLBI) (Schuh and Schmitz-Hübsch 2000;Nilsson et al. 2010Nilsson et al. , 2011Nilsson et al. , 2014. Accurate real-time PM is needed for high-precision satellite navigation and positioning and spacecraft tracking (Kalarus et al. 2010;Stamatakos 2017). However, the PM is not provided in real time due to the complexity of the measurement model and data processing; PM coordinates are available with a delay of hours to days (Bizouard and Gambis 2009;Schuh and Behrend 2012). Therefore, it is essential to predict the PM parameters precisely.
Different methods and means have been investigated and applied for PM prediction such as least squares (LS) collocation (Włodzimierz 1990), spectral analysis and LS extrapolation method (Akulenko et al. 2002), LS extrapolation of a harmonic model and autoregressive (AR) prediction (Kosek et al. 1998(Kosek et al. , 2007Xu et al. 2012), wavelets and fuzzy inference systems (Akyilmaz and Kutterer 2004;Akyilmaz et al. 2011) modeling and forecasting excitation functions (Chin et al. 2004), Kalman filter with atmospheric angular momentum forecasts (Freedman et al. 1994), and artificial neural network (ANN) (Schuh et al. 2002;Kalarus and Kosek 2004). The Earth orientation parameters prediction comparison campaign (EOP PCC) took place within (2005)(2006)(2007)(2008)(2009), and the results demonstrate that there is no particular method superior to other for all prediction intervals (Kalarus et al. 2010). Among these methods, the combination of LS and AR process is considered to be one of the most effective for PM prediction (Kalarus et al. 2010). The mentioned combination method achieved reasonable results for short-term forecasting. However, due to the complexity of the PM excitation model, it is not able to reproduce the time variation of the periodic terms that influence the long-term predictive accuracy of PM. Consequently, a new prediction method is required that could bring us significantly closer to meeting the accuracy goals pursued by the Global Geodetic Observing System (GGOS) of the International Association of Geodesy (IAG), i.e., 1 mm accuracy and 0.1 mm/ year stability on global scales in terms of the ITRF defining parameters (Plag and Pearlman 2009). Therefore new techniques or a combination of the existing methods need to be investigated for improving the efficiency of the predicted PM considering the time variation of the periodic terms and the trend. In this study, we examined the combination of singular spectrum analysis (SSA) and Copula-based analysis to predict PM. SSA is not constrained by the assumptions of using predetermined functions such as sine wave as the base; it rather exploits data-driven base functions for extracting fundamental components of the time series and applies a classification method to explore the relationship between the derived elements (Broomhead and King 1986;Vautard et al. 1992;Zotov 2005). The Copula method operates linear and nonlinear dependency between variables, and it is a potent and efficient tool for dealing with multidimensional data and modeling the relationship between parameters (Joe 1997). The combination of SSA and Copula-based methods will be applied for the first time as a novel stochastic tool for PM determination.
PM is the sum of two statistically independent parts: trend and undulation. This hybrid model consists of a deterministic annual and the Chandler component as well as long-term lower-frequency parts which are estimated by SSA. The difference between the deterministic solution and the PM data is then used in a Copula-based model to predict stochastic processes. Then, the final PM prediction is a combination of the deterministic prediction (derived from the SSA solution) and the stochastic prediction (obtained from the Copula solution). To this end, first, the time series of PM from EOP 08 C04 were analyzed, and the trend is modeled and separated by SSA. Then, a Copula prediction model is made based on the SSA-separated time series.
Finally, the accuracy of the proposed combined method is verified through different sets of PM prediction tests.

Methodology
In this study, we developed and explored the integration of Copula-based analysis and SSA for precisely predicting PM.

Singular spectrum analysis
To maximize the prediction performance, we need a mathematical tool to retrieve all time-correlated information from the time series. As a matter of fact, the existence of excitations of PM can profoundly affect the forecasting procedure, particularly in longer intervals. Therefore, the exploitation of efficient techniques is crucial to minimize the risk of having gross errors.
SSA is a nonparametric spectral estimation method which can be used for decomposing a time series into the sum of interpretable components, e.g., trend, periodic components, and noise, without a priori assumption about the constituent components (Golyandina et al. 2001).
SSA is able to remove redundancies and groups uncorrelated information into informative empirical functions which can reveal main aspects of the time series. The mentioned functions are used as bases of a subspace in which the time series is a member of and can be exploited for modeling the time series in a desired level of details. Therefore, the model can simulate the future entries of the time series using these base functions.
The SSA method for trend extraction can be succinctly expressed as two stages: Decomposition First the time series is embedded in an L-dimensional vector space. The outcome of this stage will be a trajectory matrix ( X ) which consists of L rows. The matrix has been simply formed using L-element lagged vectors taken from the time series by sliding a window of size L. (1) Having the trajectory matrix formed, singular value decomposition (SVD) is applied to factorize X in the form of U V T in order to retrieve its principal components (PC).
where U and V are the left and right singular vectors, respectively, and is a diagonal matrix consisting of singular values of X which reflect the importance of each corresponding pair of left-right singular vectors. The decomposition step can be performed using calculation of eigenvalues and eigenvectors of the matrix S = XX T .
Let 1 ≥ 2 ≥ · · · ≥ L ≥ 0 denote diagonal entries of (the eigenvalues of S ) and U 1 , U 2 , . . . , U L indicate the corresponding eigenvectors of S which are also called empirical orthogonal functions (EOF) of X . The right singular vectors of X are eigenvectors of X T X calculated by: Now, the trajectory matrix can be written as: Reconstruction This stage aims to rebuild the time series using the reconstructed version of trajectory matrix. So, a subset of A = {X 1 , X 2 , . . . , X d } can be chosen for reconstruction of the trajectory matrix. The choice of PCs of X defines how smooth would be the reconstructed version of the time series and how much detail of the original time series would be captured. Having a proper selection of PCs, a representative trend is extracted by applying diagonal averaging to the reconstructed trajectory matrix ( X ) . Let L < K , and then, the trend of the time series G = (g 1 , g 2 , . . . , g N ) is:

Copula-based analysis
There is a well-introduced method called Copula that can be applied for polar motion modeling, estimation, and prediction. The word of Copula is a Latin noun that means a link or tie. The Copula method exploits linear and nonlinear dependency between variables. It is a potent and efficient tool for dealing with multi-dimensional data and modeling the relation between parameters based on the marginal distribution functions of the variables (Embrechts et al. 2002). Copula appeared in the context for the first time by Sklar (1959). Sklar's theorem indicates that a Copula function C connects a given multivariate distribution function with its univariate marginal. For bivariate distribution, there is a bivariate Copula C which models the joint cumulative probability distribution function of two variables X and Y based on the marginal cumulative distribution functions F X (x) and F Y (y).
where C describes the joint distribution function F X,Y (x, y).u and v are transformed of X and Y to uniform distribution, respectively. Then, Joe (1997) and Nelsen (2007) developed the idea of the Copula. For many years, the Copula method has been used for modeling the dependence structure between random variables in different types of studies, e.g., economics (Rachev and Mittnik 2000;Patton 2006Patton , 2009, biomedicine (Wang and Wells 2000;Escarela and Carriere 2003), hydrology (Bárdossy and Li 2008;Bárdossy and Pegram 2009;Verhoest et al. 2015), meteorology (Laux et al. 2011;Vogl et al. 2012), hydro-geodesy (Modiri et al. 2015). A brief introduction to the concept of copula function is given in the next subsections.

Characteristic of Copulas
In the bivariate case, a Copula is represented as a function C from [0, 1] 2 to [0, 1] so that ∀u, v ∈ [0, 1] (Genest and Rivest 1993;Jaworski et al. 2010): Copula is an increasing function. It implies that Copula is a continuous function: The Copula density is computed by differentiating Copula cumulative distribution function.

Empirical Copula
The empirical Copula is an estimator for the unknown theoretical Copula distribution, and it is defined in the rank space as follows (Genest and Rivest 1993;Genest and Favre 2007;Laux et al. 2011): where, (r 1 ), (r 2 ) . . . , (r n ) denote the pairs of ranks of the variable (x 1 ), (x 2 ), . . . , (x n ), (s 1 ), (s 2 ) . . . , (s n ) denote the pairs of ranks of the variable (y 1 ), (y 2 ), . . . , (y n ), n is the length of the data vector, 1(...) is the indicator function. If the condition is true, the indicator function is equal to 1. Otherwise, the indicator function is equal to 0.

Archimedean Copula
A number of Copulas can be estimated directly with the simple form. They are named Archimedean Copulas. An Archimedean Copula can be presented in the following form: where θ is the Copula parameter and the function φ is the generator of the Copula with the following characteristics (Nelsen 2007): and φ −1 is defined by There are three commonly used Archimedean Copula which are explained as follows and will be investigated in this study (see Table 1).
(1) Clayton Copulas The generator of the Clayton Copula (see Fig. 1) is given by Therefore, the cumulative distribution function (CDF) for Clayton Copula is defined as (Clayton 1978): where θ is restricted on the interval [−1, ∞) . If θ = 0 , it shows the independence case and when θ → ∞ , indicate high dependency in the lower rank space.
(2) Frank Copula The generator of the Frank Copula (see Fig. 2) is given by The CDF for Frank Copula is given by (Joe 1997;Lee and Long 2009) Frank Copula allows to model data with positive and negative dependency. The large positive and negative θ indicate high dependency, and θ = 1 implies total independence. The Frank Copula is a suitable method for two data sets with the same dynamic characteristics (Rodriguez 2007). (

3) Gumbel Copulas
Gumbel Copula (see Fig. 3) is famous for its ability to capture strong upper tail dependence and weak lower tail dependence. Gumbel Copula is used to model asymmetric relationship in the data (Jawor- The Gumbel Copula generator is written as: The CDF for Gumbel Copula is given by (Nelsen 2007) The Copula parameter θ is on the interval [1, +∞) . If θ is equal 1, Copula shows independence. When θ → ∞ , the Gumbel Copula indicates high dependence between the random variables.

Copula parameter estimation
The widely used estimation method for the Copula parameter is the maximum likelihood (ML) estimation methodology (Joe 1997). The Copula parameters in this study are derived from ML estimation. The canonical maximum likelihood estimation (CLME) and inference for margins estimation (IFME) are two methods for estimation of the Copula parameter (Joe and Xu 1996). For both methods, the first step is marginal distribution estimation. Then, a pseudo-sample of the transformed observation is used to estimate the Copula parameter. In the IFME method, the theoretical marginal distribution parameters are estimated, and in the CMLE the univariate marginals are the empirical distribution functions (Giacomini et al. 2009). It is assumed that the sample data (X 1 , X 2 , X 3 , . . . , X n ) are n independent and identically (20)  distributed (iid) random variables. These data are transformed into uniform variates (r 1 , r 2 , r 3 , . . . , r n ).
Let c(r 1 , r 2 , r 3 , . . . , r n ) be the density function of Copula C(r 1 , r 2 , r 3 , . . . , r n ; θ) , and let θ be the Copula parameter which is estimated by maximizing the ML equation:

Computation of conditional CDF for Archimedean Copula
In this subsection, the conditional CDF of Clayton, Frank, and Gumbel Copula are computed (Yue 1999;Zhang and Singh 2007;Trivedi et al. 2007). The conditional CDF for Clayton Copula is given by (Joe 1997): and for Frank Copula: The conditional CDF of Gumbel Copula is:

Simulating from Copula-based conditional random data
This subsection provides the essential steps for data simulation using Copula-based conditional random data. The following steps are taken to fit the suitable theoretical Copula function and simulation data (Laux et al. 2011;Vogl et al. 2012).
(2) Compute the marginal distribution F X (x) and F Y (y) of the input data x and y.
(3) Transform data to rank space using the estimated marginal distributions of data with u i and v i in rank space. (4) Compute the empirical Copula to the dependence structure of random variables using the rank-transformed data. (5) Fit a theoretical Copula function C θ (u, v). (6) Compute the conditional Copula function. (7) Sample random data from the conditional Copula CDF.
log c(r 1 , r 2 , r 3 , . . . , r n ; θ) (8) Transfer the sample back to the data space using the inverse marginal.

Error analysis
The mean absolute error (MAE) standard is used in order to evaluate the prediction accuracy. It can be shown as follows: where P i is the predicted value of the i-th prediction, O i is the corresponding observation value, E i is the error, and n is the total prediction number (Willmott and Matsuura 2005).

Data description
In this paper, the PM x and PM y time series (see Fig. 4) are from the International Earth Rotation and Reference Systems Service (IERS) combined earth orientation parameter (EOP) solutions 08 C04 (available at http://hpier s.obspm .fr/eop-pc/analy sis/excit activ e.html). The EOP 08 C04 series is derived from different geodetic techniques, and it is consistent with ITRF 2008. The EOP 08 C04 time series cover the period 1962 to the present. The sampling interval is one day.

Data processing and analysis
In this study, we defined an algorithm for PM prediction which is shown in Fig. 5. The observed PM time series can be split up into two parts. The first part is dealing with periodic effects such as Chandler wobble, annual variation, and influences of solid Earth tides and ocean tides on PM. The SSA is used to model the periodic terms of the PM. Then, the difference between the observed PM and SSA estimated data is modeled by using the Copula-based analysis method. After that, the periodic terms of PM are extrapolated using the SSA a priori model. Also, the anomaly part is predicted using the Copula-based model. Finally, the anomaly solution is added to the SSA-forecasted time series. Therefore, the analysis of the data is divided into two main steps: (1) SSA Periodic Terms Estimation • Singular value decomposition of X, • Selecting a proper group of singular values and corresponding singular vectors, • Reconstruction of X, • Calculation of the trend by applying diagonal averaging to X.
(2) Copula Anomaly Modeling  Therefore, the final PM predicted data is the sum of the results of predicted periodic terms using SSA and predicted anomaly using the Copula-based model.

SSA periodic terms estimation
Window length selection is a crucial step in SSA which has a significant impact on the decomposition of the time series. The appropriate choice for L in a periodic time series with a period T is a window length proportional to the period, meaning that the L / T is an integer. Figure 6 depicts the main periods of PM time series (Golyandina and Zhigljavsky 2013). So, the Chandler period as the main period of both time series would be a reasonable choice. Making the closest choice to the half of the length of the time series (if possible, least common multiple of the Chandler and annual periods) is recommended by Golyandina and Zhigljavsky (2013), but is avoided due to the processing time.
After selection of the window length, the number of singular vectors or empirical functions for reconstruction of the time series should be determined. The goal of this procedure is to find and apply a proper set of constructive components. Most significant periodicities as well as excitation mechanisms are rather low-frequency components and reveal their impact in the first few singular vectors while high-frequency components fall in later  Figure 7 suggests that in order to achieve an accuracy of about 1 mas in polar motion modeling, we need to utilize at least first 70 singular vectors which correspond to using all components with periods more than or equal to 14 days.
Having the window length and the number of singular values determined, we construct the trajectory matrix. As it can be seen in Fig. 8 the data between the year of 1997 and 2003 is used as the training period. The cyan curve is the SSA-reconstructed PM x time series. Prediction of the future entries starts by adding initial guess of future entries to the end of the time series. Then, iteration of the SSA process is done until the result of two successive iterations has a difference less than a certain threshold. This will map the initial values to the estimated periodic terms of the time series. The residual part of the difference between original PM x time series and SSA estimated time series is named anomaly of PM x which has a stochastic behavior. Therefore, the anomaly part will be investigated by Copula-based technique.

Copula anomaly modeling
The anomaly part which is shown in Fig. 8 (lower panel) with dark violet is formed into a matrix with the same window length L. Then, the dependency structure between the rmcolumn i and rmcolumn i+1 is investigated for the whole dataset. Modeling the joint dependence structure with Copulas requires fitting marginal distribution to data. In this study, three univariate distribution functions are considered: extreme value, generalized extreme value, and generalized Pareto distribution (see Table 2). To identify which univariate distribution is the best suitable for both PM x and PM y , the root-meansquare error is estimated and the goodness of fit is examined with the Akaike and the Bayesian information criteria (AIC and BIC). and (28) where k denotes the number of the free parameters in the model. n is the sample size, and B is the maximized value of the likelihood function of the estimated model. The smallest amount of AIC or BIC, respectively, suggests the best fitting model or distribution. After estimation of the parameters by maximum likelihood approach, the AIC, BIC, and RMSE values are calculated for both PM x and PM y distribution. As it can be seen in Fig. 9, the generalized extreme value (black) provides the best fit in comparison with the generalized Pareto distribution function (blue) and extreme value distribution function (green). Furthermore, according to Tables 3 and 4, the result of the AIC, BIC, and RMSE confirmed that the generalized extreme value provides the best fit in both PM x and PM y distribution. Therefore, generalized extreme value distribution was selected in this study.

Estimating empirical Copula
Once the univariate marginal distribution is fitted, the dependence structure between the time series has to be investigated. The first step is to calculate the empirical Copula using Eq. (14). As it can be seen in Fig. 10, there is a scatter plot of two adjacent columns, and it shows a scatter linear dependency structure with the heavy tail. This kind of dependency structure can be correctly modeled using the Archimedean Copula.

Fitting a theoretical Copula function
The next step is fitting a theoretical bivariate Archimedean Copula function with its parameters estimated by maximum likelihood approach. In this study, three different theoretical Copula functions are tested (Fig. 11): Clayton, Frank, and Gumbel Copula. For the three different Copula functions, the goodness-of-fit test, which is based on the Cramer-von Mises statistics, is applied. To evaluate the performance of the Copulas, 1000 values of the test statistics are sampled, and the proportion of values larger than S n is estimated by calculating the corresponding p values. The results based on S n show that the performance of Frank Copula is slightly better than Gumbel and Clayton Copula with less error (Table 5).

365-day-ahead prediction
We utilized 6 years of observed PM time series, from January 1997 to December 2002, for the 365-day-ahead prediction. To verify the reliability of this method, the results were compared with the IERS Bulletin A predictions (https ://datac enter .iers.org/web/guest /bulle tins/-/somos /5Rgv/versi on/6). The IERS Bulletin A contains the PM parameters and the predicted PM for one year into the future, and they are released every seven days by IERS In the current prediction method, the PM prediction was the sum of the LS extrapolation model (including Fig. 9 Marginal distribution's goodness-of-fit test for PM x (left) and PM y (right). Generalized extreme value distribution is the black curve, green shows the extreme value distribution, and the blue curve is generalized Pareto distribution Extreme value (Kotz and Nadarajah 2000) f Generalized extreme value (Hosking et al. 1985) f Location µ scale σ shape ξ Generalized Pareto (Hosking and Wallis 1987) f Location µ scale σ shape ξ the Chandler period, annual, semiannual, terannual, and quarter annual terms), and the AR predictions of the LS extrapolation residuals (Kosek et al. 2007).

Discussion of results
In this study, we demonstrated the PM prediction by combination of SSA and Copula-based analysis method.
Our method is tested based on the hindcast experiments using data from the past. Hence, we have calculated the results of our methods yearly for seven years of the test     Fig. 12 shows the mean value of MAE for each year. In Fig. 12, the Bulletin A solution is shown in black and the SSA predicted data in red. Also, the results of SSA+Copula are displayed by green, blue, and pink for Clayton, Frank, and Gumbel Copula, respectively. Compared to the results from the IERS Bulletin A, the MAE of the predictions produced by the proposed method was smaller in different short-, mid-, and long-term intervals for different cases (e.g., between 1 and 5 mas progression of PM x prediction for different time intervals in 2003). The better prediction performance of the SSA + Copula prediction may have been due to the modeling of the linear change of the Chandler and the annual oscillation amplitudes. Besides, the combination of SSA+Copula improves the SSA solution because of its ability to model the stochastic behavior of the anomaly part of the PM time series. However, the proposed method did not always perform better, especially in cases of long-term prediction where the quality of the results was not as good as we expected (see Fig. 12). This may have been caused by changes of the amplitudes of the periodic terms in this six-year time span where the SSA was not able to capture all features in order to predict more precisely we would have to increase the interval of training time. Figure 13

Conclusions
The improvement in the Earth rotation prediction is a relevant, timely problem, as confirmed by the fact that the International Astronomical Union (IAU) Commission A2, the International Association of Geodesy (IAG), and the IERS have at present two Joint Working Groups on Prediction (JWG-P) and on Theory of Earth rotation and validation (JWG-ThER). According to the United Nations (UN) resolution in 2015, the primary objective of these JWGs is to assess and ensure the level of consistency of earth orientation parameter (EOP) predictions derived from theories with the corresponding EOP determined from analyses of the observational data provided by the various geodetic techniques. Therefore, accurate EOP predictions are essential to avoid any systematic drifts and/or biases between the international celestial and terrestrial reference frames (ICRF and ITRF). The results illustrate that the proposed method could efficiently and precisely predict the PM parameters. As clearly demonstrated, the SSA + Copula algorithm shows better performance for PM x prediction in comparison with the SSA prediction. The Copula-based analysis is fully successful in its aim to increase the accuracy of PM prediction by modeling the stochastic part of the PM and subtracting PM by SSA-reconstructed time series. We suspect the main error contributions come from SSA extrapolation part. So, further investigations about the SSA training time will be required to clarify this issue. Also, SSA + Copula prediction method shows periodic errors, and these errors have a significant impact on the mean (30) Success rate of PM prediction = Number of improvement in the predicted PM Total number of PM prediction × 100   absolute error. Therefore, these occasional errors should be further investigated to have a noticeable progression in the PM prediction accuracy.