Testing for Serial Correlation in Autoregressive Exogenous Models with Possible GARCH Errors

Autoregressive exogenous, hereafter ARX, models are widely adopted in time series-related domains as they can be regarded as the combination of an autoregressive process and a predictive regression. Within a more complex structure, extant diagnostic checking methods face difficulties in remaining validity in many conditions existing in real applications, such as heteroscedasticity and error correlations exhibited between the ARX model itself and its exogenous processes. For these reasons, we propose a new serial correlation test method based on the profile empirical likelihood. Simulation results, as well as two real data examples, show that our method has a good performance in all mentioned conditions.


Introduction
Time series data are frequently encountered in fields like weather forecasting, earthquake prediction, and electroencephalography. How to analyze time series data has been of great interest in statistics for a long time. Many famous models have been developed to study the relationships in time series data. Univariate time-series models include autoregressive models (AR), moving average models (MA), autoregressive moving average (ARMA), predictive regression models (PRM), and so on. Extensions of these models to handle vector-valued data contain, e.g., vector autoregression models; see, for example, ref. [1] for detailed discussions of these models.
Among others, the AR model, firstly proposed by [2], takes the simplest form. It specifies that the output variable only depends linearly on its own previous values and the random error term. But in many situations, auxiliary information is available and can be treated as covariates for modeling. Furthermore, the inclusion of auxiliary variables can improve the estimation efficiency. Hence, the so-called autoregressive exogenous (ARX) model was developed in previous literature to incorporate this benefit. The ARX model takes the following form: where Y t denotes the response, I(B) = 1 − ∑ p i=1 a i B i , and B denotes the backshift operator with B l Y t = Y t−l for l ≥ 1. Suppose X t is the covariate that contains useful auxiliary information. U t s are model errors with means of 0 and finite variances. Without confusion, we assume that {X t } follows the first-order vector autoregressive (VAR) model: where V t denotes the vector-valued random errors. α, β, µ and A are unknown parameters.
Moreover, the ARX model (1) has been widely adopted to analyze time series data in many fields. For example, in finance, as narrated in [3], an illustration of (1) can be that Y t reflects a change in an asset's price, and X t ∈ R d are lagged variables related to asset prices. A more specific example is that [4] employed the VARX model to investigate the relationship between the closing price of HRUM energy and PTBA (endogenous variable) and the exchange rate (exogenous variable). Ref. [5] applied the ARX model in the generalized space and used the international crude palm oil (CPO) prices as exogenous variables to predict the export volume of CPO. Additionally, in the environmental area, Ref. [6] analyzed hourly ozone data collected routinely at several monitoring sites in Austria using different ARX models to complete a pollution assessment. Furthermore, Ref. [7] used ARX-GARCH models to forecast air quality levels using daily data from the monitoring stations of 16 cities/counties in southeast China.
Some pre-test procedures are needed when specifying time series like (1), such as testing the existence of unit root in series Y t . To achieve this, there is a useful technique of replacing I(B) in (1) by a specific linear filter P(B) = 1 − φB − ∑ p−1 i=1 ψ i ∆ i , in which ∆ denotes the difference operator. Note that φ = ∑ p i=1 a i and ψ i = − ∑ p−1 j=i a j+1 for i = 1, · · · , p − 1. Hence, we can rewrite (1) in the following form (3). This trick is similar to the so-called Dickey-Fuller reparameterization [8,9], which is proposed for constructing a Dickey-Fuller unit root test.
The reparameterization trick given in (3) is most notable for its convenience of using φ instead of ∑ p i=1 a i to represent the stationarity property of the endogenous structure of Y t . Hence, we choose (3) as the expression of our ARX model in the following discussions.
Owing to its wide applications, many researchers have focused on the theory of model specification of ARX models, including both parametric and nonparametric methods. To name but a few, ref. [10] proposed a least squares method to estimate the parameters in linear and nonlinear ARMAX models and discussed the asymptotic properties of these estimators. Ref. [11] employed a local polynomial fitting scheme incorporated with projections to obtain nonparametric estimation of additive nonlinear ARX time series. Ref. [12] introduced a quasi-likelihood method in estimating a censored ARX model and proved the quasi-likelihood estimation computationally efficient; see, e.g., [12] and references therein for more details on this topic.
Despite the fact that numerous studies have investigated different methods to specify and estimate the ARX model, little literature has discussed the issue of testing its serial correlation. Most studies assume that errors are of no autocorrelations; however, this assumption can be often violated, and it may negatively affect our further inference. For instance, the endogeneity problem, as a common cause of serial correlation, may destroy the consistency of the regression coefficients, leading to misspecification of the model, and damage the explanatory properties in real applications. Therefore, it is of importance to check for serial correlation of the random model errors.
There are many methods of testing serial correlations based on the least squares (LS) residuals. Along this direction, Ref. [13] proposed the first corresponding test procedure, the Durbin-Watson (DW) test. Unfortunately, it only works in testing the first-order autocorrelation. To deal with this limitation, the portmanteau test constructed by Q statistics gained great popularity; this test is also known as the Q test. The two most famous forms of Q tests are the Box-Pierce (BP) and Ljung-Box (LB) tests, proposed in [14,15], respectively. However, much evidence shows that Q tests suffer from the size distortion issue for models with AR structures; see, e.g., [16], who argued against the validity of Q tests by showing that the asymptotic chi-square distribution under the null hypothesis of no autocorrelations relies on the strictly exogenous nature of all regressors to the error terms. Although some further works like [17,18] expanded the adoption of portmanteau tests to ARMA models by developing the asymptotic results of Q statistics, their techniques cannot apply to the ARX models since the combination of AR and exogenous structure leads to an increase in complexity.
Another alternative is to use the plug-in (PI) empirical likelihood method, which follows the approach in [19], or the Breusch-Godfrey (BG) test from [20,21]. However, it is worth noting that in practice, such as the financial data discussed in [3], the U t in (3) and V t in (2) might be correlated. This correlation could result in a biased LS estimator when facing a finite sample size. This means that the plug-in empirical likelihood method fails to work in such a situation. For the BG test, although it seems to perform stably against the possible correlation existing between U t and V t (as will be reported in our simulations), there are still some more scenarios for which the BG test cannot be applied. For instance, the heteroscedasticity mentioned above can violate its disturbance assumptions. Consequently, it not only leads to inefficient parameter estimations, but also breaks the conditions for the application of the BG testing method that relies on LS residuals.
Note that the ARX models are often used in financial data analysis, in which volatility clustering often occurs. Ref. [22] firstly noted this phenomenon. From his description, volatility clustering is defined by the idea that large changes tend to be followed by large changes of either sign, and small changes tend to be followed by small changes. A number of following studies focused on the impacts of its existence in real applications; see, e.g., discussions in [23][24][25]. To address such a feature, a family of widely adopted models named autoregressive conditional heteroscedasticity (ARCH) and generalized autoregressive conditional heteroscedasticity (GARCH) have been proposed, respectively, by [26,27]. With the conditional serial dependence structures, known as conditional heteroscedasticity, economists have found a more appropriate approach to modeling data than following the obviously false assumption of homoscedasticity.
It is interesting to find that when dealing with GARCH errors, many existing serial correlation tests, including the LB, BG, and plug-in empirical likelihood method, suffer from size distortion; see our simulation in the sequel.
To deal with this problem, we construct a new test for the ARX model based on the profile empirical likelihood method motivated by [28]. It turns out that the proposed test for serial correlations performs robustly in the heteroscedastic as well as correlated U t and V t situations.
The rest of this paper is organized as follows. In Section 2, we present the serial correlation test and the main results. Section 3 reports the finite sample performance of the proposed testing statistic. Section 4 further applies the test to financial data and environmental data. Section 5 concludes the whole paper. Detailed proofs of the main results can be found in Appendix A.

Methodologies and Main Results
Assume that {Y t , X t−1 } n t=1 are generated from (2) and (3). Our aim is to construct a test method to check whether serial correlations exist in {U t } n t=1 . The hypothesis of interest is as follows: i.e., to check whether there exist q-th serial correlations in the model errors, where γ k = E(U t U t−k ) and γ = (γ 1 , · · · , γ q ) for some positive integer q.
Since the definition of γ involves the expectation, a straightforward idea is to construct some testing statistics for H 0 by using the plug-in (PI) empirical likelihood method; that is, plugging the LS residuals of U t in the constraints of empirical likelihood ratio function. The details are as follows: define θ = (α, φ, ψ , β ) , and denote the LS estimate of θ as θ := ( α, φ, ψ , β ) ; then, the residuals { U t } n t=p+q+1 can be easily obtained. Write t 0 = p + q + 1 and N = n − p − q. Then, the PI method is obtained by calculating the following EL ratio function as Following a similar proof to that of [19], it is possible to show that −2 log L(0) is asymptotically chi-square distributed once {U t } is a sequence of independent and identically distributed errors under some mild conditions. Hence, for given observations {Y t , X t−1 } n t=1 , one may reject the null hypothesis once −2 log L(0) > χ 2 q (1 − τ) at the significance level τ ∈ (0, 1), where χ 2 q (1 − τ) denotes the (1 − τ)-th quantile of the chi-squared distributed variable with q degree of freedoms.
However, although this method is computationally fast and can avoid the variance estimation, it still suffers from a strong undersized distortion. As an improvement, we may treat θ as redundant parameters, and, similar to [28], construct a series correlation test by using the profile empirical likelihood (PEL) method. In Then, similar to [19], an empirical likelihood function for the unknown θ and γ can be defined as follows: where Z t (θ, γ) = ( Z t,1 (θ, γ), Z t,2 (θ, γ), · · · , Z t,p+d+q+1 (θ, γ)) with Here, X t−1,l denotes the l-th component of X t−1 for l = 1, 2, · · · , d.
Since we are interested in γ, we may eliminate the effect of θ by the profile method, as did [28]. The resulting profile empirical likelihood function L(γ) for γ can be defined as follows: To facilitate studying the limiting distribution of L(γ), we need to first specify the following assumptions: • A1. {X t } are a strictly stationary sequence with the initial value X 0 being a constant vector.
Hereafter, F t denotes the information set available at time t ∈ {1, 2, · · · , n}. Based on these assumptions, we can obtain our main result of Theorem 1.
as n → ∞, where χ 2 q denotes the chi-squared distributed variable with q degrees of freedom, and d −→ denotes the convergence in distribution.
Theorem 1 is desirable because there is no need to estimate the asymptotic variance, which is difficult especially when the model errors follow the GARCH process. In practice, we may reject the null hypothesis H 0 if at the significance level τ ∈ (0, 1), analyzing the random observations {Y t , X t−1 } n t=1 .

Simulation Studies
In this section, we conduct simulations to investigate the performance of our proposed testing method, i.e., the PEL test for the ARX model. Meanwhile, we compare the testing with the other two extant testings, including the PI and BG tests. Firstly, we set two processes in generating U t , and let {e t } n t=1 be an independent identically distributed standard normal series with e 0 = 0. The two forms of U t are, respectively, The autocorrelation of U t in (9) is from the simple AR(1) structure generated by nonzero λ, and (9) as a GARCH(1,1) process depicts a relatively complex heteroscedastic structure. When λ = 0, the results refer to the test sizes. To examine the local powers, we then set λ = 1/ √ n, · · · , 4/ √ n. The correlations between U t and V t are designed by ρ = (ρ 1 , · · · , ρ d ) . This vector determines the variance matrix of the joint normal distribution of U t and V t in which we assume For convenience, we consider the ARX model in the following form: The intercept and slope coefficient of X t were set as constants; that is, µ = 0.01 and β = 2.14. The slope of difference term ∆Y t−1 is ψ = 0.13. There are two settings of α, 0 and 0.01, to imitate the presence or absence of an intercept. For both φ and A, we considered two conditions, 0.6 and 1, in which 0.6 refers to the stationary condition, and 1 refers to the unit root process. The covariance of U t and V t , i.e. ρ, was set to 0, 0.4, 0.8, representing no correlation, light correlation and heavy correlation, respectively. The sample size n was 200 and 800. For each case, we generated 10,000 random samples.
Within the programming of profile empirical likelihood function in the PEL and PI tests, we used R functions "el.test" in the package "emplik". Meanwhile, the R function "nlm" was applied to complete the profile optimization in PEL. We used the R function "bg.test" to obtain the results of the BG test. Then, we reported the comparisons of empirical sizes and powers of PEL, PI and BG in Tables 1 and 2 and Figures 1-6.
The results of the empirical sizes and powers of the AR structural U t are shown in Table 1. As we expected, our proposed PEL test performs well in all assumed conditions. While in some settings, the simulation with smaller sample size n = 200 shows a slight oversized distortion, e.g., when both φ and A are 1 (nonstationary AR structures of X t and Y t ), as the sample size increases to n = 800, the sizes all converge to the presumed τ. Note that there are some methods in literature that can be used to adjust the size of the empirical likelihood method when the sample size is relatively small. One may improve the size performance by the Bartlett correction; see, e.g., [29] for details. Furthermore, from Table 1, we can see that with the AR structural U t , the BG test also works well, and it is efficient regardless of the presence of a strong correlation between U t and V t . However, for the PI tests, the results indicate that the corresponding sizes suffer from obvious undersized distortions due to the autoregressive structure. Similar to portmanteau Q tests such as BP and LB, the PI test has the constraint that the regressors should be exogenous to avoid the disturbance from the plug-in estimators of the autoregressive coefficients, so PI is not suitable for ARX models.
Results of empirical sizes and powers of the GARCH structural U t are shown in Table 2. Different from the cases with AR structural U t , the sizes of the BG test are highly deviated from the presumed τ, which means the heteroscedastic structure of random disturbances completely destroys the validity of the BG test. However, the PI test works well in some settings under GARCH-type errors. This is also shown in Figure 4. For instance, when U t is not or weakly correlated with V t , such as ρ = 0 and 0.4, the results of the PI test are similar to the size values τ given in advance. However, when the correlations increase or even close to 1, which often occurs in practice, the PI test fails and shows a strong undersized distortion. Compared to both the BG and PI tests, the PEL method shows its robustness. Note that only with the smaller sample size n = 200, the sizes of PEL tests show notable overrejections, but these distortions mostly disappear with the increase in sample sizes; see the results of n = 800. Except in a special condition when X t and Y t are both nonstationary(φ = A = 1), and U t is extremely strongly associated with V t , we can find a non-negligible overrejection by comparing the sizes of PEL with the predetermined sizes τ, even at a sample size of 800.
Part of the simulation results are plotted in Figures 1-6. Comparing the first two figures which show the results of the PEL test in AR and GARCH errors, the increasing speeds of rejection rates of GARCH errors vary with λ and are slower than those of AR errors. Nevertheless, the PEL test remains effective with GARCH errors. From the next two figures, we can see a good performance of PI test in limited endogenous GARCH errors, but it fails when the endogenous errors are too strong. The results of PI testing in AR-type errors also indicate apparent distortions. Finally, in the last two figures, we can find that the BG test performs well in AR errors but becomes invalid with respect to GARCH errors, as discussed above.

Real Applications
This section is devoted to applying the proposed PEL test to examining the error autocorrelations in two real applications from the financial and environmental fields. We then compare the findings with the results from the other two testing methods, the PI and BG methods.

A Financial Example
In last few decades, numerous studies have investigated whether stock returns can be predicted by financial and economic variables, such as the dividend-price ratio, the earnings-price and other measures of the interest rate. For instance, ref. [30] implemented a test of predictability on U.S. equity data. They indeed found reliable evidence for predictability with the earnings-price ratio, the T-bill rate, and the yield spread on stock returns, and weaker evidence for predictability with the dividend-price ratio. In addition, ref. [17,18] concluded that an ARMA(1,1) model of the stock returns cannot be rejected. Therefore, the stock returns at lag 1 might be a potential explanatory variable.
In this subsection, we apply the new proposed PEL method to revisit U.S. equity data analyzed by [30]. More specifically, two series of stock returns are involved in the empirical study, monthly S&P 500 index data from Global Financial Data and monthly CRSP value-weighted index data from the Center for Research in Security Prices (CRSP). See Figure 7.
Similar to [30], stock return indexes were regarded as the predicted variable, and three different sample periods were considered : 1927-2002, 1927-1994 and 1952-2002. The reasons for splitting the sample period are that around 1994, valuation ratios drifted to historical lows, making the process more nonstationary, and the nature of U.S. interest rates was changed by the Fed's pegging rate policy after 1952. For the first two periods, 1927-2002 and 1927-1994, we only included two exogenous variables, log earnings-price (lep) ratio and log dividend-price (ldp) ratio. The earnings-price ratio was computed as a moving average of earnings over the past ten years divided by the current price, and the dividend-price ratio is computed as dividends over the past year divided by the current price. For the sample period 1952-2002, this allowed us to add two additional exogenous variables, the treasury-bill (tbl) rates and the long yield (lty) spread. The short rate was the one-month treasury-bill rate, and the long yield spread was Moody's seasoned Aaa corporate bond yield. Finally, the stock returns at lag 1 were treated as another predictor among all sample periods.  4 1927 1936 1945 1954 1963 1972 1981 1990 1999 Time Return Index CRSP SP500 Figure 7. Tendencies of two stock return indexes.
Firstly, we fitted the same predictive regression models of the monthly data as [30] for the three different sample periods and applied PEL testings to examine serial correlation. The p-values are reported in Table 3. Even though [30] provided some evidence about the predictability of log earnings-price ratio on stock returns for the subsample periods 1927-2002 and 1927-1994, our PEL testings indicate that we should reject the null hypothesis of no higher-order serial correlations at the 10% level of significance. Therefore, it might be necessary to add further model correction steps when forecasting the stock returns. Next, besides the multiple econometric variables mentioned above, we also considered the stock returns at lag 1 as predicting variables in ARX models. Then, PEL and BG testings were applied for diagnostic checking, respectively. All results are reported in Table 4. From Table 4, we cannot reject the null hypothesis at the 5% significance level based on the results from PEL testings for all sample periods. The PI testings indicate the same conclusions as the PEL testings at the 5% level of significance. However, the corresponding p-values seem to be extremely huge. This could be because the high correlations between U t and V t make the PI testings invalid: the Pearson's correlation coefficients from all three periods are greater than 0.95. For the BG testings, the results provide the opposite conclusions when looking at sample periods 1927-2002 and 1927-1994. The conclusions show that the ARX models are not adequate. However, for the subsample period 1952-2002, BG testing gives evidence that we should not reject the null hypothesis of no serial correlation at the 5% level of significance. The reason for the different evidence is that BG testing could be invalid when the dataset suffers from heteroscedasticity problems, like the sample periods 1927-2002 and 1927-1994 (see Figure 8). This is consistent with our simulation studies. Furthermore, compared to the PEL testing proposed in this paper, the BG testing might be more vulnerable to over-fitting. Significance levels: * p < 0.1, ** p < 0.05, *** p < 0.01.

An Environmental Example
As mentioned before, ARX models are also popular in the area of environmetrics. A typical usage is in pollution assessment. In this subsection, we use the ARX model to analyze daily ozone data in 2016, which are collected by a monitoring station in Texas, USA. See Figure 9. The ozone dataset can be downloaded from United States Environmental Protection Agency (EPA). Meteorological conditions can strongly affect the rates and completeness of the reactions producing ozone, as well as the subsequent transport and depletion; see [31]. Therefore, in line with [6], the set of exogenous variables includes air temperature, wind speed, extraterrestrial radiation and relative humidity. The meteorological data were collected by the USBA-ARS Conservation and Production Research Laboratory located in Bushland, Texas.
The ozone dataset contains some missing values from 27 June 2016 to 10 July 2016. Rather than deleting these dates or simply using means to impute, we split the dataset into two series: 1 January 2016-26 June 2016 and 11 July 2016-31 December 2016. In [6], they analyzed hourly data points and suggested lags greater than 48 for the autoregressive part. In the current study, the dataset is in daily frequencies, and so we chose ozone at lag 1 and lag 2 as predictors as well. Then, we applied the three different testings, PEL, PI and BG testings, on the specified ARX models separately. All results are reported in Table 5. From Table 5, at least for the lower-order serial correlations, there is no reliable evidence for the presence of serial correlation for the subsample period 1 January-26 June at the 5% significance level. This means the ARX model is adequate. However, for the subsample period 11 July-31 December, the results of PEL testing indicate that we should reject the null hypothesis of no existence of error autocorrelations at the 5% significance level. The findings from PEL testings are consistent with the findings of [6] that global ARX models might be not adequate and it is necessary to consider other advanced models, such as ARXd and ARXr models. However, for the two subsample periods, all results of PI testings provide evidence that we should not reject the null hypothesis at the 5% significance level. This indicates the global ARX models are adequate, which contradicts the findings of previous studies. We then conducted the Pearson's correlation tests and found that we should reject the null hypothesis of the zero correlations between U t and V t at the 5% significance level. Therefore, it is confirmed that the PI method are prone to the underfitting problem when U t and V t are correlated. When looking at the results from BG testings, we can see more sophisticated results. When q = 2, it indicates we should reject the null hypothesis of no serial correlation at the 5% level of significance. However, as q increases to greater than 4, the serial correlation seems to disappear at the 5% level of significance. Significance levels: * p < 0.1, ** p < 0.05, *** p < 0.01.

Conclusions
In this paper, we develop a new test procedure for serial correlation in ARX models based on the profile empirical likelihood method in [28] and propose that the testing procedure holds even when the innovations of disturbances are heteroscedastic or the error series of ARX models and the exogenous processes are correlated. In the simulation studies, the PEL testing performs well in these conditions. Furthermore, we analyze two real data in financial and environmental domains by specifying ARX models and then use the PEL testing method for model diagnosis. The first application shows that when using the PEL testing method, the first-order autoregressive models with exogenous variables, such as log earning-price and log dividend-price, are useful and appropriate in interpreting stock returns. In the second application, we examined whether the first-order or second-order autoregressive model with four exogenous variables, including air temperature, wind speed, extraterrestrial radiation and relative humidity, is appropriate for forecasting the ozone concentrations. The results indicate the pollution assessment could more complex and we should consider more advanced models.
In addition, our simulations and real applications mentioned possible nonstationary models, including unit root or nearly unit root processes. The discussions about the related theory and real examples can be seen in [32][33][34], etc. However, since (A4) in Lemma A2 does not hold, as (A5) does not converge in probability when φ = 1, the rigorous proofs of the chi-squared distribution of proposed EL statistic in nonstationary situations are not given. We seek to address this issue in the future.
Proof of Lemma A1. We define U t = U t (θ 0 ) for short, where t = t 0 , · · · , n. To show the asymptotic result of (A1), we firstly expand the first element of Z t (θ, γ 0 ) as the following form: so we can get that as n → ∞, By assumption A4 and Markov inequality, max t 0 ≤t≤n |U t | = o p ( √ n) is obtained. As for the items of max t 0 ≤t≤n |Y t−1 |, max t 0 ≤t≤n ∑ p−1 i=1 |∆Y t−i | and max t 0 ≤t≤n X t−1 , we can also prove that they are of o p ( √ n). We take max t 0 ≤t≤n |Y t−1 | as an example to display. For arbitrary > 0 and some δ > 0, as n → ∞, The results of two other items are obtained similarly. Thus, max t 0 ≤t≤n sup B | Z t,1 (θ, With this result, the further proofs of max t 0 ≤t≤n sup B | Z t,j (θ, γ 0 )| = o p ( √ n) for j = 2, · · · , p + d + 1 are trivial and similar; we omitted them. Meanwhile, those of the last q components are shown below. That is, when k = 1, · · · , q, The details of proving max t 0 ≤t≤n |B i | = o p ( √ n) for i = 1, 2, 3 are included in that for i = 4, so we only expand B 4 ; that is, By the aforementioned techniques, we can find that the maximum of each term from t in t 0 to n in the expansion of B 4 is o p ( √ n). This finishes the proof of (A1). Next, for j = 1, 2, · · · , p + d + 1 as n → ∞, the proofs of are trivial. Thus, we only show the details of the last q components of (A2); that is, when k = 1, · · · , q, we have Moreover, we expand C 1 as Similarly, we can show that C 2 , C 3 = O p (1) as n → ∞. By then, we have finished the proof of (A2).
The proof of (A3) contains no new questionable part that need to be discussed, so we omit it as an analogue of (A2).
where d −→ denotes convergence in distribution, and p −→ denotes convergence in probability. Matrix Σ is defined as E Z t (θ 0 , 0) Z t (θ 0 , 0) . Note that Σ is a positive definite matrix.
Proof of Lemma A2. According to assumption A2, it is easy to check that { Z t (θ, γ)} n t=t 0 is a Martingale difference series with respect to F t−1 . Hence, based on assumptions A1-A5, we can further show (A4) by the Martingale central limit theorem in [35] and the Cramér-Wold device.
The proof of (A5) is trivial based on (A4); we omit the details.
Proof of Theorem 1. Based on the results obtained in Lemma A1 and Lemma A2, we can finally prove the asymptotic chi-squared distribution of our proposed test statistics. As shown in [36], the solution of p t by Lagrange multipliers is , t = t 0 , · · · , n.