Partially Specified Space Time Autoregressive Model with Artificial Neural Network

TL;DR

This paper extends the classical space-time autoregressive model by integrating exogenous variables, non-Gaussian errors, and a neural network component, enhancing flexibility and modeling complex spatial-temporal relationships.

Contribution

It introduces a partially specified space-time autoregressive model with neural networks, providing estimation methods and theoretical properties like consistency and asymptotic normality.

Findings

Neural network component improves model flexibility.

Simulation studies validate finite sample approximations.

Real-world application demonstrates practical utility.

Abstract

The space time autoregressive model has been widely applied in science, in areas such as economics, public finance, political science, agricultural economics, environmental studies and transportation analyses. The classical space time autoregressive model is a linear model for describing spatial correlation. In this work, we expand the classical model to include related exogenous variables, possibly non-Gaussian, high volatility errors, and a nonlinear neural network component. The nonlinear neural network component allows for more model flexibility, the ability to learn and model nonlinear and complex relationships. We use a maximum likelihood approach for model parameter estimation. We establish consistency and asymptotic normality for these estimators under some standard conditions on the space time model and neural network component. We investigate the quality of the asymptotic approximations for finite samples by means of numerical simulation studies. For illustration, we include a real world application.

Tables

Table 1. Table 1 : Empirical means and standard errors (in parentheses) of parameter estimates when ε 𝜀 \varepsilon is sampled from a standard normal, standardized student t distribution and a Laplace distribution. The asymptotic standard errors are displayed for reference in square brackets.

Model 1: $Y_t={\varphi }_0{W}_n{Y}_t+{\varphi }_1{W}_n{Y}_{t-1}+𝑭\left(X_t{𝜸}^{\prime }\right)\lambda +{\epsilon }_t$
		${\epsilon }_t$		${\hat{\varphi }}_0$	${\hat{\varphi }}_1$	$\hat{\lambda }$	${\hat{\gamma }}_1$	${\hat{\gamma }}_2$
		true value		$0.6$	$-0.274$	$1.50$	$0.75$	$-0.35$
		$N(0,1)$		0.5997	-0.2743	1.5025	0.7485	-0.3476
				(0.0065)	(0.0079)	(0.0274)	(0.0269)	(0.0134)
				[0.0079]	[0.0085]	[0.0308]	[0.0310]	[0.0147]
		$t\left(4\right)$		0.5994	-0.2737	1.5000	0.7531	-0.3507
				(0.0059)	(0.0069)	(0.0236)	(0.0249)	(0.0112)
				[0.0068]	[0.0071]	[0.0259]	[0.0258]	[0.0122]
		$Laplace$		0.5999	-0.2736	1.4992	0.7501	-0.3504
		$(0,\frac{\sqrt{2}}{2})$		(0.0048)	(0.0058)	(0.0199)	(0.0196)	(0.0097)
Model 2: $Y_t={\varphi }_0{W}_n{Y}_t+{\varphi }_1{W}_n{Y}_{t-1}+X_t\beta +𝑭\left(X_t{𝜸}_1^{\prime }\right){\lambda }_1+𝑭\left(X_t{𝜸}_2^{\prime }\right){\lambda }_2+{\epsilon }_t$
${\epsilon }_t$	${\hat{\varphi }}_0$	${\hat{\varphi }}_1$	$\hat{\beta }$		${\hat{\lambda }}_1$	${\widehat{𝜸}}_1$		${\hat{\lambda }}_2$	${\widehat{𝜸}}_2$
	$0.6$	$-0.274$	$0.24$	$-0.70$	$2$	$0.75$	$0.7$	$0.8$	$0.35$	$-1$
$N(0,1)$	0.6000	-0.2748	0.2402	-0.6985	1.9927	0.7503	0.7030	0.8076	0.3577	-1.0159
	(0.0039)	(0.0046)	(0.0137)	(0.0140)	(0.0928)	(0.0962)	(0.0369)	(0.0450)	(0.0899)	(0.1209)
	[0.0040]	[0.0044]	[0.0135]	[0.0141]	[0.0921]	[0.0920]	[0.0390]	[0.0449]	[0.0835]	[0.1243]
$t\left(4\right)$	0.5999	-0.2740	0.2402	-0.7005	2.0008	0.7496	0.7016	0.7989	0.3521	-1.0078
	(0.0036)	(0.0034)	(0.0130)	(0.0106)	(0.0727)	(0.0714)	(0.0332)	(0.0392)	(0.0749)	(0.0972)
	[0.0035]	[0.0036]	[0.0116]	[0.0113]	[0.0759]	[0.0324]	[0.0371]	[0.0758]	[0.0697]	[0.1026]
$Laplace$	0.6006	-0.2743	0.2408	-0.6997	1.9983	0.7509	0.7026	0.8034	0.3477	-1.0089
$(0,\frac{\sqrt{2}}{2})$	(0.0030)	(0.0030)	(0.0100)	(0.00995)	(0.0638)	(0.0624)	(0.0256)	(0.0293)	(0.0621)	(0.0873)

Table 2. Table 2 : Model Comparisons: PSAR-ANN model with one neuron ( 6.2 ), PSTAR-ANN models with one ( 6.2 ) and two neurons ( 6.2 )

Models	PSAR-ANN	PSTAR-ANN	PSTAR-ANN
	(one neuron)	(one neuron)	(two neurons)
$\#$ Parameters	$9$	$10$	$14$
Moran’s Test	$0.0745$	$0.2336$	$0.3368$
	$\left(1.7836\right)$	$\left(-1.1910\right)$	$\left(-0.9604\right)$
$-\ln L$	$1879.35$	$1734.5$	$1710.33$
$AIC$	$3776.17$	$3489$	$3448.669$

Table 3. Table 3 : Parameter estimates of PSTAR-ANN model ( 6.2 ) parameters with 95 % percent 95 95\% confidence intervals ( ∗ ∗ \ast indicates the insignificance)

Parameter	Estimate	Std.	$95\%$ C.I.
${\varphi }_0$	0.425	0.0086	$(0.4081,0.4419)$
${\varphi }_1$	0.464	0.0182	$(0.4283,0.4997)$
${\beta }_0$	-1.173	0.3283	$(-1.8165,-0.5295)$
${\beta }_1$	0.148	0.0697	$(0.0114,0.2846)$
${\beta }_2$	-1.177	0.1638	$(-1.4980,-0.8560)$
${\beta }_3$	-0.153	0.1079	${(-0.3645,0.0585)}^{\ast }$
$\lambda$	3.056	0.6397	$(1.8022,4.3098)$
${\gamma }_1$	-0.722	0.1278	$(-0.9725,-0.4715)$
${\gamma }_2$	1.689	0.1762	$(1.3436,2.0344)$
${\gamma }_3$	0.248	0.1890	${(-0.1224,0.6184)}^{\ast }$

Table 4. Table 4 : Parameter estimates of PSTAR-ANN model ( 6.2 ) parameters with 95 % percent 95 95\% confidence intervals ( ∗ ∗ \ast indicates the insignificance)

Parameter	Estimate	Std.	$95\%$ C.I.
${\varphi }_0$	0.417	0.0086	$(0.4000,0.4339)$
${\varphi }_1$	0.467	0.0178	$(0.4321,0.5019)$
${\beta }_0$	-1.576	0.3063	$(-2.1764,-0.9756)$
${\beta }_1$	0.203	0.0731	$(0.0598,0.3462)$
${\beta }_2$	-1.222	0.1507	$(-1.5174,-0.9266)$
${\beta }_3$	-0.057	0.0926	${(-0.2385,0.1245)}^{\ast }$
${\lambda }_1$	3.180	1.2624	$(0.7057,5.6543)$
${\gamma }_{11}$	-0.621	0.1193	$(-0.8548,-0.3872)$
${\gamma }_{12}$	1.621	0.1521	$(1.3230,1.9190)$
${\gamma }_{13}$	0.495	0.1063	$(0.2866,0.7034)$
${\lambda }_2$	0.699	0.2397	$(0.2291,1.1689)$
${\gamma }_{21}$	-1.294	0.6060	$(-2.4368,-0.0612)$
${\gamma }_{22}$	0.084	0.4859	${(-0.8683,1.0363)}^{\ast }$
${\gamma }_{23}$	-3.247	0.3570	$(-3.9469,-2.5471)$

Equations

$$Y_t=\sum _{i=1}^p\sum _{k=0}^{{\lambda }_i}{\varphi }_{ik}{W}^{(k)}{Y}_{t-i}+{\epsilon }_t$$

$$Y_t=\sum _{i=1}^p\sum _{k=0}^{{\lambda }_i}{\varphi }_{ik}{W}^{(k)}{Y}_{t-i}+\sum _{i=0}^{p^{\prime }}{X}_{t-i}{\beta }_i+{\epsilon }_t$$

$$Y_t=\sum _{i=0}^p{\varphi }_i{W}_n{Y}_{t-i}+X_t\beta +\mathbf{F}(X_t{\mathbf{\gamma }}^{\prime })\lambda +{\mathbf{\epsilon }}_t,T=1,\dots ,T$$

$$\left[\begin{array}{cccc}F(x_{1,t}^{\prime }{\mathbf{\gamma }}_1)& F(x_{1,t}^{\prime }{\mathbf{\gamma }}_2)& \dots & F(x_{1,t}^{\prime }{\mathbf{\gamma }}_h)\\ F(x_{2,t}^{\prime }{\mathbf{\gamma }}_1)& F(x_{2,t}^{\prime }{\mathbf{\gamma }}_2)& \dots & F(x_{2,t}^{\prime }{\mathbf{\gamma }}_h)\\ ⋮& ⋮& \cdots & ⋮\\ F(x_{n,t}^{\prime }{\mathbf{\gamma }}_1)& F(x_{n,t}^{\prime }{\mathbf{\gamma }}_2)& \dots & F(x_{n,t}^{\prime }{\mathbf{\gamma }}_h)\end{array}\right]\left[\begin{array}{c}{\lambda }_1\\ {\lambda }_2\\ ⋮\\ {\lambda }_h\end{array}\right]\in {\mathbb{R}}^n$$

$$F(x_{s,t}^{\prime }{\mathbf{\gamma }}_i)=(1+e^{-x_{s,t}^{\prime }{\mathbf{\gamma }}_i}{)}^{-1},s=1,2,\dots ,n,i=1,2,\dots ,h$$

$$Y_t=0.6W_n{Y}_t-0.274W_n{Y}_{t-1}+X_t\left(\begin{array}{c}0.24\\ -0.7\end{array}\right)+1.5\mathbf{F}(X_t{\mathbf{\gamma }}^{\prime })+{\mathbf{\epsilon }}_t,X_t={\left(\begin{array}{ccc}{x}_{11,t}& \dots & x_{1n,t}\\ x_{21,t}& \dots & x_{2n,t}\end{array}\right)}^{\prime }$$

$$A_0=I_n-{\varphi }_0{W}_n,A_i={\varphi }_i{W}_{n}i=1,\dots ,p$$

$$A_0{Y}_t=\sum _{i=1}^p{A}_i{Y}_{t-i}+X_t\beta +\mathbf{F}(X_t{\mathbf{\gamma }}^{\prime })\lambda +{\mathbf{\epsilon }}_t$$

$$Y_t=A_0^{-1}\sum _{i=1}^p{A}_i{Y}_{t-i}+A_0^{-1}{X}_t\beta +A_0^{-1}\mathbf{F}(X_t{\mathbf{\gamma }}^{\prime })\lambda +A_0^{-1}{\mathbf{\epsilon }}_t$$

$$Y_t=A(L{)}^{-1}(X_t\beta +\mathbf{F}(X_t{\mathbf{\gamma }}^{\prime })\lambda +{\mathbf{\epsilon }}_t)$$

$${\Psi }_j={\Theta }_j+\sum _{k=1}^{\infty }{A}_0^{-1}{A}_k{\Psi }_{j-k}$$

$${\Psi }_0$$

$${\Psi }_1$$

$${\Psi }_2$$

$$\cdots$$

$$Y_t=\sum _{j=0}^{\infty }{\Psi }_j{A}_0^{-1}(X_{t-j}\beta +\mathbf{F}(X_{t-j}{\mathbf{\gamma }}^{\prime })+{\mathbf{\epsilon }}_{t-j})$$

$$f_{Y_T,Y_{T-1},\dots ,Y_1}(\mathbf{\theta }|Y_0,\dots ,Y_{1-p},\{X_t\})=$$

$$f_{Y_t}(\mathbf{\theta }|Y_{t-1},\dots ,Y_{1-p},\{X_t\})=$$

$$f_{Y_T,Y_{T-1},\dots ,Y_1}(\mathbf{\theta }|Y_0,\dots ,Y_{1-p},\{X_t\})=$$

$${\mathcal{L}}_{n,T}(\mathbf{\theta })=T\ln |A_0|+\sum _{t=1}^T\sum _{s=1}^n\ln f({\epsilon }_{s,t}(\mathbf{\theta }))$$

$$\det \left[z^p{A}_0-\sum _{i=1}^p{\varphi }_i{W}_n{z}^{p-i}\right]=0$$

$$W=P\left(\begin{array}{cccc}{\tau }_1& 0& \cdots & 0\\ 0& {\tau }_2& \cdots & 0\\ 0& 0& \ddots & 0\\ 0& 0& \cdots & {\tau }_n\end{array}\right)P^{-1},A_0^{-1}=P\left(\begin{array}{cccc}\frac{1}{1-{\varphi }_0{\tau }_1}& 0& \cdots & 0\\ 0& \frac{1}{1-{\varphi }_0{\tau }_2}& \cdots & 0\\ 0& 0& \ddots & 0\\ 0& 0& \cdots & \frac{1}{1-{\varphi }_0{\tau }_n}\end{array}\right)P^{-1}$$

$${\hat{\mathbf{\theta }}}_{n,T}=\arg \underset{\mathbf{\theta }\in \mathbf{\Theta }}{\max }{\mathcal{L}}_{n,T}(\mathbf{\theta })$$

$$\frac{1}{nT}{\mathcal{L}}_{n,T}(\mathbf{\theta })$$

$${\hat{\mathbf{\theta }}}_{n,T}$$

$${\epsilon }_{s,t}(\mathbf{\theta })$$

$${\mathcal{L}}_{n,T}(\mathbf{\theta })=$$

$$l(\mathbf{\theta }|z_{s,t})=$$

$$\underset{n,T\to \infty }{\lim }(nT{)}^{-1}{\mathcal{L}}_{n,T}(\mathbf{\theta })$$

$$=\underset{T\to \infty }{\lim }\frac{1}{T}\sum _{i=1}^T\underset{n_1,n_2\to \infty }{\lim }\frac{1}{n_1{n}_2}\sum _{s\in \mathcal{S}}l(\mathbf{\theta }|z_{s,t})a.s.$$

$$(nT{)}^{-1}{\mathcal{L}}_{n,T}(\mathbf{\theta })\approx \frac{1}{nT}\sum _{i=1}^T\sum _{s\in \mathcal{S}}l(\mathbf{\theta }|z_{s,t})\text{for }{n}_1,n_2,T\to \infty$$

$${\hat{\mathbf{\theta }}}_{n,T}\approx \arg \underset{\mathbf{\theta }\in \mathbf{\Theta }}{\max }\underset{\begin{array}{c}T\to \infty \\ n_1,n_2\to \infty \end{array}}{\lim }\frac{1}{nT}\sum _{i=1}^T\sum _{s\in \mathcal{S}}l(\mathbf{\theta }|z_{s,t})$$

$${\mathbf{\theta }}_0=\underset{\mathbf{\theta }\in \mathbf{\Theta }}{\max }\mathbb{E}{\mathcal{L}}_{n,T}(\mathbf{\theta })\equiv \underset{\mathbf{\theta }\in \mathbf{\Theta }}{\max }\mathbb{E}\frac{1}{nT}{\mathcal{L}}_{n,T}(\mathbf{\theta })$$

$$\mathbb{E}{\mathcal{L}}_{n,T}(\mathbf{\theta })-\mathbb{E}{\mathcal{L}}_{n,T}({\mathbf{\theta }}_0)=\mathbb{E}\ln \frac{L_{n,T}(\mathbf{\theta })}{L_{n,T}({\mathbf{\theta }}_0)}$$

Full text

Partially Specified Space Time Autoregressive Model with Artificial Neural Network

Wenqian Wang

Beth Andrews

Department of Statistics

Northwestern University

1 Introduction

In Chapter 1, we proposed a PSAR model enhanced by a neural network component which aims at explaining the spatial dependence through a nonlinear approach. However, sometimes we may collect data across time as well as space. For this type of data, we want to construct a model with dependence over time taken into consideration which has a broad application especially in environmental sciences. One interesting application is forecasting the weather. For example, in a fixed location, the everyday temperature will change from time to time but in the meanwhile, it would also be affected by temperatures in the neighboring locations.

A class of such linear models known as space-time autoregressive (STAR) and space-time autoregressive moving average (STARMA) models was introduced by Cliff and Ord (1973) and Martin and Oeppen (1975) in 1970s. In general, STAR models contain a hierarchical ordering of “neighbors” of each site. For instance, on a regular grid, one can categorize neighbors of a site as first-order and second-order neighborhoods and so on. An observation at each site is then modeled as a linear function of the previous time observations at the same site and of the weighted previous observations at the neighboring sites of each order. Let $\{Y_{t}:t=0,\pm 1,\pm 2,\ldots\}$ be a multivariate time series of $n$ location components. Weights are incorporated in weight matrices $W^{(k)}$ for order $k$. An STAR model with autoregressive order $p$ and spatial order $(\lambda_{1},\ldots,\lambda_{p})$ considerded in Borovkova et.al (2008) is defined as

[TABLE]

where $\lambda_{i}$ is the spatial order of the $i$th autoregressive term, $\phi_{ik}$ is the autoregressive parameter at time lag $i$ and spatial lag $k$. Similarly an STAR model with $n$ space locations and $q$ exogenous variables is given by Stoffer (1985) as, for $Y_{t}\in\mathbb{R}^{n}$,

[TABLE]

where values of the exogenous variables $\{X_{t}:t=0,\pm 1,\pm 2,\ldots\}$ are $n\times q$ covariate matrices containing $q$ values of exogenous variables for all $n$ locations at time $t$. $X_{t}=(x_{1,t},\ldots,x_{n,t})^{\prime}$ and $x_{s,t}\in\mathbb{R}^{q}$. $p^{\prime}$ is the autoregressive order for $\{X_{t}\}$ and $\beta_{i}$ is a $q\times 1$ model parameter.

STAR models have been widely applied in many areas of science. In genomics, Epperson (1993) analyzed population gene frequencies using STAR models where he assumed genes may vary over space and time. This model is also well known in economics (Giacomini and Granger, 2004) and has been applied to forecasting regional employment (Hernandez and Owyang, 2004) as well as traffic flow (Garrido 2000; Kamarianakis and Prastacos, 2004). For instance, the traffic flow of a road network observed at different fixed locations can be simultaneously modelled as a linear combination of past observations and current observations at neighboring sites. Through weight matrices, an STAR model assumes that near sites exert more influence on each other than distant ones.

In this chapter, we want to extend an STAR model to a semi-parametric model such that this new model can capture nonlinear dependence between covariates and the spatial observations of interest.

2 PSTAR-ANN$(p)$ model

We define a Partially Specified Space-Time Autoregressive model with Artificial Neural Network (PSTAR-ANN$(p)$) as follows.

[TABLE]

where $Y_{t}=\{y_{s,t}\}_{s=1}^{n}$ contains observations of dependent variables at $n$ locations and at time $t$. The independent variable matrix $X_{t}=(x_{1,t},\ldots,x_{n,t})^{\prime}$ is the covariate matrix at time $t$, where $x_{s,t}\in\mathbb{R}^{q\times 1}$ is a vector containing exogenous regressors at location $s$ and time $t$, $s=1,\ldots,n$. $\boldsymbol{\varepsilon}_{t}=\{\varepsilon_{s,t}\}_{s=1}^{n}$ denote a vector of $n$ noise terms which are independent identically distributed across $s$ and $t$ with density function $f$, mean [math] and variance $\sigma^{2}=1$.

Exogenous parameters $\beta=(\beta_{1},\ldots,\beta_{q})^{\prime}\in\mathbb{R}^{q}$ and scalars $\phi_{i},i=0,1,\ldots,p$, the spatial/space-time autoregressive parameters, are assumed to be the same over all regions. $W_{n}=\{w_{ij}\}\in\mathbb{R}^{n\times n}$ is a known spatial weight matrix which characterizes the connection between neighboring regions. For the ease of illustration, we define some notations. Given a function $f\in C^{1}(R^{1})$ continuous in $\mathbb{R}$, we define a new matrix map $R^{n}\rightarrow R^{n}$ as $\boldsymbol{f}$ s.t. $\boldsymbol{f}(x_{1},\ldots,x_{n})=(f(x_{1}),\ldots,f(x_{n}))^{\prime}$.

Using the notation defined above, the artificial neural network component (Medeiros et al. [14]) can be written as $\boldsymbol{F}(X_{t}\boldsymbol{\gamma}^{\prime})\lambda=$

[TABLE]

$\boldsymbol{F}(X_{t}\boldsymbol{\gamma}^{\prime})\lambda$ represents two layer NN component where the first layer has $h$-neurons with the sigmoid activation function and the second layer has only one neuron with an identity activation function. In the first layer, the input is $X_{t}$ and weights are $\boldsymbol{\gamma}=(\boldsymbol{\gamma}^{\prime}_{1},\ldots,\boldsymbol{\gamma}^{\prime}_{h})\in\mathbb{R}^{h\times q}$ where $\boldsymbol{\gamma}_{i}=(\gamma_{i1},\ldots,\gamma_{iq})^{\prime}$ is the weights in the $i$th neuron. $F(\cdot)$ is the sigmoid activation function in this layer.

[TABLE]

In the second layer, the inputs are $F(x_{s,t}^{\prime}\boldsymbol{\gamma}_{i}),i=1,\ldots,h$ and the weights are $\lambda_{1},\ldots,\lambda_{h}$. So final output is $\sum_{i=1}^{h}\lambda_{i}F(x_{s,t}^{\prime}\boldsymbol{\gamma}_{i})$ for each $x_{s,t}$.

The weight matrix $W_{n}$ is a measure of distance between the spatial units, and in our application, we begin by using a square symmetric matrix with $(i,j)$ element equals to 1 if regions $i$ and $j$ are neighbors and [math] otherwise. The diagonal elements of the matrix are set to zero. Then we row standardize this matrix denoted by $W_{n}$. For more details on construction of the weight matrix, you can refer to the previous chapter or LeSage [12]. The following plot provides a preview of the data we are working with. This data is generated from a PSTAR-ANN$(2)$ model in a 10 by 10 lattice. The model equation is shown below:

[TABLE]

$\boldsymbol{\gamma}=(0.75,-0.35)$, with $\{x_{1i,t}\}_{i=1}^{n},\{x_{2i,t}\}_{i=1}^{n}$ are generated i.i.d from $N(0,1.5^{2})$, $N(0,3^{2})$ and the error $\boldsymbol{\varepsilon}_{t}$ is from $N(0,1)$. Figure 1 shows the heatmaps of $Y_{t}$ simulated at $t=30,29,28$ using (2).

The color scale represents the value in each cell. We can observe colors in cells changing gradually with the spatial and time dependence ($\phi_{1}=-0.274$, there is a little flip in cell color comparing the left figure with the middle one).

3 The Model and the Likelihood Function

3.1 The Model

Let

[TABLE]

Suppose $A_{0}$ is invertible, then model (1) can be rewritten as:

[TABLE]

Let $L$ be the usual backshift operator such that $L^{i}Y_{t}=Y_{t-i}$, $A(L)=A_{0}-\sum_{i=1}^{p}A_{i}L^{i}$. Assuming that $A^{-1}(L)$ exists, we can rewrite $Y_{t}$ as

[TABLE]

In order to derive asymptotic properties, we also need $Y_{t}$ to be a causal spatial temporal process. Referring to the definition in Brockwell and Davis [5], the process $Y_{t}$ is causal if there exists matrices $\{\Psi_{j}\}$ with absolutely summable components such that $A^{-1}(L)=\sum_{j=0}^{\infty}\Psi_{j}A_{0}^{-1}L^{j}$. Let $A(z)=A_{0}-A_{1}z-A_{2}z^{2}-\cdots-A_{p}z^{p}=A_{0}(I_{n}-A_{0}^{-1}A_{1}z-A_{0}^{-1}A_{2}z^{2}-\cdots-A_{0}^{-1}A_{p}z^{p})$ be a matrix-valued polynomial. Causality is equivalent to the condition $det(A(z))\neq 0$ for all $z\in\mathbb{C}$ such that $|z|\leq 1$.

The matrices $\Psi_{j}$ can be found recursively from the equations

[TABLE]

where we define $\Theta_{0}=I_{n}$, $\Theta_{j}=0_{n}$ for $j>0$, $A_{j}=0_{n}$ for $j>p$ and $\Psi_{j}=0_{n}$ for $j<0$. Therefore, this gives us

[TABLE]

Then

[TABLE]

With this expansion, we need few assumptions on $\sum_{j=0}^{\infty}\Psi_{j}A_{0}^{-1}$ and will be discussed later.

3.2 Likelihood Function

Denote $\boldsymbol{\theta}=(\phi_{0},\phi_{1},\ldots,\phi_{p},\beta_{1},\ldots,\beta_{q},\lambda,\boldsymbol{\gamma}_{1}^{\prime},\ldots,\boldsymbol{\gamma}_{h}^{\prime})^{\prime}\in\boldsymbol{\Theta}$. Since $\varepsilon_{s,t}$ has an identical density function $f$, the conditional joint density of $Y_{T},Y_{T-1},\ldots,Y_{1}$ conditioned on a finite number of past values $\{Y_{0},\ldots,Y_{1-p}\}$ and $\{X_{t}\}_{t=1}^{T}$ is

[TABLE]

Since

[TABLE]

we have

[TABLE]

Hence, the log-likelihood function of $\boldsymbol{\theta}$ is given by [2, p. 63],

[TABLE]

where $\boldsymbol{\varepsilon}_{t}(\boldsymbol{\theta})=\{\varepsilon_{s,t}(\boldsymbol{\theta})\}_{s=1}^{n}=A(L)Y_{t}-X_{t}\beta-\boldsymbol{F}(X_{t}\boldsymbol{\gamma})\lambda$ for $t=1,\ldots,T$.

For the analysis of identification and estimation of the PSTAR-ANN$(p)$ model, we adopt the following assumptions:

Assumption 1**.**

The $p+(q+1)(h+1)$ parameter vector $\boldsymbol{\theta}=(\phi_{0},\phi_{1},\ldots,\phi_{p},\beta^{\prime},\lambda^{\prime},\boldsymbol{\gamma}_{1}^{\prime},\ldots,\boldsymbol{\gamma}_{h}^{\prime})^{\prime}\in\boldsymbol{\Theta}$, where $\boldsymbol{\Theta}$ is a subset of the $p+(q+1)(h+1)$ dimensional Euclidean space, $\mathbb{R}^{p+(q+1)(h+1)}$. $\boldsymbol{\Theta}$ is a closed and bounded compact set and contains the true parameter value $\boldsymbol{\theta}_{0}$ as an interior point.

Assumption 2**.**

The spatial correlation coefficient $\phi_{0}$ satisfies $|\phi_{0}|<1$ and $\phi_{0}\in(-1/\tau,1/\tau)$, where $\tau=max\{|\tau_{1}|,|\tau_{2}|,\ldots,|\tau_{n}|\}$, $\tau_{1},\ldots,\tau_{n}$ are eigenvalues of spatial weight matrix $W_{n}$. To avoid the non-stationarity issue when $\phi_{0}$ approaches to 1, we assume $\sup_{\phi_{0}\in\boldsymbol{\Theta}}|\phi_{0}|<1$.

Assumption 3**.**

We assume $W_{n}$ is defined by queen contiguity and is uniformly bounded in row and column sums in absolute value as $n\rightarrow\infty$ so $A_{0}^{-1}$ is also uniformly bounded in both column and row sums as $n\rightarrow\infty$.

Assumption 4**.**

We assume a causal spatial process $Y_{t}$ which means that every $z$ which solves

[TABLE]

lie inside a unit circle. So the operator $A(L)$ is causal [16].

Assumption 5**.**

$X_{t}$ is stationary, ergodic satisfying $\mathbb{E}\,|x_{s,t}|^{2}<\infty$ and $X_{t}$ is full column rank for $t=1,2\ldots,T$.

Assumption 6**.**

The error terms $\varepsilon_{s,t}$, $s=1,2,\ldots,n$, $t=1,2\ldots,T$ are independent and identically distributed with density function $f(\cdot)$, zero mean and unit variance $\sigma^{2}=1$. The moment $\mathbb{E}(|\varepsilon_{s,t}|^{2+r})$ exists for some $r>0$ and $\mathbb{E}|\ln f(\varepsilon_{s,t})|<\infty$.

Assumption 2 defines the parameter space for $\phi_{0}$ such that $A_{0}$ is strictly diagonally dominant. By the Levy-Desplanques theorem [19], it follows that $A_{0}^{-1}$ exists for any values $\phi_{0}$ in $(-1/\tau,1/\tau)$. In real applications, since $W_{n}$ is row standardized, one just searches $\hat{\phi}_{0}$ over a parameter space on $(-1,1)$ to find the optimizer [7, p. 749-754].

It is natural to consider the neighborhood by connections and in many practical studies, since entries scaled to sum up to 1, each row of $W_{n}$ sums up to 1, which guarantees that all nonzero weights are in $(0,1]$. For simplicity, we define the weight matrix $W_{n}$ using the queen criterion and do row standardization. Assumption 3 is originated by Kelejian and Prucha [9, 10] and is also used in Lee [11]. With $W_{n}$ to be uniformly bounded, we can prove that $(I_{n}-\phi_{0}W_{n})^{-1}$ is also uniformly bounded in row and column sums for $\phi_{0}\in(-1/\tau,1/\tau)$ and $\sup_{\phi_{0}\in\boldsymbol{\Theta}}|\phi_{0}|<1$, by Lemma A.4 in Lee[11]. This result is a necessary condition for Assumption 4.

From Assumption 2 and 3, we can decompose $W_{n}$ by its eigenvalue and eigenvector pairs $\tau_{i},v_{i}$: $W_{n}=P\Lambda P^{-1}$, where $\Lambda$ is a diagonal matrix with eigenvalues $\tau_{i}$ on its diagonals and $P=[v_{1},v_{2},\ldots,v_{n}]$ (we assume $v_{i}$’s are normalized eigenvectors). So

[TABLE]

It is trivial that $A_{0}^{-1}W_{n}=W_{n}A_{0}^{-1}$.

Assumption 4 guarantees that $A(L)$ is a causal operator and there exists a casual solution $\{Y_{t}\}$ to the system of the model equation (1). Then $\sum_{j=0}^{\infty}\Psi_{j}A_{0}^{-1}$ is absolutely summable. This requirement serves to determine a region of possible $\phi_{i}$ values that will result in a stationary process $\{Y_{t}\}$.

Assumption 5 is a trivial one when exogenous variables are included in a space time model. Similar to previous chapter, the stationarity of $\{x_{s,t}\}$ is necessary in the ergodic theorem in later proofs.

Assumption 6 imposes restrictions for the random error. In this paper we mainly consider the heavy tailed density functions such scaled $t$ distributions and Laplace distributions. When the degrees of freedom goes to infinity, the scaled $t$ distribution would approximate a standard normal distribution. So we would like to concentrate more on the scaled $t$ distribution with lower degrees of freedom.

4 Model Identification

In the previous section, we have some restrictions on the weight matrices $W_{n}$ and $A_{i}$’s to guarantee the identification of a classical spatial time autoregressive model. We now investigate the conditions under which PSTAR(p)-ANN model is identified. By Rothenberg [17], a parameter $\theta_{0}\in\boldsymbol{\Theta}$ is globally identified if there is no other $\theta$ in $\boldsymbol{\Theta}$ that observationally equivalent to $\theta_{0}$ such that $f(y,\theta)=f(y,\theta_{0})$; or the parameter $\theta_{0}$ is locally identified if there is no such $\theta$ in an open neighborhood of $\theta_{0}$ in $\boldsymbol{\Theta}$. The model (1), in principle, is neither globally nor locally identified due to the neural network component. The lack of identification of neural network models has been discussed in many papers (Hwang and Ding [8]; Medeiros et al. [14]). Here we extend the discussion to our proposed PSTAR(p)-ANN model. Three characteristics imply non-identification of our model: (a) the interchangeable property: the value of the likelihood function may remain unchanged if we permute the hidden units. For a model with $h$ neurons, this will result in $h!$ different models that are indistinguishable from each other and have equal local maximums of the log-likelihood function; (b) the “symmetry” property: for a logistic function, $F(x)=1-F(-x)$ allows two equivalent parametrization for each hidden unit; (c) the reducible property: the presence of irrelevant neurons in model (1) happens when $\lambda_{i}=0$ for at least one $i$ and parameters $\boldsymbol{\gamma}_{i}$ remain unidentified. Conversely, if $\boldsymbol{\gamma}_{i}=\mathbf{0}$, $F(X_{t}\boldsymbol{\gamma}_{i})$ is a constant and $\lambda_{i}$ can take any value without affecting the value of likelihood functions.

The problem of interchangeability (as mentioned in (a)) can be solved by imposing the following restriction, as in Medeiros et al. [14]:

Restriction 1. *parameters $\lambda_{1},\ldots,\lambda_{h}$ are restricted such that: $\lambda_{1}\geq\cdots\geq\lambda_{h}$.

*And to tackle (b) and (c), we can apply another restriction:

Restriction 2. *The parameters $\lambda_{i}$ and $\gamma_{i1}$ should satisfy:

(1) $\lambda_{i}\neq 0$, $\forall i\in\{1,2,\ldots,h\}$; and

(2) $\gamma_{i1}>0$, $\forall i\in\{1,2,\ldots,h\}$.

*To guarantee the non-singularity of model matrices and the uniqueness of parameters, we impose the following basic assumption:

Assumption 7**.**

The true parameter vector $\boldsymbol{\theta}_{0}$ satisfies Restrictions 1-2.

Referring to the section 4.3 by Medeiros et al. [14], we can conclude the identifiability of the PSAR-ANN model.

Lemma 1**.**

Under the Assumptions 1-7, this PSTAR-ANN$(p)$ model (1) is globally identified.

5 Asymptotic Results

Let the true parameter vector as $\boldsymbol{\theta}_{0}$ and the solution which maximizes the log-likelihood function (6) as $\hat{\boldsymbol{\theta}}_{n,T}$. Hence, $\hat{\boldsymbol{\theta}}_{n,T}$ should satisfy

[TABLE]

Suppose as $n$ is large enough, $T$ goes to infinity, $\hat{\boldsymbol{\theta}}_{n,T}$ is equivalent to maximizing the average of the likelihood function $\mathcal{L}_{n,T}(\boldsymbol{\theta})$ shown as follows:

[TABLE]

At specific time $t$, suppose we have a $n_{1}\times n_{2}$ lattice where we consider asymptotic properties of $\hat{\boldsymbol{\theta}}_{n,T}$ when $n=n_{1}n_{2}\rightarrow\infty$. Write the location $s$ as the coordinate $(s_{x},s_{y})$ in the $[1,n_{1}]\times[1,n_{2}]$ lattice space. The distance between two locations $s,j$ is defined as $d(s,j)=\max(|s_{x}-j_{x}|,|s_{y}-j_{y}|)$. So if observations at $s,j$ locations are neighbors (by queen criterion), their coordinates should satisfy $(s_{x}-j_{x})^{2}+(s_{y}-j_{y})^{2}\leq 2$ or $d(s,j)=1$.

In a spatial context, we should notice that the functional form of $y_{s,t}$ is not identical for all the locations due to values of the weights $\{w_{si}\}_{i=1}^{n}$. For example, in a lattice, units at edges, vertexes or in the interior have different density functions due to different neighborhood structures (Figure 2). Denote $\mathcal{N}_{s}$ as a neighborhood set for location $s$. For an interior point (Figure 2(c)), its neighborhood set $\mathcal{N}_{s}$ contains eight neighbors where $w_{sj}=1/8$ if $d(s,j)=1$ otherwise $w_{sj}=0$, for $j=1,2,\ldots,n$. Similarly, an edge point (Figure 2(b)) has five neighboring units with $w_{sj}=1/5$ for $j\in\mathcal{N}_{s}$ and the weight of a vertex neighborhood is $1/3$ because a vertex unit has only three neighbors. This is known as an edge effect in spatial problems.

To deal with this, referring to Yao and Brockwell [23], we construct an edge effect correction scheme based on the way that the sample size tends to infinity. In a space $[1,n_{1}]\times[1,n_{2}]$, we consider its interior area as $\mathcal{S}=\{(s_{x},s_{y}):b_{1}\leq s_{x}\leq n_{1}-b_{1},b_{2}\leq s_{y}\leq n_{2}-b_{2}\}$, where $b_{1},b_{2},n_{1},n_{2}\rightarrow\infty$ satisfying that $b_{1}/n_{1},b_{2}/n_{2}\rightarrow 0$ and other locations belong to the boundary areas $\mathcal{M}$. Therefore the set $\mathcal{S}$ contains $n^{\ast}=(n_{1}-2b_{1})(n_{2}-2b_{2})$ interior locations while the set $\mathcal{M}$ contains $n-n^{\ast}$ boundary locations. Then $n^{\ast}/n\rightarrow 1$ and $\mathcal{L}_{n,T}(\boldsymbol{\theta})$ can be split into a sum of two parts (interior $\mathcal{S}$ and boundary $\mathcal{M}$ parts):

[TABLE]

where $Z_{t}=(W_{n}Y_{t},W_{n}Y_{t-1},\ldots,W_{n}Y_{t-p},X_{t})$ and $z_{s,t}$ is the $s$ row of $Z_{t}$.

Therefore, given that $\lim_{n_{1},n_{2}\rightarrow\infty}\frac{|\mathcal{M}|}{n}=0$, $n^{-1}\sum_{s\in\mathcal{M}}l(\boldsymbol{\theta}|z_{s,t})$ vanishes a.s. as $n$ tends to infinity for any $\boldsymbol{\theta}\in\boldsymbol{\Theta}$. Therefore,

[TABLE]

In this equation, every location $s\in\mathcal{S}$ has eight neighboring units under the queen criterion with nonzero weights $w_{sj}=1/8$. Hence for an interior unit $s\in\mathcal{S}$, $\sum_{i=1}^{n}w_{si}y_{i}=\sum_{j=1}^{n}\frac{1}{8}y_{j}I_{\{d(s,j)=1\}}$. And the log likelihood function $\mathcal{L}_{n,T}(\boldsymbol{\theta})$ is approximately

[TABLE]

So the maximum likelihood estimator $\hat{\boldsymbol{\theta}}_{n,T}$ approximately maximizes

[TABLE]

5.1 Consistency Results

To establish the consistency of $\hat{\boldsymbol{\theta}}_{n,T}$, the heuristic insight is that because $\hat{\boldsymbol{\theta}}_{n,T}$ maximizes $\frac{1}{nT}\mathcal{L}_{n,T}(\boldsymbol{\theta})$, it approximately maximizes $\frac{1}{nT}\sum_{i=1}^{T}\sum_{s\in\mathcal{S}}l(\boldsymbol{\theta}|z_{s,t})$. By equation (8), $\frac{1}{nT}\mathcal{L}_{n,T}(\boldsymbol{\theta})$ can generally be shown tending to a real function $\mathcal{L}:\boldsymbol{\Theta}\rightarrow\mathbb{R}$ with maximizer $\boldsymbol{\theta}_{0}$ as $n,T\rightarrow\infty$ under mild conditions on the data generating process, then $\hat{\boldsymbol{\theta}}_{n,T}$ should tend to $\boldsymbol{\theta}_{0}$ almost surely. Before the formal proof of the consistency, we need the following assumptions on the density function $f(\cdot)$ satisfied (similar assumptions are made in White [22], Andrews, Davis and Breidt [1], Lii and Rosenblatt [13]).

Assumption 8**.**

For all $s\in\mathbb{R}$, $f(s)>0$ and $f(s)$ is twice continuously differentiable with respect to $s$.

Assumption 9**.**

The density should satisfy the following equations:

$\bullet$ $\int sf^{\prime}(s)\,ds=sf(s)|^{\infty}_{-\infty}-\int f(s)\,ds=-1$

$\bullet$ $\int f^{\prime\prime}(s)\,ds=f^{\prime}(s)|^{\infty}_{-\infty}=0$

$\bullet$ $\int s^{2}f^{\prime\prime}(s)\,ds=s^{2}f^{\prime}(s)|^{\infty}_{-\infty}-2\int sf^{\prime}(s)\,ds=2$

Assumption 10**.**

The density should follow the following dominance conditions:

$\left|\frac{f^{\prime}(s)}{f(s)}\right|$, $\left|\frac{f^{\prime}(s)}{f^{(}s)}\right|^{2}$, $\left|\frac{f^{\prime}(s)}{f(s)}\right|^{4}$, $\frac{f^{\prime\prime}(s)}{f(s)}$, and $\frac{f^{\prime\prime}(s)f^{\prime 2}(s)}{f^{3}(s)}$ are dominated by $a_{1}+a_{2}\left|s\right|^{c_{1}}$, where $a_{1}$, $a_{2}$, $c_{1}$ are non-negative constants and $\int\left|s\right|^{c_{1}+2}f(s)\,ds<\infty$.

Assumption 11**.**

If $c_{1}>2$ in previous assumption, we further assume $\mathbb{E}\,|x_{s,t}|^{c_{1}}<\infty$.

Discussed in Breidt, Davis, Lii and Rosenblatt [4] and Andrews, Davis and Breidt [1, p. 1642-1645], these assumptions on the density $f(\cdot)$ are satisfied by the t-distribution case when $\nu>2$ and by a mixture of Gaussian distributions. The assumption $\mathbb{E}|\ln f(s)|<\infty$ (see Assumption 6) is also checked satisfied by the normal and t distributions ($\nu>2$). The Laplace distribution does not strictly satisfy the Assumptions 8-10, since it is not differentiable at 0 but it satisfies these boundedness conditions almost everywhere so we believe the consistency and asymptotic normality results remain valid for parameter estimates. This will be shown in the simulation section. Assumption 11 is a necessary to boundedness conditions in later proof.

Lemma 2**.**

Given Assumptions 1-10,

[TABLE]

Proof.

$L_{n,T}$ is the joint density function of $Y_{t},X_{t}$ for $t=1,\ldots,T$.

[TABLE]

Denote $Z=(Y_{T},X_{T},\ldots,Y_{1},X_{1})$. By Jensen’s inequality,

[TABLE]

So $\mathbb{E}\,\mathcal{L}_{n,T}(\boldsymbol{\theta})<\mathbb{E}\,\mathcal{L}_{n,T}(\boldsymbol{\theta}_{0})$. By Lemma 1, the PSTAR(p)-ANN model is globally identified and therefore $\mathbb{E}\,\mathcal{L}_{n,T}(\boldsymbol{\theta})$ is uniquely maximized at $\boldsymbol{\theta}_{0}$ for all $n,T$. Since the parameter vector $\boldsymbol{\theta}$ does not depend on $n$ and $T$, it is equivalent to say that $\boldsymbol{\theta}_{0}=\max_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\mathbb{E}\,\frac{1}{nT}\mathcal{L}_{n,T}(\boldsymbol{\theta})$. ∎

We define a Hadamard product denoted by $\circ$, s.t. for vectors $a,b_{1},\ldots,b_{n}\in\mathbb{R}^{n}$, a matrix $B=(b_{1},\ldots,b_{n})\in\mathbb{R}^{n\times n}$,

[TABLE]

And let

[TABLE]

To facilitate the proof later on, we provide a lemma as follows.

Lemma 3**.**

Given Assumptions 1-11,

[TABLE]

Proof.

As illustrated in equation (8), in a lattice with size $n_{1}\times n_{2}$,

[TABLE]

Therefore, to prove (9) is equivalent to show that

[TABLE]

where $\mathcal{S}$ denotes the interior units mentioned before. Since the interior units have the same neighboring structure, the space process for them is stationary when $n_{1},n_{2}$ go to infinity. We first show $\left|\frac{1}{nT}\sum_{t=1}^{T}\left(\sum_{s\in\mathcal{S}}\ln f(\varepsilon_{s}(\boldsymbol{\theta}))-\mathbb{E}\,\frac{1}{n}\sum_{s\in\mathcal{S}}\ln f(\varepsilon_{s}(\boldsymbol{\theta}))\right)\right|\xrightarrow{p}0$ for fixed $\boldsymbol{\theta}$.

To prove this, we want to show that $\mathbb{E}\,|\ln f(\varepsilon_{s,t}(\boldsymbol{\theta}))|<\infty$. Expanding $\ln f(\varepsilon_{s,t}(\boldsymbol{\theta}))$ around $\boldsymbol{\theta}_{0}$ with respect to $\boldsymbol{\theta}$,

[TABLE]

where $\tilde{\boldsymbol{\theta}}_{n,T}$ is between $\boldsymbol{\theta}$ and $\boldsymbol{\theta}_{0}$. Under the true parameter values, $\varepsilon_{s,t}(\boldsymbol{\theta}_{0})$ (denoted as $\varepsilon_{s,t}$ or $\boldsymbol{\varepsilon}_{t}$ as its vector form in the following) is independent and identically distributed. From Assumption 6, $\mathbb{E}\,|\ln f(\varepsilon_{s,t})|<\infty$. For $\mathbb{E}\,\left|\frac{f^{\prime}(\varepsilon_{s,t}(\tilde{\boldsymbol{\theta}}_{n,T}))}{f(\varepsilon_{s,t}(\tilde{\boldsymbol{\theta}}_{n,T}))}\frac{\partial\varepsilon_{s,t}(\tilde{\boldsymbol{\theta}}_{n,T})}{\partial\boldsymbol{\theta}^{\prime}}\right|$, $\left|\frac{\partial\varepsilon_{s,t}(\tilde{\boldsymbol{\theta}}_{n,T})}{\partial\boldsymbol{\theta}}\right|$ can be expressed as

[TABLE]

where $A^{-1}(L)(\boldsymbol{g}(X_{t-i},\boldsymbol{\theta}_{0})+\boldsymbol{\varepsilon}_{t-i}(\boldsymbol{\theta}_{0})=\sum_{j=0}^{\infty}\Psi_{j}A_{0}^{-1}(\boldsymbol{g}(X_{t-i-j},\boldsymbol{\theta}_{0})+\boldsymbol{\varepsilon}_{t-i-j}(\boldsymbol{\theta}_{0}))$. Function $g(x_{s,t},\boldsymbol{\theta})=x_{s,t}^{\prime}\beta+\boldsymbol{F}(x_{s,t}^{\prime}\gamma)\lambda$. Consider $\boldsymbol{\varepsilon}_{t}(\tilde{\boldsymbol{\theta}}_{n,T})$,

[TABLE]

Denote $P(x^{c})$ is a polynomial about $x$ with highest order $c$. Since we have assumed that $A^{-1}(L)$ existed and the expansion $\sum_{j=0}^{\infty}\Psi_{j}A_{0}^{-1}$ is absolutely summable so $W_{n}\sum_{j=0}^{\infty}\Psi_{j}A_{0}^{-1}$ is finite. By Assumption 10-11, $\left|\frac{f^{\prime}(\varepsilon_{s,t}(\tilde{\boldsymbol{\theta}}_{n,T}))}{f(\varepsilon_{s,t}(\tilde{\boldsymbol{\theta}}_{n,T}))}\right|<a_{1}+a_{2}|\varepsilon_{s,t}(\tilde{\boldsymbol{\theta}}_{n,T})|^{c_{1}}$ and $\mathbb{E}\left|\frac{f^{\prime}(\varepsilon_{s,t})}{f(\varepsilon_{s,t})}\right|,\mathbb{E}\left|\frac{f^{\prime}(\varepsilon_{s,t})}{f(\varepsilon_{s,t})}\right|^{2}$ are dominated by $a_{1}+a_{2}|\varepsilon_{s,t}|^{c_{1}}$, $\mathbb{E}|\varepsilon_{s,t}|^{c_{1}}<\infty,\,\mathbb{E}|x_{s,t}|^{c_{1}}<\infty$. Let $c^{\ast}=\max(1,c_{1})$, then,

[TABLE]

So also $\mathbb{E}\left|\frac{f^{\prime}(\varepsilon_{s,t}(\tilde{\boldsymbol{\theta}}_{n,T}))}{f(\varepsilon_{s,t}(\tilde{\boldsymbol{\theta}}_{n,T}))}\right|<\infty$. With Cauchy–Schwarz inequality [18] and the finite second moment of $x_{s,t}$, we can have,

[TABLE]

Because $W_{n}A^{-1}(L)$ is well defined and $X_{t}$ is stationary with finite second moment, so component (17) is finite. (16) is dominated by $P(\mathbb{E}\,|\varepsilon_{s,t}|^{c^{\ast}+1})$ so with the dominance assumption, (16) is finite. Hence, with (12)-(17) finite, $\mathbb{E}\,|\ln f(\varepsilon_{s,t}(\boldsymbol{\theta}_{0}))|<\infty$, we can conclude that $\mathbb{E}\,|\ln f(\varepsilon_{s,t}(\boldsymbol{\theta}))|<\infty$. Then by ergodic theorem [3],

[TABLE]

To complete the proof of uniform convergence, we also need to show $\frac{1}{nT}\sum_{t=1}^{T}\sum_{s\in\mathcal{S}}\ln f(\varepsilon_{s,t}(\boldsymbol{\theta}))$ is equicontinuous for $\boldsymbol{\theta}\in\boldsymbol{\Theta}$, i.e., for all $\boldsymbol{\theta}_{1},\boldsymbol{\theta}_{2}\in\boldsymbol{\Theta}$,

[TABLE]

Applying the mean value theorem to the left side in (18):

[TABLE]

where $\tilde{\boldsymbol{\theta}}_{n,T}$ is some value between $\boldsymbol{\theta}_{1}$ and $\boldsymbol{\theta}_{2}$. Since $\boldsymbol{\theta}$ is in a compact set $\boldsymbol{\Theta}$, we show in (5.1) that, for all $s,t$, $\varepsilon_{s,t}(\boldsymbol{\theta})$ is bounded by some function of $Z_{t}$ not depending on $\boldsymbol{\theta}$.

[TABLE]

Similarly, referring to (5.1), it is easy to show that $\left|\frac{\partial\varepsilon_{s,t}(\boldsymbol{\theta})}{\partial\boldsymbol{\theta}}\right|$ is bounded by some function about $Y_{t}$ and $X_{t}$. Therefore, due to the dominance of $\left|\frac{f^{\prime}(s)}{f(s)}\right|$ (see Assumption 10) and stationarity of $X_{t},Y_{t}$, for $\tilde{\boldsymbol{\theta}}_{n,T}$ between $\boldsymbol{\theta}_{1}$ and $\boldsymbol{\theta}_{2}$, there exists a constant $M$ such that

[TABLE]

Hence, for $\boldsymbol{\theta}_{1},\boldsymbol{\theta}_{2}\in\boldsymbol{\Theta}$

[TABLE]

So $\frac{1}{nT}\sum_{t=1}^{T}\sum_{s\in\mathcal{S}}\ln f(\varepsilon_{s,t}(\boldsymbol{\theta}))$ is equicontinuous for $\boldsymbol{\theta}\in\boldsymbol{\Theta}$. With the pointwise convergence and equicontinuity, we can conclude the uniform convergence in (10) and furthermore (9) follows. ∎

Similar to Chapter 1, we now give a formal statement of the consistency results.

Theorem 1**.**

Given Assumptions 1-11, $\hat{\boldsymbol{\theta}}_{n,T}\xrightarrow{p}\boldsymbol{\theta}_{0}$ as $n,T\rightarrow\infty$.

Proof.

Similar to the proof by Lung-fei Lee [11], we need to show the stochastic equicontinuity of $\frac{1}{n}\ln|A_{0}|$ to have the uniform convergence of the log likelihood function $\mathcal{L}_{n,T}(\boldsymbol{\theta})$. Applying the mean value theorem,

[TABLE]

where $\phi_{0_{n,T}}^{\ast}$ is between $\phi_{0}^{\dagger}$ and $\phi_{0}^{\ddagger}$. By Assumption 2 and 3, $\sup_{\phi_{0}\in\boldsymbol{\Theta}}|\phi_{0}|<1$, $W_{n}$ is bounded in both rows and column sums uniformly and using (7),

[TABLE]

where $C_{1}$ is a constant not depending on $n$. So $\left|\frac{1}{n}(\ln|I_{n}-\phi_{0}^{\dagger}W_{n}|-\ln|I_{n}-\phi_{0}^{\ddagger}W_{n}|)\right|\leq C_{1}|\phi_{0}^{\dagger}-\phi_{0}^{\ddagger}|$ and with Lemma 3 we can conclude the uniform convergence that

[TABLE]

With the assumptions 1-10, the parameter space $\boldsymbol{\Theta}$ is compact; $\frac{1}{nT}\mathcal{L}_{n,T}(\boldsymbol{\theta})$ is continuous in $\boldsymbol{\theta}\in\boldsymbol{\Theta}$ and is a measurable of $Y_{t},X_{t},t=1,\ldots,T$ for all $\boldsymbol{\theta}\in\boldsymbol{\Theta}$. $\mathbb{E}\,\frac{1}{nT}\mathcal{L}_{n,T}(\boldsymbol{\theta})$ is continuous on $\boldsymbol{\Theta}$ and by Lemma 2, $\mathbb{E}\,\frac{1}{nT}\mathcal{L}_{n,T}(\boldsymbol{\theta})$ has a unique maximum at $\boldsymbol{\theta}_{0}$. Referring to Theorem 3.5 in White [21] with the uniform convergence in (9), we can conclude that $\hat{\boldsymbol{\theta}}_{n,T}\xrightarrow{p}\boldsymbol{\theta}_{0}$ as $n,T\rightarrow\infty$. ∎

5.2 Asymptotic Distribution

Assumption 12**.**

The limit $A(\boldsymbol{\theta}_{0})=-\lim_{n,T\rightarrow\infty}\mathbb{E}\,\frac{1}{nT}\frac{\partial^{2}\mathcal{L}_{n,T}(\boldsymbol{\theta}_{0})}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{\prime}}$ is nonsingular.

Assumption 13**.**

The limit $B(\boldsymbol{\theta}_{0})=\lim_{n,T\rightarrow\infty}\mathbb{E}\,\frac{1}{nT}\frac{\partial\mathcal{L}_{n,T}(\boldsymbol{\theta}_{0})}{\partial\boldsymbol{\theta}}\frac{\partial\mathcal{L}_{n,T}(\boldsymbol{\theta}_{0})}{\partial\boldsymbol{\theta}^{\prime}}$ is nonsingular.

These assumptions are to guarantee the existence of the covariance matrix of the limiting distribution of parameters in a PSTAR(p)-ANN model. We now give the asymptotic distribution of the maximum likelihood estimator $\hat{\boldsymbol{\theta}}_{n,T}$.

Theorem 2**.**

Under Assumptions 1-13,

[TABLE]

where $\boldsymbol{\Omega}_{0}=A(\boldsymbol{\theta}_{0})^{-1}B(\boldsymbol{\theta}_{0})A(\boldsymbol{\theta}_{0})^{-1}=A(\boldsymbol{\theta}_{0})^{-1}$

Proof.

Since $\hat{\boldsymbol{\theta}}_{n,T}$ maximizes $\mathcal{L}_{n,T}(\boldsymbol{\theta})$, $\frac{\partial\mathcal{L}_{n,T}(\hat{\boldsymbol{\theta}}_{n,T})}{\partial\boldsymbol{\theta}}=0$. By the mean value theorem, expand $\frac{\partial\mathcal{L}_{n,T}(\hat{\boldsymbol{\theta}}_{n,T})}{\partial\boldsymbol{\theta}}$ around $\boldsymbol{\theta}_{0}$ with respect to $\boldsymbol{\theta}$,

[TABLE]

where $\tilde{\boldsymbol{\theta}}_{n,T}$ is between $\hat{\boldsymbol{\theta}}_{n,T}$ and $\boldsymbol{\theta}_{0}$. Therefore, we can have the following equation:

[TABLE]

From (5.1), denote $\frac{\boldsymbol{f}^{\prime}(\varepsilon_{t}\boldsymbol{\theta})}{\boldsymbol{f}(\varepsilon_{t}(\boldsymbol{\theta}))}$ as $V_{t}(\boldsymbol{\theta})\in\mathbb{R}^{n}$ and $\frac{\boldsymbol{f}^{\prime}(\varepsilon_{t}\boldsymbol{\theta}_{0})}{\boldsymbol{f}(\varepsilon_{t}(\boldsymbol{\theta}_{0}))}=V_{t}$.

Recall that $Z_{t}=(W_{n}Y_{t},W_{n}Y_{t-1},\ldots,W_{n}Y_{t-p},X_{t})$ so the first order derivatives can be expressed as

[TABLE]

By Lemma 2, the true parameter values maximize $\frac{1}{nT}\mathbb{E}\,\mathcal{L}_{n,T}(\boldsymbol{\theta})$, so $\frac{1}{nT}\frac{\partial\mathbb{E}\,\mathcal{L}_{n,T}(\boldsymbol{\theta}_{0})}{\partial\boldsymbol{\theta}}=\boldsymbol{0}$. In (12)-(17) and (5.1), we showed that $\mathbb{E}\left|\frac{\partial\ln f(\varepsilon_{s,t}(\boldsymbol{\theta}))}{\partial\boldsymbol{\theta}}\right|$ is dominated by some function not related to $\boldsymbol{\theta}$ and (20) indicates that $\mathbb{E}\left|\frac{\partial\ln f(\varepsilon_{s,t}(\boldsymbol{\theta}))}{\partial\boldsymbol{\theta}}\right|$ is bounded for interior units in $\mathcal{S}$. Hence,$\mathbb{E}\,\frac{\partial\ln f(\varepsilon_{s,t}(\boldsymbol{\theta}))}{\partial\boldsymbol{\theta}}=\frac{\partial}{\partial\boldsymbol{\theta}}\mathbb{E}\ln f(\varepsilon_{s,t}(\boldsymbol{\theta}))$, it follows that, with $\frac{1}{nT}\mathcal{L}_{n,T}(\boldsymbol{\theta})=\frac{1}{n}\ln|A_{0}|+\frac{1}{nT}\sum_{s=1}^{n}\sum_{t=1}^{T}\ln f(\varepsilon_{s,t}(\boldsymbol{\theta}))$, we can have,

[TABLE]

Therefore, with Assumption 13,

[TABLE]

And under this $A(\boldsymbol{\theta}_{0})=B(\boldsymbol{\theta}_{0})$. Since $\frac{\partial\mathcal{L}_{n,T}(\boldsymbol{\theta}_{0})}{\partial\boldsymbol{\theta}}$ is the sum of $T$ identical and ergodic random variables, by the central limit theorem for stationary ergodic processes [15], the limiting distribution of $\frac{1}{\sqrt{nT}}\frac{\partial\mathcal{L}_{n,T}(\boldsymbol{\theta}_{0})}{\partial\boldsymbol{\theta}}$ is $N(\boldsymbol{0},B(\boldsymbol{\theta}_{0}))$.

Next we would like to show that $\frac{1}{nT}\frac{\partial^{2}\mathcal{L}_{n,T}(\tilde{\boldsymbol{\theta}}_{n,T})}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{\prime}}-\frac{1}{nT}\frac{\partial^{2}\mathcal{L}_{n,T}(\boldsymbol{\theta}_{0})}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{\prime}}\xrightarrow{p}0$. Following the results in (23), define $U_{t}(\boldsymbol{\theta})=\frac{\boldsymbol{f}^{\prime\prime}(\varepsilon_{t}(\boldsymbol{\theta}))}{\boldsymbol{f}(\varepsilon_{t}(\boldsymbol{\theta}))}-\frac{\boldsymbol{f}^{\prime 2}(\varepsilon_{t}(\boldsymbol{\theta}))}{\boldsymbol{f}^{2}(\varepsilon_{t}(\boldsymbol{\theta}))}\in\mathbb{R}^{n}$, and write $U_{t}=U_{t}(\boldsymbol{\theta}_{0})$ so the second order derivatives are given below $-\frac{1}{nT}\frac{\partial^{2}\mathcal{L}_{n,T}(\boldsymbol{\theta})}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{\prime}}=$

[TABLE]

[TABLE]

Since $\tilde{\boldsymbol{\theta}}_{n,T}$ is between $\hat{\boldsymbol{\theta}}_{n,T}$ and $\boldsymbol{\theta}_{0}$, $\hat{\boldsymbol{\theta}}_{n,T}\xrightarrow{p}\boldsymbol{\theta}_{0}$ so $\tilde{\boldsymbol{\theta}}_{n,T}$ also converges to $\boldsymbol{\theta}_{0}$ in probability as $n\rightarrow\infty$. By Assumption 10, $\left|\frac{f^{\prime}(s)}{f(s)}\right|,\left|\frac{f^{\prime\prime}(s)}{f(s)}\right|$ and $\left|\frac{f^{\prime 2}(s)}{f^{2}(s)}\right|$ are continuous and are bounded by $a_{1}+a_{2}\left|s\right|^{c_{1}}$ so $U_{t}(\boldsymbol{\theta}),V_{t}(\boldsymbol{\theta})$ are continuous. With $\phi_{0}\in(-\frac{1}{\tau},\frac{1}{\tau})$, $tr((W_{n}A_{0}^{-1})^{2})=\sum_{i=1}^{n}\frac{\tau_{i}^{2}}{(1-\phi_{0}\tau_{i})^{2}}$ is also a continuous function of $\phi_{0}$.

Therefore elements in $\frac{1}{nT}\frac{\partial^{2}\mathcal{L}_{n,T}(\boldsymbol{\theta})}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{\prime}}$ are continuous functions for $\boldsymbol{\theta}$ in $\boldsymbol{\Theta}$. Then by the continuity,

[TABLE]

Finally we will prove that $\left|\frac{1}{nT}\frac{\partial^{2}\mathcal{L}_{n,T}(\boldsymbol{\theta}_{0})}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{\prime}}-\mathbb{E}\frac{1}{nT}\frac{\partial^{2}\mathcal{L}_{n,T}(\boldsymbol{\theta}_{0})}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{\prime}}\right|\xrightarrow{p}0$. Since $\ln|A_{0}|$ can be decomposed as $\sum_{i=1}^{n}\ln(1-\phi_{0}\tau_{i})$, to show $\mathbb{E}\frac{1}{nT}\frac{\partial^{2}\mathcal{L}_{n,T}(\boldsymbol{\theta}_{0})}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{\prime}}<\infty$ is equivalent to show

[TABLE]

We first discuss the second derivative with respect to $\phi_{0}$ component in (26). By triangular inequality, \mathbb{E}\left|\frac{\partial^{2}}{\partial\phi_{0}\partial\phi_{0}}\frac{1}{nT}\sum_{t=1}^{T}\sum_{s=1}^{n}\Big{(}\ln(1-\phi_{00}\tau_{s})+\ln f(\varepsilon_{s,t}(\boldsymbol{\theta}_{0}))\Big{)}\right|<\mathbb{E}\left|\frac{1}{n}\sum_{s=1}^{n}\frac{\partial^{2}\ln(1-\phi_{00}\tau_{s})}{\partial\phi_{0}\partial\phi_{0}}\right|+\mathbb{E}\left|\frac{\partial^{2}\ln f(\varepsilon_{s,t}(\boldsymbol{\theta}_{0}))}{\partial\phi_{0}\partial\phi_{0}}\right| where $\phi_{00}$ is the true value of $\phi_{0}$. Consider $\mathbb{E}\left|\frac{1}{n}\sum_{s=1}^{n}\frac{\partial^{2}\ln(1-\phi_{00}\tau_{i})}{\partial\phi_{0}\partial\phi_{0}}\right|+\mathbb{E}\left|\frac{\partial^{2}\ln f(\varepsilon_{s,t}(\boldsymbol{\theta}_{0}))}{\partial\phi_{0}\partial\phi_{0}}\right|$, under stationarity, it can be simplified as

[TABLE]

Define $M_{n}=\{m_{i,j}\}=W_{n}A_{0}^{-1}$ and by assumptions, $M_{n}$ is uniformly bounded in row and column. Suppose the row sum or column sum of $M_{n}$ is bounded by a constant $b$. We know $\frac{1}{n}tr(W_{n}A_{0}^{-1})^{2}<\infty$. So we only need to show $\mathbb{E}\left|\left(\frac{f^{\prime^{2}}(\varepsilon_{s,t})}{f^{2}(\varepsilon_{s,t})}-\frac{f^{\prime\prime}(\varepsilon_{s,t})}{f(\varepsilon_{s,t})}\right)\left(\sum_{k=1}^{n}w_{sk}y_{k,t}\right)^{2}\right|<\infty$.

By simple linear algebra,

[TABLE]

So $W_{n}Y_{t}=W_{n}\sum_{j=0}^{\infty}\Psi_{j}A_{0}^{-1}(\boldsymbol{g}(X_{t-j},\boldsymbol{\theta}_{0})+\boldsymbol{\varepsilon}_{t-j})$. Therefore $(\sum_{k=1}^{n}w_{sk}y_{k,t})^{2}$ is the $s^{th}$ component of $(W_{n}Y_{t}\circ W_{n}Y_{t})$ and we expand $(W_{n}Y_{t}\circ W_{n}Y_{t})_{s}=$

[TABLE]

From assumptions 3 and 4, we know that $W_{n}$ is uniformly bounded and $\sum_{j=0}^{\infty}\Psi_{j}A_{0}^{-1}$ is absolute summable so $(\ref{sq-y-1})<\infty$ under the stationary condition of $X_{t}$. Hence, $\mathbb{E}\left|\left(\frac{f^{\prime^{2}}(\varepsilon_{s,t})}{f^{2}(\varepsilon_{s,t})}-\frac{f^{\prime\prime}(\varepsilon_{s,t})}{f(\varepsilon_{s,t})}\right)\cdot\text{(\ref{sq-y-1})}\right|<\infty$.

For (29), when $j>0$, $\boldsymbol{\varepsilon}_{t-j}$ is independent from $\boldsymbol{\varepsilon}_{t}$. So for all $k$ when $j>0$,

[TABLE]

when $j=0$, this reduces to $\mathbb{E}\left|\left(\frac{f^{\prime^{2}}(\varepsilon_{s,t})}{f^{2}(\varepsilon_{s,t})}-\frac{f^{\prime\prime}(\varepsilon_{s,t})}{f(\varepsilon_{s,t})}\right)\cdot\left[W_{n}A_{0}^{-1}\boldsymbol{\varepsilon}_{t}\right]_{s}\right|<k_{1}|b-m_{ss}|+k_{2}|m_{ss}|$. So $\mathbb{E}\left|\left(\frac{f^{\prime^{2}}(\varepsilon_{s,t})}{f^{2}(\varepsilon_{s,t})}-\frac{f^{\prime\prime}(\varepsilon_{s,t})}{f(\varepsilon_{s,t})}\right)\cdot\text{(\ref{sq-y-2})}\right|<\infty$.

For (30), similar to (29), we can have $\mathbb{E}\left|\left(\frac{f^{\prime^{2}}(\varepsilon_{s,t})}{f^{2}(\varepsilon_{s,t})}-\frac{f^{\prime\prime}(\varepsilon_{s,t})}{f(\varepsilon_{s,t})}\right)\cdot\text{(\ref{sq-y-3})}\right|<Constant\cdot(k_{2}+k_{3}+\mathbb{E}|\varepsilon_{s,t}|)<\infty$.

Therefore combining all these components together, \mathbb{E}\left|\left(\frac{f^{\prime^{2}}(\varepsilon_{s,t})}{f^{2}(\varepsilon_{s,t})}-\frac{f^{\prime\prime}(\varepsilon_{s,t})}{f(\varepsilon_{s,t})}\right)\left(\sum_{k=1}^{n}w_{sk}y_{k,t}\right)^{2}\right|=\mathbb{E}\left|\left(\frac{f^{\prime^{2}}(\varepsilon_{s,t})}{f^{2}(\varepsilon_{s,t})}-\frac{f^{\prime\prime}(\varepsilon_{s,t})}{f(\varepsilon_{s,t})}\right)\Big{(}(\ref{sq-y-1})+(\ref{sq-y-2})+(\ref{sq-y-3})\Big{)}\right|<\infty. So equation (27) is finite.

Because $\sum_{t=1}^{T}\sum_{s=1}^{n}\ln(1-\phi_{0}\tau_{s})$ in (26) only relates to $\phi_{0}$, this term goes away when taken second derivative with respect to other parameters. Similar to the proof of $\mathbb{E}\left|\frac{\partial^{2}\ln f(\varepsilon_{s,t}(\boldsymbol{\theta}_{0}))}{\partial\phi_{0}\partial\phi_{0}}\right|<\infty$, we can show that $\mathbb{E}\left|\frac{\partial^{2}\ln f(\varepsilon_{s,t}(\boldsymbol{\theta}_{0}))}{\partial\phi_{i}\partial\phi_{j}}\right|<\infty$ for $i=0,1,\ldots,p$ and $j=1,\ldots,p$, i.e.,

[TABLE]

Other elements in the matrix (26) equal to those in $\mathbb{E}\left|\frac{\partial^{2}\ln f(\varepsilon_{s,t}(\boldsymbol{\theta}_{0}))}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{\prime}}\right|$ and they are also finite.

[TABLE]

Then we can apply the ergodic theorem [3] and conclude that

[TABLE]

Recall the equation (22), we have proved that $\frac{1}{\sqrt{nT}}\frac{\partial\mathcal{L}_{n,T}(\boldsymbol{\theta}_{0})}{\partial\boldsymbol{\theta}}$ has the limiting distribution $N(\boldsymbol{0},B(\boldsymbol{\theta}_{0}))$. With (41), for $\tilde{\boldsymbol{\theta}}_{n,T}$ between $\hat{\boldsymbol{\theta}}_{n,T}$ and $\boldsymbol{\theta}_{0}$, $-\frac{1}{nT}\frac{\partial^{2}\mathcal{L}_{n,T}(\tilde{\boldsymbol{\theta}}_{n,T})}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{\prime}}\xrightarrow{p}A(\boldsymbol{\theta}_{0})$ so we can conclude that $\sqrt{nT}(\hat{\boldsymbol{\theta}}_{n,T}-\boldsymbol{\theta}_{0})\xrightarrow{d}N(\boldsymbol{0},\boldsymbol{\Omega}_{0})$, where $\boldsymbol{\Omega}_{0}=A^{-1}(\boldsymbol{\theta}_{0})B(\boldsymbol{\theta}_{0})A^{-1}(\boldsymbol{\theta}_{0})$. ∎

6 Numerical Results

6.1 Simulation Study

In this section, we conduct simulation experiments to examine the estimators’ behavior for finite samples. We look at two PSTAR-ANN$(1)$ models with one and two neurons with model parameters specified below:

[TABLE]

Simulations are conducted in a 30 by 30 lattice grid, so $n=900$ and $p=1$, $T=30$. Random errors are sampled respectively from three distributions (standard normal, rescaled t-distribution and Laplace distribution) with variance 1. We generated data for two exogenous variables, observed at different time points $t$ and location $s$. Let $X_{t}=\begin{pmatrix}x_{11,t}&\ldots&x_{1n,t}\\ x_{21,t}&\ldots&x_{2n,t}\end{pmatrix}^{\prime}$. Usually we would like to normalize predictors before fitting a neural network model to avoid the computation overflow [14] so values of $x_{s,t}[i],i=1,2$, were generated independently from normal distributions $N(0,1.5^{2})$ and $N(0,3^{2})$ respectively. The log-likelihood function $\mathcal{L}_{n,T}(\boldsymbol{\theta})$ is given in (44) and we use L-BFGS-B method [6, 24] (recommended for bound constrained optimization) to find the parameter estimates $\hat{\boldsymbol{\theta}}$ which maximize (44).

[TABLE]

For the models under consideration, we estimated the covariance of the asymptotic normal distribution equation (21). Since matrices $A(\boldsymbol{\theta}_{0})$ and $B(\boldsymbol{\theta}_{0})$ involve expected values with respect to the true parameter $\boldsymbol{\theta}_{0}$, given merely observations, in practice they can be estimated as follows:

[TABLE]

where

[TABLE]

Using (23) and (24), we can calculate $\hat{A}(\boldsymbol{\theta}_{0}),\hat{B}(\boldsymbol{\theta}_{0})$ to assess the asymptotic properties of parameter estimates. Note that the derivative of the log-likelihood with respect to $\phi_{0}$ cannot be calculated directly because it requires taking derivative with respect to a log-determinant of $I_{n}-\phi_{0}W_{n}$. For small sample sizes, we can compute the determinant directly and get the corresponding derivatives; but for large sample sizes, for example a dataset with $n=900$ observations, $W_{n}$ is a $900\times 900$ weight matrix which makes it impossible to calculate the derivative directly. Since $W_{n}$ is a square matrix, we can apply the spectral decomposition such that $W_{n}$ can be expressed in terms of its $n$ eigenvalue-eigenvector pairs in (7). So we can apply the following approach to calculate the derivative of $\ln|I_{n}-\phi_{0}W_{n}|$, which greatly reduces the burden of computations (Viton [20]).

[TABLE]

Further the derivatives of the log-likelihood function with respect to $\phi_{0}$ is

[TABLE]

Finally we can estimate the covariance matrix by equation (45).

[TABLE]

In each simulation study, we compute $\hat{\boldsymbol{\theta}}$ for each of 200 replicates. The estimated $\hat{\boldsymbol{\Omega}}$ of the asymptotic covariance matrix $\hat{\boldsymbol{\Omega}}$ is computed based on a sample with $n=10000,T=100$ simulated observations. Table 1 compares the empirical mean and standard errors (in parentheses) of $\hat{\boldsymbol{\theta}}$ with the true value and their estimated asymptotic standard deviations. From simulation results of the two models, the empirical standard deviations of $\hat{\boldsymbol{\theta}}$ are close to the asymptotic standard deviations, which implies that the estimators’ large finite sample behavior roughly matches their asymptotic distributions. Note that when $\varepsilon_{t}$ is sampled from a Laplace distribution, this covariance matrix cannot be computed because its second order derivative is not differentiable at [math]. But the simulated $\hat{\boldsymbol{\theta}}$’s still exhibit normal properties. Normal plots for parameter estimates are shown in Figure 3 and give a strong indication of normality.

6.2 Real Data Example

Spatial models have a lot of applications in understanding spatial interactions in cross-sectional data. In our first chapter we applied a partially specified spatial autoregressive model to understand the relationships between vote choices and social factors. In this chapter, we want to use a partially specified space time autoregressive model to further analyze the time influence in the electoral dynamics.

We focus on the proportion of votes cast for U.S. presidential candidates at the county level in 2004. Counties are grouped by state, and let $Y_{t},Y_{t-1}$ (so $t=1,2$, i.e., observe $Y_{1}$ and $Y_{2}$) be the corresponding fraction of votes (vote-share) in a county for the Democratic candidate in 2004 and 2000. Predictors $X_{t}$ are chosen from economic and social factors covering the living standard, economy development and racial distribution.

Figure 4 shows the observed values of $Y_{2}$ for 2004 and $Y_{1}$ for 2000 in a US map. Despite the strong spatial correlation (by Moran’s Test on $Y_{t}$ test statistic = $52.4$, P-value $<2.2\times 10^{-16}$), these heat maps also exhibit the correlation across time since the two heat maps look rather similar. This indicates that $Y_{t}$, the fraction of vote-share for Democratic candidate, is not independently distributed across the space or time. Therefore we consider fitting a space time model to the data.

In our analysis, we exclude the four U.S. counties with no neighbors (San Juan, Dukes, Nantucket, Richmond) to avoid the non-singularity of our spatial weight matrix $W_{n}$ in the modeling, so the total number of observations is $n=3107$. Continuing our analysis in the first chapter, the selected explanatory variables are percent residents under 18 years in 2004 $X_{1,t}$ (UNDER18), percent white residents in 2004 $X_{2,t}$ (WHITE), percent residents below poverty line in 2004 $X_{3,t}$ (pctpoor).

We also assume the random error follows a scaled $t(8)$ distribution and, similar to previous chapter, perform variable transformations as follows:

[TABLE]

Figure 5 illustrates histograms of $Y_{t}^{\ast}$ (first row) and histograms of exogenous variables $X_{1,t}^{\ast},X_{2,t}^{\ast},X_{3,t}^{\ast}$ when $t=1,2,3$ ($t=3$ represents the year 2008) respectively. Comparing their histograms at different years, we can observe that the distributions look similar so we may consider $X_{t}$ and $Y_{t}$ as stationary processes across time.

The estimated PSAR-ANN model in chapter 1 is:

[TABLE]

In this chapter we would like to add time into the model and we fit two PSTAR-ANN$(1)$ models with one and two neurons respectively. Similarly we find the parameter estimates by maximizing the corresponding log-likelihood functions and use the L-BFGS-B algorithm to search for the optimum. Detailed optimization steps are similar to those in chapter 1. The model fits are shown below. One is the PSTAR-ANN$(1)$ with one neuron:

[TABLE]

Another is the PSTAR-ANN$(1)$ with two neurons:

[TABLE]

Comparing the three models (6.2), (6.2) and (6.2), the coefficients estimates are all positive so it is apparent that there exist a positive space correlation, between $y_{s,t}$ and its neighbors, and also a positive time correlation between $Y_{t}$ and $Y_{t-1}$. The P-values of Moran’s test statistic of PSTAR-ANN$(1)$ model residuals (residuals of model (6.2) and (6.2)) are higher than that of model (6.2), which indicates that PSTAR-ANN$(1)$ models are able to describe more spatial correlations than the PSAR-ANN model.

For the preliminary comparison purpose, we compare the AICs (AIC $=2\#\text{parameters}-2\ln L_{n,T}(\hat{\boldsymbol{\theta}})$) of the three models (See table 2). For likelihood ratio test ($\mathcal{H}_{0}$: Model (6.2) is adequate, $\mathcal{H}_{1}$: model (6.2) is adequate), the test statistic $-2\ln L_{\text{Model (\ref{psar-ann-model})}}+2\ln L_{\text{Model (\ref{pstar-oneneuron-model})}}=287.17$ with $df=6$, P-value $<0.05$, so we rejected $\mathcal{H}_{0}$ and conclude that the PSTAR-ANN$(1)$ model with one neuron is a better fit. Similarly we apply the same method to compare the two PSTAR-ANN$(1)$ models and conclude that the model with two neurons is better (the test statistic $-2\ln L_{\text{Model (\ref{pstar-oneneuron-model})}}+2\ln L_{\text{Model (\ref{pstar-twoneuron-model})}}=40.33$ with $df=4$, P-value $<0.05$).

The covariance matrices for the parameter estimates of model (6.2) and (6.2) are calculated and the $95\%$ confidence intervals for the model parameters are shown in Tables 3 and 4.

From Table 3, all the parameters, except $X_{3,t}$ (pctpoor), are significant at 0.05 significance level. Table 4 shows the $95\%$ level of parameter estimates in model (6.2).

From Table 3 and 4, we can see that values of $y_{s,t}$ are positively spatially correlated in both space and time. Looking at the signs of parameter estimates of coefficients, we can see that the sign of variable UNDER18 in model (6.2) is negative while positive in model (6.2) and (6.2). Considering its parameter estimate significant in all models, this indicates that age and vote-shares for Democratic candidates can be dependent but the percent residents under 18 may not be a good measurement for this social factor. We should consider using other age related variables to predict $Y_{t}$ such as the percent young voters between 18 and 30 years old. Variable WHITE is negatively correlated with $Y_{t}$ in all three fitted models and this negative correlation accords with our common sense that white voters tend to support the Republican candidate. The last variable pctpoor is bit tricky because it is not significant in model (6.2) but is significant in the neural network component in model (6.2). Regarding to this, it needs further assessment to decide if pctpoor should be included in the model. In chapter 3, we will further discuss the model selection in detail. To conclude, our proposed model PSTAR-ANN appears to successfully capture some presidential election dynamics over both space and time. It allows for non-Gaussian random errors and is flexible in learning nonlinear relationships between the response and exogenous variables.