# How to compare 2 non-stationary time series to determine a correlation?

I have two data series that plot median age at death over time. Both series demonstrate an increased age at death over time, but one much lower than another. I want to determine if the increase in the age at death of the lower sample is significantly different than that of the upper sample.

Here are the data, ordered by year (from 1972 through 2009 inclusive) as rounded to three decimal places:

Cohort A    70.257  70.424  70.650  70.938  71.207  71.263  71.467  71.763  71.982  72.270  72.617  72.798  72.964  73.397  73.518  73.606  73.905  74.343  74.330  74.565  74.558  74.813  74.773  75.178  75.406  75.708  75.900  76.152  76.312  76.558  76.796  77.057  77.125  77.328  77.431  77.656  77.884  77.983
Cohort B    5.139   8.261   6.094   12.353  11.974  11.364  12.639  11.667  14.286  12.794  12.250  14.079  17.917  16.250  17.321  18.182  17.500  20.000  18.824  21.522  21.500  21.167  21.818  22.895  23.214  24.167  26.250  24.375  27.143  24.500  23.676  25.179  24.861  26.875  27.143  27.045  28.500  29.318


Both series are non-stationary - how can I compare the two please? I'm using STATA. Any advice would be gratefully received.

• If you provide a link to your data, Matt, we can edit your question to include those data. – whuber Nov 17 '11 at 16:48
• Many thanks for your interest in my plight - link to data added. Any help would be appreciated.Matt – Matt Hurley Nov 17 '11 at 17:20
• @ Matt: Glancing at your data, it looks like they are both upward trends. So are you essentially interested in the hypothesis that one cohort is increasing more quickly than the other? – Andrew Nov 17 '11 at 17:41
• Yes Andrew - the upper cohort is the general population, while the cohort with the poorer age of death is a group dying of the same condition. The null hypothesis being that if they are closely correlated any improvement in survival is potentially due to common factors (and not improved care of said condition). – Matt Hurley Nov 17 '11 at 17:47
• The increases, however measured, are so obviously different that no formal test is necessary. (You will get p-values of $10^{-10}$ or less almost no matter how you assess and compare the slopes, no matter how you model the variation.) The difference in life expectancies exponentially decreased at a rate of 0.83% per year. The interesting thing is the sudden setback in Cohort B at 2001; this change--equivalent to an instantaneous loss of six years of progress--is statistically significant. – whuber Nov 17 '11 at 18:06

This is a simple situation; let's keep it so. The key is to focus on what matters:

• Obtaining a useful description of the data.

• Assessing individual deviations from that description.

• Assessing the possible role and influence of chance in the interpretation.

• Maintaining intellectual integrity and transparency.

There are still many choices and many forms of analysis will be valid and effective. Let's illustrate one approach here that can be recommended for its adherence to these key principles.

To maintain integrity, let's split the data into halves: the observations from 1972 through 1990 and those from 1991 through 2009 (19 years in each). We will fit models to the first half and then see how well the fits work in projecting the second half. This has the added advantage of detecting significant changes that may have occurred during the second half.

To obtain a useful description, we need to (a) find a way to measure the changes and (b) fit the simplest possible model appropriate for those changes, evaluate it, and iteratively fit more complex ones to accommodate deviations from the simple models.

(a) You have many choices: you can look at the raw data; you can look at their annual differences; you can do the same with the logarithms (to assess relative changes); you can assess years of life lost or relative life expectancy (RLE); or many other things. After some thought, I decided to consider RLE, defined as the ratio of life expectancy in Cohort B relative to that of the (reference) Cohort A. Fortunately, as the graphs show, the life expectancy in Cohort A is increasing regularly in a stable fashion over time, so that most of the random-looking variation in the RLE will be due to changes in Cohort B.

(b) The simplest possible model to start with is a linear trend. Let's see how well it works.

The dark blue points in this plot are the data retained for fitting; the light gold points are the subsequent data, not used for the fit. The black line is the fit, with a slope of .009/year. The dashed lines are prediction intervals for individual future values.

Overall, the fit looks good: examination of residuals (see below) shows no important changes in their sizes over time (during the data period 1972-1990). (There is some indication that they tended to be larger early on, when life expectancies were low. We could handle this complication by sacrificing some simplicity, but the benefits for estimating the trend are unlikely to be great.) There is just the tiniest hint of serial correlation (exhibited by some runs of positive and runs of negative residuals), but clearly this is unimportant. There are no outliers, which would be indicated by points beyond the prediction bands.

The one surprise is that in 2001 the values suddenly fell to the lower prediction band and stayed there: something rather sudden and large happened and persisted.

Here are the residuals, which are the deviations from the description mentioned previously.

Because we want to compare the residuals to 0, vertical lines are drawn to the zero level as a visual aid. Again, the blue points show data used for the fit. The light gold ones are the residuals for data falling near the lower prediction limit, post-2000.

From this figure we can estimate that the effect of the 2000-2001 change was about -0.07. This reflects a sudden drop of 0.07 (7%) of a full lifetime within Cohort B. After that drop, the horizontal pattern of residuals shows that the previous trend continued, but at the new lower level. This part of the analysis should be considered exploratory: it was not specifically planned, but came about due to a surprising comparison between the held-out data (1991-2009) and the fit to the rest of the data.

One other thing--even using just the 19 earliest years of data, the standard error of the slope is small: it's only .0009, just one-tenth of the estimated value of .009. The corresponding t-statistic of 10, with 17 degrees of freedom, is extremely significant (the p-value is less than $10^{-7}$); that is, we can be confident the trend is not due to chance. This is one part of our assessment of the role of chance in the analysis. The other parts are the examinations of the residuals.

There appears to be no reason to fit a more complicated model to these data, at least not for the purpose of estimating whether there's a genuine trend in RLE over time: there is one. We could go further and split the data into pre-2001 values and post-2000 values in order to refine our estimates of the trends, but it wouldn't be completely honest to conduct hypothesis tests. The p-values would be artificially low, because the splitting testing were not planned in advance. But as an exploratory exercise, such estimation is fine. Learn all you can from your data! Just be careful not to deceive yourself with overfitting (which is almost sure to happen if you use more than a half dozen parameters or so or use automated fitting techniques), or data snooping: stay alert to the difference between formal confirmation and informal (but valuable) data exploration.

Let's summarize:

• By selecting an appropriate measure of life expectancy (the RLE), holding out half the data, fitting a simple model, and testing that model against the remaining data, we have established with high confidence that: there was a consistent trend; it has been close to linear over a long period of time; and there was a sudden persistent drop in RLE in 2001.

• Our model is strikingly parsimonious: it requires just two numbers (a slope and intercept) to describe the early data accurately. It needs a third (the date of the break, 2001) to describe an obvious but unexpected departure from this description. There are no outliers relative to this three-parameter description. The model is not going to be substantially improved by characterizing serial correlation (the focus of time-series techniques generally), attempting to describe the small individual deviations (residuals) exhibited, or introducing more complicated fits (such as adding in a quadratic time component or modeling changes in the sizes of the residuals over time).

• The trend has been 0.009 RLE per year. This means that with each passing year, the life expectancy within Cohort B has had 0.009 (almost 1%) of a full expected normal lifetime added to it. Over the course of the study (37 years), that would amount to 37*0.009 = 0.34 = one-third of a full lifetime improvement. The setback in 2001 reduced that gain to about 0.28 of a full lifetime from 1972 to 2009 (even though during that period overall life expectancy increased 10%).

• Although this model could be improved, it would likely need more parameters and the improvement is unlikely to be great (as the near-random behavior of the residuals attests). On whole, then, we should be content to arrive at such a compact, useful, simple description of the data for so little analytical work.

• :whuber the identified one time pulses play no role in the forecast to wit: – IrishStat Nov 18 '11 at 1:26

I think that whuber's answer is straightforward and a simple one for a non-time series person like me to understand. I base mine on his. My answer is in R not Stata as I don't know stata that well.

I wonder if the question is actually asking us to look at whether the absolute year on year increase is the same in the two cohorts (rather than relative). I think this is important and illustrate it as follows. Consider the following toy example:

a <- 21:40
b <- 41:60
x <- 1:20
plot(y = a, x = x, ylim = c(0, 60))
points(y = b, x = x, pch = 2)


Here we have 2 cohorts, each of which have a steady 1 year per year increase in median survival. So each year both cohorts in this example increase by the same absolute amount, but the RLE gives the following:

rle <-  a / b
plot(rle)


Which obviously has an upward trend, and the p value to test the hypothesis that the gradient of the line 0 is 2.2e-16. The fitted straight line (let's ignore that this line looks curved) has a gradient of 0.008. So even though both cohorts have the same absolute increase in a year, the RLE has an upward slope.

So if you use RLE when you want to look for absolute increases, then you'll inappropriately reject the null hypothesis.

Using the supplied data, calculating the absolute difference between the cohorts we get:

Which implies that the absolute difference between median survival is gradually decreasing (i.e. the cohort with the poor survival is gradually getting closer to the cohort with the better survival).

• :Andrew note the two cluster of residuals at the end of your graph. This suggests a potential deficiency in your analysis. Unfortunately even skilled Mathematicians who are also skilled Statisticians sometimes are not time series people . What I suggested is what is standard operating procedure for time series analysis. – IrishStat Nov 18 '11 at 15:10
• @andrew Nice response. I upvoted, I believe in your qualifications! – Adam Nov 19 '11 at 8:38
• :Adam Thanks for your nice words. You will note that I only opine on time series problems/questions where I have some expertise as I have specialized in that area for the last 40 years. – IrishStat Nov 19 '11 at 18:03
• @IrishStat a touch discursive. Nice experience, you've been putting punks like me into place since before I was born. – Adam Nov 20 '11 at 2:20
• :Adam The whole idea is to help other people. I truly hope that has been the case. I am confused by the word "punks" as I don't have that feeling at all. Just trying to help ! – IrishStat Nov 20 '11 at 14:27

These two time series seems to have a deterministic trend. This is one relation that you obviously want to remove before further analysis. Personally, i would proceed as follows:

1) I would run a regression for each time series against a constant and a time, and compute the residual for each time series.

2) Taking the two residuals series, computed in the step above, i would run a simple linear regression (without a constant term) and look at the t-statistic, p-value, and decided on whether or not there is further dependence between the two series.

This analysis assumes the same set of assumption you make in a linear regression.

• :user3544 Running a regression against a constant in time is a form of detrending which is one form of pre-whitening ; differencing is another form of pre-whitening : Both are presumptive as there can be multiple trends or different forms of differencing operators. Note that a differencing operator is a particular case of an ARIMA filter which converts a series to white noise. IN general one wants to filter X to make it noise (x) and then apply that filter to Y to create y ( not necessarily white noise )for purposes of Identifying structure or the transfer between Y and X. – IrishStat Nov 17 '11 at 20:40
• :User3544 I should have applauded your using simple uncomplicated single trends but I think that sometimes one shouldn't be assumptive.A simple trend often is useless if there are level shifts in the series or there are a number of trends. Testing hypothesis regarding constancy of parameters needs to be robust and conducted where one searches for the point in time where the parameters may have changed rather than arbitrarily visually selecting the point using the Chow Test. It is well known that unusual values downward bias tests for auto-correlation which is why one needs to detect them. – IrishStat Nov 18 '11 at 0:49
• IrishStat: your applaud is well received and heard.. :) I totally agree with you comments, however, given the the plots of the two time series, i thought "Let's keep it simple".. :) – Lalas Nov 19 '11 at 3:37
• :user3544 My favorite Einstein quote is "Make everything as simple as possible, but not simpler" or restated Make models as simple as possible, but not simpler because some people think simple is the objective while it can be the cause of insufficient analysis. In this case your suggestion would have sufficed to identify the contemporaneous correlation and the lag correlation between the two cohorts while illuminating the level shift. Read more: brainyquote.com/quotes/quotes/a/…. – IrishStat Nov 19 '11 at 14:00

In some cases one knows a theoretical model which can used to test your hypothesis. In my world tis "knowledge" is often absent and one must resort to statistical techniques that can be classified as exploratory data analysis which summarizes what follows.When analyzing time series data that are non-stationary i.e. has autocorrelative properties simple cross-correlation tests are often misleading insofar as false positives can be easily found. One of the earliest analysis of this is found in Yule, G.U, 1926, "Why do we sometimes get nonsense correlations between time series? A study in sampling and the nature of time series", Journal of the Royal Statistical Society 89, 1–64 . Alternatively when one or more of the series themselves have been effected by exceptional activity ( see whuber "the sudden setback in Cohort B at 2001 ) which can effectively hide significant relationships. Now detecting a relationship between time series extends to examining not only contemporaneous relationships but possible lagged relationships. Continuing, if either series has been effected by anomalies ( one-time events ) then we must robustify our analysis by adjusting for these one-time distortions. The literature of time series points out how to identify the relationship via pre-whitening in order to more clearly identify structure. Pre-whitening adjusts for intra-correlative structure prior to identifying inter-correlative structure. Notice the key word was identifying structure. This approach easily leads to the following "useful model" :

Y(T) = -194.45
+[X1(T)][(+ 1.2396+ 1.6523B** 1)] COHORTA

   +[X2(T)][(- 3.3924)]                :PULSE          3

+[X3(T)][(- 2.4760)]                :LEVEL SHIFT   30 reflecting persistant  unusal activity

+[X4(T)][(+ 1.1453)]                :PULSE         29

+[X5(T)][(- 2.7249)]                :PULSE         11

+[X6(T)][(+ 1.5248)]                :PULSE         27

+[X7(T)][(+ 2.1361)]                :PULSE          4

+[X8(T)][(+ 1.6395)]                :PULSE         13

+[X9(T)][(- 1.6936)]                :PULSE         12

+[X10(T)[(- 1.6996)]                :PULSE         19

+[X11(T)[(- 1.2749)]                :PULSE         10

+[X12(T)[(- 1.2790)]                :PULSE         17

+       [A(T)]


which suggests a contemporary relationship of 1.2936 and a lagged effect of 1.6523. Note that there were a number of years where unusual activity was identified viz. (1975,2001,1983,1999,1976,1985,1984,1991 and 1989). The adjustments for the years allows us to more clearly assess the relationship between these two series.

In terms of making a forecast

MODEL EXPRESSED AS AN XARMAX
Y[t] = a[1]Y[t-1] + ... + a[p]Y[t-p]
+ w[0]X[t-0] + ... + w[r]X[t-r]
+ b[1]a[t-1] + ... + b[q]a[t-q]
+ constant

THE RIGHT-HAND SIDE CONSTANT IS: -194.45

COHORTA 0 1.239589 X( 39 ) * 78.228616 = 96.971340

COHORTA 1 1.652332 X( 38 ) * 77.983000 = 128.853835

I~L00030 0 -2.475963 X( 39 ) * 1.000000 = -2.475963

      NET PREDICTION FOR Y(    39 )=                     28.894826


Four coefficients is all that is required to make a forecast and of course a prediction for CohortA at time period 39 (78.228616) obtained from the ARIMA model for Cohorta.

• Nine out of 38 years exhibit "unusual" activities? In a model with (apparently) 25 parameters?! Something's wrong with that interpretation. Apart from all these pulses and the level shift, have you found the nonlinear component yet in the trend of cohort b? – whuber Nov 17 '11 at 20:50
• What are the numbers in the right column (3, 30, 29, 11, etc.)? They appear to be part of your description of the data and therefore are parameters as well. Even if we don't count them, using 14 parameters to describe 38 values, especially when the question is simply "is there a trend?", seems excessive. For that matter, what exactly is the trend? Where in all these estimates does one dig it out? If a physician were to approach you and ask, "ok, what has been accomplished for the patients in Cohort B since 1972," could you tell them in one clear sentence? – whuber Nov 17 '11 at 22:13
• Re "very powerful prediction": I may be misunderstanding what your model amounts to, but in general, better descriptions of a method (time series or otherwise) that identifies a quarter of one's data as "unusual" and requiring "adjustment" would be "overfitted" and "needlessly complex." Your assertion of no trend in Cohort B is simply incredible. – whuber Nov 17 '11 at 22:27
• @Adam, because this analysis essentially ignores the information about variability contained in the 10 "pulses," any prediction bands it places around forecasts will be highly over-optimistic (too tight). Moreover, a deeper analysis that includes all the data (unlike my illustrative analysis that includes only the first half) will detect a nonlinear component consistent with a slight decline in the trend, and that too is not detected here. Of greater importance than forecasting is understanding the 2000-2001 effect: if it could be repeated, then all forecasts are likely wrong. – whuber Nov 18 '11 at 15:11
• @whuber I admittedly don't know all the technical jargon, but your explanation makes a good deal of sense. Thanks a lot. – Adam Nov 19 '11 at 1:45