# An $L^{2}$-bound for the Barban–Vehov weights

Olivier Ramaré  and  Sebastian Zuniga-Alterman

###### Abstract.

Let $\lambda$ the Barban–Vehov weights, defined in (1). Let $X\geq z_{1}\geq
100$ and $z_{2}=z_{1}^{\tau}$ for some $\tau>1$. We prove that

$\sum_{n\leq
X}\frac{1}{n}\Bigl{(}\sum_{\begin{subarray}{c}d|n\end{subarray}}\lambda_{d}\Bigr{)}^{2}\leq
f(\tau)\frac{\log X}{\log(z_{2}/z_{1})},$

for a completely determined function $f:(1,\infty)\to\mathbb{R}_{>0}$. In
particular, we may take $f(2)=30$, saving more than a factor of $5$ on what
was the best known result for $\tau=2$. Two related estimates are also
provided for general $\tau>1$.

∗ The second author has been supported by the Finnish Centre of Excellence in
Randomness and Structures grant 346307 of the Academy of Finland,

##  1\. Introduction and results

In [1], M. Barban and P. Vehov introduced the weights

(1)
$\lambda_{d}=\lambda_{d}(z_{2},z_{1})=\mu(d)\frac{\log^{+}(z_{2}/d)-\log^{+}(z_{1}/d)}{\log(z_{2}/z_{1})},$

defined for fixed $z_{2}>z_{1}>0$, where $\log^{+}=\max\\{\log,0\\}$. They
discovered that the $L^{2}$-mean of $\sum_{d|n}\lambda_{d}$ is remarkably
small, even when $n$ is smaller than $z_{2}^{2}$.

The Barban–Vehov weights were further studied for instance by S. Graham [2],
Y. Motohashi [5, $\S$ 1.3] and M. Haye Betah [4].

In particular, Y. Motohashi found that, in order to use such weights to obtain
log-free density estimates, it is enough to upper bound their harmonic
$L^{2}$-average and that, in fact, it is sufficient to bound $\sum_{n\geq
1}\bigl{(}\sum_{\begin{subarray}{c}d|n\end{subarray}}\lambda_{d}\bigr{)}^{2}/n^{1+\varepsilon}$,
for any $\varepsilon>0$. Here is our first result.

######  Theorem 1.1.

Let $\varepsilon>0$ and $z_{1}\geq 100$. Consider $z_{2}=z_{1}^{\tau}$ for
some $\tau>1$. Then,

$\sum_{n\geq
1}\frac{1}{n^{1+\varepsilon}}\biggl{(}\sum_{\begin{subarray}{c}d|n\end{subarray}}\lambda_{d}\biggr{)}^{2}\leq\frac{e^{\gamma\varepsilon}}{\varepsilon\log(z_{2}/z_{1})}\frac{A(\tau+1)+[B-\tfrac{C}{z_{1}^{\varepsilon}}+(B-\tfrac{C}{z_{2}^{\varepsilon}})\tau^{2}]\varepsilon\log
z_{1}}{\tau-1},$

where $A=1.084$, $B=1.301$, $C=0.318$.

The next is a result that enables one to compare Theorem 1.1 to previous work.

######  Corollary 1.2.

When $X\geq z_{1}\geq 100$ and $z_{2}=z_{1}^{\tau}$, for some $\tau>1$, we
have

$\sum_{n\leq
X}\frac{1}{n}\biggl{(}\sum_{\begin{subarray}{c}d|n\end{subarray}}\lambda_{d}\biggr{)}^{2}\leq
3.09\frac{\log
X}{\log(z_{2}/z_{1})}\frac{1.084(\tau+1)+1.301(1+\tau^{2})-0.116}{\tau-1}.$

In particular, when $\tau=2$, the right hand side bound reads $30\log X/\log
z_{1}$.

The above estimate not only generalises the result of [4], but also saves more
than a factor of $5$ on the special case considered therein, namely when
$\tau=2$. Note that a more precise explicit study of these weights when $X\geq
z_{2}^{2}$ has been the subject of [9] by the second author.

On the other hand, while studying zeros of the Riemann zeta function, S.
Graham also used the Barban–Vehov weights. Here is a corollary that may be
compared to [3, Lemma 9].

######  Corollary 1.3.

Let $X\geq z_{1}\geq 100$ and $z_{2}=z_{1}^{\tau}$, for some $\tau>1$. Then,
for any $\alpha\in[1/2,1]$, we have

$\sum_{n\leq
X}\frac{1}{n^{2\alpha-1}}\biggl{(}\sum_{\begin{subarray}{c}d|n\end{subarray}}\lambda_{d}\biggr{)}^{2}\leq
3.09X^{2-2\alpha}\frac{\log
X}{\log(z_{2}/z_{1})}\frac{1.084(\tau+1)+1.301(1+\tau^{2})-0.116}{\tau-1}.$

As the reader will see, the proofs of the above results are remarkably simple,
a fact that is reflected thanks to the very moderate constants that appear
throughout our bounds. It is noteworthy that, giving our historical remarks,
it has taken quite a while to reach the reasonable bounds we present in
Theorem 1.1, Corollary 1.2 and Corollary 1.3.

Notation. We shall use throughout the notation $\sigma=1+\varepsilon$ — we
will use both notations interchangeably. Further, and in accordance with
earlier work on the subject, we use the family of functions

(2) $\check{m}_{q}(X,\sigma)=\sum_{\begin{subarray}{c}n\leq X\\\
(n,q)=1\end{subarray}}\frac{\mu(n)}{n^{\sigma}}\log\biggl{(}\frac{X}{n}\biggr{)}.$

We shall also require the Euler $\varphi_{s}$-function defined, for any
complex number $s$, by
$\varphi_{s}(q)=q^{s}\prod_{p|q}\left(1-{1}/{p^{s}}\right)$. Finally, the
variable $p$ always corresponds to a prime number.

##  2\. Auxiliary results

The following result may be found in [6, Lemma 5.1].

######  Lemma 2.1.

Let $\varepsilon>0$. Then $\zeta(1+\varepsilon)\leq
e^{\gamma\varepsilon}/\varepsilon$.

On the other hand, we have

######  Lemma 2.2.

Suppose that $Y\geq 100$. Then

$\sum_{p\leq Y}\frac{\log p}{p}\log\biggl{(}\frac{Y}{p}\biggr{)}\geq
0.318\log^{2}Y.$

The above sum is in fact asymptotic to $\frac{1}{2}\log^{2}Y$, converging very
slowly.

###### Proof.

We verify the claimed inequality for $Y\in[100,10^{8}]$ by using the Pari-GP
script GP/MertensAux.gp.

Suppose now that $Y\geq 10^{8}$. By summation by parts, $\sum_{p\leq
Y}\frac{\log p}{p}\log(\tfrac{Y}{p})$ equals

$\displaystyle\log Y\sum_{p\leq 100}\frac{\log p}{p}-\sum_{p\leq
100}\frac{\log^{2}(p)}{p}+\int_{100}^{Y}(\vartheta(t)-\vartheta(100))\frac{1+\log(Y/t)}{t^{2}}dt$
(3) $\displaystyle\qquad\geq 3.369\log
Y-8.739+\int_{100}^{Y}(\vartheta(t)-\vartheta(100))\frac{1+\log(Y/t)}{t^{2}}dt,$

where $\vartheta(t)=\sum_{p\leq t}\log p$ is the Chebyshev function. Now, by
[8, Thm. 6], we have that $\vartheta(t)>0.998697\ t$ for any $t\geq 1155901$.
On the other hand, one can verify that $\vartheta(t)>0.835\ t$ provided that
$100\leq t\leq 1155901$. Hence,

$\displaystyle\int_{100}^{Y}(\vartheta(t)-\vartheta(100))\frac{1+\log(Y/t)}{t^{2}}dt$
$\displaystyle\geq\int_{100}^{Y}(0.835\
t-\vartheta(100))\frac{1+\log(Y/t)}{t^{2}}dt$ (4)
$\displaystyle\geq\log^{2}(Y)\left[\frac{1}{2}+\frac{0.845}{\log
Y}-\frac{5.058}{\log^{2}Y}\right].$

By combining (2) and (2), and then using that $Y\geq 10^{8}$, we arrive at

$\sum_{p\leq Y}\frac{\log
p}{p}\log\biggl{(}\frac{Y}{p}\biggr{)}\geq\log^{2}(Y)\left[\frac{1}{2}-\frac{13.8}{\log^{2}(10^{8})}\right]\geq
0.376\log^{2}Y.$

As $0.376>0.318$, we derive the inequality of the statement.

Finally, thanks to the prime number theorem, $\vartheta(t)\sim t$, so we
indeed have that $\sum_{p\leq Y}\frac{\log
p}{p}\log(\frac{Y}{p})\sim\frac{1}{2}\log^{2}Y$. ∎

##  3\. Some lemmas involving $\check{m}_{q}(X,\sigma)$

The following result is crucial for our improvement on [4]. It can be found in
[7, Theorem 1.1].

######  Lemma 3.1.

Let $X>0$. When $k\geq 1$ and $\sigma\geq 1$, we have

$0\leq\sum_{\begin{subarray}{c}n\leq X\\\
(n,q)=1\end{subarray}}\frac{\mu(n)}{n^{\sigma}}\log^{k}\biggl{(}\frac{X}{n}\biggr{)}\leq
1.00303\frac{q}{\varphi(q)}\bigl{(}k+(\sigma-1)\log X\bigr{)}\log^{k-1}X.$

When $k\geq 2$, we may replace $1.00303$ by 1.

We also have the following result.

######  Lemma 3.2.

Let $Y\geq 100$. Then

$U(Y)=\sum_{\delta\leq
Y}\frac{\mu^{2}(\delta)\varphi_{\sigma}(\delta)}{\delta^{2\sigma}}\check{m}_{\delta}\Bigl{(}\frac{Y}{\delta};\sigma\Bigr{)}^{2}\leq\left(A+(\sigma-1)\left(B-\tfrac{C}{Y^{\varepsilon}}\log
Y\right)\right)\log Y.$

where $A=1.084$, $B=1.301$, $C=0.318$.

###### Proof.

Each $\check{m}_{\delta}(y/\delta;\sigma)$ is non-negative, thanks to Lemma
3.1 with $k=1$. Thus, we have the following estimation

(5) $U(Y)\leq 1.00303\sum_{\delta\leq
Y}\frac{\mu^{2}(\delta)\varphi_{\sigma}(\delta)}{\delta^{2\sigma}}\frac{\delta}{\varphi(\delta)}\check{m}_{\delta}\Bigl{(}\frac{Y}{\delta};\sigma\Bigr{)}\biggl{(}1+\varepsilon\log\bigg{(}\frac{Y}{\delta}\bigg{)}\biggr{)}.$

Let us define the auxiliary non-negative multiplicative function $f$ on primes
$p$ by $f(p)=\frac{1-p^{-\sigma}}{1-1/p}>0$ and on prime powers by
$f(p^{k})=0$ for $k\in\mathbb{Z}_{>1}$. We observe that

(6)
$\frac{\mu^{2}(\delta)\varphi_{\sigma}(\delta)}{\delta^{\sigma}}\frac{\delta}{\varphi(\delta)}=\mu^{2}(\delta)\sum_{a|\delta}f(a).$

On using identity (6) on estimation (5), and recalling definition (2) , we
derive

$\displaystyle U(Y)$ $\displaystyle\leq 1.00303\sum_{a\leq
Y}\mu^{2}(a)\frac{f(a)}{a^{\sigma}}\sum_{\begin{subarray}{c}b\leq Y/a\\\
(b,a)=1\end{subarray}}\frac{\mu^{2}(b)}{b^{\sigma}}\check{m}_{ab}\Bigl{(}\frac{Y}{ab};\sigma\Bigr{)}\biggl{(}1+\varepsilon\log\bigg{(}\frac{Y}{ab}\bigg{)}\biggr{)}$
$\displaystyle\leq 1.00303\sum_{a\leq
Y}\mu^{2}(a)\frac{f(a)}{a^{\sigma}}\biggl{(}1+\varepsilon\log\bigg{(}\frac{Y}{a}\bigg{)}\biggr{)}\sum_{\begin{subarray}{c}d\leq
Y/a\\\
(d,a)=1\end{subarray}}\frac{\mu^{2}(d)}{d^{\sigma}}\log\biggl{(}\frac{Y/a}{d}\biggr{)}\sum_{bc=d}\mu(c)$
$\displaystyle\qquad-1.00303\varepsilon\sum_{a\leq
Y}\mu^{2}(a)\frac{f(a)}{a^{\sigma}}\sum_{\begin{subarray}{c}d\leq Y/a\\\
(d,a)=1\end{subarray}}\frac{\mu^{2}(d)}{d^{\sigma}}\log\biggl{(}\frac{Y/a}{d}\biggr{)}\sum_{bc=d}\mu(c)\log
b.$

By Möbius inversion, the innermost sum of the above first term reduces to
$d=1$. Moreover, as $\sum_{bc=d}\mu(c)\log(b)=\Lambda(d)$, where $\Lambda$ is
the von Mangoldt function, the innermost sum of the above second term is
supported on prime numbers. Hence,

$\displaystyle\frac{U(Y)}{1.00303}\leq$ $\displaystyle\sum_{a\leq
Y}\mu^{2}(a)\frac{f(a)}{a^{\sigma}}\biggl{(}1+\varepsilon\log\bigg{(}\frac{Y}{a}\bigg{)}\biggr{)}\log\bigg{(}\frac{Y}{a}\bigg{)}$
$\displaystyle\qquad-\varepsilon\,\sum_{a\leq
Y}\mu^{2}(a)\frac{f(a)}{a^{\sigma}}\sum_{\begin{subarray}{c}p\leq Y/a\\\
(p,a)=1\end{subarray}}\frac{\log
p}{p^{\sigma}}\log\bigg{(}\frac{Y/a}{p}\bigg{)}.$

In the above inequality, we consider the two terms coming from $a=1$. Then, by
using the inequality $p^{-\sigma}\geq y^{-\varepsilon}p^{-1}$, valid whenever
$p\leq y$, we derive

$\displaystyle\frac{U(Y)}{1.00303}\leq$ $\displaystyle(1+\varepsilon\log
Y)\log Y+\sum_{2\leq a\leq
Y}\mu^{2}(a)\frac{f(a)}{a^{\sigma}}\biggl{(}1+\varepsilon\log\bigg{(}\frac{Y}{a}\bigg{)}\biggr{)}\log\bigg{(}\frac{Y}{a}\bigg{)}$
$\displaystyle\qquad-\varepsilon\,\sum_{\begin{subarray}{c}p\leq
Y\end{subarray}}\frac{\log
p}{p^{\sigma}}\log\bigg{(}\frac{Y}{p}\bigg{)}-\varepsilon\,\sum_{2\leq a\leq
Y}\mu^{2}(a)\frac{f(a)}{a^{\sigma}}\sum_{\begin{subarray}{c}p\leq Y/a\\\
(p,a)=1\end{subarray}}\frac{\log
p}{p^{\sigma}}\log\bigg{(}\frac{Y/a}{p}\bigg{)}.$

Next, by trivially bounding the above last term and using Lemma 2.2, we have

$\displaystyle\frac{U(Y)}{1.00303}$ $\displaystyle\leq(1+\varepsilon\log
Y)\log Y-\frac{0.318\varepsilon}{Y^{\varepsilon}}\log^{2}Y$
$\displaystyle\quad+\biggl{[}-1+\prod_{p\geq
2}\biggl{(}1+\frac{p^{-\varepsilon}(1-p^{-\varepsilon})}{p(p-1)}\biggr{)}\biggr{]}\biggl{(}1+\varepsilon\log\left(\frac{Y}{2}\right)\biggr{)}\log\left(\frac{Y}{2}\right).$

Furthermore, by using that $t\in\mathbb{R}_{>0}\mapsto t^{-1}(1-t^{-1})$ has a
maximum $1/4$, we deduce that

$\displaystyle-1+\prod_{p\geq
2}\biggl{(}1+\frac{p^{-\varepsilon}(1-p^{-\varepsilon})}{p(p-1)}\biggr{)}\leq-1+\biggl{(}1+\frac{2^{-\varepsilon}(1-2^{-\varepsilon})}{2}\biggr{)}\prod_{p\geq
3}\biggl{(}1+\frac{1}{4p(p-1)}\biggr{)}$
$\displaystyle\qquad\leq-1+\biggl{(}1+\frac{2^{-\varepsilon}(1-2^{-\varepsilon})}{2}\biggr{)}1.08\leq
0.08+0.54\cdot 2^{-\varepsilon}(1-2^{-\varepsilon}),$

so that

$\displaystyle\bigl{(}0.08$ $\displaystyle+0.54\cdot
2^{-\varepsilon}(1-2^{-\varepsilon})\bigr{)}(1+\varepsilon\log Y)$
$\displaystyle\leq 0.08+\varepsilon\biggl{(}0.08+0.54\cdot
2^{-\varepsilon}(1-2^{-\varepsilon})+0.54\frac{2^{-\varepsilon}(1-2^{-\varepsilon})}{\varepsilon\log
Y}\biggr{)}\log Y$ $\displaystyle\leq 0.08+0.297\,\varepsilon\log Y,$

where we have used that $t\in\mathbb{R}_{\geq 0}\mapsto(1-2^{-t})/t$ is
decreasing and has as maximum $\log 2$ at $t=0$. We conclude the result by
observing that

$\displaystyle\frac{U(Y)}{1.00303}$ $\displaystyle\leq(1+\varepsilon\log
Y)\log
Y-\frac{0.318\varepsilon}{Y^{\varepsilon}}\log^{2}Y+\biggl{(}0.08+0.297\varepsilon\log
Y\biggr{)}\log\left(\frac{Y}{2}\right)$ $\displaystyle\leq 1.08\log
Y+\varepsilon\log^{2}Y\bigg{(}1.297-\tfrac{0.318}{Y^{\varepsilon}}\bigg{)}.$

∎

##  4\. Main result

Proof of Theorem 1.1. Recall definition (1). Let us define

(7)
$\Sigma(\sigma)=\Sigma(\sigma,z_{2},z_{1})=\sum_{n}\frac{1}{n^{\sigma}}\Bigl{(}\sum_{\begin{subarray}{c}d|n\end{subarray}}\lambda_{d}\Bigr{)}^{2}.$

By expanding the square, we observe that

(8)
$\Sigma(\sigma)=\sum_{d,e}\lambda_{d}\lambda_{e}\sum_{[d,e]|n}\frac{1}{n^{\sigma}}=\zeta(\sigma)\sum_{d,e}\frac{\lambda_{d}\lambda_{e}}{[d,e]^{\sigma}},$

where the variables $d$ and $e$ may be restricted to $d,e\leq z_{2}$.

Let us use now the Selberg diagonalisation process. For any pair of square-
free numbers $d,e$, we have

(9)
$\frac{d^{\sigma}e^{\sigma}}{[d,e]^{\sigma}}=(d,e)^{\sigma}=\sum_{\delta|(d,e)}\varphi_{\sigma}(\delta).$

As the support of $\lambda$ is contained in the set of square-free numbers, we
may derive from (9) the following expression

(10)
$\sum_{d,e}\frac{\lambda_{d}\lambda_{e}}{[d,e]^{\sigma}}=\sum_{\delta}\frac{\varphi_{\sigma}(\delta)}{\delta^{2\sigma}}\biggl{(}\sum_{\begin{subarray}{c}\ell\\\
(\ell,\delta)=1\end{subarray}}\frac{\lambda_{\delta\ell}}{\ell^{\sigma}}\biggr{)}^{2}.$

Thus, by definition (1) and equation (8), we obtain

(11)
$\Sigma(\sigma)\leq\frac{\zeta(\sigma)}{\log^{2}(z_{2}/z_{1})}\sum_{\delta\leq
z_{2}}\frac{\mu^{2}(\delta)\varphi_{\sigma}(\delta)}{\delta^{2\sigma}}\biggl{[}\check{m}_{\delta}\Bigl{(}\frac{z_{2}}{\delta};\sigma\Bigr{)}-\check{m}_{\delta}\Bigl{(}\frac{z_{1}}{\delta};\sigma\Bigr{)}\biggr{]}^{2},$

where we have used the definition (2). Furthermore, as
$\check{m}_{q}(X;\sigma)\geq 0$, thanks to Lemma 3.1, we find that

$\bigg{[}\check{m}_{\delta}\bigg{(}\frac{z_{2}}{\delta};\sigma\bigg{)}-\check{m}_{\delta}\bigg{(}\frac{z_{1}}{\delta};\sigma\bigg{)}\bigg{]}^{2}\leq\check{m}^{2}_{\delta}\bigg{(}\frac{z_{2}}{\delta};\sigma\bigg{)}+\check{m}^{2}_{\delta}\bigg{(}\frac{z_{1}}{\delta};\sigma\bigg{)}$

with equality if and only if
$\check{m}_{\delta}({z_{2}}/{\delta};\sigma)\check{m}_{\delta}({z_{1}}/{\delta};\sigma)=0$.
Therefore, we see from (11) that

$\frac{\log^{2}(z_{2}/z_{1})}{\zeta(\sigma)}\Sigma(\sigma)\leq\sum_{\delta\leq
z_{2}}\frac{\mu^{2}(\delta)\varphi_{\sigma}(\delta)}{\delta^{2\sigma}}\check{m}^{2}_{\delta}\Bigl{(}\frac{z_{2}}{\delta};\sigma\Bigr{)}+\sum_{\delta\leq
z_{1}}\frac{\mu^{2}(\delta)\varphi_{\sigma}(\delta)}{\delta^{2\sigma}}\check{m}^{2}_{\delta}\Bigl{(}\frac{z_{1}}{\delta};\sigma\Bigr{)},$

so that

$\Sigma(\sigma)\leq\frac{\zeta(\sigma)}{\log^{2}(z_{2}/z_{1})}(U(z_{2})+U(z_{1}))=\frac{\zeta(\sigma)}{\log(z_{2}/z_{1})}\frac{U(z_{2})+U(z_{1})}{(\tau-1)\log
z_{1}}.$

Theorem 1.1 follows by applying Lemma 3.2 twice, recalling that
$z_{2}=z_{1}^{\tau}$ and using then Lemma 2.1.

$\square$

With the help of (10), it is worth noticing that, for any $\varepsilon>0$ and
any weight $\lambda$, the sum
$\sum_{d,e}\frac{\lambda_{d}\lambda_{e}}{[d,e]^{\sigma}}$ is always non-
negative.

##  5\. Corollaries

Proof of Corollary 1.2. By using Rankin’s trick, we have

$\Sigma(0)=\sum_{n\leq
X}\frac{1}{n}\biggl{(}\sum_{\begin{subarray}{c}d|n\end{subarray}}\lambda_{d}\biggr{)}^{2}\leq\sum_{n}\frac{X^{\varepsilon}}{n^{1+\varepsilon}}\Bigl{(}\sum_{\begin{subarray}{c}d|n\end{subarray}}\lambda_{d}\Bigr{)}^{2}$

for any parameter $\varepsilon$ that we shall soon choose. Now, by using
Theorem 1.1, we obtain the bound

(12)
$\Sigma(0)\leq\frac{X^{\varepsilon}e^{\gamma\varepsilon}}{\varepsilon\log(z_{2}/z_{1})}\frac{A(\tau+1)+[B-\tfrac{C}{z_{1}^{\varepsilon}}+(B-\tfrac{C}{z_{2}^{\varepsilon}})\tau^{2}]\varepsilon\log
z_{1}}{\tau-1},$

where $A=1.084$, $B=1.301$, $C=0.318$.

Set $\varepsilon=1/\log X$. Since there is no further assumptions on $z_{2}$,
we use the bound $-Cz_{2}^{-\varepsilon}\leq 0$. Thus, by writing
$u=\varepsilon\log z_{1}$, we obtain

$\displaystyle\Sigma(0)$ $\displaystyle\leq\frac{e^{1+\gamma/\log X}\log
X}{\log(z_{2}/z_{1})}\frac{A(\tau+1)+[B(1+\tau^{2})-\tfrac{C}{e^{u}}]u}{\tau-1}$
$\displaystyle\leq\frac{e^{1+\gamma/\log 100}\log
X}{\log(z_{2}/z_{1})}\frac{A(\tau+1)+B(1+\tau^{2})-\tfrac{C}{e}}{\tau-1}$

where we have used that $X\geq 100$ and $0<u\leq 1$.

In the specific and common situation when $\tau=2$, unlike the general case,
we may take advantage of the factor $Cz_{2}^{-\varepsilon}=Ce^{-2u}$ appearing
in (12). Hence, by using again that $X\geq 100$ and $0<u\leq 1$, we finally
obtain

$\Sigma(0,z_{1}^{2},z_{1})\leq e^{1+\gamma/\log
100}[A(\tau+1)+B(1+\tau^{2})-\tfrac{C}{e}-\tfrac{C\tau^{2}}{e}]\frac{\log
X}{\log(z_{2}/z_{1})}\leq\frac{29.18\log X}{\log(z_{2}/z_{1})}.$

$\square$

Proof of Corollary 1.3. We simply write

$\displaystyle\sum_{n\leq
X}\frac{1}{n^{2\alpha-1}}\biggl{(}\sum_{\begin{subarray}{c}d|n\end{subarray}}\lambda_{d}\biggr{)}^{2}$
$\displaystyle=\sum_{n\leq
X}\frac{(\sum_{\begin{subarray}{c}d|n\end{subarray}}\lambda_{d})^{2}}{n}n^{2-2\alpha}\leq
X^{2-2\alpha}\sum_{n\leq
X}\frac{1}{n}\bigg{(}\sum_{\begin{subarray}{c}d|n\end{subarray}}\lambda_{d}\bigg{)}^{2},$

since $2-2\alpha\geq 0$. Then the result follows by using Corollary 1.2.
$\square$

## References

  * [1] M.B. Barban and P.P. Vehov.  An extremal problem.  Trudy Moskov. Mat., 18:83–90, 1968.
  * [2] S.W. Graham.  An asymptotic estimate related to Selberg’s sieve.  J. Number Theory, 10:83–94, 1978.
  * [3] S.W. Graham.  The distribution of squarefree numbers.  J. London Math. Soc. (2), 24(1):54–64, 1981.
  * [4] M. Haye Betah.  Explicit expression of a Barban & Vehov theorem.  Funct. Approx. Comment. Math., 60(2):177–193, 2019.
  * [5] Y. Motohashi.  Sieve Methods and Prime Number Theory.  Tata Lectures Notes, page 205, 1983.
  * [6] O. Ramaré.  An explicit density estimate for Dirichlet ${L}$-series.  Math. Comp., 85(297):335–356, 2016.
  * [7] O. Ramaré and S. Zuniga-Alterman.  Möbius function: an identity factory with applications.  arXiv:2312.05138, 2023.
  * [8] L. Schoenfeld.  Sharper bounds for the Chebyshev Functions $\vartheta(x)$ and $\psi(x)$ ii.  Math. Comp., 30(134):337–360, 1976.
  * [9] S Zuniga-Alterman.  Smoothing and cancellation: the Barban-Vehov sieve made explicit.  PhD thesis, Mathematics, Université Paris-Diderot, 2020.  81 pp.

O. RAMARÉ - CNRS / Institut de Mathématiques de Marseille, Aix Marseille
Université, U.M.R. 7373, Campus de Luminy, Case 907, 13288 MARSEILLE Cedex 9,
France.
<EMAIL_ADDRESS>
S. ZUNIGA ALTERMAN - Department of Mathematics and Statistics, University of
Turku, 20014 TURKU, Finland.
<EMAIL_ADDRESS>


# A Dynamic Spatiotemporal Stochastic Volatility Model with an Application to
Environmental Risks

Philipp Otto Leibniz University, Hannover, Germany, email:
<EMAIL_ADDRESS>Osman Doğan Economics Department, Istanbul
Technical University, Istanbul, Turkey, email<EMAIL_ADDRESS>Süleyman
Taşpınar Department of Economics, Queens College, The City University of New
York, United States, email<EMAIL_ADDRESS>
###### Abstract

This article introduces a dynamic spatiotemporal stochastic volatility (SV)
model with explicit terms for the spatial, temporal, and spatiotemporal
spillover effects. Moreover, the model includes time-invariant site-specific
constant log-volatility terms. Thus, this formulation allows to distinguish
between spatial and temporal interactions, while each location may have a
different volatility level. We study the statistical properties of an outcome
variable under this process and show that it introduces spatial dependence in
the outcome variable. Further, we present a Bayesian estimation procedure
based on the Markov Chain Monte Carlo (MCMC) approach using a suitable data
transformation. After providing simulation evidence on the proposed Bayesian
estimator’s performance, we apply the model in a highly relevant field, namely
environmental risk modeling. Even though there are only a few empirical
studies on environmental risks, previous literature undoubtedly demonstrated
the importance of climate variation studies. For example, for local air
quality in Northern Italy in 2021, we show pronounced spatial and temporal
spillovers and larger uncertainties/risks during the winter season compared to
the summer season.

JEL-Classification: C13, C21, C31.
Keywords: Environmental risk, MCMC, spatial dependence, spatiotemporal
stochastic volatility, air quality

##  1 Introduction

When analyzing geo-referenced data, statistical models have to account for
instantaneous spatial correlations due to the geographical proximity between
the observations. This is commonly known as Tobler’s first law of geography
“everything is related to everything else, but near things are more related
than distant things” (Tobler,, 1970). This observation was already noted by
Ronald A. Fisher in 1935 as follows, “the widely verified fact that patches in
close proximity are commonly more alike, as judged by the yield of crops, than
those which are further apart” (Fisher,, 1935). Even though the similarity is
typically considered to be in the (conditional) mean level at each location,
there might also be spatial correlations in the (conditional) variance or
variation of the random process. In particular, for small-scale spatial units,
the variance of the process is increased (known as Arbia’s law of geography,
Arbia and Espa,, 1996). In addition to the instantaneous spatial correlations,
we also have to account for the natural temporal correlations, which usually
occur if we repeatedly observe a random process over time. The closer two
observations are to each other in time, the more strongly they can correlate
in general. In this paper, we introduce a new model for spatial, temporal, and
spatiotemporal correlations in the log-volatilities allowing for additional
random errors in the mean and volatility equation. Furthermore, we apply this
model to environmental data and show for the first time how it can be used to
analyze environmental risk factors such as air pollution.

There are generally two ways to account for spatial/cross-sectional
correlations in spatial statistics. Firstly, it can be modeled in the
covariance matrix of the process, where each entry is supposed to follow a
certain (non-)parametric covariance function depending on the distance between
their locations. This idea is typically known as the geostatistical approach
(Cressie and Wikle,, 2015; Zimmerman,, 2019). The selection of a suitable
parametric covariance function with certain properties such as stationarity,
separability and full symmetry is one of the main modeling issues of this
approach (Porcu et al.,, 2016; Huang et al.,, 2011). See Gneiting et al.,
(2007) for a review on the spatiotemporal covariance functions suggested in
the literature. Secondly, the observations on an outcome variable can be
explicitly correlated with the adjacent observations, where the adjacency is
defined fairly generally by a spatial weights matrix. This second approach is
referred to as spatial autoregression in spatial econometrics, where the
spatial lags of variables are used to model spatial correlations. See LeSage
and Pace, (2009); Anselin, (1988); Elhorst, (2014); Lee, (2004); Kelejian and
Prucha, (2010) on the specification and estimation issues in spatial
econometrics. Both approaches can be equivalent under certain conditions (see
e.g. Ver Hoef et al., 2018a, , for simultaneous and conditionally
autoregressive models).

In this paper, we consider the second approach. Our model consists of an
outcome and log-volatility equation with separate independent error terms,
whereby the log-volatility process introduces spatial dependence in the
outcome variable. Specifically, the log-volatility equation allows for
spatial, temporal, and spatiotemporal correlations, as well as time-invariant
site-specific effects (unobserved heterogeneity). Also, assuming that the
error terms in both equations have normal distributions, it is possible to
show that the outcome variable has a leptokurtic symmetric distribution under
our suggested model. To introduce a Bayesian estimation approach, we use a
transformation approach such that the outcome equation becomes linear in the
log-volatility terms. We use a Gaussian mixture distribution to approximate
the distribution of the transformed error terms in the outcome equation. This
approximation turns our model into a linear state-space model, where the log-
volatility equation becomes the state equation. Following recent developments
in the precision-based algorithms (Chan and Jeliazkov,, 2009; Chan,, 2017), we
suggest a Gibbs sampler that consists of five steps for the estimation. We
provide simulation evidence showing that the suggested sampler can perform
satisfactorily.

Theoretically, our paper is related to the spatial econometric literature that
addresses the presence of cross-sectional correlations in higher moments of
spatial data. This strand of the literature considers the spatial extensions
of generalized autoregressive conditional heteroskedasticity (GARCH) and
stochastic volatility (SV) models to account for volatility clustering
patterns observed over space (Otto et al.,, 2018; Hølleland and Karlsen,,
2020; Sato and Matsuda,, 2021; Taşpınar et al.,, 2021; Robinson,, 2009; Yan,,
2007). Our model can be considered as the longitudinal data (panel data)
extension of the cross-sectional spatial SV models suggested by Yan, (2007);
Taşpınar et al., (2021); Robinson, (2009). While these studies allow for the
presence of spatial dependence in the log-volatility equations, they do not
include temporal, spatiotemporal and unobserved heterogeneity terms in the
log-volatility equations. Our suggested process is also related to the
separable and non-separable space-time filters considered in the spatial panel
data models for modeling spatiotemporal interactions. For example, Parent and
LeSage, (2012, 2011) consider a separable space-time filter for the outcome
variable, while Lee and Yu, (2015) and Wang and Lee, (2018) consider non-
separable space-time filters for the disturbance terms of spatial panel data
models. In contrast to these studies, we consider a general space-time filter
that also allows for unobserved heterogeneity, i.e., the site-specific
effects, for the log-volatility of an outcome variable.

In an empirical application, we use our suggested model to assess
environmental risk stemming from the variations in the log-volatility of air
quality predictions. Spatial and spatiotemporal interactions in the local
climate and environmental risks were addressed in comparably few studies in
previous empirical research, even though it has been shown that an increased
variation in environmental processes can be harmful (Tewksbury et al.,, 2008;
Paaijmans et al.,, 2013; Vasseur et al.,, 2014; Iaco et al.,, 2012). Most
previous studies focused on correlations in the (conditional) mean levels of
ecological processes (Ver Hoef et al., 2018b, ; Wilby et al.,, 2009, e.g.). In
our case, we model the log-volatility of fine dust concentrations of particles
having a diameter less than $10\mu m$ in Lombardy, Northern Italy. Following
the literature on the ecological processes (Ver Hoef et al., 2018b, ; Wilby et
al.,, 2009, e.g.), we first model the variation in the conditional mean of our
outcome variable through a conventional spatial panel data model that allows
for the unobserved site and time heterogeneity. Our results from this initial
model indicate that there is strong and moderate spatial correlations in the
outcome variable in the model with only site fixed effects and the model with
both site and time fixed effects, respectively. In the next step, we use the
errors from this initial model as an outcome variable in our suggested
specification and aim to model the variations in its log-volatility terms. The
estimation results from our suggested Bayesian approach indicate that the
spatial and temporal effects in the log-volatilities are moderate while the
spatiotemporal effects appear to be of minor importance. We were also able to
detect a noticeable variation in air pollution risk across the year and
identify measurement stations that are associated with higher risks. These are
mostly located in valleys in the Alpine regions.

The rest of this paper proceeds in the following way. In Section 2, we
introduce our suggested model specification, including conditions ensuring the
stability of the model and prior specifications. In Section 3, we investigate
the statistical properties of the suggested model. In Section 4, we provide
the details on the posterior analysis of our model, and state an algorithm for
the estimation. Section 5 provides a simulation study on the performance of
the suggested Gibbs sampler. Then, in Section 6, such stochastic volatility
model is applied to environmental risks for the first time. More precisely,
our focus in on local air quality modeling. Finally, in Section 7, we provide
our concluding remarks.

##  2 Model Specification

Suppose that we observe the spatiotemporal process across a constant set of
$n$ locations in a geographical domain at $T$ equidistant time points. These
locations can be measurement stations (i.e., marked point data), atmospheric
and remotely-sensed data, or sets of municipalities, counties, states (i.e.,
areal data). Moreover, these locations does not have to seen in a strict
geographical sense, but can also be vertices in a network. Let
$\mathbf{Y}_{t}=(y_{1t},y_{2t},\ldots,y_{nt})^{{}^{\prime}}$ be the $n\times
1$ vector of the outcome variable at time $t$ for $t=1,2,\ldots,T$. We assume
the following data generating process (DGP) for $\mathbf{Y}_{t}$:

$\displaystyle\mathbf{Y}_{t}=\mathbf{H}^{1/2}_{t}\mathbf{V}_{t},$ (2.1)

for $t=1,2,\ldots T$, where
$\mathbf{H}^{1/2}_{t}=\operatorname{Diag}(e^{\frac{1}{2}h_{1t}},\ldots,e^{\frac{1}{2}h_{nt}})$
is the $n\times n$ diagonal matrix containing stochastic volatility terms
$h_{it}$’s, which are specified subsequently, and
$\mathbf{V}_{t}=(v_{1t},\ldots,v_{nt})^{{}^{\prime}}$ is the $n\times 1$
vector of disturbance terms. We assume that $v_{it}$’s are i.i.d standard
normal random variables. Let
$\mathbf{h}_{t}=(h_{1t},\ldots,h_{nt})^{{}^{\prime}}$ be the $n\times 1$
vector of stochastic volatility at time $t$. We assume the following process
for $\mathbf{h}_{t}$:

$\displaystyle\mathbf{h}_{t}-\boldsymbol{\mu}=\rho_{1}\mathbf{W}(\mathbf{h}_{t}-\boldsymbol{\mu})+\rho_{2}(\mathbf{h}_{t-1}-\boldsymbol{\mu})+\rho_{3}\mathbf{W}(\mathbf{h}_{t-1}-\boldsymbol{\mu})+\mathbf{U}_{t},$
(2.2)

for $t=1,2,\ldots T$, where
$\boldsymbol{\mu}=(\mu_{1},\ldots,\mu_{n})^{{}^{\prime}}$ is the $n\times 1$
vector of constant means, i.e., the time-invariant site-specific effects, and
$\mathbf{U}_{t}=(u_{1t},\ldots,u_{nt})^{{}^{\prime}}$ is the $n\times 1$
vector of i.i.d. disturbance terms such that $u_{it}\sim N(0,\sigma^{2})$ for
all $i$ and $t$. In (2.2), $\mathbf{W}$ is the $n\times n$ spatial weights
matrix that has zero diagonal elements. This matrix specifies how volatility
terms are related over space. In a network setting, this matrix is
equivalently specified as an adjacency matrix. The scalar parameter $\rho_{1}$
captures contemporaneous spatial correlation, $\rho_{2}$ measures the temporal
effect, i.e., the time dynamic effect, and $\rho_{3}$ represents the
spatiotemporal effect, i.e., the spatial diffusion effect. Let $\mathbf{L}$ be
the matrix of time-lag operator such that
$\mathbf{L}\mathbf{h}_{t}=\mathbf{h}_{t-1}$. Then, (2.2) can be written as

$\displaystyle\left((\mathbf{I}_{n}-\rho_{1}\mathbf{W})-(\rho_{2}\mathbf{L}+\rho_{3}\mathbf{W}\mathbf{L})\right)(\mathbf{h}_{t}-\boldsymbol{\mu})=\mathbf{U}_{t},$
(2.3)

where
$\left((\mathbf{I}_{n}-\rho_{1}\mathbf{W})-(\rho_{2}\mathbf{L}+\rho_{3}\mathbf{W}\mathbf{L})\right)$
is called the general space-time filter (Parent and LeSage,, 2011, 2012; Lee
and Yu,, 2015). Under the assumption that $\rho_{3}=-\rho_{1}\rho_{2}$, this
general filter is separable and decomposes into a product of the space filter
$(\mathbf{I}_{n}-\rho_{1}\mathbf{W})$ and the time filter
$(\mathbf{I}_{n}-\rho_{2}\mathbf{L})$. In our analysis, we do not impose this
restrictive assumption. Let
$\mathbf{S}(\rho_{1})=(\mathbf{I}_{n}-\rho_{1}\mathbf{W})$. Then, under the
assumption that $\mathbf{S}(\rho_{1})$ is invertible, the reduced form of
volatility equation is

$\displaystyle\mathbf{h}_{t}-\boldsymbol{\mu}=\mathbf{S}^{-1}(\rho_{1})\mathbf{A}(\rho_{2},\rho_{3})(\mathbf{h}_{t-1}-\boldsymbol{\mu})+\mathbf{S}^{-1}(\rho_{1})\mathbf{U}_{t},$
(2.4)

where
$\mathbf{A}(\rho_{2},\rho_{3})=(\rho_{2}\mathbf{I}_{n}+\rho_{3}\mathbf{W})$.
When the cross-sectional dimension is fixed, the process for the log-
volatility is stable if all eigenvalues of
$\mathbf{S}^{-1}(\rho_{1})\mathbf{A}(\rho_{2},\rho_{3})$ lie inside the unit
ball (Hamilton,, 1994, Proposition 10.1). Let $\vartheta_{i}(\mathbf{W})$ be
the $i$th eigenvalue of $\mathbf{W}$ for $i=1,\ldots,n$. We assume that the
parameter space of $\rho_{1}$, $\rho_{2}$ and $\rho_{3}$ are chosen such that
the following conditions hold:111In the spatial econometric literature, there
are alternative ways to specify the parameter spaces for spatial
autoregressive parameters. Among others, see Anselin, (1988), LeSage and Pace,
(2009), Lee, (2004), Kelejian and Prucha, (2010) and Elhorst, (2014).

$\displaystyle(i)\max_{1\leq i\leq
n}\left|\vartheta_{i}(\rho_{1}\mathbf{W})\right|<1,\quad\text{and}\;(ii)\,\max_{1\leq
i\leq
n}\left|\vartheta_{i}\left(\mathbf{S}^{-1}(\rho_{1})\mathbf{A}(\rho_{2},\rho_{3})\right)\right|<1.$
(2.5)

The first condition is the sufficient condition for the invertibility of
$\mathbf{S}(\rho_{1})$ (Kelejian and Prucha,, 2010, Lemma 1). By the spectral
radius theorem (Horn and Johnson,, 2012, Theorem 5.6.9), we can use any matrix
norm to define relatively restrictive conditions that ensure the conditions in
(2.5). Let $\|\cdot\|$ be any matrix norm. Then, the sufficient conditions for
(2.5) are (i) $\|\rho_{1}\mathbf{W}\|<1$ and (ii)
$\|\mathbf{S}^{-1}(\rho_{1})\mathbf{A}(\rho_{2},\rho_{3})\|<1$. Note that

$\displaystyle\left\|\mathbf{S}^{-1}(\rho_{1})\mathbf{A}(\rho_{2},\rho_{3})\right\|$
$\displaystyle\leq\left\|\mathbf{S}^{-1}(\rho_{1})\right\|\times\left\|\rho_{2}\mathbf{I}_{n}+\rho_{3}\mathbf{W}\right\|$
$\displaystyle=\|\mathbf{I}_{n}+\rho_{1}\mathbf{W}+\rho^{2}_{1}\mathbf{W}^{2}+\rho^{3}_{1}\mathbf{W}^{3}+\ldots\|\times\|\rho_{2}\mathbf{I}_{n}+\rho_{3}\mathbf{W}\|$
$\displaystyle\leq\left(\|\mathbf{I}_{n}\|+\|\rho_{1}\mathbf{W}\|+\|\rho_{1}\mathbf{W}\|^{2}+\ldots\right)\times\left(|\rho_{2}|+|\rho_{3}|\cdot\|\mathbf{W}\|\right)$
$\displaystyle=\frac{1}{1-\|\rho_{1}\mathbf{W}\|}\times\left(|\rho_{2}|+|\rho_{3}|\cdot\|\mathbf{W}\|\right),$
(2.6)

where last equality follows since we assume that $\|\rho_{1}\mathbf{W}\|<1$.
If we choose the matrix row sum norm $\|\cdot\|_{\infty}$ and assume that
$\mathbf{W}$ is row normalized, then (2) reduces to
$(|\rho_{2}|+|\rho_{3}|)/(1-|\rho_{1}|)$. Thus, a further restrictive
sufficient condition for the stability of (2.4) is
$|\rho_{1}|+|\rho_{2}|+|\rho_{3}|<1$. We will impose these restrictions during
the sampling steps for $\rho_{1}$, $\rho_{2}$ and $\rho_{3}$ in our suggested
Gibbs sampler.

Finally, to complete the model in (2.1), we assume the following prior
distributions for the posterior analysis:

$\displaystyle\rho_{1}\sim\text{Uniform}(-1,1),\quad\rho_{2}\sim\text{Uniform}(-1,1),\,\rho_{3}\sim\text{Uniform}(-1,1),$
(2.7) $\displaystyle\boldsymbol{\mu}|\mathbf{b}_{\mu},\mathbf{B}_{\mu}\sim
N(\mathbf{b}_{\mu},\mathbf{B}_{\mu}),\,\sigma^{2}|a,b\sim\text{IG}(a,b),$

where $\text{Uniform}(c_{1},\,c_{2})$ denotes the uniform distribution over
the interval $(c_{1},\,c_{2})$ and $\text{IG}(a,b)$ denotes the inverse gamma
distribution with the shape parameter $a$ and the scale parameter $b$. The
priors for $\rho_{1}$, $\rho_{2}$ and $\rho_{3}$ are subject to the stability
conditions sated in (2.5).

##  3 Statistical Properties

The outcome equation of our model can be written as

$\displaystyle y_{it}=e^{\frac{1}{2}h_{it}}v_{it}\qquad\text{for all
$i=1,\ldots,n$ and $t=1,\ldots,T$.}$ (3.1)

Thus, the conditional variance of $y_{it}$ given $h_{it}$ is
$\text{Var}(y_{it}|h_{it})=e^{h_{it}}$, indicating that the conditional
variance is both time and space varying. Following the time series literature,
we refer to $h_{it}$ as the log-volatility since
$h_{it}=\log\left(\text{Var}(y_{it}|h_{it})\right)$. In order to determine the
unconditional moments of $y_{it}$, we need to determine the distribution of
$\mathbf{h}=(\mathbf{h}^{{}^{\prime}}_{1},\ldots,\mathbf{h}^{{}^{\prime}}_{T})^{{}^{\prime}}$.
Let $\boldsymbol{\rho}=(\rho_{1},\rho_{2},\rho_{3})^{{}^{\prime}}$, and define
the $nT\times nT$ matrix $\mathbf{J}(\boldsymbol{\rho})$ as

$\displaystyle\mathbf{J}(\boldsymbol{\rho})=\begin{pmatrix}\mathbf{S}(\rho_{1})&\mathbf{0}&\ldots&\mathbf{0}&\mathbf{0}\\\
-\mathbf{A}(\rho_{2},\rho_{3})&\mathbf{S}(\rho_{1})&\ldots&\mathbf{0}&\mathbf{0}\\\
\vdots&\ddots&\ddots&\vdots&\vdots\\\
\mathbf{0}&\mathbf{0}&\ldots&-\mathbf{A}(\rho_{2},\rho_{3})&\mathbf{S}(\rho_{1})\end{pmatrix}.$
(3.2)

Then, we can express the log-volatility equation as

$\displaystyle\mathbf{J}(\boldsymbol{\rho})(\mathbf{h}-\boldsymbol{l}_{T}\otimes\boldsymbol{\mu})=\begin{pmatrix}\mathbf{S}(\rho_{1})(\mathbf{h}_{1}-\boldsymbol{\mu})\\\
\mathbf{U}_{2}\\\ \vdots\\\ \mathbf{U}_{T}\\\ \end{pmatrix}.$ (3.3)

We assume that the spatial dynamic process for $\mathbf{h}_{t}$ has been
operating for a long time so that we can express
$\mathbf{S}(\rho_{1})(\mathbf{h}_{1}-\boldsymbol{\mu})$ in the following way:

$\displaystyle\mathbf{S}(\rho_{1})(\mathbf{h}_{1}-\boldsymbol{\mu})=\sum_{j=0}^{\infty}\left(\mathbf{A}(\rho_{2},\rho_{3})\mathbf{S}^{-1}(\rho_{1})\right)^{j}\mathbf{U}_{1-j},$
(3.4)

which implies that

$\displaystyle\mathrm{Var}(\mathbf{S}(\rho_{1})(\mathbf{h}_{1}-\boldsymbol{\mu}))$
$\displaystyle=\sigma^{2}\sum_{j=0}^{\infty}\left(\mathbf{A}(\rho_{2},\rho_{3})\mathbf{S}^{-1}(\rho_{1})\right)^{j}\left(\mathbf{A}(\rho_{2},\rho_{3})\mathbf{S}^{-1}(\rho_{1})\right)^{{}^{\prime}j}=\sigma^{2}\mathbf{K}(\boldsymbol{\rho}),$

where
$\mathbf{K}(\boldsymbol{\rho})=\sum_{j=0}^{\infty}\left(\mathbf{A}(\rho_{2},\rho_{3})\mathbf{S}^{-1}(\rho_{1})\right)^{j}\left(\mathbf{A}(\rho_{2},\rho_{3})\mathbf{S}^{-1}(\rho_{1})\right)^{{}^{\prime}j}$.
Then, from (3.3), we obtain

$\displaystyle\mathrm{Var}\left((\mathbf{h}-\boldsymbol{l}_{T}\otimes\boldsymbol{\mu})\right)=\sigma^{2}\mathbf{J}^{-1}(\boldsymbol{\rho})\mathbf{P}(\boldsymbol{\rho})\mathbf{J}^{{}^{\prime}-1}(\boldsymbol{\rho}),$
(3.5)

where

$\displaystyle\mathbf{P}(\boldsymbol{\rho})=\begin{pmatrix}\mathbf{K}(\boldsymbol{\rho})&\mathbf{0}&\ldots&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&\mathbf{I}_{n}&\ldots&\mathbf{0}&\mathbf{0}\\\
\vdots&\vdots&\ddots&\vdots&\vdots\\\
\mathbf{0}&\mathbf{0}&\ldots&\mathbf{0}&\mathbf{I}_{n}\end{pmatrix}.$ (3.6)

Let
$\boldsymbol{\Omega}=\sigma^{2}\mathbf{J}^{-1}(\boldsymbol{\rho})\mathbf{P}(\boldsymbol{\rho})\mathbf{J}^{{}^{\prime}-1}(\boldsymbol{\rho})$.
Then, the distribution of $\mathbf{h}$ is

$\displaystyle\mathbf{h}|\boldsymbol{\rho},\boldsymbol{\mu},\sigma^{2}\sim
N(\boldsymbol{l}_{T}\otimes\boldsymbol{\mu},\,\boldsymbol{\Omega}).$ (3.7)

Note that when $\boldsymbol{\rho}=\mathbf{0}$, $\boldsymbol{\Omega}$ reduces
to $\sigma^{2}\mathbf{I}_{nT}$, and thus the result in (3.7) becomes
$\mathbf{h}|\boldsymbol{\mu},\sigma^{2}\sim
N(\boldsymbol{l}_{T}\otimes\boldsymbol{\mu},\,\sigma^{2}\mathbf{I}_{nT})$.
Consider the following partition of $\boldsymbol{\Omega}$:

$\displaystyle\boldsymbol{\Omega}=\begin{pmatrix}\boldsymbol{\Omega}_{11}&\boldsymbol{\Omega}_{12}&\ldots&\boldsymbol{\Omega}_{1,T-1}&\boldsymbol{\Omega}_{1T}\\\
\boldsymbol{\Omega}_{21}&\boldsymbol{\Omega}_{22}&\ldots&\boldsymbol{\Omega}_{2,T-1}&\boldsymbol{\Omega}_{2T}\\\
\vdots&\vdots&\ddots&\vdots&\vdots\\\
\boldsymbol{\Omega}_{T1}&\boldsymbol{\Omega}_{T2}&\ldots&\boldsymbol{\Omega}_{T,T-1}&\boldsymbol{\Omega}_{TT}\\\
\end{pmatrix},$ (3.8)

where each $\boldsymbol{\Omega}_{st}$ for $s,t=1,2\ldots,T$ is an $n\times n$
sub-matrix of $\boldsymbol{\Omega}$. Let $\Omega_{ij,st}$ be the $(i,j)$th
element of $\boldsymbol{\Omega}_{st}$ for $i,j=1,2,\ldots,n$. Let
$r\in\mathbb{N}$ be a natural even number. Then, the even moments of $y_{it}$
can be expressed as

$\displaystyle\mathbb{E}\left(y^{r}_{it}\right)=\mathbb{E}\left(e^{\frac{r}{2}h_{it}}\right)\mathbb{E}(v^{r}_{it})=\exp\left(\frac{\mu_{i}r}{2}+\frac{r^{2}}{8}\Omega_{ii,tt}\right)\gamma(r)$
(3.9)

where $\gamma(r)=\frac{r!}{2^{r/2}(r/2)!}$. Then, it follows that
$\mathbb{E}(y^{4}_{it})/\left(\mathbb{E}(y^{2}_{it})\right)^{2}-3=3\left(\exp(\Omega_{ii,tt})-1\right)>0$.
Thus, our specification suggests that $y_{it}$ has a leptokurtic symmetric
distribution. Next, we consider the covariance between $y^{r}_{it}$ and
$y^{r}_{js}$:

$\displaystyle\text{Cov}\left(y^{r}_{it},y^{r}_{js}\right)$
$\displaystyle=\mathbb{E}\left(e^{\frac{r}{2}(h_{it}+h_{js})}v^{r}_{it}v^{r}_{js}\right)-\mathbb{E}\left(e^{\frac{r}{2}h_{it}}v^{r}_{it}\right)\mathbb{E}\left(e^{\frac{r}{2}h_{js}}v^{r}_{js}\right)$
$\displaystyle=\gamma^{2}(r)\exp\left(\frac{r(\mu_{i}+\mu_{j})}{2}+\frac{r^{2}}{8}\left(\Omega_{ii,tt}+\Omega_{jj,ss}+2\Omega_{ij,ts}\right)\right)$
$\displaystyle-\gamma^{2}(r)\exp\left(\frac{\mu_{i}r}{2}+\frac{r^{2}}{8}\Omega_{ii,tt}\right)\exp\left(\frac{\mu_{j}r}{2}+\frac{r^{2}}{8}\Omega_{jj,ss}\right)$
$\displaystyle=\gamma^{2}(r)\exp\left(\frac{r(\mu_{i}+\mu_{j})}{2}+\frac{r^{2}}{8}\left(\Omega_{ii,tt}+\Omega_{jj,ss}\right)\right)\left(\exp\left(\frac{r^{2}}{4}\Omega_{ij,ts}\right)-1\right).$
(3.10)

This result indicates that our specification introduces spatial dependence in
the outcome variable, since $\text{Cov}\left(y^{r}_{it},y^{r}_{js}\right)\neq
0$ in general. Note that $\text{Cov}\left(y^{r}_{it},y^{r}_{js}\right)=0$ when
$\boldsymbol{\rho}=\mathbf{0}$ holds, because
$\left(\exp\left(\frac{r^{2}}{4}\Omega_{ij,ts}\right)-1\right)=0$.

##  4 Posterior Analysis

To introduce a Bayesian MCMC estimation approach, we first transform our model
such that the resulting outcome equation is linear in $\mathbf{h}_{t}$. We
then determine the conditional likelihood function of the transformed model by
approximating to the distribution of transformed disturbance term with a
Gaussian mixture distribution (Kim et al.,, 1998; Chib et al.,, 2002; Omori et
al.,, 2007). The conditional likelihood function of the transformed model
facilitates the sampling steps for $\mathbf{h}_{t}$ and the auxiliary mixture
component indicator defined subsequently. We also provide the conditional
likelihood function of the original model, which we use to determine the
sampling steps of other parameters in our model.

We square both sides of (3.1) and then take the logarithm to obtain

$\displaystyle y^{*}_{it}=h_{it}+v^{*}_{it},$ (4.1)

where $y^{*}_{it}=\log y^{2}_{it}$ and $v^{*}_{it}=\log v^{2}_{it}$. The
density of $v^{*}_{it}$ is highly skewed with a long tail on the left and can
be expressed as

$\displaystyle
p(v^{*}_{it})=\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{1}{2}(e^{v^{*}_{it}}-v^{*}_{it})\right),\quad-\infty<v^{*}_{it}<\infty,\quad
i=1,2,\ldots,n,\,t=1,\ldots,n.$ (4.2)

It can be shown that $\mathbb{E}(v^{*}_{it})\approx-1.2704$ and
$\text{Var}(v^{*}_{it})=\pi^{2}/2\approx 4.9348$. Define
$\mathbf{Y}^{*}_{t}=(y^{*}_{1t},y^{*}_{2t},\ldots,y^{*}_{nt})^{{}^{\prime}}$
and
$\mathbf{V}^{*}_{t}=(v^{*}_{1t},v^{*}_{2t},\ldots,v^{*}_{nt})^{{}^{\prime}}$.
Then, in vector form, we have

$\displaystyle\mathbf{Y}^{*}_{t}=\mathbf{h}_{t}+\mathbf{V}^{*}_{t},$ (4.3)

In order to convert (4.3) into a linear Gaussian state-space model, we
approximate $p(v^{*}_{it})$ with an $m$-component Gaussian mixture
distribution:

$\displaystyle
p(v^{*}_{it})\approx\sum_{j=1}^{m}p_{j}\times\phi(v^{*}_{it}|\mu_{j},\,\sigma^{2}_{j}),$
(4.4)

where $\phi(v^{*}_{it}|\mu_{j},\,\sigma^{2}_{j})$ denotes the Gaussian density
function with mean $\mu_{j}$ and variance $\sigma^{2}_{j}$, $p_{j}$ is the
probability of $j$th mixture component and $m$ is the number of components. In
particular, we use the ten-component Gaussian mixture distribution suggested
by Omori et al., (2007) to approximate $p(v^{*}_{it})$. We provide the
parameter values of the ten-component Gaussian mixture distribution in Table
1. The parameters in this table are chosen by matching the first four moments
of the ten component Gaussian mixture distribution with that of
$p(v^{*}_{it})$. This approach has two advantages. First, the Gaussian mixture
distribution with the pre-determined parameter values in Table 1 provides a
well enough approximation to $p(v^{*}_{it})$ (Omori et al.,, 2007). Second,
this approach does not pose any estimation difficulties since the mixture
parameters in Table 1 are pre-determined.

We can equivalently write (4.4) in terms of an auxiliary discrete random
variable $z_{it}\in\\{1,2,\ldots,m\\}$ that serves as the mixture component
indicator:

$\displaystyle v^{*}_{it}|(z_{it}=j)\sim
N(\mu_{j},\,\sigma^{2}_{j}),\quad\text{and}\quad\mathbb{P}(z_{it}=j)=p_{j},\quad
j=1,2,\ldots,m,$ (4.5)

where $\mathbb{P}(z_{it}=j)=p_{j}$ is the probability that $z_{it}$ takes the
$j$th value. Let $\mathbf{Z}_{t}=(z_{1t},\ldots,z_{nt})^{{}^{\prime}}$,
$\mathbf{d}_{t}=(\mu_{z_{1t}},\ldots,\mu_{z_{nt}})^{{}^{\prime}}$ and
$\boldsymbol{\Sigma}_{t}=\operatorname{Diag}(\sigma^{2}_{z_{1t}},\ldots,\sigma^{2}_{z_{nt}})$.
Then, from (4.5), we have $\mathbf{V}^{*}_{t}|\mathbf{Z}_{t}\sim
N(\mathbf{d}_{t},\,\boldsymbol{\Sigma}_{t})$, which indicates that our model
in (4.3) is now conditionally linear Gaussian given the component indicator
variable. Thus, from (4.3), we have

$\displaystyle\mathbf{Y}^{*}_{t}|\mathbf{Z}_{t},\,\mathbf{h}_{t}\sim
N\left(\mathbf{h}_{t}+\mathbf{d}_{t},\,\boldsymbol{\Sigma}_{t}\right),$ (4.6)

which facilitates the sampling steps for $\mathbf{h}_{t}$ and $\mathbf{Z}_{t}$
in our suggested Gibbs sampler given in Algorithm 1. The sampling steps for
the remaining parameters requires the following conditional distribution:

$\displaystyle\mathbf{Y}_{t}|\mathbf{h}_{t}\sim
N\left(\mathbf{0},\,\mathbf{H}_{t}\right).$ (4.7)

Table 1: The ten-component Gaussian mixture for $p(v^{*}_{it})$ Components | $p_{j}$ | $\mu_{j}$ | $\sigma^{2}_{j}$
---|---|---|---
1 | 0.00609 | 1.92677 | 0.11265
2 | 0.04775 | 1.34744 | 0.17788
3 | 0.13057 | 0.73504 | 0.26768
4 | 0.20674 | 0.02266 | 0.40611
5 | 0.22715 | -0.85173 | 0.62699
6 | 0.18842 | -1.97278 | 0.98583
7 | 0.12047 | -3.46788 | 1.57469
8 | 0.05591 | -5.55246 | 2.54498
9 | 0.01575 | -8.68384 | 4.16591
10 | 0.00115 | -14.65000 | 7.33342

We are now in a position to design a Gibbs sampler by using our results on (i)
the mixture component indicators in (4.5), (ii) the conditional likelihood
function of transformed model in (4.6), (iii) the conditional likelihood
function of $\mathbf{Y}_{t}$ in (4.7) and (iv) the distribution of
$\mathbf{h}$ in (3.7). Let
$\mathbf{Y}=(\mathbf{Y}^{{}^{\prime}}_{1},\ldots,\mathbf{Y}^{{}^{\prime}}_{T})^{{}^{\prime}}$
and
$\mathbf{Z}=(\mathbf{Z}^{{}^{\prime}}_{1},\ldots,\mathbf{Z}^{{}^{\prime}}_{T})^{{}^{\prime}}$.
The joint posterior distribution
$p(\mathbf{h},\mathbf{Z},\boldsymbol{\mu},\boldsymbol{\rho},\sigma^{2}|\mathbf{Y})$
can be expressed as

$\displaystyle
p(\mathbf{h},\mathbf{Z},\boldsymbol{\mu},\boldsymbol{\rho},\sigma^{2}|\mathbf{Y})$
$\displaystyle\propto p(\mathbf{Y}^{*}|\mathbf{Z},\mathbf{h})\times
p(\mathbf{h}|\boldsymbol{\rho},\boldsymbol{\mu},\sigma^{2})\times
p(\mathbf{Z})\times p(\boldsymbol{\mu})\times p(\boldsymbol{\rho})\times
p(\sigma^{2}).$ (4.8)

Then, our suggested Gibbs sampler for generating draws from
$p(\mathbf{h},\mathbf{Z},\boldsymbol{\mu},\boldsymbol{\rho},\sigma^{2}|\mathbf{Y})$
consists of the steps given in Algorithm 1.

######  Algorithm 1 (Estimation Algorithm).

  1. 1.

Sampling step for $\mathbf{Z}$: Note that $\mathbf{Z}_{t}$ is a discrete
random variable, and its conditional posterior probability mass function is

$\displaystyle
p(\mathbf{Z}_{t}|\mathbf{Y}_{t},\mathbf{h}_{t},\boldsymbol{\mu},\boldsymbol{\rho},\sigma^{2})\propto
p(\mathbf{Z}_{t})p(\mathbf{Y}^{*}_{t}|\mathbf{Z}_{t},\mathbf{h}_{t})=\prod_{i=1}^{n}p(y^{*}_{it}|z_{it},h_{it})p(z_{it}),$
(4.9)

for $t=1,\ldots,T$. Thus,

$\displaystyle\mathbb{P}(z_{it}=j|y^{*}_{it})=\frac{p_{j}\phi(y^{*}_{it}|h_{it}+\mu_{j},\,\sigma^{2}_{j})}{\sum_{k=1}^{10}p_{k}\phi(y^{*}_{it}|h_{it}+\mu_{k},\,\sigma^{2}_{k})},\,j=1,\ldots,10,\,i=1,\ldots,n,$
(4.10)

for $t=1,\ldots,T$, where the denominator is the normalization constant.

  2. 2.

Sampling step for $\mathbf{h}$: Let
$\boldsymbol{\Sigma}=\operatorname{Diag}\left(\boldsymbol{\Sigma}_{1},\ldots,\boldsymbol{\Sigma}_{T}\right)$
and
$\mathbf{d}=(\mathbf{d}^{{}^{\prime}}_{1},\ldots,\mathbf{d}^{{}^{\prime}}_{T})^{{}^{\prime}}$.
Using standard regression results on
$p(\mathbf{h}|\mathbf{Y},\mathbf{Z},\boldsymbol{\mu},\boldsymbol{\rho},\sigma^{2})\propto\prod_{t=1}^{T}p(\mathbf{Y}^{*}_{t}|\mathbf{Z},\mathbf{h}_{t})p(\mathbf{h}|\boldsymbol{\rho},\boldsymbol{\mu},\sigma^{2})$,
we obtain

$\displaystyle\mathbf{h}|\mathbf{Y},\mathbf{Z},\boldsymbol{\mu},\boldsymbol{\rho},\sigma^{2}\sim
N(\widehat{\text{$\mathbf{b}$}}_{h},\widehat{\text{$\mathbf{B}$}}_{h}),$
(4.11)

where

$\displaystyle\widehat{\text{$\mathbf{B}$}}_{h}=\left(\boldsymbol{\Omega}^{-1}+\boldsymbol{\Sigma}^{-1}\right)^{-1},\quad\widehat{\text{$\boldsymbol{b}$}}_{h}=\widehat{\text{$\mathbf{B}$}}_{h}\left(\boldsymbol{\Omega}^{-1}(\boldsymbol{l}_{T}\otimes\boldsymbol{\mu})+\boldsymbol{\Sigma}^{-1}(\mathbf{Y}^{*}-\mathbf{d})\right).$

  3. 3.

Sampling step for $\boldsymbol{\mu}$: Using (3.2) and (3.6), we can find that

$\displaystyle\boldsymbol{\Omega}^{-1}=\begin{pmatrix}\boldsymbol{\Omega}^{*}_{11}&\boldsymbol{\Omega}^{*}_{12}&\mathbf{0}&\ldots&\mathbf{0}&\mathbf{0}\\\
\boldsymbol{\Omega}^{*}_{21}&\boldsymbol{\Omega}^{*}_{22}&\boldsymbol{\Omega}^{*}_{23}&\ldots&\mathbf{0}&\mathbf{0}\\\
\vdots&\ddots&\ddots&\ddots&\vdots&\vdots\\\
\mathbf{0}&\mathbf{0}&\mathbf{0}&\ldots&\boldsymbol{\Omega}^{*}_{T-1,T-1}&\boldsymbol{\Omega}^{*}_{T-1,T}\\\
\mathbf{0}&\mathbf{0}&\mathbf{0}&\ldots&\boldsymbol{\Omega}^{*}_{T,T-1}&\boldsymbol{\Omega}^{*}_{TT}\end{pmatrix},$

where

$\displaystyle\boldsymbol{\Omega}^{*}_{11}=\sigma^{-2}\mathbf{S}^{{}^{\prime}}(\rho_{1})\mathbf{K}^{-1}(\boldsymbol{\rho})\mathbf{S}(\rho_{1})+\sigma^{-2}\mathbf{A}^{{}^{\prime}}(\rho_{2},\rho_{3})\mathbf{A}(\rho_{2},\rho_{3}),\quad\boldsymbol{\Omega}^{*}_{TT}=\sigma^{-2}\mathbf{S}^{{}^{\prime}}(\rho_{1})\mathbf{S}(\rho_{1}),$
$\displaystyle\boldsymbol{\Omega}^{*}_{ii}=\sigma^{-2}\mathbf{S}^{{}^{\prime}}(\rho_{1})\mathbf{S}(\rho_{1})+\sigma^{-2}\mathbf{A}^{{}^{\prime}}(\rho_{2},\rho_{3})\mathbf{A}(\rho_{2},\rho_{3}),\quad
i=2,\ldots,T-1,$
$\displaystyle\boldsymbol{\Omega}^{*}_{i,i+1}=\boldsymbol{\Omega}^{*^{\prime}}_{i+1,i}=-\sigma^{-2}\mathbf{A}^{{}^{\prime}}(\rho_{2},\rho_{3})\mathbf{S}(\rho_{1}),\quad
i=1,\ldots,T-1.$

Then, from
$p(\boldsymbol{\mu}|\mathbf{Y},\mathbf{h},\mathbf{Z},\boldsymbol{\mu},\boldsymbol{\rho},\sigma^{2})\propto
p(\mathbf{h}|\boldsymbol{\rho},\boldsymbol{\mu},\sigma^{2})p(\boldsymbol{\mu})$,
we obtain

$\displaystyle\boldsymbol{\mu}|\mathbf{Y},\mathbf{h},\mathbf{Z},\boldsymbol{\mu},\boldsymbol{\rho},\sigma^{2}\sim
N(\widehat{\text{$\mathbf{b}$}}_{\mu},\widehat{\text{$\mathbf{B}$}}_{\mu}),$
(4.12)

where

$\widehat{\text{$\mathbf{B}$}}_{\mu}=\left(\mathbf{B}^{-1}_{\mu}+\sum_{j=1}^{T}\sum_{i=1}^{T}\boldsymbol{\Omega}^{*}_{ij}\right)^{-1},\quad\widehat{\text{$\mathbf{b}$}}_{\mu}=\widehat{\text{$\mathbf{B}$}}_{\mu}\left(\mathbf{B}^{-1}_{\mu}\boldsymbol{b}_{\mu}+\sum_{j=1}^{T}\sum_{i=1}^{T}\boldsymbol{\Omega}^{*}_{ij}\mathbf{h}_{i}\right)$

  4. 4.

Sampling step for $\sigma^{2}$: From
$p(\sigma^{2}|\mathbf{Y},\mathbf{h},\boldsymbol{\mu},\boldsymbol{\rho})\propto
p(\mathbf{h}|\boldsymbol{\rho},\boldsymbol{\mu},\sigma^{2})p(\sigma^{2})$, we
obtain

$\displaystyle\sigma^{2}|\mathbf{Y},\mathbf{h},\mathbf{Z},\boldsymbol{\mu},\boldsymbol{\rho}\sim\text{IG}(\widehat{\text{$a$}},\,\widehat{\text{$b$}}),$

where

$\widehat{\text{$a$}}=a+nT/2,\quad\widehat{\text{$b$}}=b+\frac{1}{2}\left(\mathbf{h}-(\boldsymbol{l}_{T}\otimes\boldsymbol{\mu})\right)^{{}^{\prime}}\left(\mathbf{J}^{{}^{\prime}}(\boldsymbol{\rho})\mathbf{P}^{-1}(\boldsymbol{\rho},\boldsymbol{\lambda})\mathbf{J}(\boldsymbol{\rho})\right)\left(\mathbf{h}-(\boldsymbol{l}_{T}\otimes\boldsymbol{\mu})\right).$

  5. 5.

Sampling step for $\boldsymbol{\rho}$: The conditional posterior density of
$\boldsymbol{\rho}$ does not take any known form in our model. We use the
adaptive Metropolis (AM) algorithm suggested in Haario et al., (2001) and
Roberts and Rosenthal, (2009) to generate draws from
$p(\boldsymbol{\rho}|\mathbf{Y},\mathbf{h},\mathbf{Z},\sigma^{2})$.222Han and
Lee, (2016) and Han et al., (2017) use this algorithm to generate draws for
the spatial parameters in spatial panel data models. Their results show that
this algorithm can perform satisfactorily. At the iteration $g$, we use the
following proposal distribution to generate the candidate value
$\tilde{\boldsymbol{\rho}}$:

$\displaystyle
f_{g}\left(\boldsymbol{\rho}|\boldsymbol{\rho}^{(0)},\ldots,\boldsymbol{\rho}^{(g-1)}\right)=\begin{cases}N\left(\boldsymbol{\rho}^{(g-1)},\,\frac{(0.1)^{2}}{3}\times\mathbf{I}_{3}\right),\quad\text{for}\quad
g\leq g_{0},\\\ 0.95\times
N\left(\boldsymbol{\rho}^{(g-1)},\,\frac{c(2.38)^{2}}{3}\times\text{Cov}\left(\boldsymbol{\rho}^{(0)},\ldots,\boldsymbol{\rho}^{(g-1)}\right)\right)\\\
\quad\quad+0.05\times
N\left(\boldsymbol{\rho}^{(g-1)},\,\frac{(0.1)^{2}}{3}\times\mathbf{I}_{3}\right),\quad\text{for}\quad
g>g_{0},\end{cases}$

where $g_{0}$ is the length of initial sampling period,
$\text{Cov}\left(\boldsymbol{\rho}^{(0)},\ldots,\boldsymbol{\rho}^{(g-1)}\right)$
is the empirical covariance matrix of historical draws given by
$\text{Cov}\left(\boldsymbol{\rho}^{(0)},\ldots,\boldsymbol{\rho}^{(g-1)}\right)=\frac{1}{g}\sum_{j=0}^{g-1}\boldsymbol{\rho}^{(j)}\boldsymbol{\rho}^{(j)^{\prime}}-\hbox
to0.0pt{\hskip
0.08507pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$\boldsymbol{\rho}$}}^{(g-1)}\hbox
to0.0pt{\hskip
0.08507pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$\boldsymbol{\rho}$}}^{(g-1)^{\prime}}$
with $\hbox to0.0pt{\hskip
0.08507pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$\boldsymbol{\rho}$}}^{(g-1)}=\frac{1}{g}\sum_{j=0}^{g-1}\boldsymbol{\rho}^{(j)}$,
and $c$ is a scalar tuning parameter used to achieve a reasonable acceptance
rate. We then check whether $\tilde{\boldsymbol{\rho}}$ satisfies the
stability conditions in Section 2. If not, we regenerate
$\tilde{\boldsymbol{\rho}}$ until it meets the stability conditions. We
compute the following acceptance probability:

$\displaystyle\mathbb{P}(\boldsymbol{\rho}^{(g-1)},\tilde{\boldsymbol{\rho}})=\min\left(\frac{p\left(\mathbf{h}|\tilde{\boldsymbol{\rho}},\boldsymbol{\mu}^{(g)},\sigma^{2(g)}\right)}{p\left(\mathbf{h}|\boldsymbol{\rho}^{(g-1)},\boldsymbol{\mu}^{(g)},\sigma^{2(g)}\right)},\,1\right),$

where $p(\mathbf{h}|\boldsymbol{\rho},\boldsymbol{\mu},\sigma^{2})$ is given
in (3.7). Finally, we return $\tilde{\boldsymbol{\rho}}$ with probability
$\mathbb{P}(\boldsymbol{\rho}^{(g-1)},\tilde{\boldsymbol{\rho}})$; otherwise
return $\boldsymbol{\rho}^{(g-1)}$.

######  Remark 1.

In Step 1, the result in (4.10) indicates that the mixture components are
conditionally independent given $y^{*}_{it}$. Thus, each component is a
discrete random variable taking integer values in the interval $[1,10]$ with
the conditional posterior probability $\mathbb{P}(z_{it}=j|y^{*}_{it})$. The
conditional posterior results in Steps 2, 3 and 4 are obtained from a standard
Bayesian analysis as in a linear regression model. In the AM algorithm
described in Step 5, the proposal distribution
$f_{g}\left(\boldsymbol{\rho}|\boldsymbol{\rho}^{(0)},\ldots,\boldsymbol{\rho}^{(g-1)}\right)$
has two parts. The first part is
$N\left(\boldsymbol{\rho}^{(g-1)},\,\frac{(0.1)^{2}}{3}\times\mathbf{I}_{3}\right)$,
and is used when the number of iterations is less than or equal to $g_{0}$.
The second part consists of two normal distributions. The first component is
specified as
$N\left(\boldsymbol{\rho}^{(g-1)},\,\frac{c(2.38)^{2}}{3}\times\text{Cov}\left(\boldsymbol{\rho}^{(0)},\ldots,\boldsymbol{\rho}^{(g-1)}\right)\right)$,
where the covariance matrix is determined from the historical MCMC draws of
$\boldsymbol{\rho}$. The second component is
$N\left(\boldsymbol{\rho}^{(g-1)},\,\frac{(0.1)^{2}}{3}\times\mathbf{I}_{3}\right)$.
The candidate values generated from
$f_{g}\left(\boldsymbol{\rho}|\boldsymbol{\rho}^{(0)},\ldots,\boldsymbol{\rho}^{(g-1)}\right)$
are subject to the stability conditions given in (2.4). Finally, we adjust the
tuning parameter $c$ during the estimation to achieve an acceptance rate that
falls between 40 percent and 60 percent.

######  Remark 2.

In the sampling step for $\boldsymbol{\rho}$,
$\mathbb{P}(\boldsymbol{\rho}^{(g-1)},\tilde{\boldsymbol{\rho}})$ is
calculated at each pass of the sampler, and therefore,
$p(\mathbf{h}|\boldsymbol{\rho},\boldsymbol{\mu},\sigma^{2})$ is evaluated
twice at each pass of the sampler. In other words,
$p(\mathbf{h}|\boldsymbol{\rho},\boldsymbol{\mu},\sigma^{2})=(2\pi)^{-nT/2}|\boldsymbol{\Omega}|^{-1/2}\exp\left(-\frac{1}{2}(\mathbf{h}-\boldsymbol{l}_{T}\otimes\boldsymbol{\mu})^{{}^{\prime}}\boldsymbol{\Omega}^{-1}(\mathbf{h}-\boldsymbol{l}_{T}\otimes\boldsymbol{\mu})\right)$
must be calculated twice. Since
$\boldsymbol{\Omega}=\sigma^{2}\mathbf{J}^{-1}(\boldsymbol{\rho})\mathbf{P}(\boldsymbol{\rho})\mathbf{J}^{{}^{\prime}-1}(\boldsymbol{\rho})$,
we have
$|\boldsymbol{\Omega}|^{-1/2}=(\sigma^{2})^{-nT/2}|\mathbf{J}(\boldsymbol{\rho})||\mathbf{P}(\boldsymbol{\rho})|^{-1/2}$.
From (3.2), since $\mathbf{J}(\boldsymbol{\rho})$ is a triangular matrix, we
have
$|\mathbf{J}(\boldsymbol{\rho})|=\prod_{t=1}^{T}|\mathbf{S}(\rho_{1})|=|\mathbf{S}(\rho_{1})|^{T}$.
Also, from (3.6), since $|\mathbf{P}(\boldsymbol{\rho})|$ is a block-diagonal
matrix, we have
$|\mathbf{P}(\boldsymbol{\rho})|=|\mathbf{K}(\boldsymbol{\rho})|$.

##  5 Simulations

In this section, we provide simulation evidence to assess sampling properties
of the suggested Bayesian algorithm. The data generating process follows (2.1)
and (2.2). More specifically, the elements of $\mathbf{V}_{t}$ and
$\mathbf{U}_{t}$ are drawn independently from the standard normal distribution
for $t=1,2,\ldots,T$. Therefore, the value of $\sigma^{2}$ is set to $0.25$ in
all experiments. The elements of $\boldsymbol{\mu}$ are drawn independently
from the normal distribution with mean $3.3$ and standard deviation $0.35$. To
initialize the process, we use (3.7), and the series expression for
$\mathbf{K}(\boldsymbol{\rho})$ is truncated at 15. We consider two sets of
values for $\boldsymbol{\rho}$, $\\{(0.6,0.35,-0.025),(0.3,0.65,-0.025)\\}$.
These parameter values are chosen to ensure that the data generating process
mimics the findings from our empirical application in the next section. The
number of spatial units $n$ is set to $98$, and the number of time periods $T$
is fixed at $50$.

For the spatial weights matrix $\mathbf{W}$, we consider row-standardized rook
and queen contiguity weights matrices. To this end, we first generate a vector
containing a random permutation of the integers from $1$ to $n$ without
repeating elements. Then, we reshape this vector into an $k\times m$
rectangular lattice, where $m=n/k$. In the case of rook contiguity, we set
$w_{ij}=1$ if the $j$th observation is adjacent (left/right/above or below) to
the $i$th observation on the lattice. In the case of queen contiguity, we set
$w_{ij}=1$ if the $j$th observation is adjacent to, or shares a border with
the $i$th observation. We set $k=7$, and row-normalize all spatial weights
matrices. For the prior distributions, we consider the following:
$\sigma^{2}\sim\text{IG}(3,\,2)$ and $\boldsymbol{\mu}\sim
N(\mathbf{0},\,10\mathbf{I}_{n})$. The length of the Markov chain is $22000$
draws, and the first $2000$ draws are discarded to dissipate the effects of
the initial values.

To determine the adequacy of the length the chains and their mixing
properties, some exemplary trace plots are provided in Figure 1 and 3. For the
sake of brevity, we only present the results for the queen contiguity case. In
these trace plots, the red solid lines correspond to the estimated posterior
means. We observe that the Bayesian estimator performs satisfactorily and
seems to mix well in all cases. Note also that for both $\boldsymbol{\rho}$
and $\sigma^{2}$, the $95\%$ credible intervals contain the true values chosen
in the experiments.

For the $n\times 1$ vector $\boldsymbol{\mu}$ and the $nT\times 1$ vector
$\mathbf{h}$, we provide evidence on the performance of our Bayesian estimator
in Figure 2 and 4. In these plots, the true values are represented with solid
lines and the estimates are presented with dashed lines. In the first panel,
we observe that the estimated posterior means for the components of
$\boldsymbol{\mu}$ are in general close to the true values. Here, the shaded
region refers to the $95\%$ credible interval. For $\mathbf{h}$, we calculate
the average of true values over $n$ and over $T$ respectively, and plot them
against the average posterior means over $n$ and over $T$, respectively. The
second panel presents the case where the average is taken over $n$, and the
last panel is the case where the average is taken over $T$. We observe that
the Bayesian estimator performs satisfactorily in terms of capturing the log-
volatility over cross-sections as well as over time.

Figure 1: Trace plots for $\boldsymbol{\rho}=(0.6,0.35,-0.025)$ and
$\mathbf{W}$ is a queen contiguity matrix

Figure 2: Estimates of $\boldsymbol{\mu}$ and $\mathbf{h}$ (dashed lines) and
their corresponding true data-generating values (solid lines) for
$\boldsymbol{\rho}=(0.6,0.35,-0.025)$ and $\mathbf{W}$ is a queen contiguity
matrix

Figure 3: Trace plots for $\boldsymbol{\rho}=(0.30,0.65,-0.025)$ and
$\mathbf{W}$ is a queen contiguity matrix

Figure 4: Estimates of $\boldsymbol{\mu}$ and $\mathbf{h}$ (dashed lines) and
their corresponding true data-generating values (solid lines) for
$\boldsymbol{\rho}=(0.30,0.65,-0.025)$ and $\mathbf{W}$ is a queen contiguity
matrix

##  6 Air Quality Modeling and Environmental Risks

In this section, we demonstrate the usage of the dynamic stochastic volatility
model using an empirical example from environmental science. Like for
applications in financial economics, the volatility of the process can be
interpreted as risks (i.e., environmental risks in our case). Previous
ecological studies mostly focused on changes in the mean behavior (cf. Wilby
et al.,, 2009), but also an increased variation might be harmful, e.g.,
Vasseur et al., (2014) showed that an increased temperature variation (i.e.,
changes in the variance) poses a greater risk than global warming (i.e.,
changes in the mean). Similar results on this topic were found by Paaijmans et
al., (2013) and Tewksbury et al., (2008). Moreover, there is also a connection
between climate changes and the financial market, such that climate variations
might have an impact on the risk of financial markets (see also Giglio et al.,
2021; Hong et al., 2020).

In this article, we analyze the variation of fine dust concentrations of
particles having a diameter less than 10 $\mu m$, $PM_{10}$, in Lombardy,
Northern Italy. The region is surrounded by the Alps from the west, so that
the wind circulation is reduced and Lombardy becomes one of the regions in
Europe with the lowest air quality (see Fassò et al.,, 2022). In the following
empirical analysis, we use the daily $PM_{10}$ concentrations from 1.1.2021 to
31.12.2021 from the official monitoring stations of the regional environmental
authority, ARPA Lombardia (Maranzano,, 2022). The data are open-source
provided by the Agrimonia project (Fassò et al., 2022a, ; Fassò et al., 2022b,
). In total, there are $n=103$ measurement stations and $T=365$ daily
observations. To provide a first overview of the dataset, we depicted the
median concentrations in $\mu g/m^{3}$ across space and time as a time-series
plot (Figure 5, left) and displayed on map (Figure 5, right), respectively.

Figure 5: Median $PM_{10}$ concentrations ($\mu g/m^{3}$) and the 5% and 95%
quantiles across space displayed as time series (left) and median
concentrations across time shown on a map (right).

For our analysis, we first estimated a spatial panel model to describe the
mean variations. The spatial correlations have been modeled in an
autoregressive manner. To be precise, the mean model is given by

$\mathbf{C}_{t}=\psi\mathbf{W}\mathbf{C}_{t}+\mathbf{X}_{t}\boldsymbol{\beta}+\mathbf{W}\mathbf{X}_{t}\boldsymbol{\gamma}+\mathbf{s}+\mathbf{E}_{t}\,.$
(6.1)

The outcome variable $\mathbf{C}_{t}$ is the $n$-dimensional vector of
$PM_{10}$ concentrations at time point $t$, $\mathbf{W}$ is spatial weights
matrix which we specify below in more detail, $\mathbf{X}_{t}$ is an $n\times
p$ matrix of exogenous regressors, and $\mathbf{E}_{t}$ denotes the model
error terms. Spatial interactions are included via the spatial autoregressive
term $\psi\mathbf{W}\mathbf{C}_{t}$ with an unknown autoregressive parameter
$\psi$. Moreover, spatial fixed effects $\mathbf{s}$ are present for each
station (even if non-significant), because there are different types of
stations included in the data set, e.g., urban traffic stations located at
major roads in the cities, or rural background stations located in the Alps.
The spatial fixed effects describe the station-specific $PM_{10}$
concentrations and thereby serve as the model intercept. We refer to this
model as Model A. Moreover, we consider the same model with spatial and
temporal fixed effects as an alternative model B, that is,

$\mathbf{C}_{t}=\psi\mathbf{W}\mathbf{C}_{t}+\mathbf{X}_{t}\boldsymbol{\beta}+\mathbf{W}\mathbf{X}_{t}\boldsymbol{\gamma}+\mathbf{s}+a_{t}\mathbf{1}_{n}+\mathbf{E}_{t}\,.$
(6.2)

The temporal fixed effects $a_{t}\mathbf{1}_{n}$ will remove any additional
(station invariant) seasonal variation across time, which is not explained by
the exogenous regressors. More precisely, we included $p=4$ covariates that
are known to influence PM concentrations, namely the maximum height of the
planetary boundary layer (PBL), the relative humidity, the air temperature,
and the pressure level. All covariates are available as daily observations for
each measurement station. Moreover, they were standardized to compare the size
of the effects. The spatial weights matrix has been chosen as row-standardized
binary contiguity matrix, where all locations within 15 miles are considered
as adjacent stations. On average, each stations has 5.2621 neighbors leading
to a sparsity level of $\mathbf{W}$ of 94.95 %. We estimated all parameters
$\\{\mathbf{s},\boldsymbol{\beta},\boldsymbol{\gamma},\psi,\sigma_{E}^{2}\\}$
using the maximum likelihood approach implemented in the spatial econometrics
MATLAB toolbox (cf. Bivand and Piras, 2015; LeSage, 1999). The resulting
estimates including their asymptotic 95% confidence intervals are summarized
in Table 2.

Table 2: Estimated parameters of the mean model |  | Model A (spatial fixed effects) | Model B (spatiotemporal fixed effects)
---|---|---|---
|  | Estimate | Standard error | 95% confidence | Estimate | Standard error | 95% confidence
|  |  |  | interval (asymptotic) |  |  | interval (asymptotic)
Average fixed effects |  |  |  |  |  |
| Spatial fixed effects | 0.0000 | 16.3766 | (-32.0975, 32.0975) | 0.0000 | 13.4717 | (-26.4040, 26.4040)
| Temporal fixed effects | / | / | / | 0.0000 | 5.8340 | (-11.4343, 11.4343)
Regressors |  |  |  |  |  |
|  $\beta_{1}$ (max PBL height) | -2.6270 | 0.1989 | (-3.0169, -2.2371) | -0.7158 | 0.2038 | (-1.1152, -0.3164)
|  $\beta_{2}$ (temperature) | -3.0972 | 0.2205 | (-3.5295, -2.6650) | -2.8055 | 0.4430 | (-3.6739, -1.9372)
|  $\beta_{3}$ (relative humidity) | -0.2110 | 0.1787 | (-0.5612, 0.1392) | 0.1272 | 0.1785 | (-0.2226, 0.4771)
|  $\beta_{4}$ (air pressure) | 20.7476 | 1.0724 | (18.6457, 22.8494) | 16.5707 | 2.0569 | (12.5392, 20.6021)
Spatially lagged regressors |  |  |  |  |  |
|  $\gamma_{1}$ (max PBL height) | 1.3812 | 0.2104 | (0.9688, 1.7935) | 0.6939 | 0.2057 | (0.2908, 1.0970)
|  $\gamma_{2}$ (temperature) | 2.3897 | 0.2270 | (1.9449, 2.8346) | 1.7405 | 0.2334 | (1.2831, 2.1980)
|  $\gamma_{3}$ (relative humidity) | 0.2534 | 0.1925 | (-0.1238, 0.6307) | 0.2444 | 0.1916 | (-0.1310, 0.6199)
|  $\gamma_{4}$ (air pressure) | -12.9406 | 1.1049 | (-15.1061, -10.7751) | -6.3344 | 1.0464 | (-8.3852, -4.2835)
Spatial dependence |  |  |  |  |  |
| $\psi$ | 0.7240 | 0.0028 | (0.7185, 0.7294) | 0.5050 | 0.0045 | (0.4962, 0.5138)
Error Variance |  |  |  |  |  |
| $\sigma_{E}^{2}$ | 62.8830 |  |  | 56.4973 |  |
| Coefficient of |  |  |  |  |  |
| determination $R^{2}$ | 0.7829 |  |  | 0.8049 |  |

The estimated parameters are in line with our expectations and we see a good
general fit of the model ($R^{2}=0.7829$ and $0.8049$ with only spatial and
spatiotemporal fixed effects, respectively). However, for this paper, our
focus is not on the mean variations but on the model errors $\mathbf{E}_{t}$,
i.e., the model uncertainties. Hence, we do not go into further detail on the
interpretation of the mean models. The errors of the mean model represent the
unexplained variations and, therefore, the environmental risks. When
considering the standard deviation of $\mathbf{E}_{t}$ at each time point
and/or for each measurement station, we see that the error variance is varying
across space and time (see Figure 6). There are time periods of increased
variations, so-called volatility clusters, e.g. at the end of the year.
Moreover, we observe a similar clustering across space. Measurement stations
with a higher volatility are located in close proximity, e.g., North-East of
Milan or around Brescia. This provides motivation for estimating a dynamic
stochastic volatility model for $\mathbf{E}_{t}$ in the second step, that is,

$\displaystyle\mathbf{E}_{t}=\mathbf{H}^{1/2}_{t}\mathbf{V}_{t}$ (6.3)

as defined in (2.1) and (2.2).

Figure 6: Standard deviation of residuals $\mathbf{E}_{t}$

Table 3 reports the median of the posterior draws for each parameter including
the corresponding 95% credible intervals (500,000 posterior draws, and 100,000
burn-in draws). The space-time interaction parameters $\rho_{1}$, $\rho_{2}$,
and $\rho_{3}$ represent the uncertainty/risk spillovers. In general, we
observe moderate instantaneous spatial interactions (i.e.,
$\widehat{\text{$\rho$}}_{1}=0.5943$ for Model A and
$\widehat{\text{$\rho$}}_{1}=0.3037$ for Model B) and temporal autoregressive
interactions (i.e., $\widehat{\text{$\rho$}}_{2}=0.3781$ for Model A and
$\widehat{\text{$\rho$}}_{2}=0.6699$ for Model B), while spatiotemporal
effects appear to be of minor importance (i.e.,
$\widehat{\text{$\rho$}}_{3}=-0.0278$ and
$\widehat{\text{$\rho$}}_{3}=-0.0236$ for Models A and B, respectively). The
posterior mean estimates for both the spatial and temporal autoregressive
parameters are positive. That is, if the environmental risk (log-volatility
$h_{it}$) is high it is likely to spillover to the neighboring regions and to
future time points. For this reason, spatial and temporal volatility clusters
are formed. The negative signs of the posterior mean estimate for the
spatiotemporal lag are opposing this behavior, but they are close to zero.
Moreover, we observe that the spatial spillovers are dominating the temporal
ones for Model A without temporal fixed effects. That is, the missing temporal
variation of mean model could be picked up within the uncertainty model. When
including both temporal and spatial fixed effects (i.e., Model B), the spatial
dependence in the volatility of the errors is reduced, while the temporal
dependence (i.e., $\rho_{2}$) is significantly higher. When looking at the
overall volatility level represented by the average $\mathbf{h}_{t}$, we see
that the model uncertainty is reduced because of the additional temporal fixed
effects in the mean equation. Nevertheless, all space-time interaction
parameters are significantly different from zero indicating a space-time
dependence in the local uncertainties of the environmental process (here: the
particulate matter concentrations).

Moreover, we analyzed the conditional volatilities for each station across the
time horizon. For Model B, the results are displayed in Figure 7 as a time-
series plot (left) and on a map (right). More precisely, we computed the
posterior medians of $h_{it}$ for all locations and time points leading to an
estimate of the $n\times T$ matrix of the logarithmic conditional
volatilities. Then, for visualization, we computed the median and the 5% and
95% quantiles of these estimates across all spatial locations (see time series
plot) and across all time points (see map). From the time series plot, we see
that there is a clear annual pattern – the model uncertainties are lower in
summer and highest in winter. These varying conditional volatility levels can
be interpreted as environmental risks, which are the highest the winter
season. In the same manner, measurement stations with a greater
uncertainty/risk can be identified from the map. They are mostly located in
valleys in the mountain areas.

Table 3: Estimated parameters of the dynamic stochastic volatility model, medians and 95% posterior credible intervals. |  | Model A | Model B
---|---|---|---
|  | Median | 95% credible interval | Median | 95% credible interval
Intercept $\boldsymbol{\mu}$ |  |  |  |
| Average constant term | 3.4326 | (3.3386, 3.5274) | 3.2774 | (3.1645, 3.3931)
Space-time interactions |  |  |  |
|  $\rho_{1}$ (spatial interaction) | 0.5943 | (0.5547, 0.6372) | 0.3037 | (0.2806, 0.3325)
|  $\rho_{2}$ (temporal interaction) | 0.3781 | (0.3353, 0.4181) | 0.6699 | (0.6406, 0.6939)
|  $\rho_{3}$ (spatiotemporal interaction) | -0.0242 | (-0.0282, -0.0189) | -0.0236 | (-0.0279, -0.0185)
Stochastic volatilities $\mathbf{h}_{t}$ |  |  |  |
| Average log-volatility | 3.4140 | (3.3972, 3.4313) | 3.2299 | (3.2134, 3.2467)
Error variance $\sigma^{2}$ |  |  |  |
|  $\sigma^{2}$ (variance of $\mathbf{U}_{t}$) | 0.2267 | (0.1996, 0.2510) | 0.3312 | (0.3110, 0.3521)

Figure 7: Median posterior draws of $\mathbf{h}_{t}$ (Model B) and the 5% and
95% quantiles across space displayed as time series (left) and median
posterior draws of $\mathbf{h}_{t}$ across time shown on a map (right).

##  7 Conclusion

We have introduced a novel spatiotemporal statistical model for stochastic
volatilities, which are spatially, temporally, and spatiotemporally correlated
and have different levels for each location reflecting the heterogeneity of
the process across space. The model includes two different error terms for the
mean and the log-volatility equations. To estimate the parameters, we
suggested a Bayesian MCMC approach. To this end, we applied a log-square
transformation to transform our non-linear state-space model into a linear
one. Further, we used a Gaussian mixture distribution to approximate the
distribution of the transformed error terms in the space equation,
transforming our model into a linear Gaussian state-space model. We analyzed
the estimation performance for different parameter settings and spatial
interactions in a simulation study, and showed that the suggested Bayesian
sampler performs satisfactorily.

Moreover, we applied the dynamic spatiotemporal stochastic volatility model in
a completely new empirical framework, namely, in the field of modeling
environmental and climate risks. First, statistical modeling of the volatility
process of climate variables have not been done extensively yet, even though
there is a large scientific consensus that an increased variability of
environmental processes, e.g., the temperature variability, is harmful for the
environment. Second, stochastic volatility models were predominantly applied
to financial data, because of the straightforward interpretation of the log
conditional volatilities as the (return) risk of financial assets. We have
transferred to this idea to environmental risks and showed how the volatility
of $PM_{10}$ predictions are correlated across space and time. In particular,
measurement stations with an increased model uncertainty could be identified
in this way. In addition, we showed in our application that there are
significant temporal and spatial spillovers in these environmental risks.
Also, the temporally lagged spatial spillovers, i.e., the spatiotemporal
correlations, appear to play a minor role.

Both from a theoretical and applied perspective, spatiotemporal stochastic
volatility models and environmental risk modeling are important topics for
future research. The current model does not allow for temporally varying
constant terms in the conditional volatilities, which are constant across
space. Regarding the latter case, other environmental processes, like
temperature variations, soil droughts, or atmospheric ozone concentration and
optical depths are important processes, where a deep understanding of the
spatiotemporal interactions in the variabilities is essential, both for
obtaining accurate prediction intervals and the planning of interventions (and
their impact on the environment).

## References

  * Anselin, (1988) Anselin, L. (1988).  Spatial econometrics: Methods and Models.  Springer, New York.
  * Arbia and Espa, (1996) Arbia, G. and Espa, G. (1996).  Effects of the maup on image classification.  Journal of Geographical Systems.
  * Bivand and Piras, (2015) Bivand, R. and Piras, G. (2015).  Comparing implementations of estimation methods for spatial econometrics.  Journal of Statistical Software, 63:1–36.
  * Chan and Jeliazkov, (2009) Chan, J. C. and Jeliazkov, I. (2009).  Efficient simulation and integrated likelihood estimation in state space models.  International Journal of Mathematical Modelling and Numerical Optimisation, 1(1-2):101–120.
  * Chan, (2017) Chan, J. C. C. (2017).  The stochastic volatility in mean model with time-varying parameters: An application to inflation modeling.  Journal of Business & Economic Statistics, 35(1):17–28.
  * Chib et al., (2002) Chib, S., Nardari, F., and Shephard, N. (2002).  Markov chain monte carlo methods for stochastic volatility models.  Journal of Econometrics, 108(2):281 – 316.
  * Cressie and Wikle, (2015) Cressie, N. and Wikle, C. K. (2015).  Statistics for spatio-temporal data.  John Wiley & Sons.
  * Elhorst, (2014) Elhorst, J. P. (2014).  Spatial Econometrics: From Cross-Sectional Data to Spatial Panels.  Springer Briefs in Regional Science. Springer Berlin Heidelberg, New York.
  * Fassò et al., (2022) Fassò, A., Maranzano, P., and Otto, P. (2022).  Spatiotemporal variable selection and air quality impact assessment of covid-19 lockdown.  Spatial Statistics, 49.
  * (10) Fassò, A., Rodeschini, J., Fusta Moro, A., Shaboviq, Q., Maranzano, P., Cameletti, M., Finazzi, F., Golini, N., Ignaccolo, R., and Otto, P. (2022a).  Agrimonia: Open access dataset correlating livestock and air quality in the lombardy region, italy.
  * (11) Fassò, A., Rodeschini, J., Moro, A. F., Shaboviq, Q., Maranzano, P., Cameletti, M., Finazzi, F., Golini, N., Ignaccolo, R., and Otto, P. (2022b).  Agrimonia: a dataset on livestock, meteorology and air quality in the lombardy region, italy.
  * Lee and Yu, (2015) Lee, L.-F. and Yu, J. (2015).  Estimation of fixed effects panel regression models with separable and nonseparable space-time filters.  Journal of Econometrics, 184(1):174 – 192.
  * Fisher, (1935) Fisher, R. A. (1935).  The design of experiments.  Oliver and Boyd, Edinburgh.
  * Giglio et al., (2021) Giglio, S., Kelly, B., and Stroebel, J. (2021).  Climate finance.  Annual Review of Financial Economics, 13:15–36.
  * Gneiting et al., (2007) Gneiting, T., Genton, M. G., and Guttorp, P. (2007).  Geostatistical space-time models, stationarity, separability and full symmetry.  In Finkenstadt, B., Held, L., and Isham, V., editors, Statistical Methods for Spatio-Temporal Systems, Monographs on Statistics and Applied Probability 107, pages 151–175. Chapman and Hall/CRC.
  * Haario et al., (2001) Haario, H., Saksman, E., and Tamminen, J. (2001).  An adaptive metropolis algorithm.  Bernoulli, 7(2):223–242.
  * Hamilton, (1994) Hamilton, J. (1994).  Time Series Analysis.  Princeton University Press.
  * Han et al., (2017) Han, X., Hsieh, C.-S., and Lee, L.-F. (2017).  Estimation and model selection of higher-order spatial autoregressive model: An efficient bayesian approach.  Regional Science and Urban Economics, 63:97 – 120.
  * Han and Lee, (2016) Han, X. and Lee, L.-F. (2016).  Bayesian analysis of spatial panel autoregressive models with time-varying endogenous spatial weight matrices, common factors, and random coefficients.  Journal of Business & Economic Statistics, 34(4):642–660.
  * Hølleland and Karlsen, (2020) Hølleland, S. and Karlsen, H. A. (2020).  A stationary spatio-temporal garch model.  Journal of Time Series Analysis, 41(2):177–209.
  * Hong et al., (2020) Hong, H., Karolyi, G. A., and Scheinkman, J. A. (2020).  Climate finance.  The Review of Financial Studies, 33(3):1011–1023.
  * Horn and Johnson, (2012) Horn, R. and Johnson, C. (2012).  Matrix Analysis.  Matrix Analysis. Cambridge University Press.
  * Huang et al., (2011) Huang, W., Wang, K., Jay Breidt, F., and Davis, R. A. (2011).  A class of stochastic volatility models for environmental applications.  Journal of Time Series Analysis, 32(4):364–377.
  * Iaco et al., (2012) Iaco, S. D., Maggio, S., Palma, M., and Posa, D. (2012).  Advances in spatio-temporal modeling and prediction for environmental risk assessment.  In Haryanto, B., editor, Air Pollution, chapter 14. IntechOpen, Rijeka.
  * Kelejian and Prucha, (2010) Kelejian, H. H. and Prucha, I. R. (2010).  Specification and estimation of spatial autoregressive models with autoregressive and heteroskedastic disturbances.  Journal of Econometrics, 157:53–67.
  * Kim et al., (1998) Kim, S., Shephard, N., and Chib, S. (1998).  Stochastic volatility: Likelihood inference and comparison with arch models.  The Review of Economic Studies, 65(3):361–393.
  * Lee, (2004) Lee, L.-f. (2004).  Asymptotic distributions of quasi-maximum likelihood estimators for spatial autoregressive models.  Econometrica, 72(6):1899–1925.
  * LeSage and Pace, (2009) LeSage, J. and Pace, R. K. (2009).  Introduction to Spatial Econometrics (Statistics: A Series of Textbooks and Monographs.  Chapman and Hall/CRC, London.
  * LeSage, (1999) LeSage, J. P. (1999).  Spatial econometrics using MATLAB.  The Web: www.econ.utoledo.edu.
  * Maranzano, (2022) Maranzano, P. (2022).  Air quality in lombardy, italy: An overview of the environmental monitoring system of arpa lombardia.  Earth, 3(1):172–203.
  * Omori et al., (2007) Omori, Y., Chib, S., Shephard, N., and Nakajima, J. (2007).  Stochastic volatility with leverage: Fast and efficient likelihood inference.  Journal of Econometrics, 140(2):425 – 449.
  * Otto et al., (2018) Otto, P., Schmid, W., and Garthoff, R. (2018).  Generalised spatial and spatiotemporal autoregressive conditional heteroscedasticity.  Spatial Statistics, 26:125–145.
  * Paaijmans et al., (2013) Paaijmans, K. P., Heinig, R. L., Seliga, R. A., Blanford, J. I., Blanford, S., Murdock, C. C., and Thomas, M. B. (2013).  Temperature variation makes ectotherms more sensitive to climate change.  Global change biology, 19(8):2373–2380.
  * Parent and LeSage, (2011) Parent, O. and LeSage, J. P. (2011).  A space-time filter for panel data models containing random effects.  Computational Statistics & Data Analysis, 55(1):475 – 490.
  * Parent and LeSage, (2012) Parent, O. and LeSage, J. P. (2012).  Spatial dynamic panel data models with random effects.  Regional Science and Urban Economics, 42(4):727–738.
  * Porcu et al., (2016) Porcu, E., Bevilacqua, M., and Genton, M. G. (2016).  Spatio-temporal covariance and cross-covariance functions of the great circle distance on a sphere.  Journal of the American Statistical Association, 111(514):888–898.
  * Roberts and Rosenthal, (2009) Roberts, G. O. and Rosenthal, J. S. (2009).  Examples of adaptive mcmc.  Journal of Computational and Graphical Statistics, 18(2):349–367.
  * Robinson, (2009) Robinson, P. M. (2009).  Large-sample inference on spatial dependence.  Econometrics Journal, 12.
  * Sato and Matsuda, (2021) Sato, T. and Matsuda, Y. (2021).  Spatial extension of generalized autoregressive conditional heteroskedasticity models.  Spatial Economic Analysis, 16(2):148–160.
  * Taşpınar et al., (2021) Taşpınar, S., Doğan, O., Chae, J., and Bera, A. K. (2021).  Bayesian inference in spatial stochastic volatility models: An application to house price returns in chicago.  Oxford Bulletin of Economics and Statistics, 83(5):1243–1272.
  * Tewksbury et al., (2008) Tewksbury, J. J., Huey, R. B., and Deutsch, C. A. (2008).  Putting the heat on tropical animals.  Science, 320(5881):1296–1297.
  * Tobler, (1970) Tobler, W. R. (1970).  A computer movie simulating urban growth in the Detroit region.  Economic geography, 46(sup1):234–240.
  * Vasseur et al., (2014) Vasseur, D. A., DeLong, J. P., Gilbert, B., Greig, H. S., Harley, C. D., McCann, K. S., Savage, V., Tunney, T. D., and O’Connor, M. I. (2014).  Increased temperature variation poses a greater risk to species than climate warming.  Proceedings of the Royal Society B: Biological Sciences, 281(1779):20132612.
  * (44) Ver Hoef, J. M., Hanks, E. M., and Hooten, M. B. (2018a).  On the relationship between conditional (CAR) and simultaneous (SAR) autoregressive models.  Spatial statistics, 25:68–85.
  * (45) Ver Hoef, J. M., Peterson, E. E., Hooten, M. B., Hanks, E. M., and Fortin, M.-J. (2018b).  Spatial autoregressive models for statistical inference from ecological data.  Ecological Monographs, 88(1):36–59.
  * Wang and Lee, (2018) Wang, W. and Lee, L.-F. (2018).  Gmm estimation of spatial panel data models with common factors and a general space–time filter.  Spatial Economic Analysis, 13(2):247–269.
  * Wilby et al., (2009) Wilby, R. L., Troni, J., Biot, Y., Tedd, L., Hewitson, B. C., Smith, D. M., and Sutton, R. T. (2009).  A review of climate risk information for adaptation and development planning.  International Journal of Climatology: A Journal of the Royal Meteorological Society, 29(9):1193–1215.
  * Yan, (2007) Yan, J. (2007).  Spatial stochastic volatility for lattice data.  Journal of Agricultural, Biological, and Environmental Statistics, 12(1):25–40.
  * Zimmerman, (2019) Zimmerman, D. L. (2019).  Geostatistics.  Wiley StatsRef: Statistics Reference Online, pages 1–10.