WOLFRAM|DEMONSTRATIONS PROJECT

Three Alternatives to the Likelihood Maximization for Estimating a Centered Matérn (3/2) Process

The problem of estimating the parameters of a stationary Gaussian process whose autocorrelation function belongs to the Matérn class appears in many contexts (e.g. [1, 2]). This Demonstration considers the Matérn subclass on



, which has its "differentiability" parameter fixed to 3/2, an often-used value (e.g. [3], [1]). Also, assume that the constant mean of the process is zero. To be precise, such a process then satisfies the stochastic differential equation (SDE)

dx'(t)(-2θx'(t)-

x(t))dt+

dw(t)

, where

x'(·)

denotes the derivative of

x(·)

(in the mean-square sense),

follows a standard Wiener process, and

is a positive parameter. In geostatistics,

is called the inverse-range parameter, because

-1

is the "range" (or "support") of the corresponding autocorrelation function

ρ(·,θ)

(see the Details section). The process variance, denoted by

, is known to be given by

τ=

c/(4

)

, and the parameter

c=4

is called the "microergodic coefficient". Assume that

observations of one realization

x(·)

are given; more precisely, the values

x(t)

for

in the set

{0,δ,2δ,⋯,(n-2)δ,(n-1)δ=1}

. This Demonstration thus has some common points with the Demonstration "Estimating a Centered Ornstein-Uhlenbeck Process under Measurement Errors".

This Demonstration is restricted to the case of no measurement error. In the Demonstration mentioned in the last paragraph, the exact maximum-likelihood (ML) method is implemented for this case (indeed the exact minimizer is then explicit for the Ornstein–Uhlenbeck case). In this Demonstration, an exact ML method is not implemented because, in spite of the existence of a fast order-

algorithm to compute the likelihood criterion, its global maximization (even for its "profile" version) is not an easy task, unless one is content to find an approximate solution by a rough grid search.

A classic approach (especially in geostatistics) to the problem of estimating the two parameters

and

is "variogram fitting." Here this consists of trying to minimize the differences between the so-called variogram, that is, the list

n-j

n-1-j

∑

i=0

[x(iδ+jδ)-x(iδ)]

(where the lags

jδ

describe the set

{δ,2δ,⋯,(n-1)δ}

of possible lags) and the theoretical model

[1-ρ(·,θ)]

(evaluated at the same lags). This is known as Cressie's weighted least-squares method [4], when these differences are summarized by an appropriately weighted least-squares distance (which is minimized with respect to

and

; it is relatively easy to see that this can be reduced to a one-dimensional minimization in

, for which a global search over a grid is very fast). In this Demonstration, a grid of 313 values is chosen for

equispaced, in logarithmic scale, between 1.5 and 2000.

A second method, called a "hybrid" method, which was first advocated in [5], is to retain the above least-squares estimate of

, denoted by



, but to replace the associated least-squares estimate of

by the constrained-ML estimate defined as the maximizer of the likelihood when the unknown

is constrained to be



. Indeed, it is well-known that this maximizer, say



, can be expressed (with

〈x,y〉

denoting the scalar product of

and

) in closed form by



x,

-1

R





x

, where

is the

n×1

column vector

{x(0),x(δ),⋯,x((n-2)δ),x(1)}

and

R(θ)

is the correlation matrix of

when the inverse-range parameter is

, that is, the

n×n

matrix whose element in the

row and

column is

[R(θ)]

i,j

=ρ(i-jδ,θ)

. It has been shown (see [3] and the survey in [2]) that there is no loss of asymptotic ("infill" framework) statistical efficiency for the resulting hybrid estimate



of the microergodic parameter, when compared to its ML estimate.

For any

θ>0

, the quadratic form

x,

-1

R(θ)

x

is called "the candidate Gibbs energy of

associated with

." It is relatively easy to see that computing such a Gibbs energy has only a cost of order

, in the setting of this Demo. In fact, to simplify the code, this computation is directly constructed via two calls of the Mathematica's built-in LogLikelihood[ARMAProcess[{a1,a2},{b1},1],.] function.

Note that the grid search (used to approximate



) is restricted to values of

greater than 1.5 to prevent ill-conditioned computations that may occur when computing Gibbs energies for

less than this bound. Such an "ill-conditioning" phenomenon is well-known in geostatistics for cases with strong correlation and large data size.

The third alternative to ML considered in this Demonstration is the recently proposed "CGEM-EV" method [6], whose definition is quite simple in the case of no measurement errors: first

is simply estimated by the (unbiased) empirical variance



n-1

∑

i=0

[x(iδ)]

; second, the closest match between model Gibbs energy

x,

-1

R(θ)

x

and

is invoked to estimate

. It is easy to show that

x,

-1

R(θ)

x=

is an unbiased estimating equation; indeed, ensemble averaging of the left-hand term coincides with

when

is set to its true value. Stronger properties are studied in [6]. A simple fixed-point algorithm, always starting from the "small" value 1.5 as a first guess for

(see the above remark on ill-conditioning in computing candidate Gibbs energies) is used here to search for an approximate root greater than 1.5 (provided the root exists) for this estimating equation, and it proves to be successful with quite fast convergence (for example, for

n=512

the number of iterations, and thus of candidate Gibbs energies to be computed, is always between 4 and 11) in all the settings that can be chosen in this Demonstration.