raw data stationary but still can see trend and seaonality is stl

Question

So I am looking at unit sales data. I am doing a univariate time series analysis. My data is weekly sales numbers figures, spanning 2012- 2014 (obviously no till end 2014). I first ploted my response var against time. I suspected trend and seasonality which would make sense. So I did a seasonal trend decompose with loess (STL). Here is the result. However when I did a stationarity test on my original dataset, it gives me stationarity. And I tried Box-Ljung Philips-Perron unit root test augmented dicky fuller test with stationary and explosive as alternative both and kpss. Stationary. So why is it that my trend and sesonality not significant? I assigned to every observation a month index and did boxplot(sales~month, data=???) Here is the result: Then asked myself if differences in mean were significant, they were not. I tried aov pairwise.t.test and TuckeyHD. even if difference in mean for some means was important, overall not. Could I say then that no significant trend and seasonality appears in the dataset and model with an ARIMA? Also I have only 156 observations in my dataset. Here is the acf and pacf what ARMA would you specify? if ARMA at all. I am really at loss it is the first time I see a staionary raw sales figure.

Answer 1

This is what I get using the R package forecast for ARIMA model selection and tsoutliers for detection of possible additive outliers, level shifts or temporary changes.

These are the data:

x <- structure(c(1291, 2497, 2643, 2633, 2637, 3246, 2764, 2552, 2190, 2423, 2393, 2278, 1573, 792, 1826, 1819, 1610, 1344, 1762, 1917, 1920, 1804, 1936, 1912, 2000, 1742, 2120, 2483, 2475, 2252, 2061, 2145, 1901, 2223, 1934, 1439, 1963, 1879, 1587, 2467, 2852, 3811, 3154, 2428, 2881, 2625, 2416, 2538, 1947, 2526, 2421, 1277, 1334, 
2009, 2693, 2624, 2714, 2882, 3044, 2847, 2439, 2690, 2677, 2461, 1495, 2501, 2614, 2486, 2390, 1901, 2315, 2710, 2436, 2184, 2079, 2234, 2382, 2213, 2506, 2673, 2579, 2341, 2246, 2331, 1670, 2263, 2453, 2411, 2653, 2754, 2818, 2635, 2577, 2897, 2796, 2039, 2188, 
2749, 2790, 2510, 2128, 2839, 2270, 1791, 1373, 2078, 2878, 3198, 3095, 2964, 3230, 3122, 2809, 2902, 3048, 2984, 2927, 2740, 2732, 1572, 2829, 1831, 2821, 2907, 2973, 2690, 2511, 2632, 2502, 2590, 2907, 3046, 2953, 2663, 2613, 2637, 2061, 2812, 2660, 2653, 2708,  2921, 2862, 2526, 2863, 3002, 3227, 2254, 2606, 2959, 2963, 3302), 
.Tsp = c(2012, 2014.90384615385, 52), class = "ts")

Fit an ARIMA model:

require(tsoutliers)
res <- tso(x, args.tsmethod = list(ic="bic"))
res
# output
# Series: x 
# ARIMA(2,0,0) with non-zero mean 
# Coefficients:
#          ar1      ar2  intercept       TC6       AO14      LS40       AO65
#       0.6172  -0.2038  1960.4325  1024.713  -885.4593  494.7448  -932.6690
# s.e.  0.0836   0.0846    83.7532   296.303   252.3412  102.8566   250.3566
#          LS107       AO120      AO122
#       370.7175  -1219.8245  -987.1700
# s.e.   95.7657    252.4097   252.4363    
# sigma^2 estimated as 88411:  log likelihood=-1081.49
# AIC=2184.99   AICc=2186.87   BIC=2218.25    
# Outliers:
#   type ind    time coefhat  tstat
# 1   TC   6 2012:06  1024.7  3.458
# 2   AO  14 2012:14  -885.5 -3.509
# 3   LS  40 2012:40   494.7  4.810
# 4   AO  65 2013:13  -932.7 -3.725
# 5   LS 107 2014:03   370.7  3.871
# 6   AO 120 2014:16 -1219.8 -4.833
# 7   AO 122 2014:18  -987.2 -3.911

The data are fitted by a stationary AR(2) model with some outliers summarized in the output. According to this model, a cyclic behaviour rather than seasonality is observed in the data.

These are the residuals from the fitted model.

-562.694; 861.332; 215.01; 360.635; 400.555; -19.663; 63.132; 89.367; -159.542; 310.001; 101.865; 
80.454; -540.365; -10.707; -60.819; -131.584; -301.359; -436.549; 105.306; -50.281; -56.658; -142.147; 
62.602; -66.119; 63.857; -253.158; 302.135; 379.365; 224.416; 80.363; 25.394; 181.858; -152.887; 
336.83; -200.61; -451.631; 318.987; -189.269; -322.656; 225.7; 313.434; 1113.322; -57.1; -182.213; 
584.986; -98.521; -57.223; 141.604; -567.276; 401.332; -181.432; -1142.653; -401.002; 5.72; 284.743; 
-68.867; 203.088; 301.483; 378.136; 115.386; -138.022; 324.644; 73.6; -83.233; 14.096; 63.987; 124.937;
-57.862; -51.839; -507.671; 188.566; 228.418; -205.009; -207.419; -212.721; -44.265; -13.322; -242.081;
185.378; 137.111; -0.256; -146.214; -113.48; -18.343; -751.16; 267.111; -43.557; -81.991; 224.644; 
167.73; 218.705; 16.786; 84.769; 403.277; 92.964; -536.499; 59.121; 373.918; 99.044; -91.952; -292.79; 
596.918; -488.727; -471.684; -709.996; 155.382; 64.387; 263.095; 50.069; 47.84; 373.703; 74.842; 
-117.303; 146.866; 171.693; 36.535; 48.783; -116.079; -20.282; 6.377; 5.001; -16.583; 0.506; 82.552; 
96.052; -210.158; -201.05; -27.24; -268.39; -75.503; 160.697; 121.984; 7.788; -196.493; -86.462; 
-90.693; -691.693; 419.689; -313.173; -73.34; -44.991; 132.638; -46.613; -302.799; 229.549; 92.099; 
299.978; -783.563; 214.792; 152.29; 10.151; 418.609

$\quad$

plot(residuals(res$fit))

Autocorrelations of the residuals:

par(mfrow = c(2,1), mar = c(5,4,4,3))
tmp <- acf(residuals(res$fit), lag.max = 53, plot = FALSE)
tmp$acf[1] <- NA
plot(tmp, ylim = c(-0.4, 0.4))
pacf(residuals(res$fit), lag.max = 53, ylim = c(-0.4, 0.4))

The autocorrelations of the residuals reveal the significance of the lag of order 52 (one year lag). One way to incorporate this lag is to define it as follows:

xlag52 <- ts(c(rep(NA, 52), window(x, start = c(2013))))
tsp(xlag52) <- tsp(x)

and pass xlag52 through argument xreg in function tso or arima. However, probably because the sample is small for such lag, the algorithm run into difficulties when fitting an ARIMA model.

Edit

The idea of seasonal dummies mentioned by @IrishStat is much more appropriate than trying to add a lag of order 52. You can use these functions of the forecast package: seasonaldummies to generate seasonal dummies or fourier to generate the trigonometric representation of seasonal dummies. For illustration and as a complement to the results already shown by @IrishStat, below I show what I get using the first 10 components of the seasonal cycles.

seas <- fourier(x, K=10)
trend <- seq_along(x)
mo <- outliers(c("TC", rep("AO", 5)), c(2, 14, 42, 65, 120, 148))
outl <- outliers.effects(mo, length(x))
seas <- fourier(x, K=10)
trend <- seq_along(x)
res2 <- arima(x, order = c(1,0,0), xreg = cbind(seas, trend, outl))

This yields the following residuals

plot(residuals(res2))

par(mfrow = c(2,1), mar = c(5,4,4,3))
acf(residuals(res2), lag.max = 52)
pacf(residuals(res2), lag.max = 52)

Now, we can see that the autocorrelations of order 52 are less prominent in the residuals.

The following estimate of the seasonal pattern is obtained:

fitted.seas <- ts(seas %*% coef(res2)[3:22])
tsp(fitted.seas) <- tsp(x)
plot(fitted.seas, type = "l", main = "estimated seasonal component")

Finally, as you used R for the results shown in you question, if you want to work with the model proposed by @IrishStat, you can do the following:

mo <- outliers(rep("AO", 5), c(14, 65, 42, 39, 120))
outl <- outliers.effects(mo, length(x))
seas <- matrix(0, nrow = length(x), ncol = 52)
for (i in seq(52))
  seas[seq(i, length(x), 52),i] <- 1
trend <- seq_along(x)
res3 <- arima(x, order = c(1,0,0), xreg = cbind(seas[,-c(1,18,33)], trend, outl))

Answer 2

In the spirit of Christmas and a sincere respect for the persistence of the OP for an answer to his reasonable question I have taken the 152 weekly observations ( nearly 3 years of data) and used AUTOBOX in a totally automatic hands-free mode, a piece of software that I have helped develop , to give what I think is a more complete and comprehensive answer. One could use any flexible analysis system to duplicate my results with some requisite programming required to accomplish some technical advances. @javlacalle (previously reported results) concluded that although the residual ACF ( at period ) suggested a cyclic effect but he failed to provide a complete solution although he had some clues to what was wrong with his well-intentioned approach. He I believe didn’t have the facility to seamlessly integrate a set of weekly dummies into his model. The problem here is very straight-forward : the identification of a significant lag at period 52 in the acf can be evidence of either the need for an ARIMA augmentation OR the need for a deterministic seasonal dummy augmentation . Before Box and Jenkins and their army of zealots (of which I was a leader!) limited there approach to ARIMA (memory based) solutions a rich history of deterministic models (i.e. seasonal dummies ) was in vogue and apparently has been ignored . Developments in integrating deterministic seasonal structure in the presence of. ARIMA structure and empirically identified interventions in the presence of time varying coefficients and time varying error variance are now routinely available. In this case the toolkit used by @javlacalle whose scholarship I generally respect was not robust enough in my opinion to correctly and fully solve the problem suggested by the data. The simplistic approach to use initial ARIMA model identification assuming among other things that there were no anomalies, no deterministic seasonal dummies, no deterministic trends prior to Intervention Detection , which assumes a “good ARIMA model”, seems to have failed as the series has a significant upwards trend in the data and contains strong deterministic weekly effects. His Intervention Detection schemes yielded a mixed bag of results as the tentative set of residuals data was strongly affected by the unspecified i.e. explicitly omitted/ needed weekly dummies and a time trend.

I have posted my complete results on the web at http://www.autobox.com/stack/weekly.zip . The final model included an AR(1) coefficient value .606 curiously close to the reported .617 and 48 weekly dummies (missing weeks 1,18 and 33) , a trend value of 5.54 and some 5 Additive Outliers (AO) at periods 12,14,65,42, and 39. The final model statistics showing an MSE of 44800 a reduction of 50% from the previously reported 88411 emphasizes the differences in the models. Following is a plot of the residuals suggesting randomness . The final ACF of the residuals , noticeably free of any lag 52 effect is available in weekly.zip as finalacf.txt . A plot of the final model residuals is shown here. It is important to notice that the range of the errors is dramatically reduced as compared to the the previously reported results. The fit/forecast graph nicely presents the essence of the model and is shown here. The Cleansed Data graph nicely highlights the parsimonious identification of 5 anomalies. The incorporation of weekly dummies requires additional parameters (needed in my opinion as they are statistically significant) although I can already hear the gasps of some readers who look askance at the utilization of so many parameters.

The Actual/Fit and Forecast graph is also useful in the understanding of these powerful results. While all models are wrong, this model seems particularly useful in explaining the past and predicting the next 52 weeks.

In summary while the ARIMA structures are the same, the models differ significantly in terms of trend and seasonal components with a reduction of nearly 50% in the estimated error variance. Restated the reported results have a 100% larger error variance than my suggested model. I feel sure that if it were somehow possible to automatically integrate the weekly dummies the results would still suggest a deficit. The trend and seasonality visually evident to the OP are now safely reflected in the model and the forecast. Finally a cyclical model is not necessarily an auto-regressive (memory ) model it could be a seasonal dummy model as in this case. QED