How to analyze trend in non-periodic time series

Question

Suppose I have following non-periodic time series. Obviously the trend is decreasing and I would like to prove it by some test (with p-value). I am unable to use classic linear regression due to strong temporal (serial) auto-correlation among values.

library(forecast)
my.ts <- ts(c(10,11,11.5,10,10.1,9,11,10,8,9,9,
               6,5,5,4,3,3,2,1,2,4,4,2,1,1,0.5,1),
            start = 1, end = 27,frequency = 1)
plot(my.ts, col = "black", type = "p",
     pch = 20, cex = 1.2, ylim = c(0,13))
# line of moving averages 
lines(ma(my.ts,3),col="red", lty = 2, lwd = 2)

What are my options?

Answer 1

As you said, the trend in your example data is obvious. If you want just to justify this fact by hypothesis test, than besides using linear regression (the obvious parametric choice), you can use non-parametric Mann-Kendall test for monotonic trend. The test is used to

assess if there is a monotonic upward or downward trend of the variable of interest over time. A monotonic upward (downward) trend means that the variable consistently increases (decreases) through time, but the trend may or may not be linear. (http://vsp.pnnl.gov/help/Vsample/Design_Trend_Mann_Kendall.htm)

moreover, as noted by Gilbert (1987), the test

is particularly useful since missing values are allowed and the data need not conform to any particular distribution

The test statistic is the difference between negative and positive $x_j-x_i$ differences among all the $n(n-1)/2$ possible pairs, i.e.

$$ S = \displaystyle\sum_{i=1}^{n-1}\displaystyle\sum_{j=i+1}^{n}\mathrm{sgn}(x_j-x_i) $$

where $\mathrm{sgn}(\cdot)$ is a sign function. $S$ can be used to calculate $\tau$ statistics that is similar to correlation as it ranges from $-1$ to $+1$, where the sign suggests negative, or positive trend and value of $\tau$ is proportional to slope of the trend.

$$ \tau = \frac{S}{n(n-1)/2} $$

Finally, you can compute $p$-values. For samples of size $n \le 10$ you can use tables of precomputed $p$-values for different values of $S$ and different sample sizes (see Gilbert, 1987). With larger samples, first you need to compute variance of $S$

$$ \mathrm{var}(S) = \frac{1}{18}\Big[n(n-1)(2n+5) - \displaystyle\sum_{p=1}^{g}t_p(t_p-1)(2t_p+5)\Big] $$

and then compute $Z_{MK}$ test statistic

$$ Z_{MK} = \begin{cases} \frac{S-1}{\mathrm{var}(S)} & \text{if} ~ S > 0 \\ 0 & \text{if} ~ S = 0 \\ \frac{S+1}{\mathrm{var}(S)} & \text{if} ~ S < 0 \end{cases} $$

the value of $Z_{MK}$ is compared to standard normal values

• $Z_{MK} \ge Z_{1-\alpha}$ for upward trend,
• $Z_{MK} \le -Z_{1-\alpha}$ for downward trend,
• $|Z_{MK}| \ge Z_{1-\alpha/2}$ for upward or downward trend.

In this thread you can find R code implementing this test.

Since the $S$ statistic is compared to all possible pairs of observations then, instead of using normal approximation for $p$-value you can use permutation test that is obvious for this case. First, you compute $S$ statistic from your data and then you randomly shuffle your data multiple times and compute it for each of the samples. $p$ is simply the proportion of cases when $S_\text{data} \ge S_\text{permutation}$ for upward trend or $S_\text{data} \le S_\text{permutation}$ for downward trend.

---

Gilbert, R.O. (1987). Statistical Methods for Environmental Pollution Monitoring. Wiley, NY.

Önöz, B., & Bayazit, M. (2003). The power of statistical tests for trend detection. Turkish Journal of Engineering and Environmental Sciences, 27(4), 247-251.

Answer 2

You can use OLS because there is no serial autocorrelation (at least in the sample you supplied); note the Durbin-Watson test statistic of 1.966 (≈2).

So, the significantly negative coefficient estimate for x1 is all you need to say something like

The observed count of [certain species] is decreasing by about 1,000 per year.

or

The observed count of [certain species] is decreasing by between 628 and 1,408 per year (at 95% confidence).

This assumes that the methodology for counting the species has good coverage and is consistent over the years in your sample.

This was produced with this Python code (sorry; don't have R handy):

import numpy as np
import statsmodels.api as sm

y = [10,12,10,11,8,9,6,4,2,4]
x = np.arange(len(y))
x = sm.add_constant(x)

mod = sm.OLS(y, x)
result = mod.fit()
print(result.summary())

Answer 3

The problem that you have "I am unable to use classic linear regression due to strong temporal (serial) auto-correlation among values." is in reality an opportunity. I took your 27 values and used AUTOBOX a piece of software (which I have helped develop) that can (optionally) automatically determine a possible model. Here is the actual/fit and forecast graph . The ACF of the residuals is here with residual plot here . The model is here and here and here . Two coefficients aptly describe the data with estimated "trend" aka "drift" i.e. period to period differential of -.596. Note that this is one kind of trend where your model used the counting numbers 1,2,...27 as a predictor variable. If your data suggested that kind of trend then the software would have found it to be more applicable. I will try and find an earlier post of mine which fully detailed/contrasted these two kinds of trends. Here Identifying a stochastic trend model and Detecting initial trend or outliers

Answer 4

You can use Spearman's rank correlation coefficient to determine the degree to which you data is monotonic. It returns positive values for monotonic increasing data and negative values for monotonic decreasing data (between -1 and +1). Following the link above, there is also a section dealing significance testing, although I am sure most software packages will have a p-value done for you when computing correlation coefficients (e.g. in Matlab: [RHO,PVAL] = corr(...); in R: cor.test(x,...))

Answer 5

Knowing the source of data would be very helpful, and also the information if the values of my.ts could get negative or not.

However, taking a quick look on the plot, rather than seeing a constant linear trend, I rather suggest that the time series is not stationary, hence integrated. As an example, stock prices are also integrated, but stock returns not anymore (they fluctuate near 0).

This hypothesis can also be tested using the Augmented Dickey Fuller Test:

require(tseries)
adf.test(my.ts)

Augmented Dickey-Fuller Test
Dickey-Fuller = -2.9557, Lag order = 2, p-value = 0.7727
alternative hypothesis: stationary

Given the p-value is not lower 0.05, there is no evidence that the process is stationary.

To get the data stationary, you have to difference it:

diff.ts <- diff(my.ts)
plot(diff.ts)

Now the data shows no trend anymore, and the only thing you will find is an autoregressive term of order 2 (using acf(diff.ts)).