Efficient Bayesian analysis of long-memory processes with applications in climate.

Summary

Long-memory processes are stochastic processes where correlations decay slowly over time. Unlike short-memory models such as autoregressive (AR) processes, long memory processes retain statistical dependence across long time horizons. Long memory properties have been observed in many environmental and geophysical time series, including global mean surface temperature records. However, Bayesian inference for long-memory processes is computationally expensive The main aim of my PhD thesis was to develop an efficient approximation of long-memory processes that allows for fast Bayesian inference, and to use this methodology to analyze global mean surface temperature records.

1. Introduction

Long-range dependence (LRD) or long memory is a property relating to the rate of decay of statistical dependence of two observations as the distance between them increases.A common formulation for long-memory behavior is given by a hyperbolic (as opposed to an exponential) decay of the autocovariance function.

Here, $H$ measures the degree of long-range dependence and is commonly referred to as the Hurst exponent. The exponent is named after H.E. Hurst, who, while studying the storage requirements of reservoirs on the Nile River, observed that the rescaled range statistic of hydrological records followed a power-law (Hurst, 1951), which could not be explained using contemporary short memory models. Hurst's results lead to increased interest in the topic, and both hydrologists and mathematicians tried to find theoretical explanations of the so-called Hurst phenomenon.

A major breakthrough came with the introduction of fractional Brownian motion and its increment process, fractional Gaussian noise (fGn), by Mandelbrot and van Ness (1968)This is a process with an autocovariance function that decays hyperbolically

The fractional Gaussian noise has since become a widely used model for long-memory behaviour, with applications in fields such as hydrology, dendrochronologyn, finance, among others. For my thesis I focused on applications in climate science, specifically on global mean surface temperature records, where long memory has been observed. (CITECITECITE)However, long memory processes are difficult to work with due to

2. Efficient Bayesian inference

It can be shown (Theorem XXX, Rue and Held, 2005) that if two variables $x_i$ and $x_j$, in a Gaussian random field $\bm x$, are conditionally independent given the other variables in the field, $\bm{x}_{-ij}$, the corresponding term in the precision matrix is zero, i.e. $Q_{ij}=0$. Thus, if a $\bm x$ has a high degree of conditional independence we say that it is a Gaussian Markov random field (GMRF) and the precision matrix (inverse covariance matrix) will then be sparse. For example, since each node of an AR(1) process only depends on its immediate predecessor it is a GMRF with a tridiagonal precision matrix.

The sparse precision matrix associated with the Markov property allows for many efficient algorithms for key operations to be utilized. This includes computing the log-likelihood (Rue, 2001), performing Cholesky factorization and producing samples. This is particularly important for Bayesian inference, where the the posterior distributions are typically evaluated using time-consuming sampling-based methods such as Markov chain Monte Carlo (MCMC). The Markov property allows for the use of more efficient algorithms such as integrated nested Laplace approximations (INLA, Rue et al., 2009), which utilizes the sparse structure of the precision matrix to compute the posterior marginal distributions numerically.

Unlike Markov processes like the AR(1), long-memory processes exhibit persistent correlations over extended time periods, where each variable is conditionally dependent on all other variables. This means that the precision matrix is dense, and key inference algorithms have a computational cost that scales cubically with the number of observations, making them impractical for large datasets. In simple cases it is possible to exploit the Toeplitz structure of the covariance matrix to achieve $\mathcal{O}(n^2)$ cost, but in many situations, such as when we have inhomogeneous variance, the Toeplitz structure breaks and the cost increases to $\mathcal{O}(n^3)$. This motivates the need for an approximate model that can capture the long range dependence of the fGn process while retaining the Markov property, allowing for efficient inference algorithms to be utilized.

Granger (1980) showed that if the lag-one autocorrelation parameter $\phi$ of an AR(1) process is sampled from a beta distribution the aggregation of AR(1) processes would exhibit long-range dependence according to a fractionally integrated process. It should therefore be possible to express an fGn process through a mixture of short memory processes, which would grant the Markov property and possibly reducing computational cost.

Unfortunately, as pointed out by Haldrup and Valdés (2017), the number of AR(1) processes (with $\phi$ drawn from a beta distribution) required to exhibit long-range dependence is so large that any gains in the Markov property would be lost to the increase in variables associated with the AR(1) processes. Making a sampling-based approach impractical for inference.

3. AR(1) mixture approximation

In Sørbye et al. (2019) we propose a numerical approach to approximating long-range dependence by aggregating AR(1) process. Specifically, we approximate the fGn process $\bm z$ by:

where $\bm x_i$ denote independent AR(1) processes with unit variance and lag-one correlation $\phi_i$. To avoid a singular precision matrix for the joint GMRF $(\tilde{\bm{z}},\bm x_1,...,\bm x_m)$ we add a small noise term $\tau\sim \mathcal{N}(0,\exp(-15))$. Instead of drawing the lag-one parameters from a random distribution we instead used numerical optimization to find the weights $\bm w$ and lag-one parameters $\bm \phi$ which results in the autocorrelation function (acf), $\gamma_{\tilde{\bm z}}$, which (up to a selected lag $k_{\text{max}}$) best approximates the target acf, $\gamma_H$, of an fGn with some Hurst exponent $H$,

Formally, the optimization problem can be expressed as:

and is repeated for all $H$ in a grid covering the persistent range $H \in (0.5, 1)$. The weighting term $1/k$ is used to give more importance to the fit of the acf at small lags, which improves estimation accuracy. The optimization yields a mapping from $H$ to a corresponding set of weights $\bm w$ and lag-one autocorrelation parameters $\bm \phi$, which in turn gives us an approximate fGn process expressed as a mixture of $m$ AR(1) processes. We find that only using $m=4$ AR(1) processes in the aggregation provides excellent approximation of the autocorrelation function. We also find good estimation accuracy even for time series exceeding $k_{\text{max}}=1000$. See Sørbye et al. (2019) or Myrvoll-Nilsen (2020) for more details.

To fully exploit the computational advantages of a GMRF approximation of an fGn we formulate a latent Gaussian model, and perform full Bayesian inference using integrated nested Laplace approximations (INLA, Rue et al., 2009). Being able to incorporate the model into the R-INLA framework allows us to easily add other model components such as linear or seasonal trends. The model has since been incorporated into R-INLA as a standard model.The approximate fGn model laid the foundation for my PhD thesis, and was subsequently applied to a number of real-world problems in climate science.

4. Incorporating radiative forcing

Reliable quantification of the global mean surface temperature (GMST) response to radiative forcing is essential for assessing the risk of dangerous anthropogenic climate change. For time scales ranging from months to centuries this response can be expressed as scale-invariant (Rypdal and Rypdal, 2016; Rybski et al., 2006; Lovejoy and Schertzer, 2013; Huybers and Curry, 2006; Franzke, 2010; Fredriksen and Rypdal, 2016), implying long memory properties.

To understand how

where $\Delta Q = C\Delta T$. The forcing can be decomposed into a known component $F_{\text{known}}$, and an unknown, stochastic component $dB(t)$. Solving the energy balance model yields solution:

Where $\Delta T_{\text{known}}$ is the GMST reponse to the known forcing component and $\bm \varepsilon$ is an Ornstein-Uhlenbeck process. When discretized, the stochastic component results in an AR(1) process. This model is, however, not consistent with observations (Rypdal and Rypdal, 2014) due to the slow climate response associated with the energy exchange with the deep ocean. This simple energy balance model can be extended to an $m$-box model where heat is exchanged between $m$ 'boxes', representing different layers of the deep ocean. Mathematically, the $m$-box model can be express as

Here, dawdawdawdawdawdawdAs the number of boxes increase towards infinity, $g$ can be expressed as a scale invariant response function (Fredriksen and Rypdal, 2017)

which means that the response to the stochastic component can be expressed as a fractional Gaussian noise process with Hurst exponent $H$.

Since none of the default R-INLA models support this specific formulation we need to construct the latent model component using the custom modeling framework of R-INLA called rgeneric. To make the temperature response model easily accessible to other researchers we have developed a user-friendly R-package called INLA.climate, which takes care of the technical implementation of the model and allows users to easily fit the model to their data using R-INLA. This package is available at my GitHub repository: INLA.climate. A description and tutorial of this package can be found in Appendix A in my PhD thesis (Myrvoll-Nilsen, 2020). A demonstration and short description is also available under the 'Software' section of this website.

5. Climate applications

There are several benefits of having the model specified through the R-INLA framework. For example, any missing values in the GMST will be imputed automatically using the posterior predictive distribution (but the known forcing has to be provided as input for all points). This also means that we can very easily use the model for prediction. In Myrvoll-Nilsen et al. (2020) we predict the GMST under the four Representation Concentration Pathways (RCP) provided by the Coupled Model Intercomparison Project (CMIP) for the period 2006-2100, and find that the model provides good predictive performance even for long-term predictions up to 100 years into the future.

In other applications we use the model to estimate both the transient climate sensitivity (Myrvoll-Nilsen et al., 2020) and the equilibrium climate sensitivity (Rypdal et al., 2018), which are key metrics for understanding the risk of dangerous anthropogenic climate change. In Myrvoll-Nilsen et al. (2019) we exploit the reduced computational cost in order to fit the model to all local temperature records in the spatio-tepmoral dataset, GISS Surface Temperature Analysis version 4 (GISTEMP v) (GISTEMP Team, 2018; Lenssen et al., 2019).

6. Relevant publications in this field

• Sørbye, S. H., Myrvoll-Nilsen, E., & Rue, H. (2019). A computationally efficient approximation of the fractional Gaussian noise process. Statistics and Computing, 29(6), 1221–1233.
• Myrvoll-Nilsen, E., ...
• Myrvoll-Nilsen, E. (2020). Efficient Bayesian inference for long-memory processes applied to climate. PhD thesis, UiT The Arctic University of Norway.

References