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

Initializing live version
Download to Desktop

Requires a Wolfram Notebook System

Interact on desktop, mobile and cloud with the free Wolfram Player or other Wolfram Language products.

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 (where the lags describe the set of possible lags) and the theoretical model (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.

For any , the quadratic form 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 ; second, the closest match between model Gibbs energy and is invoked to estimate . It is easy to show that 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 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.


Contributed by: Didier A. Girard (January 2015)
(CNRS-LJK and Université Grenoble Alpes)
Open content licensed under CC BY-NC-SA



It can be shown (e.g. in the Papoulis book cited in [1]) that the stationary solution of the above second-order SDE is a centered Gaussian process whose autocorrelation function has the simple expression , which is the more classical definition of a Matérn (3/2) process (here denotes the ensemble averaging, i.e. from infinitely repeated simulations of the process under the true model).

By choosing equispaced abscissae, the simulations are greatly facilitated, since the observations then coincide to a particular ARMA(2,1) time series, the transformation giving the ARMA parameters from and requiring only simple and well-known computations [7].

For all the simulated time series, the variance is fixed to 1, called "true (unknown) value of variance" (this choice is without loss of generality), and various can be tried for the true value of the inverse-range parameter. The resulting true value for the microergodic parameter is displayed at the bottom-left (third line) of the table of the numerical results.

Snapshot 1: Selecting 12 as the true inverse-range parameter and , you can first compare the estimates produced by the classical variogram least-squares fitting (first column) versus Zhang and Zimmerman hybrid estimates (second column; note that the estimate of is also displayed in this second column, even though it is unchanged from the first column): for this seed (827), the variance estimate (second line) is observed to be degraded by the hybrid method, but the estimate of is much improved. Looking now at the third column, which displays the results produced by the third estimation method, CGEM-EV, you observe an improvement for the estimate of and a little change for the estimate of , as compared to the hybrid method.

By changing the seed (and thus a new is used to generate the data), you can be convinced that these remarks about the three estimates of are not an accident. But you can also observe a rather large variability (i.e. from seed to seed) for each of the three estimates of the variance (or of ); this is much larger than the variability of the estimates of the microergodic coefficient , which can be observed for the hybrid or CGEM-EV methods. Notice that this observation is in very good agreement with the known theory about the estimation of the variance, the inverse-range parameter, and their product (see [6], [3], [2] and the references therein).

Snapshot 2: Now selecting 3 as true , a clear improvement (in the estimation of or of ) by using CGEM-EV in place of the hybrid method can be observed: indeed, for such "large range" settings, the least-squares fitting method produces even more variable estimates (and it often produces estimates at the boundary of the grid search, here 1.5). Increasing (for example, setting ), you can observe that the estimate of is also improved, although to a lesser extent.

Snapshot 3: Staying at and now selecting 100 as true (this is a "rather small range" setting), a ranking of CGEM-EV and the hybrid method based on their statistical performances becomes much more difficult. One also observes that the relative accuracy of both these estimates of or of is clearly increased as compared to the above "large range" settings.


[1] C. E. Rasmussen and C. K. I. Williams, Gaussian Processes for Machine Learning, Cambridge, MA: MIT Press, 2006.

[2] H. Zhang, "Asymptotics and Computation for Spatial Statistics," in Advances and Challenges in Space-time Modelling of Natural Events (Lecture Notes in Statistics, Vol. 207) (E. Porcu, J. M. Montero, and M. Schlather, eds.), New York: Springer, 2012 pp. 239–252. doi:10.1007/978-3-642-17086-7_10.

[3] C. G. Kaufman and B. A. Shaby, "The Role of the Range Parameter for Estimation and Prediction in Geostatistics," Biometrika, 100(2), 2013 pp. 473–484. doi:10.1093/biomet/ass079.

[4] N. Cressie, "Fitting Variogram Models by Weighted Least Squares," Journal of the International Association for Mathematical Geology, 17(5), 1985 pp. 563–586. doi:10.1007/BF01032109.

[5] H. Zhang and D. L. Zimmerman, "Hybrid Estimation of Semivariogram Parameters," Mathematical Geology, 39(2), 2007 pp. 247–260. doi:10.1007/s11004-006-9070-8.

[6] D. A. Girard, "Asymptotic Near-Efficiency of the 'Gibbs-Energy and Empirical-Variance' Estimating Functions for Fitting Matérn Models to a Dense (Noisy) Series."

[7] M. S. Phadke and S. M. Wu, "Modeling of Continuous Stochastic Processes from Discrete Observations with Application to Sunspots Data," Journal of the American Statistical Association, 69(346), 1974 pp. 325–329.

Feedback (field required)
Email (field required) Name
Occupation Organization
Note: Your message & contact information may be shared with the author of any specific Demonstration for which you give feedback.