PQMLE of a Partially Linear Varying Coefficient Spatial Autoregressive Panel Model with Random Effects

Abstract

This article deals with asymmetrical spatial data which can be modeled by a partially linear varying coefficient spatial autoregressive panel model (PLVCSARPM) with random effects. We constructed its profile quasi-maximum likelihood estimators (PQMLE). The consistency and asymptotic normality of the estimators were proved under some regular conditions. Monte Carlo simulations implied our estimators have good finite sample performance. Finally, a set of asymmetric real data applications was analyzed for illustrating the performance of the provided method.

1. Introduction

Spatial regression models are used to deal with spatially dependent data which widely exist in many fields, such as economics, environmental science and geography. According to the different types of spatial interaction effects, spatial regression models can be sorted into three basic categories (see [1]). The first category is spatial autoregressive (SAR) models, which include endogenous interaction effects among observations of response variables at different spatial scales (see [2]). The second category is spatial durbin models (SDM), which include exogenous interaction effects among observations of covariates and endogenous interaction effects among observations of response variables at different spatial scales (see [3,4]). The third category is spatial error models (SEM), which include interaction effects among disturbance terms of different spatial scales (see [5]). Among them, the SAR models proposed by [6] may be the most popular. The developments in the testing and estimation of SAR models for sectional data were summarized in the books by [7,8,9] and the surveys by [10,11,12,13,14,15,16], among others. Compared to sectional data models, panel data models exhibit an excellent capacity for capturing the complex situations by using abundant datasets built up over time and adding individual-specific or time-specific fixed or random effects. Their theories, methods and applications can be found in the books by [17,18,19,20] and the surveys by [21,22,23,24,25,26], among others.

The above mentioned research literature mainly focuses on linear parametric models. Although the estimations and properties of these models have been well established, they are often unrealistic in application, for the reason that they are unable to accommodate sufficient flexibility to accommodate complex structures (e.g., nonlinearity). Moreover, mis-specification of the data generation mechanism by a linear parametric model could lead to excessive modeling biases or even erroneous conclusions.

Ref. [27] pointed out that the relationship between variables in space usually exhibits highly complexity in reality. Therefore, researches have proposed a number of solution methods. In video coding systems, transmission problems are usually dealt with wavelet-based methods. Ref. [28] proposed a low-band-shift method which avoids the shift-variant property of the wavelet transform and performs the motion compensation more precisely and efficiently. Ref. [29] obtained a wavelet-based lossless video coding scheme that switches between two operations based on the different amount of motion activities between two consecutive frames. More theories and applications can be found in [30,31,32]. In econometrics, some nonparametric and semiparametric spatial regression models have been developed to relax the linear parametric model settings. For sectional spatial data, ref. [33] studied the GMM estimation of a nonparametric SAR model; Ref. [34] investigated the PQMLE of a partially linear nonparametric SAR model. However, such nonparametric SAR models may cause the “curse of dimensionality” when the dimension of covariates is higher. In order to overcome this drawback, ref. [35] studied two-step SGMM estimators and their asymptotic properties for spatial models with space-varying coefficients. Ref. [36] proposed a semiparametric series-based least squares estimating procedure for a semiparametric varying coefficient mixed regressive SAR model and derived the asymptotical normality of the estimators. Some other related research works can be found in [37,38,39,40,41,42]. For panel spatial data, ref. [43] obtained PMLE and its asymptotical normality of varying the coefficient SAR panel model with random effects; Ref. [44] applied instrumental variable estimation to a semiparametric varying coefficient spatial panel data model with random effects and the investigated asymptotical normality of the estimators.

In this paper, we extend the varying coefficient spatial panel model with random effects given in [43] to a PLVCSARPM with random effects. By adding a linear part to the model of [43], we can simultaneously capture linearity, non-linearity and the spatial correlation relationships of exogenous variables in a response variable. By using the profile quasimaximum likelihood method to estimate PLVCSARPM with random effects, we proved the consistency and asymptotic normality of the estimators under some regular conditions. Monte Carlo simulations and real data analysis show that our estimators perform well.

This paper is organized as follows: Section 2 introduces the PLVCSARPM with random effects and constructs its PQMLE. Section 3 proves the asymptotic properties of estimators. Section 4 presents the small sample estimates using Monte Carlo simulations. Section 5 analyzes a set of asymmetric real data applications for illustrating the performance of the proposed method. A summary is given in Section 6. The proofs of some important theorems and lemmas are given in Appendix A.

2. The Model and Estimators

Consider the following partially linear varying coefficient spatial autoregressive panel model (PLVCSARPM) with random effects:

$$y_{it}=\rho {\left(W_0{Y}_t\right)}_i+x_{it}^{{}^{\prime }}\beta \left(u_{it}\right)+z_{it}^{{}^{\prime }}\alpha +b_i+{\epsilon }_{it},1\le i\le N,1\le t\le T,$$

(1)

where i refers to a spatial unit; t refers to a given time period; $y_{it}$ are observations of a response variable; $Y_t={(y_{1t},y_{2t},\cdots ,y_{Nt})}^{{}^{\prime }}$; $x_{it}={(x_{it1},x_{it2},\cdots ,x_{itp})}^{{}^{\prime }}$ and $z_{it}={(z_{it1},z_{it2},\cdots ,z_{itq})}^{{}^{\prime }}$ are observations of p-dimensional and q-dimensional covariates, respectively; $\beta \left(u_{it}\right)={({\beta }_1\left(u_{it}\right),{\beta }_2\left(u_{it}\right),\cdots ,{\beta }_p\left(u_{it}\right))}^{\prime }$ is an unknown univariate varying coefficient function vector; ${\beta }_s\left(u\right)$$(s=1,2,\cdots ,p)$ are unknown smoothing functions of u; $\rho$$\left(\right|\rho |<1)$ is an unknown spatial correlation coefficient; $\alpha ={({\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_q)}^{\prime }$ is a regression coefficient vector of $z_{it}$; $W_0$ is an $N\times N$ predetermined spatial weight matrix; ${\left(W_0{Y}_t\right)}_i$ is the ith component of $W_0{Y}_t$; ${\epsilon }_{it}$ are i.i.d. error terms with zero means and variance ${\sigma }_{\epsilon }^2$; $b_i$ are i.i.d. variables with zero means and variance ${\sigma }_b^2$, ${\epsilon }_{it}$ are independent of $b_i$. Let ${\theta }_0={({\rho }_0,{\alpha }_0^{\prime },{\sigma }_{b0}^2,{\sigma }_{\epsilon 0}^2)}^{\prime }$ be the true parameter vector of $\theta ={(\rho ,{\alpha }^{\prime },{\sigma }_b^2,{\sigma }_{\epsilon }^2)}^{\prime }$ and ${\beta }_0\left(u_{it}\right)$ be the true varying coefficient function of $\beta \left(u_{it}\right)$.

The model (1) can be simplified as the following matrix form:

$$Y=\rho WY+X\beta \left(u\right)+Z\alpha +Ub+\epsilon ,$$

(2)

where $Y={(y_{11},y_{12},\cdots ,y_{NT})}^{{}^{\prime }}$; $X={(x_{11}^{{}^{\prime }},x_{12}^{{}^{\prime }},\cdots ,x_{NT}^{{}^{\prime }})}^{\prime }$; $Z={(z_{11}^{{}^{\prime }},z_{12}^{{}^{\prime }},\cdots ,z_{NT}^{{}^{\prime }})}^{\prime }$; $\epsilon ={({\epsilon }_{11},{\epsilon }_{12},\cdots ,{\epsilon }_{NT})}^{{}^{\prime }}$; $W=W_0\otimes I_T$; $b={(b_1,b_2,\cdots ,b_N)}^{{}^{\prime }}$; $U=I_N\otimes e_T$; $I_N$ is an $N\times N$ identity matrix; $I_T$ is a $T\times T$ identity matrix; ⊗ denotes the Kronecker product; $e_T$ is a $T\times 1$ vector consisting of 1.

Define $A\left(\rho \right)=I-\rho W$; then the model (2) can be rewritten as:

$$A\left(\rho \right)Y=X\beta \left(u\right)+Z\alpha +Ub+\epsilon ,$$

(3)

where I is an $NT\times NT$ unit matrix. For the model (3), it is easy to get the following facts:

$$\partial (Ub+\epsilon )/\partial \left(Y\right)=A\left(\rho \right),\mathsf{\Sigma }=Var(Ub+\epsilon )={\sigma }_{\epsilon }^{2}I+{\sigma }_b^2{I}_N\otimes \left(e_T{e}_T^{{}^{\prime }}\right),$$

$$\left|\mathsf{\Sigma }\right|={\sigma }_{\epsilon }^{2N(T-1)}{({\sigma }_b^2+T{\sigma }_{\epsilon }^2)}^N,{\mathsf{\Sigma }}^{-1}=\frac{1}{{\sigma }_{\epsilon }^2}I+(\frac{1}{{\sigma }_{\epsilon }^2+T{\sigma }_b^2}-\frac{1}{T{\sigma }_b^2})I_N\otimes \left(\frac{1}{T}{e}_T{e}_T^{{}^{\prime }}\right).$$

According to [12], the quasi-log-likelihood function of the model (3) can be written as follows:

$$\begin{array}{cc}\hfill lnL\left(\theta \right)=& -\frac{N(T-1)}{2}ln{\sigma }_{\epsilon }^2-\frac{N}{2}ln({\sigma }_{\epsilon }^2+T{\sigma }_b^2)\hfill \\ \hfill & -\frac{1}{2({\sigma }_{\epsilon }^2+T{\sigma }_b^2)}{[A\left(\rho \right)Y-X\beta \left(u\right)-Z\alpha ]}^{\prime }H[A\left(\rho \right)Y-X\beta \left(u\right)-Z\alpha ]\hfill \\ \hfill & -\frac{1}{2{\sigma }_{\epsilon }^2}{[A\left(\rho \right)Y-X\beta \left(u\right)-Z\alpha ]}^{\prime }(I-H)[A\left(\rho \right)Y-X\beta \left(u\right)-Z\alpha ]+ln\left|A\left(\rho \right)\right|+c,\hfill \end{array}$$

(4)

where $H=I_N\otimes \left(T^{-1}{e}_T{e}_T^{\prime }\right)$, and c is a constant.

By maximizing the above quasi-log-likelihood function with respect to $\theta$, the quasi-maximum likelihood estimators of $\alpha$, ${\sigma }_b^2$ and ${\sigma }_{\epsilon }^2$ can be easily obtained as

$$\begin{array}{cc}\hfill & \tilde{\alpha }={\left[Z^{\prime }(I-S)(\frac{H}{{\sigma }_{\epsilon }^2+T{\sigma }_b^2}+\frac{I-H}{{\sigma }_{\epsilon }^2})(I-S)Z\right]}^{-1}{Z}^{\prime }(I-S)(\frac{H}{{\sigma }_{\epsilon }^2+T{\sigma }_b^2}+\frac{I-H}{{\sigma }_{\epsilon }^2})A\left(\rho \right)Y,\hfill \end{array}$$

(5)

$${\tilde{\sigma }}_{\epsilon }^2={\left[N(T-1)\right]}^{-1}{[A\left(\rho \right)Y-X\beta \left(u\right)-Z\alpha ]}^{\prime }(I-H)[A\left(\rho \right)Y-X\beta \left(u\right)-Z\alpha ],$$

(6)

$$\begin{array}{cc}\hfill & {\tilde{\sigma }}_b^2={\left(NT\right)}^{-1}{[A\left(\rho \right)Y-X\beta \left(u\right)-Z\alpha ]}^{\prime }H[A\left(\rho \right)Y-X\beta \left(u\right)-Z\alpha ]-T^{-1}{\tilde{\sigma }}_{\epsilon }^2.\hfill \end{array}$$

(7)

By substituting (5)–(7) into (4), we have the concentrated quasi-log-likelihood function of $\rho$ as

$$\begin{array}{cc}\hfill & lnL\left(\rho \right)=\frac{-N(T-1)}{2}ln{\tilde{\sigma }}_{\epsilon }^2-\frac{N}{2}ln({\tilde{\sigma }}_{\epsilon }^2+T{\tilde{\sigma }}_b^2)-\frac{NT}{2}+ln\left|A\left(\rho \right)\right|.\hfill \end{array}$$

It is obvious that we cannot directly obtain the quasi-maximum likelihood estimator of $\rho$ by maximizing the above formula because $\beta \left(u_{it}\right)$ is an unknown function. In order to overcome this problem, we use the PQMLE method and working independence theory ([45,46]) to estimate the unknown parameters and varying coefficient functions of the model (1).

The main steps are as follows:

Step 1 Suppose that $\theta$ is known. ${\beta }_s(·)$ can be approximated by the first-order Taylor expansion

$${\beta }_s\left(u_{it}\right)\approx {\beta }_s\left(u\right)+(u_{it}-u){\dot{\beta }}_s\left(u\right)\widehat{=}{a}_{1s}+a_{2s}(u_{it}-u),s=1,2,\cdots ,p$$

for $u_{it}$ in a neighborhood of u, where $a_{2s}$ is the first order derivative of $\beta \left(u\right)$. Let $a_1={(a_{11},a_{12},\cdots ,a_{1p})}^{\prime }$ and $a_2={(a_{21},a_{22},\cdots ,a_{2p})}^{\prime }$; then estimators of $a_1$ and $a_2$ can be obtained by

$$({\widehat{a}}_1,{\widehat{a}}_2)=arg\underset{a_1,a_2}{min}\frac{1}{NT}\sum _{i=1}^N\sum _{t=1}^T{[{\tilde{y}}_{it}-\sum _{s=1}^p(a_{1s}+a_{2s}(u_{it}-u))x_{its}]}^2{K}_h(u_{it}-u),$$

where ${\tilde{y}}_{it}=y_{it}-\rho {\left(W_0{Y}_t\right)}_i-z_{it}^{\prime }\alpha$, $K_h(u_{it}-u)=h^{-1}K((u_{it}-u)/h)$, $K(·)$ is a kernel function and h is the bandwidth. Therefore, the feasible initial estimator of $\beta \left(u\right)$ can be obtained by ${\widehat{\beta }}_{IN}\left(u\right)={\widehat{a}}_1={({\widehat{a}}_{11},{\widehat{a}}_{12},\cdots ,{\widehat{a}}_{1p})}^{\prime }$.

Denote $W_u=diag(K_h(u_{11}-u),K_h(u_{12}-u),\cdots ,K_h(u_{NT}-u))$, $\tilde{Y} = ({\tilde{y}}_{11},{\tilde{y}}_{12},\cdots ,$${\tilde{y}}_{NT}{)}^{\prime }$, $\delta ={(a_{11},\cdots ,a_{1p},h{a}_{21},\cdots ,h{a}_{2p})}^{\prime }$ and $D_u={\left(\begin{array}{ccc}{x}_{11}& \cdots & x_{NT}\\ \frac{u_{11}-u}{h}{x}_{11}& \cdots & \frac{u_{NT}-u}{h}{x}_{NT}\end{array}\right)}^{\prime },$ then we have

$$\widehat{\delta }=arg\underset{\delta }{min}\frac{1}{NT}{(\tilde{Y}-D_u\delta )}^{\prime }{W}_u(\tilde{Y}-D_u\delta ).$$

Therefore, we obtain

$$\widehat{\delta }={\left(D_u^{\prime }{W}_u{D}_u\right)}^{-1}{D}_u^{\prime }{W}_u\tilde{Y}.$$

Let $S_u=\frac{1}{NT}{D}_u^{\prime }{W}_u{D}_u$, $T_u=\frac{1}{NT}{D}_u^{\prime }{W}_u\tilde{Y}$ and $e_0=(I_p,0_p)$. It is easy to know that

$${\widehat{\beta }}_{IN}\left(u\right)={({\widehat{a}}_{11},{\widehat{a}}_{12},\cdots ,{\widehat{a}}_{1p})}^{\prime }=e_0^{{}^{\prime }}{S}_u^{-1}{T}_u,$$

(8)

where $e_0^{{}^{\prime }}{S}_u=s\left(u\right)$. Consequently, the initial estimator of $X\beta \left(u\right)$ is given by

$$X{\widehat{\beta }}_{IN}\left(u\right)=\left(\begin{array}{c}\left(x_{11}^{\prime }0\right){\left(D_{u_{11}}^{\prime }{W}_{u_{11}}{D}_{u_{11}}\right)}^{-1}{D}_{u_{11}}^{\prime }{W}_{u_{11}}\\ \left(x_{12}^{\prime }0\right){\left(D_{u_{12}}^{\prime }{W}_{u_{12}}{D}_{u_{12}}\right)}^{-1}{D}_{u_{12}}^{\prime }{W}_{u_{12}}\\ ⋮\\ \left(x_{NT}^{\prime }0\right){\left(D_{u_{NT}}^{\prime }{W}_{u_{NT}}{D}_{u_{NT}}\right)}^{-1}{D}_{u_{NT}}^{\prime }{W}_{u_{NT}}\end{array}\right)\tilde{Y}\widehat{=}S\tilde{Y},$$

where

$$S=\left(\begin{array}{c}\left(x_{11}^{\prime }0\right){\left(D_{u_{11}}^{\prime }{W}_{u_{11}}{D}_{u_{11}}\right)}^{-1}{D}_{u_{11}}^{\prime }{W}_{u_{11}}\\ \left(x_{12}^{\prime }0\right){\left(D_{u_{12}}^{\prime }{W}_{u_{12}}{D}_{u_{12}}\right)}^{-1}{D}_{u_{12}}^{\prime }{W}_{u_{12}}\\ ⋮\\ \left(x_{NT}^{\prime }0\right){\left(D_{u_{NT}}^{\prime }{W}_{u_{NT}}{D}_{u_{NT}}\right)}^{-1}{D}_{u_{NT}}^{\prime }{W}_{u_{NT}}\end{array}\right).$$

Step 2 Replacing $\beta \left(u\right)$ of (4) with ${\widehat{\beta }}_{IN}\left(u\right)$, the estimator $\widehat{\theta }$ of $\theta$ can be obtained by maximizing approximate quasi-log-likelihood function:

$$\begin{array}{cc}\hfill ln\widehat{L}\left(\theta \right)=& -\frac{N(T-1)}{2}ln{\sigma }_{\epsilon }^2-\frac{1}{2({\sigma }_{\epsilon }^2+T{\sigma }_b^2)}{[A\left(\rho \right)Y-Z\alpha ]}^{\prime }{(I-S)}^{\prime }H(I-S)[A\left(\rho \right)Y-Z\alpha ]\hfill \\ \hfill & -\frac{N}{2}ln({\sigma }_{\epsilon }^2+T{\sigma }_b^2)-\frac{1}{2{\sigma }_{\epsilon }^2}{[A\left(\rho \right)Y-Z\alpha ]}^{\prime }{(I-S)}^{\prime }(I-H)(I-S)[A\left(\rho \right)Y-Z\alpha ]\hfill \\ \hfill & +ln\left|A\right(\rho \left)\right|+c.\hfill \end{array}$$

(9)

In the real estimation of $\theta$, the procedure is realized by following steps:

Firstly, assume $\rho$ is known. The initial estimators of ${\sigma }_b^2$, ${\sigma }_{\epsilon }^2$ and $\alpha$ are obtained by maximizing (9). Then, we find

$${\widehat{\sigma }}_{\epsilon IN}^2={\left[N(T-1)\right]}^{-1}{[\tilde{Y}-X{\widehat{\beta }}_{IN}\left(u\right)]}^{\prime }(I-H)[\tilde{Y}-X{\widehat{\beta }}_{IN}\left(u\right)],$$

(10)

$${\widehat{\sigma }}_{bIN}^2={\left(NT\right)}^{-1}{[\tilde{Y}-X{\widehat{\beta }}_{IN}\left(u\right)]}^{\prime }H[\tilde{Y}-X{\widehat{\beta }}_{IN}\left(u\right)]-T^{-1}{\widehat{\sigma }}_{\epsilon IN}^2,$$

(11)

$${\widehat{\alpha }}_{IN}={\left[Z^{\prime }(I-S)(\frac{H}{{\widehat{\sigma }}_{\epsilon IN}^2+T{\widehat{\sigma }}_{bIN}^2}+\frac{I-H}{{\widehat{\sigma }}_{\epsilon IN}^2})(I-S)Z\right]}^{-1}{Z}^{\prime }(I-S)(\frac{H}{{\widehat{\sigma }}_{\epsilon IN}^2+T{\widehat{\sigma }}_{bIN}^2}+\frac{I-H}{{\widehat{\sigma }}_{\epsilon IN}^2})A\left(\rho \right)Y.$$

(12)

Secondly, with the estimated ${\widehat{\sigma }}_{\epsilon IN}^2$, ${\widehat{\sigma }}_{bIN}^2$ and ${\widehat{\alpha }}_{IN}$, update $\widehat{\rho }$ by maximizing the concentrated quasi-log-likelihood function of $\rho$:

$$ln\tilde{L}\left(\rho \right)=-\frac{N(T-1)}{2}ln{\widehat{\sigma }}_{\epsilon IN}^2-\frac{N}{2}ln({\widehat{\sigma }}_{\epsilon IN}^2+T{\widehat{\sigma }}_{bIN}^2)+ln\left|A\left(\rho \right)\right|+c.$$

Therefore, the estimator of $\rho$ is obtained by:

$$\widehat{\rho }=arg\underset{\rho }{max}ln\tilde{L}\left(\rho \right).$$

Step 3 By substituting $\widehat{\rho }$ into (10), (11) and (12), respectively, the final estimators of ${\sigma }_b^2$, ${\sigma }_{\epsilon }^2$ and $\alpha$ are computed as follows:

$${\widehat{\sigma }}_{\epsilon }^2={\left[N(T-1)\right]}^{-1}{\widehat{\tilde{Y}}}^{\prime }{(I-S)}^{\prime }(I-H)(I-S)\widehat{\tilde{Y}},$$

(13)

$${\widehat{\sigma }}_b^2={\left(NT\right)}^{-1}{\widehat{\tilde{Y}}}^{\prime }(I-S)H(I-S)\widehat{\tilde{Y}}-T^{-1}{\widehat{\sigma }}_{\epsilon }^2,$$

$$\widehat{\alpha }={\left[Z^{\prime }(I-S)(\frac{H}{{\widehat{\sigma }}_{\epsilon }^2+T{\widehat{\sigma }}_b^2}+\frac{I-H}{{\widehat{\sigma }}_{\epsilon }^2})(I-S)Z\right]}^{-1}{Z}^{\prime }(I-S)(\frac{H}{{\widehat{\sigma }}_{\epsilon }^2+T{\widehat{\sigma }}_b^2}+\frac{I-H}{{\widehat{\sigma }}_{\epsilon }^2})A\left(\widehat{\rho }\right)Y.$$

Step 4 By replacing $\rho$ with $\widehat{\rho }$ in (8), we get the ultimate estimator of $\beta \left(u\right)$:

$$\widehat{\beta }\left(u\right)=e_0^{{}^{\prime }}{\left[D_u^{\prime }{W}_u{D}_u\right]}^{-1}{D}_u^{\prime }{W}_u\widehat{\tilde{Y}}.$$

where $\widehat{\tilde{Y}}=A\left(\widehat{\rho }\right)Y-Z\widehat{\alpha }$, $A\left(\widehat{\rho }\right)=I-\widehat{\rho }W$.

3. Asymptotic Properties for the Estimators

In this section, we focus on studying consistency and asymptotic normality of the PQMLEs given in Section 2. To prove these asymptotic properties, we need the following assumptions to hold.

To provide a rigorous analysis, we make the following assumptions.

Assumption 1.

(i)

${\epsilon }_{it}$ and $b_i$ are uncorrelated to $x_{it}$ and $z_{it}$, and satisfy $E\left(b_i\right|x_{it},z_{it})=0$, $Var\left(b_i\right|x_{it},z_{it})={\sigma }_b^2<\infty$, $E(\parallel b_i{x}_{it}^{\prime }\parallel )<\infty$, $E(\parallel b_i{z}_{it}^{\prime }\parallel )<\infty$, $E\left({\epsilon }_i\right|x_{it},z_{it})=0$, $E(\parallel {\epsilon }_i{x}_{it}^{\prime }\parallel )<\infty$, $Var\left({\epsilon }_i\right|x_{it},z_{it})={\sigma }_{\epsilon }^2<\infty$ and $E(\parallel {\epsilon }_i{z}_{it}^{\prime }\parallel )<\infty$, where $\parallel ·\parallel$ represents the Euclidean norm. Moreover, $E\left(\right|\epsilon {|}^r)<\infty$ and $E\left(\right|b{|}^r)<\infty$ for some $r>4$.

(ii)

$\left\{\right(u_{it}$, ${\epsilon }_{it}{\left)\right\}}_{i=1,t=1}^{N,T}$ are $i.i.d$ random sequences. The density function $f_t\left(u\right)$ of $u_{it}$ is non-zero, uniformly bounded and second-order continuously differentiable on $\mathbb{U}$, where $\mathbb{U}$ is the supporting set of $K\left(u\right)$. ${\mathsf{\Omega }}_{11}\left(u\right)=E\left(x_{it}{x}_{it}^{\prime }\right|u_{it}=u)<\infty$ exists, and it has a second-order continuous derivative; ${\mathsf{\Omega }}_{12}\left(u\right)=E\left(x_{it}{z}_{it}^{\prime }\right|u_{it}=u)$ and ${\mathsf{\Omega }}_{22}\left(u\right)=E\left(z_{it}{z}_{it}^{\prime }\right|u_{it}=u)$ are second-order continuously differentiable on $\mathbb{U}$; $E\left(x_{it}{x}_{is}^{\prime }\right|u_{it}=u)=0$, $E\left(z_{it}{z}_{is}^{\prime }\right|u_{it}=u)=0$ and $E\left(x_{it}{z}_{is}^{\prime }\right|u_{it}=u)=0$ for $t\ne s$; $E\left(x_{it}{x}_{it}^{\prime }\right|u_{it}=u)\ne 0$, $E\left(x_{it}{\epsilon }_{it}\right|u_{it})=0$ and $E\left(x_{it}{b}_{it}\right|u_{it})=0$.

(iii)

Real valued function ${\beta }_s\left(u\right)(s=1,2,\cdots ,p)$ is second-order, continuously differentiable and satisfies the first order Lipschitz condition, $|x_{it}^{\prime }{\beta }_s\left(u\right)|\le m_{\beta }$ at any u, $|z_{it}^{\prime }\alpha |\le m_{\alpha }$, where $m_{\beta }$ and $m_{\alpha }$ are positive and constant.

(iv)

${\{x_{it},z_{it}\}}_{i=1,t=1}^{N,T}$ are i.i.d. random variables from the population. Moreover, $\underset{i,t}{max}E\left(\right|x_{it}{|}^r)<\infty$ and $\underset{i,t}{max}E\left(\right|z_{it}{|}^r)<\infty$ for $r\ge 4$.

Assumption 2.

(i)

As a normalization, the diagonal elements ${\varpi }_{ii}$ of $W_0$ are 0 for all i and ${\varpi }_{ij}$ are at most of order $l_N^{-1}$, denoted by $O(1/l_N)$.

(ii)

The ratio $l_N/N\to 0$ as N goes to infinity.

(iii)

$A\left(\rho \right)$ is nonsingular for any $\left|\rho \right|<1$.

(iv)

The matrices $W_0$ and $A^{-1}\left(\rho \right)$ are uniformly bounded in both row and column sums in absolute value.

Assumption 3.

$K(·)$ is a nonnegative continuous even function. Let ${\mu }_l=\int k\left(\upsilon \right){\upsilon }^{l}d\upsilon$, ${\upsilon }_l=\int k^2\left(\upsilon \right){\upsilon }^{l}d\upsilon$; then ${\mu }_l={\upsilon }_l=0$ for any positive odd number. Meanwhile, ${\mu }_0=1$, ${\mu }_2\ne 0$.

Assumption 4.

If $N\to \infty$, $T\to \infty$ and $h\to 0$, then $NTh\to \infty$.

Assumption 5.

There is an unique ${\theta }_0$ to make the model (1) tenable.

Assumption 6.

${\varphi }_{Q_{0}Z}{\varphi }_{ZZ}^{-1}{\varphi }_{Q_{0}Z}^{\prime }+{\varphi }_{Q_0{Q}_0}\ge 0$, where

$${\varphi }_{Q_{0}Z}=\underset{N\to \infty }{lim}\underset{T\to \infty }{lim}\frac{1}{NT}{Q}_0^{\prime }{(I-S)}^{\prime }(I-S)Z,{\varphi }_{ZZ}=\underset{N\to \infty }{lim}\underset{T\to \infty }{lim}\frac{1}{NT}{Z}^{\prime }{(I-S)}^{\prime }(I-S)Z,$$

$${\varphi }_{Q_0{Q}_0}=\underset{N\to \infty }{lim}\underset{T\to \infty }{lim}\frac{1}{NT}{Q}_0^{\prime }{(I-S)}^{\prime }(I-S)Q_0,Q_0=G_0(Z{\alpha }_0+X{\beta }_0),G_0=W{A}^{-1}\left({\rho }_0\right).$$

Remark 1.

Assumption 1 provides the essential features of the regressors and disturbances for the model (see [47]). Assumption 2 concerns the basic features of the spatial weights matrix and the parallel Assumptions 2–5 of [12]. Assumption 2(i) is always satisfied if $l_N$ is a bounded sequence. We allow $l_N$ to be divergent, but at a rate smaller than N, as specified in Assumption 2(ii). Assumption 2(iii) guarantees that model (2) has an equilibrium given by (3). Assumption 2(iv) is also assumed in [12], and it limits the spatial correlation to some degree but facilitates the study of the asymptotic properties of the spatial parameter estimators. Assumptions 3 and 4 concern the kernel function and bandwidth sequence. Assumption 5 offers a unique identification condition, and Assumption 6 is necessary for proof of asymptotic normality.

In order to prove large sample properties of estimators, we need introduce the following useful lemmas. Before that, we simplify the model (2) and obtain reduced form equation of Y as follows:

$$Y=X{\beta }_0\left(u\right)+Z{\alpha }_0+\rho G(Z{\alpha }_0+X{\beta }_0\left(u\right))+A^{-1}\left({\rho }_0\right)(Ub+\epsilon ),$$

where $A^{-1}\left({\rho }_0\right)=I+{\rho }_0{G}_0$; $G_0=W{A}^{-1}\left({\rho }_0\right)$. The above equations are frequently used in a later derivation.

Lemma 1.

Let $\{(X_1,Y_1),(X_2,Y_2),\cdots ,(X_n,Y_n)\}$ be i.i.d. random vectors, where $Y_i$ are scalar random variables. Further, assume that ${E\left|y\right|}^{r}f(x,y)dy<\infty$, where f denotes the joint density of $(X,Y)$. Let K be a bounded positive function with a bounded support, satisfying a Lipschitz condition. Given that $n^{2\tau -1}h\to \infty$ for some $\tau <1-r^{-1}$,

$$\underset{x}{sup}|\frac{1}{n}\sum _{i=1}^n[K_h(X_i-x)Y_i-E{K}_h(X_i-x)Y_i]|=O_p{\left(\{\frac{log(1/h)}{nh}\right\}}^{1/2}).$$

The proof can be found in [48].

Lemma 2.

Under Assumptions 1–4, we have

$$S_u\stackrel{p}{⟶}F\left(u\right)\left(\begin{array}{cc}{\mathsf{\Omega }}_{11}\left(u\right)& 0\\ 0& {\mu }_2{\mathsf{\Omega }}_{11}\left(u\right)\end{array}\right),$$

where $F\left(u\right)=\underset{T\to \infty }{lim}\frac{1}{T}{\sum }_{i=1}^T{f}_t\left(u\right)$ and ${\mathsf{\Omega }}_{11}=E\left(x_{it}{x}_{it}^{{}^{\prime }}\right|u_{it}=u)$, $i=1,2,\cdots ,N,t=1,2,\cdots ,T$.

Proof See the Appendix A.

Lemma 3.

Under Assumptions 1–4, we have${\widehat{\beta }}_{IN}\left(u\right)\stackrel{p}{⟶}{\beta }_0\left(u\right)$.

Proof See the Appendix A.

Lemma 4.

Under Assumptions 1–4, we have

(i)

$\frac{1}{NT}{Z}^{\prime }{(I-S)}^{\prime }(I-H)(I-S)Z\stackrel{p}{⟶}{\mathsf{\Omega }}_{22}\left(u\right){\mathsf{\Omega }}_{11}^{-1}\left(u\right){\mathsf{\Omega }}_{12}\left(u\right)$;

(ii)

$(I-S)X\beta \left(u\right)=op\left(1\right)$, $(I-S)Z\alpha =o_p\left(1\right)$;

(iii)

$(I-S)GX\beta \left(u\right)=O_p\left(1\right)$, $(I-S)GZ\alpha =O_p\left(1\right)$;

(iv)

$S(Ub+\epsilon )=o_p\left(1\right)$, $SG(Ub+\epsilon )=o_p\left(1\right).$

The proof can be found in [37].

Lemma 5.

Under Assumptions 1–6, ${\mathsf{\Sigma }}_{{\theta }_0}=E\left[\frac{1}{\sqrt{NT}}\frac{{\partial }^{2}lnL\left({\theta }_0\right)}{\partial \theta \partial {\theta }^{\prime }}\right]$ is a positive definite matrix and

$$\begin{array}{ccc}\hfill & \frac{1}{\sqrt{NT}}\frac{{\partial }^{2}lnL\left({\theta }_0\right)}{\partial \theta \partial {\theta }^{\prime }}=\hfill & \\ \hfill & \left(\begin{array}{cccc}\frac{{\sigma }_b^{2}E\left[Q_0^{\prime }P{Q}_0\right]}{N({\sigma }_{\epsilon 0}^2+T{\sigma }_{b0}^2){\sigma }_{\epsilon 0}^2}+\frac{{\sigma }_{b0}^2[tr\left(G_0^{\prime }{G}_0\right)+tr\left(G_0^{\prime }H{G}_0\right)]}{N{\sigma }_{\epsilon 0}^2}+\frac{tr\left(G_0^{\prime }{G}_0\right)-tr\left(G_0^2\right)}{NT}& *& *& *\\ \frac{tr\left(G_0^{\prime }H\right)}{NT}+\frac{({\sigma }_{\epsilon 0}^2+T{\sigma }_{b0}^2)[tr\left(G_0^{\prime }\right)-tr\left(G_0^{\prime }H\right)]}{NT{\sigma }_{\epsilon 0}^4}& -\frac{1}{2{\sigma }_{\epsilon 0}^4}& *& *\\ \frac{tr\left(G_0^{\prime }H\right)}{NT({\sigma }_{\epsilon 0}^2+T{\sigma }_{b0}^2)}& -\frac{1}{2{({\sigma }_{\epsilon 0}^2+T{\sigma }_{b0}^2)}^2}& -\frac{T}{2{({\sigma }_{\epsilon 0}^2+T{\sigma }_{b0}^2)}^2}& *\\ \frac{{\sigma }_{b0}^{2}E\left[Q_0^{\prime }PZ\right]}{N({\sigma }_{\epsilon 0}^2+T{\sigma }_{b0}^2){\sigma }_{\epsilon 0}^2}& 0& 0& -\frac{E\left[Z^{\prime }PZ\right]}{NT{\sigma }_{\epsilon 0}^2({\sigma }_{\epsilon 0}^2+T{\sigma }_{b0}^2)}\end{array}\right)\hfill \\ \hfill & +o_p\left(1\right).\hfill \end{array}$$

Proof See the Appendix A.

With the above lemmas, we state main results as follows. Their detailed proofs are given in the Appendix A.

Theorem 1.

Under Assumptions 1–5, we have $\widehat{\rho }-{\rho }_0=op\left(1\right)$ and $\widehat{\alpha }-{\alpha }_0=op\left(1\right)$.

Theorem 2.

Under Assumptions 1–5, we have $\widehat{\beta }\left(u\right)-{\beta }_0\left(u\right)=op\left(1\right)$.

Theorem 3.

Under Assumptions 1–5, we have ${\widehat{\sigma }}_{\epsilon }^2-{\sigma }_{\epsilon 0}^2=op\left(1\right)$ and ${\widehat{\sigma }}_b^2-{\sigma }_{b0}^2=op\left(1\right)$.

Theorem 4.

Under Assumptions 1–6, we have

$$\sqrt{NT}(\widehat{\theta }-{\theta }_0)\stackrel{D}{⟶}N(0,{\mathsf{\Sigma }}_{{\theta }_0}^{-1}+{\mathsf{\Sigma }}_{{\theta }_0}^{-1}{\mathsf{\Omega }}_{\theta }{\mathsf{\Sigma }}_{{\theta }_0}^{-1}).$$

where “$\stackrel{D}{⟶}$” represents convergence in distribution, ${\mathsf{\Sigma }}_{{\theta }_0}=-E\left[\frac{1}{\sqrt{NT}}\frac{{\partial }^{2}lnL\left({\theta }_0\right)}{\partial \theta \partial {\theta }^{\prime }}\right]$ is the average Hessian matrix (information matrix when ε and b are normally distributed, respectively), ${\mathsf{\Omega }}_{\theta }=E\left({\mathsf{\Omega }}_{{\theta }_{NT}}\right)$ and ${\mathsf{\Omega }}_{{\theta }_{NT}}=\frac{1}{\sqrt{NT}}\frac{\partial lnL\left({\theta }_0\right)}{\partial \theta }·\frac{1}{\sqrt{NT}}\frac{\partial lnL\left({\theta }_0\right)}{\partial {\theta }^{\prime }}+\frac{1}{NT}\frac{{\partial }^{2}lnL\left({\theta }_0\right)}{\partial \theta \partial {\theta }^{\prime }}$.

Theorem 5.

Under Assumptions 1–5, we have

$$\sqrt{NTh}(\widehat{\beta }\left(u\right)-{\beta }_0\left(u\right)-\phi \left(u\right))\stackrel{D}{⟶}N(0,{\gamma }^2\left(u\right)),$$

where ${\gamma }^2\left(u\right)$ and $\phi \left(u\right)$ satisfy $\phi \left(u\right)=\frac{1}{2}{h}^2{\mu }_2\ddot{\beta }\left(u\right)$ and ${\gamma }^2\left(u\right)={\upsilon }_0({\sigma }_{\epsilon ,0}^2+{\sigma }_{b,0}^2){\left[F\left(u\right){\mathsf{\Omega }}_{11}\left(u\right)\right]}^{-1},$ $F\left(u\right)=\underset{T\to \infty }{lim}\frac{1}{T}{\sum }_{t=1}^T{f}_t\left(u\right),$ ${\mathsf{\Omega }}_{11}\left(u\right)=E\left(x_{it}{x}_{it}^{{}^{\prime }}\right|u_{it}=u)$ and $\ddot{\beta }\left(u\right)$ is the second-order derivative of $\beta \left(u\right)$, $i=1,2,\cdots ,N,t=1,2,\cdots ,T$. Furthermore, if $NT{h}^5\to 0$, then we have

$$\sqrt{NTh}(\widehat{\beta }\left(u\right)-{\beta }_0\left(u\right))\stackrel{D}{⟶}N(0,{\gamma }^2\left(u\right)).$$

4. Monte Carlo Simulation

In this section, Monte Carlo simulations are presented, which were carried out to investigate the finite sample performance of PQMLEs. The sample standard deviation ($SD$) and two root mean square errors ($RMSE1$ and $RMSE2$) were used to measure the estimation performance. $RMSE1$ and $RMSE2$ are defined as:

$$\mathrm{RMSE}1={\left[\frac{1}{mcn}\sum _{i=1}^{mcn}{({\widehat{\theta }}_i-{\theta }_0)}^2\right]}^{\frac{1}{2}},\mathrm{RMSE}2={({\widehat{\theta }}_{0.5}-{\theta }_0)}^2+\frac{{\widehat{\theta }}_{0.75}-{\widehat{\theta }}_{0.25}}{1.35},$$

where $mcn$ is the number of iterations; ${\widehat{\theta }}_i(i=1,2,\cdots ,mcn$) are estimates of $\theta$ for each iteration, ${\theta }_0$ is the true value of $\theta$; ${\widehat{\theta }}_{0.25}$, ${\widehat{\theta }}_{0.5}$ and ${\widehat{\theta }}_{0.75}$ are the upper quartile, median and lower quartile of parametric estimates, respectively. For the nonparametric estimates, we took the mean absolute deviation error ($MADE$) as the evaluation criterion:

$${\mathrm{MADE}}_j=Q^{-1}\sum _{q=1}^Q|{\widehat{\beta }}_j\left(u_q\right)-{\beta }_j\left(u_q\right)|,j=1,2,\cdots ,mcn,$$

where ${\left\{u_q\right\}}_{q=1}^Q$ are Q fixed grid points in the support set of u. In the simulations, we applied the rule of thumb method of [48] to choose the optimal width and let kernel function be an Epanechnikov kernel $K\left(u\right)=\left(\frac{3}{4\sqrt{5}}\right)(1-\frac{1}{5}{u}^2)1(u^2\le 5)$ (see [33]).

We ran a small simulation experiment with $mcn=500$ and generated the simulation data from following model:

$$y_{it}=\rho {\left(W_0{Y}_t\right)}_i+x_{it}^{{}^{\prime }}\beta \left(u_{it}\right)+z_{it}^{{}^{\prime }}\alpha +b_i+{\epsilon }_{it},i=1,2,\cdots ,N,t=1,2,\cdots ,T,$$

where we assume that $x_{it}\sim U[-2,2]$, $z_{it}\sim U[-2,2]$, $u_{it}\sim U[-3,3]$, ${\epsilon }_{it}\sim i.i.d.N(0,0.5)$, $b_i\sim i.i.d.N(0,1)$, ${\beta }_1\left(u_{it}\right)=0.5u_{it}+sin\left(1.5u_{it}\right)$, ${\beta }_2\left(u_{it}\right)=u_{it}^2+0.5u_{it}$, $({\alpha }_1,{\alpha }_2)=(1,1.5)$ and $\rho =0.25,0.5$ and 0.75, respectively. Furthermore, we chose the Rook weight matrix (see [7]) to investigate the influence of the spatial weight matrix on the estimates. Our simulation results for both cases $T=10$ and $T=15$ are presented in

Table 1. The medians and SDs of MADE values for ${\beta }_1\left(u\right)$.

Parameter		T = 10				T = 15
		N = 49	N = 64	N = 81	N = 100	N = 49	N = 64	N = 81	N = 100
	Median	0.1681	0.1519	0.1415	0.1326	0.1472	0.1340	0.1337	0.1301
$\rho =0.25$	SD	0.0352	0.0337	0.0283	0.0247	0.0289	0.0248	0.0240	0.0203
	Median	0.1769	0.1617	0.1498	0.1424	0.1587	0.1475	0.1400	0.1385
$\rho =0.5$	SD	0.0377	0.0338	0.0295	0.0256	0.0298	0.0296	0.0243	0.0220
	Median	0.1553	0.1421	0.1316	0.1219	0.1384	0.1205	0.1137	0.1048
$\rho =0.75$	SD	0.0335	0.0280	0.0254	0.0239	0.0256	0.0226	0.0211	0.0200

,

Table 2. The medians and SDs of MADE values for ${\beta }_2\left(u\right)$.

Parameter		T = 10				T = 15
		N = 49	N = 64	N = 81	N = 100	N = 49	N = 64	N = 81	N = 100
	Median	0.0944	0.0864	0.0777	0.0677	0.084	0.0724	0.0704	0.0677
$\rho =0.25$	SD	0.0326	0.0313	0.0249	0.0224	0.0274	0.0230	0.0201	0.0178
	Median	0.0985	0.0891	0.0836	0.0768	0.0868	0.0785	0.074	0.0724
$\rho =0.5$	SD	0.0322	0.0304	0.0258	0.0238	0.0276	0.0234	0.0206	0.0186
	Median	0.0957	0.0863	0.0795	0.0729	0.0807	0.0743	0.0648	0.0608
$\rho =0.75$	SD	0.0314	0.0257	0.0211	0.0213	0.0254	0.0223	0.0191	0.0173

,

Table 3. The results of parametric estimates with T = 10(1).

Parameter	True	N = 49				N = 64
	Value	MEAN	SD	RMSE1	RMSE2	MEAN	SD	RMSE1	RMSE2
$\rho$	0.2500	0.2836	0.0685	0.0859	0.0657	0.2720	0.0662	0.0851	0.0622
${\sigma }_{\epsilon }^2$	0.5000	0.5949	0.0625	0.1136	0.0608	0.5844	0.0619	0.1046	0.0592
${\sigma }_b^2$	1.0000	0.9070	0.1971	0.2176	0.1979	0.9226	0.1786	0.1944	0.1874
${\alpha }_1$	1.0000	0.8755	0.0511	0.0072	0.0057	0.8879	0.0439	0.0065	0.0042
${\alpha }_2$	1.5000	1.5721	0.0111	0.0155	0.0097	1.5983	0.0039	0.0126	0.0052
$\rho$	0.5000	0.5236	0.0532	0.0896	0.0599	0.5222	0.0480	0.0878	0.0544
${\sigma }_{\epsilon }^2$	0.5000	0.6026	0.0729	0.1258	0.0688	0.5913	0.0644	0.1116	0.0639
${\sigma }_b^2$	1.0000	0.8785	0.2026	0.2360	0.2107	0.9059	0.1819	0.2046	0.1999
${\alpha }_1$	1.0000	0.8690	0.0058	0.0076	0.0098	0.9145	0.0036	0.0049	0.0044
${\alpha }_2$	1.5000	1.5764	0.0079	0.0172	0.0286	1.6692	0.0056	0.0073	0.0058
$\rho$	0.7500	0.8083	0.0253	0.0635	0.0295	0.8110	0.0246	0.0657	0.0294
${\sigma }_{\epsilon }^2$	0.5000	0.6059	0.0702	0.1270	0.0734	0.6048	0.0659	0.1237	0.0679
${\sigma }_b^2$	1.0000	0.8864	0.2063	0.2352	0.2212	0.8623	0.1694	0.2181	0.1944
${\alpha }_1$	1.0000	0.8421	0.0472	0.0091	0.0009	0.8935	0.0442	0.0061	0.0007
${\alpha }_2$	1.5000	1.5260	0.0023	0.0249	0.0002	1.5152	0.0012	0.0113	0.0001

,

Table 4. The results of parametric estimates with T = 10(2).

Parameter	True	N = 81				N = 100
	Value	MEAN	SD	RMSE1	RMSE2	MEAN	SD	RMSE1	RMSE2
$\rho$	0.2500	0.3033	0.0578	0.0785	0.0620	0.3051	0.0527	0.0762	0.0531
${\sigma }_{\epsilon }^2$	0.5000	0.5585	0.0456	0.0741	0.0432	0.5530	0.0389	0.0657	0.0369
${\sigma }_b^2$	1.0000	0.9289	0.1596	0.1745	0.1691	0.9616	0.1524	0.1569	0.1477
${\alpha }_1$	1.0000	0.9468	0.0364	0.0032	0.0075	1.0757	0.0363	0.0044	0.0036
${\alpha }_2$	1.5000	1.5410	0.0402	0.0476	0.0089	1.5246	0.0373	0.0057	0.0039
$\rho$	0.5000	0.5112	0.0396	0.0814	0.0465	0.5033	0.0339	0.0808	0.0414
${\sigma }_{\epsilon }^2$	0.5000	0.5765	0.0505	0.0916	0.0537	0.5679	0.0436	0.0806	0.0457
${\sigma }_b^2$	1.0000	0.9203	0.1575	0.1763	0.1481	0.9352	0.1480	0.1613	0.1565
${\alpha }_1$	1.0000	0.8403	0.0082	0.0092	0.0017	0.9422	0.0014	0.0033	0.0015
${\alpha }_2$	1.5000	1.4713	0.0034	0.0255	0.0008	1.4743	0.0032	0.0033	0.0007
$\rho$	0.7500	0.8094	0.0208	0.0629	0.0218	0.8101	0.0178	0.0627	0.0215
${\sigma }_{\epsilon }^2$	0.5000	0.5823	0.0548	0.1057	0.0600	0.5821	0.0522	0.0933	0.052
${\sigma }_b^2$	1.0000	0.8966	0.1520	0.1837	0.1669	0.9054	0.1464	0.1741	0.1610
${\alpha }_1$	1.0000	0.8568	0.0061	0.0083	0.0088	0.9889	0.0035	0.0006	0.0005
${\alpha }_2$	1.5000	1.6519	0.0022	0.0205	0.0231	1.4913	0.0017	0.0001	0.0001

,

Table 5. The results of parametric estimates with T = 15(1).

Parameter	True	N = 49				N = 64
	Value	MEAN	SD	RMSE1	RMSE2	MEAN	SD	RMSE1	RMSE2
$\rho$	0.2500	0.2751	0.0670	0.0801	0.0608	0.2740	0.0637	0.0780	0.0614
${\sigma }_{\epsilon }^2$	0.5000	0.5761	0.0556	0.0942	0.0548	0.5626	0.0523	0.0815	0.0385
${\sigma }_b^2$	1.0000	0.9374	0.2005	0.1989	0.1046	0.9363	0.1752	0.1861	0.1872
${\alpha }_1$	1.0000	1.1449	0.0356	0.0054	0.0049	0.9441	0.0335	0.0032	0.0025
${\alpha }_2$	1.5000	1.4568	0.0052	0.0210	0.0073	1.5855	0.0034	0.0031	0.0019
$\rho$	0.5000	0.5749	0.0517	0.0910	0.0613	0.5766	0.0433	0.0879	0.0439
${\sigma }_{\epsilon }^2$	0.5000	0.5922	0.0596	0.1098	0.0614	0.5779	0.0598	0.0982	0.0515
${\sigma }_b^2$	1.0000	0.9193	0.2160	0.2303	0.2313	0.9126	0.1729	0.1935	0.1814
${\alpha }_1$	1.0000	1.0623	0.0202	0.0036	0.0101	1.0652	0.0021	0.0038	0.0038
${\alpha }_2$	1.5000	1.3248	0.0026	0.0039	0.0307	1.3651	0.0023	0.0043	0.0182
$\rho$	0.7500	0.8118	0.0270	0.0674	0.0331	0.8146	0.0226	0.0684	0.0264
${\sigma }_{\epsilon }^2$	0.5000	0.6058	0.0707	0.1272	0.0671	0.5929	0.0611	0.1112	0.0622
${\sigma }_b^2$	1.0000	0.8746	0.2018	0.2173	0.2024	0.8788	0.1688	0.2076	0.1708
${\alpha }_1$	1.0000	1.0408	0.0093	0.0024	0.0008	1.0420	0.0046	0.0024	0.0016
${\alpha }_2$	1.5000	1.4716	0.0006	0.0017	0.0002	1.4956	0.0004	0.0015	0.0008

and

Table 6. The results of parametric estimates with T = 15(2).

Parameter	True	N = 81				N = 100
	Value	MEAN	SD	RMSE1	RMSE2	MEAN	SD	RMSE1	RMSE2
$\rho$	0.2500	0.2662	0.0545	0.0714	0.0576	0.2646	0.0493	0.0664	0.0525
${\sigma }_{\epsilon }^2$	0.5000	0.5566	0.0441	0.0717	0.0397	0.5490	0.0392	0.0627	0.0330
${\sigma }_b^2$	1.0000	0.9508	0.1555	0.1629	0.1524	0.9772	0.1507	0.1522	0.1565
${\alpha }_1$	1.0000	0.9450	0.0202	0.0031	0.0055	0.9862	0.0008	0.0106	0.0024
${\alpha }_2$	1.5000	1.5622	0.0413	0.0030	0.0169	1.3702	0.0332	0.0028	0.0017
$\rho$	0.5000	0.5750	0.0343	0.0733	0.0448	0.4930	0.0297	0.0063	0.0204
${\sigma }_{\epsilon }^2$	0.5000	0.5633	0.0488	0.0800	0.0476	0.5364	0.0366	0.0602	0.0433
${\sigma }_b^2$	1.0000	0.9311	0.1517	0.1634	0.1525	0.9532	0.0880	0.1337	0.0935
${\alpha }_1$	1.0000	0.9822	0.0351	0.0026	0.0018	0.9673	0.0283	0.0013	0.0013
${\alpha }_2$	1.5000	1.4033	0.0029	0.0020	0.0005	1.4980	0.0203	0.0011	0.0007
$\rho$	0.7500	0.8156	0.0205	0.0687	0.0235	0.8066	0.0153	0.0543	0.0137
${\sigma }_{\epsilon }^2$	0.5000	0.5904	0.0522	0.0974	0.0564	0.5806	0.0501	0.0949	0.0512
${\sigma }_b^2$	1.0000	0.8806	0.1486	0.1783	0.1680	0.9254	0.1369	0.1635	0.1597
${\alpha }_1$	1.0000	1.0338	0.0049	0.0023	0.0077	0.9856	0.0026	0.0002	0.0005
${\alpha }_2$	1.5000	1.3659	0.0326	0.0018	0.0180	1.4977	0.0008	0.0001	0.0001

.

By observing the simulation results in Table 1, Table 2, Table 3, Table 4, Table 5 and Table 6, one can obtain the following findings: (1) The $RMSE1s$ and $RMSE2s$ for $\widehat{\rho }$, ${\widehat{\sigma }}_{\epsilon }^2$, ${\widehat{\alpha }}_1$ and ${\widehat{\alpha }}_2$ were fairly small for almost all cases, and they decreased as N increased. The $SDs$ and $RMSEs$ for ${\widehat{\sigma }}_b^2$ are not negligible for small sample sizes, but they decreased as N increased. (2) For fixed T, as N increased, the $SDs$ for $\widehat{\rho }$, ${\widehat{\sigma }}_{\epsilon }^2$, ${\widehat{\sigma }}_b^2$, ${\widehat{\alpha }}_1$ and ${\widehat{\alpha }}_2$ decreased for all cases. $RMSE1s$ and $RMSE2s$ for ${\widehat{\alpha }}_1$ and ${\widehat{\alpha }}_2$ decreased rapidly, whereas the $RMSEs$ for estimates of the others parameters did not change much. For fixed N, as T increased, the behavior of the estimates of parameters was similar to the case where N changed under the fixed T. (3) The $SDs$ for $MADEs$ of varying coefficient functions ${\beta }_1\left(u\right)$ and ${\beta }_2\left(u\right)$ decreased as T or N increased. Combined with the above three findings, we conclude that the estimates of the parameter and unknown varying coefficient functions were convergent.

Figure 1 and Figure 2 present the fitting results and 95% confidence intervals of ${\beta }_1\left(u\right)$ and ${\beta }_2\left(u\right)$ under $\rho =0.5$, respectively, where the short dashed curves are the average fits over 500 simulations ${\widehat{\beta }}_s\left(u\right)$$(s=1,2)$ by PQMLE, the solid curves are the true values of ${\beta }_s\left(u\right)$$(s=1,2)$ and the two long dashed curves are the corresponding 95% confidence bands. By observing every subgraph in Figure 1 and Figure 2, we can see that the short dashed curve is fairly close to solid curve and the corresponding confidence bandwidth is narrow. This illustrates that the nonparametric estimation procedure works well for small samples. To save space, we do not present the cases $\rho =0.25$ and $\rho =0.75$ because they had similar results as the case $\rho =0.5$.

5. Real Data Analysis

We applied the housing prices of Chinese city data to analyze the proposed model with real data. The data were obtained from the China Statistical Yearbook, the China City Statistical Yearbook and the China Statistical Yearbook for Regional Economies. Based on the panel data of related variables of 287 cities at/above the prefecture level (except the cities in Taiwan, Hong Kong and Macau) in China from 2011 to 2018, we explored the influencing factors of housing prices of Chinese cities by PLVCSARPM with random effects.

Taking [49,50,51] as references, we collected nine variables related to housing prices of China cities, including each city’s average selling price of residential houses (denoted by HP, yuan/sq.m), the expectation of housing price trends (EHP, %), population density (POD, person), annual per capita disposable income of urban households (ADI, yuan), loan–to–GDP ratio (MON, %), natural growth rate of population (NGR, %), sulphur dioxide emission (SDE, 10,000 tons), area of grassland (AOG, 10,000 hectares) and the value of buildings completed (VBC, yuan/sq.m). According to result of non-linear regression, the established model is given by

$$y_{it}=\rho {\left(W_0{Y}_t\right)}_i+x_{it}^{{}^{\prime }}\beta \left(u_{it}\right)+z_{it}^{{}^{\prime }}\alpha +b_i+{\epsilon }_{it},1\le i\le 287,1\le t\le 8,$$

(14)

where $y_{it}$ represents the ith observation of ln(HP) at time t; $z_{it}$ represents the ith observation of EHP at time t; $x_{it}$ means the ith observation of POD, ln(ADI), MON, NGR, SDE and AOG at time t, respectively, $u_{it}$ represents the ith observation of ln(VBC) at time t.

In order to transfer the asymmetric distribution of POD to nearly uniform distribution on (0,1), we set

$$POD=\left(PO{D}^{1/3}-min\left(PO{D}^{1/3}\right)\right)/\left(max\left(PO{D}^{1/3}\right)-min\left(PO{D}^{1/3}\right)\right).$$

The spatial weight matrix $W_0$ we adopted is calculated as follows:

$$w_{ij}=exp\left(-\parallel s_i-s_j\parallel \right)/\sum _{k=1}^{264}exp\left(-\parallel s_i-s_k\parallel \right),$$

where $s_i=(longitud{e}_i,latitud{e}_i).$

Figure 3 presents the estimation results and corresponding 95% confidence intervals of varying coefficient functions ${\beta }_s(s=1,2,\cdots ,6)$, where the black dashed line curves are the average fits over 500 simulations and the red solid lines are the corresponding 95% confidence bands. It can be seen from Figure 3 that all covariates variables have obvious non-linear effects on housing prices of Chinese cities.

The estimation results of parameters in the model (14) are reported in

Table 7. The estimation results of unknown parameters in the model (14).

	$\widehat{\mathit{\rho }}$	$\widehat{\mathit{\alpha }}$	${\widehat{\mathit{\sigma }}}_{\mathit{b}}^2$	${\widehat{\mathit{\sigma }}}_{\mathit{\epsilon }}^2$
Estimators	$0.3820{}^{***}$	$0.9869{}^{***}$	$0.1375{}^{**}$	$0.0163{}^{**}$
SD	0.0243	0.0397	5.0344	0.3687

Notes: ** and *** are significant at the significance levels of 5% and 1%, respectively.

. It can be seen from Table 7: (1) all estimates of parameters are all significant; (2) spatial correlation coefficient $\widehat{\rho }=0.3820>0$, which means there exists a spatial spill over effect for housing prices of Chinese cities; (3) $\widehat{\alpha }=0.9869>0$ indicates that the expectation of housing price trends (EHP) has a promotional effect on housing prices of Chinese cities; (4) ${\widehat{\sigma }}_{\epsilon }^2=0.0163$ shows that the growth of housing prices in different regions is relatively stable and is less affected by external fluctuations.

6. Conclusions

In this paper, we proposed PQMLE of PLVCSARPM with random effects. Our model has the following advantages: (1) It can overcome the “curse of dimensionality” in the nonparametric spatial regression model effectively. (2) It can simultaneously study the linear and non-linear effects of coveriates. (3) It can investigate the spatial correlation effect of response variables. Under some regular conditions, consistency and asymptotic normality of the estimators for parameters and varying coefficient functions were derived. Monte Carlo simulations showed the proposed estimators are well behaved with finite samples. Furthermore, the performance of the proposed method was also assessed on a set of asymmetric real data.

This paper only focused on the PQMLE of PLVCSARPM with random effects. In future research, we may try to extend our method to more general models, such as a partially linear varying coefficient spatial autoregressive model with autoregressive disturbances. In addition, we also need study the issues of Bayesian analysis, variable selection and quantile regression in these models.