# Holt-Winters with a Quantile Loss Function

Exponential Smoothing is an old technique, but it can perform extremely well on real time series, as discussed in Hyndman, Koehler, Ord & Snyder (2008)),

when Gardner (2005) appeared, many believed that exponential smoothing should be disregarded because it was either a special case of ARIMA modeling or an ad hoc procedure with no statistical rationale. As McKenzie (1985) observed, this opinion was expressed in numerous references to my paper. Since 1985, the special case argument has been turned on its head, and today we know that exponential smoothing methods are optimal for a very general class of state-space models that is in fact broader than the ARIMA class.

Furthermore, I like it because I think it has nice pedagogical features. Consider simple exponential smoothing, $$L_{t}=\alpha Y_{t}+(1-\alpha)L_{t-1}$$ where $\alpha\in(0,1)$ is the smoothing weight. It is locally constant, in the sense that ${}_{t}\hat Y_{t+h} = L_{t}$

 library(datasets) X=as.numeric(Nile) SimpleSmooth = function(a){ T=length(X) L=rep(NA,T) L[1]=X[1] for(t in 2:T){L[t]=a*X[t]+(1-a)*L[t-1]} return(L) } plot(X,type="b",cex=.6) lines(SimpleSmooth(.2),col="red")

When using the standard R function, we get

hw=HoltWinters(X,beta=FALSE,gamma=FALSE, l.start=X[1]) hw$alpha [1] 0.2465579 Of course, one can replicate that optimal value V=function(a){ T=length(X) L=erreur=rep(NA,T) erreur[1]=0 L[1]=X[1] for(t in 2:T){ L[t]=a*X[t]+(1-a)*L[t-1] erreur[t]=X[t]-L[t-1] } return(sum(erreur^2)) } optim(.5,V)$par [1] 0.2464844

Here, the optimal value for $\alpha$ is the one that minimizes the one-step prediction, for the $\ell_2$ loss function, i.e. $$\sum_{t=2}^n(Y_t-{}_{t-1}\hat Y_t)^2$$ where here ${}_{t-1}\hat Y_t = L_{t-1}$. But one can consider another loss function, for instance the quantile loss function, $$\ell_{\tau}(\varepsilon)=\varepsilon(\tau-\mathbb{I}_{\varepsilon\leq 0})$$. The optimal coefficient is then obtained using

HWtau=function(tau){ loss=function(e) e*(tau-(e&lt;=0)*1) V=function(a){ T=length(X) L=erreur=rep(NA,T) erreur[1]=0 L[1]=X[1] for(t in 2:T){ L[t]=a*X[t]+(1-a)*L[t-1] erreur[t]=X[t]-L[t-1] } return(sum(loss(erreur))) } optim(.5,V)$par } Here is the evolution of $\alpha^\star_\tau$ as a function of $\tau$ (the level of the quantile considered). T=(1:49)/50 HW=Vectorize(HWtau)(T) plot(T,HW,type="l") abline(h= hw$alpha,lty=2,col="red")

Note that the optimal $\alpha$ is decreasing with $\tau$. I wonder how general this result can be…

Of course, one can consider more general exponential smoothing, for instance the double one, with $$L_t=\alpha Y_t+(1-\alpha)[L_{t-1}+B_{t-1}]$$and$$B_t=\beta[L_t-L_{t-1}]+(1-\beta)B_{t-1}$$so that the prediction is now ${}_{t}\hat Y_{t+h} = L_{t}+hB_t$ (it is now locally linear – and no longer constant).

hw=HoltWinters(X,gamma=FALSE,l.start=X[1]) hw$alpha alpha 0.4200241 hw$beta beta 0.05973389

The code to compute the smoothed series is the following

DoubleSmooth = function(a,b){ T=length(X) L=B=rep(NA,T) L[1]=X[1]; B[1]=0 for(t in 2:T){ L[t]=a*X[t]+(1-a)*(L[t-1]+B[t-1]) B[t]=b*(L[t]-L[t-1])+(1-b)*B[t-1] } return(L+B) }

Here also it is possible to replicate R using the $\ell_2$ loss function

V=function(A){ a=A[1] b=A[2] T=length(X) L=B=erreur=rep(NA,T) erreur[1]=0 L[1]=X[1]; B[1]=X[2]-X[1] for(t in 2:T){ L[t]=a*X[t]+(1-a)*(L[t-1]+B[t-1]) B[t]=b*(L[t]-L[t-1])+(1-b)*B[t-1] erreur[t]=X[t]-(L[t-1]+B[t-1]) } return(sum(erreur^2)) } optim(c(.5,.05),V)$par [1] 0.41904510 0.05988304 (up to numerical optimization approximation, I guess). But here also, a quantile loss function can be considered HWtau=function(tau){ loss=function(e) e*(tau-(e&lt;=0)*1) V=function(A){ a=A[1] b=A[2] T=length(X) L=B=erreur=rep(NA,T) erreur[1]=0 L[1]=X[1]; B[1]=X[2]-X[1] for(t in 2:T){ L[t]=a*X[t]+(1-a)*(L[t-1]+B[t-1]) B[t]=b*(L[t]-L[t-1])+(1-b)*B[t-1] erreur[t]=X[t]-(L[t-1]+B[t-1]) } return(sum(loss(erreur))) } optim(c(.5,.05),V)$par }

and we can plot those values on a graph

T=(1:49)/50 HW=Vectorize(HWtau)(T) plot(HW[1,],HW[2,],type="l") abline(v= hw$alpha,lwd=.4,lty=2,col="red") abline(h= hw$beta,lwd=.4,lty=2,col="red") points(hw$alpha,hw$beta,pch=19,col="red")

(with $\alpha$ on the $x$-axis, and $\beta$ on the $y$-axis). So here, it is extremely simple to change the loss function, but so far, it should be done manually. Of course, one do it also for the seasonal exponential smoothing model.

# Forecast, Automatic Routines vs. Experience

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

> urls = "http://freakonometrics.free.fr/report-headphones-2015.csv"
> tail(report)
580 2015-02-08 - 2015-02-14         53
581 2015-02-15 - 2015-02-21         52
582 2015-02-22 - 2015-02-28         51
583 2015-03-01 - 2015-03-07         50
584 2015-03-08 - 2015-03-14         49
585 2015-03-15 - 2015-03-21         49

If we plot that weekly time series, we have

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

> source(
+   "http://freakonometrics.blog.free.fr/public/code/H2M.R")
> plot(headphones)