{"word_count": 465, "character_count": 4791}
Exponential Smoothing is an old technique, but it can perform extremely well on real time series, as discussed in Hyndman, Koehler, Ord &#38; 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, [latex display="true"]L_{t}=\alpha Y_{t}+(1-\alpha)L_{t-1}[/latex] where [latex]\alpha\in(0,1)[/latex] is the smoothing weight. It is locally constant, in the sense that [latex]{}_{t}\hat Y_{t+h} = L_{t}[/latex] 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 gethw=HoltWinters(X,beta=FALSE,gamma=FALSE, l.start=X[1])hw$alpha[1] 0.2465579Of course, one can replicate that optimal valueV=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.2464844Here, the optimal value for [latex]\alpha[/latex] is the one that minimizes the one-step prediction, for the [latex]\ell_2[/latex] loss function, i.e. [latex display="true"]\sum_{t=2}^n(Y_t-{}_{t-1}\hat Y_t)^2[/latex] where here [latex]{}_{t-1}\hat Y_t = L_{t-1}[/latex]. But one can consider another loss function, for instance the quantile loss function, [latex display="true"]\ell_{\tau}(\varepsilon)=\varepsilon(\tau-\mathbb{I}_{\varepsilon\leq 0})[/latex]. The optimal coefficient is then obtained usingHWtau=function(tau){loss=function(e) e*(tau-(e&#60;=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 [latex]\alpha^\star_\tau[/latex] as a function of [latex]\tau[/latex] (the level of the quantile considered).T=(1:49)/50HW=Vectorize(HWtau)(T)plot(T,HW,type="l")abline(h= hw$alpha,lty=2,col="red")Note that the optimal [latex]\alpha[/latex] is decreasing with [latex]\tau[/latex]. I wonder how general this result can be...Of course, one can consider more general exponential smoothing, for instance the double one, with [latex display="true"]L_t=\alpha Y_t+(1-\alpha)[L_{t-1}+B_{t-1}][/latex]and[latex display="true"]B_t=\beta[L_t-L_{t-1}]+(1-\beta)B_{t-1}[/latex]so that the prediction is now [latex]{}_{t}\hat Y_{t+h} = L_{t}+hB_t[/latex] (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.05973389The code to compute the smoothed series is the followingDoubleSmooth = 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 [latex]\ell_2[/latex] loss functionV=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 consideredHWtau=function(tau){loss=function(e) e*(tau-(e&#60;=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 graphT=(1:49)/50HW=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 [latex]\alpha[/latex] on the [latex]x[/latex]-axis, and [latex]\beta[/latex] on the [latex]y[/latex]-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.