Efficient Bayesian inference for natural time series using ARFIMA processes

Abstract

Many geophysical quantities, such as atmospheric temperature, water levels in rivers, and wind speeds, have shown evidence of long memory (LM). LM implies that these quantities experience non-trivial temporal memory, which potentially not only enhances their predictability, but also hampers the detection of externally forced trends. Thus, it is important to reliably identify whether or not a system exhibits LM. In this paper we present a modern and systematic approach to the inference of LM. We use the flexible autoregressive fractional integrated moving average (ARFIMA) model, which is widely used in time series analysis, and of increasing interest in climate science. Unlike most previous work on the inference of LM, which is frequentist in nature, we provide a systematic treatment of Bayesian inference. In particular, we provide a new approximate likelihood for efficient parameter inference, and show how nuisance parameters (e.g., short-memory effects) can be integrated over in order to focus on long-memory parameters and hypothesis testing more directly. We illustrate our new methodology on the Nile water level data and the central England temperature (CET) time series, with favorable comparison to the standard estimators. For CET we also extend our method to seasonal long memory.

Introduction

Many natural processes are sufficiently complex that a stochastic model is essential, or at the very least an efficient description . Such a process will be specified by several properties, of which a particularly important one is the degree of memory in a time series, often expressed through a characteristic autocorrelation time over which fluctuations will decay in magnitude. In this paper, however, we are concerned with specific types of stochastic processes that are capable of possessing long memory (LM) . Long memory is the notion of there being correlation between the present and all points in the past. A standard definition is that a (finite variance, stationary) process has long memory if its autocorrelation function (ACF) has power-law decay: $\mathit{\rho }(\cdot )$ such that $\mathit{\rho }\left(k\right)$ $\sim$ $c_{\mathit{\rho }}{k}^{\mathrm{2}d-\mathrm{1}}$ as $k$ $\to$ $\mathrm{\infty }$, for some non-zero constant $c_{\mathit{\rho }}$, and where 0 $<$ $d$ $<$ $\frac{\mathrm{1}}{\mathrm{2}}$. The parameter $d$ is the memory parameter; if $d$ $=$ 0 the process does not exhibit long memory, while if $-\frac{\mathrm{1}}{\mathrm{2}}$ $<$ $d$ $<$ 0 the process is said to be anti-persistent.

The asymptotic power-law form of the ACF corresponds to an absence of a characteristic decay timescale, in striking contrast to many standard (stationary) stochastic processes where the effect of each data point decays so fast that it rapidly becomes indistinguishable from noise. An example of the latter is the exponential ACF, where the e-folding timescale sets a characteristic correlation time. The study of processes that do possess long memory is important because they exhibit unusual properties, because many familiar mathematical results fail to hold, and because of the numerous examples of data sets where LM is seen.

The study of long memory originated in the 1950s in the field of hydrology, where studies of the levels of the Nile demonstrated anomalously fast growth of the rescaled range of the time series. After protracted debates

For a detailed exposition of this period of mathematical history, see .

about whether this was a transient (finite time) effect, the mathematical pioneer Benoît B. Mandelbrot showed that if one retained the assumption of stationarity, novel mathematics would then be essential to sufficiently explain the Hurst effect. In doing so he rigorously defined the concept of long memory .

Most research into long memory and its properties has been based on classical statistical methods, spanning parametric, semi-parametric, and non-parametric modeling seefor a review. Very few Bayesian methods have been studied, most probably due to computational difficulties. The earliest works are parametric and include , , and . If computational challenges could be mitigated, the Bayesian paradigm would offer advantages over classical methods including flexibility in specification of priors (i.e., physical expertise could be used to elicit an informative prior). It would offer the ability to marginalize out aspects of a model apparatus and data, such as short-memory or seasonal effects and missing observations, so that statements about long-memory effects can be made unconditionally.

Towards easing the computational burden, we focus on the autoregressive fractional integrated moving average (ARFIMA) class of processes as the basis of developing a systematic and unifying Bayesian framework for modeling a variety of common time series phenomena, with particular emphasis on (marginally) detecting potential long-memory effects (i.e., averaging over short-memory and seasonal effects). ARFIMA has become very popular in statistics and econometrics because it is generalizable and its connection to the autoregressive moving average (ARMA) family and to fractional Gaussian noise is relatively transparent. A key property of ARFIMA is its ability to simultaneously yet separately model long and short memory.

Here we present a Bayesian framework for the efficient and systematic estimation of the ARFIMA parameters. We provide a new approximate likelihood for ARFIMA processes that can be computed quickly for repeated evaluation on large time series, and which underpins an efficient Markov chain Monte Carlo (MCMC) scheme for Bayesian inference. Our sampling scheme can be best described as a modernization of a blocked MCMC scheme proposed by – adapting it to the approximate likelihood and extending it to handle a richer form of (known) short-memory effects. We then further extend the analysis to the case where the short-memory form is unknown, which requires trans-dimensional MCMC, in which the model order (the $p$ and $q$ parameters in the ARFIMA model) varies and, thus, so does the dimension of the problem. This aspect is similar to the work of , who considered the simpler autoregressive-integrated moving average (ARIMA) model class, and to , who worked with a non-parametric long-memory process. Our contribution has aspects in common with , who presented a more limited method focused on model selection rather than averaging. The advantage of averaging is that the unknown form of short-memory effects can be integrated out, focusing on long memory without conditioning on nuisance parameters.

The aim of this paper is to introduce an efficient Bayesian algorithm for the inference of the parameters of the ARFIMA($p$,$d$,$q$) model, with particular emphasis on the LM parameter $d$. Our Bayesian inference algorithm has been designed in a flexible fashion so that, for instance, the innovations can come from a wide class of different distributions, e.g., $\mathit{\alpha }$ stable or $t$ distribution (to be published in a companion paper). The remainder of the paper is organized as follows. Section  discusses the important numerical calculation of likelihoods, representing a hybrid between earlier classical statistical methods and our new contributions towards a full-Bayesian approach. Section  describes our proposed Bayesian framework and methodology in detail, focusing on long memory only. Then, in Sect. , we consider extensions for additional short memory and the computational techniques required to integrate them out. Empirical illustration and comparison of all methods is provided in Sect.  via synthetic and real data including the Nile water level data and the central England temperature (CET) time series, with favorable comparison to the standard estimators. In the case of the Nile data, we find strong evidence for long memory. The CET analysis requires a slight extension to handle seasonal long memory, and we find that the situation here is more nuanced in terms of evidence for long memory. The paper concludes with a discussion in Sect.  focused on the potential for further extension.

Likelihood evaluation for Bayesian inference

ARFIMA model

We provide here a brief review of the ARFIMA model. More details are given in Appendix .

An ARFIMA model is given by $$\mathrm{\Phi }\left(\mathcal{B}\right)(\mathrm{1}-\mathcal{B}{)}^d{X}_t=\mathrm{\Theta }(\mathcal{B}){\mathit{\epsilon }}_t.$$ We define the backshift operator $\mathcal{B}$, where $\mathcal{B}{X}_t$ $=$ $X_{t-\mathrm{1}}$, and powers of $\mathcal{B}$ are defined iteratively: ${\mathcal{B}}^k{X}_t$ $=$ ${\mathcal{B}}^{k-\mathrm{1}}{X}_{t-\mathrm{1}}$ $=$ $\mathrm{\cdots }$ $=$ $X_{t-k}$. $\mathrm{\Phi }$ is the autoregressive component and $\mathrm{\Theta }$ is the moving average component and constitute the short-memory components of the ARFIMA model. These are defined in more detail in Appendix  and in .

Likelihood function

For now, we restrict our attention to a Bayesian analysis of an ARFIMA(0,$d$,0) process, having no short-ranged ARMA components ($p$ $=$ $q$ $=$ 0), placing emphasis squarely on the memory parameter $d$. As we explain in our Appendix, the resulting process is identical to a fractionally integrated processes with memory parameter $d$.

Here we develop an efficient and new scheme for evaluating the (log) likelihood, via approximation. Throughout, the reader should suppose that we have observed the vector $\mathit{x}$ $=$ ($x_{\mathrm{1}}$, …, $x_n{)}^{\top }$ as a realization of a stationary, causal and invertible ARFIMA(0,$d$,0) process $\mathit{\left\{}{X}_t\mathit{\right\}}$ with mean $\mathit{\mu }$ $\in$ $\mathbb{R}$. The innovations will be assumed to be independent, and taken from a zero-mean location-scale probability density $f(\cdot$; 0, $\mathit{\sigma }$, $\mathit{\lambda }$), which means the density can be written as $f(x$; $\mathit{\delta }$, $\mathit{\sigma }$, $\mathit{\lambda }$) $\equiv$ $\frac{\mathrm{1}}{\mathit{\sigma }}f(\frac{x-\mathit{\delta }}{\mathit{\sigma }}$; 0, 1, $\mathit{\lambda }$). The parameters $\mathit{\delta }$ and $\mathit{\sigma }$ are called the location and scale parameters, respectively. The $m$-dimensional $\mathit{\lambda }$ is a shape parameter (if it exists, i.e., $m$ $>$ 0). A common example is the Gaussian $\mathcal{N}(\mathit{\mu }$, ${\mathit{\sigma }}^{\mathrm{2}}$), where $\mathit{\delta }$ $\equiv$ $\mathit{\mu }$ and there is no $\mathit{\lambda }$. We classify the four parameters $\mathit{\mu }$, $\mathit{\sigma }$, $\mathit{\lambda }$, and $d$ into three distinct classes: (1) the mean of process, $\mathit{\mu }$; (2) innovation distribution parameters, $\mathit{\upsilon }$ $=$ ($\mathit{\sigma }$, $\mathit{\lambda }$); and (3) memory structure, $d$. Together, $\mathit{\psi }$ $=$ ($\mathit{\mu }$, $\mathit{\upsilon }$, $\mathit{\omega }$), where $\mathit{\omega }$ will later encompass the short-range parameters $p$ and $q$.

Our proposed likelihood approximation uses a truncated autoregressive model (AR) ($\mathrm{\infty }$) approximation (cf. ). We first re-write the AR($\mathrm{\infty }$) approximation of ARFIMA(0,$d$,0) to incorporate the unknown parameter $\mathit{\mu }$, and drop the ($d$) superscript for convenience: $X_t$ $-$ $\mathit{\mu }$ $=$ ${\mathit{\epsilon }}_t$ $-$ $\sum _{k=\mathrm{1}}^{\mathrm{\infty }}{\mathit{\pi }}_k(X_{t-k}$ $-$ $\mathit{\mu }$). Then we truncate this AR($\mathrm{\infty }$) representation to obtain an AR($P$) one, with $P$ large enough to retain low frequency effects, e.g., $P$ $=$ $n$. We denote ${\mathrm{\Pi }}_P$ $=$ $\sum _{k=\mathrm{0}}^P{\mathit{\pi }}_k$ and, with ${\mathit{\pi }}_{\mathrm{0}}$ $=$ 1, rearrange terms to obtain the following modified model: $$X_t={\mathit{\epsilon }}_t+{\mathrm{\Pi }}_P\mathit{\mu }-\sum _{k=\mathrm{1}}^P{\mathit{\pi }}_k{X}_{t-k}.$$

It is now possible to write down a conditional likelihood. For convenience the notation ${\mathit{x}}_k$ $=$ ($x_{\mathrm{1}}$, …, $x_k{)}^{\top }$ for $k$ $=$ 1, …, $n$ will be used (and ${\mathit{x}}_{\mathrm{0}}$ is interpreted as appropriate where necessary). Denote the unobserved $\mathit{P}$ vector of random variables ($x_{\mathrm{1}-P}$, …, $x_{-\mathrm{1}}$, $x_{\mathrm{0}}{)}^{\top }$ by ${\mathit{x}}_{\text{A}}$ (in the Bayesian context these will be auxiliary, hence “$A$”). Consider the likelihood $L\left(\mathit{x}\right|\mathit{\psi })$ as a joint density, which can be factorized as a product of conditionals. Writing $g_t\left(x_t\right|{\mathit{x}}_{t-\mathrm{1}}$, $\mathit{\psi })$ for the density of $X_t$ conditional on ${\mathit{x}}_{t-\mathrm{1}}$, we obtain $L\left(\mathit{x}\right|\mathit{\psi })$ $=$ ${\prod }_{t=\mathrm{1}}^n{g}_t\left(x_t\right|{\mathit{x}}_{t-\mathrm{1}}$, $\mathit{\psi })$.

This is still of little use because the $g_t$ may have a complicated form. However, by further conditioning on ${\mathit{x}}_{\text{A}}$, and writing $h_t\left(x_t\right|{\mathit{x}}_{\text{A}}$, ${\mathit{x}}_{t-\mathrm{1}}$, $\mathit{\psi })$ for the density of $X_t$ conditional on ${\mathit{x}}_{t-\mathrm{1}}$ and ${\mathit{x}}_{\text{A}}$, we obtain $L\left(\mathit{x}\right|\mathit{\psi }$, ${\mathit{x}}_{\text{A}}$) $=$ ${\prod }_{t=\mathrm{1}}^n{h}_t\left(x_t\right|{\mathit{x}}_{\text{A}}$, ${\mathit{x}}_{t-\mathrm{1}}$, $\mathit{\psi }$). Returning to Eq. () observe that, conditional on both the observed and unobserved past values, $X_t$ is simply distributed according to the innovations' density $f$ with a suitable change in location: $X_t|{\mathit{x}}_{t-\mathrm{1}}$, ${\mathit{x}}_{\text{A}}$ $\sim$ $f(\cdot$; [${\mathrm{\Pi }}_P\mathit{\mu }$ $-$ $\sum _{k=\mathrm{1}}^P{\mathit{\pi }}_k{x}_{t-k}$], $\mathit{\sigma }$, $\mathit{\lambda }$). Then using location-scale representation: $$\begin{array}{cc}{h}_t\left(x_t,|,{\mathit{x}}_{\text{A}},,,{\mathit{x}}_{t-\mathrm{1}},,,\mathit{\psi }\right)& \approx f\left(x_t,;,\left[{\mathrm{\Pi }}_P,\mathit{\mu },-,\sum _{k=\mathrm{1}}^P,{\mathit{\pi }}_k,x_{t-k}\right],,,\mathit{\sigma },,,\mathit{\lambda }\right)\\ & \equiv {\displaystyle \frac{\mathrm{1}}{\mathit{\sigma }}}f\left({\displaystyle \frac{c_t-{\mathrm{\Pi }}_P\mathit{\mu }}{\mathit{\sigma }}},;,\mathrm{0},,,\mathrm{1},,,\mathit{\lambda }\right),\\ & & \text{where}{c}_t=\sum _{k=\mathrm{0}}^P{\mathit{\pi }}_k{x}_{t-k},t=\mathrm{1},\mathrm{\dots },n.\end{array}$$ Therefore, $L\left(\mathit{x}\right|\mathit{\psi }$, ${\mathit{x}}_{\text{A}}$) $\approx$ ${\mathit{\sigma }}^{-n}{\prod }_{t=\mathrm{1}}^{n}f\left(\frac{c_t-{\mathrm{\Pi }}_P\mathit{\mu }}{\mathit{\sigma }},;,,\mathit{\lambda }\right)$, or equivalently $$\mathrm{\ell }\left(\mathit{x},|,\mathit{\psi },,,{\mathit{x}}_{\text{A}}\right)\approx -n\log \mathit{\sigma }+\sum _{t=\mathrm{1}}^n\log \left\{f,\left({\displaystyle \frac{c_t-{\mathrm{\Pi }}_P\mathit{\mu }}{\mathit{\sigma }}},;,\mathit{\lambda }\right)\right\}.$$

Evaluating this expression efficiently depends upon efficient calculation of $\mathbf{c}$ $=$ ($c_{\mathrm{1}}$, …, $c_n{)}^t$ and $\log f(\cdot )$. From Eq. (), $\mathbf{c}$ is a convolution of the augmented data, $(\mathit{x}$, ${\mathit{x}}_{\text{A}})$, and coefficients depending on $d$, which can be evaluated quickly in the R language for statistical computing via convolve via fast Fourier transform (FFT). Consequently, evaluation of the conditional on ${\mathit{x}}_A$ likelihood in the Gaussian case costs only $\mathcal{O}\left(n\log n\right)$ – a clear improvement over the exact method. Obtaining the unconditional likelihood requires marginalization over ${\mathit{x}}_{\text{A}}$, which is analytically infeasible. However, this conditional form will suffice in the context of our Bayesian inferential scheme, presented below.

A Bayesian approach to long-memory inference

We are now ready to consider Bayesian inference for ARFIMA(0,$d$,0) processes. Our method can be succinctly described as a modernization of the blocked MCMC method of . Isolating parameters by blocking provides significant scope for modularization, which helps to accommodate our extensions for short memory. Pairing with efficient likelihood evaluations allows much longer time series to be entertained than ever before. Our description begins with the appropriate specification of priors, which are more general than previous choices, yet still encourages tractable inference. We then provide the relevant updating calculations for all parameters, including those for auxiliary parameters ${\mathit{x}}_{\text{A}}$.

We follow earlier work and assume a priori independence for components of $\mathit{\psi }$. Each component will leverage familiar prior forms with diffuse versions as limiting cases. Specifically, we use a diffuse Gaussian prior on $\mathit{\mu }$: $\mathit{\mu }$ $\sim$ $\mathcal{N}({\mathit{\mu }}_{\mathrm{0}}$, ${\mathit{\sigma }}_{\mathrm{0}}^{\mathrm{2}}$), with ${\mathit{\sigma }}_{\mathrm{0}}$ large. The improper flat prior is obtained as the limiting distribution when ${\mathit{\sigma }}_{\mathrm{0}}$ $\to$ $\mathrm{\infty }$: $p_{\mathit{\mu }}\left(\mathit{\mu }\right)$ $\propto$ 1. We place a gamma prior on the precision $\mathit{\tau }$ $=$ ${\mathit{\sigma }}^{-\mathrm{2}}$ implying a root-inverse gamma distribution $\mathcal{R}({\mathit{\alpha }}_{\mathrm{0}}$, ${\mathit{\beta }}_{\mathrm{0}})$ for $\mathit{\sigma }$, with density $f\left(\mathit{\sigma }\right)$ $=$ $\frac{\mathrm{2}}{\mathrm{\Gamma }\left(\mathit{\alpha }\right)}{{\mathit{\beta }}_{\mathrm{0}}}^{{\mathit{\alpha }}_{\mathrm{0}}}{\mathit{\sigma }}^{-(\mathrm{2}{\mathit{\alpha }}_{\mathrm{0}}+\mathrm{1})}\exp \left(-,\frac{{\mathit{\beta }}_{\mathrm{0}}}{y^{\mathrm{2}}}\right)$, $\mathit{\sigma }$ $>$ 0. A diffuse/improper prior is obtained as the limiting distribution when ${\mathit{\alpha }}_{\mathrm{0}}$, ${\mathit{\beta }}_{\mathrm{0}}$ $\to$ 0: $p_{\mathit{\sigma }}\left(\mathit{\sigma }\right)$ $\propto$ ${\mathit{\sigma }}^{-\mathrm{1}}$, which, in the asymptotic limit, is equivalent to a log uniform prior. Finally, we specify $d$ $\sim$ $\mathcal{U}(-\frac{\mathrm{1}}{\mathrm{2}}$, $\frac{\mathrm{1}}{\mathrm{2}}$).

Updating $\mathit{\mu }$: following , we use a symmetric random walk (RW) Metropolis–Hastings (MH) update with proposals ${\mathit{\xi }}_{\mathit{\mu }}$ $\sim$ $\mathcal{N}(\mathit{\mu }$, ${\mathit{\sigma }}_{\mathit{\mu }}^{\mathrm{2}})$, for some ${\mathit{\sigma }}_{\mathit{\mu }}^{\mathrm{2}}$. The acceptance ratio is $$\begin{array}{cc}{A}_{\mathit{\mu }}& \left(\mathit{\mu },,,{\mathit{\xi }}_{\mathit{\mu }}\right)=\sum _{t=\mathrm{1}}^n\log \left\{,f,\left({\displaystyle \frac{c_t-{\mathrm{\Pi }}_P{\mathit{\xi }}_{\mathit{\mu }}}{\mathit{\sigma }}},;,\mathit{\lambda }\right),\right\}\\ & & -\sum _{t=\mathrm{1}}^n\log \left\{,f,\left({\displaystyle \frac{c_t-{\mathrm{\Pi }}_P\mathit{\mu }}{\mathit{\sigma }}},;,\mathit{\lambda }\right),\right\}+\log \left[{\displaystyle \frac{p_{\mathit{\mu }}\left({\mathit{\xi }}_{\mathit{\mu }}\right)}{p_{\mathit{\mu }}\left(\mathit{\mu }\right)}}\right]\end{array}$$ under the approximate likelihood.

Updating $\mathit{\sigma }$: we diverge from here, who suggest independent MH with moment-matched inverse gamma proposals, finding poor performance under poor moment estimates. We instead prefer a RW MH approach, which we conduct in log space since the domain is ${\mathbb{R}}^{+}$. Specifically, we set: $\log {\mathit{\xi }}_{\mathit{\sigma }}$ $=$ $\log \mathit{\sigma }$ $+$ $\mathit{\upsilon }$, where $\mathit{\upsilon }$ $\sim$ $\mathcal{N}$(0, ${\mathit{\sigma }}_{\mathit{\sigma }}^{\mathrm{2}}$) for some ${\mathit{\sigma }}_{\mathit{\sigma }}^{\mathrm{2}}$. ${\mathit{\xi }}_{\mathit{\sigma }}|\mathit{\sigma }$ is log-normal and we obtain $\frac{q(\mathit{\sigma };{\mathit{\xi }}_{\mathit{\sigma }})}{q({\mathit{\xi }}_{\mathit{\sigma }};\mathit{\sigma })}$ $=$ $\frac{{\mathit{\xi }}_{\mathit{\sigma }}}{\mathit{\sigma }}$. Recalling Eq. (), the MH acceptance ratio under the approximate likelihood is $$\begin{array}{cc}{A}_{\mathit{\sigma }}\left(\mathit{\sigma },,,{\mathit{\xi }}_{\mathit{\sigma }}\right)& =\sum _{t=\mathrm{1}}^n\log \left\{f,\left({\displaystyle \frac{c_t-{\mathrm{\Pi }}_P\mathit{\mu }}{{\mathit{\xi }}_{\mathit{\sigma }}}},;,\mathit{\lambda }\right)\right\}\\ & -\sum _{t=\mathrm{1}}^n\log \left\{f,\left({\displaystyle \frac{c_t-{\mathrm{\Pi }}_P\mathit{\mu }}{\mathit{\sigma }}},;,\mathit{\lambda }\right)\right\}\\ & +\log \left[{\displaystyle \frac{p_{\mathit{\sigma }}\left({\mathit{\xi }}_{\mathit{\sigma }}\right)}{p_{\mathit{\sigma }}\left(\mathit{\sigma }\right)}}\right]+(n-\mathrm{1})\log \left[{\displaystyle \frac{\mathit{\sigma }}{{\mathit{\xi }}_{\mathit{\sigma }}}}\right].\end{array}$$

The MH algorithm, applied alternately in a Metropolis-within-Gibbs fashion to the parameters $\mathit{\mu }$ and $\mathit{\sigma }$, works well. However, actual Gibbs sampling is an efficient alternative in this two-parameter case (i.e., for known $d$; see ).

Update of $d$: updating the memory parameter $d$ is far less straightforward than either $\mathit{\mu }$ or $\mathit{\sigma }$. Regardless of the innovations' distribution, the conditional posterior ${\mathit{\pi }}_{d|{\mathit{\psi }}_{-d}}\left(d\right|{\mathit{\psi }}_{-d},\mathit{x})$ is not amenable to Gibbs sampling. We use RW proposals from truncated Gaussian ${\mathit{\xi }}_d$ $\sim$ ${\mathcal{N}}^{(a,b)}(d$, ${\mathit{\sigma }}_d^{\mathrm{2}})$, with density $$\begin{array}{cc}f(x;\mathit{\mu },\mathit{\sigma },a,b)=& {\displaystyle \frac{\mathrm{1}}{\mathit{\sigma }}}{\displaystyle \frac{{\mathit{\varphi }}^{\left(\mathcal{N}\right)}\left[\right(x-\mathit{\mu })/\mathit{\sigma }]}{{\mathrm{\Phi }}^{\left(\mathcal{N}\right)}\left[\right(b-\mathit{\mu })/\mathit{\sigma }]-{\mathrm{\Phi }}^{\left(\mathcal{N}\right)}\left[\right(a-\mathit{\mu })/\mathit{\sigma }],}}\\ & & a<x<b,\end{array}$$ where ${\mathrm{\Phi }}^{\left(\mathcal{N}\right)}$ and ${\mathit{\varphi }}^{\left(\mathcal{N}\right)}$ are the standard normal cumulative density function (CDF) and probability density function (PDF), respectively. In particular, we use ${\mathit{\xi }}_d|d$ $\sim$ ${\mathcal{N}}^{(-\mathrm{1}/\mathrm{2},\mathrm{1}/\mathrm{2})}(d$, ${\mathit{\sigma }}_d^{\mathrm{2}})$ via rejection sampling from $\mathcal{N}(d$, ${\mathit{\sigma }}_d^{\mathrm{2}})$ until ${\mathit{\xi }}_d$,$\in$ ($-\frac{\mathrm{1}}{\mathrm{2}}$, $\frac{\mathrm{1}}{\mathrm{2}}$). Although this may seem inefficient, it is perfectly acceptable; for example, if ${\mathit{\sigma }}_d$ $=$ 0.5 the expected number of required variates is still less than 2, regardless of $d$. More refined methods of directly sampling from truncated normal distributions exist – see for example – but we find little added benefit in our context.

A useful cancellation in $q(d$; ${\mathit{\xi }}_d)/q({\mathit{\xi }}_d$; $d)$ obtained from Eq. () yields $$\begin{array}{cc}{A}_d& =\mathrm{\ell }\left(\mathit{x},|,{\mathit{\xi }}_d,,,{\mathit{\psi }}_{-d}\right)-\mathrm{\ell }\left(\mathit{x},|,d,,,{\mathit{\psi }}_{-d}\right)+\log \left[{\displaystyle \frac{p_d\left({\mathit{\xi }}_d\right)}{p_d\left(d\right)}}\right]\\ & +\log \left\{,{\displaystyle \frac{{\mathrm{\Phi }}^{\left(\mathcal{N}\right)}\left[\left(\frac{\mathrm{1}}{\mathrm{2}},-,d\right),/,{\mathit{\sigma }}_d\right]-{\mathrm{\Phi }}^{\left(\mathcal{N}\right)}\left[\left(-,\frac{\mathrm{1}}{\mathrm{2}},-,d\right),/,{\mathit{\sigma }}_d\right]}{{\mathrm{\Phi }}^{\left(\mathcal{N}\right)}\left[\left(\frac{\mathrm{1}}{\mathrm{2}},-,{\mathit{\xi }}_d\right),/,{\mathit{\sigma }}_d\right]-{\mathrm{\Phi }}^{\left(\mathcal{N}\right)}\left[\left(-,\frac{\mathrm{1}}{\mathrm{2}},-,{\mathit{\xi }}_d\right),/,{\mathit{\sigma }}_d\right]}},\right\}.\end{array}$$ Denote ${\mathit{\xi }}_{c_t}$ $=$ $\sum _{k=\mathrm{0}}^P{\mathit{\xi }}_{{\mathit{\pi }}_k}{x}_{t-k}$ for $t$ $=$ 1, …, $n$, where $\mathit{\left\{}{\mathit{\xi }}_{{\mathit{\pi }}_k}\mathit{\right\}}$ are the proposed coefficients $\mathit{\left\{}{\mathit{\pi }}_k^{\left({\mathit{\xi }}_d\right)}\mathit{\right\}}$; ${\mathit{\pi }}_k^{\left(d\right)}$ $=$ $\frac{\mathrm{1}}{\mathrm{\Gamma }(k+\mathrm{1})}\frac{\mathrm{\Gamma }(k-d)}{\mathrm{\Gamma }(-d)}$; furthermore, ${\mathit{\xi }}_{{\mathrm{\Pi }}_P}$ $-$ $\sum _{k=\mathrm{0}}^P{\mathit{\xi }}_{{\mathit{\pi }}_k}$. Then in the approximate case $$\begin{array}{cc}{A}_d& =\sum _{t=\mathrm{1}}^n\log \left\{f,\left({\displaystyle \frac{{\mathit{\xi }}_{c_t}-{\mathit{\xi }}_{{\mathrm{\Pi }}_P}\mathit{\mu }}{\mathit{\sigma }}},;,\mathit{\lambda }\right)\right\}\\ & -\sum _{t=\mathrm{1}}^n\log \left\{f,\left({\displaystyle \frac{c_t-{\mathrm{\Pi }}_P\mathit{\mu }}{\mathit{\sigma }}},;,\mathit{\lambda }\right)\right\}+\log \left[{\displaystyle \frac{p_d\left({\mathit{\xi }}_d\right)}{p_d\left(d\right)}}\right]\\ & & +\log \left\{{\displaystyle \frac{{\mathrm{\Phi }}^{\left(\mathcal{N}\right)}\left[\left(\frac{\mathrm{1}}{\mathrm{2}},-,d\right),/,{\mathit{\sigma }}_d\right]-{\mathrm{\Phi }}^{\left(\mathcal{N}\right)}\left[\left(-,\frac{\mathrm{1}}{\mathrm{2}},-,d\right),/,{\mathit{\sigma }}_d\right]}{{\mathrm{\Phi }}^{\left(\mathcal{N}\right)}\left[\left(\frac{\mathrm{1}}{\mathrm{2}},-,{\mathit{\xi }}_d\right),/,{\mathit{\sigma }}_d\right]-{\mathrm{\Phi }}^{\left(\mathcal{N}\right)}\left[\left(-,\frac{\mathrm{1}}{\mathrm{2}},-,{\mathit{\xi }}_d\right),/,{\mathit{\sigma }}_d\right]}}\right\}.\end{array}$$

Optional update of ${\mathit{x}}_{\text{A}}$: when using the approximate likelihood method, one must account for the auxiliary variables ${\mathit{x}}_{\text{A}}$, a $\mathit{P}$ vector (e.g., $P$ $=$ $n$). We find that, in practice, it is not necessary to update all the auxiliary parameters at each iteration. In fact the method can be shown to work perfectly well, empirically, if we never update them, provided they are given a sensible initial value (such as the sample mean of the observed data $\stackrel{\mathrm{‾}}{x}$). This is not an uncommon tactic in the long-memory (big-$n$) context e.g.,; for further discussion refer to Appendix C.

For a full-MH approach, we recommend an independence sampler to backward project the observed time series. Specifically, first relabel the observed data: $y_{-i}$ $=$ $x_{i+\mathrm{1}}$, $i$ $=$ 0, … $n$ $-$ 1; furthermore, use the vector ($y_{-(n-\mathrm{1})}$, …, $y_{-\mathrm{1}}$, $y_{\mathrm{0}}{)}^t$ to generate a new vector of length $n$, ($Y_{\mathrm{1}}$, …, $Y_n{)}^t$ where $Y_t$ via Eq. (): $Y_t$ $=$ ${\mathit{\epsilon }}_t$ $+$ ${\mathrm{\Pi }}_P\mathit{\mu }$ $-$ $\sum _{k=\mathrm{1}}^n{\mathit{\pi }}_k{Y}_{t-k}$, where the coefficients $\mathit{\left\{}\mathit{\pi }\mathit{\right\}}$ are determined by the current value of the memory parameter(s). Then take the proposed ${\mathit{x}}_{\text{A}}$, denoted ${\mathit{\xi }}_{x_{\text{A}}}$, as the reverse sequence: ${\mathit{\xi }}_{x_{-i}}$ $=$ $y_{i+\mathrm{1}}$, $i$ $=$ 0, …, $n$ $-$ 1. Since this is an independence sampler, calculation of the acceptance probability is straightforward. It is only necessary to evaluate the proposal density $q\left({\mathit{\xi }}_{x_{\text{A}}}\right|\mathit{x}$, $\mathit{\psi })$. But this is easy using the results from Sect. . For simplicity, we prefer uniform prior for ${\mathit{x}}_{\text{A}}$.

Besides simplicity, justification for this approach lies primarily in is preservation of the autocorrelation structure – this is clear since the ACF is symmetric in time. The proposed vector has a low acceptance rate, and the potential remedies (e.g., multiple-try methods) seem unnecessarily complicated given the success of the simpler method.

Extensions to accommodate short memory

Simple ARFIMA(0,$d$,0) models are mathematically convenient but have limited practical applicability because the entire memory structure is determined by just one parameter, $d$. Although $d$ is often of primary interest, it may be unrealistic to assume no short-memory effects. This issue is often implicitly acknowledged since semi-parametric estimation methods, such as those used as comparators in Sect. , are motivated by a desire to circumvent the problem of specifying precisely (and inferring) the form of short memory (i.e., the values of $p$ and $q$ in an ARIMA model). Full parametric Bayesian modeling of ARFIMA($p$,$d$,$q$) processes represents an essentially untried alternative, primarily due to computational challenges. Related, more discrete, alternatives show potential. considered all four models with $p$, $q$ $\le$ 1, whereas considered 16 with $p$, $q$ $\le$ 3.

Such approaches, especially ones allowing larger $p$, $q$, can be computationally burdensome as much effort is spent modeling unsuitable processes towards a goal (inferring $p$, $q$), which is not of primary interest ($d$ is). To develop an efficient, fully parametric, Bayesian method of inference that properly accounts for varying models, and to marginalize out these nuisance quantities, we use reversible-jump (RJ) MCMC . We extend the parameter space to include the set of models ($p$ and $q$), with chains moving between (i.e., changing $p$ and/or $q$) and within (sampling $\mathit{\varphi }$ and $\mathit{\theta }$ given particular fixed $p$ and $q$) models, and focus on the marginal posterior distribution of $d$ obtained by (Monte Carlo) integration over all models and parameters therein. RJ methods, which mixes so-called trans-dimensional, between-model moves with the conventional within-model ones, have previously been applied to both autoregressive models , and full-ARMA models . In the long-memory context, applied RJ to fractional exponential processes (FEXP). However for ARFIMA, the only related work we are aware of is by who demonstrated a promising if limited alternative.

Below we show how the likelihood may be calculated with extra short-memory components when $p$ and $q$ are known, and subsequently how Bayesian inference can be applied in this case. Then, the more general case of unknown $p$ and $q$ via RJ is described. The result is a Monte Carlo inferential scheme that allows short-memory effects to be marginalized out when summarizing inferences for the main parameter of interest: $d$, for long memory.

Likelihood derivation and inference for known short memory

Recall that short-memory components of an ARFIMA process are defined by the AR and moving average (MA) polynomials, $\mathrm{\Phi }$ and $\mathrm{\Theta }$, respectively (see Sect. ). Here, we distinguish between the polynomial, $\mathrm{\Phi }$, and the vector of its coefficients, $\mathit{\varphi }$ $=$ (${\mathit{\varphi }}_{\mathrm{1}}$, …, ${\mathit{\varphi }}_p$). When the polynomial degree is required explicitly, bracketed superscripts will be used: ${\mathrm{\Phi }}^{\left(p\right)}$, ${\mathit{\varphi }}^{\left(p\right)}$, ${\mathrm{\Theta }}^{\left(p\right)}$, ${\mathit{\theta }}^{\left(p\right)}$.

We combine the short-memory parameters $\mathit{\varphi }$ and $\mathit{\theta }$ with $d$ to create a single memory parameter, $\mathit{\omega }$ $=$ ($\mathit{\varphi }$, $\mathit{\theta }$, $d$). For a given unit-variance ARFIMA($p$,$d$,$q$) process, we denote its autocovariance (ACV) by ${\mathit{\gamma }}_{\mathit{\omega }}(\cdot )$, with ${\mathit{\gamma }}_d(\cdot )$ and ${\mathit{\gamma }}_{\mathit{\varphi }}$, $\mathit{\theta }(\cdot )$ those of the relevant unit-variance ARFIMA(0,$d$,0) and ARMA($p$,$q$) processes, respectively. The spectral density function (SDF) of the unit-variance ARFIMA($p$,$d$,$q$) process is written as $f_{\mathit{\omega }}(\cdot )$, and its covariance matrix is ${\mathbf{\Sigma }}_{\mathit{\omega }}$.

An exact likelihood evaluation requires an explicit calculation of the ACV ${\mathit{\gamma }}_{\mathit{\omega }}(\cdot )$; however, there is no simple closed form for arbitrary ARFIMA processes. Fortunately, our proposed approximate likelihood method of section can be ported over directly. Given the coefficients $\mathit{\left\{}{\mathit{\pi }}_k^{\left(d\right)}\mathit{\right\}}$ and polynomials $\mathrm{\Phi }$ and $\mathrm{\Theta }$, it is straightforward to calculate the $\mathit{\left\{}{\mathit{\pi }}_k^{\left(\mathit{\omega }\right)}\mathit{\right\}}$ coefficients required by again applying the numerical methods of Sect. 3.3.

To focus the exposition, consider the simple, yet useful, ARFIMA(1,$d$,0) model where the full memory parameter is $\mathit{\omega }$ $=$ ($d$, ${\mathit{\varphi }}_{\mathrm{1}}$). Because the parameter spaces of $d$ and ${\mathit{\varphi }}_{\mathrm{1}}$ are independent, it is simplest to update each of these parameters separately; $d$ with the methods of Sect.  and ${\mathit{\varphi }}_{\mathrm{1}}$ similarly: ${\mathit{\xi }}_{{\mathit{\varphi }}_{\mathrm{1}}}|{\mathit{\varphi }}_{\mathrm{1}}$ $\sim$ ${\mathcal{N}}^{(-\mathrm{1},\mathrm{1})}({\mathit{\varphi }}_{\mathrm{1}}$, ${\mathit{\sigma }}_{{\mathit{\varphi }}_{\mathrm{1}}}^{\mathrm{2}})$, for some ${\mathit{\sigma }}_{{\mathit{\varphi }}_{\mathrm{1}}}^{\mathrm{2}}$. In practice however, the posteriors of $d$ and ${\mathit{\varphi }}_{\mathrm{1}}$ typically exhibit significant correlation so independent proposals are inefficient. One solution would be to parametrize to some $d^{*}$ and orthogonal ${\mathit{\varphi }}_{\mathrm{2}}^{*}$, but the interpretation of $d^{*}$ would not be clear. An alternative to explicit reparametrization is to update the parameters jointly, but in such a way that proposals are aligned with the correlation structure. This will ensure a reasonable acceptance rate and mixing.

To propose parameters in the manner described above, a two-dimensional, suitably truncated Gaussian random walk, with covariance matrix aligned with the posterior covariance, is required. To make proposals of this sort, and indeed for arbitrary $\mathit{\omega }$ in larger $p$ and $q$ cases, requires sampling from a hypercuboid-truncated multivariate normal (MVN) ${\mathcal{N}}_r^{(a,b)}(\mathit{\omega }$, ${\mathrm{\Sigma }}_{\mathit{\omega }}$), where ($\mathbf{a}$, $\mathbf{b}$) describe the coordinates of the hypercube. We find that rejection sampling-based unconstrained similarly parametrized MVN samples (e.g., using mvtnorm, ) works well, because in the RW setup the mode of the distribution always lies inside the hypercuboid. Returning to the specific ARFIMA(1,$d$,0) case, $r$ $=$ 2, $\mathbf{b}$ $=$ (0.5, 1), and $\mathbf{a}$ $=$ $-\mathbf{b}$ are appropriate choices. Calculation of the MH acceptance ratio $A_{\mathit{\omega }}(\mathit{\omega }$, ${\mathit{\xi }}_{\mathit{\omega }}$) is trivial; it simply requires numerical evaluation of ${\mathrm{\Phi }}_r^{\left(\mathcal{N}\right)}(\cdot$; $\cdot$, ${\mathbf{\Sigma }}_{\mathit{\omega }}$), e.g., via mvtnorm, since the ratios of hypercuboid normalization terms would cancel. We find that initial values ${\mathit{\varphi }}^{\left[\mathrm{0}\right]}$ chosen uniformly in ${\mathcal{C}}_{\mathrm{1}}$; i.e., the interval ($-$1, 1), and $d^{\left[\mathrm{0}\right]}$ from $\mathit{\{}-$0.4, $-$0.2, 0, 0.2, 0.4$\mathit{\}}$ work well. Any choice of prior for $\mathit{\omega }$ can be made, although we prefer flat (proper) priors.

The only technical difficulty is the choice of proposal covariance matrix ${\mathbf{\Sigma }}_{\mathit{\omega }}$. Ideally, it would be aligned with the posterior covariance; however, this is not a priori known. We find that running a pilot chain with independent proposals via ${\mathcal{N}}_r^{(a,b)}(\mathit{\omega }$, ${\mathit{\sigma }}_{\mathit{\omega }}^{\mathrm{2}}{\mathbf{I}}_r$) can help choose a ${\mathbf{\Sigma }}_{\mathit{\omega }}$. A rescaled version of the sample covariance matrix from the pilot posterior chain, following , works well (see Sect. ).

Unknown short-memory form

We now expand the parameter space to include models $M$ $\in$ $\mathcal{M}$, the set of ARFIMA models with $p$ and $q$ short-memory parameters, indexing the size of the parameter space ${\mathrm{\Psi }}^{\left(M\right)}$. For our trans-dimensional moves, we only consider adjacent models, on which we will be more specific later. For now, note that the choice of bijective function mapping between model spaces (whose Jacobian term appears in the acceptance ratio) is crucial to the success of the sampler. To illustrate, consider transforming from ${\mathrm{\Phi }}^{(p+\mathrm{1})}$ $\in$ ${\mathcal{C}}_{p+\mathrm{1}}$ down to ${\mathrm{\Phi }}^{\left(p\right)}$ $\in$ ${\mathcal{C}}_p$. This turns out to be a non-trivial problem, however, because (for $p$ $>$ 1) ${\mathcal{C}}_p$ has a very complicated shape. The most natural map would be (${\mathit{\varphi }}_{\mathrm{1}}$, …, ${\mathit{\varphi }}_p$, ${\mathit{\varphi }}_{p+\mathrm{1}}$) $⟼$ (${\mathit{\varphi }}_{\mathrm{1}}$, …, ${\mathit{\varphi }}_p$). However, there is no guarantee that the image will lie in ${\mathcal{C}}_p$. Even if the model dimension is fixed, difficulties are still encountered; a natural proposal method would be to update each component of $\mathit{\varphi }$ separately but, because of the awkward shape of ${\mathcal{C}}_p$, the allowable values for each component are a complicated function of the others. Non-trivial proposals are required.

A potential approach is to parametrize in terms of the inverse roots (poles) of $\mathrm{\Phi }$, as advocated by : by writing $\mathrm{\Phi }\left(z\right)$ $=$ ${\prod }_{i=\mathrm{1}}^p$(1 $-$ ${\mathit{\alpha }}_{i}z$), we have ${\mathit{\varphi }}^{\left(p\right)}$ $\in$ ${\mathcal{C}}_p$ $⟺$ $\left|{\mathit{\alpha }}_i\right|$ $<$ 1 for all $i$. This looks attractive because it transforms ${\mathcal{C}}_p$ into $D^p$ $=$ $D$ $\times$ $\mathrm{\cdots }$ $\times$ $D$ ($p$ times) where $D$ is the open unit disc, which is easy to sample from. But this method has serious drawbacks when we consider the RJ step. To decrease dimension, the natural map would be to remove one of the roots from the polynomial. But because it is assumed that $\mathrm{\Phi }$ has real coefficients (otherwise the model has no realistic interpretation), any complex ${\mathit{\alpha }}_i$ must appear as conjugate pairs. There is then no obvious way to remove a root; a contrived method might be to remove the conjugate pair and replace it with a real root with the same modulus; however, it is unclear how this new polynomial is related to the original, and to other aspects of the process, like ACV.

Reparametrization of $\mathrm{\Phi }$ and $\mathrm{\Theta }$

We therefore propose reparametrization $\mathrm{\Phi }$ (and $\mathrm{\Theta }$) using the bijection between ${\mathcal{C}}_p$ and ($-$1, 1)${}^p$ advocated by various authors, e.g.,  and . To our knowledge, these methods have not previously been deployed towards integrating out short-memory components in Bayesian analysis of ARFIMA processes.

defined a mapping ${\mathit{\varphi }}^{\left(p\right)}$ $⟷$ ${\mathit{\phi }}^{\left(p\right)}$ recursively as follows: $${\mathit{\varphi }}_i^{(k-\mathrm{1})}={\displaystyle \frac{{\mathit{\varphi }}_i^{\left(k\right)}-{\mathit{\varphi }}_k^{\left(k\right)}{\mathit{\varphi }}_{k-i}^{\left(k\right)}}{\mathrm{1}-{\left({\mathit{\varphi }}_k^{\left(k\right)}\right)}^{\mathrm{2}}}},k=p,\mathrm{\dots },\mathrm{2},i=\mathrm{1},\mathrm{\dots },k-1.$$ Then set ${\mathit{\phi }}_k^{\left(p\right)}$ $=$ ${\mathit{\varphi }}_k^{\left(k\right)}$ for $k$ $=$ 1, …, $p$. The reverse recursion is given by $${\mathit{\varphi }}_i^{\left(k\right)}=\left\{\begin{array}{lcll}{\mathit{\phi }}_k^{\left(p\right)}& \text{for}& i=k& k=\mathrm{1},\mathrm{\dots },p\\ {\mathit{\varphi }}_i^{(k-\mathrm{1})}+{\mathit{\phi }}_k^{\left(p\right)}{\mathit{\varphi }}_{k-i}^{(k-\mathrm{1})}& \text{for}& i=\mathrm{1},\mathrm{\dots },k-\mathrm{1}& k=\mathrm{2},\mathrm{\dots },p\end{array}\right..$$ Note that ${\mathit{\varphi }}_p^{\left(p\right)}$ $=$ ${\mathit{\phi }}_p^{\left(p\right)}$. Moreover, if $p$ $=$ 1, the two parametrizations are the same, i.e., ${\mathit{\varphi }}_{\mathrm{1}}$ $=$ ${\mathit{\phi }}_{\mathrm{1}}$ (consequently the brief study of ARFIMA(1,$d$,0) in Sect.  fits in this framework). The equivalent parametrized form for $\mathit{\theta }$ is $\mathit{\vartheta }$. The full memory parameter $\mathit{\omega }$ is parametrized as $\stackrel{\mathrm{‾}}{\mathrm{\Omega }}$ $=$ ($-\mathrm{1}/\mathrm{2}$, $\mathrm{1}/\mathrm{2}$) $\times$ (the image of ${\mathcal{C}}_{p,q}$). However, recall that in practice, ${\mathcal{C}}_{p,q}$ will be assumed equivalent to ${\mathcal{C}}_p$ $\times$ ${\mathcal{C}}_q$, so the parameter space is effectively $\stackrel{\mathrm{‾}}{\mathrm{\Omega }}$ $=$ ($-\mathrm{1}/\mathrm{2}$, $\mathrm{1}/\mathrm{2}$) $\times$ ($-$1, 1)${}^{p+q}$.

Besides mathematical convenience, this bijection has a very useful property cf., which helps motivate its use in defining RJ maps. If $d$ $=$ $q$ $=$ 0, using this parametrization for $\mathit{\phi }$ when moving between different values of $p$ allows one to automatically choose processes that have very closely matching ACFs at low lags. In the MCMC context this is useful because it allows the chain to propose models that have a similar correlation structure to the current one. Although this property is nice, it may be of limited value for full-ARFIMA models, since the proof of the main result does not easily lend itself to the inclusion of either a MA or long-memory component. Nevertheless, our empirical results similarly indicate a near match for a full-ARFIMA($p$,$d$,$q$) model.

Application of RJ MCMC to ARFIMA($p$,$d$,$q$) processes

We now use this reparametrization to efficiently propose new parameter values. Firstly, it is necessary to propose a new memory parameter $\mathit{\varpi }$ while keeping the model fixed. Attempts at updating each component individually suffer from the same problems of excessive posterior correlation that were encountered in Sect. . Therefore, the simultaneous update of the entire $r$ $=$ ($p$ $+$ $q$ $+$ 1)-dimensional parameter $\mathit{\varpi }$ is performed using the hypercuboid-truncated Gaussian distribution from definition ${\mathit{\xi }}_{\mathit{\varpi }}|\mathit{\varpi }$ $\sim$ ${\mathcal{N}}_r^{{\mathcal{H}}_r}(\mathit{\varpi }$, ${\mathbf{\Sigma }}_{\mathit{\varpi }})$, where ${\mathcal{H}}_r$ defines the $r$-dimensional rectangle. The covariance matrix ${\mathbf{\Sigma }}_{\mathit{\varpi }}$ is discussed in some detail below. The choice of prior $p_{\mathit{\varpi }}(\cdot )$ is arbitrary. used a uniform prior for $\mathit{\omega }$, which has an explicit expression in the $\mathit{\varpi }$ parametrization . However, their expression is unnecessarily complicated since a uniform prior over $\mathrm{\Omega }$ holds no special interpretation. We therefore prefer uniform prior over $\stackrel{\mathrm{‾}}{\mathrm{\Omega }}$: $p_{\mathit{\varpi }}\left(\mathit{\varpi }\right)$ $\propto$ 1, $\mathit{\varpi }$ $\in$ $\stackrel{\mathrm{‾}}{\mathrm{\Omega }}$.

Now consider the between-model transition. We must first choose a model prior $p_{\mathcal{M}}(\cdot )$. A variety of priors are possible; the simplest option would be to have a uniform prior over $\mathcal{M}$, but this would of course be improper. We may in practice want to restrict the possible values of $p$, $q$ to 0 $\le$ $p$ $\le$ $P$ and 0 $\le$ $q$ $\le$ $Q$ for some $P$, $Q$ (say 5), which would render the uniform prior proper. However, even in this formulation, a lot of prior weight is being put onto (larger) more complicated models that, in the interests of parsimony, might be undesired. As a simple representative of potential priors that give greater weight to smaller models, we prefer a truncated joint Poisson distribution with parameter $\mathit{\lambda }$: $p_{\mathcal{M}}(p$, $q)$ $\propto$ $\frac{{\mathit{\lambda }}^{p+q}}{p\mathrm{!}q\mathrm{!}}\mathbb{I}(p$ $\le$ $P$, $q$ $\le$ $Q)$.

Now, denote the probability of jumping from model $M_{p,q}$ to model $M_{p^{\prime },q^{\prime }}$ by $U_{(p,q),(p^{\prime },q^{\prime })}$. $U$ could allocate non-zero probability for every model pair, but for convenience we severely restrict the possible jumps (while retaining irreducibility) using a two-dimensional bounded birth and death process. Consider the subgraph of ${\mathbb{Z}}^{\mathrm{2}}$: $G$ $=$ $\mathit{\left\{}\right(p,q)$: 0 $\le$ $p$ $\le$ $P$, 0 $\le$ $q$ $\le$ $Q\mathit{\}}$, and allocate uniform non-zero probability only to neighboring values, i.e., if and only if $|p$ $-$ $p^{\prime }|$ $+$ $|q$ $-$ $q^{\prime }|$ $=$ 1. Each point in the body of $G$ has four neighbors, each point on the line boundaries has three, and each of the four corner points has only two neighbors. Therefore, the model transition probabilities $U_{(p,q),(p^{\prime },q^{\prime })}$ are either $\mathrm{1}/\mathrm{4}$, $\mathrm{1}/\mathrm{3}$, $\mathrm{1}/\mathrm{2}$, or 0.

Now suppose the current ($p$ $+$ $q$ $+$ 3)-dimensional parameter is ${\mathit{\psi }}^{(p,q)}$, given by ${\mathit{\psi }}^{(p,q)}$ $=$ ($\mathit{\mu }$, $\mathit{\sigma }$, $d$, ${\mathit{\phi }}^{\left(p\right)}$, ${\mathit{\vartheta }}^{\left(q\right)}$), using a slight abuse of notation. Because the mathematical detail of the AR and MA components are almost identical, we consider only the case of decreasing/increasing $p$ by 1 here; all of the following remains valid if $p$ is replaced by $q$, and $\mathit{\phi }$ replaced by $\mathit{\vartheta }$. We therefore seek to propose a parameter ${\mathit{\xi }}^{(p+\mathrm{1},q)}$ $=$ (${\mathit{\xi }}_{\mathit{\mu }}$, ${\mathit{\xi }}_{\mathit{\sigma }}$, ${\mathit{\xi }}_d$, ${\mathit{\xi }}_{\mathit{\phi }}^{(p+\mathrm{1})}$, ${\mathit{\xi }}_{\mathit{\vartheta }}^{\left(q\right)}$), that is somehow based on ${\mathit{\psi }}^{(p,q)}$. We further simplify by regarding the other three parameters ($\mathit{\mu }$, $\mathit{\sigma }$, and $d$) as having the same interpretation in every model, choosing ${\mathit{\xi }}_{\mathit{\mu }}$ $=$ $\mathit{\mu }$, ${\mathit{\xi }}_{\mathit{\sigma }}$ $=$ $\mathit{\sigma }$ and ${\mathit{\xi }}_d$ $=$ $d$. For simplicity we also set ${\mathit{\xi }}_{\mathit{\vartheta }}^{\left(q\right)}$ $=$ ${\mathit{\vartheta }}^{\left(q\right)}$. Now consider the map ${\mathit{\phi }}^{\left(p\right)}$ $\to$ ${\mathit{\xi }}_{\mathit{\phi }}^{(p+\mathrm{1})}$. To specify a bijection, we match dimensions by adding in a random scalar $u$. The most obvious map is to specify $u$, so that its support is the interval ($-$1, 1) and then set: ${\mathit{\xi }}_{\mathit{\phi }}^{(p+\mathrm{1})}$ $=$ (${\mathit{\phi }}^{\left(p\right)}$, $u$). The corresponding map for decreasing the dimension ${\mathit{\phi }}^{(p+\mathrm{1})}$ $\to$ ${\mathit{\xi }}_{\mathit{\phi }}^{\left(p\right)}$ is ${\mathit{\xi }}_{\mathit{\phi }}^{\left(p\right)}$ $=$ (${\mathit{\phi }}_{\mathrm{1}}^{(p+\mathrm{1})}$, …, ${\mathit{\phi }}_p^{(p+\mathrm{1})}$). We either add, or remove the final parameter, while keeping all others fixed with the identity map, so the Jacobian is unity. The proposal $q\left(u\right|{\mathit{\psi }}^{(p,q)})$ can be made in many ways – we prefer the simple $\mathcal{U}$($-$1,1). With these choices the RJ acceptance ratio is $$\begin{array}{cc}A& ={\mathrm{\ell }}_{(p^{\prime },q^{\prime })}\left(\mathit{x},|,{\mathit{\xi }}^{(p^{\prime },q^{\prime })}\right)-{\mathrm{\ell }}_{(p,q)}\left(\mathit{x},|,{\mathit{\psi }}^{(p,q)}\right)\\ & +\log \left\{{\displaystyle \frac{p_{\mathcal{M}}(p^{\prime },q^{\prime })}{p_{\mathcal{M}}(p,q)}},{\displaystyle \frac{U_{(p^{\prime },q^{\prime }),(p,q)}}{U_{(p,q),(p^{\prime },q^{\prime })}}}\right\},\end{array}$$ which applies to both increasing and decreasing dimensional moves.

Posterior outputs; (a) Bayesian estimate ${\widehat{d}}^{\left(\text{B}\right)}$ values on the $y$ axis against the true $d_I$ on the $x$ axis, (b) residuals ${\widehat{d_{\text{R}}}}^{\left(\text{B}\right)}$ from the Bayesian estimate from the truth against that truth, $d_I$. Each “x” plotted represts one estimate or residual.

[Figure omitted. See PDF]

Construction of ${\mathbf{\Sigma }}_{\mathit{\varpi }}$: much of the efficiency of the above scheme, including within- and between-model moves, depends on the choice of ${\mathbf{\Sigma }}_{\mathit{\varpi }}$ $\equiv$ ${\mathrm{\Sigma }}^{(p,q)}$, the within-model move RW proposal covariance matrix. We first seek an appropriate ${\mathrm{\Sigma }}^{(\mathrm{1},\mathrm{1})}$, as in Sect. , with a pilot tuning scheme. That matrix is shown on the left below, where we have blocked it out $$\begin{array}{cc}& {\mathrm{\Sigma }}^{(\mathrm{1},\mathrm{1})}=\left(\begin{array}{ccc}{\mathit{\sigma }}_d^{\mathrm{2}}& {\mathit{\sigma }}_{d,{\mathit{\phi }}_{\mathrm{1}}}& {\mathit{\sigma }}_{d,{\mathit{\vartheta }}_{\mathrm{1}}}\\ & {\mathit{\sigma }}_{{\mathit{\phi }}_{\mathrm{1}}}^{\mathrm{2}}& {\mathit{\sigma }}_{{\mathit{\phi }}_{\mathrm{1}},{\mathit{\vartheta }}_{\mathrm{1}}}\\ & & {\mathit{\sigma }}_{{\mathit{\vartheta }}_{\mathrm{1}}}^{\mathrm{2}}\end{array}\right),\\ & & {\mathrm{\Sigma }}^{(p,q)}=\left(\begin{array}{ccc}{\mathit{\sigma }}_d^{\mathrm{2}}& {\mathrm{\Sigma }}_{d,\mathit{\phi }}^{\left(p\right)}& {\mathrm{\Sigma }}_{d,\mathit{\vartheta }}^{\left(q\right)}\\ & {\mathbf{\Sigma }}_{\mathit{\phi }}^{\left(p\right)},{\mathit{\phi }}^{\left(p\right)}& {\mathbf{\Sigma }}_{\mathit{\phi }}^{\left(p\right)},{\mathit{\vartheta }}^{\left(q\right)}\\ & & {\mathrm{\Sigma }}_{\mathit{\vartheta }}^{\left(q\right)},{\mathit{\vartheta }}^{\left(q\right)}\end{array}\right)\end{array}$$ (where each block is a scalar), so that we can extend this idea to the ($p$, $q$) case in the obvious way – on the right above – where ${\mathbf{\Sigma }}_{\mathit{\phi }}^{\left(p\right)}$, ${\mathit{\phi }}^{\left(p\right)}$ is a $p$ $\times$ $p$ matrix, ${\mathrm{\Sigma }}_{\mathit{\vartheta }}^{\left(q\right)}$, ${\mathit{\vartheta }}^{\left(q\right)}$ is a $q$ $\times$ $q$ matrix, etc. If either (or both) $p$, $q$ $=$ 0 then the relevant blocks are simply omitted. To specify the various sub-matrices, we propose ${\mathit{\phi }}_{\mathrm{2}}$, …, ${\mathit{\phi }}_p$ with equal variances, and independently of $d$, ${\mathit{\phi }}_{\mathrm{1}}$, ${\mathit{\vartheta }}_{\mathrm{1}}$, (and similarly for ${\mathit{\vartheta }}_{\mathrm{2}}$, …, ${\mathit{\vartheta }}_q$). In the context of Eq. (), the following holds true: $$\begin{array}{cc}& {\mathrm{\Sigma }}_{d,\mathit{\phi }}^{\left(p\right)}=\left(\begin{array}{cc}{\mathit{\sigma }}_{d,{\mathit{\phi }}_{\mathrm{1}}}& \mathbf{0}\end{array}\right),\\ & {\mathrm{\Sigma }}_{d,\mathit{\vartheta }}^{\left(q\right)}=\left(\begin{array}{cc}{\mathit{\sigma }}_{d,{\mathit{\vartheta }}_{\mathrm{1}}}& \mathbf{0}\end{array}\right),\end{array}$$ $$\begin{array}{cc}& {\mathbf{\Sigma }}_{\mathit{\phi }}^{\left(p\right)},{\mathit{\phi }}^{\left(p\right)}=\left(\begin{array}{cc}{\mathit{\sigma }}_{{\mathit{\phi }}_{\mathrm{1}}}^{\mathrm{2}}& \mathbf{0}\\ \mathbf{0}& {\mathit{\sigma }}_{\mathit{\phi }}^{\mathrm{2}}{I}_{p-\mathrm{1}}\end{array}\right),\\ & {\mathrm{\Sigma }}_{\mathit{\vartheta }}^{\left(q\right)},{\mathit{\vartheta }}^{\left(q\right)}=\left(\begin{array}{cc}{\mathit{\sigma }}_{{\mathit{\vartheta }}_{\mathrm{1}}}^{\mathrm{2}}& \mathbf{0}\\ \mathbf{0}& {\mathit{\sigma }}_{\mathit{\vartheta }}^{\mathrm{2}}{I}_{q-\mathrm{1}}\end{array}\right),\\ & {\mathbf{\Sigma }}_{\mathit{\phi }}^{\left(p\right)},{\mathit{\vartheta }}^{\left(q\right)}=\left(\begin{array}{cc}{\mathit{\sigma }}_{{\mathit{\phi }}_{\mathrm{1}},{\mathit{\vartheta }}_{\mathrm{1}}}& \mathbf{0}\\ \mathbf{.0}& \mathbf{O}\end{array}\right),\end{array}$$ where the dotted lines indicate further blocking, $\mathbf{0}$ is a row vector of zeros, and $\mathbf{O}$ is a zero matrix. This choice of ${\mathbf{\Sigma }}_{\mathit{\phi }}$ is conceptually simple, computationally easy and preserves the positive definiteness as required (see ).

Empirical illustration and comparison

Here we provide empirical illustrations for the methods above: for classical and Bayesian analysis of long-memory models, and extensions for short memory. To ensure consistency throughout, the location and scale parameters will always be chosen as ${\mathit{\mu }}_I$ $=$ 0 and ${\mathit{\sigma }}_I$ $=$ 1. Furthermore, unless stated otherwise, the simulated series will be of length $n$ $=$ 2${}^{10}$ $=$ 1024. This is a reasonable size for many applications; it is equivalent to 85 years of monthly observations. When using the approximate likelihood method we set $P$ $=$ $n$.

Long memory

Standard MCMC diagnostics were used throughout to ensure, and tune for, good mixing. Because $d$ is the parameter of primary interest, the initial values $d^{\left[\mathrm{0}\right]}$ will be chosen to systematically cover its parameter space, usually starting five chains at the regularly spaced points $\mathit{\{}-$0.4, $-$0.2, 0, 0.2, 0.4$\mathit{\}}$. Initial values for other parameters are not varied: $\mathit{\mu }$ will start at the sample mean $\stackrel{\mathrm{‾}}{x}$; $\mathit{\sigma }$ at the sample standard deviation of the observed series $\mathit{x}$.

Efficacy of approximate likelihood method

We start with the null case; i.e., how does the algorithm perform when the data are not from a long-memory process? One hundred independent ARFIMA(0,0,0), or Gaussian white noise, processes are simulated, from which marginal posterior means, standard deviations, and credibility interval end points are extracted. Table  shows averages over the runs.

The average estimate for each of the three parameters is less than a quarter of a standard deviation away from the truth. Credibility intervals are nearly symmetric about the estimate and the marginal posteriors are, to a good approximation, locally Gaussian (not shown). Upon, applying a proxy credible-interval-based hypothesis test, one would conclude in 98 of the cases that $d$ $=$ 0 could not be ruled out. A similar analysis for $\mathit{\mu }$ and $\mathit{\sigma }$ shows that hypotheses $\mathit{\mu }$ $=$ 0 and $\mathit{\sigma }$ $=$ 1 would each have been accepted 96 times. These results indicate that the 95 % credibility intervals are approximately correctly sized.

Next, consider the more interesting case of $d_I$ $\ne$ 0. We repeat the above experiment except that 10 processes are generated with $d_I$ set to each of $\mathit{\{}-$0.45, $-$0.35, …, 0.45$\mathit{\}}$, giving 100 series total. Figure  shows a graphical analog of results from this experiment. The plot axes involve a Bayesian residual estimate of $d$, ${\widehat{d_{\text{R}}}}^{\left(\text{B}\right)}$, defined as ${\widehat{d_{\text{R}}}}^{\left(\text{B}\right)}$ $=$ ${\widehat{d}}^{\left(\text{B}\right)}$ $-$ $d_I$, where ${\widehat{d}}^{\left(\text{B}\right)}$ is the Bayesian estimate of $d$.

Posterior outputs: (a) Bayesian estimated standard deviation ${\widehat{\mathit{\sigma }}}^{\left(\text{B}\right)}$ against true $d_I$ values; (b) Bayesian estimated mean ${\widehat{\mathit{\mu }}}^{\left(\text{B}\right)}$ against $d_I$; and (c) uncertainty in the posterior for $\mathit{\mu }$, the standard deviation ${\widehat{{\mathit{\sigma }}_{\mathit{\mu }}}}^{\left(\text{B}\right)}$ against $d_I$ (semi-log scale). Each “x” plotted corresponds to an estimate.

[Figure omitted. See PDF]

Posterior summary statistics for an ARFIMA(0,0,0) process. Results are based on averaging over 100 independent ARFIMA(0,0,0) simulations for the long-memory parameter $d$, mean $\mathit{\mu }$ and noise variance $\mathit{\sigma }$.

	Mean	SD	95 % CI
$d$	0.006	0.025	$-$0.042	0.055
$\mathit{\mu }$	$-$0.004	0.035	$-$0.073	0.063
$\mathit{\sigma }$	1.002	0.022	0.956	1.041

From the figure is clear that the estimator for $d$ is performing well. Figure a shows how tight the estimates of $d$ are around the input value – recall that the parameter space for $d$ is the whole interval ($-\frac{\mathrm{1}}{\mathrm{2}}$, $\frac{\mathrm{1}}{\mathrm{2}}$). Moreover, Fig. b indicates that there is no significant change of posterior bias or variance as $d_I$ is varied.

Next, the corresponding plots for the parameters $\mathit{\sigma }$ and $\mathit{\mu }$ are shown in Fig. . We see from Fig. a that the estimate of $\mathit{\sigma }$ also appears to be unaffected by the input value $d_I$. The situation is different however in Fig. b for the location parameter $\mathit{\mu }$. Although the bias appears to be roughly zero for all $d_I$, the posterior variance clearly is affected by $d_I$. To ascertain the precise functional dependence, consider Fig. c, which shows, on a semi-log scale, the marginal posterior standard deviation of $\mathit{\mu }$, ${\widehat{{\mathit{\sigma }}_{\mathit{\mu }}}}^{\left(\text{B}\right)}$ against $d_I$.

It appears that the marginal posterior standard deviation ${\widehat{{\mathit{\sigma }}_{\mathit{\mu }}}}^{\left(\text{B}\right)}$ is a function of $d_I$; specifically, ${\widehat{{\mathit{\sigma }}_{\mathit{\mu }}}}^{\left(\text{B}\right)}$ $\propto$ $A^{d_I}$, for some $A$. The constant $A$ could be estimated via least-squares regression. Instead however, inspired by asymptotic results in literature concerning classical estimation of long-memory processes , we set $A=n$ and plotted the best-fitting such line (shown in Fig. c). Observe that, although not fitting exactly, the relation ${\widehat{{\mathit{\sigma }}_{\mathit{\mu }}}}^{\left(\text{B}\right)}$ $\propto$ $n^{d_I}$ holds reasonably well for $d_I$ $\in$ ($-\frac{\mathrm{1}}{\mathrm{2}}$, $\frac{\mathrm{1}}{\mathrm{2}}$). Indeed, motivated long memory in this way, and derived asymptotic consistency results for optimum (likelihood-based) estimators and found indeed that the standard error for $\mathit{\mu }$ is proportional to $n^{d-\mathrm{1}/\mathrm{2}}$ but the standard errors of all other parameters are proportional to $n^{-\mathrm{1}/\mathrm{2}}$.

Posterior outputs from an ARFIMA(0,0,0) series: (a) the posterior standard deviation in $d$, ${\widehat{{\mathit{\sigma }}_{\text{d}}}}^{\left(\text{B}\right)}$ against the sample size $n$; (b) posterior standard deviation in $\mathit{\mu }$, ${\widehat{{\mathit{\sigma }}_{\mathit{\mu }}}}^{\left(\text{B}\right)}$ against $n$; and (c) ${\widehat{{\mathit{\sigma }}_{\mathit{\sigma }}}}^{\left(\text{B}\right)}$ against $n$ (log–log scale).

[Figure omitted. See PDF]

Table: mean difference of estimates ${\widehat{d}}^{\left(\text{B}\right)}$ under alternative prior assumption. Plots: comparison of posteriors (solid lines) obtained under different priors (dotted lines). Time series used: ARFIMA(0,0.25,0) – (a) $n$ $=$ 2${}^{\mathrm{7}}$ $=$ 128, (b) $n$ $=$ 2${}^{10}$ $=$ 1024.

[Figure omitted. See PDF]

Effect of varying time series length

We now analyze the effect of changing the time series length. For this we conduct a similar experiment but fix $d_I$ $=$ 0 and vary $n$. The posterior statistics of interest are the posterior standard deviations ${\widehat{{\mathit{\sigma }}_d}}^{\left(\text{B}\right)}$, ${\widehat{{\mathit{\sigma }}_{\mathit{\mu }}}}^{\left(\text{B}\right)}$, and ${\widehat{{\mathit{\sigma }}_{\mathit{\sigma }}}}^{\left(\text{B}\right)}$. For each $n$ $\in$ $\mathit{\{}$128 $=$ 2${}^{\mathrm{7}}$, 2${}^{\mathrm{8}}$, …, 2${}^{14}$ $=$ 16 384$\mathit{\}}$, 10 independent ARFIMA(0,0,0) time series are generated. The resulting posterior standard deviations are plotted against $n$ (on log–log scale) in Fig. .

Observe that all three marginal posterior standard deviations are proportional to $\frac{\mathrm{1}}{\sqrt{n}}$, although the posterior of $\mathit{\mu }$ is less reliable. Combining these observations with our earlier deduction that ${\mathit{\sigma }}_{\mathit{\mu }}^{\left(\text{B}\right)}$ $\propto$ $n^{d_I}$, we conclude that for an ARFIMA(0,$d_I$,0) process of length $n$, the marginal posterior standard deviations follow those of .

Comparison with common estimators

In many practical applications, the long-memory parameter is estimated using non-/semi-parametric methods. These may be appropriate in many situations, where the exact form of the underlying process is unknown. However, when a specific model form is known (or at least assumed) they tend to perform poorly compared with fully parametric alternatives . Our aim here is to demonstrate, via a short Monte Carlo study involving ARFIMA(0,$d$,0) data, that our Bayesian likelihood-based method significantly outperforms other common methods in that case. We consider the following comparators: (i) rescaled adjusted range, or $R/S$ – we use the R implementation in the FGN package; (ii) semi-parametric Geweke–Porter–Hudak (GPH) method – implemented in R package fracdiff ; (iii) detrended fluctuation analysis (DFA), originally devised by – in the R package PowerSpectrum ; and (iv) wavelet-based semi-parametric estimators available in R package fARMA .

Each of these four methods will be applied to the same 100 time series with varying $d_I$ as were used earlier experiments above. We extend the idea of a residual, ${\widehat{d_{\text{R}}}}^{\left(\text{R}\right)}$, ${\widehat{d_{\text{R}}}}^{\left(G\right)}$, ${\widehat{d_{\text{R}}}}^{\left(D\right)}$, and ${\widehat{d_{\text{R}}}}^{\left(W\right)}$, to accommodate the new comparators, respectively, and plot them against ${\widehat{d_{\text{R}}}}^{\left(\text{B}\right)}$ in Fig. .

Observe that all four methods have a much larger variance than our Bayesian method, and moreover the $R/S$ is positively biased. Actually, the bias in some cases would seem to depend on $d_I$: $R/S$ is significantly (i.e., $>$ 0.25) biased for $d_I$ $<$ $-$0.3 but slightly negatively biased for $d$ $>$ 0.3 (not shown); DFA is only unbiased for $d_I$ $>$ 0; both the GPH and wavelet methods are unbiased for all $d$ $\in$ ($-\frac{\mathrm{1}}{\mathrm{2}}$, $\frac{\mathrm{1}}{\mathrm{2}}$).

Extensions for short memory

Known form: we first consider the MCMC algorithm from Sect.  for sampling under an ARFIMA(1,$d$,0) model where the full memory parameter is $\mathit{\omega }$ $=$ ($d$, ${\mathit{\varphi }}_{\mathrm{1}}$). Recall that method involved proposals from a hypercuboid MVN using a pilot-tuned covariance matrix. Also recall that it is a special case of the reparametrized method from Sect. .

In general, this method works very well; two example outputs are presented in Fig. , under two similar data-generating mechanisms.

Comparison of Bayesian estimator with common classical estimators: (a) $R/S$, (b) GPH, (c) DFA, and (d) wavelet.

[Figure omitted. See PDF]

Posterior samples of ($d$, $\mathit{\varphi }$): input time series (a) (1 $+$ 0.92$\mathcal{B}$)(1 $-$ $\mathcal{B}{)}^{0.25}{X}_t$ $=$ ${\mathit{\epsilon }}_t$, (b) (1 $-$ 0.83$\mathcal{B}$)(1 $-$ $\mathcal{B}{)}^{-0.35}{X}_t$ $=$ ${\mathit{\epsilon }}_t$.

[Figure omitted. See PDF]

Spectra for processes in Fig. . Green line is relevant ARMA(1,0) process, red line is relevant ARFIMA(0,$d$,0) process, black line is ARFIMA(1,$d$,0) process: (a) (1 $+$ 0.92$\mathcal{B}$)(1 $-$ $\mathcal{B}{)}^{0.25}{X}_t$ $=$ ${\mathit{\epsilon }}_t$; (b) (1 $-$ 0.83$\mathcal{B}$)(1 $-$ $\mathcal{B}{)}^{-0.35}{X}_t$ $=$ ${\mathit{\epsilon }}_t$. (c) Shows posterior samples of ($d$, $\mathit{\varphi }$) from series considered in (b) with credibility sets: red is 95 % credibility set for ($d$, $\mathit{\varphi }$), green is 95 % credibility interval for $d$, blue is 95 % credibility interval for $\mathit{\varphi }$.

[Figure omitted. See PDF]

Figure a shows relatively mild correlation ($\mathit{\rho }$ $=$ 0.21) compared with Fig. b, which shows strong correlation ($\mathit{\rho }$ $=$ 0.91). This differential behavior can be explained heuristically by considering the differing data-generating values. For the process in Fig. a the short-memory and long-memory components exhibit their effects at opposite ends of the spectrum; see Fig. a. The resulting ARFIMA spectrum, with peaks at either end, makes it easy to distinguish between short and long-memory effects, and consequently the posteriors of $d$ and $\mathit{\varphi }$ are largely uncorrelated. In contrast, the parameters of the process in Fig. b express their behavior at the same end of the spectrum. With negative $d$ these effects partially cancel each other out, except very near the origin where the negative memory effect dominates; see Fig. b. Distinguishing between the effects of $\mathit{\varphi }$ and $d$ is much more difficult in this case; consequently the posteriors are much more dependent.

Marginal posterior density of $d$ from series in Fig. , (a, b). Solid line is density obtained using reversible-jump algorithm. Dotted line is density obtained using fixed $p$ $=$ 1 and $q$ $=$ 0. The true values are $d_l$ $=$ 0.25 and $-$0.35, respectively. (c, d) Shows the posterior densities for $\mathit{\mu }$ and $\mathit{\sigma }$, respectively, corresponding to the series in Fig. a; those for Fig. b look similar. The true values are $\mathit{\mu }$ $=$ 0 and $\mathit{\sigma }$ $=$ 1. True values are marked by an X.

[Figure omitted. See PDF]

In cases where there is significant correlation between $d$ and $\mathit{\varphi }$, it arguably makes little sense to consider only the marginal posterior distribution of $d$. For example the 95 % credibility interval for $d$ from Fig. b is ($-$0.473, $-$0.247), and the corresponding interval for $\mathit{\varphi }$ is ($-$0.910, $-$0.753), yet these clearly give a rather pessimistic view of our joint knowledge about $d$ and $\mathit{\varphi }$; see Fig. c. In theory an ellipsoidal credibility set could be constructed although this is clearly less practical when $\sim$ $\mathit{\omega }$ $>$ 2.

Unknown form: the RJ scheme outlined in Sect.  works well for data simulated with $p$ and $q$ up to 3. The marginal posteriors for $d$ are generally roughly centered around $d_I$ (the data-generating value) and the modal posterior model probability is usually the correct one. To illustrate, consider again the two example data-generating contexts used above.

For both series, kernel density for the marginal posterior for $d$ are plotted in Fig. a and b, together with the equivalent density estimated assuming unknown model orders.

Notice how the densities obtained via the RJ method are very close to those obtained assuming $p$ $=$ 1 and $q$ $=$ 0. The former are slightly more heavy tailed, reflecting a greater level of uncertainty about $d$. Interestingly, the corresponding plots for the posteriors of $\mathit{\mu }$ and $\mathit{\sigma }$ do not appear to exhibit this effect; see Fig. c and d. The posterior model probabilities are presented in Table , showing that the correct modes are being picked up consistently.

Posterior model probabilities for time series from Figs. a, b and a, b for the autoregressive parameter $p$ and moving average parameter $q$. Bold numbers denote the true model.

$p\backslash q$	0	1	2	3	4	5	Marginal
(a)
0	0.000	0.000	0.000	0.000	0.000	0.000	0.000
1	0.805	0.101	0.003	0.000	0.000	0.000	0.908
2	0.038	0.043	0.001	0.000	0.000	0.000	0.082
3	0.005	0.004	0.000	0.000	0.000	0.000	0.009
4	0.000	0.001	0.000	0.000	0.000	0.000	0.001
5	0.000	0.000	0.000	0.000	0.000	0.000	0.000
Marginal	0.848	0.148	0.004	0.000	0.000	0.000
(b)
0	0.000	0.000	0.000	0.000	0.000	0.000	0.000
1	0.829	0.125	0.002	0.000	0.000	0.000	0.956
2	0.031	0.013	0.000	0.000	0.000	0.000	0.044
3	0.000	0.000	0.000	0.000	0.000	0.000	0.000
4	0.000	0.000	0.000	0.000	0.000	0.000	0.000
5	0.000	0.000	0.000	0.000	0.000	0.000	0.000
Marginal	0.860	0.138	0.002	0.000	0.000	0.000

Marginal posterior densities (a) $d$, (b) $\mathit{\alpha }$ from the model Eq. ().

[Figure omitted. See PDF]

As a test of the robustness of the method, consider a complicated short-memory input combined with a heavy-tailed $\mathit{\alpha }$-stable innovations distribution. Specifically, the time series that will be used is the following ARFIMA(2,$d$,1) process $$\begin{array}{cc}\left(\mathrm{1},-,{\displaystyle \frac{\mathrm{9}}{16}},{\mathcal{B}}^{\mathrm{2}}\right)(\mathrm{1}-\mathcal{B}{)}^{0.25}{X}_t=& \left(\mathrm{1},+,{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}},\mathcal{B}\right){\mathit{\epsilon }}_t,\\ & & \text{where}{\mathit{\epsilon }}_t\sim {\mathcal{S}}_{1.75,\mathrm{0}}.\end{array}$$ For more details, see Sect. 7.1. The marginal posterior densities of $d$ and $\mathit{\alpha }$ are presented in Fig. .

Performance looks good despite the complicated structure. The posterior estimate for $d$ is ${\widehat{d}}^{\left(\text{B}\right)}$ $=$ 0.22, with 95 % CI (0.04, 0.41). Although this interval is admittedly rather wide, it is reasonably clear that long memory is present in the signal. The corresponding interval for $\mathit{\alpha }$ is (1.71, 1.88) with estimate ${\widehat{\mathit{\alpha }}}^{\left(\text{B}\right)}$ $=$ 1.79. Finally, we see from Table  that the algorithm is very rarely in the wrong model.

Observational data analysis

We conclude with the application of our method to two long data sets: the Nile water level minima data and the CET. The Nile data are part of the R package “longmemo” and the CET time series can be downloaded from http://www.metoffice.gov.uk/hadobs/hadcet/.

Posterior model probabilities based on simulations of model Eq. () for the autoregressive parameter $p$ and moving average parameter $q$. Bold numbers denote the true model.

$p\backslash q$	0	1	2	3	4	5	Marginal
0	0.000	0.000	0.000	0.000	0.000	0.000	0.000
1	0.000	0.000	0.000	0.000	0.000	0.000	0.000
2	0.000	0.822	0.098	0.001	0.000	0.000	0.921
3	0.014	0.056	0.004	0.000	0.000	0.000	0.075
4	0.003	0.001	0.000	0.000	0.000	0.000	0.004
5	0.000	0.000	0.000	0.000	0.000	0.000	0.000
Marginal	0.017	0.880	0.102	0.002	0.000	0.000

Annual Nile minima time series.

[Figure omitted. See PDF]

Marginal posterior densities for Nile minima; (a) $d$, (b) $\mathit{\mu }$.

[Figure omitted. See PDF]

The Nile data

Because of the fundamental importance of the Nile river to the civilizations it has supported, local rulers kept measurements of the annual maximal and minimal heights obtained by the river at certain points (called gauges). The longest uninterrupted sequence of recordings is from the Roda gauge (near Cairo), between AD 622 and 1284 ($n$ $=$ 663).

There is evidence e.g., that the sequence is not actually homogeneous.

These data are plotted in Fig. .

We immediately observe the apparent low frequency component of the data. The data appear to be on the “verge” of being stationary; however, the general consensus amongst the statistical community is that the series is stationary. The posterior summary statistics are presented in Table , density estimates of the marginal posteriors of $d$ and $\mathit{\mu }$ are presented in Fig. , and the posterior model probabilities are presented in Table .

Table: summary posterior statistics for Nile minima. Plots: marginal posterior densities for Nile minima – (a) $d$, (b) $\mathit{\mu }$.

[Figure omitted. See PDF]

Posterior model probabilities for Nile minima time series for the autoregressive parameter $p$ and moving average parameter $q$. Bold numbers denote the best fit model.

$p\backslash q$	0	1	2	3	4	5	Marginal
0	0.638	0.101	0.010	0.000	0.000	0.000	0.750
1	0.097	0.124	0.011	0.000	0.000	0.000	0.232
2	0.007	0.010	0.000	0.000	0.000	0.000	0.018
3	0.000	0.000	0.000	0.000	0.000	0.000	0.000
4	0.000	0.000	0.000	0.000	0.000	0.000	0.000
5	0.000	0.000	0.000	0.000	0.000	0.000	0.000
Marginal	0.742	0.236	0.022	0.000	0.000	0.000

Summary posterior statistics for Nile minima time series for the long-memory parameter $d$, mean $\mathit{\mu }$, and noise variance $\mathit{\sigma }$.

	Mean	SD	95 % CI end points
$d$	0.402	0.039	0.336	0.482
$\mathit{\mu }$	1158	62	1037	1284
$\mathit{\sigma }$	70.15	1.91	66.46	73.97

The posterior summary statistics and marginal densities of $d$ and $\mathit{\mu }$ for the Nile data are presented in Fig. . Posterior model probabilities are presented in Table . We see that the model with the highest posterior probability is the ARFIMA(0,$d$,0) model with $d$ $\approx$ 0.4. This suggests a strong, pure, long-memory feature. Our results compare favorably with other studies .

It is interesting to compare these findings with other literature. used a semi-parametric Bayesian method on the first 512 observations of the sequence and obtained an estimate for $d$ of 0.278. used a similar method to to estimate $d$ (within an ARFIMA(0,$d$,0) model) at 0.416 with an approximate credibility interval of (0.315, 0.463). similarly found using wavelets $\widehat{d^{\text{B}}}$ $=$ 0.379 with a credibility interval of (0.327, 0.427). obtained $\widehat{d^{\text{B}}}$ $=$ 0.420. obtained $\widehat{d^{\text{B}}}$ $=$ 0.387 with a credibility interval of (0.316, 0.475) using their Bayesian FEXP method.

We note that the interpretation as persistence of the $d$ $\approx$ 0.4 ($H$ $\approx$ 0.9) value that we and others have obtained has been challenged by . In his view the analysis should be applied to the increments of the level heights rather than the level heights themselves, giving an anti-persistent time series with a negative $d$ value. The need for a short-range-dependent component that he argues for is, however, automatically included in the use of an ARFIMA model. Although ARFIMA was originally introduced in econometrics as a phenomenological model of LM, very recent progress is being made in statistics and physics on building a bridge between it and continuous time linear dynamical systems (see e.g., ).

CET time series (deseasonalized).

[Figure omitted. See PDF]

CET time series; (a) assumed deterministic seasonal component $S\left(t\right)$, (b) spectrum of deseasonalized index.

[Figure omitted. See PDF]

Posterior model probabilities for Nile minima time series for the autoregressive parameter $p$ and moving average parameter $q$.

$p\backslash q$	0	1	2	3	4	5	Marginal
0	0.638	0.101	0.010	0.000	0.000	0.000	0.750
1	0.097	0.124	0.011	0.000	0.000	0.000	0.232
2	0.007	0.010	0.000	0.000	0.000	0.000	0.018
3	0.000	0.000	0.000	0.000	0.000	0.000	0.000
4	0.000	0.000	0.000	0.000	0.000	0.000	0.000
5	0.000	0.000	0.000	0.000	0.000	0.000	0.000
Marginal	0.742	0.236	0.022	0.000	0.000	0.000

In conclusion, our findings agree with all published Bayesian long-memory results (except for the anomalous finding of ). Moreover, these findings agree with numerous classical methods of analysis e.g., that have found the best model fit is an ARFIMA(0,$d$,0) model with $d$ $\approx$ 0.4. We note that it is a result of our data analysis method that short-memory can be neglected, rather than being an a priori assumption.

Central England temperature

There is increasing evidence that surface air temperatures posses long memory but long time series are needed to get robust results. The CET index is a famous measure of the monthly mean temperature in an area of southern-central England dating back to 1659 . Given to a precision of 0.5 ${}^{\circ }$C prior to 1699 and 0.1 ${}^{\circ }$C thereafter, the index is considered to be the longest reliable known temperature record from station data. As expected, the CET exhibits a significant seasonal signal, at least some of which must be considered as deterministic. Following the approach of , the index is first deseasonalized using the additive “STL” method . This deseasonalized CET index is shown in Fig. .

The estimated seasonal function $S\left(t\right)$ that was removed is shown in Fig. a. The spectrum of the deseasonalized process is shown in Fig. b. $D$ denotes the seasonal long-memory parameter. Notice that, in addition to the obvious spectral peak at the origin, there still remains a noticeable peak at the monthly frequency $\mathit{\omega }$ $=$ $\frac{\mathit{\pi }}{\mathrm{6}}$. However, there are no further peaks in the spectrum, which would appear to rule out a seasonal ARFIMA (SARFIMA) model. These observations therefore suggest that a simple two-frequency Gegenbauer($d$,$D$;$\frac{\mathit{\pi }}{\mathrm{6}}{)}_{\mathrm{2}}$ process might be an appropriate model. See Appendix  for more details about seasonal long memory.

Applying this model, the marginal posterior statistics are presented in Table  and the joint posterior samples of ($d$, $D$) from this model are plotted in Fig. . These clearly indicate that both $d$ and $D$ are non-zero (albeit small in the case of $D$) suggesting the presence of long memory in both the conventional and seasonal sense.

Joint posterior samples of ($d$, $D$) with 95 % credibility set in red for CET time series.

[Figure omitted. See PDF]

CET time series; posterior estimate (solid line) and 95 % credibility interval (dotted line) for four blocks (black) and whole index (red) for (a) $d$, (b) $\mathit{\mu }$.

[Figure omitted. See PDF]

Posterior summary statistics for CET index for the long-memory parameter $d$, seasonal long-memory parameter $D$, mean $\mathit{\mu }$, and noise variance $\mathit{\sigma }$.

	Mean	SD	95 % CI end points
$d$	0.209	0.013	0.186	0.235
$D$	0.040	0.011	0.018	0.062
$\mathit{\mu }$	9.266	0.144	9.010	9.576
$\mathit{\sigma }$	1.322	0.015	1.294	1.353

In order to compare these results with other publications', it is important to note that to remove annual seasonality from the CET, the series of annual means is often used instead of the monthly series. This of course reduces the fidelity of the analysis. found (using rather crude estimation procedures) that the best-fitting model for the annual means of the CET was the ARFIMA(1,0.33,0) model with $\mathit{\varphi }$ $=$ 0.16. used the same series as test data for their Bayesian method; they fitted each of the ARFIMA models with $p$, $q$ $\le$ 1 and found that all models were suitable. Their estimates of $d$ ranged from 0.24 for $p$ $=$ $q$ $=$ 0 to 0.34 for $p$ $=$ 0, $q$ $=$ 1.

Of course all these studies assume the time series is stationary and, in particular, has a constant mean. The validity of this assumption was considered by who used formal hypothesis testing to consider models: $$Y_t={\mathit{\beta }}_{\mathrm{0}}+{\mathit{\beta }}_{\mathrm{1}}t+X_t,$$ where $\mathit{\left\{}{X}_t\mathit{\right\}}$ is an ARFIMA(0,$d$,0) process. For values of $d$ $=$ 0, 0.05, 0.10, 0.15, ${\mathit{\beta }}_{\mathrm{1}}$ was found to be significantly non-zero (at about 0.23 ${}^{\circ }$C per century) but for $d$ $\ge$ 0.20, statistical significance was not found. later extended this work by replacing the ARFIMA(0,$d$,0) process in Eq. () with a Gegenbauer($d$;$\mathit{\omega }$) process to obtain similar results. However, choice of $\mathit{\omega }$ was rather ad hoc, likely influencing the results.

In order to consider the stationarity of the time series, we divided the series up into four blocks of length 1024 months (chosen to maximize efficiency of the fast Fourier transform) and analyzed each block independently. The posterior statistics for each block are presented in Table  with some results presented graphically in Fig. .

Posterior summary statistics for four blocks of CET index for the long-memory parameter $d$, seasonal long-memory parameter $D$, mean $\mathit{\mu }$, and noise variance $\mathit{\sigma }$.

		Mean	SD	95 % CI end points
1659–1744	$d$	0.277	0.026	0.231	0.332
	$D$	0.054	0.022	0.013	0.097
	$\mathit{\mu }$	9.036	0.347	8.332	9.702
	$\mathit{\sigma }$	1.217	0.027	1.167	1.271
1744–1829	$d$	0.204	0.028	0.151	0.259
	$D$	0.017	0.023	$-$0.028	0.063
	$\mathit{\mu }$	9.107	0.216	8.671	9.533
	$\mathit{\sigma }$	1.348	0.031	1.290	1.409
1829–1914	$d$	0.172	0.027	0.118	0.223
	$D$	0.036	0.022	$-$0.010	0.076
	$\mathit{\mu }$	9.172	0.168	8.859	9.517
	$\mathit{\sigma }$	1.364	0.030	1.312	1.429
1914–2000	$d$	0.163	0.027	0.108	0.213
	$D$	0.063	0.022	0.023	0.109
	$\mathit{\mu }$	9.591	0.152	9.314	9.906
	$\mathit{\sigma }$	1.348	0.030	1.291	1.406

It is interesting to note that the degree of (conventional) long memory is roughly constant over the last three blocks but appears to be larger in the first block. Of particular concern is that there is no value of $d$ that is included in all four 95 % credibility intervals; this would suggest non-stationarity. Although this phenomenon may indeed have a physical explanation, it is more likely caused by the inhomogeneity of the time series. Recall that the first 50 years of the index are only given to an accuracy of 0.5 ${}^{\circ }$C compared to 0.1 ${}^{\circ }$C afterwards; this lack of resolution clearly has the potential to bias in favor of strong autocorrelation when compared with later periods.

Interestingly, the seasonal long-memory parameter $D$ has 95 % credibility intervals that include zero for the both the second and third blocks. Finally, note that the 95 % credibility intervals for $\mathit{\mu }$ all include the range (9.314, 9.517), in other words it is entirely credible that the mean is non-varying over the time period.

Conclusions

We have provided a systematic treatment of efficient Bayesian inference for ARFIMA models, the most popular parametric model combining long- and short-memory effects. Through a mixture of theoretical and empirical work we have demonstrated that our method can handle the sorts of time series data with possible long memory that we are typically confronted with.

Many of the choices made throughout, but in particular those leading to our likelihood approximation, stem from a need to accommodate further extension. For example, in future work we intend to extend them to cope with heavy-tailed innovation distributions. For more evidence of potential in this context, see Sect. 7.

Finally, an advantage of the Bayesian approach is that it provides a natural mechanism for dealing with missing data, via data augmentation. This is particularly relevant for long historical time series, which may, for a myriad of reasons, have recording gaps. For example, some of the data recorded at other gauges along the Nile have missing observations although otherwise span a similarly long time frame. For a demonstration of how this might fit within our framework, see Sect. 5.6 of .

ARFIMA model

We define an autocovariance ACV $\mathit{\gamma }(\cdot )$ of a weakly stationary process as $\mathit{\gamma }\left(k\right)$ $=$ Cov($X_t$ $\left(E\right)\left[X\right]\left)\right(X$ $-$ $\left(E\right)\left[X\right]\left)\right]$ is the lag-covariance matrix. The (normalized) ACF $\mathit{\rho }(\cdot )$ is defined as $\mathit{\rho }\left(k\right)$ $=$ $\frac{\mathit{\gamma }\left(k\right)}{\mathit{\gamma }\left(\mathrm{0}\right)}$. A stationary process $\mathit{\left\{}{X}_t\mathit{\right\}}$ is said to be causal if there exists a sequence of coefficients $\mathit{\left\{}{\mathit{\psi }}_k\mathit{\right\}}$, with finite total mean square $\sum _{k=\mathrm{0}}^{\mathrm{\infty }}{\mathit{\psi }}_k^{\mathrm{2}}$ $<$ $\mathrm{\infty }$ such that for all $t$, a given member of the process can be expanded as a power series in the backshift operator acting on the innovations, $\mathit{\left\{}{\mathit{\epsilon }}_t\mathit{\right\}}$: $$X_t=\mathrm{\Psi }\left(\mathcal{B}\right){\mathit{\epsilon }}_t,\text{where}\mathrm{\Psi }\left(z\right)=\sum _{k=\mathrm{0}}^{\mathrm{\infty }}{\mathit{\psi }}_k{z}^k.$$ The innovations are a white (i.e., stationary, zero mean, iid) noise process with variance ${\mathit{\sigma }}^{\mathrm{2}}$. Causality specifies that for every $t$, $X_t$ can only depend on the past and present values of the innovations $\mathit{\left\{}{\mathit{\epsilon }}_t\mathit{\right\}}$.

A process $\mathit{\left\{}{X}_t\mathit{\right\}}$ is said to be an autoregressive process of order $p$, AR($p$), if for all $t$: $$\begin{array}{cc}\mathrm{\Phi }\left(\mathcal{B}\right)X_t=& {\mathit{\epsilon }}_t,\text{where}\mathrm{\Phi }\left(z\right)=\mathrm{1}+\sum _{k=\mathrm{1}}^p{\mathit{\varphi }}_k{z}^k,\\ & & \text{and}\left({\mathit{\varphi }}_{\mathrm{1}},,,\mathrm{\dots },,,{\mathit{\varphi }}_p\right)\in {\mathbb{R}}^p.\end{array}$$ AR($p$) processes are invertible, stationary and causal if and only if $\mathrm{\Phi }\left(z\right)$ $\ne$ 0 for all $z$ $\in$ $\mathbb{C}$ such that $\left|z\right|$ $\le$ 1. $\mathit{\left\{}{X}_t\mathit{\right\}}$ is said to be a moving average process on the order of $q$, MA($q$), if $$\begin{array}{cc}{X}_t=& \mathrm{\Theta }\left(\mathcal{B}\right){\mathit{\epsilon }}_t,\text{where}\mathrm{\Theta }\left(z\right)=\mathrm{1}+\sum _{k=\mathrm{1}}^q{\mathit{\theta }}_k{z}^k,\\ & & \text{and}\left({\mathit{\theta }}_{\mathrm{1}},,,\mathrm{\dots },,,{\mathit{\theta }}_p\right)\in {\mathbb{R}}^q,\end{array}$$ for all $t$.

Many authors define $\mathrm{\Phi }\left(z\right)$ $=$ 1 $-$ $\sum {\mathit{\varphi }}_k{z}^k$. Our version emphasizes connections between $\mathrm{\Phi }$ and Eqs. ()–().

MA($q$) processes are stationary and causal, and are invertible if and only if $\mathrm{\Theta }\left(z\right)$ $\ne$ 0 for all $z$ $\in$ $\mathbb{C}$ such that $\left|z\right|$ $\le$ 1. A natural extension of the AR and MA classes arises by combining them .

The process $\mathit{\left\{}{X}_t\mathit{\right\}}$ is said to be an ARMA process of orders $p$ and $q$, ARMA($p$,$q$), if for all $t$: $$\mathrm{\Phi }\left(\mathcal{B}\right)X_t=\mathrm{\Theta }\left(\mathcal{B}\right){\mathit{\epsilon }}_t.$$ Although there is no simple closed form for the ACV of an ARMA process with arbitrary $p$ and $q$, so long as the process is causal and invertible, then $\left|\mathit{\rho }\right(k\left)\right|$ $\le$ $C{r}^k$, for $k$ $>$ 0; i.e., it decays exponentially fast. In other words, although correlation between nearby points may be high, dependence between distant points is negligible.

Before turning to long memory, we require one further result. Under some extra conditions, stationary processes with ACV $\mathit{\gamma }(\cdot )$ possess a SDF $f_{\text{sd}}(\cdot )$ defined such that $\mathit{\gamma }\left(k\right)$ $=$ $\underset{-\mathit{\pi }}{\overset{\mathit{\pi }}{\int }}{e}^{ik\mathit{\lambda }}{f}_{\text{sd}}\left(\mathit{\lambda }\right)$ $d\mathit{\lambda }$, $\forall k$ $\in$ $\mathbb{Z}$. This can be inverted to obtain an explicit expression for the SDF e.g.,§4.3: $f_{\text{sd}}\left(\mathit{\lambda }\right)$ $=$ $\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}\sum _{k=-\mathrm{\infty }}^{\mathrm{\infty }}\mathit{\gamma }\left(k\right)e^{-ik\mathit{\lambda }}$, where $-\mathit{\pi }$ $\le$ $\mathit{\lambda }$ $\le$ $\mathit{\pi }$.

Since ACV of a stationary process is an even function of lag, the above equation implies that the associated SDF is an even function. One therefore only needs to be interested positive arguments: 0 $\le$ $\mathit{\lambda }$ $\le$ $\mathit{\pi }$.

Finally, the SDF of an ARMA process is $$f_{\text{sd}}\left(\mathit{\lambda }\right)={\displaystyle \frac{{\mathit{\sigma }}^{\mathrm{2}}}{\mathrm{2}\mathit{\pi }}}{\displaystyle \frac{|\mathrm{\Theta }\left(e^{-i\mathit{\lambda }}\right){|}^{\mathrm{2}}}{|\mathrm{\Phi }\left(e^{-i\mathit{\lambda }}\right){|}^{\mathrm{2}}}},\mathrm{0}\le \mathit{\lambda }\le \mathit{\pi }.$$

For an ARFIMA process (Eq. ) the restriction $\left|d\right|$ $<$ $\frac{\mathrm{1}}{\mathrm{2}}$ is necessary to ensure stationarity; clearly if $\left|d\right|$ $\ge$ $\frac{\mathrm{1}}{\mathrm{2}}$ the ACF would not decay. The continuity between stationary and non-stationary processes around $\left|d\right|$ $=$ $\frac{\mathrm{1}}{\mathrm{2}}$ is similar to those that occur for the AR(1) process with $\left|{\mathit{\varphi }}_{\mathrm{1}}\right|$ $\to$ 1 (such processes are stationary for $\left|{\mathit{\varphi }}_{\mathrm{1}}\right|$ $<$ 1, but the case $\left|{\mathit{\varphi }}_{\mathrm{1}}\right|$ $=$ 1 is the non-stationary random walk).

There are a number of alternative definitions of LM, one of which is particularly useful, as it considers the frequency domain: a stationary process has long memory when its SDF follows $f_{\text{sd}}\left(\mathit{\lambda }\right)$ $\sim$ $c_f{\mathit{\lambda }}^{-\mathrm{2}d}$, as $\mathit{\lambda }$ $\to$ 0${}^{+}$ for some positive constant $c_f$, and where 0 $<$ $d$ $<$ $\frac{\mathrm{1}}{\mathrm{2}}$.

The simplest way of creating a process that exhibits long memory is through the SDF. Consider $f_{\text{sd}}\left(\mathit{\lambda }\right)$ $=$ $|$1 $-$ $e^{i\mathit{\lambda }}{|}^{-\mathrm{2}d}$, where 0 $<$ $\left|d\right|$ $<$ $\frac{\mathrm{1}}{\mathrm{2}}$. By simple algebraic manipulation, this is equivalently $f_{\text{sd}}\left(\mathit{\lambda }\right)$ $=$ ${\left(\mathrm{2},,\sin ,,\frac{\mathit{\lambda }}{\mathrm{2}}\right)}^{-\mathrm{2}d}$, from which we deduce that $f\left(\mathit{\lambda }\right)$ $\sim$ ${\mathit{\lambda }}^{-\mathrm{2}d}$ as $\mathit{\lambda }$ $\to$ 0${}^{+}$. Therefore, assuming stationarity, the process that has this SDF (or any scalar multiple of it) is a long-memory process. More generally, a process having spectral density $$f_{\text{sd}}\left(\mathit{\lambda }\right)={\displaystyle \frac{{\mathit{\sigma }}^{\mathrm{2}}}{\mathrm{2}\mathit{\pi }}}{\left|\mathrm{1},-,e^{i\mathit{\lambda }}\right|}^{-\mathrm{2}d},\mathrm{0}<\mathit{\lambda }\le \mathit{\pi },$$ is called fractionally integrated with memory parameter $d$, Fractionally Integrated FI($d$) with memory parameter $d$ . The full trichotomy of negative, short, and long memory is determined solely by $d$.

In practice this model is of limited appeal to time series analysts because the entire memory structure is determined by just one parameter, $d$. One often therefore generalizes it by taking any short-memory SDF $f_{\text{sd}}^{*}(\cdot )$, and defining a new SDF: $f_{\text{sd}}\left(\mathit{\lambda }\right)$ $=$ $f_{\text{sd}}^{*}\left(\mathit{\lambda }\right){\left|\mathrm{1},,-,,e^{i\mathit{\lambda }}\right|}^{-\mathrm{2}d}$, 0 $\le$ $\mathit{\lambda }$ $\le$ $\mathit{\pi }$. An obvious class of short-memory processes to use this way is ARMA. Taking $f^{*}$ from Eq. () yields so-called autoregressive fractionally integrated moving average process with parameter $d$, and orders $p$ and $q$ (ARFIMA($p$,$d$,$q$)), having SDF: $$f\left(\mathit{\lambda }\right)={\displaystyle \frac{{\mathit{\sigma }}^{\mathrm{2}}}{\mathrm{2}\mathit{\pi }}}{\displaystyle \frac{|\mathrm{\Theta }\left(e^{-i\mathit{\lambda }}\right){|}^{\mathrm{2}}}{\left|\mathrm{\Phi }\right(e^{-i\mathit{\lambda }}){|}^{\mathrm{2}}}}|\mathrm{1}-e^{i\mathit{\lambda }}{|}^{-\mathrm{2}d},\mathrm{0}\le \mathit{\lambda }\le \mathit{\pi }.$$ Choosing $p$ $=$ $q$ $=$ 0 recovers FI($d$) $\equiv$ ARFIMA(0,$d$,0).

Practical utility from the perspective of (Bayesian) inference demands finding a representation in the temporal domain. To obtain this, consider the operator (1 $-$ $\mathcal{B}{)}^d$ for real $d$ $>$ $-$1, which is formally defined using the generalized form of the binomial expansion Eq. 13.2.2: $$\begin{array}{cc}(\mathrm{1}-\mathcal{B}{)}^d=:& \sum _{k=\mathrm{0}}^{\mathrm{\infty }}{\mathit{\pi }}_k^{\left(d\right)}{\mathcal{B}}^k,\text{where}{\mathit{\pi }}_k^{\left(d\right)}=(-\mathrm{1}{)}^k\\ & & {\displaystyle \frac{\mathrm{1}}{\mathrm{\Gamma }(k+\mathrm{1})}}{\displaystyle \frac{\mathrm{\Gamma }(d+\mathrm{1})}{\mathrm{\Gamma }(d-k+\mathrm{1})}}.\end{array}$$ From this observation, one can show that $X_t$ $=$ (1 $-$ $\mathcal{B}{)}^{-d}{Z}_t$, where $\mathit{\left\{}{Z}_t\mathit{\right\}}$ is an ARMA process, has SDF (Eq. ). The operator (1 $-$ $\mathcal{B}{)}^d$ is called the fractional differencing operator since it allows a degree of differencing between the zeroth and first order. The process $\mathit{\left\{}{X}_t\mathit{\right\}}$ is fractionally inverse differenced; i.e., it is an integrated process. The operator is used to re-define both the ARFIMA(0,$d$,0) and more general ARFIMA($p$,$d$,$q$) processes in the time domain. A process $\mathit{\left\{}{X}_t\mathit{\right\}}$ is an ARFIMA(0,$d$,0) process if for all $t$: (1 $-$ $\mathcal{B}{)}^d{X}_t$ $=$ ${\mathit{\epsilon }}_t$. Likewise, a process $\mathit{\left\{}{X}_t\mathit{\right\}}$ is an ARFIMA($p$,$d$,$q$) process if for all $t$: $\mathrm{\Phi }\left(\mathcal{B}\right)$(1 $-$ $\mathcal{B}{)}^d{X}_t$ $=$ $\mathrm{\Theta }\left(\mathcal{B}\right){\mathit{\epsilon }}_t$, where $\mathrm{\Phi }$ and $\mathrm{\Theta }$ are given in Eqs. () and (), respectively.

Finally, to connect back to our first definition of long memory, consider the ACV of the ARFIMA(0,$d$,0) process. By using the definition of spectral density to directly integrate Eq. (), and an alternative expression for ${\mathit{\pi }}_k^{\left(d\right)}$ in Eq. () $${\mathit{\pi }}_k^{\left(d\right)}={\displaystyle \frac{\mathrm{1}}{\mathrm{\Gamma }(k+\mathrm{1})}}{\displaystyle \frac{\mathrm{\Gamma }(k-d)}{\mathrm{\Gamma }(-d)}},$$ one can obtain the following representation of the ACV of the ARFIMA(0,$d$,0) process: $${\mathit{\gamma }}_d(k;\mathit{\sigma })={\mathit{\sigma }}^{\mathrm{2}}{\displaystyle \frac{\mathrm{\Gamma }(\mathrm{1}-\mathrm{2}d)}{\mathrm{\Gamma }(\mathrm{1}-d)\mathrm{\Gamma }\left(d\right)}}{\displaystyle \frac{\mathrm{\Gamma }(k+d)}{\mathrm{\Gamma }(\mathrm{1}+k-d)}}.$$ Because the parameter ${\mathit{\sigma }}^{\mathrm{2}}$ is just a scalar multiplier, we may simplify notation by defining ${\mathit{\gamma }}_d\left(k\right)$ $=$ ${\mathit{\gamma }}_d(k;\mathit{\sigma })/{\mathit{\sigma }}^{\mathrm{2}}$, whereby ${\mathit{\gamma }}_d(\cdot )$ $\equiv$ ${\mathit{\gamma }}_d(\cdot$; 1). Then the ACF is $${\mathit{\rho }}_d\left(k\right)={\displaystyle \frac{\mathrm{\Gamma }(\mathrm{1}-d)}{\mathrm{\Gamma }\left(d\right)}}{\displaystyle \frac{\mathrm{\Gamma }(k+d)}{\mathrm{\Gamma }(\mathrm{1}+k-d)}},$$ from which Stirling's approximation gives ${\mathit{\rho }}_d\left(k\right)$ $\sim$ $\frac{\mathrm{\Gamma }(\mathrm{1}-d)}{\mathrm{\Gamma }\left(d\right)}{k}^{\mathrm{2}d-\mathrm{1}}$, confirming a power-law relationship for the ACF. Finally, note that Eq. () can be used to represent ARFIMA(0,$d$,0) as an AR($\mathrm{\infty }$) process, as $X_t$ $+$ $\sum _{k=\mathrm{1}}^{\mathrm{\infty }}{\mathit{\pi }}_k^{\left(d\right)}{X}_{t-k}$ $=$ ${\mathit{\epsilon }}_t$. Furthermore, noting that in this case ${\mathit{\psi }}_k^{\left(d\right)}$ $=$ ${\mathit{\pi }}_k^{(-d)}$ leads to the following MA($\mathrm{\infty }$) analog: $X_t$ $=$ $\sum _{k=\mathrm{0}}^{\mathrm{\infty }}\frac{\mathrm{1}}{\mathrm{\Gamma }(k+\mathrm{1})}\frac{\mathrm{\Gamma }(k+d)}{\mathrm{\Gamma }\left(d\right)}{\mathit{\epsilon }}_{t-k}$.

Seasonal long-memory models

We define a seasonal differencing operator (1 $-$ ${\mathcal{B}}^{\text{s}}$), as a natural extension to a SARFIMA processes by combining seasonal and non-seasonal fractional differencing operators : $$(\mathrm{1}-\mathcal{B}{)}^d{\left(\mathrm{1},-,{\mathcal{B}}^{\text{s}}\right)}^D{X}_t={\mathit{\epsilon }}_t.$$

The generalization to include both seasonal and non-seasonal short-memory components is obvious : $$\begin{array}{cc}{\mathrm{\Phi }}^{\left(p\right)}\left(\mathcal{B}\right){\mathrm{\Phi }}_{\text{s}}^{\left(P\right)}& \left({\mathcal{B}}^{\text{s}}\right)(\mathrm{1}-\mathcal{B}{)}^d{\left(\mathrm{1},-,{\mathcal{B}}^{\text{s}}\right)}^D{X}_t\\ & ={\mathrm{\Theta }}^{\left(q\right)}\left(\mathcal{B}\right){\mathrm{\Theta }}_{\text{s}}^{\left(Q\right)}\left({\mathcal{B}}^{\text{s}}\right){\mathit{\epsilon }}_t.\end{array}$$

Focusing on the first of these issues, considered generalising the ARFIMA(0,$d$,0) process in a different manner by retaining only one pole but at any given frequency in [0, $\mathit{\pi }$]. The model he suggested was later studied and popularized by and , and became known as the “Gegenbauer process”.

A process $\mathit{\left\{}{X}_t\mathit{\right\}}$ is a Gegenbauer ($d$; $\mathit{\omega }$) process if for all $t$: $${\left(\mathrm{1},-,\mathrm{2},u,\mathcal{B},+,{\mathcal{B}}^{\mathrm{2}}\right)}^d{X}_t={\mathit{\epsilon }}_t,$$ where $\mathit{\omega }$ $=$ ${\cos }^{-\mathrm{1}}u$ is called the Gegenbauer frequency. The obvious extension to include short-memory components ${\mathrm{\Phi }}^{\left(p\right)}$ and ${\mathrm{\Theta }}^{\left(p\right)}$ is denoted GARMA($p$,$d$,$q$;$\mathit{\omega }$).

The term “Gegenbauer” derives from the close relationship to the Gegenbauer polynomials, a set of orthogonal polynomials useful in applied mathematics. The Gegenbauer polynomials are most usefully defined in terms of their generating function. The Gegenbauer polynomial on the order of $k$ with parameter $d$, $G_k^{\left(d\right)}(\cdot )$ satisfies $${\left(\mathrm{1},-,\mathrm{2},u,z,+,z^{\mathrm{2}}\right)}^{-d}\equiv \sum _{k=\mathrm{0}}^{\mathrm{\infty }}{G}_k^{\left(d\right)}\left(u\right)z^k.$$

The spectral density function of the Gegenbauer($d$;$\mathit{\omega }$) process is $$f\left(\mathit{\lambda }\right)={\displaystyle \frac{{\mathit{\sigma }}^{\mathrm{2}}}{\mathrm{2}\mathit{\pi }}}{\left|\mathrm{2},(,\cos ,\mathit{\lambda },-,\cos ,\mathit{\omega },)\right|}^{-\mathrm{2}d},\mathrm{0}\le \mathit{\lambda }\le \mathit{\pi }.$$

Note that Gegenbauer($d$;$\mathit{\omega }$) processes possess a pole at the Gegenbauer frequency $\mathit{\omega }$. Gegenbauer processes may be considered to be somewhat ambiguous in terms of long memory. Non-trivial (i.e., $\mathit{\omega }$ $\ne$ 0) Gegenbauer processes have bounded spectral density functions at the origin, and therefore do not have long memory according to our strict definition. Consequently a more general Gegenbauer process was developed: let $\mathbf{d}$ $=$ ($d_{\mathrm{1}}$, …, $d_k$) and $\mathit{\omega }$ $=$ (${\mathit{\omega }}_{\mathrm{1}}$, …, ${\mathit{\omega }}_k$), and for all $j$, $u_j$ $=$ $\cos {\mathit{\omega }}_j$ (assumed distinct). Then a process $\mathit{\left\{}{X}_t\mathit{\right\}}$ is a $k$-factor Gegenbauer($\mathbf{d}$;$\mathit{\omega }$) process if for all $t$ : $$\prod _{j=\mathrm{1}}^k{\left(\mathrm{1},-,\mathrm{2},u_j,\mathcal{B},+,{\mathcal{B}}^{\mathrm{2}}\right)}^{d_j}{X}_t={\mathit{\epsilon }}_t.$$

The spectral density function of the $k$-factor Gegenbauer($\mathbf{d}$;$\mathit{\omega }$) process is $$f\left(\mathit{\lambda }\right)={\displaystyle \frac{{\mathit{\sigma }}^{\mathrm{2}}}{\mathrm{2}\mathit{\pi }}}\prod _{j=\mathrm{1}}^k{\left|\mathrm{2},\left(\cos ,\mathit{\lambda },-,\cos ,{\mathit{\omega }}_j\right)\right|}^{-\mathrm{2}{d}_j},\mathrm{0}\le \mathit{\lambda }\le \mathit{\pi }.$$

Indeed, $k$-factor Gegenbauer models are very flexible, and include nearly all other seasonal variants of ARFIMA processes such as the flexible-seasonal ARFIMA and fractional ARUMA processes. Importantly, they also includes SARFIMA processes : a SARFIMA(0,$d$,0) $\times$ (0,$D$,0)${}_{\text{s}}$ process is equivalent to a $⌊\frac{s+\mathrm{2}}{\mathrm{2}}⌋$ factor Gegenbauer($\mathbf{d}$;$\mathit{\omega }$) process where: $${\mathit{\omega }}_j={\displaystyle \frac{\mathrm{2}\mathit{\pi }(j-\mathrm{1})}{s}},j=\mathrm{1},\mathrm{\dots },k,$$ and $d_{\mathrm{1}}$ $=$ $\frac{d+D}{\mathrm{2}}$, $d_j$ $=$ $D$ for $j$ $=$ 2, …, $k$, unless $s$ is even in which case $d_k$ $=$ $\frac{D}{\mathrm{2}}$.

Although $k$-factor Gegenbauer models are very general, one particular sub-model is potentially very appealing. This is the two-factor model, with one pole at the origin and one at a non-zero frequency. In order to conform with notation for ARFIMA(0,$d$,0) processes, we will slightly re-define this model: a process $\mathit{\left\{}{X}_t\mathit{\right\}}$ is a simple two-frequency Gegenbauer process with parameters $d$, $D$, and $\mathit{\omega }$, denoted Gegenbauer($d$,$D$;$\mathit{\omega }{)}_{\mathrm{2}}$ if for all $t$: $${\left(\mathrm{1},-,\left(,\mathrm{2},\cos ,\mathit{\omega },\right),\mathcal{B},+,{\mathcal{B}}^{\mathrm{2}}\right)}^D(\mathrm{1}-\mathcal{B}{)}^d{X}_t={\mathit{\epsilon }}_t.$$ The Bayesian MCMC methodology developed here is easily extended to incorporate these seasonal fractional models. It is assumed that the frequency $\mathit{\omega }$, or seasonal period $s$, is a priori known.

Acknowledgements

We thank one anonymous reviewer and M. Crucifix for their comments, which helped to improve this manuscript. C. L. E. Franzke is supported by the German Research Foundation (DFG) through the cluster of excellence CliSAP (EXC177), N. W. Watkins is supported by ONR NICOP grant N62909-15-1-N143, and both are supported by the Norwegian Research Council KLIMAFORSK project 229754. N. W. Watkins thanks the University of Potsdam for hospitality. Edited by: Z. Toth Reviewed by: M. Crucifix and another anonymous referee