WOLFRAM|DEMONSTRATIONS PROJECT

Estimating a Centered Matérn (1) Process: Three Alternatives to Maximum Likelihood via Conjugate Gradient Linear Solvers

​
inverse-range parameter θ
3
6
12
24
48
96
192
seed for data generation
843
dataset size
256
512
1024
2048
4096
8192
The setting is the same as in the Demonstration "Three Alternatives to the Likelihood Maximization for Estimating a Centered Matérn (3/2) Process", except that the "differentiability" parameter (often denoted by
ν
) is fixed here to 1. This value 1 for
ν
is classically used in two dimensions following the seminal work of Whittle [1]. The corresponding autocorrelation function has the simple expression
ρ(t,θ)=
[x(s+t)x(s)]
2
τ
=(θt)
K
1
(θt)
(the process variance

2
x(s)

being denoted by
2
τ
), where
K
1
(·)
is a modified Bessel function of the second kind (see the Details section in the help page for the BesselK function) and
θ>0
is still called the inverse-range parameter,
-1
θ
being the "range" (or "support") of
ρ(·,θ)
.
One important difference with the case
ν=3/2
is that the considered process, even regularly sampled (here observed at
t
in the set
{0,δ,2δ,…,(n-2)δ,(n-1)δ=1}
), no longer coincides with any ARMA time series. Fast Fourier transforms (FFT) are then invoked for the simulation of such a series.
The same three methods as in the case
ν=3/2
are implemented here for estimating
τ
and
θ
, namely: Cressie's weighted least-squares method [4], the "hybrid" Zhang–Zimmerman method [5], and the recently proposed "CGEM-EV" method [6], here in the case of no measurement error. For each of the three methods, the estimate of the so-called "microergodic coefficient", which is now
c=
2
τ
2
θ
(the general expression for
c
being
c=
2
τ
2ν
θ
), is naturally defined as the square of the product of the estimates of
τ
and
θ
.
For all the simulated time series, the true variance
2
τ
is fixed to 1 (this choice is without loss of generality), and various
θ
can be tried for the true value of the inverse-range parameter. The legends for the table of the numerical results are similar to the ones in the above-mentioned Demonstration for the case
ν=3/2
.
Let
R(θ)
denote the candidate correlation matrix of
x={x(0),x(δ),…,x((n-2)δ),x(1)}
, with inverse-range
θ
, that is, the
n×n
matrix whose element in the
th
i
row and
th
j
column is
[R(θ)]
i,j
=ρ(i-jδ,θ)
. The only difficulty regarding implementation in the case of no measurement error, of both the hybrid method and CGEM-EV, is the computation of the quadratic form
1
n
x,
-1
R(θ)
x
, the so-called "Gibbs energy of
x
associated with a given
θ
".
The computation of
-1
R(θ)
x
uses a conjugate-gradient (CG) solver preconditioned by a classical factored sparse approximate inverse (FSAI) preconditioning (see [7] for a recent survey), each product by
R(θ)
being obtained via FFT from the standard embedding of
R(θ)
in a circulant
2n×2n
matrix.
It is important to notice that
R(θ)
becomes ill-conditioned for large
n
and
θ
decreasing toward the lower-bound
θ=1.5
used here (as in the Demonstration for
ν=3/2
); so using a preconditioner is a mandatory requirement to accelerate the CG convergence.
It is observed here that this implementation is quite fast, even for
n=8196
.