More accurate, calibrated bootstrap confidence intervals for estimating the correlation between two time series

Abstract

Estimation of Pearson’s correlation coefficient between two time series, in the evaluation of the influences of one time-dependent variable on another, is an often used statistical method in climate sciences. Data properties common to climate time series, namely non-normal distributional shape, serial correlation, and small data sizes, call for advanced, robust methods to estimate accurate confidence intervals to support the correlation point estimate. Bootstrap confidence intervals are estimated in the Fortran 90 program PearsonT (Mudelsee, Math Geol 35(6):651–665, 2003), where the main intention is to obtain accurate confidence intervals for correlation coefficients between two time series by taking the serial dependence of the data-generating process into account. However, Monte Carlo experiments show that the coverage accuracy of the confidence intervals for smaller data sizes can be substantially improved. In the present paper, the existing program is adapted into a new version, called PearsonT3, by calibrating the confidence interval to increase the coverage accuracy. Calibration is a bootstrap resampling technique that performs a second bootstrap loop (it resamples from the bootstrap resamples). It offers, like the non-calibrated bootstrap confidence intervals, robustness against the data distribution. Pairwise moving block bootstrap resampling is used to preserve the serial dependence of both time series. The calibration is applied to standard error-based bootstrap Student’s $t$ confidence intervals. The performance of the calibrated confidence interval is examined with Monte Carlo simulations and compared with the performance of confidence intervals without calibration. The coverage accuracy is evidently better for the calibrated confidence intervals where the coverage error is acceptably small already (i.e., within a few percentage points) for data sizes as small as 20.

Appendices

Appendix A: Software

The calibration method explained in the paper was adapted into the Fortran 90 software PearsonT3. The software is freely available at http://www.climate-risk-analysis.com. The installation requires copying the PearsonT3 executable file into an appropriate directory and installing the free graphic program Gnuplot (http://www.gnuplot.info/). The gnuplot executable file, gnuplot.exe, needs to be in the same directory as PearsonT3.exe. The software is command line driven and can be run from the Windows command prompt or simply by double clicking the executable file. After starting PearsonT3 the program asks for a name and path of input data file. The data file should be a simple text file and in the format

$$\begin{aligned} \begin{array}{l@{\quad }l@{\quad }l} t_{1} &{} x_{1} &{} y_{1} \\ t_{2} &{} x_{2} &{} y_{2} \\ \vdots &{} \vdots &{} \vdots \\ t_{n} &{} x_{n} &{} y_{n}, \end{array} \end{aligned}$$

where $t$ is sampling times and $x$ and $y$ are two equally long time series with data size $\ge 10$. The time series are automatically mean detrended (the mean of the data is subtracted from the data). A linear detrend option was included in the old version of the software but is not included in PearsonT3. If the time series samples contain linear or more complex trend, we recommend some detrending prior to the analysis in PearsonT3 to fulfill the weakly stationary assumptions.

The persistence times are estimated with the least-squares algorithm TAUEST Mudelsee (2002) with automatic bias correction. If the bias-corrected equivalent autocorrelation coefficient becomes $>1$ (Eq. 6), then the bias correction is not performed, which can occur if $n$ is small and the autocorrelation coefficient is large. After the estimation the time series are plotted up on the screen along with an $x-y$ scatterplot to test for the linear relationship. The results are printed on the screen, which informs about the data file name, the time interval $[t(1);\,t(n)]$, the number of data points ($n$), the persistence times ($\tau _{X}$ and $\tau _{Y}$), and the estimated correlation coefficient ($r_{XY}$) with 95  % calibrated confidence interval. The results are also written into a result file, along with the data, means, and mean detrended data. The final result file, named PearsonT3.dat, is a plain ASCII file, which is saved in the same directory as the executable file PearsonT3.exe.

Appendix B: Bivariate AR(1) process

The bivariate AR(1) process for uneven spacing is given by ((Mudelsee 2010, Ch. 7.6))

$$\begin{aligned} X(1)&= \fancyscript{E}_{\mathrm{N}(0,\, 1)}^X(1),\nonumber \\ Y(1)&= \fancyscript{E} _{\mathrm{N}(0,\, 1)}^Y(1),\nonumber \\ X(i)&= \exp \left\{ -\left[ T(i) - T(i-1) \right] / \tau _X \right\} \cdot X(i-1)\nonumber \\&+\, \fancyscript{E}_{\mathrm{N}(0,\, 1-\exp \left\{ -2\left[ T(i) - T(i-1) \right] / \tau _X \right\} )}^X(i), \qquad i = 2,\ldots ,n,\nonumber \\ Y(i)&= \exp \left\{ -\left[ T(i) - T(i-1) \right] / \tau _Y \right\} \cdot Y(i-1)\nonumber \\&+\, \fancyscript{E}_{\mathrm{N}(0,\, 1-\exp \left\{ -2\left[ T(i) - T(i-1) \right] / \tau _Y \right\} )}^Y(i), \qquad i = 2,\ldots ,n, \end{aligned}$$

(12)

where the white-noise terms are correlated as

$$\begin{aligned} \textit{CORR}\left[ \fancyscript{E}_{\mathrm{N}(0, 1)}^X(1), \fancyscript{E}_{\mathrm{N}(0, 1)}^Y(1)\right] =\rho _\fancyscript{E}, \end{aligned}$$

(13)

$$\begin{aligned}&\textit{CORR}\Bigl [ \fancyscript{E}_{\mathrm{N}(0,\, 1-\exp \left\{ -2\left[ T(i) - T(i-1) \right] / \tau _X \right\} )}^X(i), \fancyscript{E}_{\mathrm{N}(0,\, 1-\exp \left\{ -2\left[ T(j) - T(j-1) \right] / \tau _Y \right\} )}^Y(j) \Bigr ] = 0,\\&\qquad i, j = 2,\ldots ,n,\qquad i \ne j, \end{aligned}$$

$$\begin{aligned}&\textit{CORR}\Bigl [ \fancyscript{E}_{\mathrm{N}(0,\, 1-\exp \left\{ -2\left[ T(i) - T(i-1) \right] / \tau _X \right\} )}^X(i), \fancyscript{E}_{\mathrm{N}(0,\, 1)}^Y(1)\Bigr ] = 0,\\&\qquad \qquad \, i = 2,\ldots ,n, \end{aligned}$$

$$\begin{aligned}&\textit{CORR}\Bigl [\fancyscript{E}_{\mathrm{N}(0,\, 1)}^X(1), \fancyscript{E}_{\mathrm{N}(0,\, 1-\exp \left\{ -2\left[ T(i) - T(i-1) \right] / \tau _Y \right\} )}^Y(i), \Bigr ] = 0,\\&\qquad \qquad \, i = 2,\ldots ,n. \end{aligned}$$

The process is strictly stationary with the properties

$$\begin{aligned}&\hbox {E}[X(i)]= \hbox {E}[Y(i)]= 0,\end{aligned}$$

(14)

$$\begin{aligned}&\hbox {VAR}[X(i)]=\hbox {VAR}[Y(i)]=1,\end{aligned}$$

(15)

$$\begin{aligned}&\hbox {CORR}[X(i),Y(i)]=\rho _{XY}= \rho _\fancyscript{E}. \end{aligned}$$

(16)