Long-range dependence

1. Background

After finishing my master's thesis at NTNU in September 2016, I started my PhD at UiT The Arctic University of Norway, under the supervision of Sigrunn Holbek Sørbye (UiT), and co-supervision of Håvard Rue (KAUST) and Martin Rypdal (UiT). 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. In addition to my supervisors, I also collaborated with Hege-Beate Fredriksen (UiT) and Kristoffer Rypdal (UiT).

2. Introduction

2.1. Fractional Gaussian noise

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 definition 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 led 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, dendrochronology, 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 (Rypdal and Rypdal, 2016; Rybski et al., 2006; Lovejoy and Schertzer, 2013; Huybers and Curry, 2006; Franzke, 2010; Fredriksen and Rypdal, 2016).

2.2. Efficient Bayesian inference

It can be shown (see e.g. Theorem 2.2 in 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 $\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 is only conditionally dependent 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 posterior distributions are typically evaluated using time-consuming sampling-based methods such as Markov chain Monte Carlo (MCMC). The Markov property allows for these methods to be implemented more efficiently, but also allows for more efficient alternatives such as integrated nested Laplace approximations (INLA, Rue et al., 2009). INLA utilizes the sparse structure of the precision matrix to compute the posterior marginal distributions numerically, foregoing the need for sampling.

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, in theory, be possible to express an fGn process through a mixture of short memory processes, which would grant the Markov property and possibly reduce computational cost.

Unfortunately, as pointed out by Haldrup and Vera-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 high number of 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) processes. 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 small lags in the optimization procedure, which are more important for 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 INLA through its implementation in R (r-inla.org). Being able to incorporate the model into the R-INLA framework allows us to easily add other model components such as seasonal trends and random effects. 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 radiative forcing affects the GMST we can use a simple energy balance model, which describes the change in heat content of the climate system as a function of the radiative forcing $F$ and the heat loss to space. This can be expressed as

and $\Delta Q = C\Delta T$, where $\Delta Q$ is the change in the system's heat content, $C$ is the heat capacity, $\Delta T$ is the temperature anomaly and $\lambda$ is the feedback parameter. The forcing can be decomposed into a known component $F_{\text{known}}$, and an unknown, stochastic component $dB(t)$. By solving the energy balance model we find that the GMST can be expressed as the sum of two components. First, the response to known forcing

where $g(t-s) = e^{-(t-s)\lambda/C}$ is the Green's function. And second, an Ornstein-Uhlenbeck process $\bm \varepsilon$ which describes the response to the stochastic forcing component

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 expressed as

Here, $\bm K$ is a matrix representing the heat exchange between the boxes. As the number of boxes increases towards infinity, the Green's function 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$. Moreover, letting $F_0$ denote an unknown shift parameter in the known forcing, and $\sigma$ and $\sigma_f$ be scaling parameters for the stochastic and known forcing, respectively, we can express the GMST as

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, and a demonstration and short description is also available under the 'Software' section of this website. A more detailed description and tutorial of this package can be found in Appendix A in my PhD thesis (Myrvoll-Nilsen, 2020).

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-temporal dataset, GISS Surface Temperature Analysis version 4 (GISTEMP v4) (GISTEMP Team, 2018; Lenssen et al., 2019).

6. Relevant publications in this field

References