# First-order planar autoregressive model

Sergiy Shklyar

###### Abstract

This paper establishes the conditions of existence of a stationary solution to
the first order autoregressive equation on a plane as well as properties of
the stationarity solution. The first-order autoregressive model on a plane is
defined by the equation

$X_{i,j}=aX_{i-1,j}+bX_{i,j-1}+cX_{i-1,j-1}+\epsilon_{i,j}.$

A stationary solution $X$ to the equation exists if and only if
$(1-a-b-c)\\*(1-a+b+c)\\*(1+a-b+c)\\*(1+a+b-c)>0$. The stationary solution $X$
satisfies the causality condition with respect to the white noise $\epsilon$
if and only if $1-a-b-c>0$, $1-a+b+c>0$, $1+a-b+c>0$ and $1+a+b-c>0$. A
sufficient condition for $X$ to be purely nondeterministic is provided.

An explicit expression for the autocovariance function of $X$ at some points
is provided. With Yule–Walker equations, this allows to compute the
autocovariance function everywhere. In addition, all situations are described
where different parameters determine the same autocovariance function of $X$.

##### Keywords:

autoregressive models, causality, discrete random fields, purely
nondeterministic random fields, stationary random fields

##### MSC2020:

60G60, 62M10

##  1 Introduction

##### The model.

Let $\epsilon_{i,j}$ be a uncorrelated zero-mean equal-variance variables,
$\operatorname{\mathsf{E}}\epsilon_{i,j}=0$ and
$\operatorname{\mathsf{E}}\epsilon_{i,j}^{2}=\sigma_{\epsilon}^{2}>0$.
Consider the equation

$X_{i,j}=aX_{i-1,j}+bX_{i,j-1}+cX_{i-1,j-1}+\epsilon_{i,j}$ (1)

with $a$, $b$ and $c$ fixed coefficients, and $X$ an unknown array defined on
an integer lattice — for all integer indices.

Eq. (1) can be used for modelling and simulating random images.

##### Subject of the study.

We consider the following questions about the equation 1.

  * •

For what coefficients $a$, $b$ and $c$ does a stationary solution of the (1)
exist?

  * •

How to compute the autocovariance function of the stationary solution?

  * •

For what coefficients $a$, $b$ and $c$ does the stationary solution satisfy
the causality condition with respect to $\epsilon$, that is, can be
represented in the form

$X_{i,j}=\sum_{k=0}^{\infty}\sum_{l=0}^{\infty}\psi_{k,l}\epsilon_{i-k,j-l}$

with summable coefficients $\psi_{k,l}$?

  * •

How the causality condition is related to the stability of equation (1) taken
nondeterministically?

  * •

Is the stationary solution purely nondeterministic? This property is related
to the causality condition with respect to some uncorrelated random field,
called innovations: does the field of innovations exist?

  * •

Which different parameters $(a,b,c,\sigma_{\epsilon}^{2})$ determine the same
autocovariance function of the field $X$?

##### Previous art.

The review of statistical models for representing a discretely observed field
on a plane, including regression models and trend-stationary models, in given
in [13]. The study of planar models, their estimation and applications is
presented in the monograph [2].

A theory of stationary random field on $n$-dimensional integer lattice was
constructed by Tjøstheim [10]. For the stationary fields, he considers the
property of being “purely nondeterministic,” which is related to the existence
of appropriate field of innovations. He considers ARMA models, for which he
provides conditions for stability of the equation and existence of a
stationary solution. For AR models, he establishes Yule–Walker equations, and
also consistency and asymptotic normality of the Yule–Walker estimator.

Autoregressive models in a $n$-dimensional space had been studied before
Tjøstheim. Whittle [14] had been studied such topics of autoregressive model
in the plane as the spectral density of the stationary solution,
nonidentifiability, and estimation of the regression coefficients. In
theoretic part, the order of the model was not limited; in the example, he
models data with the second-order model of the form

$\xi_{s,t}=\alpha\xi_{s+1,t}+\beta\xi_{s-1,t}+\gamma\xi_{s,t+1}+\delta\xi_{s,t-1}+\epsilon_{s,t}.$

The model equivalent to (1) is considered in [8, 12]. One of simpler planar
models considered in the paper is defined by equation
$\xi_{s,t}=\alpha\xi_{s+1,t}+\beta\xi_{s,t+1}+\epsilon_{s,t},$ which, up to
different notation, is also a spectral case of (1).

The autocovariance function of a stationary field is even-symmetric, that is
$\gamma_{X}(h_{1},h_{2})=\gamma_{X}(-h_{1},-h_{2})$; however,
$\gamma_{X}(h_{1},h_{2})$ need not be an even function in each of its
arguments. The autocovariance function of the stationary stationary solution
to (1) is even each argument if the coefficients satisfy $c=-ab$. This special
case, its estimation and applications are considered in [1, 5, 3, 6, 7].

The property of random field of being purely nondeterministic is introduced in
[10]. Further discussion of this property is found in [4] and [11].

##### Methods.

We have obtained most results by manipulation with the spectral density.

##### Structure of the paper.

In Section 2, we provide necessary and sufficient conditions for existence of
the stationary solution to equation (1). We also prove that the stationary
solution must be centered and must have a spectral density. In Section 3, we
study the properties of the autocovariance function $\gamma_{X}(h_{1},h_{2})$
of the stationary solution. We find $\gamma_{X}(h_{1},h_{2})$ for $h_{1}=0$
and for $h_{2}=0$, and we prove Yule–Walker equations. This allows us to
evaluate $\gamma_{X}(h_{1},h_{2})$ everywhere. Section 4 is dedicated to the
case where the equation (1) is stable. In Section 5 we show
nonidentifiability, there the same autocovariance function (and, in Gaussian
case, the same distribution) of the random field $X$ corresponds to two
different values of the parameter vectors. In Section 6 we study Tjøstheim’s
pure nondeterminism. Section 7 concludes. Some auxiliary statements are
presented in the appendix, with is in Section 8.

##  2 Existence of a stationary solution

###  2.1 Spectral density

Denote the shift operators

$L_{x}X_{i,j}=X_{i-1,j}\quad\mbox{and}\quad L_{y}X_{i,j}=X_{i,j-1}.$ (2)

Eq. (1) rewrites as follows

$X-aL_{x}X-bL_{y}X-cL_{x}L_{y}X=\epsilon.$ (3)

We use the following definition of a spectral density. An integrable function
$f_{X}(\nu_{1},\nu_{2})$ is a spectral density for the wide-sense stationary
field $X$ if

$\displaystyle\gamma_{X}(h_{1},\,h_{2})$
$\displaystyle=\operatorname{cov}(X_{i,j},X_{i-h_{1},j-h_{2}})=$
$\displaystyle=\int_{-1/2}^{1/2}\int_{-1/2}^{1/2}\exp(2\pi{\mathrm{i}}\,(h_{1}\nu_{1}+h_{2}\nu_{2}))f_{X}(\nu_{1},\nu_{2})\,d\nu_{1}\,d\nu_{2},$
(4)

where the upright $\mathrm{i}$ is the imaginary unit. This can be rewritten
without complex numbers as

$\gamma_{X}(h_{1},h_{2})=\int_{-1/2}^{1/2}\int_{-1/2}^{1/2}\cos(2\pi\,(h_{1}\nu_{1}+h_{2}\nu_{2}))f_{X}(\nu_{1},\nu_{2})\,d\nu_{1}\,d\nu_{2}$
(5)

with further restraint that the spectral density should be an even function,
$f_{X}(\nu_{1},\nu_{2})=f_{X}(-\nu_{1},-\nu_{2})$.

Due to (3), the spectral densities of the random fields $\epsilon$ and $X$, if
they exist, are related as follows

$f_{\epsilon}(\nu_{1},\nu_{2})=|1-ae^{2\pi\mathrm{i}\nu_{1}}-be^{2\pi\mathrm{i}\nu_{2}}-ce^{2\pi\mathrm{i}(\nu_{1}+\nu_{2})}|^{2}\,f_{X}(\nu_{1},\nu_{2}).$
(6)

The random field $\epsilon$ is a white noise with constant spectral density
$f_{\epsilon}(\nu_{1},\nu_{2})=\sigma_{\epsilon}^{2}$. Thus, the spectral
density of the stationary field $X$, if it exists, is equal to

$f_{X}(\nu_{1},\nu_{2})=\frac{\sigma_{\epsilon}^{2}}{|1-ae^{2\pi\mathrm{i}\nu_{1}}-be^{2\pi\mathrm{i}\nu_{2}}-ce^{2\pi\mathrm{i}(\nu_{1}+\nu_{2})}|^{2}}\,.$
(7)

######  Lemma 1.

Denote

$D=(1-a-b-c)(1-a+b+c)(1+a-b+c)(1+a+b-c).$

Consider the denominator in (7)

$g(\nu_{1},\nu_{2})=|1-ae^{2\pi\mathrm{i}\nu_{1}}-be^{2\pi\mathrm{i}\nu_{2}}-ce^{2\pi\mathrm{i}(\nu_{1}+\nu_{2})}|^{2},\qquad\nu_{1},\nu_{2}\in\left(-{\textstyle\frac{1}{2}},\,{\textstyle\frac{1}{2}}\right].$

It attains zero if and only if $D\leq 0$. It can attain zero not more than at
$2$ points on the half-open box $\left(-\frac{1}{2},\frac{1}{2}\right]^{2}$.

###### Proof.

The following lines are equivalent:

$\displaystyle\exists\nu_{1}\in\left(-{\textstyle\frac{1}{2}},{\textstyle\frac{1}{2}}\right]\;\exists\nu_{2}\in\left(-{\textstyle\frac{1}{2}},{\textstyle\frac{1}{2}}\right]\>:\>|1-ae^{2\pi\mathrm{i}\nu_{1}}-be^{2\pi\mathrm{i}\nu_{2}}-ce^{2\pi\mathrm{i}(\nu_{1}+\nu_{2})}|^{2}=0,$
$\displaystyle\exists\nu_{1}\in\left(-{\textstyle\frac{1}{2}},{\textstyle\frac{1}{2}}\right]\;\exists\nu_{2}\in\left(-{\textstyle\frac{1}{2}},{\textstyle\frac{1}{2}}\right]\>:\>1-ae^{2\pi\mathrm{i}\nu_{1}}=e^{2\pi\mathrm{i}\nu_{2}}(b-ce^{2\pi\mathrm{i}\nu_{1}}),$
(8)
$\displaystyle\exists\nu_{1}\in\left(-{\textstyle\frac{1}{2}},{\textstyle\frac{1}{2}}\right]\>:\>|1-ae^{2\pi\mathrm{i}\nu_{1}}|=|(b-ce^{2\pi\mathrm{i}\nu_{1}})|,$
$\displaystyle\exists\nu_{1}\in\left(-{\textstyle\frac{1}{2}},{\textstyle\frac{1}{2}}\right]\>:\>|1-ae^{2\pi\mathrm{i}\nu_{1}}|^{2}=|(b-ce^{2\pi\mathrm{i}\nu_{1}})|^{2},$
$\displaystyle\exists\nu_{1}\in\left(-{\textstyle\frac{1}{2}},{\textstyle\frac{1}{2}}\right]\>:\>1-2a\cos(2\pi\nu_{1})+a^{2}=b^{2}-2bc\cos(2\pi\nu_{1})+c^{2},$
$\displaystyle\exists\nu_{1}\in\left(-{\textstyle\frac{1}{2}},{\textstyle\frac{1}{2}}\right]\>:\>1+a^{2}-b^{2}-c^{2}=2(a-bc)\cos(2\pi\nu_{1}),$
(9) $\displaystyle|1+a^{2}-b^{2}-c^{2}|\leq 2\,|a-bc|,$
$\displaystyle(1+a^{2}-b^{2}-c^{2})^{2}-4(a-bc)^{2}\leq 0,$ $\displaystyle
D\leq 0.$

The sufficient and necessary condition for the denominator
$g(\nu_{1},\nu_{2})$ to attain $0$ is proved.

Now treat $a$, $b$ and $c$ as fixed, and assume that $D\leq 0$. Equality in
(9) can be attained for not more than two
$\nu_{1}\in\left(-\frac{1}{2},\frac{1}{2}\right]$. For each $\nu_{1}$,
equality in (8) can be attained for no more than one
$\nu_{2}\in\left(-\frac{1}{2},\frac{1}{2}\right]$. Thus, the function
$1-ae^{2\pi\mathrm{i}\nu_{1}}-be^{2\pi\mathrm{i}\nu_{2}}-ce^{2\pi\mathrm{i}(\nu_{1}+\nu_{2})}$
has not more than 2 zeros in the half-open box
$\left(-\frac{1}{2},\frac{1}{2}\right]^{2}$. The denominator
$g(\nu_{1},\nu_{2})$ cannot attain zero at more than two points either. ∎

###  2.2 Necessary and sufficient conditions for the existence of a stationary
solution

######  Proposition 1.

Let $\sigma^{2}>0$, $a$, $b$ and $c$ be fixed real numbers. A collection of
$\epsilon=\\{\epsilon_{i,j},\;i,j\in\mathbb{Z}\\}$ of uncorrelated random
variables with zero mean and equal variance $\sigma^{2}$ and a wide-sense
stationary random field $X=\\{X_{i,j},\;i,j\in\mathbb{Z}\\}$ that satisfy
equation (1) exist if and only if $D>0$, where $D$ comes from Lemma 1.

###### Proof.

Denote

$g_{1}(\nu_{1},\nu_{2})=1-ae^{2\pi\mathrm{i}\nu_{1}}-be^{2\pi\mathrm{i}\nu_{2}}-ce^{2\pi\mathrm{i}(\nu_{1}+\nu_{2})};$

then $g(\nu_{1},\nu_{2})$ defined in Lemma 1 equals
$g(\nu_{1},\nu_{2})=|g_{1}(\nu_{1},\nu_{2})|^{2}$. Necessity. Let random
fields $\epsilon$ and $X$ that satisfy equation (1) and the other conditions
of Proposition 1. Both $\epsilon$ and $X$ are wide-sense stationary. Hence,
they have spectral measures — denote them $\lambda_{\epsilon}$ and
$\lambda_{X}$ — and $\lambda_{\epsilon}$ is a constant-weight Lebesgue measure
in $\left(-\frac{1}{2},\frac{1}{2}\right]^{2}$. Because of (3),
$\lambda_{\epsilon}$ has a Radon–Nikodym density $g(\nu_{1},\nu_{2})$ with
respect to (w.r.t.) $\lambda_{X}$:

$\frac{d\lambda_{\epsilon}}{d\lambda_{X}}(\nu_{1},\nu_{2})=|1-|1-ae^{2\pi\mathrm{i}\nu_{1}}-be^{2\pi\mathrm{i}\nu_{2}}-ce^{2\pi\mathrm{i}(\nu_{1}+\nu_{2})}|.$

Now assume that $D\leq 0$. Then, according to Lemma 1, $g(\nu_{1},\nu_{2})$
attains 0 at one or two points on a half-open box
$\left(-\frac{1}{2},\frac{1}{2}\right]^{2}$. Denote

$\displaystyle
A_{1}=\left\\{(\nu_{1},\nu_{2})\mathbin{\in}\left(-{\textstyle\frac{1}{2}},\,{\textstyle\frac{1}{2}}\right]^{2}:g_{1}(\nu_{1},\nu_{2})\neq
0\right\\}=\left\\{(\nu_{1},\nu_{2})\mathbin{\in}\left(-{\textstyle\frac{1}{2}},\,{\textstyle\frac{1}{2}}\right]^{2}:g(\nu_{1},\nu_{2})>0\right\\},$
$\displaystyle
A_{2}=\left\\{(\nu_{1},\nu_{2})\mathbin{\in}\left(-{\textstyle\frac{1}{2}},\,{\textstyle\frac{1}{2}}\right]^{2}:g_{1}(\nu_{1},\nu_{2})=0\right\\}=\left\\{(\nu_{1},\nu_{2})\mathbin{\in}\left(-{\textstyle\frac{1}{2}},\,{\textstyle\frac{1}{2}}\right]^{2}:g(\nu_{1},\nu_{2})=0\right\\};$
$\displaystyle\lambda_{1}(A)=\lambda_{X}(A\cap
A_{1}),\qquad\lambda_{2}(A)=\lambda_{X}(A\cap A_{2})$

for all measurable sets
$A\subset\left(-{\textstyle\frac{1}{2}},\,{\textstyle\frac{1}{2}}\right]^{2}$.

The set $A_{2}$ contains one or two points. As $\lambda_{\epsilon}$ is
absolutely continuous, $\lambda_{\epsilon}(A)=\lambda_{\epsilon}(A\setminus
A_{2})=\lambda_{\epsilon}(A\cap A_{1})$. Measures $\lambda_{1}$ and
$\lambda_{\epsilon}$ are absolutely continuous with respect of each other:
$d\lambda_{1}/d\lambda_{\epsilon}\,(\nu_{1},\nu_{2})=g(\mu_{1},\nu_{2})^{-1}$.
Thus, $\lambda_{1}$ is absolutely continuous with respect to Lebesgue measure:

$\lambda_{1}(A)=\iint_{A\cup(-1/2,1/2]^{2}}\frac{\sigma_{\epsilon}^{2}}{g(\nu_{1},\nu_{2})}\,d\nu_{1}\,d\nu_{2}.$

The measure $\lambda_{2}$ is concentrated at not more than 2 points. Thus,
$\lambda_{1}$ and $\lambda_{2}$ are absolutely continuous and discrete
components of the measure $\lambda_{X}$, and non-discrete singular component
is zero.

Let $(\nu_{1}^{(0)},\nu_{2}^{(0)})$ be one of the points of $A_{2}$. As
$g_{1}(\nu_{1},\nu_{2})$ is differentiable and $g(\nu_{1},\nu_{2})\geq 0$, at
the neighbourhood of $(\nu_{1}^{(0)},\nu_{2}^{(0)})$

$\displaystyle
g_{1}(\nu_{1},\nu_{2})=o(|\nu_{1}-\nu_{1}^{(0)}|+|\nu_{2}-\nu_{2}^{(0)}|),$
$\displaystyle
g(\nu_{1},\nu_{2})=|g_{1}(\nu_{1},\nu_{2})|^{2}=o((\nu_{1}-\nu_{1}^{(0)})^{2}+(\nu_{2}-\nu_{2}^{(0)})^{2}),$
$\displaystyle\frac{1}{g(\nu_{1},\nu_{2})}>\frac{\textrm{const}}{(\nu_{1}-\nu_{1}^{(0)})^{2}+(\nu_{2}-\nu_{2}^{(0)})^{2}}$

for some $\textrm{const}>0$;

$\int_{-1/2}^{1/2}\int_{-1/2}^{1/2}\frac{1}{g(\nu_{1},\nu_{2})}\,d\nu_{1}\,d\nu_{2}=+\infty.$

Hence

$\lambda_{X}\left(\left(-{\textstyle\frac{1}{2}},{\textstyle\frac{1}{2}}\right]^{2}\right)\geq\lambda_{1}\left(\left(-{\textstyle\frac{1}{2}},{\textstyle\frac{1}{2}}\right]^{2}\right)=\iint_{(-1/2,1/2]^{2}}\frac{\sigma^{2}}{g(\nu_{1},\nu_{2})}\,d\nu_{1}\,d\nu_{2}=\infty.$

Thus, the spectral measure of the random field $X$ is infinite, which is
impossible. The assumption $D\leq 0$ brings the contradiction. Thus, $D>0$.

##### Sufficiency.

Let $D<0$. Then $\sigma_{\epsilon}^{2}/g(\nu_{1},\nu_{2})$ is an even
integrable function on $\left[-\frac{1}{2},\frac{1}{2}\right]^{2}$ that
attains only positive values. Then there exists a zero-mean Gaussian
stationary random field $X$ with spectral density
$\sigma_{\epsilon}^{2}/g(\nu_{1},\nu_{2})$, see Lemma 6. Define $\epsilon$ by
formula (1). Then $\epsilon$ is zero-mean and stationary; $\epsilon$ has a
constant spectral density $\sigma_{\epsilon}^{2}$. Thus, $\epsilon$ is a
collection of uncorrelated random variables with zero mean and variance
$\sigma_{\epsilon}^{2}$. Thus, the random fields $\epsilon$ and $X$ satisfy
the desired conditions. ∎

######  Remark 1.

In the sufficiency part of Proposition 1, the probability space is specially
constructed. The random fields $X$ and $\epsilon$ constructed are jointly
Gaussian and jointly stationary, which means that
$\\{(X_{i,j},\epsilon_{i,j})\;i,j\mathbin{\in}\mathbb{Z}\\}$ is a two-
dimensional Gaussian stationary random field on a two-dimensional lattice.

######  Corollary 1.

Let $X$ be a wide-sense stationary field that satisfies (1). Then $X$ is
centered and has a spectral density (7).

###### Proof.

Due to Proposition 1, $D>0$ (where $D$ comes from Lemma 1); Hence,
$1-a-b-c\neq 0$. Taking expectation in (1) immediately imply that the mean of
the stationary random field $X$ is zero.

Demonstration that $X$ has a spectral density partially repeat the Necessity
part of the proof of Proposition 1; in what follows, we use notation
$g(\nu_{1},\nu_{2})$, $\lambda_{X}$, $A_{2}$, $\lambda_{1}$ and $\lambda_{2}$
from that proof. According to Lemma 1, $D>0$ imply that $A_{2}=\emptyset$. The
discrete component $\lambda_{2}$ of the spectral measure $\lambda_{X}$ of the
random filed $X$ is concentrated on the set $A_{2}$; thus, $\lambda_{2}=0$.
The non-discrete singular component of $\lambda_{X}$ is also zero. The
spectral measure $\lambda_{X}=\lambda_{1}+\lambda_{2}$ is absolutely
continuous, that is the random field $X$ has a spectral density. ∎

###  2.3 Restatement of the existence theorem

######  Proposition 2.

Let $\sigma^{2}>0$, $a$, $b$ and $c$ be fixed real numbers, and
$\epsilon_{i,j}$ be a collection of uncorrelated random variables with zero
mean and equal variance $\sigma^{2}$. A wide-sense stationary random field $X$
that satisfy equation (1) exists if and only if $D>0$, where $D$ comes from
Lemma 1.

###### Proof.

Necessity is logically equivalent to one in Proposition 1.

Sufficiency. Denote

$\psi_{k,l}=\int_{-1/2}^{1/2}\int_{-1/2}^{1/2}\frac{\exp(-2\pi\mathrm{i}(k\nu_{1}+l\nu_{2}))}{1-ae^{2\pi\mathrm{i}\nu_{1}}-be^{2\pi\mathrm{i}\nu_{2}}-ce^{2\pi\mathrm{i}(\nu_{1}+\nu_{2})}}\,d\nu_{1}\,d\nu_{2}\,.$
(10)

The integrand is well-defined due to Lemma 1. It is infinitely times
differentiable periodic function. Due to Lemma 7,

$\sum_{k=-\infty}^{\infty}\sum_{l=-\infty}^{\infty}|\psi_{k,l}|<\infty.$ (11)

Now prove that

$X_{i,j}=\sum_{k=-\infty}^{\infty}\sum_{l=-\infty}^{\infty}\psi_{k,l}\epsilon_{i-k,j-l}$
(12)

is a solution. The series in the is convergent in mean squares as well as
almost surely due to Proposition 15. The resulting field $X$ is stationary as
a linear transformation of a stationary field $\epsilon$.

By changing the indices,

$\displaystyle X_{i,j}$ $\displaystyle-aX_{i-1,j}-bX_{i,j-1}-cX_{i-1,j-1}=$
$\displaystyle=\sum_{k=-\infty}^{\infty}\sum_{l=-\infty}^{\infty}\psi_{k,l}(\epsilon_{i-k,j-l}-a\epsilon_{i-1-k,j-l}-b\epsilon_{i-k,j-1-l}-c\epsilon_{i-1-k,j-1-l})=$
$\displaystyle=\sum_{k=-\infty}^{\infty}\sum_{l=-\infty}^{\infty}(\psi_{k,l}-a\psi_{k-1,l}-b\psi_{k,l-1}-c\psi_{k-1,l-1})\epsilon_{i-l,j-1}.$

This transformation of the Fourier coefficients corresponds to multiplication
of the integrand by a certain function. Thus,

$\displaystyle\begin{aligned}
\psi_{k,l}-a\psi_{k-1,l}&-b\psi_{k,l-1}-c\psi_{k-1,l-1}=\\\
&=\int_{-1/2}^{1/2}\int_{-1/2}^{1/2}\exp(-2\pi\mathrm{i}(k\nu_{1}+l\nu_{2}))\,d\nu_{1}\,d\nu_{2}=\begin{cases}1\mbox{~{}if~{}$k=l=0$},\\\
0\mbox{~{}otherwise};\end{cases}\end{aligned}$ $\displaystyle
X_{i,j}-aX_{i-1,j}-bX_{i,j-1}-cX_{i-1,j-1}=\epsilon_{i,j},$

and the random field $X$ is indeed a solution to (1). Thus, $X$ is a desired
random field. ∎

##  3 The autocovariance function

Assume that $\\{X_{i,j},i,j\in\mathbb{Z}\\}$ is a wide-sense stationary field
that satisfies (1). According to Corollary 1, $X$ has a spectral density,
which is defined by (7). The autocovariance function can be evaluated with (4)
or (5). The explicit formula is

$\displaystyle\gamma_{X}(h_{1},h_{2})$
$\displaystyle=\int_{-1/2}^{1/2}\int_{-1/2}^{1/2}\frac{\sigma_{\epsilon}^{2}\exp(2\pi\mathrm{i}(h_{1}\nu_{1}+h_{2}\nu_{2}))}{|1-ae^{2\pi\mathrm{i}\nu_{1}}-be^{2\pi\mathrm{i}\nu_{2}}-ce^{2\pi\mathrm{i}(\nu_{1}+\nu_{2})}|^{2}}\,d\nu_{1}\,d\nu_{2}=$
$\displaystyle=\int_{-1/2}^{1/2}\int_{-1/2}^{1/2}\frac{\sigma_{\epsilon}^{2}\exp(2\pi\mathrm{i}(h_{1}\nu_{1}+h_{2}\nu_{2}))}{g(\nu_{1},\nu_{2})}\,d\nu_{1}\,d\nu_{2},$
(13)

with $g(\nu_{1},\nu_{2})$ defined in Lemma 1.

For fixed $h_{1}$ and $h_{2}$, it is possible to compute the integral in (13)
as a function of $a$, $b$ and $c$.

###  3.1 The autocovariance function for $h_{1}=0$.

The denominator $g(\nu_{1},\nu_{2})$ in (7) is equal to

$\displaystyle g(\nu_{1},\nu_{2})$
$\displaystyle=1+a^{2}+b^{2}+c^{2}+2(ac-b)\cos(2\pi\nu_{2})+\mbox{}$
$\displaystyle\quad+2((bc-a)+(ab-c)\cos(2\pi\nu_{2}))\cos(2\pi\nu_{1})+\mbox{}$
$\displaystyle\quad+2(ab+c)\sin(2\pi\nu_{1})\sin(2\pi\nu_{2})=$
$\displaystyle=A+B\cos(2\pi\nu_{1})+C\sin(2\pi\nu_{1}),$

with

$\begin{gathered}A=1+a^{2}+b^{2}+c^{2}+2(ac-b)\cos(2\pi\nu_{2}),\\\
B=2((bc-a)+(ab-c)\cos(2\pi\nu_{2})),\qquad
C=2(ab+c)\sin(2\pi\nu_{2})\end{gathered}$ (14)

We are going to use Lemma 8 to calculate the integral
$\int_{-1/2}^{1/2}g(\nu_{1},\nu_{2})^{-1}\,d\nu_{1}$. We have to verify the
conditions $A>0$ and $A^{2}-B^{2}-C^{2}>0$:

$(1+a^{2}+b^{2}+c^{2})^{2}-4(ac-b)^{2}=((a-c)^{2}+(b+1)^{2})((a+c)^{2}+(b-1)^{2})\geq
0.$ (15)

The zero is attained if either $a=c$ and $b=-1$, or $a+c=0$ and $b=1$. In both
cases $D=0$, while $D>0$ according to Proposition 1 ($D$ is defined in Lemma
1). Thus the inequality in (15) is strict,

$\displaystyle(1+a^{2}+b^{2}+c^{2})^{2}>4(ac-b)^{2},$ $\displaystyle
1+a^{2}+b^{2}+c^{2}>|2(ac-b)|\geq-2(ac-b)\cos(2\pi\nu_{2}),$ $\displaystyle
A=1+a^{2}+b^{2}+c^{2}+2(ac-b)\cos(2\pi\nu_{2})>0.$

With some trigonometric transformations,

$A^{2}-B^{2}-C^{2}=(1-a^{2}+b^{2}-c^{2}-2(ac+b)\cos(2\pi\nu_{2}))^{2}.$

Again, with some transformations and notion that $D>0$,

$\displaystyle(1-a^{2}+b^{2}-c^{2})^{2}-4(ac+b)^{2}=D>0,$
$\displaystyle|1-a^{2}+b^{2}-c^{2}|>|2(ac+b)|>2(ac+b)\cos(2\pi\nu_{2}),$
$\displaystyle 1-a^{2}+b^{2}-c^{2}\neq 2(ac+b)\cos(2\pi\nu_{2}),$
$\displaystyle
A^{2}-B^{2}-C^{2}=(1-a^{2}+b^{2}-c^{2}-2(ac+b)\cos(2\pi\nu_{2}))^{2}>0.$

According to Lemma (8), eq. (42),

$\int_{-1/2}^{1/2}\frac{d\nu_{1}}{g(\nu_{1},\nu_{2})}=\frac{1}{\sqrt{A^{2}-B^{2}-C^{2}}}=\frac{1}{|1-a^{2}+b^{2}-c^{2}-2(ac+b)\cos(2\pi\nu_{2})|}.$

Denote

$g_{2}(\nu_{2})=1-a^{2}+b^{2}-c^{2}-2(ac+b)\cos(2\pi\nu_{2})$

and compute
$\int_{-1/2}^{1/2}e^{2\pi\mathrm{i}\nu_{2}h_{2}}\,|g_{2}(\nu_{2})|^{-1}\,d\nu_{2}$
using formula (42). The expression $g_{2}(\nu_{2})$ does not attain $0$;
hence, either $g_{2}(\nu_{2})>0$ for all $\nu_{2}$, or $g_{2}(\nu_{2})<0$ for
all $\nu_{2}$. Then

$|g_{2}(\nu_{2})|=A_{2}+B_{2}\cos(2\pi\nu_{2}),$

where $A_{2}=1-a^{2}+b^{2}-c^{2}$ and $B_{2}=-2(ac+b)$ if $g_{2}(\nu_{2})>0$
for all $\nu_{2}$; otherwise, $A_{2}=-(1-a^{2}+b^{2}-c^{2})$ and
$B_{2}=2(ac+b)$ if $g_{2}(\nu_{2})<0$ for all $\nu_{2}$. In both cases,
$A_{2}=|g_{2}(\pi/2)|>0$ and

$A_{2}^{2}-B_{2}^{2}=(1-a^{2}+b^{2}-c^{2})^{2}-4(ac+b)^{2}=D>0.$

According to Lemma 8, eq. (43),

$\int_{-1/2}^{1/2}\frac{\exp(2\pi\mathrm{i}h_{2}\nu_{2})}{|g_{2}(\nu_{2})|}\,d\nu_{2}=\frac{\beta^{|h_{2}|}}{\sqrt{A_{2}^{2}-B_{2}^{2}}}=\frac{1}{\sqrt{D}}$

where

$\beta=\frac{B_{2}}{A_{2}+\sqrt{A_{2}^{2}-B_{2}^{2}}}=\frac{B_{2}}{A_{2}+\sqrt{D}}.$

The value of $\beta$ in

$\displaystyle\beta$ $\displaystyle=0\quad\mbox{if $ac+b=0$},$
$\displaystyle\beta$
$\displaystyle=\frac{1-a^{2}+b^{2}-c^{2}-\sqrt{D}}{2(ac+b)}\quad\mbox{if
$ac+b\neq 0$ and $g_{2}(\nu_{2})>0$ for all $\nu_{2}$},$ $\displaystyle\beta$
$\displaystyle=\frac{1-a^{2}+b^{2}-c^{2}+\sqrt{D}}{2(ac+b)}\quad\mbox{if
$ac+b\neq 0$ and $g_{2}(\nu_{2})<0$ for all $\nu_{2}$};$

the formula that if valid in all cases is

$\beta=\frac{2(ac+b)}{1-a^{2}+b^{2}-c^{2}+\operatorname{sign}(1-a^{2}+b^{2}-c^{2})\sqrt{D}}.$
(16)

Finally,

$\displaystyle\gamma_{X}(0,h_{2})$
$\displaystyle=\int_{-1/2}^{1/2}\int_{-1/2}^{1/2}\frac{\sigma_{\epsilon}^{2}\exp(2\pi\mathrm{i}h_{2}\nu_{2})}{g(\nu_{1},\nu_{2})}\,d\nu_{1}\,d\nu_{2}=$
$\displaystyle=\int_{-1/2}^{1/2}\frac{\sigma_{\epsilon}^{2}\exp(2\pi\mathrm{i}h_{2}\nu_{2})}{|g_{2}(\nu_{2})|}\,d\nu_{2}=\frac{\beta^{|h_{2}|}\sigma_{\epsilon}^{2}}{\sqrt{D}}.$
(17)

In particular,

$\operatorname{var}X_{i,j}=\gamma_{X}(0,0)=\frac{\sigma_{\epsilon}^{2}}{\sqrt{D}}.$
(18)

Due to symmetry,

$\gamma_{X}(h_{1},0)=\frac{\alpha^{|h_{1}|}\sigma_{\epsilon}^{2}}{\sqrt{D}},$
(19)

where

$\alpha=\frac{2(a+bc)}{1+a^{2}-b^{2}-c^{2}+\operatorname{sign}(1+a^{2}-b^{2}-c^{2})\sqrt{D}}.$
(20)

Notice that $|\alpha|<1$ and $|\beta|<1$, whence

$|\gamma_{X}(1,0)|<\gamma_{X}(0,0)\quad\mbox{and}\quad|\gamma_{X}(0,1)|<\gamma_{X}(0,0).$
(21)

###  3.2 Yule–Walker equations

We formally write down Yule–Walker equations for the autocovariance function
of $X$. These equations are particular cases of ones given in [10, Section 6]:

$\displaystyle\gamma_{X}(0,\,0)$
$\displaystyle=a\gamma_{X}(1,\,0)+b\gamma_{X}(0,\,1)+c\gamma_{X}(1,\,1)+\sigma_{\epsilon}^{2},$
(22) $\displaystyle\gamma_{X}(h_{1},\,h_{2})$
$\displaystyle=a\gamma_{X}(h_{1}-1,\,h_{2})+b\gamma_{X}(h_{1},\,h_{2}-1)+c\gamma_{X}(h_{1}-1,\,h_{2}-1).$
(23)

In Lemma 2 we obtain conditions for (23) to hold true. As (22) requires
specific conditions on coefficients $a$, $b$ and $c$, we postpone the
consideration of (22) to Section 4.

######  Lemma 2.

Let $X$ be a stationary field satisfying equation (1), and let
$\gamma_{X}(h_{1},h_{2})$ be the covariance function of the process $X$. Then
equality (23) holds true if in any of the following conditions hold true:

  1. (1)

$h_{1}\geq 1$, $(1-a-b-c)(1+a-b+c)>0$ and $(1-a+b+c)(1+a+b-c)>0$,

  2. (2)

$h_{1}\leq 0$, $(1-a-b-c)(1+a-b+c)<0$ and $(1-a+b+c)(1+a+b-c)<0$,

  3. (3)

$h_{2}\geq 1$, $(1-a-b-c)(1-a+b+c)>0$ and $(1+a-b+c)(1+a+b-c)>0$, or

  4. (4)

$h_{2}\leq 0$, $(1-a-b-c)(1-a+b+c)<0$ and $(1+a-b+c)(1+a+b-c)<0$.

###### Proof.

According to Corollary 1, the stationary field $X$ has a spectral density,
that is (4) with $f_{X}(h_{1},h_{2})$ defined in (7). Hence, as to rules how
the Fourier transform changes when the function is linearly transformed,

$\gamma_{X}(h_{1},\,h_{2})-a\gamma_{X}(h_{1}-1,\,h_{2})-b\gamma(h_{1},\,h_{2}-1)-c\gamma(h_{1}-1,\,h_{2}-1)=\\\
\begin{aligned}
&=\int_{-1/2}^{1/2}\int_{-1/2}^{1/2}e^{2\pi\mathrm{i}(h_{1}\nu_{1}+h_{2}\nu_{2})}(1-ae^{-2\pi\mathrm{i}\nu_{1}}-be^{-2\pi\mathrm{i}\nu_{2}}-ce^{-2\pi\mathrm{i}(\nu_{1}+\nu_{2})})\times\mbox{}\\\
&\qquad\qquad\qquad\times f_{X}(h_{1},h_{2})\,d\nu_{1}\,d\nu_{2}=\\\
&=\int_{-1/2}^{1/2}\int_{-1/2}^{1/2}\frac{\exp(2\pi{\mathrm{i}}(h_{1}\nu_{1}+h_{2}\nu_{2}))\sigma_{\epsilon}^{2}}{1-ae^{2\pi\mathrm{i}\nu_{1}}-be^{2\pi\mathrm{i}\nu_{2}}-ce^{2\pi\mathrm{i}(\nu_{1}+\nu_{2})}}\,d\nu_{1}\,d\nu_{2}.\end{aligned}$
(24)

##### Case (1):

$h_{1}\geq 1$, $(1-a-b-c)(1+a-b+c)>0$ and $(1-a+b+c)(1+a+b-c)>0$. In this
case,

$\displaystyle(1-b)^{2}>(a+c)^{2}\quad\mbox{and}\quad(1+b)^{2}>(a-c)^{2},$
$\displaystyle 1-a^{2}+b^{2}-c^{2}>2b+2ac\quad\mbox{and}\quad
1-a^{2}+b^{2}-c^{2}>-2b-2ac,$ $\displaystyle 1-a^{2}+b^{2}-c^{2}>2\,|b+ac|.$

Then, for every $\nu_{2}\in[-1/2,\>1/2]$, since $|\cos(2\pi\nu_{2})|\leq 1$,

$\displaystyle 1-a^{2}+b^{2}-c^{2}>2(b+ac)\cos(2\pi\nu_{2}),$ $\displaystyle
1-2b\cos(2\pi\nu_{2})+b^{2}>a^{2}+2ac\cos(2\pi\nu_{2})+c^{2},$
$\displaystyle|1-be^{2\pi\mathrm{i}\nu_{2}}|^{2}>|a+ce^{2\pi\mathrm{i}\nu_{2}}|^{2},$
$\displaystyle|1-be^{2\pi\mathrm{i}\nu_{2}}|>|a+ce^{2\pi\mathrm{i}\nu_{2}}|.$
(25)

According to Lemma 9,

$\int_{-1/2}^{1/2}\frac{\exp(2\pi{\mathrm{i}}h_{1}\nu_{1})}{1-ae^{2\pi\mathrm{i}\nu_{1}}-be^{2\pi\mathrm{i}\nu_{2}}-ce^{2\pi\mathrm{i}(\nu_{1}+\nu_{2})}}\,d\nu_{1}=0,$
(26)

which, with (24), implies (23).

##### Case (2):

$h_{1}\leq 0$, $(1-a-b-c)(1+a-b+c)<0$ and $(1-a+b+c)(1+a+b-c)<0$. In this
case,

$\displaystyle(1-b)^{2}<(a+c)^{2}\quad\mbox{and}\quad(1+b)^{2}<(a-c)^{2},$
$\displaystyle-2b-2ac<a^{2}-b^{2}+c^{2}-1\quad\mbox{and}\quad
2b+2ac<a^{2}-b^{2}+c^{2}-1,$ $\displaystyle 2\,|b+ac|<a^{2}-b^{2}+c^{2}-1.$

Then, for every $\nu_{2}\in[-1/2,\>1/2]$, since $|\cos(2\pi\nu_{2})|\leq 1$,

$\displaystyle-2(b+ac)\cos(2\pi\nu_{2})\leq 2\,|b+ac|<a^{2}-b^{2}+c^{2}-1,$
$\displaystyle 1-2b\cos(2\pi\nu_{2})+b^{2}<a^{2}+2ac\cos(2\pi\nu_{2})+c^{2},$
$\displaystyle|1-be^{2\pi\mathrm{i}\nu_{2}}|^{2}<|a+ce^{2\pi\mathrm{i}\nu_{2}}|^{2},$
$\displaystyle|1-be^{2\pi\mathrm{i}\nu_{2}}|<|a+ce^{2\pi\mathrm{i}\nu_{2}}|.$

Again, according to another case of Lemma 9, (26) still holds true. With (24),
it implies (23).

##### Cases (3) and (4)

are symmetric to cases (1) and (2), respectively. ∎

###  3.3 Uniqueness of the solution

We prove that under Sufficient conditions of Proposition 2, the stationary
solution to (1) is unique.

######  Proposition 3.

Let $\sigma^{2}>0$, $a$, $b$ and $c$ be fixed real numbers, and
$\epsilon_{i,j}$ be a collection of uncorrelated random variables with zero
mean and equal variance $\sigma^{2}$. If $D>0$, then equation (1) has only one
stationary solution $X$, namely one defined by (12) with coefficients
$\psi_{k,l}$ defined in (10). The coefficients $\psi_{k,l}$ satisfy (11);
hence, the double series in (12) converges in least squares and almost surely.

###### Proof.

It remains to prove the uniqueness of the solution, as other assertions follow
from the proof of Proposition 2.

Let $X$ be a solution to (1), which can be rewritten as (3) Let $\mathcal{H}$
be the Hilbert space spanned by $L_{x}^{k}L_{y}^{l}$,
$k,l=\ldots,-1,0,1,\ldots$ with scalar product

$\langle
Y^{(1)},Y^{(2)}\rangle=\operatorname{\mathsf{E}}Y^{(1)}_{0,0}Y^{(2)}_{0,0}.$

Formally, $\mathcal{H}$ can be defined as the closure of the set of random
fields of form
$\\{\sum_{k=-n}^{n}\sum_{l=-n}^{n}\phi_{k,l}X_{i-k,j-l},\;\allowbreak
i,j=\ldots,-1,0,1,\ldots\\}$ in a Banach space of bounded-second-moment random
fields with norm

$\|Y\|=\sup_{i,j}(\operatorname{\mathsf{E}}Y_{i,j}^{2})^{1/2}.$

All random fields within $\mathcal{H}$ are jointly stationary, which imply
that the scalar product indeed corresponds to the norm of the Banach space:

$\langle Y,Y\rangle=\|Y\|^{2}\quad\mbox{for all}\quad Y\in\mathcal{H}.$

Since (3), $\epsilon\in\mathcal{H}$. Let us construct a similar space for
$\epsilon$. Let

$\displaystyle\mathcal{H}_{1}$
$\displaystyle=\left\\{\sum_{k=-\infty}^{\infty}\sum_{l=-\infty}^{\infty}\theta_{k,l}L_{x}^{k}L_{y}^{k}\epsilon:\sum_{k=-\infty}^{\infty}\sum_{l=-\infty}^{\infty}\theta_{k,l}^{2}<\infty\right\\}=$
$\displaystyle=\left\\{\left\\{\sum_{k=-\infty}^{\infty}\sum_{l=-\infty}^{\infty}\theta_{k,l}\epsilon_{i-k,j-l},\,i,j=\ldots,-1,0,1,\ldots\\!\right\\}:\sum_{k=-\infty}^{\infty}\sum_{l=-\infty}^{\infty}\theta_{k,l}^{2}<\infty\\!\right\\}$

The series is convergent, see Proposition 16.

The random fields $\sigma_{\epsilon}^{-1}L_{x}^{k}L_{y}^{k}\epsilon$ make an
orthonormal basis of the subspace $\mathcal{H}_{1}$. The orthogonal projector
onto $\mathcal{H}_{1}$ can be defined as follows

$\displaystyle
P_{1}Y=\frac{1}{\sigma_{\epsilon}^{2}}\sum_{k=-\infty}^{\infty}\sum_{l=-\infty}^{\infty}(\operatorname{\mathsf{E}}Y_{k,l}\epsilon_{0,0})L^{k}L^{l}\epsilon,$
$\displaystyle(P_{1}Y)_{i,j}=\frac{1}{\sigma_{\epsilon}^{2}}\sum_{k=-\infty}^{\infty}\sum_{l=-\infty}^{\infty}\operatorname{\mathsf{E}}(Y_{k,l}\epsilon_{0,0})\epsilon_{i-k,j-l}$

for all $Y\in\mathcal{H}$.

It is possible to verify that

$P_{1}L_{x}=L_{x}P_{1},\qquad P_{1}L_{y}=L_{y}P_{1},\qquad
P_{1}\epsilon=\epsilon.$

Applying the operator $P_{1}$ to both sides of (3), we get that the random
field $P_{1}X$ is also a stationary solution to (1). Both random fields must
be centered, and variances of both the fields satisfy (18):

$\|X\|^{2}=\|P_{1}X\|=\frac{\sigma_{\epsilon}^{2}}{\sqrt{D}}.$

This implies that $\|X-P_{1}X\|=0$, which means $X=P_{1}X$ almost surely.

Combining (11) and (24), we get

$\displaystyle\gamma_{X}(-k,-l)-a\gamma_{X}(-k{-}1,-l)-b\gamma_{X}(-k,-l{-}1)-c\gamma_{X}(-k{-}1,-l{-}1)=\sigma_{\epsilon}^{2}\psi_{k,l},$
$\displaystyle\operatorname{\mathsf{E}}X_{0,0}X_{k,l}-a\operatorname{\mathsf{E}}X_{-1,0}X_{k,l}-b\operatorname{\mathsf{E}}X_{0,-1}X_{k,l}-c\operatorname{\mathsf{E}}X_{-1,-1}X_{k,l}=\sigma_{\epsilon}^{2}\psi_{k,l},$
$\displaystyle\operatorname{\mathsf{E}}(X_{0,0}-aX_{-1,0}-bX_{0,-1}-cX_{-1,-1})X_{k,l}=\sigma_{\epsilon}^{2}\psi_{k,l},$
$\displaystyle\operatorname{\mathsf{E}}\epsilon_{0,0}X_{k,l}=\sigma_{\epsilon}^{2}\psi_{k,l},$
(27) $\displaystyle
X_{i,j}=PX_{i,j}=\frac{1}{\sigma_{\epsilon}^{2}}\sum_{k=-\infty}^{\infty}\sum_{l=-\infty}^{\infty}(\operatorname{\mathsf{E}}X_{k,l}\epsilon_{0,0})\epsilon_{i-k,j-l}=\sum_{k=-\infty}^{\infty}\sum_{l=-\infty}^{\infty}\psi_{k,l}\epsilon_{i-k,j-l}.$

Thus, the solution to (1) must satisfy (12) almost surely, whence the
uniqueness follows. ∎

##  4 Stability and causality

###  4.1 Causality

We generalize the causality condition for random fields on a two-dimensional
lattice.

######  Definition 1.

Let $\epsilon_{i,j}$ ($i$ and $j$ integers) be random variables of zero mean
and constant, finite, non-zero variance. A random field
$\\{X_{i,j},\;i,j=\ldots,-1,0,1,\ldots\\}$ is said to be causal with respect
to (w.r.t.) the white noise $\epsilon$ if there exists an array of
coefficients $\\{\psi_{k,l},\;k,l=0,1,2,\ldots\\}$ such that

$\sum_{k=0}^{\infty}\sum_{l=0}^{\infty}|\psi_{k,l}|<\infty$ (28)

and the process $X$ allows the representation

$X_{i,j}=\sum_{k=0}^{\infty}\sum_{l=0}^{\infty}\psi_{k,l}\epsilon_{i-k,j-l}.$
(29)

A causal random field must be wide-sense stationary.

######  Proposition 4.

Let $\epsilon_{i,j}$ be a uncorrelated zero-mean equal-variance variables,
$\operatorname{\mathsf{E}}\epsilon_{i,j}=0$ and
$\operatorname{\mathsf{E}}\epsilon_{i,j}^{2}=\sigma_{\epsilon}^{2}$,
$0<\sigma_{\epsilon}^{2}<\infty$ Let $X$ be a stationary field that satisfies
(1). The random field $X$ is causal w.r.t. the white noise $\epsilon$ if and
only if these four inequalities hold true:

$\displaystyle 1-a-b-c$ $\displaystyle>0,$ $\displaystyle 1-a+b+c$
$\displaystyle>0,$ $\displaystyle 1+a-b+c$ $\displaystyle>0,$ $\displaystyle
1+a+b-c$ $\displaystyle>0.$

If the field $X$ is causal, then the coefficients in representation (29) equal

$\psi_{k,l}=\sum_{m=0}^{\min(k,l)}\binom{k}{m}\binom{l}{m}a^{k-m}b^{l-m}(ab+c)^{m}.$
(30)

###### Proof.

Sufficiency. According to Proposition 3, the field $X$ admits a representation
(12) with coefficients $\psi_{k,l}$ satisfying (28). Now prove that
$\psi_{k,l}=0$ if $k<0$ or $l<0$.

Let $k<0$. Let us look into the proof of Lemma 2, case (1), for $h_{1}=-k$.
Eq. (26) holds true,

$\int_{-1/2}^{1/2}\frac{\exp(-2\pi{\mathrm{i}}k\nu_{1})}{1-ae^{2\pi\mathrm{i}\nu_{1}}-be^{2\pi\mathrm{i}\nu_{2}}-ce^{2\pi\mathrm{i}(\nu_{1}+\nu_{2})}}\,d\nu_{1}=0,$

whence $\psi_{k,l}$ defined in (10) is zero,

$\psi_{k,l}=\int_{-1/2}^{1/2}e^{-2\pi\mathrm{i}l\nu_{2}}\biggl{(}\int_{-1/2}^{1/2}\frac{\exp(-2\pi\mathrm{i}k\nu_{1})}{1-ae^{2\pi\mathrm{i}\nu_{1}}-be^{2\pi\mathrm{i}\nu_{2}}-ce^{2\pi\mathrm{i}(\nu_{1}+\nu_{2})}}\,d\nu_{1}\biggr{)}\,d\nu_{2}=0\,.$

The case $l<0$ is symmetric to $k<0$; if $l<0$, then $\psi_{k,l}=0$ also holds
true. Thus, (12) rewrites as (29).

Necessity. Denote

$f_{1}\mathbin{=}1-a-b-c,\quad f_{2}\mathbin{=}1-a+b+c,\quad
f_{3}\mathbin{=}1+a-b+c,\quad f_{4}\mathbin{=}1+a+b-c.$ (31)

Obviously, $f_{1}+f_{2}+f_{3}+f_{4}=4>0$, and due to Proposition 1,
$D=f_{1}f_{2}f_{3}f_{4}>0$.

Now prove that $f_{1}f_{3}>0$. Equality $f_{1}f_{3}=0$ is impossible. Let
$f_{1}f_{3}<0$. Then $f_{2}f_{4}<0$. The representation (12) is unique: the
coefficients $\psi_{k,l}$ must be equal to
$\sigma_{\epsilon}^{-2}\operatorname{cov}(X_{0,0},\epsilon_{-k,-l})$. Thus,
coefficients in (29) must be defined by (10). Prove that $\psi_{k,l}=0$ for
all integers $k>0$ and $l>0$. Look into the proof of Lemma 2, case (2), for
$h_{1}=-k$. Eq. (26) still holds true and implies $\psi_{k,l}=0$ the same way
as in the proof of sufficiency. The representation (29) rewrites as
$X_{i,j}=0$. The assumption $f_{1}f_{3}<0$ have led to a contradiction. Thus,
$f_{1}f_{3}>0$.

The inequality $f_{1}f_{2}>0$ holds true due to symmetry. Inequalities
$f_{1}f_{2}>0$, $f_{1}f_{3}>0$, $f_{1}f_{2}f_{3}f_{4}>0$ and
$f_{1}+f_{2}+f_{3}+f_{4}>0$ imply that $f>0$, $f_{2}>0$, $f_{3}>0$ and
$f_{4}>0$.

The expression for $\psi_{k,l}$. Again, if $f_{1}>0$, $f_{2}>0$, $f_{3}>0$ and
$f_{4}>0$, in the proof of Lemma 2, case (1), inequality (25) was
demonstrated. Let $k$ and $l$ be nonnegative integers. By (10) and Lemmas 9
and 10,

$\displaystyle\psi_{k,l}$
$\displaystyle=\int_{-1/2}^{1/2}e^{-2\pi\mathrm{i}l\nu_{2}}\biggl{(}\int_{-1/2}^{1/2}\frac{\exp(-2\pi\mathrm{i}k\nu_{1})}{1-ae^{2\pi\mathrm{i}\nu_{1}}-be^{2\pi\mathrm{i}\nu_{2}}-ce^{2\pi\mathrm{i}(\nu_{1}+\nu_{2})}}\,d\nu_{1}\biggr{)}\,d\nu_{2}=$
$\displaystyle=\int_{-1/2}^{1/2}\frac{e^{-2\pi\mathrm{i}l\nu_{2}}(a+ce^{2\pi\mathrm{i}\nu_{2}})^{k}}{(1-be^{2\pi\mathrm{i}\nu_{2}})^{k+1}}\,d\nu_{2}=$
$\displaystyle=\sum_{n=0}^{\min(k,l)}\binom{k}{n}\binom{l}{n}a^{k-n}b^{l-n}(ab+c)^{n}.\qed$

###  4.2 Autocovariance function and Yule–Walker equations under the causality
condition

######  Proposition 5.

Let coefficients $a$, $b$ and $c$ be such that $f_{1}>0$, $f_{2}>0$, $f_{3}>0$
and $f_{4}>0$ for $f_{1},\ldots,f_{4}$ defined in (31). Let $X$ be a
stationary field and $\epsilon$ be a collection of zero-mean equal-variance
random variables, $\operatorname{var}\epsilon_{i,j}=\sigma_{\epsilon}^{2}$,
that satisfy (1). Then $\operatorname{cov}(X_{i,j},\epsilon_{k,l})=0$ if $k>i$
or $l>j$, and the following Yule–Walker equations hold true:

$\displaystyle\gamma_{X}(0,0)$
$\displaystyle=a\gamma_{X}(1,0)+b\gamma_{X}(0,1)+c\gamma_{X}(1,1)+\sigma_{\epsilon}^{2},$
$\displaystyle\gamma_{X}(h_{1},h_{2})$
$\displaystyle=a\gamma_{X}(h_{1}{-}1,h_{2})+b\gamma_{X}(h_{1},h_{2}{-}1)+c\gamma_{X}(h_{1}{-}1,h_{2}{-}1)$

if $\max(h_{1},h_{2})>0$.

###### Proof.

Because of causality and representation (29),

$\operatorname{cov}(X_{i,j},\epsilon_{k,l})=\sum_{s=0}^{\infty}\sum_{r=0}^{\infty}\psi_{r,s}\operatorname{cov}(\epsilon_{i-r,j-s}\epsilon_{k,l}).$

If $k>i$, then $i-r<k$ in all summands; if $l>j$, then $j-s<l$. In ether case,
$\operatorname{cov}(\epsilon_{i-r,j-s}\epsilon_{k,l})=0$ is all terms. Hence,
$\operatorname{cov}(X_{i,j},\epsilon_{k,l})=0$.

Eq. (23) for $h_{1}>0$ or $h_{2}>0$ follows from Lemma 2, case (1) and
case(3), respectively.

Because of (1),

$\displaystyle\operatorname{cov}(X_{i,j},\epsilon_{i,j})=\operatorname{cov}(aX_{i-1,j}+bX_{i,j-1}+cX_{i-1,j-1}+\epsilon_{i,j},\>\epsilon_{i,j}),$
$\displaystyle\begin{aligned}
\operatorname{cov}(X_{i,j},\epsilon_{i,j})&=a\operatorname{cov}(X_{i-1,j},\epsilon_{i,j})+b\operatorname{cov}(X_{i,j-1},\epsilon_{i,j})+\mbox{}\\\
&\quad+c\operatorname{cov}(X_{i-1,j-1},\epsilon_{i-1,j-1})+\operatorname{cov}(\epsilon_{i,j},\epsilon_{i,j}),\end{aligned}$
$\displaystyle\operatorname{cov}(X_{i,j},\epsilon_{i,j})=\sigma_{\epsilon}^{2}.$

Then

$\displaystyle\operatorname{cov}(X_{i,j},X_{i,j})=\operatorname{cov}(X_{i,j},\>aX_{i-1,j}+bX_{i,j-1}+cX_{i-1,j-1}+\epsilon_{i,j}),$
$\displaystyle\begin{aligned}
\operatorname{cov}(X_{i,j},X_{i,j})&=a\operatorname{cov}(X_{i,j},X_{i-1,j})+b\operatorname{cov}(X_{i,j},X_{i,j-1})+\mbox{}\\\
&\quad+c\operatorname{cov}(X_{i,j}X_{i-1,j-1})+\operatorname{cov}(X_{i,j},\epsilon_{i,j}),\end{aligned}$
$\displaystyle\gamma_{X}(0,0)=a\gamma_{X}(1,0)+b\gamma_{X}(0,1)+c\gamma_{X}(1,1)+\sigma_{\epsilon}^{2}.\qed$

######  Proposition 6.

Let coefficients $a$, $b$ and $c$ be such that $f_{1}>0,\>\ldots,\>f_{4}>0$,
with $f_{1},\ldots,f_{4}$ defined in (31). Let $X$ be a stationary random
field, and $\epsilon$ be a collection of zero-mean random variables with equal
variance $\operatorname{\mathsf{E}}\epsilon_{i,j}^{2}$ that satisfy (1). Then
random fields $X$ and $\epsilon$ are jointly stationary with cross-
autocovariance function

$\displaystyle\operatorname{cov}(X_{i+h_{1},j+h_{2}},\epsilon_{i,j})=0\qquad\mbox{if
$h_{1}<0$ or $h_{2}<0$},$
$\displaystyle\operatorname{cov}(X_{i+h_{1},j+h_{2}},\epsilon_{i,j})=\sum_{k=0}^{\min(h_{1},h_{2})}\binom{h_{1}}{k}\binom{h_{2}}{k}a^{h_{1}-k}b^{h_{2}-k}(ab+c)^{k}\sigma_{\epsilon}^{2}$

if $h_{1}\geq 0$ and $h_{2}\geq 0$.

###### Proof.

The joint stationarity follows from (3) and the stationarity of $X$ (this was
already used in the proof of Proposition 3). The formula for the
autocovariance function follows from (27) and Proposition 4. ∎

######  Proposition 7.

Under conditions of Proposition 6, the autocovariance function of the field
$X$ satisfies

$\gamma_{X}(h_{1},h_{2})=\frac{|\alpha|^{|h_{1}|}|\beta|^{|h_{2}|}\sigma_{\epsilon}^{2}}{\sqrt{D}}\qquad\mbox{if
$h_{1}h_{2}\leq 0$,}$ (32)

where $D$ is defined in Lemma 1, and

$\alpha=\frac{2(a+bc)}{1+a^{2}-b^{2}-c^{2}+\sqrt{D}},\qquad\beta=\frac{2(ac+b)}{1-a^{2}+b^{2}-c^{2}+\sqrt{D}}.$

###### Proof.

First, notice that

$\displaystyle
1+a^{2}-b^{2}-c^{2}=\frac{(1{-}a{-}b{-}c)(1{-}a{+}b{+}c)+(1{+}a{-}b{+}c)(1{+}a{+}b{-}c)}{2}>0,$
$\displaystyle
1-a^{2}+b^{2}-c^{2}=\frac{(1{-}a{-}b{-}c)(1{+}a{-}b{+}c)+(1{-}a{+}b{+}c)(1{+}a{+}b{-}c)}{2}>0.$

Thus, $\alpha$ and $\beta$ are the same as in (16) and (20), and equality in
(32) holds true for $h_{1}h_{2}=0$.

Consider three cases.

##### Case 1:

$b\neq 0$. Eq. (32) hold true for $h_{1}=0$ and all $h_{2}$. Assume it holds
for $h_{1}=n_{1}$ and all $h_{2}$, for $n_{1}$ a nonnegative integer. Eq. (32)
hold true for $h_{2}=0$ (and all $h_{1}$). Assume it holds true for
$h_{2}=-n_{2}$ and $h_{1}=n_{1}+1$, for $n_{2}$ a nonnegative integer. Now
prove that (32) holds true for $h_{2}=-n_{2}-1$ and all $h_{1}=n_{1}+1$. Apply
Lemma 2, conditions (1), for $h_{1}=n_{1}+1\geq 1$ and $h_{2}=-n_{2}$:

$\gamma_{X}(n_{1}{+}1,\,{-}n_{2})=a\gamma_{X}(n_{1},\,-n_{2})+b\gamma_{X}(n_{1}{+}1,\,-n_{2}{-}1)+c\gamma_{X}(n_{1},\,-n_{2}-1).$

Hence

$\displaystyle\gamma_{X}(n_{1}{+}1,{-}n_{2}{-}1)$
$\displaystyle=\frac{\gamma_{X}(n_{1}{+}1,\,-n_{2})-a\gamma_{X}(n_{1},\,-n_{2})-c\gamma_{X}(n_{1},\,-n_{2}-1)}{b}=$
$\displaystyle=\frac{1}{b}\left(\frac{\alpha^{n_{1}+1}\beta^{n_{2}}\sigma_{\epsilon}^{2}}{\sqrt{D}}-\frac{a\alpha^{n_{1}}\beta^{n_{2}}\sigma_{\epsilon}^{2}}{\sqrt{D}}-\frac{c\alpha^{n_{1}}\beta^{n_{2}+1}\sigma_{\epsilon}^{2}}{\sqrt{D}}\right)=$
$\displaystyle=\frac{\alpha^{n_{1}}\beta^{n_{2}}\sigma_{\epsilon}^{2}\,(\alpha-
a-c\beta)}{b\,\sqrt{D}}=\frac{\alpha^{n_{1}+1}\beta^{n_{2}+1}\sigma_{\epsilon}^{2}}{\sqrt{D}};$

equality $\alpha-a-c\beta=b\alpha\beta$ can be verified by easy but tedious
computations. Thus, (32) holds true for $h_{2}=-n_{2}-1$ and all
$h_{1}=n_{1}+1$. By induction, it holds true for $h_{1}=n_{1}+1$ and all
$h_{2}\leq 0$, and then for all $h_{1}\geq 0$ and $h_{2}\leq 2$.

Since both sides of (32) are even functions in $(h_{1},h_{2})$, that is
$\gamma_{X}(h_{1},h_{2})=-\gamma_{X}(-h_{1},-h_{2})$, (32) also holds true for
$h_{1}\leq 0$ and $h_{2}\geq 0$. Thus, (32) holds true for all integer $h_{1}$
and $h_{2}$ such that $h_{1}h_{2}\leq 0$.

##### Case 2:

$b=0$ but $a\neq 0$. This case is symmetrical to case $a=0$ and $b\neq 0$,
which is covered by case 1.

##### Case 3:

$a=b=0$. Then $\alpha=\beta=0$ and

$\gamma_{X}(h_{1},h_{2})=\int_{-1/2}^{1/2}\int_{-1/2}^{1/2}\frac{\exp(2\pi\mathrm{i}(h_{1}\nu_{1}+h_{2}\nu_{2}))}{|1-c\exp(2\pi\mathrm{i}(\nu_{1}+\nu_{2}))|^{2}}\,d\nu_{1}\,d\nu_{2}.$

Substituting $\nu_{1}=\nu_{3}-\nu_{2}$ and taking in account the periodicity,
we get

$\displaystyle\gamma_{X}(h_{1},h_{2})$
$\displaystyle=\int_{\nu_{2}=-1/2}^{1/2}\int_{\nu_{3}=\nu_{2}-\frac{1}{2}}^{\nu_{2}+\frac{1}{2}}\frac{\exp(2\pi\mathrm{i}(h_{1}\nu_{3}-h_{1}\nu_{2}+h_{2}\nu_{2}))}{|1-c\exp(2\pi\mathrm{i}(\nu_{3}))|^{2}}\,d\nu_{3}\,d\nu_{2}=$
$\displaystyle=\int_{\nu_{2}=-1/2}^{1/2}\int_{\nu_{3}=-\frac{1}{2}}^{\frac{1}{2}}\frac{\exp(2\pi\mathrm{i}(h_{1}\nu_{3}-h_{1}\nu_{2}+h_{2}\nu_{2}))}{|1-c\exp(2\pi\mathrm{i}\nu_{3})|^{2}}\,d\nu_{3}\,d\nu_{2}=$
$\displaystyle=\int_{-1/2}^{1/2}\exp(2\pi\mathrm{i}(h_{2}-h_{1})\nu_{2})\,d\nu_{2}\,\int_{\nu_{3}=-\frac{1}{2}}^{\frac{1}{2}}\frac{\exp(2\pi\mathrm{i}h_{1}\nu_{3})}{1-2c\cos(2\pi\nu_{3})+c^{2}}\,d\nu_{3}.$

If $h_{1}h_{2}<0$, then $h_{1}\neq h_{2}$ and
$\int_{-1/2}^{1/2}\exp(2\pi\mathrm{i}(h_{2}-h_{1})\nu_{2})\,d\nu_{2}=0$,
whence $\gamma_{X}(h_{1},h_{2})=0$, and (32) holds true. As noted before (32)
holds true if $h_{1}h_{2}=0$. Thus, (32) holds true for all integer $h_{1}$
and $h_{2}$ such that $h_{1}h_{2}\leq 0$. ∎ Under conditions of Proposition 7,

$\gamma_{X}(1,-1)\gamma_{X}(0,0)=\gamma_{X}(1,0)\gamma_{X}(0,1).$ (33)

###  4.3 Stability

######  Lemma 3.

Let $a$, $b$ and $c$ be such that $f_{1},\ldots,f_{4}$ defined in (31) are
positive real numbers. Let $\psi_{k,l}$ be defined in (30) if $k\geq 0$ and
$l\geq 0$, and $\psi_{k,l}=0$ if $k<0$ or $l<0$. Then $\psi_{k,l}$ satisfy
this equation:

$\psi_{k,l}-a\psi_{k-1,l}-b\psi_{k,l-1}-c\psi_{k-1,l-1}=\begin{cases}1\quad\mbox{if
$k=l=0$},\\\ 0\quad\mbox{otherwise}.\end{cases}$ (34)

.

###### Proof.

Under conditions $f_{1}>0,\,\ldots,\,f_{4}>0$, the coefficients $\psi_{k,l}$
satisfy (10) for all integer $k$ and $l$. Respective transforms of the Fourier
coefficients and the integrand in (10) yield

$\psi_{k,l}-a\psi_{k-1,l}-b\psi_{k,l-1}-c\psi_{k-1,l-1}=\\\
=\int_{-1/2}^{1/2}\int_{-1/2}^{1/2}\exp(-2\pi\mathrm{i}(k\nu_{1}+l\nu_{2}))\,d\nu_{1}\,d\nu_{2}=\begin{cases}1\quad\mbox{if
$k=l=0$},\\\ 0\quad\mbox{otherwise}.\end{cases}\qed$

######  Remark 2.

In Lemma 3 the condition on $a$, $b$ and $c$ is unnecessary. Indeed, the
coefficients $\psi_{k,l}$ are polynomials in $a$, $b$ and $c$; thus, the left-
hand side of (34) is also a polynomial in $a$, $b$ and $c$ and the fact that
(34) holds true for small $a$, $b$ and $c$ implies that (34) holds true for
all real $a$, $b$ and $c$.

######  Proposition 8.

Let $a$, $b$ and $c$ be such that $f_{1},\ldots,f_{4}$ defined in (31) are
positive real numbers. Consider the recurrence equation

$x_{i,j}=ax_{i-1,j}+bx_{i,j-1}+cx_{i-1,j-1}+v_{i,j},$ (35)

for all positive integers $i$ and $j$, where

  * •

$v_{i,j}$ are known variables, $i>0$ and $j>0$;

  * •

$x_{0,0}$, $x_{i,0}$ and $x_{0,j}$ are preset, so setting their values makes
initial/boundary conditions;

  * •

$x_{i,j}$ ($i>0$ and $j>0$) are unknown variables.

Then the explicit formula for the solution is the following:

$\displaystyle x_{i,j}$
$\displaystyle=\psi_{i,j}x_{0,0}+\sum_{k=0}^{i-1}\psi_{k,j}\,(x_{i-k,0}-ax_{i-k-1,0})+\mbox{}$
$\displaystyle\quad+\sum_{l=0}^{j-1}\psi_{i,l}\,(x_{0,j-l}-bx_{0,j-l-1})+\sum_{k=0}^{i-1}\sum_{l=0}^{j-1}\psi_{k,l}v_{i-k,j-l},$
(36)

with $\psi_{k,l}$ defined in (30).

If, in addition, the sequences of $x_{i,0}$ and $x_{0,j}$ and the array of
$v_{i,j}$ are bounded, then the solution $x$ is also bounded.

###### Proof.

First, show that under convention that “an empty sum is zero”, namely
$\sum_{i=0}^{-1}\ldots=0$ and $\sum_{j=0}^{-1}\ldots=0$, the right side of
(36) satisfies the bounding conditions. For one or both indices being zero,
the coefficients $\psi_{k,l}$ equal

$\psi_{0,0}=1,\qquad\psi_{k,0}=a^{k},\qquad\psi_{0,l}=b^{l}$

for all integer $k\geq 0$ and $l\geq 0$. Hence,

$\displaystyle\psi_{0,0}x_{0,0}$ $\displaystyle=x_{i,j},$
$\displaystyle\psi_{0,j}x_{0,0}$
$\displaystyle+\sum_{l=0}^{j-1}\psi_{0,l}(x_{0,j-l}-bx_{0,j-l-1})=b^{l}x_{0,0}+\sum_{l=0}^{j-1}b^{l}(x_{0,j-k}-bx_{0,j-l-1})=$
$\displaystyle=b^{l}x_{0,0}+\sum_{l=0}^{j-1}b^{l}x_{0,j-l}-\sum_{l=0}^{j-1}b^{l+1}x_{0,j-l-1}=$
$\displaystyle=b^{l}x_{0,0}+\sum_{l=0}^{j-1}b^{l}x_{0,j-l}-\sum_{l=1}^{j}b^{l}x_{0,j-l}=b^{l}x_{0,0}+x_{0,j}-b^{l}x_{0,0}=x_{0,j}$

for integer $j>0$, and due to symmetry for all integers $i>0$

$\psi_{i,0}x_{0,0}+\sum_{k=0}^{i-1}\psi_{k,0}\,(x_{i-k,0}-ax_{i-k-1,0})=x_{i,0}.$

Known variables $v_{i,j}$ are defined for all integer $i>0$ and $j>0$. Extend
their domain for all integer $i\geq 0$ and $j\geq 0$ by assigning

$v_{0,0}=0;\qquad v_{i,0}=x_{i,0}-ax_{i-1,0},\quad i>0;\qquad
v_{0,j}=x_{0,j}-bx_{0,j-1},\quad j>0.$

Then the right-hand side of (36) is

$x^{\rm RHS}_{i,j}=\sum_{k=0}^{i}\sum_{l=0}^{j}\psi_{k,l}v_{i-k,j-l},\qquad
i\geq 0,\quad j\geq 0.$

In particular, $x^{\rm RHS}_{i,j}=x_{i,j}$ if either $i=0$ or $j=0$.

Now prove that $x^{\rm RHS}$ satisfies (35). To that and, let $\psi_{k,l}=0$
if $k<0$ or $l<0$. Then, for $i$ and $j$ positive integers,

$x^{\rm RHS}_{i,j}-ax^{\rm RHS}_{i-1,j}-bx^{\rm RHS}_{i,j-1}-cx^{\rm
RHS}_{i-1,j-1}=\\\ \begin{aligned}
&=\sum_{k=0}^{i}\sum_{l=0}^{j}\psi_{k,l}v_{i-k,j-l}-a\sum_{k=0}^{i-1}\sum_{l=0}^{j}\psi_{k,l}v_{i-1-k,j-l}-\mbox{}\\\
&\quad-b\sum_{k=0}^{i}\sum_{l=0}^{j-1}\psi_{k,l}v_{i-k,j-1-l}-c\sum_{k=0}^{i-1}\sum_{l=0}^{j-1}\psi_{k,l}v_{i-1-k,j-1-l}=\\\
&=\sum_{k=0}^{i}\sum_{l=0}^{j}\psi_{k,l}v_{i-k,j-l}-a\sum_{k=1}^{i}\sum_{l=0}^{j}\psi_{k-1,l}v_{i-k,j-l}-\mbox{}\\\
&\quad-b\sum_{k=0}^{i}\sum_{l=1}^{j}\psi_{k,l-1}v_{i-k,j-l}-c\sum_{k=1}^{i}\sum_{l=1}^{j}\psi_{k-1,l-1}v_{i-k,j-l}=\\\
&=\sum_{k=0}^{i}\sum_{l=0}^{j}(\psi_{k,l}-a\psi_{k-1,l}-b\psi_{k,l-1}-c\psi_{k-1,l-1})v_{i-k,j-l}=v_{i,j}\,;\end{aligned}$

here we used Lemma 3.

Finally, equality $x_{i,j}=x^{\rm RHS}_{i,j}$ for all integer $i\geq 0$ and
$j\geq 0$ can be proved by induction.

If initial/boundary values and $v_{i,j}$ are bounded,

$\sup_{i>0}|x_{i,0}|<\infty,\qquad\sup_{j>0}|x_{0,j}|<\infty,\quad\mbox{and}\quad\sup_{i,j>00}|v_{i,j}|<\infty,$

then

$\displaystyle|x_{i,j}|$
$\displaystyle=\left|\sum_{k=0}^{i}\sum_{k=0}^{j}\psi_{k,l}v_{i-k,j-l}\right|\leq\sum_{k=0}^{i}\sum_{j=0}^{l}|\psi_{i,j}|\max_{0\leq
r\leq 1}\max{0\leq s\leq 1}|v_{r,s}|\leq$
$\displaystyle\leq\sum_{k=0}^{i}\sum_{j=0}^{l}|\psi_{i,j}|\times\begin{aligned}
\max\Bigl{(}&|x_{0,0}|,\;(1+|a|)\max_{1\leq r\leq i}|x_{r,0}|,\\\
&(1+|b|)\max_{1\leq s\leq i}|x_{0,s}|,\;\max_{1\leq r\leq i}\max_{1\leq s\leq
i}|v_{i,j}|\Bigr{)}\leq\end{aligned}$
$\displaystyle\leq\sum_{k=0}^{\infty}\sum_{j=0}^{\infty}|\psi_{i,j}|\times\begin{aligned}
\max\Bigl{(}&|x_{0,0}|,\;(1+|a|)\sup_{r\geq 1}|x_{r,0}|,\\\
&(1+|b|)\sup_{s\geq 1}|x_{0,s}|,\;\sup_{r,s\geq
1}|v_{i,j}|\Bigr{)}<\infty.\end{aligned}$

Thus, the solution $x$ is bounded. Here we used (28). ∎

##  5 Symmetry, nonidentifiability, and special cases

In this section we consider the question how parameters change when the random
field $X$ is flipped, and whether two or more sets of parameters correspond to
the same distribution (or the same autocovariance function) of $X$. Here
parameters are the coefficients $a$, $b$ and $c$ and the variance of the error
term $\sigma_{\epsilon}^{2}$. Usually, two sets sets of parameters correspond
to the same autocovariance function. In special cases, up to four sets sets of
parameters correspond to the same autocovariance function. We start with
special cases.

###  5.1 Special cases

####  5.1.1 Symmetric case: $ab+c=0$.

######  Proposition 9.

Let a stationary field $X$ and a collection of uncorrelated zero-mean equal-
variance random variables $\epsilon$ be a solution to (1). The autocovariance
function $\gamma_{X}$ is even with respect to each variable,
$\gamma_{X}(h_{1},h_{2})=\gamma_{X}(-h_{1},h_{2})=\gamma_{X}(h_{1},-h_{2})$ if
and only if $ab+c=0$.

###### Proof.

Eq. (4) can be rewritten as

$\gamma_{X}(h_{1},\,h_{2})=\int_{-1/2}^{1/2}e^{2\pi{\mathrm{i}}h_{2}\nu_{2}}\int_{-1/2}^{1/2}e^{2\pi{\mathrm{i}}h_{1}\nu_{1}}f_{X}(\nu_{1},\nu_{2})\,d\nu_{1}\,d\nu_{2}.$

The even symmetry with respect to $h_{2}$,
$\gamma_{X}(h_{1},h_{2})=\gamma_{X}(h_{1},-h_{2})$, is equivalent to the inner
integral being a real function,

$\displaystyle\int_{-1/2}^{1/2}e^{2\pi{\mathrm{i}}h_{1}\nu_{1}}f_{X}(\nu_{1},\nu_{2})\,d\nu_{1}=$
$\displaystyle=\int_{-1/2}^{1/2}\frac{\exp(2\pi\,{\mathrm{i}}\,h_{1}\nu_{1})\sigma_{\epsilon}^{2}}{|1-ae^{2\pi{\mathrm{i}}\nu_{1}}-be^{2\pi{\mathrm{i}}\nu_{2}}-ce^{2\pi{\mathrm{i}}(\nu_{1}+\nu_{2})}|^{2}}\,d\nu_{1}=$
$\displaystyle=\int_{-1/2}^{1/2}\frac{\exp(2\pi\,{\mathrm{i}}\,h_{1}\nu_{1})\sigma_{\epsilon}^{2}}{A+B\cos(2\pi\nu_{1})+C\sin(2\pi\nu_{2})}\in\mathbb{R}\mbox{\quad
for almost all $\nu_{2}\in\left[-\frac{1}{2},\frac{1}{2}\right]$},$

with $A$, $B$ and $C$ defined in (14) in Section 3.1. In turn, this is
equivalent to the function Fourier-transformed being an even function,

$\frac{\sigma_{\epsilon}^{2}}{A+B\cos(2\pi\nu_{1})+C\sin(2\pi\nu_{1})}=\frac{\sigma_{\epsilon}^{2}}{A+B\cos(-2\pi\nu_{1})+C\sin(-2\pi\nu_{1})}$

for almost all $\nu_{1}$ and
$\nu_{2}\in\left[-\frac{1}{2},\frac{1}{2}\right]$, which holds true if and
only if $C=2(ab+c)\sin(2\pi\nu_{2})=0$ for all $\nu_{2}$, or $ab+c=0$. ∎

Thus, we consider the case $c=-ab$. The following proposition establishes
conditions on $a$ and $b$.

######  Proposition 10.

Let $c=-ab$. The stationary field $X$ and a collection of uncorrelated zero-
mean variables with equal variance $\sigma_{\epsilon}^{2}$ which satisfy (1)
exist if and only if $|a|\neq 1$ and $|b|\neq 1$. In that case, the
autocovariance function of $X$ is equal

$\gamma_{X}(h_{1},h_{2})=\frac{a^{\pm h_{1}}b^{\pm
h_{2}}\sigma_{\epsilon}^{2}}{|1-a^{2}|\times|1-b^{2}|},$ (37)

where the signs “$\pm$” are chosen so that $|a^{\pm h_{2}}|<1$ and $|b^{\pm
h_{2}}|<1$.

The field $X$ is causal with respect to $\epsilon$ is and only if $|a|<1$ and
$|b|<1$.

###### Proof.

Apply Proposition 1. Since $c=-ab$,

$\displaystyle D$ $\displaystyle=(1-a-b+ab)(1-a+b-ab)(1+a-b-ab)(1+a+b+ab)=$
$\displaystyle=(1-a)^{2}(1-b)^{2}(1+b)^{2}(1+a)^{2}.$

Thus, $D\geq 0$, and the necessary and sufficient condition $D>0$ for the
desired process to exist is satisfied if and only if $(1-a)(1-b)(1+b)(1+a)\neq
0$, which is equivalent to $|a|\neq 1$ and $|b|\neq 1$.

The spectral density (7) splits into two factors

$f_{X}(\nu_{1},\nu_{2})=\frac{1}{|1-ae^{2\pi\mathrm{i}\nu_{1}}|^{2}}\times\frac{\sigma_{\epsilon}^{2}}{|1-be^{2\pi\mathrm{i}\nu_{2}}|^{2}}.$

Thus, the autocovariance function is

$\gamma_{X}(h_{1},\,h_{2})=\int_{-1/2}^{1/2}\frac{\exp(2\pi
h_{1}\nu_{1})}{|1-ae^{2\pi\mathrm{i}\nu_{1}}|^{2}}\,d\nu_{1}\times\int_{-1/2}^{1/2}\frac{\exp(2\pi
h_{2}\nu_{2})}{|1-be^{2\pi\mathrm{i}\nu_{2}}|^{2}}\,d\nu_{2}\,\sigma_{\epsilon}^{2}.$

With Lemma 8,

$\displaystyle\int_{-1/2}^{1/2}\frac{\exp(2\pi
h_{1}\nu_{1})}{|1-ae^{2\pi\mathrm{i}\nu_{1}}|^{2}}\,d\nu_{1}$
$\displaystyle=\int_{-1/2}^{1/2}\frac{\exp(2\pi
h_{1}\nu_{1})}{1+a^{2}-2a\cos(2\pi\nu_{1})}\,d\nu_{1}=$
$\displaystyle=\frac{\alpha^{|h_{1}|}}{\sqrt{(1+a^{2})^{2}-4a^{2}}}=\frac{\alpha^{|h_{1}|}}{|1-a^{2}|}$

with

$\alpha=\frac{2a}{1+a^{2}+|1-a^{2}|}=\begin{cases}a&\mbox{if $|a|<1$},\\\
a^{-1}&\mbox{if $|a|>1$}.\end{cases}$

Thus,

$\int_{-1/2}^{1/2}\frac{\exp(2\pi
h_{1}\nu_{1})}{|1-ae^{2\pi\mathrm{i}\nu_{1}}|^{2}}\,d\nu_{1}=\frac{a^{\pm
h_{1}}}{|1-a^{2}|}.$

Similarly,

$\int_{-1/2}^{1/2}\frac{\exp(2\pi
h_{2}\nu_{2})}{|1-be^{2\pi\mathrm{i}\nu_{2}}|^{2}}\,d\nu_{2}=\frac{b^{\pm
h_{2}}}{|1-b^{2}|},$

and (37) holds true.

The causality conditions from Proposition 4 rewrite as follows:

$\displaystyle(1-a)(1-b)$ $\displaystyle>0,$ $\displaystyle(1-a)(1+b)$
$\displaystyle>0,$ $\displaystyle(1+a)(1-b)$ $\displaystyle>0,$
$\displaystyle(1+a)(1+b)$ $\displaystyle>0.$

They mean that $1-a$, $1-b$, $1+a$ and $1+b$ should be (nonzero and) of the
same sign. As sum of these factors $(1-a)+(1-b)+(1+a)+(1+b)=4$ is positive,
the causality condition is equivalent to

$1-a>0,\qquad 1-b>0,\qquad 1+a>0,\quad\mbox{and}\quad 1+b>0,$

which in turn is equivalent to $|a|<1$ and $|b|<1$. ∎

####  5.1.2 Case $a=-bc$: uncorrelated observations along transect line

Another special case is where the coefficients in (1) satisfy relation
$a=-bc$.

######  Proposition 11.

Let $a=-bc$. The stationary field $X$ and a collection of uncorrelated zero-
mean variables with equal variance $\sigma_{\epsilon}^{2}$ which satisfy (1)
exist if and only if $|b|\neq 1$ and $|c|\neq 1$. The autocovariance function
of $X$ satisfies relations

$\gamma_{X}(h_{1},0)=0,\quad h_{1}\neq
0\,;\qquad\gamma_{X}(0,h_{2})=\frac{b^{\pm
h_{2}}\sigma_{\epsilon}^{2}}{|1-b^{2}|\times|1-c^{2}|}$

where the sign “$\pm$” is chosen so that $|b^{\pm h_{2}}|<1$.

The field $X$ is causal with respect to $\epsilon$ is and only if $|b|<1$ and
$|c|<1$.

###### Proof.

Since $a=-ac$,

$\displaystyle D$ $\displaystyle=(1+bc-b-c)(1+bc+b+c)(1-bc-b+c)(1-bc+b-c)=$
$\displaystyle=(1-b)^{2}(1-c)^{2}(1+b)^{2}(1+c)^{2}.$

Thus, $D\geq 0$, and the necessary and sufficient condition of existence of a
stationary solution, namely $D>0$, is satisfied if and only if
$(1-b)(1-c)(1+b)(1+c)\neq 0$, which is equivalent to $|b|\neq 1$ and $|c|\neq
1$.

For the autocovariance function, we use results of Section 3.1. With notation
(16) and (20),

$\displaystyle\sqrt{D}=|1-b^{2}|\times|1-c^{2}|,\qquad\alpha=0,$
$\displaystyle\begin{aligned}
\beta&=\frac{2(-b^{2}c+b)}{1-b^{2}c^{2}+b^{2}-c^{2}+\operatorname{sign}(1-b^{2}c^{2}+b^{2}-c^{2})\sqrt{D}}=\\\
&=\frac{2(1-c^{2})b}{(1+b^{2})(1-c^{2})+\operatorname{sign}(1-c^{2})|1-b^{2}|\times|1-c^{2}|}=\\\
&=\frac{2(1-c^{2})b}{(1+b^{2})(1-c^{2})+|1-b^{2}|\,(1-c^{2})}=\\\
&=\frac{2b}{(1+b^{2})+|1-b^{2}|}=\begin{cases}b&\mbox{if $|b|<1$,}\\\
b^{-1}&\mbox{if $|b|>1$.}\end{cases}\end{aligned}$

The desired formulas for the autocovariance function follow from (17) and
(19).

The causality conditions, which are obtained in Proposition 4, can be
rewritten as follows:

$\displaystyle(1-b)(1-c)$ $\displaystyle>0,$ $\displaystyle(1+b)(1+c)$
$\displaystyle>0,$ $\displaystyle(1-b)(1+c)$ $\displaystyle>0,$
$\displaystyle(1+b)(1-c)$ $\displaystyle>0.$

These conditions are equivalent to $|b|<1$ and $|c|<1$. The proof is similar
to one in Proposition 10. ∎

The values of the autocovariance function for $a=-0.1$, $b=0.5$, $c=0.2$, and
$\sigma_{\epsilon}^{2}=0.72$ are shown in Table 1. The parameters are chosen
in such a way that $\operatorname{var}X_{i,j}=1$.

Table 1: The autocovariance function $\gamma_{X}(h_{1},h_{2})$ for $a=-0.1$, $b=0.5$, $c=0.2$, and $\sigma_{\epsilon}^{2}=0.72$ | Values of $\gamma_{X}(h_{1},h_{2})$
---|---
$h_{2}$ | $h_{1}{=}{-}2$ | $h_{1}{=}{-}1$ | $h_{1}{=}0$ | $h_{1}{=}1$ | $h_{1}{=}2$ | $h_{1}{=}3$
3 | 0 | 0 | 0.125 | 0.1125 | 0.0225 | –0.00225
2 | 0 | 0 | 0.25 | 0.15 | 0.0075 | –0.003
1 | 0 | 0 | 0.5 | 0.15 | –0.015 | 0.0015
0 | 0 | 0 | 1 | 0 | 0 | 0
–1 | –0.015 | 0.15 | 0.5 | 0 | 0 | 0
–2 | 0.0075 | 0.15 | 0.25 | 0 | 0 | 0
–3 | 0.0225 | 0.1125 | 0.125 | 0 | 0 | 0

###  5.2 Nonidentifiability. Symmetry

Let $X$ be a stationary solution to eq. (1). Consider the relation between
parameters $a$, $b$, $c$ and $\sigma_{\epsilon}^{2}$ and the distribution of
the random field $X$ — more specifically, we study the multiplicity of this
relation. The situation is complicated by the fact that the distribution of
$X$ is not determined uniquely by the parameters $a$, $b$, $c$ and
$\sigma_{\epsilon}^{2}$; it also depends on the distribution of the error
field $\epsilon$ besides the scaling coefficient $\sigma_{\epsilon}$.

######  Definition 2.

Let us have a statistical model
$(\Omega,\mathcal{F},\;\mathsf{P}_{\theta,\phi},\>(\theta,\phi)\in\Psi)$, that
is a family of probability measures $\mathsf{P}_{\theta,\psi}$ on a common
measurable space $(\Omega,\mathcal{F})$ is indexed by parameter
$(\theta,\phi)\in\Psi$. The parameter $\theta$ is called identifiable if
$\mathsf{P}_{\theta_{1},\psi_{1}}=\mathsf{P}_{\theta_{2},\psi_{2}}$ implies
$\theta_{1}=\theta_{2}$.

In the model considered, the parameter of interest $\theta$ comprises the
coefficients and the error variance,
$\theta=(a,b,c,\sigma_{\epsilon}^{2})\in\Theta$, $\sigma_{\epsilon}>0$,

$\Theta=\\{(a,b,c,v):(1{-}a{-}b{-}c)(1{-}a{+}b{+}c)(1{+}a{-}b{+}c)(1{+}a{+}b{-}c)>0\mbox{~{}and~{}}v>0\\}.$

The nuisance “non-parametric” parameter $\phi$ describes the distribution of
the normalized error field $\sigma_{\epsilon}^{-1}\epsilon$. In the Gaussian
case, $\sigma_{\epsilon}^{-1}\epsilon$ is a collection of independent random
variables all of standard normal distribution. In general case,
$\sigma_{\epsilon}^{-1}\epsilon$ is a collection of uncorrelated zero-mean
unit-variance random variables.

It is much simpler to consider the autocovariance function
$\gamma_{X}(h_{1},h_{2})$ than the distribution of $X$, as
$\gamma_{X}(h_{1},h_{2})$ is uniquely determined by $\theta$ due to (13). In
the Gaussian case, the distribution of zero-mean stationary field $X$ is
uniquely determined by its autocovariance function; and obviously the
distribution uniquely determines the autocovariance function. In general case,
if parameters $\theta_{1}$ and $\theta_{2}$ determine the same autocovariance
function, than for any “distribution of a collection of uncorrelated zero-mean
unit-variance variables $\sigma_{\epsilon}^{-1}\epsilon$” described by
parameter $\phi_{1}$ there exists a corresponding distribution described by
parameter $\phi_{2}$ such that
$P_{\theta_{1},\phi_{1}}=P_{\theta_{2},\phi_{2}}$. Thus, it is enough to find
out which parameters $\theta$ yield the same autocovariance function
$\gamma_{X}(h_{1},h_{2})$.

In this section we will abuse notation by denoting the random field $X$ as
$X_{i,j}$. Similarly, $X_{-i,-j}$ will be a random field $X^{(2)}$ defined so
that $X^{(2)}_{i,j}=X_{-i,-j}$.

By substitution of indices, (1) can be rewritten as

$\displaystyle
X_{-i,-j}=-\frac{b}{c}X_{1-i,-j}-\frac{a}{c}X_{-i,1-j}+\frac{1}{c}X_{1-i,1-j}-\frac{1}{c}\epsilon_{1-i,1-j},$
$\displaystyle
X_{-i,j}=\frac{1}{a}X_{1-i,j}-\frac{c}{a}X_{-1,j-1}-\frac{b}{a}X_{1-i,j-1}-\frac{1}{a}\epsilon_{1-i,j},$
$\displaystyle
X_{i,-j}=-\frac{c}{b}X_{i-1,-j}+\frac{1}{b}X_{i,1-j}-\frac{a}{b}X_{i-1,1-j}-\frac{1}{b}\epsilon_{i,1-j}.$

Thus, random fields $X_{-i,-j}$, $X_{-i,j}$ and $X_{i,-j}$ satisfy equations
of the structure of (1), but with different coefficients and different error
term. The correspondence of parameters is presented in Table 2.

Table 2: Correspondence of terms of the recurrence equation for fields $X_{i,j}$, $X_{-i,-j}$, $X_{-i,j}$ and $X_{i,-j}$ Line | Field of | Coefficients | Error field
---|---|---|---
no., $m$ | interest |  |  |  | variance | values
1 | $X_{i,j}$ | $a$ | $b$ | $c$ | $\sigma_{\epsilon}^{2}$ | $\epsilon_{i,j}$
2 | $X_{-i,-j}$ | $-b/c$ | $-a/c$ | $1/c$ | $c^{-2}\sigma_{\epsilon}^{2}$ | $-c^{-1}\epsilon_{1-i,1-j}$
3 | $X_{-i,j}$ | $1/a$ | $-c/a$ | $-b/a$ | $a^{-2}\sigma_{\epsilon}^{2}$ | $-a^{-1}\epsilon_{1-i,j}$
4 | $X_{i,j-1}$ | $-c/b$ | $1/b$ | $-a/b$ | $b^{-2}\sigma_{\epsilon}^{2}$ | $-b^{-1}\epsilon_{i,1-j}$

However the random fields $X_{i,j}$ and $X_{-i,-j}$ have the same
autocovariance function,
$\gamma_{X_{i,j}}(h_{1},h_{2})=\gamma_{X_{-i,-j}}(h_{1},h_{2})$. Hence, two
sets of parameters $\theta_{1}=(a,b,c,\sigma_{\epsilon}^{2})$ and
$\theta_{2}=(-b/c,-a/c,1/c,c^{-2}\sigma_{\epsilon}^{2})$ determine the same
autocovariance function of the field $X$.

The autocovariance function of random fields $X_{-i,j}$ and $X_{i,-j}$ is
flip-symmetric to the autocovariance function of the field $X_{i,j}$:

$\gamma_{X_{-i,j}}(h_{1},h_{2})=\gamma_{X_{i,-j}}(h_{1},h_{2})=\gamma_{X_{i,j}}(-h_{1},h_{2})=\gamma_{X_{i,j}}(h_{1},-h_{2}).$

In the symmetric case $ab+c=0$, where the autocovariance function
$\gamma_{X_{i,j}}(h_{1},h_{2})$ is even in each of the arguments, the
autocovariance functions of the random fields $X_{i,j}$, $X_{-i,-j}$,
$X_{-i,j}$ and $X_{i,-j}$ are all equal. Hence, four sets of parameters
presented in Table 3 determine the same autocovariance function of the random
fields $X_{i,j}$. Hence, four sets of parameters
$\theta_{1}=(a,b,-bc,\allowbreak\sigma_{\epsilon}^{2})$,
$\theta_{2}=(a^{-1}!,b^{-1}\\!,-a^{-1}b^{-1}\\!,\allowbreak
a^{-2}b^{-2}\sigma_{\epsilon}^{2})$,
$\theta_{3}=(a^{-1}\\!,b,-a^{-1}b,\allowbreak a^{-2}\sigma_{\epsilon}^{2})$
and $\theta_{4}=(a,b^{-1}\\!,\allowbreak-ab^{-1}\\!,\allowbreak
b^{-2}\sigma_{\epsilon}^{2})$ determine the same autocovariance function of
the random field $X$.

Table 3: Four sets of parameters $\theta=(a,b,c,\sigma_{\epsilon}^{2})$ that determine the same autocovariance function $\gamma_{X}(h_{1},h_{2})$ $a$ | $b$ | $c=-ab$ | $\sigma_{\epsilon}^{2}$
---|---|---|---
$a^{-1}$ | $b^{-1}$ | $-a^{-1}b^{-1}$ | $a^{-2}b^{-2}\sigma_{\epsilon}^{2}$
$a^{-1}$ | $b$ | $-a^{-1}b$ | $a^{-2}\sigma_{\epsilon}^{2}$
$a$ | $b^{-1}$ | $-ab^{-1}$ | $b^{-2}\sigma_{\epsilon}^{2}$

Now we prove that there no other spurious cases where different sets of
parameters determine the same autocovariance function.

######  Lemma 4.

Let $D>0$, with $D$ defined in Lemma 1. Then of four sets of parameters listed
in Table 2, exactly one set satisfy the causality conditions specified in
Proposition 4.

The proof involves the search of 7 cases of different signs of
$f_{1},\ldots,f_{4}$ defined in (31), and evaluating inequalities. For
brevity, we do not present the full proof; however, in Table 4 we list for
what signs of $f_{1},\ldots,f_{4}$ which parameterization is causal.

Table 4: Relation between the signs of $f_{1},\ldots,f_{4}$ and line in Table 2 where parameters satisfy the causality condition Signs of $f_{1},\ldots,f_{4}$ | The | Which | Outcome
---|---|---|---
| signs | parameter- |
| imply | ization |
| that | is causal, $m$ |
$f_{1}>0$, $f_{2}>0$, $f_{3}>0$, $f_{4}>0$ |  | 1 |  $X_{i,j}$ is causal w.r.t. $\epsilon_{i,j}$
$f_{1}<0$, $f_{2}>0$, $f_{3}>0$, $f_{4}<0$ | $c>0$ | 2 |  $X_{-i,-j}$ is causal
$f_{1}>0$, $f_{2}<0$, $f_{3}<0$, $f_{4}>0$ | $c<0$ |  | w.r.t. $\epsilon_{1-i,1-j}$
$f_{1}<0$, $f_{2}<0$, $f_{3}>0$, $f_{4}>0$ | $a>0$ | 3 |  $X_{-i,j}$ is causal
$f_{1}>0$, $f_{2}>0$, $f_{3}<0$, $f_{4}<0$ | $a<0$ |  | w.r.t. $\epsilon_{1-i,j}$
$f_{1}<0$, $f_{2}>0$, $f_{3}<0$, $f_{4}>0$ | $b>0$ | 4 |  $X_{i,-j}$ is causal
$f_{1}>0$, $f_{2}<0$, $f_{3}>0$, $f_{4}<0$ | $b<0$ |  | w.r.t. $\epsilon_{i,1-j}$

######  Lemma 5.

In the causal case, that is under causality conditions stated in Proposition
4, the autocovariance function uniquely determines coefficients $a$, $b$ and
$c$ and error variance $\sigma_{\epsilon}^{2}$.

###### Proof.

Using Yule–Walker equations, (21) and (33), we can express the parameters $a$,
$b$, $c$ and $\sigma_{\epsilon}^{2}$ in terms of $\gamma_{X}(0,0)$,
$\gamma_{X}(1,0)$, $\gamma_{X}(0,1)$ and $\gamma_{X}(1,1)$. The explicit
formulas are

$\displaystyle
a=\frac{\gamma_{X}(1,0)\gamma_{X}(0,0)-\gamma_{X}(0,1)\gamma_{X}(1,1)}{\gamma_{X}(0,0)^{2}-\gamma_{X}(0,1)^{2}},$
$\displaystyle
b=\frac{\gamma_{X}(0,1)\gamma_{X}(0,0)-\gamma_{X}(1,0)\gamma_{X}(1,1)}{\gamma_{X}(0,0)^{2}-\gamma_{X}(1,0)^{2}},$
$\displaystyle
c=\frac{\gamma_{X}(1,1)-a\gamma_{X}(0,1)-b\gamma_{X}(1,0)}{\gamma_{X}(0,0)},$
$\displaystyle\sigma_{\epsilon}^{2}=\gamma_{X}(0,0)\sqrt{D},$

where $D$ is defined in Lemma 1. ∎

The following proposition provides a necessary condition for two sets of
parameters to determine the same autocovariance function of the field $X$.

######  Proposition 12.

Let stationary fields $X^{(1)}$ and $X^{(2)}$ satisfy (1)-like equation with
parameters $\theta_{1}=(a_{1},b_{1},c_{1},\sigma_{\epsilon,1}^{2})$ and
$\theta_{2}=(a_{2},b_{2},c_{2},\sigma_{\epsilon,2}^{2})$:

$\displaystyle
X^{(1)}_{i,j}=a_{1}X^{(1)}_{i-1,j}+b_{1}X^{(1)}_{i,j-1}+c_{1}X^{(1)}_{i-1,j-1}+\epsilon^{(1)}_{i,j},$
$\displaystyle
X^{(2)}_{i,j}=a_{2}X^{(2)}_{i-1,j}+b_{2}X^{(2)}_{i,j-1}+c_{2}X^{(2)}_{i-1,j-1}+\epsilon^{(2)}_{i,j},$
$\displaystyle\operatorname{var}\epsilon^{(1)}_{i,j}=\sigma_{\epsilon,1}^{2},\qquad\operatorname{var}\epsilon^{(2)}_{i,j}=\sigma_{\epsilon,2}^{2}.$

If random fields $X^{(1)}$ and $X^{(2)}$ have the same autocovariance
function, then the parameters $\theta_{1}$ and $\theta_{2}$ relate to each
other as quadruples of parameters in Table 2.

###### Proof.

Denote $T_{m}:\Theta_{m}\to\Theta$ the operator that transforms the quadruple
of parameters in the fist line of Table 2 into one in the $m$th line,
$m=1,2,3,4$. Thus,

$\displaystyle
T_{2}(a,b,c,\sigma_{\epsilon}^{2})=\textstyle\left(-\frac{b}{c},-\frac{a}{c},\frac{1}{c},c^{-2}\sigma_{\epsilon}^{2}\right),$
$\displaystyle
T_{3}(a,b,c,\sigma_{\epsilon}^{2})=\textstyle\left(\frac{1}{a},-\frac{c}{a},-\frac{b}{a},a^{-2}\sigma_{\epsilon}^{2}\right),$
$\displaystyle
T_{4}(a,b,c,\sigma_{\epsilon}^{2})=\textstyle\left(-\frac{c}{b},\frac{1}{b},-\frac{a}{b},b^{-2}\sigma_{\epsilon}^{2}\right),$

and $T_{1}$ is the identity operator, $T_{1}(\theta)=\theta$. Here
$\Theta_{m}\subset\Theta$ is the domain where the operator $T_{m}$ is well-
defined, e.g.,
$\Theta_{2}=\\{(a,b,c,\sigma_{\epsilon}^{2})\in\Theta:\allowbreak c\neq 0\\}$.
We have to prove that $\theta_{1}=T_{m}\theta_{0}$ and
$\theta_{2}=T_{n}\theta_{0}$ for some $m,n=1,\ldots,4$ and
$\theta_{0}\in\Theta$.

According to Lemma 4, one of parameterization
$T_{1}\theta_{1},\ldots,T_{4}\theta_{1}$ satisfies the conditions for
causality; let it be $T_{m}\theta_{1}=\theta_{3}$. Denote also $X^{(3)}$ and
$\epsilon^{(3)}$ the random field defined by such parameterization, and the
respective white noise; what this means is listed in Table 2. For example, if
$m=1$, that $X^{(3)}=X^{(1)}$, and if $m=2$, then
$X^{(3)}_{i,j}=X^{(1)}_{-i,-j}$.

The random fields $X^{(1)}$ and $X^{(3)}$ are either equal or flip or turn
symmetric to each other, $X^{(3)}_{i,j}=X^{(1)}_{\pm i,\pm j}$. Hence, their
autocovariance functions are either equal or one-variable symmetric to each
other: either

$\gamma_{X^{(3)}}(h_{1},h_{2})=\gamma_{X^{(1)}}(h_{1},h_{2})\quad\mbox{for all
$h_{1}$ and $h_{2}$},$

or

$\gamma_{X^{(3)}}(h_{1},h_{2})=\gamma_{X^{(1)}}(-h_{1},h_{2})=\gamma_{X^{(1)}}(h_{1},-h_{2})\quad\mbox{for
all $h_{1}$ and $h_{2}$}.$

The random field $X^{(3)}$ is causal w.r.t. white noise $\epsilon^{(3)}$.
Hence, due to Proposition 7,

$\gamma_{X^{(3)}}(h_{1},h_{2})=\frac{\gamma_{X^{(3)}}(h_{1},0)\gamma_{X^{(3)}}(0,h_{2})}{\gamma_{X^{(3)}}(0,0)}\quad\mbox{if
$h_{1}h_{2}<0$.}$ (38)

We do the same with the field $X^{(2)}$. There is $n$ and a field $X^{(4)}$
that satisfy a (1)-like equation with parameters $\theta_{4}=T_{n}\theta_{2}$;
these parameters satisfy the conditions for causality, and either

$\displaystyle\gamma_{X^{(4)}}(h_{1},h_{2})=\gamma_{X^{(2)}}(h_{1},h_{2})\quad\mbox{for
all $h_{1}$ and $h_{2}$, or}$
$\displaystyle\gamma_{X^{(4)}}(h_{1},h_{2})=\gamma_{X^{(2)}}(-h_{1},h_{2})=\gamma_{X^{(1)}}(h_{1},-h_{2})\quad\mbox{for
all $h_{1}$ and $h_{2}$},$

and also

$\gamma_{X^{(4)}}(h_{1},h_{2})=\frac{\gamma_{X^{(4)}}(h_{1},0)\gamma_{X^{(4)}}(0,h_{2})}{\gamma_{X^{(4)}}(0,0)}\quad\mbox{if
$h_{1}h_{2}<0$.}$ (39)

Similar relations hold for autocovariance functions of $X^{(3)}$ and
$X^{(4)}$: either

$\displaystyle\gamma_{X^{(4)}}(h_{1},h_{2})=\gamma_{X^{(3)}}(h_{1},h_{2})\quad\mbox{for
all $h_{1}$ and $h_{2}$, or}$
$\displaystyle\gamma_{X^{(4)}}(h_{1},h_{2})=\gamma_{X^{(3)}}(-h_{1},h_{2})=\gamma_{X^{(1)}}(h_{1},-h_{2})\quad\mbox{for
all $h_{1}$ and $h_{2}$},$

This implies $\gamma_{X^{(4)}}(h_{1},h_{2})=\gamma_{X^{(3)}}(h_{1},h_{2})$ if
$h_{1}=0$ or $h_{2}=0$, and because of (38) and (39),

$\gamma_{X^{(4)}}(h_{1},h_{2})=\gamma_{X^{(3)}}(h_{1},h_{2})\quad\mbox{if
$h_{1}h_{2}<0$.}$

All above implies that random fields $X^{(3)}$ and $X^{(4)}$ have the same
autocovariance function. These fields are causal w.r.t. respective white
noises. According to Lemma 5, their parameters are equal,
$\theta_{3}=\theta_{4}$.

The operators $T_{1},\ldots,T_{4}$ are self-inverse; $T_{k}^{2}$ is the
identity operator. Thus,
$T_{m}\theta_{1}=\theta_{3}=\theta_{4}=T_{n}\theta_{2}$ implies
$\theta_{1}=T_{n}\theta_{3}$ and $\theta_{2}=T_{m}\theta_{4}$. Thus,
parameters $\theta_{1}$ and $\theta_{2}$ relate to each other as the $m$th and
$n$th rows of Table 2. ∎

The following corollary combines necessary and sufficient conditions for
parameters to determine the same autocovariance function; it also takes into
account where the parameters make sense.

######  Corollary 2.

The parameter set $\Theta$ splits into two-element, four-element and one-
element classes of parameters that define the same autocovariance function
$\gamma_{X}(h_{1},h_{2})$ as follows:

  1. 1.

In generic case $-ab\neq c\neq 0$, the parameter
$(a,b,c,\sigma_{\epsilon}^{2})$ belongs to a two-element class. It determines
the same autocovariance function as the parameter
$(-b/c,-a/c,1/c,c^{-2}\sigma_{\epsilon}^{2})$.

  2. 2.

In symmetric case $-ab=c\neq 0$, the parameter
$(a,b,-bc,\sigma_{\epsilon}^{2})$ belongs to a four-element class. All four
sets of parameters listed in Table 3 determine the same autocovariance
function.

  3. 3.

In symmetric case $a\neq 0$, $b=c=0$, the parameter
$(a,0,0,\sigma_{\epsilon}^{2})$ belongs to a two-element class. It determines
the same autocovariance function as the parameter
$(a^{-1},0,0,a^{-2}\sigma_{\epsilon}^{2})$.

Similarly, if $b\neq 0$, then the parameter $(0,b,0,\sigma_{\epsilon}^{2})$
belongs to two-element class and determines the same autocovariance function
as the parameter $(0,b^{-1},0,b^{-2}\sigma_{\epsilon}^{2})$.

  4. 4.

For $-ab\neq c=0$, the parameter $(a,b,0,\sigma_{\epsilon}^{2})$ makes a class
of its own. Likewise, for $a=b=c=0$, the parameter
$(0,0,0,\sigma_{\epsilon}^{2})$ makes a class of its own. These are
exceptional cases where the autocovariance function uniquely determines the
parameters.

Next proposition shows how the asymmetric case is split into two major
subcases.

######  Proposition 13.

Let $X$ be a stationary field and $\epsilon$ be a collection of zero-mean
equal-variance random variables that satisfy (1).

  1. 1.

If $ab+c=0$, then

$\gamma_{X}(h_{1},h_{2})=\frac{\gamma_{X}(h_{1},0)\gamma_{X}(0,h_{2})}{\gamma_{X}(0,0)}\quad\mbox{for
all $h_{1}$ and $h_{2}$}.$

  2. 2.

If $ab+c\neq 0$ and $1+c^{2}>a^{2}+b^{2}$, then

$\displaystyle\gamma_{X}(1,-1)=\frac{\gamma_{X}(1,0)\gamma_{X}(0,1)}{\gamma_{X}(0,0)}\neq\gamma_{X}(1,1),$
$\displaystyle\gamma_{X}(h_{1},h_{2})=\frac{\gamma_{X}(h_{1},0)\gamma_{X}(0,h_{2})}{\gamma_{X}(0,0)}\quad\mbox{for
all $h_{1}$ and $h_{2}$ such that $h_{1}h_{2}\leq 0$}.$

  3. 3.

If $ab+c\neq 0$ and $1+c^{2}<a^{2}+b^{2}$, then

$\displaystyle\gamma_{X}(1,-1)\neq\frac{\gamma_{X}(1,0)\gamma_{X}(0,1)}{\gamma_{X}(0,0)}=\gamma_{X}(1,1),$
$\displaystyle\gamma_{X}(h_{1},h_{2})=\frac{\gamma_{X}(h_{1},0)\gamma_{X}(0,h_{2})}{\gamma_{X}(0,0)}\quad\mbox{for
all $h_{1}$ and $h_{2}$ such that $h_{1}h_{2}\geq 0$}.$

###### Proof.

The symmetric case $ab+c=0$ has been studied in Section 5.1.1. The desired
equality follows from Proposition 10.

Let us borrow some notation from the proof of Proposition 12, with
$X^{(1)}=X$. Let
$\theta_{3}=(a_{3},b_{3},c_{3},\sigma_{\epsilon,3}^{2})=T_{m}(a,b,c,\sigma_{\epsilon}^{2})$
be a set of parameters obtained from Lemma 4 that satisfy the causality
conditions and let $X^{(3)}$ be a stationary random field defined by this
parameterization; $X^{(3)}$ is causal w.r.t. the respective white noise. The
sign of

$1+c^{2}-a^{2}-b^{2}=\frac{f_{1}f_{4}+f_{2}f_{3}}{2}$

determines the relation between the autocovariance functions of fields $X$ and
$X^{(3)}$. If $1+c^{2}>a^{2}+b^{2}$, then either $X^{(3)}=X$ or
$X^{(3)}_{i,j}=X_{-i,-j}$; in either case the fields $X$ and $X^{(3)}$ have
the same autocovariance function. Otherwise, if $1+c^{2}<a^{2}+b^{2}$, then
either $X^{(3)}_{i,j}=X_{-i,j}$ or $X^{(3)}_{i,j}=X_{i,-j}$; in both cases the
autocovariance functions are one-variable symmetric,
$\gamma_{X^{(3)}}(h_{1},h_{2})=\gamma_{X}(h_{1},-h_{2})$.

As $X^{(3)}$ is causal, due to Proposition 7

$\gamma_{X^{(3)}}(h_{1},h_{2})=\frac{\gamma_{X^{(3)}}(h_{1},0)\gamma_{X^{(3)}}(0,h_{2})}{\gamma_{X^{(3)}}(0,0)}\quad\mbox{if
$h_{1}h_{2}\leq 0$}.$

Hence, and from the relation between $\gamma_{X}$ and $\gamma_{X^{(3)}}$ all
the desired equalities follow.

Now prove the inequalities. Assume that the would-be inequality is actually an
equality, that is
$\gamma_{X}(1,1)=\gamma_{X}(1,0)\gamma_{X}(0,1)/\gamma_{X}(0,0)$ in case
$1+c^{2}>a^{2}+b^{2}$ or
$\gamma_{X}(1,-1)=\gamma_{X}(1,0)\gamma_{X}(0,1)/\gamma_{X}(0,0)$ in case
$1+c^{2}<a^{2}+b^{2}$. Then the autocovariance function of the field $X^{(3)}$
satisfies

$\gamma_{X^{(3)}}(1,1)=\frac{\gamma_{X^{(3)}}(1,0)\gamma_{X^{(3)}}(0,1)}{\gamma_{X^{(3)}}(0,0)}.$

The field $X^{(3)}$ is causal, and formulas in Lemma 5 yield

$a_{3}=\frac{\gamma_{X^{(3)}}(1,0)}{\gamma_{X^{(3)}}(0,0)},\qquad
b_{3}=\frac{\gamma_{X^{(3)}}(0,1)}{\gamma_{X^{(3)}}(0,0)},\qquad
c_{3}=-\frac{\gamma_{X^{(3)}}(1,0)\gamma_{X^{(3)}}(0,1)}{\gamma_{X^{(3)}}(0,0)^{2}}=-a_{3}b_{3}.$

Examining 4 cases, we can verify a similar relation for the original
parameterization, $c=-ab$. ∎

##  6 Pure nondeterminism

###  6.1 Sufficient condition

We borrow the definition of a pure nondeterministic random field from [10].
For the field on a planar lattice, the definition rewrites as follows, in term
of subspaces of the Hilbert space of finite-variance random variables
$L^{2}(\Omega,\mathcal{F},\mathsf{P})$.

######  Definition 3.

Let $X$ be a field of random variables of finite variances. $X$ is called
purely nondeterministic if and only if

$\displaystyle\operatorname*{\overline{\textrm{span}}}_{i,j}X_{i,j}=\operatorname*{\overline{\textrm{span}}}_{i,j}\Bigl{(}$
$\displaystyle\operatorname*{\overline{\textrm{span}}}\\{X_{r,s}:r\leq
i\\}\cap(\operatorname*{\overline{\textrm{span}}}\\{X_{r,s}:r<i\\})^{\bot}\cap\mbox{}$
$\displaystyle\quad\cap\operatorname*{\overline{\textrm{span}}}\\{X_{r,s}:s\leq
j\\}\cap(\operatorname*{\overline{\textrm{span}}}\\{X_{r,s}:s<j\\})^{\bot}\Bigr{)}.$
(40)

Here $\operatorname*{\overline{\textrm{span}}}_{i,j}X_{i,j}$ is the smallest
(closed) subspace of $L^{2}(\Omega,\mathcal{F},\mathsf{P})$ that contains all
observations of the field $X$, while
$\operatorname*{\overline{\textrm{span}}}\\{X_{r,s}:r\leq i\\}$ is a similar
minimal subspace which contains all observations $X_{r,s}$ of the field $X$
with first index $r\leq i$. The outer
$\operatorname*{\overline{\textrm{span}}}$ in the right-hand size of (40) is
the smallest subspace that contains all subspaces used as an argument.

Pure nondeterminism is a sufficient condition for a stationary field to
representable in form (29), where $\epsilon$ is some collection of
uncorrelated variables $\epsilon$, and coefficients $\psi_{k,l}$ satisfy
$\sum_{k=0}^{\infty}\sum_{l=0}^{\infty}\psi_{k,l}^{2}<\infty$. This follows
from [10, Theorem 2.1], while the necessity is easy to verify.

######  Proposition 14.

Let $X$ be a stationary field and $\epsilon$ be a collection of uncorrelated
zero-mean unit-variance random variables that satisfy (1). If
$1+c^{2}>a^{2}+b^{2}$, then $X$ is purely nondeterministic.

###### Proof.

We use notation $D$ from Lemma 1 and notation $f_{1},\ldots,f_{4}$ from (31).

The existence of a stationary solution to (31) imply that $D>0$. As shown in
the proof of Proposition 13, $D>0$ and $1+c^{2}>a^{2}+b^{2}$ holds true if
three cases, $f_{1}>0$, $f_{2}>0$, $f_{3}>0$, $f_{4}>0$ or $f_{1}<0$,
$f_{2}>0$, $f_{3}>0$, $f_{4}<0$ or $f_{1}>0$, $f_{2}<0$, $f_{3}<0$, $f_{3}<0$.

If $f_{1}>0$, $f_{2}>0$, $f_{3}>0$, $f_{4}>0$, then the field $X$ is causal
with respect to white noise $\epsilon$. Equations (1) and (29) imply that

$\operatorname*{\overline{\textrm{span}}}\\{X_{r,s}:r\leq
i\\}=\operatorname*{\overline{\textrm{span}}}{\epsilon_{r,s}:r\leq
i},\quad\;\operatorname*{\overline{\textrm{span}}}\\{X_{r,s}:s\leq
j\\}=\operatorname*{\overline{\textrm{span}}}{\epsilon_{r,s}:s\leq j},$

whence

$\operatorname*{\overline{\textrm{span}}}\\{X_{r,s}:r\leq
i\\}\cap(\operatorname*{\overline{\textrm{span}}}\\{X_{r,s}:r<i\\})^{\bot}\cap\mbox{}\\\
\cap\operatorname*{\overline{\textrm{span}}}\\{X_{r,s}:s\leq
j\\}\cap(\operatorname*{\overline{\textrm{span}}}\\{X_{r,s}:s<j\\})^{\bot}=\operatorname{span}\epsilon_{i,j}.$

Thus,

$\operatorname*{\overline{\textrm{span}}}_{i,j}\Bigl{(}\operatorname*{\overline{\textrm{span}}}\\{X_{r,s}:r\leq
i\\}\cap(\operatorname*{\overline{\textrm{span}}}\\{X_{r,s}:r<i\\})^{\bot}\cap\mbox{}\\\
\cap\operatorname*{\overline{\textrm{span}}}\\{X_{r,s}:s\leq
j\\}\cap(\operatorname*{\overline{\textrm{span}}}\\{X_{r,s}:s<j\\})^{\bot}\Bigr{)}.=\operatorname*{\overline{\textrm{span}}}_{i,j}\epsilon_{i,j}=\operatorname*{\overline{\textrm{span}}}_{i,j}X_{i,j},$

and (40) holds true.

Now consider two other cases, where $-f_{1}$, $f_{2}$, $f_{3}$ and $-f_{4}$
are nonzero and of the same case. In these cases $c\neq 0$, and

$\tilde{\epsilon}_{i,j}=X_{i,j}+\frac{b}{c}X_{i-1,j}+\frac{a}{c}X_{i,j-1}-\frac{1}{c}X_{i-1,j-1}$

is also a collection of uncorrelated zero-mean equal-variance random variables
(to verify this, one can compute the spectral density of $\tilde{\epsilon}$).
Causality condition in Proposition 4 can be easily verified; the stationary
field $X$ is causal w.r.t. white noise $\tilde{\epsilon}$. The field $X$ is
causal in these cases likewise. ∎

###  6.2 Counterexample to Tjøstheim

It would be tempting to use the sufficient condition for pure nondeterminism
from [10, Theorem 3.1]; however, that condition is not correct. We construct a
counterexample.

Let $0<|\theta|<1$. Let $X$ be a stationary field that satisfies the equation

$X_{i-1,j}=\theta X_{i,j-1}+\epsilon_{i,j},$

where $\epsilon$ is a collection of independent random variables with standard
normal distribution, $\epsilon_{i,j}\sim\mathcal{N}(0,1)$.

The spectral density of the field $X$ is

$f_{X}(\nu_{1},\nu_{2})=\frac{1}{|e^{2\pi\mathrm{i}\nu_{1}}-\theta
e^{2\pi\mathrm{i}\nu_{2}}|^{2}}\,.$

The denominator $|e^{2\pi\mathrm{i}\nu_{1}}-\theta
e^{2\pi\mathrm{i}\nu_{2}}|^{2}$ attains only positive values and is a
continuous function. Thus, the field $X$ satisfies the sufficient condition
stated in [10, Theorem 3.1].

In the filed $X$ the diagonals $\\{X_{i,j},\;i+j=n\\}$ are independent for
different $n$. The random variables on each diagonal are jointly distributed
as values of a centered Gaussian AR(1) process with the coefficient $\theta$.

Let

$\tilde{\epsilon}_{i,j}=X_{i,j-1}-\theta X_{i-1,j}.$

Then $\tilde{\epsilon}$ is a collection of uncorrelated zero-mean unit-
variance variables. (Since $\tilde{\epsilon}$ is a Gaussian field, it is a
collection independent variables with distribution $\mathcal{N}(0,1)$.)

The field $X$ can be represented as

$X_{i,j}=\sum_{k=0}^{\infty}\theta^{k}\epsilon_{i+k+1,j-k}=\sum_{k=0}^{\infty}\theta^{k}\tilde{\epsilon}_{i-k,j+k+1}.$

Hence,

$\tilde{\epsilon}_{i,j}=-\theta\epsilon_{i,j}+(1-\theta^{2})\sum_{k=0}^{\infty}\epsilon_{i+k+1,j-k-1}$

These representations imply that

$\operatorname*{\overline{\textrm{span}}}\\{X_{r,s}:r\leq
i\\}=\operatorname*{\overline{\textrm{span}}}\\{\tilde{\epsilon}_{r,s}:r\leq
i\\},\qquad\operatorname*{\overline{\textrm{span}}}\\{X_{r,s}:s\leq
j\\}=\operatorname*{\overline{\textrm{span}}}\\{\epsilon_{r,s}:s\leq j\\}.$

Hence

$\operatorname*{\overline{\textrm{span}}}\\{X_{r,s}:r\leq
i\\}\cap(\operatorname*{\overline{\textrm{span}}}\\{X_{r,s}:r<i\\})^{\bot}\cap\mbox{}\\\
\cap\operatorname*{\overline{\textrm{span}}}\\{X_{r,s}:s\leq
j\\}\cap(\operatorname*{\overline{\textrm{span}}}\\{X_{r,s}:s<j\\})^{\bot}=\operatorname*{\overline{\textrm{span}}}_{s}\tilde{\epsilon}_{i,s}\cap\operatorname*{\overline{\textrm{span}}}_{r}\epsilon_{r,j}.$

Now prove that
$\operatorname*{\overline{\textrm{span}}}_{s}\tilde{\epsilon}_{i,s}\cap\operatorname*{\overline{\textrm{span}}}_{r}\epsilon_{r,j}$
is a trivial subspace. Let
$\zeta\in\operatorname*{\overline{\textrm{span}}}_{s}\tilde{\epsilon}_{i,s}\cap\operatorname*{\overline{\textrm{span}}}_{r}\epsilon_{r,j}$,

$\zeta=\sum_{s=-\infty}^{\infty}k_{s}\tilde{\epsilon}_{i,s}=\sum_{r=-\infty}^{\infty}c_{r}\epsilon_{r,j}.$

Here the coefficients satisfy
$\sum_{s=-\infty}^{\infty}k_{s}^{2}=\sum_{r=-\infty}^{\infty}c_{r}^{2}<\infty$;
the series converge in mean squares. The covariance between
$\tilde{\epsilon}_{i,s}$ and $\epsilon_{r,j}$ is

$\operatorname{\mathsf{E}}\tilde{\epsilon}_{i,s}\epsilon_{r,j}=\begin{cases}-\theta&\mbox{if
$i=r$ and $s=j$,}\\\ \theta^{r-i-1}(1-\theta^{2})&\mbox{if $r-i=s-j>0$,}\\\
0&\mbox{otherwise.}\end{cases}$

Anyway, $\operatorname{\mathsf{E}}\tilde{\epsilon}_{i,s}\epsilon_{r,j}=0$ if
$r-i\neq s-j$, and
$|\operatorname{\mathsf{E}}\tilde{\epsilon}_{i,s}\epsilon_{r,j}|<1$ if
$r-i=s-j$.

Compute $\operatorname{\mathsf{E}}\tilde{\epsilon}_{i,j-n}\zeta$ and
$\operatorname{\mathsf{E}}\zeta\epsilon_{i-n,j}$ using two different
expressions:

$\displaystyle\sum_{s=-\infty}^{\infty}k_{s}\operatorname{\mathsf{E}}\tilde{\epsilon}_{i,j-n}\tilde{\epsilon}_{i,s}=\operatorname{\mathsf{E}}\tilde{\epsilon}_{i,j-n}\zeta=\sum_{r=-\infty}^{\infty}c_{r}\operatorname{\mathsf{E}}\tilde{\epsilon}_{i,j-n}\epsilon_{r,j},$
$\displaystyle\sum_{s=-\infty}^{\infty}k_{s}\operatorname{\mathsf{E}}\tilde{\epsilon}_{i,s}\epsilon_{i-n,j}=\operatorname{\mathsf{E}}\zeta\epsilon_{i-n,j}=\sum_{r=-\infty}^{\infty}c_{r}\operatorname{\mathsf{E}}\epsilon_{r,j}\epsilon_{i-n,j}.$

In the series, the only nonzero term might be where $s=j-n$ and $r=i-n$. Thus,

$k_{j-n}=c_{i-n}\operatorname{\mathsf{E}}\tilde{\epsilon}_{i,j-n}\epsilon_{i-n,j},\quad
k_{j-n}\operatorname{\mathsf{E}}\tilde{\epsilon}_{i,j-n}\epsilon_{i-n,j}=c_{i-n}.$
(41)

As $|\operatorname{\mathsf{E}}\tilde{\epsilon}_{i,j-n}\epsilon_{i-n,j}|<1$,
(41) imply $k_{j-n}=c_{i-n}=0$. As this holds true for all integer $n$,
$k_{s}=c_{r}=0$ for all $s$ and $r$, and $\zeta=0$ almost surely. Thus,

$\operatorname*{\overline{\textrm{span}}}\\{X_{r,s}:r\leq
i\\}\cap(\operatorname*{\overline{\textrm{span}}}\\{X_{r,s}:r<i\\})^{\bot}\cap\mbox{}\\\
\cap\operatorname*{\overline{\textrm{span}}}\\{X_{r,s}:s\leq
j\\}\cap(\operatorname*{\overline{\textrm{span}}}\\{X_{r,s}:s<j\\})^{\bot}=\\{0\\},$

and the right-hand side of (40) is the trivial subspace. Thus, the random
field $X$ is not purely nondeterministic.

##  7 Conclusion

We considered AR(1) model on a plane. We found conditions (in terms of the
regression coefficients) under which the autoregressive equation has a
stationary solution $X$. As for the autocovariance function of the stationary
solution $X$, we presented a simple formula for it at some points, and proved
Yule–Walker equations. These allow to compute the autocovariance function
recursively at all points.

We found conditions under which the stationary solution to the autoregressive
equation satisfies the causality condition with respect to the underlying
white noise. These conditions also appear to be sufficient conditions for
stability of the deterministic problem of solving a recursive equations in a
quadrant, with preset values on the border of the quadrant.

We described sets of parameters (the coefficients and the variance of the
underlying white noise) where different parameters determine the same
autocovariance function of the stationary solution.

We found sufficient conditions for the stationary solution $X$ to be a pure
nondeterministic random field. This condition is related to the causality
condition with respect to some (non-fixed) white noise, which is called
innovations; in particular, the innovations need not coincide with the white
noise in the autoregressive equation.

The causality condition and pure-nondeterministic property seem to be too
restrictive because the coordinate-wise order between coordinates of
$\epsilon$ and $X$ in the representation (29) is a partial order. More general
representation with lexical order, which is a total order, is suggested in
[14, Section 6]:

$X_{i,j}=\sum_{k=1}^{\infty}\sum_{l=-\infty}^{\infty}\psi_{k,l}X_{i-k,j-l}+\sum_{l=0}^{\infty}\psi_{0,l}X_{i,j-l}.$

Further discussion can be found in [4, 11].

With exception for the identifiability topic, we did not consider the
estimation at all.

Most results of this paper are necessary and sufficient conditions for some
properties of the stationary field $X$. However, for pure nondeterminism of
$X$ and for stability of the deterministic equation, we obtained only
sufficient conditions. Obtaining necessary conditions is an interesting open
problem.

##  8 Appendix: Auxiliary results

###  8.1 Existence of a random field with given spectral density

Next lemma is used on the proof of Proposition 1. It was a modification of the
similar result stated for stochastic processes [9].

######  Lemma 6.

Let $f(\nu_{1},\nu_{2})$ be an even integrable function
$\left[-\frac{1}{2},\frac{1}{2}\right]^{2}\to[0,\infty]$, that is

$\displaystyle f(\nu_{1},\nu_{2})=f(-\nu_{1},-\nu_{2})\geq 0\quad\mbox{for all
$\nu 1,\nu_{2}\mathbin{\in}[-\frac{1}{2},\frac{1}{2}]$},$
$\displaystyle\iint_{[-1/2,\,1/2]^{2}}f(\nu_{1},\nu_{2})\,d\nu_{1}\,d\nu_{2}<\infty.$

Then there exists a Gaussian stationary field
$\\{X_{i,j},\;i,j\mathbin{\in}\mathbb{Z}\\}$ on some (specially constructed)
probability space that has spectral density $f(\nu_{1},\nu_{2})$.

###### Proof.

Let us construct a probability space with two Brownian fields on a plane
$\\{W_{k}(t,s),\;t,s\mathbin{\in}[0,1]\\}$ (a zero-mean Gaussian field with
covariance function
$\operatorname{cov}(W_{k}(t_{1},s_{1}),\,W_{k}(t_{2},s_{2}))=\min(t_{1},t_{2})\min(s_{1},s_{2})$
$k=1,\,2$).

Extend the function $X$ by periodicity,
$f(\nu_{1},\nu_{2})=f(\nu_{1}-1,\nu_{2})$ if $\frac{1}{2}<\nu_{1}\leq 1$, and
then $f(\nu_{1},\nu_{2}-1)=f(\nu_{1}-1,\nu_{2})$ if $\frac{1}{2}<\nu_{2}\leq
1$.

Denote

$\displaystyle X_{i,j}$
$\displaystyle=\iint_{[0,1]^{2}}\cos(2\pi(i\nu_{1}+j\nu_{2}))\sqrt{f(\nu_{1},\nu_{2})}\,d^{2}W_{1}(\nu_{1},\nu_{2})+\mbox{}$
$\displaystyle+\iint_{[0,1]^{2}}\sin(2\pi(i\nu_{1}+j\nu_{2}))\sqrt{f(\nu_{1},\nu_{2})}\,d^{2}W_{2}(\nu_{1},\nu_{2}).$

The constructed field is zero-mean Gaussian as the integral of nonstochastic
kernel with respect Gaussian field. The autocovariance function is as
expected,

$\displaystyle\operatorname{cov}(X_{i,j},X_{i+h_{1},j+h_{2}})$
$\displaystyle=\iint_{[0,1]^{2}}\cos(2\pi(h_{1}\nu_{1}+h_{2}\nu_{2}))f(\nu_{1},\nu_{2})\,d\nu_{1}\,d\nu_{2}=$
$\displaystyle=\iint_{[-1/2,1/2]^{2}}\exp(2\pi\mathrm{i}(h_{1}\nu_{1}+h_{2}\nu_{2}))f(\nu_{1},\nu_{2})\,d\nu_{1}\,d\nu_{2}.$

∎

###  8.2 Summation

The next two propositions are used implicitly when we use double-series
notation.

######  Proposition 15.

Let $\\{\xi_{i,j},\;i,j\in\mathbb{Z}\\}$ be a collection of random variables
bounded in mean squares,
$\sup_{i,j}\operatorname{\mathsf{E}}\xi_{i,j}^{2}<\infty$. Let
$\\{\xi_{i,j},\;i,j\in\mathbb{Z}\\}$ be a collection of real numbers such that
${\sum\sum}_{i,j=-\infty}^{\infty}|a_{i,j}|<\infty$. Then the double series
${\sum\sum}_{i,j=-\infty}^{\infty}a_{i,j}\xi_{i,j}$ converges in mean squares
and almost surely. The limit it the same for both types of convergence (up to
equal-almost-surely equivalence); it also does not depend on the order the
terms of the double series are added.

######  Remark 3.

  1. 1.

Under conditions of Proposition 15, the double series
${\sum\sum}_{i,j=-\infty}^{\infty}a_{i,j}\xi_{i,j}$ converges in mean and in
probability, as well — to the same limit.

  2. 2.

The iterated series
$\sum_{i=-\infty}^{\infty}\left(\sum_{j=-\infty}^{\infty}a_{i,j}\xi_{i,j}\right)$
converges to the same limit.

  3. 3.

$\operatorname{\mathsf{P}}\left({\sum\sum}_{i,j=-\infty}^{\infty}|a_{i,j}\xi_{i,j}|<\infty\right)=1$,
and whenever ${\sum\sum}_{i,j=-\infty}^{\infty}|a_{i,j}\xi_{i,j}|<\infty$, the
double series ${\sum\sum}_{i,j=-\infty}^{\infty}a_{i,j}\xi_{i,j}$ converges,
and its limit does not depend on the order its terms are added.

######  Proposition 16.

Let $\\{\epsilon_{i,j},\;i,j\in\mathbb{Z}\\}$ be a collection of random
uncorrelated variables with zero mean and same variance
$\operatorname{var}\epsilon_{i,j}=\sigma_{\epsilon}^{2}<\infty$. Let
$\\{\xi_{i,j},\;i,j\in\mathbb{Z}\\}$ be a collection of real numbers such that
${\sum\sum}_{i,j=-\infty}^{\infty}a_{i,j}^{2}<\infty$. Then the double series
${\sum\sum}_{i,j=-\infty}^{\infty}a_{i,j}\xi_{i,j}$ converges in mean squares.
The limit does not depend on the order the terms of the double series are
added: it is the same up to equal-almost-surely equivalence.

######  Remark 4.

  1. 1.

Under conditions of Proposition 16, the double series
${\sum\sum}_{i,j=-\infty}^{\infty}a_{i,j}\xi_{i,j}$ converges in mean and in
probability, as well.

  2. 2.

The iterated series
$\sum_{i=-\infty}^{\infty}\left(\sum_{j=-\infty}^{\infty}a_{i,j}\xi_{i,j}\right)$
converges in mean squares to the same limit.

Next lemma is used in Proposition 4.

######  Lemma 7.

Let $f\in C(\mathbb{R}^{2})$ be a biperiodic function of two variables with
continuous mixed derivative:

$\displaystyle f(x,y)=f(x+1,y)=f(x,y+1)\quad\mbox{for all
$x,y\in\mathbb{R}$},$ $\displaystyle f\in
C(\mathbb{R}^{2}),\quad\frac{\partial f}{\partial x}\in
C(\mathbb{R}^{2}),\quad\frac{\partial f}{\partial y}\in
C(\mathbb{R}^{2}),\quad\frac{\partial^{2}f}{\partial x\,\partial y}\in
C(\mathbb{R}^{2}).$

Then Fourier coefficients of the function $f$ are summable:

$\sum_{k=-\infty}^{\infty}\sum_{l=-\infty}^{\infty}\left|\int_{0}^{1}\int_{0}^{1}\exp(2\pi\mathrm{i}\,(kx+ly))\,f(x,y)\,dx\,dy\right|<\infty.$

###### Proof.

Denote the Fourier coefficients $c_{k,l}$:

$c_{k,l}=\int_{0}^{1}\int_{0}^{1}\exp(2\pi\mathrm{i}\,(kx+ly))\,f(x,y)\,dx\,dy.$

The coefficients $c_{k,l}$ are also equal

$\displaystyle c_{k,l}=\frac{\mathrm{i}}{2\pi
k}\int_{0}^{t}\int_{0}^{1}\exp(2\pi\mathrm{i}\,(kx+ly))\,f^{\prime}_{1}(x,y)\,dx\,dy\qquad\mbox{if
$k\neq 0$};$ $\displaystyle c_{k,l}=\frac{\mathrm{i}}{2\pi
l}\int_{0}^{t}\int_{0}^{1}\exp(2\pi\mathrm{i}\,(kx+ly))\,f^{\prime}_{2}(x,y)\,dx\,dy\qquad\mbox{if
$l\neq 0$};$ $\displaystyle
c_{k,l}=\frac{-1}{4\pi^{2}kl}\int_{0}^{t}\int_{0}^{1}\exp(2\pi\mathrm{i}\,(kx+ly))\,f^{\prime\prime}_{12}(x,y)\,dx\,dy\qquad\mbox{if
$k\neq 0$ and $l\neq 0$},$

where $f^{\prime}_{1}(x,y)$, $f^{\prime}_{2}(x,y)$, and
$f^{\prime\prime}_{12}(x,y)$ are partial and mixed derivatives; here the
periodicity of the function $f$ is used. The coefficients $c_{k,l}$ allow
bounding:

$\displaystyle c_{k,0}\leq\frac{1}{2\pi|k|}\max_{x,y}\left|\frac{\partial
f(x,y)}{\partial x}\right|\qquad\mbox{if $k\neq 0$};$ $\displaystyle
c_{0,l}\leq\frac{1}{2\pi|l|}\max_{x,y}\left|\frac{\partial f(x,y)}{\partial
y}\right|\qquad\mbox{if $l\neq 0$};$ $\displaystyle
c_{k,l}\leq\frac{1}{4\pi^{2}|kl|}\max_{x,y}\left|\frac{\partial^{2}f(x,y)}{\partial
y\partial x}\right|\qquad\mbox{if $k\neq 0$ and $l\neq 0$},$

which implies the summability. ∎

###  8.3 Integration

######  Lemma 8.

If $A>0$, $B$ and $C$ are real numbers such that $A^{2}>B^{2}+C^{2}$, then

$\int_{-1/2}^{1/2}\frac{dt}{A+B\cos(2\pi t)+C\sin(2\pi
t)}=\frac{1}{\sqrt{A^{2}-B^{2}-C^{2}}}.$ (42)

If $A$ and $B$ are real numbers, $A>|B|$, and $n$ is integer, then

$\int_{-1/2}^{1/2}\frac{\exp(2\pi\mathrm{i}nt)}{A+B\cos(2\pi
t)}\,dt=\frac{\alpha^{|n|}}{\sqrt{A^{2}-B^{2}}},$ (43)

where

$\alpha=\frac{-B}{A+\sqrt{A^{2}-B^{2}}}=\begin{cases}0&\mbox{if $B=0$},\\\
(-A+\sqrt{A^{2}-B^{2}})/B&\mbox{if $B\neq 0$}.\end{cases}$

In (43), if $\alpha=0$ and $n=0$, then $\alpha^{|n|}=1$ by convention.

###### Proof.

That is easy to check that the antiderivative in (42) is

$\int\frac{dt}{A+B\cos(2\pi t)+C\sin(2\pi t)}=\\\
=\frac{1}{\pi\sqrt{A^{2}-B^{2}-C^{2}}}\,\arctan\\!\left(\frac{(A-B)\tan(\pi
t)+C}{\sqrt{A^{2}-B^{2}-C^{2}}}\right)+\mathrm{const},$

whence (42) follows. (Notice that $A-B>0$.)

In (43), the function

$\frac{1}{A+B\cos(2\pi t)}=\frac{1+\alpha^{2}}{A\,(1-2\alpha\cos(2\pi
t)+\alpha^{2})}=\frac{1+\alpha^{2}}{A}\,\frac{1}{|1-\alpha\exp(2\pi\mathrm{i}t)|^{2}}$

is the spectral density of AR(1) stationary autoregressive process
$X_{k}=\alpha X_{k-1}+\epsilon_{k}$ with white-noise variance
$\operatorname{var}\epsilon_{k}=(1+\alpha^{2})/A$. (Conditions imply that
$|\alpha|<1$.) The integral is the autocovariance function of the process. It
equals

$\int_{-1/2}^{1/2}\frac{\exp(2\pi\mathrm{i}nt)}{A+\cos(2\pi
t)}\,dt=\operatorname{cov}(X_{k+n},X_{k})=\frac{1+\alpha^{2}}{A}\,\frac{\alpha^{|n|}}{1-\alpha^{2}}=\frac{\alpha^{|n|}}{\sqrt{A^{2}-B^{2}}}\,;$

here equality

$\frac{1+\alpha^{2}}{1-\alpha^{2}}=\frac{A}{\sqrt{A^{2}-B^{2}}}$

is used. ∎

######  Lemma 9.

If $A$ and $B$ are complex numbers, $|A|\neq|B|$, and $n$ is an integer
number, then

$\int_{-1/2}^{1/2}\frac{e^{2\pi{\mathrm{i}}n\nu}}{A-Be^{2\pi{\mathrm{i}}\nu}}\,d\nu=\begin{cases}A^{n-1}B^{-n}&\mbox{if
$n\leq 0$ and $|A|>|B|$,}\\\ 0&\mbox{if $n\leq 0$ and $|A|<|B|$,}\\\
0&\mbox{if $n\geq 1$ and $|A|>|B|$,}\\\ -A^{n-1}B^{-n}&\mbox{if $n\geq 1$ and
$|A|<|B|$.}\\\ \end{cases}$ (44)

###### Proof.

For $n=1$, the integral can be rewritten as the contour integral,

$\int_{-1/2}^{1/2}\frac{e^{2\pi{\mathrm{i}}\nu}}{A-Be^{2\pi{\mathrm{i}}\nu}}\,d\nu=\frac{1}{2\pi\mathrm{i}}\oint_{z=\exp(2\pi\mathrm{i}\nu)}\frac{1}{A-Bz}\,dz$

and the residue formula is applicable. Other cases can be reduced to the case
$n=1$ recursively. ∎

######  Lemma 10.

Let $a$, $b$ and $c$ be real numbers such that $1-a-b-c>0$, $1-a+b+c>0$,
$1+a-b+c>0$ and $1+a+b-c>0$. Let $n_{1}$ and $n_{2}$ be nonnegative integers.
Then

$\int_{-1/2}^{1/2}\frac{(a+ce^{2\pi\mathrm{i}\nu})^{n_{1}}e^{-2\pi\mathrm{i}n_{2}\nu}}{(1-be^{2\pi\mathrm{i}\nu})^{n_{1}+1}}\,d\nu=\sum_{k=0}^{\min(n_{1},n_{2})}\binom{n_{1}}{k}\binom{n_{2}}{k}a^{n_{1}-k}b^{n_{2}-k}(ab+c)^{k}.$

###### Proof.

Conditions of the lemma imply that

$\displaystyle 1-b=\frac{(1-a-b-c)+(1+a-b+c)}{2}>0,$ $\displaystyle
1+b=\frac{(1-a+b+c)+(1+a+b-c)}{2}>0,$ $\displaystyle|b|<1.$

Let $k\geq 0$ and $n$ be integer numbers. By binomial formula,

$\frac{1}{(1-bz)^{k+1}}=\sum_{m=0}^{\infty}\binom{-k-1}{m}(-bz)^{m}=\sum_{m=0}^{\infty}\binom{k+m}{m}(bz)^{m}$
(45)

for all complex $z$ such that $|bz|<1$, in particular, for all $z$ such that
$|z|\leq 1$. Here $\binom{-k-1}{m}=\frac{(-k-1)\cdots(-k-m)}{1\cdots
m}=(-1)^{m}\binom{k+m}{m}$ is a coefficient in a binomial series. The rational
function $z^{n}(1-bz)^{-k-1}$ have a singularity a point $1/b$ and potential
singularity (if $n<0$) at point $0$. The point $1/b$ lies outside the unit
circle, and point 0 lies inside the unit circle. The expansion (45) implies
that

$\displaystyle\operatorname{Res}_{z=0}\left(\frac{z^{n}}{(1-bz)^{k+1}}\right)$
$\displaystyle=\operatorname{Res}_{z=0}\left(\sum_{m=0}^{\infty}\binom{k+m}{m}b^{m}z^{m+n}\right)=$
$\displaystyle=\begin{cases}0&\mbox{if $n\geq 0$},\\\
\binom{k-n-1}{-n-1}b^{-n-1}&\mbox{if $n\leq-1$ .}\end{cases}$

By the residue formula,

$\displaystyle\oint\frac{z^{n}}{(1-bz)^{k+1}}\,dz$
$\displaystyle=2\pi\mathrm{i}\operatorname{Res}_{z=0}\left(\frac{z^{n}}{(1-bz)^{k+1}}\right)=$
$\displaystyle=\begin{cases}0&\mbox{if $n\geq 0$},\\\
2\pi\mathrm{i}\binom{k-n-1}{-n-1}b^{-n-1}&\mbox{if $n\leq-1$,}\end{cases}$

where the integral is taken along the unit-circle contour.

Recall that $n_{2}\geq 0$ is integer. Then

$\int_{-1/2}^{1/2}\frac{\exp(2\pi\mathrm{i}(k-n_{2})\nu_{2})}{(1-be^{2\pi\mathrm{i}\nu_{2}})^{1+k}}\,d\nu_{2}=\frac{1}{2\pi\mathrm{i}}\oint\frac{z^{k-n_{2}-1}}{(1-bz)^{k+1}}\,dz=\\\
=\begin{cases}0&\mbox{if $k>n_{2}$,}\\\
\binom{n_{2}}{n_{2}-k}b^{n_{2}-k}&\mbox{if $k\leq
n_{2}$}\end{cases}=\begin{cases}0&\mbox{if $k>n_{2}$,}\\\
\binom{n_{2}}{k}b^{n_{2}-k}&\mbox{if $k\leq n_{2}$.}\end{cases}$ (46)

The constant $n_{1}\geq 0$ is also integer, and by binomial formula

$\displaystyle(a+ce^{2\pi\mathrm{i}\nu})^{n_{1}}$
$\displaystyle=((ab+c)e^{2\pi\mathrm{i}\nu}+a(1-be^{2\pi\mathrm{i}\nu}))^{n_{1}}=$
$\displaystyle=\sum_{k=0}^{n_{1}}\binom{n_{1}}{k}(ab+c)^{k}e^{2\pi\mathrm{i}k\nu}a^{n_{1}-k}(1-be^{2\pi\mathrm{i}\nu})^{n_{1}-k}.$
(47)

Finally, with use of (47) and (46),

$\displaystyle\int_{-1/2}^{1/2}\frac{e^{-2\pi\mathrm{i}n_{2}\nu}(a+ce^{2\pi\mathrm{i}\nu})^{n_{1}}}{(1-be^{2\pi\mathrm{i}\nu})^{n_{1}+1}}\,d\nu=$
$\displaystyle=\int\limits_{-1/2}^{1/2}\\!\\!\frac{e^{-2\pi\mathrm{i}n_{2}\nu}}{(1{-}be^{2\pi\mathrm{i}\nu})^{n_{1}+1}}\sum_{k=0}^{n_{1}}\\!\binom{n_{1}}{k}(ab{+}c)^{k}e^{2\pi\mathrm{i}k\nu}a^{n_{1}-k}(1{-}be^{2\pi\mathrm{i}\nu})^{n_{1}-k}d\nu=$
$\displaystyle=\sum_{k=0}^{n_{1}}\binom{n_{1}}{k}a^{n_{1}-k}(ab+c)^{k}\int_{-1/2}^{1/2}\frac{e^{2\pi\mathrm{i}(k-n_{2})\nu}}{(1-be^{2\pi\mathrm{i}\nu})^{1+k}}\,d\nu=$
$\displaystyle=\sum_{k=0}^{\min(n_{1},\,n_{2})}\binom{n_{1}}{k}a^{n_{1}-k}(ab+c)^{k}\binom{n_{2}}{k}b^{n_{2}-k}.\qed$

## References

  * [1] S. Baran, G. Pap, M. C. A. van Zuijlen. Asymptotic inference far a nearly unstable sequences of stationary spatial AR models. A Journal of Theoretical and Applied Statistics, 2004, Vol. 38, No. 6, 465–482.
doi:10.1080/02331880412331319297

  * [2] N. A. C. Cressie. Statistics for Spatial Data. Wiley, 1993.
doi:10.1002/9781119115151

  * [3] A. K. Jain. Advances in mathematical models for image processing. Proceedings of the IEEE, 1981, Vol. 69, No. 5, 502–528.
doi:10.1109/PROC.1981.12021

  * [4] G. Kallianpur. Some remarks of the purely nondeterministic property of second-order random fields. In Stochastic Differential Systems, Springer, 1981, pp. 98–109. doi:10.1007/BFb0006413
  * [5] R. J. Martin. A subclass of lattice processes applied to a problem in planar sampling. Biometrika, 1979, Vol. 66, No. 2, 209–217.
doi:10.2307/2335651

  * [6] R. J. Martin. The use of time-series models and methods in the analysis of agricultural field trials. Communications in Statistics — Theory and Methods, 1990, Vol. 19, No. 1, pp. 55–81. doi:10.1080/03610929008830187
  * [7] R. J. Martin. Some results on unilateral ARMA lattice processes. Journal of Statistical Planning and Inference, 1996, Vol. 50, No. 3, 395–411.
doi:10.1016/0378-3758(95)00066-6

  * [8] D. K. Pickard. Unilateral Markov fields. Advances in Applied Probability, 1980, Vol. 12, No. 3, 655–671.
  * [9] R. H. Shumway, D. S. Stoffer. Time Series Analysis and Its Applications: With R Examples. Springer, New York, 2011.
  * [10] D. Tjøstheim. Statistical spatial series modelling. Advances in Applied Probability, 1978, Vol. 10, No. 1, 130–154. doi:10.2307/1426722
  * [11] D. Tjøstheim. Statistical spatial series modelling II. Some further results on unilateral lattice processes. Advances in Applied Probability, 1983, Vol. 15, No. 3, 562–584. doi:10.2307/1426619
  * [12] E. M. Tory and D. K. Pickard. Unilateral Gaussian fields. Advances in Applied Probability, 1992, Vol. 24, No. 1, pp. 95–112. doi:10.2307/1427731
  * [13] D. J. Unwin, L. W. Hepple. The statistical analysis of spatial series. Journal of the Royal Statistical Society: Series D, 1974, Vol. 23, No. 3/4, 211–227. doi:10.2307/2987581
  * [14] P. Whittle. On stationary processes in the plane. Biometrika, 1954, Vol. 41, No. 3–4, 434–449. doi:10.1093/biomet/41.3-4.434