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

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.

• There's upward trend, and some evidence of seasonality. I see two options to start with: ARX(1), where time is one of the variables, and ARI(1,1) both with seasonality. I wouldn't worry about MA in the beginning, you can try it later. – Aksakal Dec 11 '14 at 19:48
• You are assuming a certain kind of model when you apply STL. There may be other models that are more appropriate such as a model with multiple trends or multiple level shifts. Furthermore the seasonal structure may be determinstic or autoprojective. Furthermore there may be anomalies in the data that left untreated distort the analysis. Furthermore there may be transience in either the model parameters or the error variance. To answer what your model might be , please post your data in an excel form. – IrishStat Dec 11 '14 at 20:02
• drive.google.com/file/d/0B7uexFVHb9uwLXhpX3hIUGctVTQ/… – kebrab67 Dec 11 '14 at 21:58
• link – kebrab67 Dec 11 '14 at 21:59
• Here is alink to the data I give access to anyone who needs to see it to download. Thanks for the help. – kebrab67 Dec 11 '14 at 22:00

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))


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