B-spline estimation in varying coefficient models with correlated errors

1.   Introduction

Nonlinear phenomena, namely nonnormality, asymmetric cycles, and nonlinear relationships between lagged variables, have been well observed in some classical data sets, such as the sunspot, Canadian lynx, and Australian blowfly data. However, linear ARMA models can not adequately approximate these nonlinear phenomena^[1]. Due to this, nonparametric regression models have found important applications in modeling nonlinear time series^[2,3]. Yet, in the multivariate setting with more than two variables, it is difficult to estimate regression function with reasonable accuracy due to the "curse of dimensionality". To solve this problem, various semiparametric models have been studied. For example, Stone^[4] proposed the additive model, and Gao^[5] investigated a class of partially linear models. In addition, one of the most popular semiparametric models is the varying coefficient model ^[6], whose regression function depends linearly on some regressors and regression coefficients vary with some threshold variables.

In this study, we consider the varying coefficient model of the following form:

where $\bf{X}_{t} = (X_{t1}, \cdots, X_{tp})^{ \top}$ and $\mathit{\pmb{\beta}}(U_{t}) = (\beta_{1}(U_{t}) \cdots, \beta_{p}(U_{t}))^{ \top}$. The functions $\beta_{j}(\cdot), j = 1, \cdots, p$ are assumed to be unknown but smooth. ${ \top}$ denotes the transpose of a vector or matrix. $U_t$ is a univariate random variable, which is called the threshold variable. Both $\bf{X}_{t}$ and $U_t$ can consist of either exogenous variables or lagged variables of $Y_{t}$. $\epsilon_{t}$ is the error term that satisfies $\mathbb{E}(\epsilon_{t}| \bf{X}_{t}, U_{t}) = 0$, and $\boldsymbol{{\epsilon}} = (\epsilon_{1}, \ldots, \epsilon_{n})^{ \top}$. As a generalization of the linear model, the Model (1.1) has attracted a great deal of attention over the past two decades. When $\bf{X}_{t}$ and $U_t$ are some lagged variables of $Y_{t}$, Chen and Tsay ^[7] proposed an arranged local regression procedure for the specification of the Model (1.1), and the consistency result and a recursive algorithm were given. Cai et al.^[1] applied a local linear technique to estimate the time-varying coefficients. The asymptotic properties of the kernel estimators were investigated under the $\alpha-$mixing condition. Nevertheless, as pointed out in ^[8,9], the local smoothing method is computationally expensive because it requires re-fitting at every point where the fitted function needs to be evaluated. In contrast, the advantage of the B-spline estimation method is its computational efficiency. However, it is difficult to establish the asymptotic normality of the spline estimators ^[10,11]. For varying coefficient models, Lai et al.^[12] considered the B-spline to estimate the time-varying coefficients.

The aforementioned work depends on the assumption of independent errors. However, model misspecification, such as omitting relevant variables and wrong function form, may result in correlated errors, as mentioned in ^[13]. From this perspective, the assumption of independent errors is inappropriate. Many authors have studied the topic of non-/semi-parametric regression with correlated errors. Under the assumption that model errors follow an invertible linear process, Xiao et al. ^[14] proposed a modification of the local polynomial estimation for nonparametric regression. Su and Ullah ^[13] considered nonparametric regression with an error process in a nonparametric autocorrelated form. Lei et al.^[15] investigated a semiparametric autoregressive model with ARMA errors. In this study, we propose a global smoothing method based on B-spline for the estimation of the time-varying coefficients. Similar to ^[16] and ^[17], we have relaxed the assumption of independence to model errors, assuming that they can be correlated. That is, when the model errors are independent, the covariance matrix of the errors is $\mathbb{E}(\boldsymbol{{\epsilon}}\boldsymbol{{\epsilon}}^{ \top}) = \sigma^{2} \boldsymbol{{I}}$. Now, we assume that $\mathbb{E}(\boldsymbol{{\epsilon}}\boldsymbol{{\epsilon}}^{ \top}| \bf{X}_{t}, U_{t}) = \boldsymbol{{V}}$. Specially, if model errors follow an AR(1) process, $\epsilon_{t} = \rho \epsilon_{t-1}+e_{t}$, $| \rho| < 1$, $e_{t}\stackrel{\mathrm{i.i.d}} {\sim} (0, \sigma^{2})$,

Our study assumes that the covariance matrix of model errors is positive definite, without requiring a specific form of autoregression, in which case, the scope of application can be further expanded. Certainly, we extend Theorem 1 of ^[12] to the case of correlated errors.

The remainder of this paper has been presented as follows. In Section 2, the estimation procedure for the varying coefficient model with correlated errors has been introduced. Section 3 presents the consistency and convergence rates of the spline estimators and provides an estimation algorithm. In Section 4, we present numerical examples to illustrate the performance of the proposed estimation method. In Section 5, we compare the results of the proposed spline method with the local linear method proposed by Cai et al. ^[1]. The conclusions are presented in Section 6 and the proofs of the main results in the Appendix.

2.   Estimation method

Let $a = \xi_{0} < \xi_{1} < \cdots < \xi_{M+1} = b$ partition the interval $[a, b]$ into subintervals $[\xi_{k}, \xi_{k+1}), k = 0, ..., M$ with $M$ internal knots. A polynomial B spline of order $r$ is a function whose restriction to each subinterval is a polynomial of degree $r-1$ and globally $r-2$ times continuously differentiable. A linear space $S_{r, M}$ with a fixed sequence of knots has a normalized B-spline basis $\{B_{1}(u), \cdots, B_{K}(u)\}$ with $K = M+r$. As in ^[18], the basis satisfies (i) $B_{k}(u)\geq0, k = 1, \cdots, K,$ (ii) $\sum\limits_{s = 1}^K B_{k}(u)\equiv 1$, and (iii) $B_{j}(u)$ is supported inside an interval of length $r/K$, and at most $r$ of the basis functions are nonzero at any given $u$.

Suppose that each coefficient function $\beta_{j}(u)$ in the model (1.1) is smooth; then, it can be well approximated by a B-spline function $\beta^{*}_{j}(u)\in S_{r, M}$ ^[19]. Thus, there is a set of constants $b^{\ast}_{js}, s = 1, \cdots, K$, such that $\beta_{j}(u)\approx \beta^{\ast}_{j}(u) = \sum\limits_{s = 1}^K b^{\ast}_{js}B_{s}(u)$. Different coefficients might be approximated by B-spline with a different number of knots in principle, but for simplicity, we have assumed the same basis for the different coefficients. Let $\boldsymbol{{b}} = (\boldsymbol{{b_{1}}}^{ \top}, \cdots, \boldsymbol{{b_{p}}}^{ \top})^{ \top} = (b_{11}, \cdots, b_{1K}, \cdots, b_{p1}, \cdots, b_{pK})^{ \top}$, $\boldsymbol{{Z}}_{t} = (X_{t1}B_{1}(U_{t}), \cdots, X_{t1}B_{K}(U_{t}), \cdots, X_{tp}B_{K}(U_{t}))^{ \top}$, $\boldsymbol{{Z}} = (\boldsymbol{{Z}}_{1}, \cdots, \boldsymbol{{Z}}_{n})^{ \top}$, $\boldsymbol{{Y}} = (Y_{1}, \cdots, Y_{n})^{ \top}$. Let $\mathbb{E}(\boldsymbol{{\epsilon}}\boldsymbol{{\epsilon}}^{ \top}| \bf{X}_{t}, U_{t}) = \boldsymbol{{V}}$, and it is initially known. Subsequently, we can write the criterion as

Denoting the minimizer by $\hat{\boldsymbol{{b}}}$, we estimate $\beta_{j}(u)$ by $\hat{\beta}_{j}(u) = \sum\limits_{k = 1}^K \hat{b}_{jk}B_{k}(u)$.

3.   Asymptotic property

We have imposed the following conditions. (C1) and (C2) are mild regularity conditions. (C3) is necessary for the identification of the coefficient functions in varying coefficient models, as mentioned in Huang and Shen ^[20]. (C4), (C6), and (C7) are the same as those used in ^[1] for stationary mixing data. (C5) imposes some smoothness conditions on the coefficient functions.

3.1.   Assumptions

(C1) The eigenvalues of $\boldsymbol{{V}}$ are bounded away from zero and infinity.

(C2) The smoothing variable $U_{t}$ has a bounded density supported on $[a, b]$.

(C3) The eigenvalues of the matrix $\mathbb{E}(\bf{X}_{t} \bf{X}_{t}^{ \top}\mid U_{t} = u)$ are uniformly bounded away from zero and infinity for all $u\in [a, b]$.

(C4) The conditional density of $(U_{1}, U_{l+1})$ given $(\bf{X}_{1}, \bf{X}_{l+1})$ is uniformly bounded on the support of $(\bf{X}_{1}, \bf{X}_{l+1})$. The conditional density of $U_{1}$ given $\bf{X}_{1}$ is uniformly bounded on the support of $\bf{X}_{1}$.

(C5) For $g = \beta_{j}, 1\leq j \leq p,$ $g$ satisfies a Lipschitz condition of order $d > 1/2:|g^{(\lfloor d \rfloor)}(t)-g^{(\lfloor d \rfloor)}(s)|\leq C|s-t|^{d-\lfloor d \rfloor}$, where $\lfloor d \rfloor$ is the biggest integer strictly smaller than $d$ and $g^{(\lfloor d \rfloor)}$ is the $\lfloor d \rfloor-$ th derivative of $g$. The order of the B-spline used satisfies $q \geq d + 1/2$.

(C6) The process $\{Y_{t}, \bf{X}_{t}, U_{t}\}_{t\in\mathbb{Z}}$ is jointly strictly stationary with $U_{t}$ taking values in $\mathbb{R}$ and $\bf{X}_{t}$ taking values in $\mathbb{R}^{p}$. The $\alpha$-mixing coefficient $\alpha(l)$ of $\{Y_{t}, \bf{X}_{t}, Z_{t}\}_{t\in\mathbb{Z}}$ satisfies $\sum\limits_{l = 1}^\infty l^{c}\alpha(l)^{1-2/\delta} < \infty$ for some $\delta > 2$ and $c > 1-2/\delta$.

(C7) $\mathbb{E}(|X_{tj}|^{2\delta}) < \infty, j = 1, \cdots, d$, where $\delta$ is given in condition (C6).

3.2.   Theoretical result

Theorem 3.1: Assume (C1)$-$(C7) and that $K\rightarrow \infty, K^{3}/n\rightarrow 0$, then we have

As noted in Lai et al.^[12], we know that by choosing $K\asymp n^{1/(2d+1)}$, the well-known optimal convergence rate $O_{p}(n^{-2d/(2d+1)})$ is achieved.

If $\boldsymbol{{V}}$ is unknown and $\hat{\boldsymbol{{V}}}$ is an estimator of $\boldsymbol{{V}}$, then an application of the substitution principle leads to

A two-stage estimator $\bar{\boldsymbol{{b}}}$ can be obtained by minimizing (3.1); then, $\bar{\beta}_{j}(u) = \sum\limits_{k = 1}^K \bar{b}_{jk}B_{k}(u)$.

Theorem 3.2: Assume (C1)$-$(C7), $K\rightarrow \infty, K^{3}/n\rightarrow 0$ and that $\hat{\boldsymbol{{V}}}$ is consistent in probability estimating $\boldsymbol{{V}}$, then we have

Remark 1. In practice, $\boldsymbol{{V}}$ is usually unknown. It is customary to replace $\boldsymbol{{V}}$ in (2.1) with its consistent estimator $\hat{\boldsymbol{{V}}}$ ^[17,21]. Then, $\bar{\boldsymbol{{b}}}$ can be derived through a two-stage estimation. If the sample size is sufficiently large, we can split the data set into the training and test parts, and provide a consistent estimator of $\boldsymbol{{V}}$ from the residuals of ordinary least squares (OLS) for fitting (1.1) by grouping the training set.

3.3.   Computational aspects

However, the method in Remark 1 is not feasible when the sample size is small. It should be noted that Montoril et al.^[17] have presented an iterative procedure for estimating $\boldsymbol{{V}}$, but the iterative procedure is computationally expensive. Hence, they have provided a more efficient method. Generally, if the errors of the model follow an autoregressive process, we can estimate the $\beta_{j}(\cdot)$ by adopting some ideas of ^[17]. The estimation algorithm is as follows.

Step 1. Estimate the coefficient vector $\boldsymbol{{b}}$ by OLS, and denote it by $\boldsymbol{{b}}^{(0)}$.

Step 2. Compute the residuals via $\epsilon_{t}^{(0)} = Y_{t}-\boldsymbol{{Z}}_{t}\boldsymbol{{b}}^{(0)}$, and fit an autoregressive model to the residuals, i.e.,

Step 3. Letting $\boldsymbol{{b}}^{(0)}$ and $(\varphi_{1}^{(0)}\cdots \varphi_{p}^{(0)})$ in step 1 and step 2 as initial values, estimate $\boldsymbol{{b}}$ by minimizing numerically

where $\varphi_{p}(L) = 1-\varphi_{1}L-\cdots-\varphi_{p}L^{P}$, with the backshift operator satisfying $L^{k}V_{t} = V_{t-k}$, $k > 0$.

4.   Numerical examples

In this section, two simulated examples are considered: the threshold variable $U_{t}$ is an exogenous variable and the lagged variable of $Y_{t}$. For a given data set, we used equally spaced knots. The values used as candidates for the number of internal knots varied from one to five. The optimal number of internal knots is selected using the Bayesian information criterion (BIC) (Schwarz ^[22]). The BIC criterion function is defined as

where $n$ denotes the sample size, RMS denotes the residual mean square, $p$ is equal to the sum of number of autoregressive coefficients assumed for the errors and number of B-spline basis. The B-spline basis $\{B_{1}(U_{t}), \cdots, B_{K}(U_{t})\}$ can be obtained from function bs using the package splines in R language ^[3]. We consider time series data with lengths $n = 200,400$ and 600, and replicate the simulation 200 times in each case. For each replication, a total of $1000 + n$ observations were generated, and only the last $n$ observations were used to ensure approximate stationarity. The performance of estimators $\{\hat{\beta_{j}}(\cdot)\}$ can be demonstrated by the square root of average squared errors (RASE):

with

and $\{u_{k}, k = 1, \cdots, n_{\mathrm{grid}}\}$ are grid points on an interval over which the functions are evaluated. Because the range of the time series data varies from simulation to simulation, we need to select a common interval to compare the RASE values. The intervals selected for Example 4.1 and Example 4.2 are $[-0.45, 0.45]$ and $[-2, 2]$, respectively.

Example 4.1. Consider the following data generating process

with $\beta_{1}(u) = 0.9\sin(\pi u)$ and $\beta_{2}(u) = 0.85\cos(\pi u)$. We study two autoregressive errors, i.e. AR(1): $\epsilon_{t} = 0.8 \epsilon_{t-1}+e_{t}$ and AR(2): $\epsilon_{t} = 0.5 \epsilon_{t-1}+ 0.45 \epsilon_{t-2} + e_{t}$, where $e_{t}\stackrel{\mathrm{i.i.d}}{\sim} N(0, 0.2^{2})$ and $U_{t}\stackrel{\mathrm{i.i.d}}{\sim} U[-0.5, 0.5]$.

Example 4.2. We now discuss an exponential autoregressive (EXPAR) model

with $\beta_{1}(u) = 0.6+0.9e^{(-u^2)}$ and $\beta_{2}(u) = -0.4-1.2e^{(-u^2)}$. In this case, the autoregressive errors we study are AR(1): $\epsilon_{t} = 0.9\epsilon_{t-1}+e_{t}$ and AR(2): $\epsilon_{t} = 0.6 \epsilon_{t-1}+ 0.35 \epsilon_{t-2} + e_{t}$, where $e_{t} \stackrel{\mathrm{i.i.d}}{\sim}N(0, 0.2^{2})$.

5.   Numerical results and discussion

In this section, the numerical results for Examples 4.1 and 4.2 are presented. We compare the performance of local linear (Local) estimators proposed by Cai et al.^[1], the spline estimators under the assumption of independent errors (Spl.ind) and spline estimators (Spl.cor) proposed by us. Tables 1 and 2 show the mean and standard deviation (in parentheses) of the RASEs for $\hat{\beta}_{1}(\cdot)$ and $\hat{\beta}_{2}(\cdot)$ with linear ($k = 1$) and cubic ($k = 3$) splines under different AR($p$) errors. It is apparent that the standard deviation of the RASEs in columns Spl.cor decrease with the increase in the sample size $n$. In addition, the results of linear splines are similar to those of the cubic splines. Moreover, it can be seen that the proposed approach performs better than the method that ignores the correlated errors. Based on cubic splines, Tables 3 and 4 provide the resulting estimates while considering different AR (1) errors: $\epsilon_{t} = \theta \epsilon_{t-1}+e_{t}$ with $\theta = 0.3, 0.6, 0.9, 0.95$ when $n = 600$. We can see that with the increase in the correlation levels $\theta$, the performance of Local and Spl.ind worsens, but Spl.cor remains relatively stable. Furthermore, Spl.cor is always better than Spl.ind and Local, which emphasizes the importance of considering the correlated errors. The computational times of Example 4.2 are reported in Table 5. The local linear method is computationally expensive when the sample size is large. Moreover, when considering the correlated errors, the algorithm needs to search for the optimal solution iteratively. Therefore, the calculation time of Spl.cor is longer than that of Spl.ind.

The estimated $\beta_{1}(\cdot)$ and $\beta_{2}(\cdot)$ from a typical sample using cubic splines for Examples 4.1 and 4.2 with AR(1) error structures are plotted in Figures 1 to 4. The solid curve is the true curve, and the dotted curve represents the typical estimated curve. The typical sample is chosen in such a way that its RASE value is equal to the median in the 200 simulations ^[23]. Clearly, the proposed estimators capture the main features of the true functions well. Although it is not shown here, the examples with AR(2) error structures show similar results.

6.   Conclusions

This study considered the B-spline estimation for varying coefficient models with correlated errors using a weighted least squares method. In terms of the theoretical results, convergence rates were derived under the $\alpha$-mixing condition. Simulation studies were conducted to illustrate the performance of the proposed estimation method which showed better performance while considering correlated errors in fitting varying coefficient models.

Acknowledgments

This work was supported in part by the National Natural Science Foundation of China under Grant No. 61973096 and by GDUPS (2019).

Conflict of interest

The authors declare that they have no any competing interests.

Appendix

A.   Proof of the main results

Notation: For two positive sequences ${a_{n}}$ and $b_{n}$, $a_{n}\lesssim b_{n}$ means that $a_{n}/ b_{n}$ is uniformly bounded, $a_{n}\asymp b_{n}$ if $a_{n}\lesssim b_{n}$ and $b_{n}\lesssim a_{n}$. Let $\|g\|_{2} = \{\int_{[a, b]}g^{2}(x)dx\}^{1/2}$ be the $L_{2}$-norm of a square integrable function $g(\cdot)$ on $[a, b]$. $\mathbb{R}$ denotes the set of real numbers, and $\mathbb{Z}$ denotes the set of integers. $|\cdot|$ denotes the Euclidean norm of a vector. $\boldsymbol{{I}}$ is a identity matrix. $\boldsymbol{{e}}_{j}$ denotes a unit vector whose $j$-th entry is 1 and all other entries of which are $0$. We denote by $\lambda_{\min}(\cdot)$ $(\lambda_{\max}(\cdot))$ the smallest (the largest) eigenvalue of a matrix. Let $C$ denote a generic constant that might assume different values at different places. Given random variables $V_{n}, n\geq 1$, let $V_{n} = O_{p}(b_{n})$ means that the random variables $V_{n}/b_{n}, n\geq 1$ are bounded in probability, that is,

And $V_{n} = o_{p}(b_{n})$ means that the random variables $V_{n}/b_{n}, n\geq 1$ converge to zero in probability, namely

In order to proof the theorems, we need the following Lemmas.

Lemma 1. (Lemma 2, Lai et al.^[12]) Under assumptions (C2) and (C3), there are constants $0 < b_{1} < b_{2} < \infty$ such that the eigenvalues of $\mathbb{E} \boldsymbol{{Z}}_{1} \boldsymbol{{Z}}_{1}^{ \top}$ fall in $[\frac{b_{1}}{K}, \frac{b_{2}}{K}]$, and under the additional conditions (C4)$-$(C7), there are constants $0 < b_{3} < b_{4} < \infty$ such that the eigenvalues of $\sum\limits_{t = 1}^n \boldsymbol{{Z}}_{t} \boldsymbol{{Z}}_{t}^{ \top}/n$ fall in $[\frac{b_{3}}{K}, \frac{b_{4}}{K}]$ with probability approaching 1 as $n\rightarrow \infty$.

Lemma 2. Under the conditions (C1)$-$(C7), we have $\boldsymbol{{\epsilon}}^{ \top}\boldsymbol{{V}}^{-1} \boldsymbol{{Z}}(\boldsymbol{{Z}}^{ \top}\boldsymbol{{V}}^{-1} \boldsymbol{{Z}})^{-1} \boldsymbol{{Z}}^{ \top}\boldsymbol{{V}}^{-1}\boldsymbol{{\epsilon}} = O_{p}(K)$.

Proof. Note that $\mathbb{E}(\boldsymbol{{Z}}^{ \top} \boldsymbol{{Z}}) = \mathbb{E}(\sum\limits_{t = 1}^n \boldsymbol{{Z}}_{t} \boldsymbol{{Z}}_{t}^{ \top}) = n \mathbb{E}(\boldsymbol{{Z}}_{1} \boldsymbol{{Z}}_{1}^{ \top}).$ Hence, by Lemma 1 we obtain

Using the Markov inequality, we can derive

On the other hand, by Lemma 1 and (C1),

Consequently,

Then, it follows from (A.2) and (A.3) that

This completes the proof of Lemma 2.

Lemma 3. Suppose $\beta_{0j}(u) = \sum\limits_{k = 1}^K b_{0jk}B_{k}(u)$ is the best approximating B-spline for $\beta_{j}(u)$ with $\|\beta_{0j}-\beta_{j}\|_{\infty} = O(K^{-d})$ (such approximation property is well-known for B-spline under smoothness assumptions (C5), see for example ^[19]). Let $\boldsymbol{{b}}_{0} = (b_{011}, \cdots, b_{01K}, \cdots, b_{0p1}, \cdots, b_{0pK})^{ \top}$, $\boldsymbol{{\eta}} = \boldsymbol{{P}}_{\boldsymbol{{V}}^{-\frac{1}{2}} \boldsymbol{{Z}}}\boldsymbol{{V}}^{-\frac{1}{2}}(\boldsymbol{{Y}}- \boldsymbol{{Z}}\boldsymbol{{b}}_{0})$, where $\boldsymbol{{P}}_{\boldsymbol{{V}}^{-\frac{1}{2}} \boldsymbol{{Z}}} = \boldsymbol{{V}}^{-\frac{1}{2}} \boldsymbol{{Z}}(\boldsymbol{{Z}}^{ \top}\boldsymbol{{V}}^{-1} \boldsymbol{{Z}})^{-1} \boldsymbol{{Z}}^{ \top}\boldsymbol{{V}}^{-\frac{1}{2}}$ is a projection matrix. Then under (C1)$-$(C7) we have

Proof. Denote $r_{t} = \sum\limits_{j = 1}^p X_{tj}\beta_{j}(U_{t})$ and $\boldsymbol{{r}} = (r_{1}, \ldots, r_{n})^{ \top},$ thus $\boldsymbol{{Y}} = \boldsymbol{{r}} + \boldsymbol{{\epsilon}}$. We have

Noting that $\boldsymbol{{P}}_{\boldsymbol{{V}}^{-\frac{1}{2}} \boldsymbol{{Z}}}$ is an idempotent matrix, by Lemma 2 we obtain that

By the approximation property of B-spline we have $\|\boldsymbol{{r}} - \boldsymbol{{Z}}\boldsymbol{{b}}_{0}\|^{2} = O_{p}(\frac{n}{K^{2d}})$, hence

Thus, $S_{2} = O_{p}(\frac{n}{K^{2d}})$, the proof is complete.

Proof of Theorem 3.1.

Proof. By the definition of $\hat{\boldsymbol{{b}}}$, we get

Applying the Cauchy-Schwartz inequality, the equation (A.4) can be continued as

Hence, by Lemma 3 we obtain

Then, it follows from Lemma 1 and (C1) that

which means that $\|\boldsymbol{{b_{0}}}-\hat{\boldsymbol{{b}}}\|^{2} = O_{p}(\frac{K^{2}}{n}+\frac{1}{K^{2d-1}})$. Then, by the property of B-spline (De Boor^[19])

for some constants $b_{5}, b_{6} > 0$, we can derive

Proof of Theorem 3.2.

Proof. The proof of this theorem is similar to that of Theorem 2.2 in ^[17].

Let $\bf{R} = \boldsymbol{{Z}}^{ \top}\boldsymbol{{V}}^{-1} \boldsymbol{{Z}}$ and $\hat{ \bf{R}} = \boldsymbol{{Z}}^{ \top}\hat{\boldsymbol{{V}}}^{-1} \boldsymbol{{Z}}$. Note that $\hat{\boldsymbol{{V}}}$ is consistent, then by the Theorem 1.4.2 in ^[24] we have

This implies that $\hat{ \bf{R}}$ is a consistent estimator of $\bf{R}$, namely $\hat{ \bf{R}}^{-1} \bf{R}-\boldsymbol{{I}} = o_{p}(1)$. It follows from (2.1) that

Therefore, by (3.1) we can derive that

This ensures that

By (A.5), (A.6) and Theorem 3.1, we conclude that

This completes the proof of Theorem 3.2.