11
$\begingroup$

I am having basically the same issue than in this thread, except one thing:

The difference, in my case, is that my data is measured weekly and not daily, so the argument of a too high seasonality (> 350) does not hold for my data, since the seasonality in my case is 52 (52 weeks in a year).

And yet, when I use auto.arima(), R returns the ARIMA model (p,d,q) = (2,1,1) and (P,D,Q) = (0,0,0), while the seasonal pattern in my data is blatant... How could you explain that R completely dismisses the seasonality in my data?

Since I'm still in a learning phase, I am using the data set cmort available in the astsa library, so everyone here can use the same data as me.

And I have done cmort <- ts(cmort,frequency=52) to be sure that the seasonality in my data is taken account of, but it didn't change anything.

Improve this question
$\endgroup$
7
  • 7
    $\begingroup$ It has been suggested to close this question as off-topic. I disagree and think there is a genuine statistical question in here. $\endgroup$
    – Stephan Kolassa
    May 18, 2016 at 8:12
  • 2
    $\begingroup$ Questions about errors or error messages in R, or questions about how to perform a well-defined operation in R or similar are often better posted in the R tag at StackOverflow. Your question may be migrated there if enough people agree with the close vote. (I disagree and think it's quite on-topic here.) Incidentally, please don't just cross-post your question if you do think it's better at SO; instead, click on the little "flag" link and request migration - a moderator will happily assist you. $\endgroup$
    – Stephan Kolassa
    May 18, 2016 at 8:16
  • 2
    $\begingroup$ 52 "seasons" is a lot, & the seasonality's only blatant because we can see (in @StephanKolassa's plot) a smooth pattern that won't be seen so easily by separately comparing Wk1 for different years, then Wk2 for different years, & so on. You might well get better forecasts over the next year by regressing on a few Fourier terms. $\endgroup$
    – Scortchi - Reinstate Monica ♦
    May 18, 2016 at 16:44
  • 1
    $\begingroup$ @Scortchi Does "regressing on a few Fourier terms" mean to perform a linear regression of $x_t$ (my data) on $a_1 \cos(\omega t) + a_2 \cos(2\omega t) +...$ (where $\omega $ is the pulsation, of course) and find the coefficients $a_1$, $a_2$, ... that minimize the mean square error? $\endgroup$
    – R. Bourgeon
    May 18, 2016 at 19:14
  • 4
    $\begingroup$ There are 365.25/7 weeks in a year, which is not 52. If you have a long time series, taking the frequency as 52 will not yield a satisfactory result. Actually, if you had some 40 years of data (OK, I guess your series is shorter, but think about it anyway), the 52-week seasonality would be off by 1 year. Certainly, that would disturb the results. $\endgroup$
    – Richard Hardy
    May 18, 2016 at 19:36

3 Answers 3

Reset to default
10
$\begingroup$

(First off, cmort is already a ts object with frequency 52, so you don't need to coerce it.)

I'd say seasonality is visible, not that it is blatant:

library(forecast)
library(astsa)
seasonplot(cmort)

seasonplot

Per the help page (?auto.arima), auto.arima() decides whether or not to take seasonal differences by using a OCSB test. It's quite possible that this test simply got it wrong in this instance; it's a statistical test, after all. You can force a seasonal model by setting D=1, although auto.arima() runs for quite some time with forced seasonality. (Note that the information criteria are not comparable between the original and the differenced series.)

Auto-fitted model:

> auto.arima(cmort)
Series: cmort 
ARIMA(2,1,1)                    

Coefficients:
         ar1     ar2      ma1
      0.0957  0.2515  -0.6435
s.e.  0.4302  0.2444   0.4155

sigma^2 estimated as 33.72:  log likelihood=-1609.89
AIC=3227.77   AICc=3227.85   BIC=3244.68

Model with forced seasonality:

> auto.arima(cmort,D=1)
Series: cmort 
ARIMA(0,0,0)(1,1,0)[52] with drift         

Coefficients:
         sar1    drift
      -0.5737  -0.0257
s.e.   0.0378   0.0041

sigma^2 estimated as 47.7:  log likelihood=-1537.6
AIC=3081.21   AICc=3081.26   BIC=3093.57
Improve this answer
$\endgroup$
5
  • $\begingroup$ Thanks Stephan, the forecasting looks indeed way better if we coerce the seasonality by forcing $D=1$ ! $\endgroup$
    – R. Bourgeon
    May 18, 2016 at 8:21
  • $\begingroup$ Stephan, what do you mean when you write "Note that the information criteria are not comparable between the original and the differenced series." ? Does it mean that I cannot use the fact that the AIC and the BIC are lower in the second model to state : "the second model is more accurate than the first one?" $\endgroup$
    – R. Bourgeon
    May 18, 2016 at 10:55
  • 3
    $\begingroup$ Yes, that is exactly what it means. Information criteria are calculated as the negative log-likelihood plus a penalty term that depends on the number of parameters. But the log-likelihood of a model applied to the original data and the log-likelihood of the same model applied to differenced data are not comparable at all, since the data may well be on completely different scales. Therefore, information criteria cannot be compared between undifferenced and differenced series, either. Much better to use holdout samples and compare accuracy of models on those. $\endgroup$
    – Stephan Kolassa
    May 18, 2016 at 11:23
  • $\begingroup$ However, here, we are not in the situation of the same model applied to one set of data and then the same model applied to differenced data... I think that the situation is: we have applied two different SARIMA models (first (2,1,1,0,0,0) and then (0,0,0,1,1,0)) to the same set of data cmort... Haven't we? So the likelihood is computed from the same dataset, isn't it? $\endgroup$
    – R. Bourgeon
    May 18, 2016 at 11:28
  • 1
    $\begingroup$ No, first we apply the model to the original data (yielding the nonseasonal model), and then we apply the same or a different model to the seasonally differenced data (yielding the seasonal model). So we cannot compare two models applied to different data (original vs. seasonally differenced) based on ICs. $\endgroup$
    – Stephan Kolassa
    May 18, 2016 at 11:32
11
$\begingroup$

A glance at the manual for auto.arima shows that an explanation of precisely why it found the solution it did in this case would be complicated: depending on the fitting algorithm (conditional least squares by default); on the details of the stepwise selection procedure, & the criteria used (approximate AICc by default); & on the particular configuration & the sequence of the stationarity tests carried out. Nevertheless, I daresay that ARIMA(2,1,1) is a perfectly good model for one-step-ahead forecasts, because the evident smoothness of the seasonality means that knowing the observed value from exactly one year ago will give you little more information about this week's likely value than just knowing the last couple of weeks' observations.

If you want to forecast over the next year or two, you will need to take seasonality into account. Rather than taking seasonal differences, or looking for seasonal autoregressive or moving average terms, or estimating 52 fixed seasonal effects; you might be better off taking advantage of the smoothness & regressing on a few Fourier terms, as explained by @Glen_b here. (This also makes it easy to allow for non-integer frequencies—as @RichardHardy ponts out might be appropriate in this case (if the extra day or two each year weren't simply dropped).) Your model will then have deterministic, but smoothed, seasonal effects, & you can still allow ARIMA errors:

$$ y_t = \alpha + \sum_{i=1}^k \beta_i \sin \left(i \cdot \frac{2 \pi}{\omega}\cdot t\right) + \sum_{i=1}^k \gamma_i \cos\left(i \cdot \frac{2 \pi}{\omega} \cdot t \right) + \frac{1+\sum_{i=1}^q \theta_i B^i}{(1-B)^d\left(1-\sum_{i=1}^p \phi_i B^i\right)}\cdot\varepsilon_t $$ $$ \varepsilon\sim\mathrm{WN}(\sigma^2) $$ The frequency $\omega$ is assumed known. The parameters you need to estimate are: $\alpha$, the intercept; the $\beta$s & $\gamma$s, the coefficients of the sine & cosine terms terms; the $\phi$s & $\theta$'s, the autoregressive & moving average coefficients; & $\sigma^2$, the variance of the white noise innovations. The orders for the $\mathrm{ARIMA}(p,d,q)$ errors & the number of harmonics $k$ can be chosen heuristically or by using significance tests/information criteria.

The forecast package has a function to create the Fourier terms; here's an illustration:

library(astsa)
library(forecast)
# make a training set & test set
window(cmort, start=c(1970,1), end=c(1977,40)) -> cmort.train
window(cmort, start=c(1977,41), end=c(1979,40)) -> cmort.test
# fit non-seasonal & seasonally-differenced models for comparison
Arima(cmort.train, order=c(2,1,1)) -> mod.ns
Arima(cmort.train, seasonal=c(1,1,0), include.drift=T) -> mod.sd
# make predictors
fourier(cmort, K=10) -> fourier.terms
fourier.terms.train <- fourier.terms[1:length(cmort.train),]
fourier.terms.test <- fourier.terms[(length(cmort.train)+1):length(cmort),]
# plot first 5 harmonics
plot(cmort.train, col="grey")
lines(fitted(Arima(cmort.train, order=c(0,0,0), xreg=fourier.terms.train[,1:2])), col="red")
lines(fitted(Arima(cmort.train, order=c(0,0,0), xreg=fourier.terms.train[,1:4])), col="darkorange")
lines(fitted(Arima(cmort.train, order=c(0,0,0), xreg=fourier.terms.train[,1:6])), col="green")
lines(fitted(Arima(cmort.train, order=c(0,0,0), xreg=fourier.terms.train[,1:8])), col="blue")
lines(fitted(Arima(cmort.train, order=c(0,0,0), xreg=fourier.terms.train[,1:10])), col="violet")

harmonics

Up to the third harmonic seems to be enough, and you still need to difference.

# fit models
Arima(cmort.train, order=c(0,1,0), xreg=(fourier.terms.train[,1:6])) -> mod.ft3.df
acf(residuals(mod.ft3.df), ci.type="ma", lag.max=52)
acf(residuals(mod.ft3.df), type="partial", lag.max=52)
Arima(cmort.train, order=c(2,1,0), xreg=(fourier.terms.train[,1:6])) -> mod.ft3.df.ar2
acf(residuals(mod.ft3.df.ar2), lag.max=52)
acf(residuals(mod.ft3.df.ar2), type="partial", lag.max=52)

Examining correlograms suggests a couple of AR terms will do to leave the residuals looking like white noise. Or you can investigate the number of harmonics to use, as well as the ARMA structure, usng auto.arima:

# auto.arima
auto.arima(cmort.train, xreg=fourier.terms.train[,1:2], seasonal=F, d=1, allowdrift, stepwise=F, approximation=F) -> auto1
auto.arima(cmort.train, xreg=fourier.terms.train[,1:4], seasonal=F, d=1, allowdrift, stepwise=F, approximation=F) -> auto2
auto.arima(cmort.train, xreg=fourier.terms.train[,1:6], seasonal=F, d=1, allowdrift, stepwise=F, approximation=F) -> auto3
auto.arima(cmort.train, xreg=fourier.terms.train[,1:8], seasonal=F, d=1, allowdrift, stepwise=F, approximation=F) -> auto4
auto.arima(cmort.train, xreg=fourier.terms.train[,1:10], seasonal=F, d=1, allowdrift=T, stepwise=F, approximation=F) -> auto5
auto1
auto2
auto3
auto4
auto5

The ARMA(2,1) structure is found by auto.arima's minimizing the AICc for the model with 4 harmonics on differenced observations, & has the lowest AICc overall.

# look at accuracy on test set
plot(cmort.test, col="grey")
lines(forecast(mod.ns,h=length(cmort.test))$mean, col="red")
lines(forecast(mod.sd,h=length(cmort.test))$mean, col="darkorange")
lines(forecast(mod.ft3.df.ar2,h=length(cmort.test), xreg=fourier.terms.test[,1:6])$mean, col="green")
lines(forecast(auto4,h=length(cmort.test), xreg=fourier.terms.test[,1:8])$mean, col="blue")

forecasts

accuracy(forecast(mod.ns,h=length(cmort.test)),cmort.test)
accuracy(forecast(mod.sd,h=length(cmort.test)),cmort.test)
accuracy(forecast(mod.ft3.df.ar2,h=length(cmort.test),xreg=fourier.terms.test[,1:6]),cmort.test)
accuracy(forecast(auto4,h=length(cmort.test),xreg=fourier.terms.test[,1:8]),cmort.test)

Summarizing the root mean square errors shows that both the harmonic models perform better on the test set, with little to choose between them:

$$ \begin{array}{l,l,l} &\text{Training} & \text{Test}\\ \mathrm{ARIMA}(2,1,1) & 5.729 & 7.657\\ \mathrm{SARIMA}(1,1,0)_{52}\text{ with drift} & 6.481 & 7.390\\ \text{3 harmonics, }\mathrm{ARIMA}(2,1,0) & 5.578 & 5.151\\ \text{4 harmonics, }\mathrm{ARIMA}(2,1,1) & 5.219 & 5.188 \end{array} $$

Improve this answer
$\endgroup$
5
  • $\begingroup$ You're welcome. (1) I suspect the SARIMA model is unfairly disadvantaged on the test set by its drift - I doubt anyone can really tell whether this series is drifting, trending, or just wandering - but you get the idea. (2) Using several test sets & training sets could give a better idea of the relative out-of-sample performance of the different model-fitting procedures with this series. (3) I'd also like to show how to let the seasonality evolve in a dynamic linear model, but probably won't get around to it / remember how - in any case this series is probably too short to need that. $\endgroup$
    – Scortchi - Reinstate Monica ♦
    May 19, 2016 at 17:39
  • $\begingroup$ I think that you've reached the end of what I can understand, your suggestions outstand my knowledge in this domain (I'm quite a newbie). I'd be glad to understand them, but at least all I can say is that your answer gave me an excellent insight on what "regressing on a few Fourier terms" means in practice. If I understood correctly, fourier(cmort, K=10) computes the 10 first Fourier coefficients of the vector cmort? And your example with the additional argument xreg in arima made me understand that you regress on other variables in addition of the MA and the AR terms. Thanks. $\endgroup$
    – R. Bourgeon
    May 19, 2016 at 19:03
  • $\begingroup$ Yes, type forecast:::...fourier to see the source code. Most of it's concerned with error-checking & other peripheral matters; there's a small loop that makes the sine & cosine terms. $\endgroup$
    – Scortchi - Reinstate Monica ♦
    May 20, 2016 at 7:58
  • $\begingroup$ Most data has some trend. If you run a simple trend model, there is a strong trend. Even if you take the worst case scenario and just use the last 230 observations the coefficient is negative. $\endgroup$
    – Tom Reilly
    May 20, 2016 at 15:04
  • $\begingroup$ @TomReilly: Well the local-level model seems to do all right, even without a drift term (stochastic trend). But if I hadn't differenced - yes, I'd've needed to include a (deterministic) trend. Perhaps I'll add a little about the two approaches, but the main point was to demonstrate using Fourier terms to model seasonality. My guess would be that the series isn't long enough to distinguish them & it'd come down to which was thought more appropriate based on subject-matter knowledge. (On the other hand the autocorrelation of the differenced process was very negative at the first lag.) $\endgroup$
    – Scortchi - Reinstate Monica ♦
    May 20, 2016 at 15:57
1
$\begingroup$

I say that the seasonality is blatant. ACF Plot

I ran this through Autobox(I am affiliated) and there is also a blatant downward trend in the data and some outliers too.

enter image description here

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.