{"word_count": 809, "character_count": 7138}
This morning, in our Time Series course, we've been playing with some data I got from google.ca/trends/. Actually, we've been playing on some old version, downloaded 18 months ago (discussed in a previous post, in French).&#62; urls = "http://freakonometrics.free.fr/report-headphones-2015.csv"&#62; report=read.table(+ urls,skip=4,header=TRUE,sep=",",nrows=585)&#62; tail(report)                    Semaine headphones580 2015-02-08 - 2015-02-14         53581 2015-02-15 - 2015-02-21         52582 2015-02-22 - 2015-02-28         51583 2015-03-01 - 2015-03-07         50584 2015-03-08 - 2015-03-14         49585 2015-03-15 - 2015-03-21         49If we plot that weekly time series, we have&#62; plot(report[,2],type="l")Working with weekly series is more complicated (at least to find a simple model, with only a few lags), so let us convert that series into a monthly one,&#62; source(+   "http://freakonometrics.blog.free.fr/public/code/H2M.R")&#62; headphones=H2M(report,lang="FR",type="ts")&#62; plot(headphones)Here, we did three things,we considered the logarithm of the series, because of the increasing variability we have on recent years,we removed (too) old observationswe kept recent observations to test our models.&#62; X=log(headphones)&#62; T=as.numeric(time(X))&#62; idx=which((T&#62;2006)&#38;(T&#60;2014.3))&#62; T=T[idx]&#62; X=as.numeric(X)[idx]&#62; plot(T,X,type="l")&#62; df=data.frame(X,T)&#62; reg=lm(X~T,data=df)&#62; abline(reg,col="red")&#62; Y=residuals(reg)The residuals obtained after fitting that linear model will be our (supposed to be) stationnary time series. Let us consider various models.model 1We obtained our first model be playing with the orders, this morning. Here, the model is slightly more complex,&#62; fit1=arima(Y,+   order = c(1, 0, 10),+   seasonal = list(order = c(1, 1, 0), +   period = 12))&#62; fit1 Call:arima(x = Y, order = c(1, 0, 10), seasonal = list(order = c(1, 1, 0), period = 12)) Coefficients:        ar1     ma1    ma2     ma3    ma4           0.930  -0.033  0.078  -0.101  0.365  s.e.  0.064   0.102  0.132   0.119  0.152   ma5    ma6     ma7    ma8    ma9   ma10-0.217  0.216  -0.137  0.134  0.048  0.465  0.154  0.187   0.166  0.173  0.153  0.154   sar1-0.578 0.118 sigma^2 estimated as 0.00114:  log likelihood = 166,  aic = -305If we look at the residuals, we have some white noiseWe have here our first model.Model 2For the second one, I was much more lazy&#62; library(forecast)&#62; auto.arima(ts(Y,start=2006,+ frequency=12),max.p=12,max.q=12)Series: ts(Y, start = 2006, frequency = 12) ARIMA(0,1,0)(0,1,1)[12]                     Coefficients:        sma1      -0.950s.e.   0.672 sigma^2 estimated as 0.00124:  log likelihood=155AIC=-306   AICc=-306   BIC=-301&#62; fit2=arima(Y,+            order = c(0, 1, 0),+            seasonal = list(order = c(0, 1, 1), +            period = 12))&#62; plot(acf(residuals(fit2),lag=36),lwd=3)Just to confirm that we actually have a white noise, here, use&#62; BP = function(h) Box.test(residuals(fit2),lag=h,type='Box-Pierce')$p.value&#62; plot(1:24,Vectorize(BP)(1:24),type='b',+ ylim=c(0,1),col="red")&#62; abline(h=.05,lty=2,col="blue")One more time, we have a white noise. So this model here is valid.Model 3&#62; fit3=arima(Y,+            order = c(12, 0, 0),+            seasonal = list(order = c(0, 1, 1), +            period = 12))&#62; fit3 Call:arima(x = Y, order = c(12, 0, 0), seasonal = list(order = c(0, 1, 1), period = 12)) Coefficients:        ar1     ar2    ar3    ar4     ar5    ar6      0.918  -0.044  0.014  0.253  -0.366  0.347s.e.  0.112   0.148  0.145  0.135   0.142  0.151         ar7     ar8    ar9    ar10   ar11      -0.205  -0.030  0.395  -0.208  0.116s.e.   0.152   0.139  0.134   0.162  0.169        ar12    sma1      -0.249  -0.767s.e.   0.110   0.164Again, we go a model similar, this morning, simply by playing with the orders. And one more time, we get a white noise,Model 4For the last model, we can consider some exponential smoothing model&#62; HW = HoltWinters(log(headphones))(again, on the logarithm of the series).ForecastingFor the first three models, use&#62; Y1=predict(fit1,n.ahead=36)&#62; Y2=predict(fit2,n.ahead=36)&#62; Y3=predict(fit3,n.ahead=36)To visualise the three forecast, use, e.g.&#62; draw=function(fit=Y1){+ straight_line=predict(reg,+        newdata=data.frame(T=Tfutur))+ mth=straight_line+fit$pred+ sth2=fit$se^2+ Xmy=exp(mth+sth2/2)+ qinf=qlnorm(.05,mth,sdlog=sqrt(sth2))+ qsup=qlnorm(.95,mth,sqrt(sth2))+ + plot(headphones,type="l",xlim=c(2006,2016),+      ylim=c(25,105))+ polygon(c(Tfutur,rev(Tfutur)),+         c(qinf,rev(qsup)),+         col="yellow",border=NA)+ lines(Tfutur,Xmy,col="red")+ return(Xmy)+ }Here, we have&#62; draw(Y1)Time Series:Start = 100 End = 135 Frequency = 1  [1]  45.6  46.3  47.6  52.1  52.7  52.9  68.3   [8]  88.8  58.5  51.2  49.5  47.4  47.8  49.4  [15]  51.6  54.6  54.4  55.2  73.1  95.9  64.2  [22]  56.8  54.5  52.3  53.4  54.8  56.9  61.5  [29]  62.0  62.6  81.8 106.9  71.0  62.6  60.4  [35]  57.9&#62; draw(Y2)&#62; draw(Y3)And finally, for the fourth model we can use&#62; Y4=predict(HW,n.ahead=36,+    prediction.interval = TRUE)&#62; sth2=((Y4[,"upr"]-Y4[,"fit"])/1.96)^2&#62; mth=Y4[,"fit"]&#62; Xmy=exp(mth+sth2/2)&#62; XHW=Xmy&#62; qinf=qlnorm(.05,mth,sdlog=sqrt(sth2))&#62; qsup=qlnorm(.95,mth,sqrt(sth2))&#62; plot(headphones,type="l",xlim=c(2006,2016),+      ylim=c(25,105))&#62; polygon(c(Tfutur,rev(Tfutur)),+         c(qinf,rev(qsup)),+         col="yellow",border=NA)&#62; lines(Tfutur,Xmy,col="red")Model comparisonHere, we can compare the forecast with actual values&#62; X=headphones&#62; T=as.numeric(time(X))&#62; idx=which((T&#62;=2014.3))&#62; T=T[idx]&#62; X=as.numeric(X)[idx]&#62; D=data.frame(T=T,X=X,X1=draw(Y1)[1:length(T)],X2=draw(Y2)[1:length(T)],X3=draw(Y3)[1:length(T)],X4=XHW[2:(1+length(T))])&#62; D      T    X   X1   X2   X3   X41  2014 44.8 45.6 43.7 44.3 49.02  2014 45.2 46.3 44.3 44.9 49.63  2014 46.8 47.6 45.5 45.6 51.24  2015 49.3 52.1 47.7 48.8 53.85  2015 49.8 52.7 47.9 49.7 53.66  2015 50.9 52.9 49.0 50.8 54.47  2015 69.8 68.3 62.3 64.9 72.38  2015 91.2 88.8 82.5 85.4 92.49  2015 57.5 58.5 58.0 59.4 59.810 2015 52.5 51.2 51.9 53.6 54.911 2015 49.3 49.5 49.3 51.4 52.6We can plot the errors,&#62; barplot(D$X1-D$X,ylim=c(-7,5))&#62; barplot(D$X2-D$X,ylim=c(-7,5))x&#62; barplot(D$X3-D$X,ylim=c(-7,5))&#62; barplot(D$X4-D$X,ylim=c(-7,5))To compare the four model, compute the sum of the squares of those errors,&#62; L2_dist=function(var){+   return(sum((D[,"X"]-D[,var])^2))+ } &#62; sapply(paste("X",1:4,sep=""),L2_dist)   X1    X2    X3    X4  33.2 145.3  68.0 133.9Model 2 and Model 4 were obtained using automatic routines. One on ARIMA models, the other one on the class of exponential smoothing models. The interesting point is thatthe we've been able to obtain two better models, just by playing a litle bit. From an initial guess, we look at the residuals, try a another one, and again another one, until the residuals series is a white noise... Somehow I like this idea, that I can still beat automatic routines, just using experience I have, now, on time series.