Discussion Paper No. 78
Spatial GARCH Models
Takaki Sato
Yasumasa Matsuda

March , 2018

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Spatial GARCH Models
Takaki Sato ∗Yasumasa Matsuda †Abstract
This study proposes a spatial extension of time series generalized autoregressive conditional heteroscedasticity (GARCH) models. We call the spatial extended GARCH models as spatial GARCH(S-GARCH) models. S-GARCH models specify conditional variances given simultaneous observations,which constitutes a good contrast with time series GARCH models that specify conditional variances given
past observations. The S-GARCH model are transformed into a spatial autoregressive moving-average(SARMA) model and the parameters of the S-GARCH model are estimated by a two step procedure.
First step estimation is the quasi maximum likelihood (QML) estimation method and consistency andasymptotic normality of the proposed QML estimators are given. Second step is estimation of an intercept term by the estimator derived from another QML to avoid bias in first step and consistency of theestimator is shown. We demonstrate empirical properties of the model by simulation studies and real
data analyses of land price data in Tokyo areas. We find the estimators have small bias regardless of
distributions of error terms from simulation studies and real data analyses show that spatial volatility
in land price has global spillover and volatility clustering, namely units with higher spatial volatility are
clustered in some specific districts like time series financial data.

Keywords: GARCH model, Spatial ARMA model, Quasi maximum likelihood, areal data, spatial volatility.

1

Introduction

Volatility models for time series financial data have developed with their application in academia and thefinancial industry. The seminal work by Engle (1982) introduces the autoregressive conditional heteroscedasticity (ARCH) model and Bollerslev (1986) proposes a extension known as the generalized ARCH (GARCH)
model. These models are widely used to model and forecast volatility of univariate time series data for
calculation of the price of options or value at risk of a financial position in risk management. Subsequently,
Multivariate extensions of univariate models are proposed by Bollerslev et al (1988), Bollerslev (1990) andEngle and Kroner (1995) for modeling dynamic relationships between volatility of multiple asset returns.
A major challenge of multivariate volatility modeling is to overcome the curse of dimensionality; there aren(n + 1)/2 variances and covariances for n-dimensional asset return series. One solution for the problem isconsider simpler structures of covariance matrices to reduce parameters.
The ideas of spatial econometrics have been applied to volatility models in recent years. Two main
objectives of the applications are to reduce parameters in covariance matrices and to extend time series
volatility models to spatial models for spatial data. Caporin and Paruolo (2008) and Borovkova and Lopuhaa(2012) have applied the ideas of spatial econometrics to time series multivariate GARCH models from theformer view point. On the other hand, Yan (2007) and Robinson (2009) have done a spatial extension ofstochastic volatility models which are another kind of volatility models and Sato and Matsuda (2017) haveextend time series ARCH models to spatial ARCH (S-ARCH) models from both view points.
This paper aims to extend S-ARCH models to spatial generalized ARCH (S-GARCH) models. The SGARCH model have two interesting features. Firstly, volatility at a point or an area in map is specified
by surrounding observations in the S-ARCH model, whereas that of the S-GARCH model is characterized
by surrounding observations and volatility. Thus, the S-GARCH model captures global spatial spillover in∗ Graduate
Schoolof
Economics
andManagement,
TohokuUniversity,
Sendai
980-8576,
Japan.
takaki.sato.s4@dc.tohoku.ac.jp† Graduate School of Economics and Management, Tohoku University, Sendai 980-8576, Japan. matsuda@econ.tohoku.ac.jp

1volatility in spatial data. Secondly, the S-GARCH model can be transformed into a spatial autoregressive
moving average (SARMA) model. This means the existence condition of the S-GARCH model are easily
established.
Parameters in the S-GARCH model are estimated by the quasi-maximum likelihood (QML) estimation
method and we show the QMLE estimators have consistency and asymptotic normality. Two estimation
methods are basically used in spatial econometrics literature. First one is the moment method. Kelejian andRobinson (1993) and Kelejian and Prucha (1997, 1998) propose two stage least squares estimation methods
and Lee (2007) propose the generalized method of moments (GMM) for spatial autoregressive (SAR) models
and spatial autoregressive models which also have a spatial autoregressive process in disturbances (SARAR).
Moreover, Dogan and Taspinar (2013) consider the GMM methodology for spatial autoregressive models
with moving average disturbances (SARMA). Second one is the QML estimation method. Lee (2004), Yu et
al (2008), Su and Yang (2015) propose the QML estimation method for SAR models and spatial dynamic
panel (SDP) models and also Yang (2015) shows M-estimator based on the QML for SDP models which havespatial autoregressive process in both dependent variables and disturbances. However, asymptotic propertiesof the QML estimator for SARMA models has not been discussed. As mentioned above, S-GARCH models
can be transformed into SARMA models. Therefore, we show asymptotic properties of the QML estimatorfor SARMA models.
This paper proceeds as follows. Section 2 introduces the S-GARCH model and discusses characteristicsof the model. Estimation methods for the model and asymptotic properties of the estimators are derived in
Section 3. Section 4 examines empirical properties of the model by applying to simulated and land price datain Tokyo area. Section 5 concludes the paper. All the proofs are collected in the Appendix.

2

Model specification

We consider the S-GARCH model of the form√hi εi ,
yi =
n
n
∑
∑
log hi = λwi,j log hj + ρ
wi,j log yj2 + α + zi′ δ,
j=1

i

=

j=1

1, . . . , n,√where yi is an areal data, hi is volatility, εi is an independent and identically distributed (i.i.d) randomvariable with zero mean and variance 1, zi is (k × 1) non-stochastic regressors, and wi,j is a spatial weight
which is predetermined and quantifies a closeness from area i to area j. Parameters in this model are(λ, ρ, α, δ ′ )′ . Scalar parameters λ and ρ characterizes the simultaneous eﬀect, α is an intercept term and δ isusual regression coeﬃcients.
The S-GARCH model is diﬀerent from the time series GARCH model is the following two points. First
one is a description of volatility. Spatial volatility in the S-GARCH model is described by observationsand volatility at all other units, on the other hands time series volatility is defined by past observationsand volatility following the flow of time. Although the descriptions of time series and spatial volatility are
diﬀerent, we have found in this paper that they have the similar properties. For instance, volatility clustering
exists, namely large changes tend to be followed by large changes and small changes tend to be followed bysmall changes. This is a stylized feature of financial time series data and land price data also has the similar
property that a large change at one area leads to large changes at surrounding areas.
Second one is the log transformation of volatility. Log transformation is used to ensure the existence of
areal data yi . If we define non logarithmic volatility, it would be diﬃcult to guarantee the existence condition
unlike that of time series GARCH models that can be derived from Markov process theories (Fan and Yao(2003)). On the other hand, the log transformation of volatility makes it much easier to show the existence
condition because the S-GARCH model can be transformed into the spatial autoregressive moving average(SARMA) model as shown below and the existence condition of the SARMA model is already known.
Let us introduce the following SARMA transformation of the S-GARCH model. Denoting log y 2 = (log y12 ,. . . , log yn2 )′ , log h = (log h1 , . . . , log hn )′ , log ε2 = (log ε21 , . . . , log ε2n )′ , Zn = (z1 , . . . , zn )′ , 1n = (1, . . . , 1)′ andIn is a n × n identity matrix, the model has the following vector form representation,2

log y 2=

log h =

log h + log ε2(1)2

λWn log h + ρWn log y + α1n + Zn δ,(2)where Wn is a spatial weight matrix whose elements are wi,j . From (2),log h = (In − λWn )−1 (ρWn log y 2 + α1n + Zn δ),
By substituting (2) into (1),log y 2(In − λWn ) log y 2
log y 2= (In − λWn )−1 (ρWn log y 2 + α1n + Zn δ) + log ε2 ,= ρWn log y 2 + α1n + Zn δ + (In − λWn ) log ε2 ,= λWn log y 2 + ρWn log y 2 + α1n + Zn δ + (In − λWn ) log ε2 .

(3)This is the SARMA model and the existence condition holds when |λ| + |ρ| < 1.

3

Estimation

We consider the estimation of the parameters (λ, ρ, α, δ ′ )′ and the asymptotic properties of the estimators.
Parameters are estimated by a two step procedure. First step is the estimation of (λ, ρ, δ ′ )′ by the QML
estimation method. The proposed QML estimator are consistent and asymptotically normal. However, log ε2
in (3) is not zero mean error terms. Thus, the estimator for α in the first step has bias, therefore we needto estimate α by another method. In second step, α is estimated with consistent estimator derived from the
QML based on the likelihood diﬀerent from the one in the first step.

3.1

First step estimation

Parameters λ, ρ and δ are estimated in first step by the QML estimation method.
To apply the QML estimation method, we need to modify the error term because the mean of log ε2 in(3) is not zero as already mentioned. From (3),α1n + (In − λWn ) log ε2= α1n + (In − λWn ){log ε2 − E(log ε21 )1n + E(log ε21 )1n },= {α + (1 − λ)E(log ε21 )}1n + (In − λWn ){log ε2 − E(log ε21 )1n }.

Noting that intercept term has a bias by (1 − λ)E(log ε21 ).
Denoting Yn = log y 2 , Xn = [1n , Z], Vn = {log ε2 − E(log ε21 )1n } and β = ({α + (1 − λ)E(log ε21 )}, δ ′ )′ ,the model has the following representation,
Yn = λWn Yn + ρWn Yn + Xn β + (In − λWn )Vn ,(4)

where Vn is already zero mean processes.
Now, let us consider the QML estimation of (4). Regarding vi′ s as independent Gaussian variables withmean zero and variance σ 2 , the likelihood function of (4) is
log Ln (ψ) = −nV ′ (θ, β)Vn (θ, β)− log |Rn (λ)| + log |Sn (θ)|,
log(2πσ 2 ) − n
2
2σ 2

(5)

where θ = (λ, ρ)′ , ψ = (β ′ , σ2 , θ′ )′ , Rn (λ) = In − λWn , Rn = In − λ0 Wn , Sn (θ) = In − λWn − ρWn ,Sn = In − λ0 Wn − ρ0 Wn and Vn (θ, β) = Rn−1 (λ)[Sn (θ)Yn − Xn β]. The QML estimator is the extreme
estimator derived form the maximization of (5).
It is convenient to work with the concentrated likelihood by concentrating β and σ 2 out for computation
and asymptotic analysis on the estimator. From the first order condition of (5), the concentrated QML
estimators of β and σ 2 isβ̂n (θ) =(Xn′ Rn′−1 (λ)Rn−1 Xn )−1 Xn′ Rn′−1 (λ)Rn−1 (λ)Sn (θ)Yn ,3

σ̂n2 (θ) =V̂n′ (θ)V̂n (θ),
nwhere V̂n (θ) = Rn−1 (λ)[Sn (θ)Yn − Xn β̂n (θ)]. The concentrated likelihood function of θ islog Ln (θ) =nn− (log(2π) + 1) − log σ̂n2 (θ) − log |Rn (λ)| + log |Sn (θ)|.
2
2

(6)

The QML estimator θ̂n maximizes the concentrated likelihood function (6) and the QML estimators of β andσ 2 are β̂n (θ̂n ) and σ̂n2 (θ̂n ), respectively.
For our analysis of the asymptotic properties of first step estimators, we need the following assumptions:
Assumption 1. The disturbances {vi }, i = 1, . . . n are i.i.d. across i with zero mean, variance σ02 andE|vi |4+η < ∞ for some η > 0.
Assumption 2. The elements wn,ij of Wn are nonnegative and row normalized and the column sums of Wn
are uniformly bounded.
Assumption 3. The space Θ is compact, and the true parameter θ0 lies in its interior.
Assumption 4. The matrix Sn , Sn (θ), Rn , and Rn (λ) are uniformly bounded both row and column sums
and nonsingular.
Assumption 5. The elements of Xn are uniformly bounded constants. The limn→∞ n1 (Xn′ Rn′−1 (λ)Rn−1 (λ)Xn )exists and is nonsingular.
Assumption 6. 0 ≤ cy ≤ inf θ∈Θ γmin (V ar(Sn (θ)Yn )) ≤ supθ∈Θ γmax (V ar(Sn (θ)Yn )) ≤ cy < ∞.
Assumption 7. 0 ≤ cr ≤ inf λ∈Λ γmin (Rn′−1 (λ)Rn−1 (λ)) ≤ supλ∈Λ γmax (Rn′−1 (λ)Rn−1 (λ)) ≤ cr < ∞Assumption 8. limn→∞ n1 β0′ Xn′ Sn′−1 Sn′ (θ)Rn′−1 (λ)Mn Rn−1 (λ)Sn (θ)Sn−1 Xn β0 ̸= 0,
where Mn = In − Rn−1 Xn (Xn′ Rn′−1 (λ)Rn−1 (λ)Xn )−1 Xn′ Rn′−1 .
To derive the consistency of the QML estimators, we need to show the identification of θ0 . Define
Qn (θ) = maxβ,σ2 E(log Ln (ψ)). The optimal solutions of this maximization problem are given byβn∗ (θ) =(Xn′ Rn′−1 (λ)Rn−1 Xn )−1 Xn′ Rn′−1 (λ)Rn−1 (λ)Sn (θ)E(Yn ),1=
E(Vn′∗ (θ)Vn∗ (θ)),
nσn2∗where Vn∗ (θ) = Rn−1 (λ)[Sn (θ)Yn − Xn βn∗ (θ)]. Therefore,
Qn (θ) =nn− (log(2π) + 1) − log σn2∗ (θ) − log |Rn (λ)| + log |Sn (θ)|,
22and identification of θ0 is based on n1 Qn (θ).
Consistency of the QML estimators θ̂ follows from the uniform convergence of n1 log Ln (θ) − n1 Qn (θ) tozero on Θ and identification of θ0 .
Theorem 1. Under Assumptions 1-8, θ0 is globally identifiable and θn is a consistent estimator of θ0 .
To derive the asymptotic distribution of the QMLE ψ̂n , we need to make the Taylor∂
log
Ln (ψ̂) = 0 at ψ0 . The first-order derivatives of the log-likelihood function at ψ0 are∂ψ
1 ∂ log Ln (ψ0 )√∂βn1 ∂ log Ln (ψ0 )√∂σ 2n
1 ∂ log Ln (ψ0 )
√
∂ρ
n

=
=
=

expansion of1√ X ′ R′−1 Vn ,σ02 n n n
1
√ (V ′ Vn − nσ02 ),2σ04 n n
11√ β ′ X ′ S ′−1 Wn′ Rn′−1 Vn + 2 √ (Vn′ Rn′ Sn′−1 Wn′ Rn′−1 Vn − σ02 tr(Sn−1 Wn )),σ02 n 0 n n
σ0 n
4

1 ∂ log Ln (ψ0 )
√
∂λ
n

=

1
√ β ′ X ′ S ′−1 Wn′ Rn′−1 Vn +σ02 n 0 n n
1√ (V ′ (R′ S ′−1 Wn′ Rn′−1 − Wn′ Rn′−1 )Vn − σ02 tr(Sn−1 Wn ) + σ02 tr(Rn−1 Wn )),σ02 n n n n

where tr(·) denote the trace of a matrix.
These involve linear and quadratic function of Vn . The asymptotic distribution of these score functions
are derived from the central limit theorems for linear-quadratic forms in Kelejian and Prucha (2001).
∂
The variance matrix of √1n ∂ψ
log Ln (ψ0 ) is())(
1 ∂ log Ln (ψ0 ) 1 ∂ log Ln (ψ0 )1 ∂ 2 log Ln (ψ0 )
√
E √
=
−E
+ Ωψ,n ,
∂ψ
∂ψ ′
n
∂ψψ ′
n
n
)
( ∂2
is the average Hessian matrix and Ωψ,n is a symmetric matrix and both are
where −E n1 ∂ψψ′ log Ln (ψ0 )given in Appendix A. When Vn is normally distributed, Ωψ,n = 0.
The score function and Hessian matrix have proper asymptotic behavior, therefore we have the following
theorem.
Theorem 2. Under Assumptions 1-8,
) d
(
)
√ (
−1
−1
n ψ̂n − ψ0 −→ N 0, Σ−1ψ + Σ ψ Ωψ Σ ψ ,( ∂2)where Σψ = − limn→∞ E n1 ∂ψψand Ωψ = limn→∞ Ωψ,n . Σψ and Ωψ assume to exist and −Σψ′ log Ln (ψ0 )
) d
(
)
√ (
to be positive definite, suﬃciently large n. When errors are normally distributed, n ψ̂n −ψ0 −→ N 0, Σ−1
ψ .

3.2

Second step estimation

Let us consider the estimation of α in the second step. As we mentioned, β1 = {α + (1 − λ)E(log ε21 )} in thefirst step has the bias. Therefore, we need to estimate α separately.
We regard εi in (3) as independent Gaussian variables, not log ε2i . Then, log ε2i follows a log chi-squared
distribution on 1 degress of freedom (Lee 2012, p379). The probability density function of log ε2i is given by
)
(
1
1
1
f (x) = √ exp − exp(x) + x .
(7)
2
2
2π
For notational purposes, we define ϕ = (λ, ρ, α, δ ′ )′ , Yn = log y 2 and U = log ε2 . Then, from (3),Rn−1 (λ)(S(θ)Y − α1n − Zn δ),α= Rn−1 (λ)(S(θ)Y − Zn δ) −
1n ,
1−λ
α
= C−
1n ,
1−λ

Un (ϕ) =where C = Rn−1 (λ)(S(θ)Y − Zn δ).
Therefore, the likelihood function based on the density (7) is
(
)}
n {
∑
1
n
α
1 exp(Ci )+log Ln (ϕ) = log 2π −Ci −− log |Rn (λ)| + log |Sn (θ)|,
−
α
) 2
2
2 exp( 1−λ1−λi=1where Ci is the i-th element of C.
Diﬀerentiating it with respect to α, the concentrated QML estimator of α given λ, ρ and δ is)n1∑
exp(Ci )αn (λ, ρ, δ) = (1 − λ) log
n i=1
(

5

Finally, substituting the proposed QML estimator (λ̂, ρ̂, δ̂ ′ ) in the first step for (λ, ρ, δ ′ )′ , we propose
( ∑
)
n
[{ −1
}]
1
α̂n = (1 − λ̂) log
exp Rn (λ̂)(S(θ̂)Yn − Zn δ̂) i ,
(8)
n i=1
as an estimator for α.
The estimator α̂n has consistency.
Theorem 3. Under Assumptions 1-8, α̂n is a consistent estimator of α0 .

4Empirical analysis

We examine the empirical properties of the S-GARCH model by applying to simulated and land price datain Tokyo areas. Monte Carlo experiments are carried out to investigate the finite sample performance of the
proposed estimators.

4.1

Simulation studies

To investigate finite sample properties of the proposed estimators, we use the following data generatingprocess:√yi =hi εi ,
nn∑∑
log hi = λwi,j log hj + ρ
wi,j log yj2 + α + xi β,
j=1

j=1

where xi ’s are randomly generated form independent normal distributions and the spatial weights matrix is
generated according to Rook contiguity and row normalizing. The error, εi , distributions are (i) standard
normal distributions, (ii) chi-squared distributions with 3 degrees of freedom and (iii) log normal distributions. Let ϕ = (λ, ρ, α, β)′ . We choose ϕ10 = (0.9, 0.05, 0.5, 1)′ , ϕ20 = (0.45, 0.45, 0.5, 1)′ , ϕ30 = (0.05, 0.9, 0.5, 1)′
and n = 100 or n = 400. Each set of Monte Carlo results is based on 1000 samples and the parameters are
estimated by the two step procedure.
The empirical means and square root of mean squared errors (RMSE) for the proposed estimators arereported in Table 1. The results show the estimators in the firs step, (λ̂, ρ̂, β̂ ′ )′ are nearly unbiased and notsensitive to the choice of the error distributions. On the other hand, the second step estimator, α̂ depends
on true parameters and the error distribution. Small λ may be attributed to the poor performance of theestimator because 1 − λ̂ in (8) eﬀects on estimated errors from true value as shown in the proof of Theorem
3. Moreover, as the error distribution is more discrepant from the Gaussian distribution, the estimator has
bigger bias and less eﬃciency. However, the empirical performances of the estimator improve as n becomeslarger.

4.2

Land price data analysis

We apply the S-GARCH model to land price data in Tokyo area
Let us introduce land price data used in this section. We use prefectural land price research as land
price data. The Japanese Ministry of Land, Infrastructure, Transport, and Tourism publishes land prices on
sampling points scattered irregularly all over Japan in the form of price per m2 in July. We focus on the land
prices over Tokyo area (Tokyo, Kanagawa, Saitama, Chiba, Tochigi, Ibaraki, Gunma) observed form 2009 to
2014. The log returns of the land prices are averaged in municipal units. Therefore, our data set consist of
344 discrete unit’s average log returns from 2010 to 2014.
Before application of the S-GARCH model, we remove spatial correlations in data with the spatial autoregressive (SAR) model year by year. This modification is similar to that we apply the ARMA model todata before fitting the GARCH model to remove correlation in time series analysis. The SAR model isyi = ζ + κ

344∑

wi,j yj + ui,t , ui,t ∼ i.i.d(0, τ 2 ).

j=1

6

where W = (wi,j ) is given the first-order contiguity relation that takes 1 when two units have a common
boarder.
We apply the S-GARCH model to the residuals obtained after fitting the SAR model year by year, where
the same spatial weight matrix as one in the SAR model was employed. Explanatory variables are intercept
term and each unit’s area. Areas of observations are included to hold Assumption 8 which is important for
the identification uniqueness. Table 2 shows the estimated values of λ, ρ, α and β. Here, α and β is intercept
and the coeﬃcient of logarithm of areas, respectively. The standard errors of λ̂ and ρ̂ are derived in Theorem
2 by replacing the population moments with the corresponding sample moments. Figure 1 express the spatial
volatility evaluated bylog ĥ = (In − λ̂Wn )−1 (ρ̂Wn log y 2 + α̂1n + xβ̂),
where x is the vector of the areas of observations.
we find estimates of spatial correlation of volatility, λ, are significant after the Great East Japan Earthquake in 2011 until 2013 from Table 2. This may show that volatility in land prices have strong correlation
when a big event occurs. The eﬀects from simultaneous returns, ρ, are not large and this is similar to empirical results of the time series GARCH model. The sum λ̂ + ρ̂ takes near 1 values between 2011 and 2013.
Thus, volatility is persistent to far areas and may generate volatility clustering. From Figure 1, not only thevolatility of costal area which hit by the Tsunami but also that of near Fukushima areas is high. This may becaused as the result of Fukushima nuclear accident. Moreover, we find the volatility clustering as explained
above. Therefore, volatility in land price takes similar behavior to that of time series financial data like stock
price. In addition, we can find volatility in land price has global spillover from figure 2. The model fittingof the S-GARCH model to the data is better than that of the S-ARCH model. The estimated volatility of
the S-ARCH model makes small clusters. On the other hand, that of the S-GARCH model generates large
clusters. This result shows that the estimated volatility is globally strongly spatially correlated.

5

Conclusion

We have proposed a spatial generalized autoregressive conditional heteroskedasticity (S-GARCH) model as
extension of a spatial autoregressive conditional heteroskedasticity (S-ARCH) model (Sato and Matsuda(2017)) to evaluate spatial volatility. The S-GARCH models can be expressed in the form of a spatial
autoregressive moving average (SARMA) model and we propose the two step estimation procedure to estimate
the parameters in the model. The quasi maximum likelihood (QML) estimators in each step have desired
asymptotic properties. Finite sample performances of the estimators are reasonably good from Monte Carloexperiments. We find volatility in land prices is similar behavior to that of time series data from real data
analysis.
For future research, we describe possible extensions. We used the first-order contiguity relations to makethe spatial weight matrix. The choice of spatial weight matrices is important matter in empirical analysis.
Thus, applying other spatial weight matrices can improve our volatility analysis using the S-GARCH model.
Moreover, spatiotemporal extension of the S-GARCH model which considers eﬀects from both space andtime would make it possible to analyze volatility structures in more detail.

7

Table 1: The empirical means and square root of mean squared errors (RMSE) of the estimators.
normaln=100ϕBias RMSE
0.9
0.029
0.082
0.05 -0.039
0.069
0.5
0.039
0.378
1.0
0.021
0.188
0.45 -0.060
0.238
0.45 -0.001
0.155
0.5
-0.014
0.292
1.0
0.034
0.232
0.05 -0.027
0.139
0.9
-0.0110.108
0.5-0.431
0.829
1.0
0.013
0.236
Note: ϕ = (λ, ρ, α, β)′

n=400
Bias
0.007
-0.0090.006
0.009-0.0150.002-0.0020.012-0.011
0.002-0.100
0.007

RMSE
0.029
0.026
0.105
0.089
0.098
0.072
0.092
0.113
0.080
0.069
0.240
0.114

chi(3)
n=100
Bias
0.032
-0.0400.015
0.023-0.0650.002-0.0540.044-0.023-0.017-0.627
0.022

RMSE
0.078
0.066
0.310
0.176
0.243
0.159
0.313
0.229
0.141
0.115
1.089
0.228

n=400
Bias
0.009
-0.010-0.0030.004-0.0150.002-0.0070.018-0.011
0.002-0.114
0.006

RMSE
0.030
0.027
0.100
0.082
0.103
0.075
0.113
0.109
0.079
0.068
0.295
0.109

log normal
n=100
Bias RMSE
0.031
0.080
-0.0380.068-0.037
0.310
0.020
0.173
-0.053
0.224
0.002
0.151
-0.277
0.595
0.041
0.216
-0.0290.132-0.0060.108-1.009
1.660
0.007
0.215

n=400
Bias
0.009
-0.011-0.0180.004-0.0160.003-0.0860.007-0.0130.003-0.313
0.004

RMSE
0.029
0.026
0.101
0.077
0.097
0.073
0.255
0.103
0.079
0.068
0.630
0.105

Table 2: Estimated values and standard errors of λ and ρ and estimated values of α and β in the S-ARCH
model and the S-GARCH model applied to the residuals by fitting the SAR model to log returns of land
priced data year by year.
S-ARCH
2010
2011
λ̂
se(λ)
ρ̂
se(ρ)
α̂
β̂
AIC

0.240
0.083
0.569
-0.022
1538.70.244
0.081-0.5180.212
1481.7

2012

2013

2014

0.274
0.082
-0.6060.232
1549.80.279
0.083-0.1930.109
1573.80.184
0.084-0.8040.225
1537.7

8

S-GARCH
2010
2011
0.772
0.845
0.206
0.139
0.110
0.076
0.077
0.055
0.162 -0.121-0.001
0.052
1536.6 1475.3

2012
0.874
0.128
0.059
0.048
-0.1300.049
1547.9

2013
0.893
0.100
0.060
0.045
-0.021
0.025
1570.4

2014
0.601
0.415
0.104
0.086
-0.4120.120
1537.9Figure 1: The identified volatility in 2010 and 2011. Notice that the great earth quake occurred in 2011.

Figure 2: A comparison between the estimated volatility of the S-ARCH model and that of the S-GARCH
model.

A

Hessian, average Hessian and symmetric matrix Ωψ,n2∂
The Hessian matrix Hn (ψ) ≡ ∂ψ∂ψ′ log Ln (ψ) has the elements:

Hββ ′

= −

1 ′ ′−1
X R (λ)Rn−1 (λ)Xn ,σ2 n n9

1 ′ ′−1
X R (λ)V (θ),σ4 n n1Hβρ = − 2 Xn′ Rn′−1 (λ)Rn−1 (λ)Wn Yn ,σ
1 ′ ′−1Hβλ =
X R (λ)(Wn′ Rn′−1 (λ)Vn (θ) + Rn−1 (λ)Wn Vn (θ) − Rn−1 (λ)Wn Yn ),σ2 n nnVn′ (θ)Vn (θ)
Hσ 2 σ 2 =
−
,
2σ 4
σ6
1
Hσ2 ρ = − 4 Yn′ Wn′ Rn′−1 (λ)V (θ),σ
1
Hσ2 λ =(V ′ (θ) − Yn′ )Wn′ Rn′−1 (λ)Vn (θ),σ4 n1Hρρ = − 2 Yn′ Wn′ Rn′−1 (λ)Rn−1 (λ)Wn Yn − tr(Sn−1 (θ)Wn Sn−1 (θ)Wn ),σ
1 ′ ′ ′−1Hρλ =Y W R (λ)(Wn′ Rn′−1 (λ)V (θ) + Rn−1 (λ)Wn Vn (θ) − Rn−1 (λ)Wn Yn ) − tr(Sn−1 (θ)Wn Sn−1 (θ)Wn ),σ2 n n n1
Hλλ =(Y ′ − Vn′ (θ))Wn′ Rn′−1 (λ)(2Wn′ Rn′−1 (λ)Vn (θ) + Rn−1 (λ)Wn Vn (θ) − Rn−1 (λ)Wn Yn )σ2 n+tr(Rn−1 (λ)Wn Rn−1 (λ)Wn ) − tr(Sn−1 (θ)Wn Sn−1 (θ)Wn ).
( ∂2)The average Hessian matrix Σψ,n ≡ −E n1 ∂ψψ
has the elements:′ log Ln (ψ0 )
Hβσ2

Σββ ′
Σβσ2
Σβρ
Σβλ
Σσ 2 σ 2
Σσ 2 ρ
Σσ 2 λ
Σρρ
Σρλ

Σλλ

= −

1
X ′ R′−1 Rn−1 Xn ,nσ02 n n= 0,
1
X ′ R′−1 Rn−1 Wn Sn−1 Xn β0 ,=nσ02 n n
1=
X ′ R′−1 Rn−1 Wn Sn−1 Xn β0 ,nσ02 n n
1
=
,
2σ04
1
=
tr(Wn Sn−1 ),nσ02
1tr(Wn Sn−1 − Wn Rn−1 ),=nσ02
1 ′ ′ ′−1 ′ ′−1 −1
1=
β X S Wn Rn Rn Wn Sn−1 Xn β0 + tr(Rn′ Sn′−1 Wn′ Rn′−1 Rn−1 Wn Sn−1 Rn + Sn−1 Wn Sn−1 Wn ),nσ02 0 n nn
1 ′ ′ ′−1 ′ ′−1 −1
1=
β0 Xn Sn Wn Rn Rn Wn Sn−1 Xn β0 + tr(Rn′ Sn′−1 Wn′ Rn′−1 Rn−1 Wn Sn−1 Rn + Sn−1 Wn Sn−1 Wn )
2
nσ0
n
1
′ ′−1
′ ′−1 −1−1−1
− tr(Rn Sn Wn Rn Rn Wn + Sn Wn Rn Wn ),n
1
1 ′ ′ ′−1 ′ ′−1 −1
β X S Wn Rn Rn Wn Sn−1 Xn β0 + tr(Rn′ Sn′−1 Wn′ Rn′−1 Rn−1 Wn Sn−1 Rn + Sn−1 Wn Sn−1 Wn )=nσ02 0 n nn2
1− tr(Rn′ Sn′−1 Wn′ Rn′−1 Rn−1 Wn + Sn−1 Wn Rn−1 Wn ) + tr(Rn−1 Wn Rn−1 Wn + Wn′ Rn′−1 Rn−1 Wn ).
nn=

The symmetric matrix Ωψ,n has the elemetns:
Ωββ ′
Ωβσ2
Ωβρ

= 0,
µ3
=
X ′ R′−1 1n ,2nσ06 n nnµ3 ∑={(Rn−1 Xn )i }′ (Rn−1 Wn Sn−1 Rn )ii ,nσ04 i10µ3 ∑{(Rn−1 Xn )i }′ (Rn−1 Wn Sn−1 Rn − Rn−1 Wn )ii ,nσ04 inΩβλ=Ωσ 2 σ 2=Ωσ 2 ρ

=Ωσ 2 λ

=

Ωρρ

=

Ωρλ

=

µ4 − 3σ04,
4σ08µ3 ′ ′ ′−1 ′ ′−1µ4 − 3σ04
β0 Xn Sn Wn Rn 1n +tr(Sn−1 Wn ),
6
2nσ0
2nσ06
µ4 − 3σ04µ3 ′ ′ ′−1 ′ ′−1
β0 Xn Sn Wn Rn 1n +tr(Sn−1 Wn − Rn−1 Wn ),
6
2nσ0
2nσ06
n
n
2µ3 ∑ −1
µ4 − 3σ04 ∑
−1
−1
−1
(R
W
S
X
β
)
(R
W
S
R
)
+
{(Rn−1 Wn Sn−1 Rn )ii }2 ,
n
n
0
i
n
n
ii
n
n
n
n
nσ04 i=1nσ04 i=1nµ3 ∑ −1(R Wn Sn−1 Xn β0 )i (2Rn−1 Wn Sn−1 Rn − Rn−1 Wn )ii
nσ04 i=1 nµ4 − 3σ04 ∑ −1
(R Wn Sn−1 Rn )ii (Rn−1 Wn Sn−1 Rn − Rn−1 Wn )ii ,nσ04 i=1 n
n

+
Ωλλ

=

n
2µ3 ∑ −1
(R Wn Sn−1 Xn β0 )i (Rn−1 Wn Sn−1 Rn − Rn−1 Wn )ii
nσ04 i=1 nµ4 − 3σ04 ∑{(Rn−1 Wn Sn−1 Rn − Rn−1 Wn )ii }2 ,nσ04 i=1n+where µ3 and µ4 are the third and fourth moments of vi s, respectively, (Rn−1 Xn )i is the i-th row of (Rn−1 Xn ),(Rn−1 Wn Sn−1 Xn β0 )i is the i-th element of (Rn−1 Wn Sn−1 Xn β0 ) and (Rn−1 Wn Sn−1 Rn )ii , (Rn−1 Wn Sn−1 Rn −Rn−1 Wn )ii
and (2Rn−1 Wn Sn−1 Rn − Rn−1 Wn )ii are the (i, j)th element of (Rn−1 Wn Sn−1 Rn ), (Rn−1 Wn Sn−1 Rn − Rn−1 Wn ) and(2Rn−1 Wn Sn−1 Rn − Rn−1 Wn ), respectively.

B

Some useful Lemmas

Lemma B.1 (Proposition 8.4.13, Bernstein (2009)). Let A and B be matrices. We use γmax and γmin to
denote the largest and smallest eigenvalues of a matrix. If A is symmetric and B is positive semi definite,thenγmin (A)tr(B) ≤ tr(AB) ≤ γmax (A)tr(B).
Lemma B.2 (Lee, 2002, p.256; Lee, 2004, p1918). Let {An } and {Bn } be two two sequences of n×n matricesthat are uniformly bounded in both row and column sums and the elements of an n × n matrix {Cn } be O(1)uniformly. Then
1. the sequence {An Bn } are uniformly bounded in both row and column sums,
2. the elements of Cn Bn have the uniform order O(1), and3. the elements of An are uniformly bounded and tr(An ) = O(n).
Lemma B.3 (Lee, 2004, p1918). The elements, the vi′ s of Vn are assumed to be i.i.d. with zero mean anda finite variance and the fourth moment of the v ′ s is assumed to exist. Suppose that An is a square matrix
with tis column sums being uniformly bounded and elements of the n × K matrix Zn are uniformly bounded.
Let {Bn } be uniformly bounded either in row or column sums and their elements bn,ij have O(1) uniformly
in i and j. Then
1. √1n Zn′ An Vn = Op (1) and
2. n1 E(Vn′ Bn Vn ) = O(1) and n1 [Vn′ Bn Vn − E(Vn′ Bn Vn )] = op (1).
11

C

Proofs of Theorems 1-3

C.1

Proof of Theorem1

The consistency of θ̂ will follow from the uniform convergence of n1 (log Ln (θ) − Qn (θ)) to zero on Θ andthe uniqueness identification condition that, for any ϵ > 0, lim supn→∞ maxθ∈Nϵc (θ0 ) n1 (Qn (θ) − Qn (θ0 )) < 0,where Nϵc (θ0 ) is the complement of an open neighborhood of θ0 in Θ of diameter ϵ (Theorem 3.4 of white(1994)).
C.1.1

Proof of the uniform convergence of n1 (log Ln (θ) − Qn (θ))

First, we shall prove the uniform convergence of n1 (log Ln (θ) − Qn (θ)) to zero on Θ. The proof follows from:(a) inf θ∈Θ σn∗2 (θ) is bounded away from zero,(b) supθ∈Θ |σ̂n2 (θ) − σn∗2 (θ)| = op (1),(c) supθ∈Θ | n1 (log Ln (θ) − Qn (θ))| = op (1).
Proof of (a)Vn∗ (θ)By the definition of Vn∗ (θ),=

Rn−1 (λ)(Sn (θ)Yn − Xn βn∗ (θ)),=

Rn−1 (λ)Sn (θ)Yn − Rn−1 Xn (Xn′ Rn′−1 (λ)Rn−1 (λ)Xn )−1 Xn′ Rn′−1 (λ)Rn−1 (λ)Sn (θ)E(Yn ),=Rn−1 (λ)Sn (θ)Yn − Pn Rn−1 (λ)Sn (θ)E(Yn ),=Mn Rn−1 (λ)Sn (θ)Yn + Pn Rn−1 (λ)Sn (θ)(Yn − E(Yn )),where, Pn = Rn−1 Xn (Xn′ Rn′−1 (λ)Rn−1 (λ)Xn )−1 Xn′ Rn′−1 and Mn = In − Pn .
From the orthogonality between the two symmetric idempotent matrices Mn and Pn , we have,σn∗2 (θ) ===1
E(Vn′∗ (θ)Vn∗ (θ)),
n1
E[Yn′ Sn′ (θ)Rn′−1 (λ)Mn Rn−1 (λ)Sn (θ)Yn + (Yn − E(Yn ))′ Sn′ (θ)Rn′−1 (λ)Pn Rn−1 (λ)Sn (θ)(Yn − E(Yn ))],n1
1
E(Yn′ )Sn′ (θ)Rn′−1 (λ)Mn Rn−1 (λ)Sn (θ)E(Yn ) + tr(Rn′−1 (λ)Rn−1 (λ)V ar(Sn (θ)Yn )).
nn

The matrix Mn is positive semi definite because Mn is a symmetric idempotent matrix (Lemma 14.2.14
of Harville (1997)). Thus, the first term is nonnegative uniformly in θ ∈ Θ.
Because the matrix V ar(Sn (θ)Yn ) is symmetric and γmin V ar(Sn (θ)Yn ) > 0 from the assumption, thematrix is positive semi definite (Theorem 3.25 of Schott (2005)). By Lemma A.1, the second term is1
tr(Rn′−1 (λ)Rn−1 (λ)V ar(Sn (θ)Yn ))n1γmin (Rn′−1 (λ)Rn−1 (λ))tr(V ar(Sn (θ)Yn )),
n
1
≥
c c ,
n r y
> 0, uniformly in θ ∈ Θ.
≥

It follow that inf θ∈Θ σn∗2 (θ) is bounded away from zero.
Proof of (b)V̂n (θ)Noting that=Rn−1 (λ)(Sn (θ)Yn − Xn β̂n (θ)),=

Rn−1 (λ)Sn (θ)Yn − Rn−1 Xn (Xn′ Rn′−1 (λ)Rn−1 (λ)Xn )−1 Xn′ Rn′−1 (λ)Rn−1 (λ)Sn (θ)Yn ,=Mn Rn−1 (λ)Sn (θ)Yn .

Hence,
σ̂n2 (θ) =1 ˆ′V n (θ)V̂n (θ),
n
12

=

1 ′ ′
Y S (θ)Rn′−1 (λ)Mn Rn−1 (λ)Sn (θ)Yn .
n n n

It follows that
)
1 ′ ′
1 (
Y S (θ)Rn′−1 (λ)Mn Rn−1 (λ)Sn (θ)Yn − E Yn′ Sn′ (θ)Rn′−1 (λ)Mn Rn−1 (λ)Sn (θ)Yn
n n n
n
)
1 (
− E (Yn − E(Yn ))′ Sn′ (θ)Rn′−1 (λ)Pn Rn−1 (λ)Sn (θ)(Yn − E(Yn )) ,n
= (Q1 − EQ1 ) − EQ2 ,
(
1 ′ ′
1
′−1
−1
′ ′
′−1
−1
where, Q) 1 = n Yn Sn (θ)Rn (λ)Mn Rn (λ)Sn (θ)Yn and EQ2 = n E (Yn −E(Yn )) Sn (θ)Rn (λ)Pn Rn (λ)Sn (θ)(Yn −E(Yn )) .
p
To show the result, it suﬃcient to show Q1 − EQ1 −→ 0 and EQ2 −→ 0, uniformly in θ ∈ Θ.
pFirst, we show that Q1 − EQ1 −→ 0 uniformly in θ ∈ Θ. By Theorem 1 of Andrews (1992), theuniform convergence of Q1 − EQ1 to zero in probability follows from the pointwise convergence for eachθ ∈ Θ and stochastic equicontinuity of Q1 , i.e., for any ϵ > 0, there exists a positive number δ such thatlim supn→∞ P (supθ∈Θ supθ′ ∈B(θ,δ) > ϵ) < ϵ, where B(θ, δ) denote a closed ball in Θ of radius δ ≥ 0 centeredat θ.
First of all, the pointwise convergence of Q1 − EQ1 will be shown. We have, by the identity: Yn =Sn−1 Xn β0 + Sn−1 Rn Vn ,σ̂n2 (θ) − σn∗2 (θ) =

Q1

=

1 −1
(S Xn β0 + Sn−1 Rn Vn )′ Sn′ (θ)Rn′−1 (λ)Mn Rn−1 (λ)Sn (θ)(Sn−1 Xn β0 + Sn−1 Rn Vn ),n n
1 ′ ′ ′−1 ′(β X S S (θ)Rn′−1 (λ)Mn Rn−1 (λ)Sn (θ)Sn−1 Xn β0 + 2β0′ Xn′ Sn′−1 Sn′ (θ)Rn′−1 (λ)Mn Rn−1 (λ)Sn (θ)Sn−1 Rn Vnn 0 n n n+Vn′ Rn′ Sn′−1 Sn′ (θ)Rn′−1 (λ)Mn Rn−1 (λ)Sn (θ)Sn−1 Rn Vn ),=Q1,1 (θ) + 2Q1,2 (θ) + Q1,3 (θ),=

where Q1,1 (θ) = n1 (β0′ Xn′ Sn′−1 Sn′ (θ)Rn′−1 (λ)Mn Rn−1 (λ)Sn (θ)Sn−1 Xn β0 ),Q1,2 (θ) = n1 (β0′ Xn′ Sn′−1 Sn′ (θ)Rn′−1 (λ)Mn Rn−1 (λ)Sn (θ)Sn−1 Rn Vn ) andQ1,3 (θ) = n1 (Vn′ Rn′ Sn′−1 Sn′ (θ)Rn′−1 (λ)Mn Rn−1 (λ)Sn (θ)Sn−1 Rn Vn ). The two terms Q1,2 (θ) and Q1,3 (θ) arestochastic.
For the second term, the column sums of Sn′−1 Sn′ (θ)Rn′−1 (λ)Mn Rn−1 (λ)Sn (θ)Sn−1 Rn are uniformly boundedfrom assumption 3 and Lemma 2 and E(Q1,2 (θ)) = 0. Thus, the pointwise convergence of Q1,2 (θ)−E(Q1,2 )(θ)
follow from Lemma 3. Similarly, the column sums of Rn′ Sn′−1 Sn′ (θ)Rn′−1 (λ)Mn Rn−1 (λ)Sn (θ)Sn−1 Rn are uniformly bounded and the pointwise convergence of Q1,3 (θ) − E(Q1,3 )(θ) follows from Lemma 3. Therefore,pQ1 − EQ1 −→ 0, for each θ ∈ Θ.
Next, we show that Q1 is stochastic equicontinuous. We have by the mean value theorem:
Q1,ℓ (θ1 ) − Q1,ℓ (θ2 )=≤∂Q1,ℓ (θ̄)(θ2 − θ1 ),
∂θ′
∂
sup
Q (θ) (θ2 − θ1 ),
′ 1,ℓ
∂θ
θ∈Θ

where ℓ = 1, 2, 3 and θ̄ lies between θ1 and θ2 . For stochastic equicontinuous, it suﬃces to show that∂′−1 ′supθ∈Θ ∂θSn (θ)Rn′−1 (λ)Mn Rn−1 (λ)Sn (θ)Sn−1 ,′ Q1,ℓ (θ) = Op (1) by Theorem 21.10 of Davidson (1994). Let Π1 be Sn′′ ′−1 ′′−1
−1
−1
′ ′−1 ′′−1
Π2 be β0 Xn Sn Sn (θ)Rn (λ)Mn Rn (λ)Sn (θ)Sn Rn and Π3 be Rn Sn Sn (θ)Rn (λ)Mn Rn−1 (λ)Sn (θ)Sn−1 Rn .
∂∂The partial derivatives ∂θ′ Π1,ℓ take simple form and consequently ∂θ ′ Π1,ℓ are also uniformly bounded in both′′ ∂row and column sums. For Q1,1 , for any θ, the elements of β0 Xn ∂θ′ Sn′−1 Sn′ (θ)Rn′−1 (λ)Mn Rn−1 (λ)Sn (θ)Sn−1 and∂′−1 ′Sn (θ)Rn′−1 (λ)Mn
Xn β0 are uniformly bounded. Thus, there exists constants c1 and c2 such that |{β0′ Xn′ ( ∂θ
′ Sn
∂
′−1 ′
Rn−1 (λ)Sn (θ)Sn−1 )}i |≤ c1 and |(Xn β0 )i | ≤ c2 where {β0′ Xn′ ( ∂θSn (θ)Rn′−1 (λ)Mn Rn−1 (λ)Sn (θ)Sn−1 )}i and′ Sn∂(Xn β0 )i are the i-th elements of each vector. It follows that ∂θ′ Q1,1 ≤ c1 c2 = O(1). For Q1,2 , for any θ,
( ∂
∂
∂
∂
≤ c3 where ∂θ′ Π1,2,i is the i-th element of ∂θ ′ Π1,2 . Therefore, from Lemma B.3, P∂θ)′ Π1,2,i∂θ ′ Q1,2 >)( −1 )( 1 ∑n∂∂2. For Q1,3 , for any θ, ∂θ′ Π1,3,ij ≤ c4 where where ∂θ′ Π1,3,ij is the
M ≤ P n i=1 c3 vi > M = O n
13

( ∂
)
( 1 ∑n ∑)∂(i, j)th element of ∂θ′ Π1,3 . Thus, from Lemma B.3, P
i=1
j=1 c4 vi vj > M =∂θ ′ Q1,3 > M ≤ P
n∂O(1). Thus, supθ∈Θ ∂θ
′ Q1,ℓ (θ) = Op (1) It follow that Q1 is stochastic equicontinuous. Hence, by Theoremp
1 of Andrews (1992), Q1 − EQ1 −→ 0 uniformly in θ ∈ Θ.
Secondly, we show( that EQ2 −→ 0,) uniformly in θ ∈ Θ. There exist cx such that0 < cx ≤ inf λ∈Λ γmin n1 Xn′ Rn′−1 Rn−1 X from assumption. By Assumption, Lemma 1 and 2 and theorem 3.4
of Schott (2005), We have,
EQ2

=
=
=
≤
=
≤
=

)
1 (
E (Yn − E(Yn ))′ Sn′ (θ)Rn′−1 (λ)Pn Rn−1 (λ)Sn (θ)(Yn − E(Yn )) ,n1tr(Rn′−1 (λ)Pn Rn−1 (λ)V ar(Sn (θ)Yn )),n1tr(Rn′−1 (λ)Rn−1 Xn (Xn′ Rn′−1 (λ)Rn−1 (λ)Xn )−1 Xn′ Rn′−1 Rn−1 (λ)V ar(Sn (θ)Yn )),
n
1 −1
2
γ
(X ′ R′−1 Rn−1 X)γmax(Rn′−1 (λ)Rn−1 (λ))γmax (V ar(Sn (θ)Y ))tr(Xn′ Xn )),n min ( n n)1
1 −1 Xn′ Rn′−1 Rn−1 X
2γγmax(Rn′−1 (λ)Rn−1 (λ))γmax (V ar(Sn (θ)Y )) tr(Xn′ Xn )),
n min
n
n
1 −1 2 1c c cy tr(Xn′ Xn )),n (x r) nO n−1

Hence, EQ2 −→ 0, uniformly in θ ∈ Θ.
Therefore, supθ∈Θ |σ̂n2 (θ) − σn∗2 (θ)| = op (1), completing the proof of (b).
Proof of (C)We show that supθ∈Θ n1 (log Ln (θ) − Qn (θ)) = op (1). Note that11(log Ln (θ) − Qn (θ)) = − (log σ̂n2 (θ) − log σn∗2 (θ)).
n2

By the Taylor expansion,log σ̂n2 (θ) − log σn∗2 (θ) =1
σ̃n2 (θ)σ̂n2 (θ) − σn∗2 (θ) ,where σ̃n2 (θ) lies between σ̂n2 (θ) and σn∗2 (θ). From the proof (a) and (b), it follow that σ̂n2 (θ) is uniformly
bounded away from zero on Θ. Moreover, σ̃n2 (θ) is also uniformly bounded away from zero on Θ becauseσ̃n2 (θ) exists between σ̂n2 (θ) and σn∗2 (θ) and thereby σ̃21(θ) is uniformly bounded. As σ̂n2 (θ) − σn∗2 (θ) covergesnin probability to zero uniformly on Θ, | log σ̂n2 (θ) − log σn∗2 (θ)| = op (1) uniformly on Θ.
Consequently, supθ∈Θ n1 (log Ln (θ) − Qn (θ)) = op (1).
C.1.2

Proof of the identification uniqueness condition

Secondly, we shall prove the identification uniqueness condition. The proof follow from:(i) n1 Qn (θ) is uniformly equicontinuous on Θ.
(ii) Show some properties of an auxiliary model.
(iii) Show that the identification uniqueness condition holds.
Proof of (i) We show that n1 Qn (θ) = 21 (log 2π + 1) − 12 log σn∗2 (θ) − n1 log |Rn (λ)| + n1 log |Sn (θ)| is uniformly
equicontinuous on Θ. It is suﬃcient to show that partial derivatives of each term are uniformly bounded. The
uniform continuity of log σn∗2 (θ) on Θ follows because σ∗21(θ) is uniformly bounded since σn∗2 (θ) is uniformly
n

d 1
1
−1
bounded away form zero on Θ. For n1 log |Rn (λ)|, dλn log |Rn (λ)| = n tr(Rn (λ)Wn ). From assumption andLemma 2, the elements of Rn−1 (λ)Wn are uniformly bounded. Thus, n1 tr(Rn−1 (λ)Wn ) = O(1) from Lemma∂ 1
1
2. Similarly, ∂θn log |Sn (θ)| = O(1). Hence, n Qn (θ) is uniformly equicontinuous on Θ.

14

Proof of (ii)

It is useful to establish an auxiliary process:Yn = λWn Yn + ρWn Yn + Rn (λ)Vn ,where Vn ∼ N (0, σ02 In ). The log-likelihood function of the above auxiliary process is given bylog Lp,n (θ, σ 2 ) =

−

n
n
1
log(2π) − log(σ 2 (θ)) − log |Rn (λ)| + log |Sn (θ)| − 2 Yn′ Sn′ (θ)Rn′−1 (λ)Rn−1 (λ)Sn (θ)Yn .
222σ

Let Ep be the expectation under this auxiliary process. Define Qp,n (θ) = maxσ2 Ep (log Lp,n (θ)). The
optimal solutions of this maximization problem is
σn2 (θ)==1
Ep (Yn′ Sn′ (θ)Rn′−1 (λ)Rn−1 (λ)Sn (θ)Yn ),nσ2tr(Rn Sn−1 Sn′ (θ)Rn′−1 (λ)Rn−1 (λ)Sn (θ)Sn−1 Rn ).
nHence,
Qp,n (θ) = −nnlog(2π + 1) + log σn2 (θ) − log |Rn (λ)| + log |Sn (θ)|.
22

By Shannon-Kolmogorov Information Inequality (Ferguson (1996), p113), Qp,n (θ) ≤ Qp,n (θ0 ) for allθ ∈ Θ. This implies that n1 (Qp,n (θ) − Qp,n (θ0 ) ≤ 0 for all θ ∈ Θ.
Proof of (iii)

We show that the identification uniqueness condition holds by contradiction.
(
)
1
1
1
(Qn (θ) − Qn (θ0 )) = − log σn∗2 (θ) − log |Rn (λ)| + log |Sn (θ)| − − log σ02 − log |Rn | + log |Sn |
n
2
2
(
)
1
1
1
=
− (log σn2 (θ) − log σ02 ) − (log |Rn (λ)| − log |Rn |) + (log |Sn (θ)| − log |Sn |)2nn1− (log σn∗2 (θ) − log σn2 (θ)),
n
) 1
1(
=
Qp,n (θ) − Qp,n (θ0 ) − (log σn∗2 (θ) − log σn2 (θ)).
n2Moreover,σn∗2 (θ) − σn2 (θ) =1 ′ ′ ′−1 ′β X S S (θ)Rn′−1 (λ)Mn Rn−1 (λ)Sn (θ)Sn−1 Xn β0 .
n 0 n n n

Mn is positive semi definite and thereby σn∗2 (θ) − σn2 (θ) ≥ 0. This implies − 21 (log σn∗2 (θ) − log σn2 (θ)) ≤ 0.
Now, suppose that the identification uniquenesscondition does
not hold. Then, there exists an ϵ > 0 and a()
sequence {θn } in Nϵc (θ0 ) such that limn→∞ n1 Qn (θ) − Qn (θ0 ) = 0. By the compactness of Nϵc (θ0 ), there exists a convergent subsequence {θnm } of {θn } with the limit (θ+ of θnm being in Nϵc)(θ0 ). This implies(that θ+ ̸=
θ0 . As n1 Qn (θ) is uniformly equicontinuous, limnm →∞ n1m Qnm (θ+ ) − Qnm (θ0 ) = 0. Because n1 Qp,n (θ) −
)
(
)
(
)
Qp,n (θ0 ) ≤ 0 and − 21 log σn∗2 (θ)−log σn2 (θ) ≤ 0, this is possible only if limnm →∞ n1m Qnm (θ+ )−Qnm (θ0 ) =()
1 ′′ ′−1 ′′−1
−1
−1
0 and − 12 log σn∗2 (θ)−log σn2 (θ) ≤ 0. However, limn→∞n β0 Xn Sn Sn (θ)R() n (λ)Mn Rn (λ)Sn (θ)Sn Xn β0 ̸=1∗2
2
0 from the assumptionin Theorem)3.1 . Thus, − 2 log σn (θ) − log σn (θ) < 0 and consequently(
1limnm →∞ nm Qnm (θ+ ) − Qnm (θ0 ) ̸= 0. This is a contradiction. Therefore, the identification uniquenesscondition must hold.
The consistency of θ̂ follow form uniform convergence and the identification uniqueness condition. This
completes the proof of the theorem.

C.2

Proof of Theorem 2

We have by the Taylor expansion,
0 =1 ∂ log Ln (ψ̂n )
√
,
∂ψ
n
15

=

1 ∂ log Ln (ψ0 )
√
+
∂ψ
n

(

)
1 ∂ 2 log Ln (ψ̄n ) √n(ψ̂n − ψ0 ),n∂ψ∂ψ ′

where ψ̄n lies between ψ̂n and ψ0 . Thus, the asymptotic normality of ψ̂n follows if()
Ln (ψ0 ) D(a) √1n ∂ log ∂ψ−→ N 0, limn→∞ Γ(ψ0 ) ,( 2 Ln (ψ0 ) ) p
2
Ln (ψ0 )(b) n1 ∂ log− E n1 ∂ log−→ 0, and∂ψ∂ψ ′∂ψ∂ψ ′22p

Ln (ψ̄n )
Ln (ψ0 )(c) n1 ∂ log−→ 0.
− n1 ∂ log∂ψ∂ψ ′∂ψ∂ψ ′Ln (ψ0 )Proof of (a) The asymptotic normality of √1n ∂ log ∂ψfollows from the central limit theorems for linearquadratic forms in Kelejian and Prucha (2001). We need to check that the score vector holds Assumptions in Kelejian and Prucha (2001). To check assumptions for asymptotic normality, it is suﬃcient toshow some matrices hold desired boundaly conditions. From assumptions of this paper and Lemma A.2,
(Rn′ Sn′−1 Wn′ Rn′−1 − Wn′ Rn′−1 ) and Rn′ Sn′−1 Wn′ Rn′−1 are uniformly bounded in column sums, and the elementsof Xn′ Sn′−1 Wn′ Rn′−1 are uniformly bounded. Thus, each score function holds the assumptions and the asymptotic normality of each score function follows. Finally, the Cramér-Wold devise (Proposition 6.3.1 of Brockwell
and Davis (1991)) leads to the joint asymptotic normality.

Proof of (b)( 2 Ln (ψ0 ) )2
Ln (ψ0 )− E n1 ∂ log. Then, Dψψ has the elements:
Let Dψψ be n1 ∂ log∂ψ∂ψ ′∂ψ∂ψ ′

Dββ ′

= 0,

Dβσ2

= −

Dβρ

=

Dβλ

=

Dσ 2 σ 2

=

Dσ2 ρ

=

Dσ2 λ

=

Dρρ

=

Dρλ

=

Dλλ

=

1
X ′ R′−1 Vn ,nσ04 n n
1− 2 Xn′ Rn′−1 Rn−1 Wn Sn−1 Rn Vn ,nσ01
X ′ (R′−1 Wn′ Rn′−1 + Rn′−1 Rn−1 W − Rn′−1 Rn−1 Wn Sn−1 Rn )Vn ,nσ02 n n
1
1−
V ′ Vn ,σ04
nσ06 n1
1− 4 β0′ Xn′ Sn′−1 Wn′ Rn′−1 Vn −(V ′ R′ S ′−1 Wn′ Rn′−1 Vn − σ02 tr(Sn′−1 Wn′ )),nσ0
nσ04 n n n
1
1(V ′ W ′ R′−1 Vn − σ02 tr(Wn′ Rn′−1 ))− 4 β0′ Xn′ Sn′−1 Wn′ Rn′−1 Vn +nσ0nσ04 n n n
1− 4 (Vn′ Rn′ Sn′−1 Wn′ Rn′−1 Vn − σ02 tr(Sn′−1 Wn′ )),nσ02 ′ ′ ′−1 ′ ′−1 −1−
β X S Wn Rn Rn Wn Sn−1 Rn Vn
nσ0 0 n n
1
− 2 (Vn′ Rn′ Sn′−1 Wn′ Rn′−1 Rn−1 Wn Sn−1 Rn Vn − σ02 tr(Rn′ Sn′−1 Wn′ Rn′−1 Rn−1 Wn Sn−1 Rn )),nσ01 ′ ′ ′−1 ′ ′−1 ′ ′−1
β X (S Wn Rn Wn Rn + Sn′−1 Wn′ Rn′−1 Rn−1 Wn − 2Sn′−1 Wn′ Rn′−1 Rn−1 Wn Sn−1 Rn )Vn
nσ02 0 n n
1+ 2 (Vn′ Rn′ Sn′−1 Wn′ Rn′−1 Wn′ Rn′−1 Vn − σ02 tr(Sn′−1 Wn′ Rn′−1 Wn′ ))nσ01+ 2 (Vn′ Rn′ Sn′−1 Wn′ Rn′−1 Rn−1 Wn Vn − σ02 tr(Rn′ Sn′−1 Wn′ Rn′−1 Rn−1 Wn ))nσ01− 2 (Vn′ Rn′ Sn′−1 Wn′ Rn′−1 Rn−1 Wn Sn−1 Rn Vn − σ02 tr(Rn′ Sn′−1 Wn′ Rn′−1 Rn−1 Wn Sn−1 Rn )),
nσ0
1 ′ ′
β X (2Sn′−1 Wn′ Rn′−1 Wn′ Rn′−1 + Sn′−1 Wn′ Rn′−1 Rn−1 Wn − 2Sn′−1 Wn′ Rn′−1 Rn−1 Wn Sn′−1 Rnnσ02 0 n16

−2Rn′−1 Wn′ Rn′−1 Wn′ Rn′−1 − Rn′−1 Wn′ Rn′−1 Rn−1 Wn + 2Rn′−1 Wn′ Rn′−1 Rn−1 Wn Sn−1 Rn+2Rn′−1 Wn′ Rn′−1 Wn′ Rn−1 + Rn′−1 Wn′ Rn′−1 Rn−1 Wn − Rn′−1 Wn′ Rn′−1 Rn−1 Wn Sn−1 Rn )Vn2+ 2 (Vn′ Rn′ Sn′−1 Wn′ Rn′−1 Wn′ Rn′−1 Vn − σ02 tr(Sn′−1 Wn′ Rn′−1 Wn′ ))nσ01+ 2 (Vn′ Rn′ Sn′−1 Wn′ Rn′−1 Rn−1 Wn Vn − σ02 tr(Rn′ Sn′−1 Wn′ Rn′−1 Rn−1 Wn ))nσ01− 2 (Vn′ Rn′ Sn′−1 Wn′ Rn′−1 Rn−1 Wn Sn−1 Rn Vn − σ02 tr(Rn′ Sn′−1 Wn′ Rn′−1 Rn−1 Wn Sn−1 Rn ))nσ02− 2 (Vn′ Wn′ Rn′−1 Wn′ Rn′−1 Vn − σ02 tr(Wn′ Rn′−1 Wn′ Rn′−1 ))nσ01+ 2 (Vn′ Wn′ Rn′−1 Rn−1 Wn Vn − σ02 tr(Wn′ Rn′−1 Rn−1 Wn ))nσ01+ 2 (Vn′ Wn′ Rn′−1 Rn−1 Wn Sn−1 Rn Vn − σ02 tr(Wn′ Rn′−1 Rn−1 Wn Sn−1 Rn )).
nσ0Thus, the elements of Dψψ are decomposed into sums of the forms: n1 Xn′ An (θ)Vn , n1 β0′ Xn′ An (θ)Vn ,
1
1
1
′
′
′
n (Vn An (θ)Vn − E(Vn An (θ)Vn )) and σ04 − nσ06 Vn Vn , where a matrix An (θ) is uniformly bounded in both row
and column sums. From Lemma A.3, n1 Xn′ An (θ)Vn , n1 β0′ Xn′ An (θ)Vn and n1 (Vn′ An (θ)Vn − E(Vn′ An (θ)Vn )) areppconvergence to zero in probability. Moreover, σ14 − nσ1 6 Vn′ Vn −→ 0 because n1 Vn Vn −→ σ02 by the law of large00( 2 Ln (ψ0 ) ) p
2
Ln (ψ0 )numbers. Therefore, it follow that n1 ∂ log− E n1 ∂ log−→ 0.
∂ψ∂ψ ′∂ψ∂ψ ′
Proof of (c)22

Ln (ψ̄n )
Ln (ψ0 )From Lemma B.2 and B.3, it is easy to show that n1 ∂ log= Op (1) and n1 ∂ log=∂ψ∂ψ ′∂ψ∂ψ ′p∂
Op (1). Here, σ̄ −r = σ0−r + op (1), r = 2, 4, 6 because σ̄ 2 −→ σ02 and σ r appears in Hn (ψ) ≡ ∂ψ∂ψ′ log Ln (ψ)
22multiplicatively, thus it results in an asymptotically negligible error to replace σ̄ by σ0 . The elements of the∂2′
Hessian matrix, Hn (ψ) ≡ ∂ψ∂ψ′ log Ln (ψ), are decomposed into sums of terms of the forms: Xn An (θ)Xn ,Xn′ An (θ)Yn , Xn′ An (θ)V (θ), Yn′ An (θ)Yn , 2σn4 − σ16 Vn′ (θ)Vn (θ), Yn′ An (θ)Vn (θ), Vn′ (θ)An (θ)Vn (θ) and tr(An (θ)),
where a matrix An (θ) is uniformly bounded in both row and column sums. Therefore, it is suﬃcient to showthat the diﬀerence between each term at ψ̄ and ψ0 converges to zero in probability and moreover this can beeasily shown. We show some examples corresponding each term of the Hessian matrix.
Noting thatRn−1 (λ) − Rn−1

2

= Rn−1 (λ)(Rn − Rn (λ))Rn−1 ,= (λ0 − λ)Rn−1 (λ)Wn Rn−1 .

For Xn′ An (θ)Xn ,1 ′ ′−1
1
X R (λ̄)Rn−1 (λ̄)Xn − Xn′ Rn′−1 Rn−1 Xnn n nn

1 ′ ′−1
1
X (R (λ̄) − Rn′−1 + Rn′−1 )Rn−1 (λ̄)Xn − Xn′ Rn′−1 Rn−1 Xn ,n n nn
1 ′ ′−1=
X (R (λ̄) − Rn′−1 )Rn−1 (λ̄)Xn +n n n
1 ′ ′−1 −1
1
X R Rn (λ̄)Xn − Xn′ Rn′−1 Rn−1 Xn ,
n n n
n
1
= (λ0 − λ̄) Xn Rn′−1 (λ)Wn′ Rn−1 Rn−1 (λ̄)Xnn1+(λ0 − λ̄) Xn′ Rn′−1 Rn−1 (λ)Wn Rn−1 Xn ,n
= op (1)O(1) + op (1)O(1),== op (1).
Moreover, the convergence of Xn′ An (θ)Yn is shown similarly.

17Noting thatRn−1 (λ)Rn (λ)Vn (θ),

Vn (θ) == Rn−1 (λ)(S(θ)Yn − Xn β),= Rn−1 (λ)((λ0 − λ)Wn Yn + (ρ0 − ρ)Wn Yn + Xn (β0 − β) + Rn Vn ).
Thus, for Xn′ An (θ)V (θ),
1 ′ ′−1
1
X R (λ̄)Vn (θ̄) − Xn′ Rn′−1 Vnn n nn) 1 ′ ′−1
1
X R (λ̄)Wn Yn + Xn′ Rn′−1 (λ̄)Xn (β0 − β)
n n n
n
1
1
+ Xn′ Rn′−1 (λ̄)Rn Vn − Xn′ Rn′−1 Vn ,nn
= op (1)Op (1) + Op (1)op (1) + op (1) + op (1),=((λ0 − λ̄) + (ρ0 − ρ̄)= op (1),where the convergence of last two terms follow from Lemma B.3.
Here,
)2 1 ′ ′ ′−1
Y W R (λ̄))Rn−1 (λ̄)Wn Ynn n n n11+(β0 − β)′ Xn′ Rn′−1 (λ̄)Rn−1 (λ̄)Xn (β0 − β) + Vn′ Rn′ Rn′−1 (λ̄)Rn−1 (λ̄)Rn Vnnn) ′ ′ ′−12(
−1
+ (λ0 − λ̄) + (ρ0 − ρ̄) Yn Wn Rn (λ̄)Rn (λ̄)Xn (β0 − β)n)2(2
+ (λ0 − λ̄) + (ρ0 − ρ̄) Yn′ Wn′ Rn′−1 (λ̄)Rn−1 (λ̄)Rn Vn + (β0 − β)′ Xn′ Rn′−1 (λ̄)Rn−1 (λ̄)Rn Vn ,nn
= op (1)Op (1) + op (1)O(1)op (1) + σ02 + op (1)Op (1)op (1) + op (1)Op (1) + op (1)op (1),1 ′V (θ̄)Vn (θ̄) =n n((λ0 − λ̄) + (ρ0 − ρ̄)= σ02 + op (1).
It follows that 2σ1 4 − nσ1 6 Vn′ (θ)Vn (θ) = op (1).
00Before next proof, we show an example. Yn′ Sn (θ)Vn = β ′ Xn′ Sn−1 S(θ)Vn + Vn′ Rn′ Sn−1 Sn (θ)Vn and1 ′ ′ −1
1
V R S Sn (θ)Vn − Vn′ Rn′ Sn−1 Sn Vnn n n nn=)1(λ0 − λ) + (ρ0 − ρ) Vn′ Rn′ Sn−1 Vn ,n
op (1)Op (1),=op (1).

=(

It follows that n1 Yn′ Sn (θ)Vn − n1 Yn′ Sn Vn = op (1) and similarly n1 Yn′ An (θ)Vn − n1 Yn′ An Vn = op (1) and1 ′1 ′n Yn An (θ)Yn − n Yn An Yn = op (1) where An is An (θ) at true value θ0 .
′Now, for Yn An (θ)Vn (θ),
1 ′ ′ ′−1
1
Yn Wn Rn (λ)Vn (θ) − Yn′ Wn′ Rn′−1 Vnnn)1(λ0 − λ̄) + (ρ0 − ρ̄) Yn′ Wn′ Rn′−1 (λ)Rn−1 (λ)Wn Ynn1 ′ ′ ′−1
−1
+ Yn Wn Rn (λ)Rn (λ)Xn (β0 − β̄)n11+ Yn′ Wn′ Rn′−1 (λ̄)Rn−1 (λ̄)Rn V − Yn′ Wn′ Rn′−1 Vnnn= op (1)Op (1) + Op (1)op (1) + op (1)=(= op (1).
Moreover, the convergence of Vn (θ)′ An (θ)Vn (θ) is also shown similary.
Finally, for tr(An (θ)), by the Taylor expansion,
11
tr(Rn−1 (λ)Wn Rn−1 (λ)Wn ) − tr(Rn−1 Wn Rn−1 Wn ) =
n
n
18

d
tr(Rn−1 (λ̃)Wn Rn−1 (λ̃)Wn )(λ̄ − λ0 ),dλ= O(1)op (1),= op (1),where λ̃ lies between λ̄ and λ0 .
2
Ln (ψ̄n )The convergence of the other elements of the Hessian matrix are shown similarly, hence n1 ∂ log−∂ψ∂ψ ′
2
1 ∂ log Ln (ψ0 ) p−→ 0.
n∂ψ∂ψ ′

This completes the proof of the theorem.

C.3

Proof of Theorem 3

The estimator for α is(α̂n = (1 − λ̂) log)n1∑exp{(Rn−1 (λ̂)[S(θ̂)Yn − Zn δ̂])i } ,
n i=1

Here,
S(θ̂)Yn − Zn δ̂= Yn − λ̂Wn Yn − ρ̂Wn Yn − Zn δ̂,=(λ0 − λ̂)Wn Yn + (ρ0 − ρ̂)Wn Yn + Zn (δ0 − δ̂) + α0 1n + Rn Vn ,= D + α0 1n + Rn Vn ,where D = (λ0 − λ̂)Wn Yn + (ρ0 − ρ̂)Wn Yn + Zn (δ0 − δ̂).
α0
Because Rn−1 (λ̂)(S(θ̂)Yn − Zn δ̂) = 1−
1 + Rn−1 (λ̂)D + Rn−1 (λ̂)Rn Vn ,
λ̂ n
) ∑
(
n
n
1∑
1
α
exp{(Rn−1 (λ̂)[S(θ̂)Yn − Zn δ̂])i } = exp
exp{(Rn−1 (λ̂)D + Rn−1 (λ̂)Rn Vn )i }.
n i=1
1 − λ n i=1Thus,(α̂ − α0 = (1 − λ̂) log)n1∑exp{(Rn−1 (λ̂)D + Rn−1 (λ̂)Rn Vn )i } .
n i=1

To prove consistency, it is suﬃcient that the right side of (9) converges to zero in probability.
By the Taylor expansion,1∑exp{(Rn−1 (λ̂)D + Rn−1 (λ̂)Rn Vn )i }
n i=1

= 1+

n
{
}
1∑
exp(bi ) (Rn−1 (λ̂)D + Rn−1 (λ̂)Rn Vn )i
n i=1

= 1+

1 ′ −1b (Rn (λ̂)D + Rn−1 (λ̂)Rn Vn ),nnwhere bi lies between 0 and (Rn−1 (λ̂)D + Rn−1 (λ̂)Rn Vn )i , and b = (b1 , . . . , bn )′ .
From Assumptions, Theorem 1 and Lemma B.2 and B.3,
1 ′ −1b (Rn (λ̂)D + Rn−1 (λ̂)Rn Vn )n11
= (λ0 − λ̂) b′ Rn−1 (λ̂)Wn Yn + (ρ0 − ρ̂) b′ Rn−1 (λ̂)Wn Ynnn
1 ′ −1
1 ′ −1+ b Rn (λ̂)Zn (δ0 − δ̂) + b Rn (λ̂)Rn Vn ,nn
= op (1)Op (1) + op (1)Op (1) + O(1)op (1) + op (1),= op (1).

∑npThus, n1 i=1 exp{(Rn−1 (λ̂)D + Rn−1 (λ̂)Rn Vn )i } −→ 1 and)( 1 ∑np(1 − λ̂) log n i=1 exp{(Rn−1 (λ̂)D + Rn−1 (λ̂)Rn Vn )i } −→ 0.

19

(9)

References[1] Andrews, D. W. (1992). Generic uniform convergence. Econometric theory, 8(2), 241-257.
[2] Bernstein, D. S. (2009). Matrix mathematics: Theory, facts, and formulas. 2nd edition. Princeton:
Princeton University Press.
[3] Bollerslev, T. (1986). Generalized autoregressive conditional heteroskedasticity. Journal of econometrics,
31(3), 307-327.
[4] Bollerslev, T. (1990). Modelling the coherence in short-run nominal exchange rates: a multivariate
generalized ARCH model. The review of economics and statistics, 498-505.
[5] Bollerslev, T., Engle, R. F., & Wooldridge, J. M. (1988). A capital asset pricing model with time-varying
covariances. Journal of political Economy, 96(1), 116-131.
[6] Borovkova, S. & Lopuhaa, R., (2012). Spatial GARCH: A Spatial Approach to Multivariate Volatility
Modeling. Available at SSRN: https://ssrn.com/abstract=2176781.
[7] Brockwell, P. J., & Davis, R. A. (1991). Time series: theory and methods. Springer-Verlag.
[8] Caporin, M., and Paruolo, P., 2009, Structured multivariate volatility models, Available at SSRN:
http://ssrn.com/abstract=1318639.
[9] Davidson, J. (1994). Stochastic limit theory. Oxford: Oxford University Press.
[10] Doan, O., & Tapnar, S. (2013). GMM estimation of spatial autoregressive models with moving average
disturbances. Regional Science and Urban Economics, 43(6), 903-926.
[11] Engle, R. F. (1982). Autoregressive conditional heteroscedasticity with estimates of the variance ofUnited Kingdom inflation. Econometrica, 50, 987-1007.
[12] Engle, R. F., & Kroner, K. F. (1995). Multivariate simultaneous generalized ARCH. Econometric theory,11(1), 122-150.
[13] Fan, J., & Yao, Q. (2003). Nonlinear time series: nonparametric and parametric methods. Springer
Science & Business Media.
[14] Ferguson, T. S. (1996). A course in large sample theory. Chapman & Hall.
[15] Harville, D. A. (1997). Matrix algebra from a statistician’s perspective. New York: Springer.
[16] Kelejian, H. H., & Prucha, I. R. (1997). Estimation of spatial regression models with autoregressive
errors by two-stage least squares procedures: a serious problem. International regional science review,20(1-2), 103-111.
[17] 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.
[18] Kelejian, H. H., & Prucha, I. R. (2001). On the asymptotic distribution of the Moran I test statisticwith applications. Journal of Econometrics, 104(2), 219-257.
[19] Kelejian, H. H., & Robinson, D. P. (1993). A suggested method of estimation for spatial interdependent
models with autocorrelated errors, and an application to a county expenditure model. Papers in regional
science, 72(3), 297-312.
[20] Lee, L. F. (2002). Consistency and eﬃciency of least squares estimation for mixed regressive, spatial
autoregressive models. Econometric theory, 18(2), 252-277.
[21] Lee, L. F. (2004). Asymptotic Distributions of QuasiMaximum Likelihood Estimators for Spatial Autoregressive Models. Econometrica, 72(6), 1899-1925.
20

[22] Lee, L. F. (2007). GMM and 2SLS estimation of mixed regressive, spatial autoregressive models. Journalof Econometrics, 137(2), 489-514.
[23] Lee, P. M. (2012). Bayesian statistics: an introduction. John Wiley & Sons.
[24] Robinson, P. M. (2009). Largesample inference on spatial dependence. The Econometrics Journal, 12(s1).
[25] Sato, T. and Matsuda, Y. (2017). Spatial Autoregressive Conditional Heteroskedasticity Models. J.Japan
Statist. Soc., Vol. 47, 2[26] Schott, J. R. (2005). Matrix analysis for statistics. 2nd edition. John Wiley & Sons.
[27] Su, L., & Yang, Z. (2015). QML estimation of dynamic panel data models with spatial errors. Journalof Econometrics, 185(1), 230-258.
[28] White, H. (1994). Estimation, inference and specification analysis. New York: Cambridge university
press.
[29] Yan, J. (2007). Spatial stochastic volatility for lattice data. Journal of agricultural, biological, andenvironmental statistics, 12(1), 25.
[30] Yang Z. (2015). Unified M-Estimation of Fixed-Eﬀects Spatial Dynamic Models with Short Panels,
Working Papers 14-2015, Singapore Management University, School of Economics.
[31] Yu, J., De Jong, R., & Lee, L. F. (2008). Quasi-maximum likelihood estimators for spatial dynamic
panel data with fixed eﬀects when both n and T are large. Journal of Econometrics, 146(1), 118-134.

21