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

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 (the process variance being denoted by ), where is a modified Bessel function of the second kind (see the Details section in the help page for the BesselK function) and is still called the inverse-range parameter, being the "range" (or "support") of .
One important difference with the case is that the considered process, even regularly sampled (here observed at in the set ), 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 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 (the general expression for being ), is naturally defined as the square of the product of the estimates of and .
For all the simulated time series, the true variance 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 .
Let denote the candidate correlation matrix of , with inverse-range , that is, the matrix whose element in the row and column is . 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 , the so-called "Gibbs energy of associated with a given ".
The computation of 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 being obtained via FFT from the standard embedding of in a circulant matrix.
It is important to notice that becomes ill-conditioned for large and decreasing toward the lower-bound used here (as in the Demonstration for ); 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 .


  • [Snapshot]
  • [Snapshot]
  • [Snapshot]


[1] P. Whittle, "On Stationary Processes in the Plane," Biometrika, 41(3–4), 1954 pp. 434–449. doi:10.2307/2332724.
[2] H. Zhang, "Asymptotics and Computation for Spatial Statistics," in Advances and Challenges in Space-time Modelling of Natural Events (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] C. Janna, M. Ferronato, F. Sartoretto, and G. Gambolati, "FSAIPACK: A Software Package for High Performance FSAI Preconditioning," ACM Transactions on Mathematical Software, 41(2), forthcoming.
    • Share:

Embed Interactive Demonstration New!

Just copy and paste this snippet of JavaScript code into your website or blog to put the live Demonstration on your site. More details »

Files require Wolfram CDF Player or Mathematica.