Discussion Paper No. 104
Estimation of Partially Linear Spatial
Autoregressive Models withAutoregressive Disturbances
Takaki Sato

October, 2019

 Data Science and Service Research
Discussion Paper

Center for Data Science and Service Research
Graduate School of Economic and Management

Tohoku University27-1 Kawauchi, Aobaku

Sendai 980-8576, JAPANEstimation of Partially Linear Spatial Autoregressive Models withAutoregressive Disturbances
Takaki Sato ∗Abstract
This study considers semiparametric partially linear spatial autoregressive models with autoregressive
disturbances that contain an unspeciﬁed nonparametric component and allow for spatial lags in both
the dependent variables and disturbances. Having the nonparametric function approximated by basis
functions, we propose a three-step estimation procedure for the proposed model. We also establish theconsistency and asymptotic normality of the proposed estimators. Then, the ﬁnite sample performancesof the proposed estimators are examined using Monte Carlo simulations. As an empirical application,
we use the proposed model and estimation method to analyze Boston housing price data to evaluate theeﬀect of air pollution on the value of owner-occupied homes.

Keywords: Partially linear models, Series estimation, Spatial econometrics, Instrumental variables.

1

Introduction

Recently, the spatial autoregressive (SAR) model proposed by Clif and Ord (1973) has received increasing
attention in both theoretical and applied econometrics research. Speciﬁcally, the data in the ﬁeld of regional,
urban, and environmental economics usually show the spatial dependency of cross-sectional units and SAR
models are used to capture this dependency. The class of SAR models is extended by considering spatial
interaction eﬀects in both the dependent variables and disturbances. We call these models SAR models withspatial autoregressive disturbances (SARAR).
Anselin (1988) and Lee (2004) propose the (quasi) maximum likelihood (ML) to estimate such parametricspatial econometric models. However, one drawback of ML estimation is the computational load when thesample size is large, because there is no closed-form expression of ML estimators; therefore, it is necessary
to calculate the determinant of a large matrix, whose size depends on the sample size. Another approach∗ Advanced Institute for Yotta Informatics, Tohoku University, Sendai 980-8577, Japan. satotaka@riec.tohoku.ac.jp

1

for the estimation of spatial econometric models consists of moment-based estimations. Kelejian and Prucha(1998, 2010) introduce generalized spatial two-stage least squares (2SLS) estimation methods, while Lee andLiu (2010) consider the generalized method of moments (GMM) estimation methods.
To avoid mis-speciﬁcation of the data generating process in parametric models, the semiparametric extensions of spatial econometric models have received signiﬁcant attention owing to the simple interpretationof parametric terms and the ﬂexibility of nonparametric terms. A popular semiparametric regression modelis the partially linear one, which contains explanatory variables nonlinearly rerated with dependent variables.
As semiparametric extensions of the SAR models, Su and Jin (2010) and Du et al. (2018) propose partiallylinear SAR (PL-SAR), while Su (2012) considers partially linear SARAR (PL-SARAR) models. Zhang and
Sun (2015) further study the spatial dynamic panel extension of PL-SAR models. Another semiparametricextension is the varying coeﬃcient model, in which the impact of some explanatory variables depends on spatial units. Zhang and Shen (2015) consider semiparametric varying coeﬃcient-speciﬁed spatial panel models
and Hoshino (2018) proposes functional coeﬃcient SAR models with endogenous regressors.
A popular method for estimating nonparametric terms in regression models is the kernel approach. Su(2012) applies kernel methods and proposes the estimation method for the PL-SARAR models in which the
nonparametric terms are proﬁled out. However, as the sample size increases, the computational load of theseestimation methods increases signiﬁcantly, making them less manageable. Another estimation method for
nonparametric terms is series estimation. One advantage of series methods is their computational simplicity. As such, we apply moment-based estimation methods for the estimation of nonparametric terms byapproximating the nonparametric terms using basis functions such as polynomials and splines.
We consider the moment-based estiamtion method for PL-SARAR models for computational simplicity.
Accordingly, we propose a three-step estimation procedure by applying the 2SLS and nonlinear least squares(NLS) methods for the parametric terms and series methods for the nonparametric term in the proposedmodel. The consistency and asymptotic normality of the proposed estimators are established and the small
sample properties of the proposed estimators are then evaluated.
As an empirical analysis, we apply the SARAR and PL-SARAR models to Boston land price data toevaluate the causal eﬀect of air pollution on housing prices. In the model, the dependent variable is the medianvalue of owner-occupied homes and the explanatory variable is the nitrogen oxide (NOX) concentration. Our
empirical ﬁndings are as follows. First, housing prices show spatial correlations even after we control for the
potential determinants of housing prices. Second, air pollution has strong negative eﬀects on housing prices
in both the parametric and semiparametric models. Finally, the eﬀect of air pollution of housing prices is
not linear and the negative eﬀect increases signiﬁcantly when the proportion of NOX in the air is above athreshold value.
2

The rest of paper proceeds as follows. We introduce PL-SARAR models and propose a three-step estimation method in section 2. The asymptotic properties of the proposed estimators are established in section 3.
Section 4 examines the small sample properties of the proposed estimators using Monte Carlo simulations.
In section 5, we apply the proposed models to Boston land price data to investigate the empirical propertiesof the proposed model. Section 6 presents the concluding remarks. The proofs of Lemmas and Theorems areprovided in the Appendix.
Notation: We use In to denote an n × n identity matrix. For matrix An , ||An || denotes its Frobenius
norm: ||An || = {tr(An An )}1/2 , where tr(·) is the trace operator. When An is a symmetric matrix, γmax (An )and γmin (An ) denote the largest and smallest eigenvalues of An , respectively.

2

Model Speciﬁcation and Estimation

Let us consider the following PL-SARAR models:

yn,i
un,i

=

=

λ0
ρ0

n

j=1
n




wn,i,j yn,j + xn,i β0 + g0 (sn,i ) + un,i ,(1)

mn,i,j un,j + εn,i ,

j=1
(1)

(d )

where n is the number of spatial units, yn,i is an observed dependent variable, xn,i = (xn,i , . . . , xn,ix ) is adx × 1 vector of exogenous regressors, sn,i is a nonparametric regressor, g0 (·) is an unknown function, εn,i isan independently and identically distributed (i.i.d.) disturbance with mean zero and variance σ02 , and wn,i,j
and mn,i,j are the (i, j)th elements of predetermined n × n spatial weight matrices Wn and Mn , respectively.
Scalar parameters λ0 and ρ0 are SAR parameters and β0 is a coeﬃcient vector.
We apply the series approximation method to estimate the nonparametric term. Let {pk (·) : k = 1, 2, . . .}be a sequence of basis functions such as polynomials, splines, and Fourier series. We assume that nonparametric function g0 (sn,i ) can be approximated by P K (sn,i ) α0 , where P K (·) = (p1 (·), . . . , pK (·)) , K is thenumber of basis functions, and α0 is a K × 1 vector of parameters. Therefore, the series approximation errorof the nonparametric function is given by:

vn,i = g0 (sn,i ) − P K (sn,i )α0 ,3and model (1) is expressed as

=

yn,i
un,i

=

λ0
ρ0

n

j=1
n




wn,i,j yn,j + xn,i β0 + P K (sn,i )α0 + vn,i + un,i ,(2)

mn,i,j un,j + εn,i .

j=1

For notational simplicity, we consider the following matrix notation of the proposed model. Let Yn =(yn,1 , . . . , yn,n ) , Xn = (xn,1 , . . . , xn,n ) , Bn = (Wn Yn , Xn ), δ0 = (ρ0 , β0 ) , Pn = (P K (sn,1 ), . . . , P K (sn,n )) ,Vn = (vn,1 , . . . , vn,n ) , and εn = (εn,1 , . . . , εn,n ) . When In − ρ0 Mn are nonsingular, model (2) is rewritten
as:

Yn

=

Bn δ0 + Pn α0 + Vn + (In − ρ0 Mn )−1 εn .

(3)

For the estimation of the parameters in model (3), we propose a three-step estimation procedure. In theﬁrst step, we apply 2SLS to model (3) to estimate δ0 because the spatial lagged dependent variable, Wn Yn , iscorrelated with the error term, (In − ρ0 Mn )−1 εn . In the second step, we estimate the coeﬃcient of the basis
function, α0 , and the unknown function, g0 (·), by ordinary least squares (OLS). In the third step, the spatial
autoregressive parameter and variance of disturbances, ρ0 and σ02 , respectively, are estimated by applying
NLS to the residuals obtained in the ﬁrst and second steps.
The ﬁrst step is the estimation of parameter δ0 by 2SLS because the correlation of the spatial lagged
dependent variable and the error term leads to the inconsistency of the OLS estimator (see, e.g., Kelejian
and Prucha (1998)). Let Zn be an n × dz matrix of instrumental variables. For example, we may use matrices(Xn , Wn Xn , Wn Wn Xn ) as instrumental variables.
Following Zhang and Sun (2015) and Du et al. (2018), we partial out the series approximation. LetΠn = Pn (Pn Pn )−1 Pn denote the projection matrix onto the space spanned by Pn . Then, we obtain:= (In − Πn )Bn δ0 + (In − Πn )Vn + (In − Πn )(In − ρMn )−1 εn .

(In − Πn )Yn(4)

Applying 2SLS to model (4) with instrument variables Zn , we propose the following 2SLS estimator forparameter δ0 :








δ̂ = (Bn (In − Πn )Hn (In − Πn )Bn )−1 Bn (In − Πn )Hn (In − Πn )Yn ,where Hn = Zn (Zn Zn )−1 Zn .

4

In the second step, we consider the estimation of the coeﬃcient on the series approximation, α0 , byapplying the OLS method and derive the estimator of the unknown function, g0 (·). Using OLS, we obtainthe following estimator for α and g0 (sn,i ):α̂=(Pn Pn )−1 Pn (Yn − Bn δ̂),ĝ(sn,i )=P K (sn,i )α̂,where δ0 is the 2SLS estimator obtained in the ﬁrst step.
The third step represents the estimation of the spatial autoregressive parameter and the variance of the
disturbances, ρ0 and σ02 , respectively by NLS. Let un = Wn un , un = Wn un and εn = Wn εn . Moreover, wedenote the i-th elements of un , un , un , and εn by un,i , un,i , un,i , and εn,i , respectively.
The spatial correlation of the disturbance term indicates the following moment condition:un − ρun

=

εn ,

(5)

un − ρun

=

εn .

(6)

Following Kelejian and Prucha (1999), we deﬁne the two matrices for the NLS estimation based on (5)and (6), respectively:
⎛
Gn

=

⎜
1⎜
⎜
n⎜
⎝

=

n

2

n

n

−

i=1 E(un,i un,i )i=1 E(un,i un,i )2
i=1 E(un,i un,i + un,i )

⎛

gn

2

⎞

−

⎞

n

2
i=1 E(un,i )n2
i=1 E(un,i )tr(Mn Mn )i=1 E(un,i un,i )

0

n

−
n

⎟
⎟
⎟,
⎟
⎠

n
2
i=1 E(un,i )⎜⎟⎟1⎜⎜ n E(u2 ) ⎟ .
⎜
⎟
n,i
i=1
n⎝
⎠
n
E(u
u
)
n,i n,i
i=1

(7)

(8)

We derive the objective function for the NLS estimation by replacing the disturbances in (7) and (8) with
ˆ = W û . Moreover, we denote the i-th
the sample moments. Let ûn = Yn − Bn δ̂ − Pn α̂, ûn = Wn ûn and unn nˆ by û , û , and uˆ , respectively. We also deﬁneelements of ûn , ûn , and u
n
n,i
n,i
n,i
⎛
Ĝn

=

⎜
1⎜
⎜
n⎜
⎝

2
2
n

n

n

2
− i=1 ûn,i2
n

i=1 ûn,i ûn,inˆ
i=1 ûn,i un,i2ˆi=1 (ûn,i un,i + ûn,i )

5

−

−
n

ˆ
i=1 un,iˆ

i=1 ûn,i un,i⎞ntr(Mn Mn )
0

⎟
⎟
⎟,
⎟
⎠

⎛
ĝn

=

⎜
1⎜
⎜
n⎜
⎝

n

2
i=1 ûn,inn2
i=1 ûn,i

⎞
⎟
⎟
⎟.
⎟
⎠

i−1 ûn,i ûn,iLet η = (ρ, ρ2 , σ 2 ) . Then, the NLS estimators for ρ and σ 2 , ρ̂ and σ̂ 2 , respectively, are deﬁned as theminimizers of (Ĝn − gn η) (Gn − gn η). Therefore, the third step estimator is deﬁned by(ρ̂, σ̂ 2 ) = argmin (Ĝn − ĝn η) (Ĝn − ĝn η) .

3

Asymptotic Properties

Here, we consider the asymptotic properties of the proposed estimators. We introduce the following assumptions.
Assumption 1
1. All the diagonal elements of Wn and Mn are zero.
2. Matrices In − λWn and In − ρMn are nonsingular for all |λ| < 1 and |ρ| < 1.
3. The row and column sums of matrices Wn , Mn , (In −λ0 Wn )−1 and (In −ρ0 Mn )−1 are uniformly boundedin absolute value.
Assumption 2

Disturbance εi,n is i.i.d. with E(εi,n ) = 0 and V (εi,n ) = σ02 . Moreover, the disturbancehas a ﬁnite fourth moment.
Assumption 3
1. Exogenous regressors Xn are non-stochastic and the elements of Xn are uniformly bounded in absolutevalue.
2. Instrumental variables Zn are non-stochastic and the elements of Zn are uniformly bounded in absolutevalue.
3. Nonparametric regressor Sn = (sn,1 , . . . , sn,n ) is non-stochastic and the set of possible values for sn,i ,S, is a compact space.

6

Assumption 4
1. There exist α0 ∈ RK and rs > 0 so that sups∈S |P K (s) α0 − f (s)| = O(K −rs ) for each K.
2. sups∈S ||P K (s)|| = O(K 1/2 ).
3.

√

2

nK −rs → 0 and Kn → 0 as n → ∞.

4. There exist constants cPn and cPn so that 0 < cPn < γmin

Pn Pnn≤ γmaxPn Pnn< cPn < ∞.

5. g0 (s) is uniformly bounded in absolute value.
Assumption 5

Let B̃n = Wn (In − λ0 Wn )−1 (Bn + g0 (Sn )). There exist constants cB̃n and cB̃n so that

0 < cB̃n < γminB̃n B̃nn≤ γmaxB̃n B̃nn< cB̃n < ∞.

We deﬁne

Σn,1

=

Σn,2

=

Assumption 6

1 B̃ (In − Πn )Hn (In − Πn )B̃n ,
n n

1 
B (In − Πn )Hn (In − Πn )(I − ρ0 Mn )−1 (I − ρ0 Mn ) −1 (In − Πn )Hn (In − Πn )Bn .
n n
Σ1 = limn→∞ Σn,1 and Σ2 = limn→∞ Σn,2 exist and are bounded away from zero andinﬁnity.
Assumption 7

There exists constant cGn so that 0 < cGn < γmin (Gn Gn ).

Assumption 1.1 leads to the normalization of the proposed model and Assumption 1.2 to the existence
condition of the model. We say that the row sums of matrix An are uniformly bounded in absolute value if
there exists constant cA so thatmax

1≤i≤n,n≥1n|an,i,j | < cA ,j=1

where an,i,j is the (i, j)th element of An . The uniform boundedness of column sums is similarly deﬁned.
Assumption 1.3 limits the spatial correlation between the elements of Yn and εn . Assumption 2 provides theessential features of the disturbances. Assumption 3 is the standard set of assumptions in spatial econometrics
literatures. Assumption 4.1 indicates the approximation error reduction at K −rs , assumption 4.2 imposesa restriction on the basis functions, assumption 4.3 ensures that the series approximation bias does not
aﬀect the limiting distribution of the proposed estimators, and assumptions 4.4 and 4.5 are required for thederivation of the asymptotic properties of the proposed estimators. Assumption 5 limits spatial correlation toa certain degree and is required to establish the asymptotic properties of the proposed estimator. Assumption
7

6 is required to derive the limiting distribution of the ﬁrst-step estimator. Assumption 7 is required for theidentiﬁability of the third-step nonlinear estimator.
First, we consider the asymptotic behaviors of the ﬁrst-step estimator, δ̂. The limiting distribution of thisestimator is centered at δ0 and is asymptotically normal.
Theorem 1.

If Assumptions 1–6 hold, then,√−1n(δ̂ − δ0 ) −→ N (0, σ02 Σ−1
1 Σ2 Σ1 ).
dSecond, we consider the asymptotic properties of the second-step estimators, α̂ and ĝ(·). We deﬁne:

σ 2 (s) = σ 2 (P K (s)(Pn Pn )−1 Pn (In − ρ0 Mn )−1 (In − ρ0 Mn ) −1 Pn (Pn Pn )−1 P K (s)).
Then, the convergence rates of |α̂ −α0 | and sups |ĝ(s)−g0 (s)| are derived. Moreover, the limiting distributionof estimator ĝ(·) is centered at g0 (·) and asymptotically normal for a given s ∈ S.
Theorem 2.

If Assumptions 1–6 hold, then,

1. α̂ = α0 + Op (


K/N + K −rs ).

√2. sups |ĝ(s) − g0 (s)| = Op (K/ n + K (1−2rs )/2 ).
d

3. (ĝ(s) − g0 (s)) −→ N (0, σ 2 (s)).
Finally, we show the consistency of the third-step estimator.
Theorem 3.

If Assumptions 1–7 hold, then,pρ̂ −→ ρ0 ,p

σ̂ 2 −→ σ02 .

4Monte Carlo Simulation

Here, we examine the small sample performances of the proposed three-step estimators through a set ofsimulation experiments. We consider the following data generating process for the Monte Carlo simulations:

yn,i

=

λ0

n


wn,i,j yn,j + xn,i β0 + g0 (sn,i ) + un,i ,

j=1

8

un,i

=

ρ0

n


wn,i,j un,j + εn,i ,

j=1

where xn,i ∼ i.i.d. N (0, 1), sn,i ∼ i.i.d. Uniform[0, 1], g0 (sn,i ) = sin(3πsn,i ) and εn,i ∼ i.i.d. N (0, σ02 ) for alli = 1, . . . , n. Spatial weight matrix Wn is deﬁned according to rook contiguity with row normalization (see,e.g., Arbia (2014)). As basis functions for the approximation of the nonparametric function, we use cubic
B-splines (see, e.g., Hastie et al. (2009)). Following a simple rule-of-thumb, we set the numbers of the basis
functions as n1/5 + 2 × 4, where n1/5 denotes the integer part of n1/5 .
We set β0 = 2 and σ02 = 1 as true values. As pairs of spatial autoregressive parameters (λ0 , ρ0 ), we
consider the following four cases: (λ0 , ρ0 ) ∈ {(0.2, 0.2), (0.8, 0.8), (0.2, 0.8), (0.8, 0.2)}. For each parametervalue, we generate a sample of size n (= 400, 900) and calculate the estimators. This step is repeated 1000
times. For the estimators of λ0 , ρ0 , β0 and σ02 , we report the bias and root mean squared errors (RMSE). To
evaluate the estimation performance of the nonparametric term, we use the average RMSE (ARMSE):
1000

ARM SE =1  1 l[ĝ (sn,i ) − g0 (sn,i )]2
1000
n i=1
n

1/2

,

l=1

where ĝ l (·) indicates the estimate from the l-th replicated dataset.
Table 1 summarizes the estimation results of λ0 , ρ0 , β0 , σ02 , and g0 (·). As the sample size of observations
increases, estimations become more accurate. The results demonstrate the consistency of the proposed
estimators. The ARMSE of the estimator for the nonparametric function, ĝ(·), is larger than the RMSEof the estimator for the parametric functions because the convergence rate of ĝ(·) is slower than root-N.
Moreover, the bias and RMSE of the third-step estimator, ρ̂ and σ02 , are larger than those of the 2SLS
estimator, λ̂ and β̂, respectively. With regard to the magnitude of the spatial autoregressive parameters, λ0and ρ0 , their degree does not aﬀect the estimation accuracy of the parametric terms. However, the bias and
RMSE of the estimators for the nonparametric function tend to increase as ρ0 increases.

5

Real Data Analysis

We apply the SARAR and PL-SARAR models to Boston housing price data collected by Harrison and
Rubinﬁeld (1978) to investigate the empirical properties of the PL-SARAR model and evaluate the eﬀect of
air pollution on house value. The data contain the median house prices in 506 Boston area census tracts,
NOX concentrations per town as an index of air pollution, and other potential determinants of house values.
The deﬁnitions of the variables are summarized in Table 2.

9Table 1: Small sample performances of the proposed estimators by biases and root mean square errors.

λ0
ρ0
β0
σ02
g0 (·)

Bias
RMSE
Bias
RMSE
Bias
RMSE
Bias
RMSE
ARMSE

λ0 = 0.2, ρ0 = 0.2
β0 = 1, σ 2 = 1
n = 400 n = 900-0.0001-0.00100.05110.0343-0.0321-0.0126
0.08930.0598-0.0006-0.00020.05140.0349-0.0242-0.0124
0.07130.04780.16220.1104

Dependent variable
Explanatory variables

λ0 = 0.8, ρ0 = 0.8
β0 = 1, σ 2 = 1
n =400n=900-0.0128 -0.0072
0.07790.0514-0.0202 -0.0108
0.0856
0.0563-0.0042 -0.0002
0.05120.0349-0.0111 -0.0023
0.07320.05150.57290.4171λ0 = 0.8, ρ0 = 0.2
β0 = 1, σ 2 = 1
n = 400 n = 900-0.0008-0.0006
0.03090.0197-0.0329-0.0160
0.08990.0572-0.0002-0.0009
0.05170.0337-0.0277-0.0102
0.07580.04890.17750.1219λ0 = 0.2, ρ0 = 0.8
β0 = 1, σ 2 = 1
n=400
n=900
0.0023 -0.0045
0.10370.0714-0.0326 -0.0111
0.08480.0519-0.0056 -0.0039
0.07080.0462-0.0063 -0.0047
0.0831
0.05660.51590.3716

Table 2: Variable deﬁnitions.
Variable
Deﬁnition
MEDVMedian value of owner-occupied homes.
CRIM
Per capita crime rate by town.
RM
Average number of rooms per dwelling.
AGE
Proportion of owner units built prior to 1940.
TAX
Full value property tax rate per USD 10,000 per town.
LSTAT
Proportion of lower status of the population.
INDUS
Proportion of non-retail business acres per town.
B
Black proportion of population.
DIS
Weighted distances from ﬁve Boston employment centers.
RAD
Index of accessibility to radial highways.
PTRATIO
Pupil-teacher ratio by town school district.
NOX
Nitrogen oxide concentration per town.

We compare the partially linear with the parametric linear models. Model 1 is deﬁned by:

M EDV

=

λWn M EDV + β1 + β2 CRIM + β3 RM + β4 AGE + β5 T AX + β6 LST AT+β7 IN DU S + β8 B + β9 DIS + β1 0RAD + β11 P T RAT IO + g(N OX) + un ,un=

ρWn un + εn ,where g(·) is an unknown function of NOX. We set the number of basis functions as 3 + 2 × 4 following
a simple rule-of-thumb. In model 2, we assume explanatory variable NOX is linearly correlated with thedependent variable. Therefore, we replace g(N OX) in model 1 with β12 N OX in model 2. According to Pace
and Gilley (1997) and Du et al. (2018), we deﬁne the (i, j)th element of the spatial weight matrix by:

wn,i,j = max 1 −

10

di,j
,0 ,
d0

Figure 1: Estimates of nonparametric function g(N OX) in model 1 and its 95% conﬁdence interval.
where di,j is the Euclidean distance calculated by the longitude and latitude coordinates of the two observations and d0 is the threshold distance, chosen as 0.025 in this analysis. Furthermore, we normalized theweight matrix so that the sums of rows are equal to one. The parameters in model 1 are estimated by the
proposed three-step estimation method and the ones in model 2 are estimated by 2SLS (see, e.g., Kelejian
and Prucha (1998)).
Table 3 shows the estimation results of the regression coeﬃcient, spatial autoregressive parameters, andvariances in innovation. The estimation results of models 1 and 2 are similar and the sign and statisticalsigniﬁcance of the regression coeﬃcients are consistent with previous empirical research on Boston housepricing data (see, e.g. Pace and Gilley (1997) and Arbia (2014)). Figure 1 shows the estimation results of thenonparametric function in model 1. The solid line corresponds to the estimates of g(·) and the dotted ones
to the 95% conﬁdence interval. Our empirical ﬁndings are as follows. First, a spatial correlation between
the dependent variables and disturbances exists even after we control for some of the potential determinantsof housing prices. This indicates that house values in surrounding areas have a positive eﬀect on housing
prices and there may exist unobserved shocks following a spatial pattern. Second, air pollution has a strong
negative eﬀect on housing prices in both the parametric and semiparametric models because the regression
coeﬃcient on NOX in model 2 and the estimates of g(N OX) in Figure 1 take negative values. Third, the11Table 3: Estimation results for the coeﬃcients in models 1 and 2.

Variable
CRIM
RM
INDUS
AGE
DIS
RAD
PTRARIO
B
LSTAT
TAX
NOX
Constantλρ
σ2Model 1
Coeﬃcient-0.1116
4.1387-0.0449-0.0012-0.80680.5106-0.98770.0091-0.5531-0.0215—
23.86480.57750.806225.0267

Std. error
0.03820.45220.06510.01550.32270.15780.16820.00370.05720.0056—7.74330.2847——

Model 2
Coeﬃcient-0.1025
3.8561-0.0126
0.0020-1.3219
0.2916-0.9638
0.0099-0.5362-0.0120-14.7740
26.0959
0.40370.851822.1511

Std. error
0.03270.41330.06160.01340.33750.06900.13510.00270.05100.00384.1372
7.23770.1892——eﬀect of air pollution of house prices is not linear and the negative eﬀect increases when the proportion of
NOX is over a threshold value. Figure 1 shows the proportion of NOX tends to negatively aﬀect house prices
and this negative eﬀect increases rapidly for values above 0.65. These results suggest that air pollution has
negative eﬀects on house values but that people are tolerant of air pollution to a certain extent.

6

Conclusions

In this study, we consider the PL-SARAR model and series estimation methods are employed to estimate
the nonparametric term of the proposed model. For model estimation, we propose a three-step estimation
procedure. The ﬁrst step is the estimation of the parametric regression coeﬃcient and spatial autoregressive
parameters for the dependent variables using 2SLS. The series approximation coeﬃcient for the nonparametricfunction is then estimated by OLS in the second step. The third step entails the estimation of variancesand spatial autoregressive parameters in disturbances using NLS. We then establish the consistency andasymptotic normality of the proposed estimators. Monte Carlo simulations indicate that the small sample
performances of the proposed estimator are reasonably good. Subsequently, we apply the proposed modeland estimators to Boston land price data. We ﬁnd that the proportion of NOX in the air tends to negatively
aﬀect house prices, the negative eﬀect rapidly increasing for values above 0.65.
In future studies, some extensions of this study could be considered as follows. First, GMM could be used
for the estimation of spatial autoregressive parameters in the proposed model instead of 2SLS and NLS. Lee
and Liu (2010) indicate that GMM estimators are more eﬃcient for the estimation of spatial autoregressive
12

parameters. Applying GMM estimation procedures to the proposed model improves the eﬃciency of estimation. Second, the extension of the proposed model to spatial dynamic panel data models could be considered.
Such models can control the dynamics of economic activities and unobserved time invariant heterogeneity
across spatial units. This spatial dynamic panel extension would be helpful to investigate dynamic spatial
spillover and causal eﬀects in the empirical analysis.

References[1] Anselin, L. (1988). Spatial econometrics: methods and models. Kluwer Academic Publishers, Boston.
[2] Arbia, G. (2014). A primer for spatial econometrics: with applications in R. Palgrave Macmillan, UK.
[3] Bernstein, D. S. (2009). Matrix mathematics: theory, facts, and formulas. Princeton university press,
Princeton.
[4] Cliﬀ, A.D., & Ord, J.K. (1973). Spatial Autocorrelation. Pion, London.
[5] Du, J., Sun X., Cao, R. & Zhang, Z., (2018). Statistical inference for partially linear additive spatial
autoregressive models. Spatial Statistics, 25, 52-67.
[6] Harrison Jr, D., & Rubinfeld, D. L. (1978). Hedonic housing prices and the demand for clean air. Journalof environmental economics and management, 5(1), 81-102.
[7] Hastie, T., Tibshirani, R., & Friedman, J. (2009). The elements of statistical learning. Springer Verlag,
New York.
[8] Hoshino, T. (2018). Semiparametric spatial autoregressive models with endogenous regressors: With anapplication to crime data. Journal of Business & Economic Statistics, 36(1), 160-172.
[9] Kelejian, H. H., & Prucha, I. R. (1998). A generalized spatial two-stage least squares procedure for
estimating a spatial autoregressive model with autoregressive disturbances. The Journal of Real Estate
Finance and Economics, 17(1), 99-121.
[10] Kelejian, H. H., & Prucha, I. R. (1999). A generalized moments estimator for the autoregressive parameter in a spatial model. International economic review, 40(2), 509-533.
[11] Kelejian, H. H., & Prucha, I. R. (2010). Speciﬁcation and estimation of spatial autoregressive models
with autoregressive and heteroskedastic disturbances. Journal of econometrics, 157(1), 53-67.

13

[12] Lee, L. F. (2004). Asymptotic Distributions of QuasiMaximum Likelihood Estimators for Spatial Autoregressive Models. Econometrica, 72(6), 1899-1925.
[13] Lee, L. F., & Liu, X. (2010). Eﬃcient GMM estimation of high order spatial autoregressive models withautoregressive disturbances. Econometric Theory, 26(1), 187-230.
[14] Pace, R. K., & Gilley, O. W. (1997). Using the spatial conﬁguration of the data to improve estimation.
The Journal of Real Estate Finance and Economics, 14(3), 333-340.
[15] Ptscher, B. M., & Prucha, I. R. (1997). Dynamic nonlinear econometric models, Asymptotic theory.
Springer Verlag, New York.
[16] Su, L. (2012). Semiparametric GMM estimation of spatial autoregressive models. Journal of Econometrics, 167(2), 543-560.
[17] Su, L., & Jin, S. (2010). Proﬁle quasi-maximum likelihood estimation of partially linear spatial autoregressive models. Journal of Econometrics, 157(1), 18-33.
[18] Zhang, Y., & Shen, D. (2015). Estimation of semi-parametric varying-coeﬃcient spatial panel data
models with random-eﬀects. Journal of Statistical Planning and Inference, 159, 64-80.
[19] Zhan, Y. & Sun, Y., (2015). Estimation of partially speciﬁed dynamic spatial panel data models withﬁxed-eﬀects, Regional Science and Urban Economics. 51, 37-46.

Appendix
The following facts summarize some basic properties on matrix algebras.
Fact 1.

If the row and column sums of n × n matrices C1 and C2 are uniformly bounded in absolute value,then the row and column sums of C1 C2 and C2 C1 are also uniformly bounded in absolute value (see, e.g.,
Kelejian and Prucha (1998)).
Fact 2.

Let C1 be a symmetric matrix and C2 be a positive semideﬁnite matrix. Then, γmin (C1 )tr(C2 ) ≤tr(C1 C2 ) ≤ γ(max)(C1 )tr(C2 ).
Fact 3.

For an n × n matrix C, its spectral radius is bounded by maxi

nj=1 |cn,i,j |, with cn,i,j being the(i, j)-th element of Cn (see, the appendix of Hoshino (2018)).
The following lemmas are essential for the proofs of the main results of this paper.

14

Lemma 1.

Let An be an n × n matrix whose row and column sums are uniformly bounded in absolutevalue, and Dn be a symmetric and idempotent matrix. Suppose that Assumptions 1-5 hold. Then,√
Bn An (In − Dn )Hn (In − Dn )An Bn = B̃n  An (In − Dn )Hn (In − Dn )An B̃n + Op ( n),where B̃n = (Wn (In − λ0 Wn )−1 (Xn β0 + g0 (Sn )), Xn ).
Proof. By the deﬁnition of the matrix Bn , we have

Bn=

(Wn Yn , Xn ),=(Wn (In − λ0 Wn )−1 (Xn β0 + g0 (Sn ), Xn ) + (Wn (In − λ0 Wn )−1 (In − ρ0 Mn )−1 εn , 0n×dx ),=B̃n + ε̃n .

where B̃n = (Wn (In − λ0 Wn )−1 (Xn β0 + g0 (Sn ), Xn ) and ε̃n = (Wn (In − λ0 Wn )−1 (In − ρ0 Mn )−1 εn , 0n×dx )and 0n×dx is an n × dx matrix whose components are zero.
Thus,




Bn An (In − Dn )Hn (In − Dn )Bn=(B̃n + ε̃n ) An (In − Dn )Hn (In − Dn )(B̃n + ε̃n ),=B̃n  An (In − Dn )Hn (In − Dn )An B̃n + B̃  An (In − Dn )Hn (In − Dn )An ε̃n+ε̃n (In − Dn )Hn (In − Dn )An B̃n + ε̃n An (In − Dn )Hn (In − Dn )An ε̃,=R11 + R12 + R13 + R14,

where R11 = B̃n  An (In − Dn )Hn (In − Dn )An B̃n , R12 = B̃  An (In − Dn )Hn (In − Dn )An ε̃n , R13 = ε̃n (In −Dn )Hn (In − Dn )An B̃n and R14 = ε̃n An (In − Dn )Hn (In − Dn )An ε̃.
Firstly, we consider R14. Let Tn = An Wn (In − λ0 Wn )−1 (In − ρ0 Mn )−1 . The row and column sums of
Tn is uniformly bounded in absolute value by Assumption 1 and Fact 1, and γmax (Tn Tn ) = O(1) by Fact 3.
Noting that the largest eigenvalue of an idempotent matrix is at most one, by Assumption 2 and Fact 2,E(εn Tn (In − Dn )Hn (In − Dn )Tn εn ) =11

σ 2 tr((Zn Zn ) 2 Zn (In − Dn )Tn Tn (In − Dn )Zn (Zn Zn ) 2 ,11≤

σ 2 γmax (Tn Tn )tr((Zn Zn ) 2 Zn (I − Dn )Zn (Zn Zn ) 2 ),≤σ 2 γmax (Tn Tn )tr((Zn Zn ) 2 Zn Zn (Zn Zn ) 2 ),

=

O(1).

1

15

1

Then, it follows by Markov’s inequality that R14 = Op (1).
Next, we consider R12. By assumption 5,
E||B̃n An (In − Dn )Hn (In − Dn )Tn εn ||2=

Etr(εn Tn (In − Dn )Hn (In − Dn )An B̃n B̃n An (In − Dn )Hn (In − Dn )Tn εn ),≤nσ 2 cB̃n γmax (An An )tr(Tn (In − Dn )Hn (In − Dn )Hn (In − Dn )Tn ),≤nσ 2 cB̃n γmax (An An )γmax (Tn Tn )tr(Hn ),

=

O(n).

√
√
Thus, R12 = Op ( n) by Jensen’s inequality and Markov’s inequality. Similarly, we have R13 = Op ( n).
By combining the convergence rate of R12, R13 and R14, we have

Bn An (In − Dn )Hn (In − Dn )An Bn

Lemma 2=√
R11 + Op ( n).

Let An be an n × n matrix whose row and column sums are uniformly bounded in absolutevalue, Dn be a symmetric and idempotent matrix. Suppose that Assumptions 1-5 hold. Then,
Bn An (In − Dn )Hn (In − Dn )An Vn = O(nK −rs ),Proof. By the deﬁnition of Bn , we have
Bn An (In − Dn )Hn (In − Dn )An Vn=

B̃n An (In − Dn )Hn (In − Dn )An Vn + ε̃n An (In − Dn )Hn (In − Dn )An Vn ,=R21 + R22,

where R21 = B̃n An (In − Dn )Hn (In − Dn )An Vn and R22 = ε̃n An (In − Dn )Hn (In − Dn )An Vn .
Firstly, we consider R21. By Assumption 4 and 5,||B̃n An (In − Dn )Hn (In − Dn )An Vn ||2=tr(Vn An (In − Dn )Hn (In − Dn )An B̃ B̃n An (In − Dn )Hn (In − Dn )An Vn ),≤ncB̃n γmax (An An )tr(Vn An (In − Dn )Hn (In − Dn )An Vn ),≤ncB̃n γmax (An An )γmax (An An )||Vn ||2 ,≤n2 cB̃n γmax (An An )γmax (An An ) sup |p(s) B0 − f (s)|2 ,

=

s∈S

2

O(n K

−2rs

).

Thus, R21 = O(nK −rs ) by Jensen’s inequality.

16

Next, we consider R22. Similarly, by assumption 4 and 5,
E||εn Tn (In − Dn )Hn (In − Dn )An Vn ||2=

Etr(Vn An (In − Dn )Hn (In − Dn )Tn εn εn Tn (In − Dn )Hn (In − Dn )An Vn ),≤σ 2 γmax (Tn Tn )γmax (An An )||Vn ||2 ,≤σ 2 γmax (Tn Tn )γmax (An An )n sup |p(s) B0 − f (s)|2 ,=O(nK −2rs ).

s∈S√Thus, R22 = Op ( nK −rs ) by Jensen’s inequality and Markov’s inequality.
By combining the convergence rate of R21 and R22, we have
Bn An (In − Dn )Hn (In − Dn )An Vn = O(nK −rs ).

Proof of Theorem 1δ̂

By the deﬁnition of δ̂,= (Bn (In − Πn )Hn (In − Πn )Bn )−1 Bn (In − Πn )Hn (In − Πn )Yn ,=δ0 + (Bn (In − Πn )Hn (In − Πn )Bn )−1 Bn (In − Πn )Hn (In − Πn )V+(Bn (In − Πn )Hn (In − Πn )Bn )−1 Bn (In − Πn )Hn (In − Πn )(I − ρ0 Mn )−1 εn .

Thus,
√

−1

1 
B (In − Πn )Hn (In − Πn )Bnn nn(δ̂ − δ0 ) =+1 
B (In − Πn )Hn (In − Πn )Bn
n n1√ Bn (In − Πn )Hn (In − Πn )Vn−1

1√ Bn (In − Πn )Hn (In − Πn )(In − ρ0 Mn )−1 εn .
n

By Lemma 1 and 2,
1 
B (In − Πn )Hn (In − Πn )Bn
n n1√ Bn (In − Πn )Hn (In − Πn )Vnp−→ Σ2 ,p−→ 0.

By Slutsky’s theorem and a central limit theorem, we have√n(δ̂ − δ)=R11+ O(n−1 )n−11√ Bn (In − Πn )Hn (In − Πn )(In − ρ0 Mn )−1 εn + O(K −rs ) ,n−1−→ N (0, σ 2 Σ−1
2 Σ1 Σ2 ).
d17

Proof of Theorem 2α̂Firstly, we consider the convergence rate of α̂. By the deﬁnition of α̂,=(Pn Pn )−1 Pn (Yn − Bn δ̂),=

α0 + (Pn Pn )−1 Pn Bn (δ0 − δ̂) + (Pn Pn )−1 Pn Vn + (Pn Pn )−1 Pn (In − ρ0 Mn )−1 εn ,=α0 + R31 + R32 + R33,

where R31 = (Pn Pn )−1 Pn Bn (δ0 − δ̂), R32 = (Pn Pn )−1 Pn Vn and R33 = (Pn Pn )−1 Pn (In − ρ0 Mn )−1 εn .
By the deﬁnition of Bn , we have

R31 =
=

(Pn Pn )−1 Pn B̃(δ0 − δ̂) + (Pn Pn )−1 Pn ε̃(δ0 − δ̂),
R41 + R42,

where R41 = (Pn Π)−1 Π B̃(δ0 − δ̂) and R42 = (Pn Pn )−1 Pn ε̃(δ0 − δ̂).
Firstly we consider R41.
||(Pn Pn )−1 Pn B̃n (δ0 − δ̂)||2=

tr((δ0 − δ̂) B̃n Pn (Pn Pn )−2 Pn B̃n (δ0 − δ̂)),≤1 −1c tr((δ0 − δ̂) B̃n P (P  P )−1 P  B̃n (δ0 − δ̂)),
n Pn

≤


c−1
Pn cB̃n tr((δ0 − δ̂) (δ0 − δ̂)),≤cΠ cB̃ tr((δ0 − δ̂) (δ0 − δ̂)),=O(n−1 ).

Thus, R41 = O(n−1/2 ) by Jensen’s inequality.
Similarly, we consider R42.
E||(Pn Pn )−1 Pn Tn εn (λ0 − λ̂)||2= (λ0 − λ̂)2 σ 2 tr(Tn Pn (Pn Pn )−2 P  Tn ),1tr(Tn Tn ),= (λ0 − λ̂)2 σ 2 c−1
n Pn=

O(n−1 )

Thus, R42 = Op (n−1/2 ) by Jensen’s inequality and Markov’s inequality.
Therefore, we have R31 = Op (n−1/2 ) by combining the convergence rate of R41 and R42.

18

Next, we consider R32.
||(Pn Pn )−1 Pn Vn ||2=tr(Vn Pn (Pn Pn )−2 Pn Vn ),≤1 −1c tr(V  V ),
n Pn

≤


2
c−1
Pn sup |p(s) B0 − f (s)| ,

=

O(K

s∈S
−2rs

).

Thus, R32 = O(K −rs ) by Jensen’s inequality.
Finally, we consider R33.
E||(Pn Pn )−1 Pn (In − ρ0 Mn )−1 εn ||2

Etr(εn (In − ρ0 Mn ) −1 Pn (Pn Pn )−2 Pn (In − ρ0 Mn )−1 εn ,=
1 2 −1σ cPn γmax ((In − ρ0 Mn ) −1 (In − ρ0 Mn )−1 )tr(Pn (Pn Pn )−1 Pn ),
n
K
.
O
n

≤
=

√ √
Thus, R33 = Op ( K/ n)by Jensen’s inequality and Markov’s inequality.
√√K + K −rsnTherefore, we obtain α̂ = α0 + Op.

Next, we consider the uniform convergence rate of ĝ(·). By the triangle inequity and Cauchy-Schwarz
inequality, we havesup ||ĝ(s) − g0 (s)|| ≤sup ||P K (s)(α̂ − α0 )|| + sup ||P K (s)α0 − g0 (s)||,

s

s

s

||α̂ − α0 || sup ||P (s)|| + O(K −rs ),

≤

K

s

=

Op

K
√ + K (1−2rs )/2 .
nFinally, we consider the limiting distribution of ĝ(·). By the deﬁntion of ĝ(s),ĝ(s) − g0 (s) ==P K (s)α̂ − (pK α0 + O(K −rs ),pK (R31 + R32 + R33) + O(K −rs ).

It follows by the above discussion that||pK R31||=||P K (s)(Pn Pn )−1 Pn B(δ − δ̂)||,≤ ||P K (s)|| ||(Pn Pn )−1 Pn B(δ − δ̂)||,

19

=

√
K
O √
,
n

and

||pK R32||=||P K (s)(Pn Pn )−1 Pn V ||,≤||P K (s)|| ||(Pn Pn )−1 Pn Vn ||,=O(K (1−2rs )/2 ).

Thus,ĝ(s) − g0 (s) =P K (s)(Pn Pn )−1 Pn (In − ρ0 Mn )−1 εn + O√
K√ + K (1−2rs )/2 .
n

Let us consider the variance of the ﬁrst term of the above equation.
σ 2 (s) =E(P K (s)(Pn Pn )−1 Pn (In − ρ0 Mn )−1 εn εn (In − ρ0 Mn ) −1 Pn (Pn Pn )−1 P K (s)),
=σ 2 (P K (s)(Pn Pn )−1 Pn (In − ρ0 Mn )−1 (In − ρ0 Mn ) −1 Pn (Pn Pn )−1 P K (s)),≤σ 2 γmax ((In − ρ0 Mn )−1 (In − ρ0 Mn ) −1 )(P K (s)(Pn Pn )−1 P K (s)),

≤
=







1
σ 2 γmax ((In − ρ0 Mn )−1 (In − ρ0 Mn ) −1 )cPn (P K (s)P K (s)),
n
K
.
O
n

Similarly, σ 2 (s) ≥ O(K/n), Thus σ 2 (s) = O(K/n).
By Slutsky’s theorem and a central limit theorem, we obtaineĝ(s) − g0 (s) −→ N (0, σ 2 (s)).
dProof of Theorem 3

Let ûn = Yn − Bn δ̂ − Pn α̂. As the ﬁrst step, we show that11 û An ûn − E un An un = op (1),n nnwhere An is a matrix whose row and column sums are uniformly bounded in absolute values.

20Note that11 ûn An ûn − E ûn An ûn
n
n

1 
1
ûn An ûn − un An un ,
n
n
1 
1
+ un An un − E un An un .
nn=Firstly, we show that n1 un An un − E n1 un An un = op (1). By the deﬁnition of un
1 u An un
n n

=
=

1 
ε (In − ρ0 Mn ) An (In − ρ0 Mn )εn ,
n n
1  ∗
ε A εn ,n n n

where A∗n = (In − ρ0 Mn ) An (In − ρ0 Mn ) and the row and column sums of A∗n are uniformly bounded inabsolute value by Fact 1. Thus it follows that n1 un An un − E n1 un An un = op (1) immediately from the basic
property of laws of large numbers in Lee (2004).
Next, we consider that n1 ûn An ûn − n1 un An un = op (1). By the deﬁnition of ûn ,ûn=Yn − Bn δ̂ − Pn α̂,=Yn − λ̂Wn Yn − Xn β̂ − Pn α̂,=un + (λ0 − λ̂)Wn Yn + Xn (β0 − β̂) + (g0 (Sn ) − Pn α̂),=un + (λ0 − λ̂)Wn (In − λ0 Wn )−1 (In − ρ0 Mn )εn+(λ0 − λ̂)Wn (In − λ0 Wn )−1 (Xn β0 + g0 (Sn )) + Xn (β0 − β̂) + (g0 (Sn ) − Pn α̂),=un + ψ1 + ψ2 + ψ3 + ψ4 ,where ψ1 = (λ0 − λ̂)Wn (In − λ0 Wn )−1 (In − ρ0 Mn )εn , ψ2 = (λ0 − λ̂)Wn (In − λ0 Wn )−1 (Xn β0 + g0 (Sn )), ψ3 =Xn (β0 − β̂) and ψ4 = (g0 (Sn ) − Pn α̂).
Thus,
1
1 
ûn An ûn − un An unnn=φ1 + φ2 + φ3 + φ4 + 2φ5 + 2φ6 + 2φ7 + 2φ8 + 2φ9 + 2φ10
+2φ11 + 2φ12 + 2φ13 + 2φ14 ,where φ1 = n1 ψ1 An ψ1 , φ2 = n1 ψ2 An ψ2 , φ3 = n1 ψ3 An ψ3 , φ4 = n1 ψ4 An ψ4 , φ5 = n1 un An ψ1 , φ6 = n1 un An ψ2 , φ7 =1 n u n An ψ 3 , φ 8=

1 n u n An ψ 4 , φ 9=1 n ψ1 An ψ2 , φ10=1 n ψ1 An ψ3 , φ11

21

=

1 
n ψ1 An ψ4 , φ12=

1 n ψ2 An ψ3 , φ13

=

1 
1 
n ψ2 An ψ4 and φ14 = n ψ3 An ψ4 . We show that φi , i = 1, . . . , 14, are or order op (1). Here, note thatρ̂ − ρ0=β̂ − β0=g0 (s) − Pn α̂

=
=
=

1
Op √ ,
n1
Op √ ,nK
Op √ + K (1−2rs )/2 ,n√K
K√nK −rs ,Op √ + √nn
op (1).

For example,

Eφ1

=
=
=
=


1
E (λ0 − λ̂)2 εn (In − ρ0 Mn ) (In − λ0 Wn ) −1 Wn An (λ0 − λ̂)Wn (In − λ0 Wn )−1 (In − ρ0 Mn )εn ,n1
E(λ0 − λ̂)2 εn Ãn εn ,
n
n

2
21
ãn,i,j ,σ0 (λ0 − λ̂)
n i=1

op (1),


where Ãn = (In − ρ0 Mn ) (In − λ0 Wn ) −1 Wn An (λ0 − λ̂)Wn (In − λ0 Wn )−1 (In − ρ0 Mn ) and ãn,i,j is the (i, j)thelement of the matrix Ãn . The remaining terms can be shown to be op (1) in the same way. Therefore,1 
1 n ûn An ûn − E n un An un = op (1)We prove the consistency of the third step estimator following Kelejian and Prucha (1999). The objectivefunction of the nonlinear least squares estimator and its corresponding counterpart are given by

Rn (θ)=[Gn − gn ] [Gn − gn ],R̂n (θ)=[Ĝn − ĝn ] [Ĝn − ĝn ],where θ = (ρ, σ 2 ) .
Let θ0 = (ρ0 , σ02 ) . By Assumption 7,
Rn (θ) − Rn (θ0 )=[ρ − ρ0 , ρ2 − ρ20 , σ 2 − σ02 ]GN Gn [ρ − ρ0 , ρ2 − ρ20 , σ 2 − σ02 ] ,≥cGn [ρ − ρ0 , σ 2 − σ02 ][ρ − ρ0 , σ 2 − σ02 ] ,=cGn ||θ − θ0 ||2 .

22

It follow that for every ε > 0 and any N ,inf

θ:||θ−θ0 ||≥ε[Rn (θ) − Rn (θ0 )]≥c Gn ε2 ,>0.

Thus, the identiﬁability of θ is proved.
Let Fn = [Gn , −gn ], F̂n = [Ĝn , −ĝn ], ρ ∈ [−a, a] and σ 2 ∈ [0, b].
|Rn (θ) − R̂n (θ)|

=





[ρ, ρ2 , σ 2 , 1][Fn Fn − F̂n F̂n ][ρ, ρ2 , σ 2 , 1] ,≤ ||Fn Fn − F̂n F̂n || [1 + a2 + a4 + b2 ].
The elements of Fn and F̂n are all of the form n1 ûn An ûn and E n1 ûn An ûn where the row and column sumsof An are uniformly bounded in absolute value. We have shown that n1 ûn An ûn − E n1 ûn An ûn = op (1). Thus,
Fn − F̂n = op (1). It follow thatsup |Rn (θ) − R̂n |ρ,σ≤||Fn Fn − F̂n F̂n ||[1 + a2 + a4 + b2 ],

p

−→ 0.
The consistency of ρ̂n and σ̂n2 follows form Lemma 3.1 in Potscher and Prucha (1997).

23