What to do if time series data remains autocorrelated?

Question

Data:

I have 92 years of monthly climate data. One of my variables is a drought index (SPEI) ranging from -2 (dry) to 2 (wet). All the data can be found here.

Data Structure:

year month  order.ID  season  temp.avg  ppt.avg    GDD    pdo     spei
1923     0        13       1    6.1612  0.23516  73.54   0.27   0.6544
1923     1        14       1    5.0967  0.24161  40.66   0.01   0.4837
1923     2        15       1     3.425  0.24428  47.19  -0.95   0.5207
1923     3        16       2    9.7870  0.46612  161.3  -1.23   0.2631
1923     4        17       2    12.753    0.304  238.1  -0.64   0.2476

note: PDO = Pacific Decadal Oscillation Index | order.ID = (monthly) time series order

Graph of SPEI across time:

I am trying to determine the long-term (92yr) trends of the SPEI data. I am doing so with a general least squares approach using gls function from the nlme package in R.

---

Issue:

After examining my best (lowest AIC) model ( gls(spei ~ I(year - 1950) + pdo) ) using ACF and PACF, it was very clear that there was underlying temporal autocorrelation and likely periodicity.

I've attempted to account for these issues in two ways:

1. Accounting for periodicity by including sin + cos parameters to the model. I did so by examining how including the following parameters,

   I(cos(2*pi/P*(order.ID))) + I(sin(2*pi/P*(order.ID)))
   

   in my model using various period lengths (P) affected the AIC of the model:

   

   • Though including the sin/cos parameters using the period (630) that created the model with the lowest AIC improved the AIC of my model, it did not correct for the issues I was seeing in my ACF/PACF plots. In fact, even including upwards of 4 combinations of SIN/COS parameters using numerous "strong" periods (e.g., 630, 238, 183 and 49) failed to eliminate the ACF/PACF issues:

   

2. Accounting for autocorrelation using corARMA argument in gls function. Specifically, I tried ARMA models using all combinations of p = {1:4} & q = {0:3} with no luck in eliminating the autocorrelation.

   • The best performing (lowest AIC) model was

     gls(spei ~ I(year - 1950) + pdo, correlation = corARMA(p=4,q=2), method = ML)
     

   • But the resulting ACF/PACF still had issues:

   

Note: I did try a combination of 1 and 2 (using sin+cos parameters and a corARMA argument in my model), but even the combination of approaches failed to "fix" the problem.

Additional note:

I tried adding months and seasons (separately and together) to the model in two ways: 1st, as periods (12 and 4 respectively) in SIN/COS parameters and, second, as categorical dummy variables in the model. Both approaches resulted in models with higher AICs - essentially suggesting the periodicity is not due to months or seasons??

---

Question:

What do I do now? How do I go about "fixing" this data - i.e., how do I properly account for the periodicity and autocorrelation so that I can accurately view the long-term annual (linear) trend (and correct p value of the coefficient)?

Obviously nothing I've tried so far has worked...

Answer 1

[EDIT: I think I'm in over my head here and would appreciate input. I originally specified a ARMA(2,1) with seasonal ARMA(1,1), assuming that any trend that would normally require differencing would be handled by the year xreg variable. But the ar1 coefficient was 1.5, which isn't good. So I've now switched to ARIMA(1,1,1) with seasonal ARIMA(0,0,1), but I'm not sure that's right either. If you look at How to know if a time series is stationary or non-stationary?, StasK makes the point that non-economic time series (climate in particular) are more complex than ARIMA can handle -- I think.]

I believe I've calculated what you're looking for, using R's arima to account for a yearly seasonality and also an ARMA (2, 1):

climate <- read.csv ("http://theforestecologist.web.unc.edu/files/2015/07/df.clim_.csv")

climate$yr <- (climate$year - 1950) + (climate$month / 12)

spei <- ts (climate$spei, start=c(1923, 1), frequency=12)
pdo  <- ts (climate$pdo,  start=c(1923, 1), frequency=12)
yr   <- ts (climate$yr,   start=c(1923, 1), frequency=12)


spei2 <- window (spei, start=c(1950, 1))
pdo2  <- window (pdo,  start=c(1950, 1))
yr2   <- window (yr,   start=c(1950, 1))

mod1 <- arima (spei2, order=c(1, 1, 1), seasonal=c(0, 0, 1),  xreg=yr2)
mod1

tsdiag (mod1, gof.lag=24)

The results are:

Coefficients:
         ar1      ma1     sma1     yr2
      0.5285  -0.3086  -0.6463  0.0292
s.e.  0.0984   0.1089   0.0255  0.0433

sigma^2 estimated as 0.03562:  log likelihood = 190.25,  aic = -370.5

I also tried with your pdo variable in, but it doesn't help:

mod2 <- arima (spei2, order=c(1, 1, 1), seasonal=c(0, 0, 1), xreg=cbind (yr2, pdo2))
mod2

tsdiag (mod2, gof.lag=24)

Does this look reasonable?

At any rate, the coefficient for year (the linear trend) is not significant with an s.e. nearly as large as the coefficient. (Basically the same if we look at 1950-2014 or 1923-2014.) Hopefully this gives you a start.

arima is a little backwards, syntax-wise, from a regression, but has flexible ARIMA and Seasonal ARIMA terms.

Answer 2

Did you try using Dicky Fuller or Augmented Dicky Fuller? I'm not sure if there is something odd about climate data, but it works for things like air pollution and inflation, once you've tried everything else keep going until you don't have a unit root