0. Introducció

La suavització exponencial es va proposar a finals de la dècada de 1950 (Brown, 1959; Holt, 1957; Winters, 1960) i ha motivat alguns dels mètodes de predicció més reeixits. Les prediccions produïdes mitjançant mètodes de suavització exponencial són mitjanes ponderades d’observacions passades, amb els pesos que disminueixen exponencialment a mesura que les observacions es fan més antigues. En altres paraules, com més recent sigui l’observació, més alt serà el pes associat. Aquest marc genera prediccions fiables i per a una àmplia gamma de sèries temporals, cosa que és un gran avantatge i de gran importància per a les aplicacions a la indústria.

Aquí presentarem la mecànica dels mètodes de suavització exponencial més importants i la seva aplicació en la predicció de sèries temporals amb diverses característiques. Això ens ajuda a desenvolupar una intuïció de com funcionen aquests mètodes. En aquest context, seleccionar i utilitzar un mètode de predicció pot semblar una mica improvisat. La selecció del mètode es basa generalment en el reconeixement dels components clau de la sèrie temporal (tendència i estacional) i la manera com aquests entren al mètode de suavització (per exemple, de manera additiva, amortiguada o multiplicativa).

Posreiorment presentem els models estadístics que fonamenten els mètodes de suavització exponencial. Aquests models generen previsions acurades idèntiques als mètodes discutits a la primera part del capítol, però també generen intervals de predicció. A més, aquest marc estadístic permet una selecció real de models entre models competidors.

Referències

G. J. Hyndman and G. Athanasopoulos, Forecasting: Principles and Practice,3rd edition, Chapter 5, 2021, OTexts: Melbourne, Australia. Disponible a otexts.com/.
També la segona edició del mateix text. Moltes de les solucions de R que apliquem estan més en relació amb aquesta edició.

Datasets

La llibreria fp3 conté tots els conjunts de dades necessaris per als exemples i exercicis del llibre “Forecasting: principles and practice, 3d ed.” de Rob J Hyndman i George Athanasopoulos. Per aixo s’utilitzarà principalment datasets d’aquest paquest.

La llibreria tsibbledata proporciona diversos conjunts de dades a l’estructura de dades “tsibble”. Aquests conjunts de dades són útils per aprendre i demostrar com es poden ordenar, visualitzar i preveure les dades temporals ordenades.

La lliberia fma conté tots els conjunts de dades de “Forecasting: methods and applications” de Makridakis, Wheelwright i Hyndman (Wiley, 3a ed., 1998).

library(fpp2)
library(fpp3)
library(tsibbledata)
library(fma)

Possiblement se n’utilitzaran d’altres, però aquestes seran la més usades.

Llibreries principals

Llibreries generals

A més de les llibreries base de R, atilitzarem algunes llibreria de l’entorn tidyverse::dplyr per a la manipulació de dades, ggplot2 per a gràfics i patchwork per a gràfics múltiples. Puntualment, n’hem utilitzat d’altres, que apareixeran indicades al codi.

library(ggplot2)
library(dplyr)
library(patchwork)

Lliberiries d’estadística

Pel que fa a les funcions estadístiques, a més de stats que ve per defecte amb R, farem servir en vàries ocasions forecast (veure) i astsa que són llibreries específiques per anàlisi de sèries temposrals. També es podran fer servir altres paquets quan no estigui disponible alguna funció o les disponibles no tinguin totes les característiques desitjades. En tot cas, en el codi sempre es trobarà indicada la llibreria que correspongui.

Advertència

Ateses le nombroses ambigüitats que sovint es poden generar, sovint les funcions i els datasets van acompanyats pel nom del paque d’origen (per exemple, ggplot2::autoplot(), stats::acf(), etc.).

1. Suavització exponencial simple

El més simple dels mètodes de suavització exponencial s’anomena naturalment suavització exponencial simple (SES). Aquest mètode és adequat per a la previsió de dades sense una tendència o un patró estacional clar. Per exemple, les dades de la Figura 8.1 no mostren cap comportament de tendència clar ni cap estacionalitat. (Hi ha un descens en els darrers anys, cosa que podria suggerir una tendència. Més endavant en aquest capítol, considerarem si un mètode de tendència seria millor per a aquesta sèrie.) Ja hem considerat el mètode ingenu i el mètode mitjà com a possibles mètodes per predir aquestes dades.

algeria_economy <- 
  tsibbledata::global_economy |>
  dplyr::filter(Country == "Algeria")

algeria_economy |>
  ggplot2::autoplot(Exports) +
  labs(y = "% del PIB", x="Any", title = "Exportacions: Algèria")+
  theme_bw()

Utilitzant el mètode naïf, totes les previsions per al futur són iguals a l’últim valor observat de la sèrie, \[ \hat{y}_{T+h|T} = y_{T}, \]

per \(h=1,2,…\). Per tant, el mètode ingenu assumeix que l’observació més recent és l’única important, i totes les observacions anteriors no proporcionen informació per al futur. Això es pot considerar com una mitjana ponderada on tot el pes es dóna a l’última observació.

Utilitzant el mètode de la mitjana, totes les previsions futures són iguals a una mitjana simple de les dades observades,

\[ \hat{y}_{T+h|T} = \frac1T \sum_{t=1}^T y_t, \] per a \(h=1,2,…\). Per tant, el mètode de la mitjana assumeix que totes les observacions tenen la mateixa importància i els dóna pesos iguals a l’hora de generar previsions.

Sovint volem alguna cosa entre aquests dos extrems. Per exemple, pot ser sensat assignar pesos més grans a les observacions més recents que a les observacions del passat llunyà. Aquest és exactament el concepte darrere del suavització exponencial simple. Les previsions es calculen mitjançant mitjanes ponderades, on els pesos disminueixen exponencialment a mesura que les observacions provenen de més enllà en el passat; els pesos més petits s’associen amb les observacions més antigues: \[ \begin{equation} \hat{y}_{T+1|T} = \alpha y_T + \alpha(1-\alpha) y_{T-1} + \alpha(1-\alpha)^2 y_{T-2}+ \cdots, \tag{8.1} \end{equation} \]

on \(0 \le \alpha \le 1\) és el paràmetre de suavització. La previsió d’un pas endavant per al temps \(T+1\) és una mitjana ponderada de totes les observacions de la sèrie \(y_1,\dots,y_T\). La velocitat a la qual disminueixen els pesos està controlada pel paràmetre \(α\).

La taula següent mostra els pesos assignats a les observacions per a quatre valors diferents de \(α\) quan es fa una previsió mitjançant un suavització exponencial simple. Cal tenir en compte que la suma dels pesos, fins i tot per a un valor petit de \(α\) serà aproximadament una per a qualsevol mida de mostra raonable

\[ \begin{array}{|l|c|c|c|c|} & α=0.2& α=0.4& α=0.6& α=0.8 \\ \hline y_T& 0.2& 0.4& 0.6& 0.8 \\ y_{T−1}& 0.16 & 0.24 &0.24 &0.16\\ y_{T−2}& 0.128 &0.144 &0.096 &0.032\\ y_{T−3}& 0.1024& 0.0864& 0.0384& 0.0064\\ y_{T−4}& 0.0819& 0.0518& 0.0154& 0.0013\\ y_{T−5}& 0.0655& 0.0311& 0.0061& 0.0003\\ \hline \end{array} \]

Per a qualsevol \(α\) entre 0 i 1, els pesos assignats a les observacions disminueixen exponencialment a mesura que retrocedim en el temps, d’aquí el nom “suavitzat exponencial”. Si \(α\) és petit (és a dir, proper a 0), es dóna més pes a les observacions del passat més llunyà. Si \(α\) és gran (és a dir, proper a 1), es dóna més pes a les observacions més recents. Per al cas extrem on \(\alpha=1,\;\hat{y}_{T+1|T}=y_T\), de manera que les previsions són iguals a les previsions ingenues.

Presentem dues formes equivalents de suavitzacio exponencial simple, cadascuna de les quals condueix a l’equació de previsió (8.1).

Suavització exponencial per mitjana ponderada

La previsió en el moment \(T+1\) és igual a una mitjana ponderada entre l’observació més recent \(y_T\) i la previsió anterior \(\hat{y}_{T|T-1}\):

\[ \hat{y}_{T+1|T} = \alpha y_T + (1-\alpha) \hat{y}_{T|T-1}, \] on \(0≤α≤1\) és el paràmetre de suavització. De manera similar, podem escriure els valors ajustats com \[ \hat{y}_{t+1|t} = \alpha y_t + (1-\alpha) \hat{y}_{t|t-1}, \] per a \(t=1,…,T\) (Recordem que els valors ajustats són simplement previsions d’un pas de les dades d’entrenament).

El procés ha de començar en algun lloc, de manera que deixem que el primer valor ajustat en el temps 1 es denoti per \(ℓ_0\) (que haurem d’estimar). Aleshores

\[ \begin{align*} \hat{y}_{2|1} &= \alpha y_1 + (1-\alpha) \ell_0\\ \hat{y}_{3|2} &= \alpha y_2 + (1-\alpha) \hat{y}_{2|1}\\ \hat{y}_{4|3} &= \alpha y_3 + (1-\alpha) \hat{y}_{3|2}\\ \vdots\\ \hat{y}_{T|T-1} &= \alpha y_{T-1} + (1-\alpha) \hat{y}_{T-1|T-2}\\ \hat{y}_{T+1|T} &= \alpha y_T + (1-\alpha) \hat{y}_{T|T-1}. \end{align*} \]

Substituint cada equació a l’equació següent, obtenim

\[ \begin{align*} \hat{y}_{3|2} & = \alpha y_2 + (1-\alpha) \left[\alpha y_1 + (1-\alpha) \ell_0\right] \\ & = \alpha y_2 + \alpha(1-\alpha) y_1 + (1-\alpha)^2 \ell_0 \\ \hat{y}_{4|3} & = \alpha y_3 + (1-\alpha) [\alpha y_2 + \alpha(1-\alpha) y_1 + (1-\alpha)^2 \ell_0]\\ & = \alpha y_3 + \alpha(1-\alpha) y_2 + \alpha(1-\alpha)^2 y_1 + (1-\alpha)^3 \ell_0 \\ & ~~\vdots \\ \hat{y}_{T+1|T} & = \sum_{j=0}^{T-1} \alpha(1-\alpha)^j y_{T-j} + (1-\alpha)^T \ell_{0}. \end{align*} \] L’últim terme esdevé petit per a \(T\) gran. Per tant, la forma de mitjana ponderada condueix a la mateixa equació de previsió (8.1).

Suavització amb components

Una representació alternativa és la forma de components. Per a la suavització exponencial simple, l’únic component inclòs és el nivell, \(ℓ_t\). (Altres mètodes que es consideren més endavant en aquest capítol també poden incloure una tendència \(b_t\) i un component estacional \(s_t\).) Les representacions en forma de components dels mètodes de suavització exponencial comprenen una equació de previsió i una equació de suavització per a cadascun dels components inclosos en el mètode. La forma de components de la suavització exponencial simple ve donada per:

\[ \begin{align*} \text{Forecast equation} && \hat{y}_{t+h|t} & = \ell_{t}\\ \text{Smoothing equation} && \ell_{t} & = \alpha y_{t} + (1 - \alpha)\ell_{t-1}, \end{align*} \]

recordant que \(ℓ_t\) és una mitjana ponderada (amb ponderació \(\alpha\) de l’observació \(y_t\) i la previsió d’entrenament amb un pas d’antelació per al temps \(t\).

Si substituïm \(\ell_tt\) per \(\hat{y}_{t+1|t}\) i \(\ell_{t-1}\) per \(\hat{y}_{t|t-1}\) a l’equació de suavització, recuperarem la forma mitjana ponderada de la suavització exponencial simple.

La forma de components de la suavització exponencial simple no és particularment útil per si sola, però serà la forma més fàcil d’utilitzar quan comencem a afegir altres components.

Previsions planes

La suavització exponencial simple té una funció de previsió “plana”: \[ \hat{y}_{T+h|T} = \hat{y}_{T+1|T}=\ell_T, \qquad h=2,3,\dots. \]

És a dir, totes les previsions prenen el mateix valor, igual a l’últim component de nivell. Recordeu que aquestes previsions només seran adequades si la sèrie temporal no té cap tendència ni component estacional.

Optimització

L’aplicació de cada mètode de suavització exponencial requereix que es triïn els paràmetres de suavització i els valors inicials. En particular, per a la suavització exponencial simple, hem de seleccionar els valors d’\(α\) i \(ℓ_0\). Totes les previsions es poden calcular a partir de les dades un cop coneixem aquests valors. Per als mètodes següents, normalment hi ha més d’un paràmetre de suavització i més d’un component inicial a triar.

En alguns casos, els paràmetres de suavització es poden triar de manera subjectiva: el pronosticador especifica el valor dels paràmetres de suavització basant-se en l’experiència prèvia. Tanmateix, una manera més fiable i objectiva d’obtenir valors per als paràmetres desconeguts és estimar-los a partir de les dades observades.

En regressió linial, els coeficients s’estimen del model de regressió minimitzant la suma dels residuals al quadrat (normalment coneguda com a SSE o “suma d’errors al quadrat”). De la mateixa manera, els paràmetres desconeguts i els valors inicials per a qualsevol mètode de suavització exponencial es poden estimar minimitzant la SSE. Els residuals s’especifiquen com a \(e_t=y_t - \hat{y}_{t|t-1}\) per a \(t=1,\dots,T\). Per tant, trobem els valors dels paràmetres desconeguts i els valors inicials que minimitzen \[ \begin{equation} \text{SSE}=\sum_{t=1}^T(y_t - \hat{y}_{t|t-1})^2=\sum_{t=1}^Te_t^2. \tag{8.2} \end{equation} \]

A diferència del cas de regressió (on tenim fórmules que retornen els valors dels coeficients de regressió que minimitzen l’SSE), això implica un problema de minimització no lineal i necessitem utilitzar una eina d’optimització per resoldre’l.

Exemple: le exportacions d’Algèria

Ena aquest exemple s’aplica la suavització exponencial a la predicció de sobre l’exportació de bens di serveis d’Algèria.

La funció `ets()`

Utilitzarem la funció ets() de forecast::. Aquesta funció té un parametre important per a definir el model de suavització que es vol dur a teme: model. És una cadena de tres carácters (tipus “ZZZ”, per defecte) que utilitza la terminologia del marc de treball de Hyndman et al. (2002) i Hyndman et al. (2008). La primera lletra indica el tipus d’error (“A”, “M” o “Z”); la segona lletra indica el tipus de tendència (“N”, “A”, “M” o “Z”); i la tercera lletra indica el tipus de estacionalitat (“N”, “A”, “M” o “Z”). En tots els casos, “N” = cap, “A” = additiu, “M” = multiplicatiu i “Z” = seleccionat automàticament. Així, per exemple, “ANN” és una suavització exponencial simple amb errors additius, “MAM” és el mètode multiplicatiu de Holt-Winters amb errors multiplicatius, etc.

També és possible que el model sigui de la classe ets() i igual a la sortida d’una crida anterior a ets(). En aquest cas, el mateix model s’ajusta a y sense reestimar cap paràmetre de suavització. Vegeu també l’argument use.initial.values.

tsalgeria <- ts(algeria_economy[,8], start = 1960, frequency = 1)

ajustalg <- ets(tsalgeria, model = "ANN")

ajustalg

## ETS(A,N,N) 
## 
## Call:
## ets(y = tsalgeria, model = "ANN")
## 
##   Smoothing parameters:
##     alpha = 0.84 
## 
##   Initial states:
##     l = 39.539 
## 
##   sigma:  5.9691
## 
##      AIC     AICc      BIC 
## 446.7154 447.1599 452.8968

Això dóna estimacions de paràmetres \(\hat α=0.84\) i \(\hat ℓ_0=39.5\), obtingudes minimitzant SSE durant períodes \(t=1, 2,…, 58\), subjectes a la restricció que \(0≤α≤1\).

A la Taula 8.1 demostrem el càlcul utilitzant aquests paràmetres. La penúltima columna mostra el nivell estimat per als temps de \(t=0\) a \(t=58\); les últimes files de l’última columna mostren les previsions per a \(h=1\) a 5 passos més endavant.

Taula 8.1: Previsió de les exportacions de béns i serveis d’Algèria mitjançant suavització exponencial simple.

Any	Temps	Observació	Nivell	Previsió
	\(t\)	\(y_t\)	\(ℓ_t\)	\(\hat y_{t\|t−1}\)
1959	0		39.54
1960	1	39.04	39.12	39.54
1961	2	46.24	45.1	39.12
1962	3	19.79	23.84	45.1
1963	4	24.68	24.55	23.84
1964	5	25.08	25	24.55
1965	6	22.6	22.99	25
1966	7	25.99	25.51	22.99
1967	8	23.43	23.77	25.51
⋮		⋮	⋮	⋮
2014	55	30.22	30.8	33.85
2015	56	23.17	24.39	30.8
2016	57	20.86	21.43	24.39
2017	58	22.64	22.44	21.43
	\(h\)			\(\hat y_{T+h\|T}\)
2018	1			22.44
2019	1			22.44
2020	1			22.44
2021	1			22.44
2022	1			22.44

plot(forecast(ajustalg, h=10),
     col = "dodgerblue3", fcol= "red", panel.first=grid(), flwd = 2.5,lwd=2.5,
     main = "Exportacions Algèria: històric i previsió", ylab = "% del PIB")
lines(ajustalg$fitted, col= "darkred",lwd= 2)
legend("bottomleft", title = "Intervals de confiança",
       legend=c("80%", "95%"), fill = c("grey65","grey85"),
       bg="white",cex = 0.7, box.lty = NULL)
legend("bottom", title = "Línees",
       legend=c("Valors bruts", "Valors ajustats", "Previsió"), 
       fill = c("dodgerblue3","darkred","red"),
       bg="white",cex = 0.7, box.lty = NULL)

Figura 8.2

Les previsions per al període 2018–2022 es representen a la Figura 8.2. També es representen els valors ajustats un pas per endavant juntament amb les dades durant el període 1960–2017. El gran valor de \(α\) en aquest exemple es reflecteix en el gran ajust que té lloc en el nivell estimat \(ℓ_t\) en cada moment. Un valor més petit d’\(α\) comportaria canvis més petits al llarg del temps, i per tant la sèrie de valors ajustats seria més suau.

Els intervals de predicció mostren que hi ha una incertesa considerable en les exportacions futures durant el període de previsió de cinc anys. Per tant, interpretar les previsions puntuals sense tenir en compte la gran incertesa pot ser molt enganyós.

2. Mètodes amb tendència

Método de tendencia lineal de Holt

Holt (1957) va ampliar el suavització exponencial simple per permetre la previsió de dades amb una tendència. Aquest mètode implica una equació de previsió i dues equacions de suavització (una per al nivell i una per a la tendència):

\[ \begin{align*} 1. &&\text{Forecast equation}&& \hat{y}_{t+h|t} &= \ell_{t} + hb_{t} \\ 2. &&\text{Level equation} && \ell_{t} &= \alpha y_{t} + (1 - \alpha)(\ell_{t-1} + b_{t-1})\\ 3.&&\text{Trend equation} && b_{t} &= \beta^*(\ell_{t} - \ell_{t-1}) + (1 -\beta^*)b_{t-1} \end{align*}\tag{8.3} \] on \(ℓ_t\) denota una estimació del nivell de la sèrie en el temps \(t\), \(b_t\) denota una estimació de la tendència (pendent) de la sèrie en el temps \(t\), \(α\) és el paràmetre de suavització per al nivell, \(0≤α≤1\), i \(β^∗\) és el paràmetre de suavització per a la tendència, \(0≤β^∗≤1\). (Ho denotem com \(aβ^∗\) en lloc de \(β\) per raons que s’explicaran en pròxims apartats.)

Igual que amb la suavització exponencial simple, l’equació de nivell mostra que \(ℓ_t\) és una mitjana ponderada de l’observació \(y_t\) i la previsió d’entrenament amb un pas d’antelació per al temps \(t\), donada aquí per \(ℓ_{t−1}+b_{t−1}\). L’equació de tendència mostra que \(b_t\) és una mitjana ponderada de la tendència estimada en el moment t basada en \(ℓ_t−ℓ_{t−1}\) i \(b_{t−1}\), l’estimació anterior de la tendència.

La funció de previsió ja no és plana sinó que té tendència. La previsió amb h passos per davant és igual a l’últim nivell estimat més h multiplicat per l’últim valor de tendència estimat. Per tant, les previsions són una funció lineal de h.

Exemple: la població d’Austràlia

aus_economy <- global_economy |>
  filter(Code == "AUS") |>
  mutate(Pop = Population / 1e6)
autoplot(aus_economy, Pop) +
  labs(y = "Millions", title = "Australian population")+
  theme_bw()

Transformarem les dades en format ts per poder-les treballar amb la funció específica de forecast:: per a suvització exponencial, executant amb el parámetre model="AAN" específic per a funció exponencial amb tendència de Holdt.

tsaus <- ts(aus_economy[,10], start = 1960, frequency = 1)

ajustaus <- ets(tsaus, model = "AAN")

ajustaus

## ETS(A,A,N) 
## 
## Call:
## ets(y = tsaus, model = "AAN")
## 
##   Smoothing parameters:
##     alpha = 0.9999 
##     beta  = 0.3266 
## 
##   Initial states:
##     l = 10.0541 
##     b = 0.2225 
## 
##   sigma:  0.0643
## 
##       AIC      AICc       BIC 
## -76.98569 -75.83184 -66.68347

El coeficient de suavització estimat per al nivell és \(\hat α=0,9999\). El valor molt alt mostra que el nivell canvia ràpidament per capturar la sèrie amb alta tendència. El coeficient de suavització estimat per al pendent és \(\hat β^∗=0,3267\). Això és relativament gran, cosa que suggereix que la tendència també canvia sovint (fins i tot si els canvis són lleus).

Any	Temps	Observació	Nivell	Pendent	Previsió
	\(t\)	\(y_t\)	\(ℓ_t\)		\(y_{t+1\|t}\)
1959	0		10.05	0.22
1960	1	10.28	10.28	0.22	10.28
1961	2	10.48	10.48	0.22	10.50
1962	3	10.74	10.74	0.23	10.70
1963	4	10.95	10.95	0.22	10.97
1964	5	11.17	11.17	0.22	11.17
1965	6	11.39	11.39	0.22	11.39
1966	7	11.65	11.65	0.23	11.61
⋮		⋮		⋮	⋮
2014	55	23.50	23.50	0.37	23.52
2015	56	23.85	23.85	0.36	23.87
2016	57	24.21	24.21	0.36	24.21
2017	58	24.60	24.60	0.37	24.57
	\(h\)				\(\hat y_{T+h∣T}\)
2018	1				24.97
2019	2				25.34
2020	3				25.71
2021	4				26.07
2022	5				26.44
2023	6				26.81
2024	7				27.18
2025	8				27.55
2026	9				27.92
2027	10				28.29

On apliquem les 3 equacions que acabem de veure (8.3). Per exemple, la primera fila serà (1960, Temps=1):

l0 = 10.05; b0 = 0.22
alpha=0.9999; beta=0.3267

l1 <- alpha*10.28+(1-alpha)*(l0+b0)
b1 <- beta*(l1-l0)+(1-beta)*b0
for1 <- l0+1*b0

i la segona fila

l2 <- alpha*10.48+(1-alpha)*(l1+b1)
b2 <- beta*(l2-l1)+(1-beta)*b1
for2 <- l1+1*b1

així que

data.frame(t= c(1,2),
           l_t= c(l1,l2),
           b_t= c(b1,b2),
           forcast = c(for1, for2))

##   t   l_t       b_t  forcast
## 1 1 10.28 0.2232667 10.27000
## 2 2 10.48 0.2156665 10.50327

… i etcétera.

Mètodes de tendència amortida

Les previsions generades pel mètode lineal de Holt mostren una tendència constant (augmentant o decreixent) indefinidament en el futur. L’evidència empírica indica que aquests mètodes tendeixen a pronosticar excessivament, especialment per a horitzons de previsió més llargs. Motivats per aquesta observació, Gardner i McKenzie (1985) van introduir un paràmetre que “esmorteeix” la tendència a una línia plana en algun moment en el futur. Els mètodes que inclouen una tendència esmorteïda han demostrat tenir molt d’èxit i, sens dubte, són els mètodes individuals més populars quan es requereixen previsions automàticament per a moltes sèries.

En conjunció amb els paràmetres de suavització \(α\) i \(β^∗\) (amb valors entre 0 i 1 com en el mètode de Holt), aquest mètode també inclou un paràmetre d’amortiment \(0<ϕ<1\) :

\[ \begin{align*} \hat{y}_{t+h|t} &= \ell_{t} + (\phi+\phi^2 + \dots + \phi^{h})b_{t} \\ \ell_{t} &= \alpha y_{t} + (1 - \alpha)(\ell_{t-1} + \phi b_{t-1})\\ b_{t} &= \beta^*(\ell_{t} - \ell_{t-1}) + (1 -\beta^*)\phi b_{t-1}. \end{align*} \]

Si ϕ=1, el mètode és idèntic al mètode lineal de Holt. Per a valors entre 0 i 1 , \(ϕ\) esmorteeix la tendència de manera que s’acosta a una constant en el futur. De fet, les previsions convergeixen a \(ℓ_T+ϕb_T/(1−ϕ)\) com a \(h→∞\) per a qualsevol valor \(0<ϕ<1\). Això vol dir que les previsions a curt termini tenen tendència mentre que les previsions a llarg termini són constants.

A la pràctica, \(ϕ\) rarament és inferior a 0,8, ja que l’amortiment té un efecte molt fort per a valors més petits. Els valors de \(ϕ\) propers a 1 significaran que un model amortit no es pot distingir d’un model no amortit. Per aquests motius, normalment restringim \(ϕ\) a un mínim de 0,8 i un màxim de 0,98.

Exemple continuació de població australiana

Ara compararem els models suavització amb trend linial i amb trend amorti.

ajustausdam <- ets(tsaus, 
                model = "AAN", 
                damped = T,
                phi=0.9)
holdt <- forecast::forecast(ajustaus,15)
damped <- forecast::forecast(ajustausdam,15)

plot(ts(aus_economy[,10], start = 1960, frequency = 1), ylim= c(10,30), xlim = c(1960,2032),
     lwd=1.8, col="dodgerblue3", panel.first=grid(), xlab="",ylab="", 
     main= "Població d'Austràlia i forecasts \namb trend linial i trend amortit (Mètode Holt)")
lines(holdt$mean, col= "darkred",lwd= 2)
lines(damped$mean, col= "red", lwd=2)
legend("topleft", title = "Models de predicció",
       legend=c("Dades brutes", "Trend linial", "Tend amortit"), 
       fill = c("dodgerblue3","darkred","red"),
       bg="white",cex = 0.7, box.lty = NULL)

Hem establert el paràmetre d’amortiment en un nombre relativament baix \((ϕ=0,90)\) per exagerar l’efecte de l’amortiment per comparar-lo. Normalment, estimaríem \(ϕ\) juntament amb els altres paràmetres. També hem utilitzat un horitzó de previsió força gran ( \(h=15\) ) per destacar la diferència entre una tendència amortida i una tendència lineal.

Exemple: ús d’Internet

En aquest exemple, comparem el rendiment de previsió dels tres mètodes de suavització exponencial que hem considerat fins ara en la previsió del nombre d’usuaris connectats a Internet mitjançant un servidor. Les dades s’observen durant 100 minuts i es mostren a la figura 8.5.

www_usage <- as_tsibble(WWWusage)
www_usage |> autoplot(value) +
  labs(x="Minut", y="Nombre d'usuaris",
       title = "Ús d'internet per minut")+
  theme_bw()

Validació creant creant un apartat train de la sèrie

WWWtrain <- WWWusage[1:80]
WWWtest <- WWWusage[81:10]

# Models de previsió

ses <- forecast(forecast::ets(WWWtrain, model = "ANN"),h=20)
trend <- forecast(forecast::ets(WWWtrain, model = "AAN"),h=20)
holt_rect <- forecast(forecast::holt(WWWtrain, damped = F,h=20),h=20)
amort <- forecast(forecast::ets(WWWtrain, model = "AAN", damped = T, phi = 0.95),h=20)

pred_www <- list(
ses, trend, holt_rect,amort
)

accur_www <- lapply(pred_www, function(x) {
  as.data.frame(forecast::accuracy(
  x, WWWusage))
  }) |>
  bind_rows() |>
  select(2,3,5,6) |>
  slice(2,4,6,8)

accur_www$Mètode <- c("SES", "Holt no linial","Holt linial", "Damped (φ=0.95)")  


accur_www[,c(5,1:4)] |>
  flextable::flextable() |>
  flextable::set_caption("Diferents mesures de la precisió de les previsions respecte a les dades reals") %>% 
  flextable::width(j=1,width = 2)

Diferents mesures de la precisió de les previsions respecte a les dades reals
Mètode	RMSE	MAE	MAPE	MASE
SES	82.81817	76.200984	38.881010	18.695273
Holt no linial	53.57516	47.316924	23.683118	11.608811
Holt linial	28.14652	20.505873	9.888785	5.030944
Damped (φ=0.95)	10.56447	9.412866	4.999946	2.309368

plot(WWWusage, xlim = c(1, 105),
     lwd=1.8, col="dodgerblue", panel.first=grid(), xlab="",ylab="", 
     main= "Població d'Austràlia i forecasts \namb trend linial i trend amortit (Mètode Holt)")
lines(ses$mean, col= "darkred",lwd= 2)
lines(trend$mean, col= "red", lwd=2)
lines(holt_rect$mean, col= "violet", lwd=2)
lines(amort$mean, col= "blue", lwd=2)

legend("topleft", title = "Models de predicció",
       legend=c("Dades brutes", "SES: 'ANN'", "Holt no linial: 'AAN'", "Holt linial: 'holt()'", "Holt amortit: 'AAN' i 'φ=0.95'"), 
       fill = c("dodgerblue","darkred","red", "violet","blue"),
       bg="white",cex = 0.7, box.lty = NULL)

Cal fixar-se que la **funció ets(model="AAN", damped=F) produeix en aquest cas una línia corbada*. Això perquè el mètode lineal de Holt (‘AAN’) produeix una línia recta, però la funció ets() reestima automàticament els paràmetres d’estat intern (com el nivell i la tendència)** basant-se en les dades d’entrenament, sovint resultant en una curvatura suau i no lineal a la trama en lloc d’una extrapolació lineal estricta.

Components additius a l’espai d’estats: mentre que el mètode lineal de Holt (error additiu, tendència additiva) utilitza tècnicament una fórmula lineal (\(\hat{y}_{t+h|t} = \ell_{t}\)), el marc d’espai d’estats (“Levele equation i Trend equation”) de ets() a R forecast:: calcula el nivell òptim (\(\ell\)) i la tendència (\(b\)) que s’actualitzen a partir de les dades. El component de tendència resultant (\(b\)) pot no ser completament pla a tot l’horitzó de previsió, donant lloc a una corba.
La funció ets() utilitza mètodes específics per inicialitzar el nivell i la tendència inicials, que poden crear una curvatura en els passos immediats de previsió del futur, en lloc d’una línia perfectament recta des de l’últim punt observat.
No és una funció lineal “simple”: si utilitzeu la funció manual holt() amb damped = F, la línia seria perfectament recta. En canvi, ets() optimitza els paràmetres de suavització (\(\alpha\) i \(\beta\)) per a tota la sèrie, fent que el “pendent” sigui una mitjana ponderada dels canvis recents, que sovint es corba per seguir els darrers punts de dades més de prop.

Validació creuada

Utilitzarem la validació creuada de sèries temporals per comparar la precisió de la previsió en un sol pas dels tres mètodes.

www_usage |>
  stretch_tsibble(.init = 10) |>
  model(
    SES = ETS(value ~ error("A") + trend("N") + season("N")),
    Holt = ETS(value ~ error("A") + trend("A") + season("N")),
    Damped = ETS(value ~ error("A") + trend("Ad") +
                   season("N"))
  ) |>
  forecast(h = 1) |>
  accuracy(www_usage) |>
  select(1,2,4,5,7,8) |>
  flextable::flextable() |>
  flextable::set_caption("Diferents mesures de la precisió de les previsions respecte a les dades reals, obtingudes per Cross Validation")

Diferents mesures de la precisió de les previsions respecte a les dades reals, obtingudes per Cross Validation
.model	.type	RMSE	MAE	MAPE	MASE
Damped	Test	3.686810	2.999724	2.264128	0.6628854
Holt	Test	3.874378	3.174028	2.382212	0.7014035
SES	Test	6.049937	4.813149	3.548391	1.0636199

S’estima que el paràmetre de suavització del pendent és gairebé 1, cosa que indica que la tendència canvia per reflectir majoritàriament el pendent entre els dos últims minuts d’ús d’Internet. El valor de \(α\) és molt proper a 1, cosa que demostra que el nivell (\(\ell\)) reacciona fortament a cada nova observació.

Les previsions resultants semblen sensates amb una tendència decreixent, que s’aplana a causa del baix valor del paràmetre d’amortiment (\(\phi=0,815\)) i dels intervals de predicció relativament amplis que reflecteixen la variació de les dades històriques. Els intervals de predicció es calculen mitjançant els mètodes descrits més endavant.

En aquest exemple, el procés de selecció d’un mètode és relativament fàcil, ja que les comparacions de MSE i MAE van suggerir el mateix mètode (Holt amortit). No obstant això, de vegades diferents mesures de precisió suggeriran mètodes de previsió diferents, i llavors cal decidir quin mètode de previsió preferim utilitzar. Com que les tasques de previsió poden variar en moltes dimensions (longitud de l’horitzó de previsió, mida del conjunt de proves, mesures d’error de previsió, freqüència de dades, etc.), és poc probable que un mètode sigui millor que tots els altres per a tots els escenaris de previsió. El que necessitem d’un mètode de previsió són previsions coherentment raonables, i aquestes s’han d’avaluar amb freqüència en funció de la tasca en qüestió.

3. Mètodes amb estacionalitat

El mètode estacional de Holt-Winters inclou l’equació de previsió i tres equacions de suavització: una per al nivell \(ℓ_t\), una per a la tendència \(b_t\) i una per al component estacional \(s_t\), amb els paràmetres de suavització corresponents \((α, β^∗, γ)\). Utilitzem \(m\) per indicar el període de l’estacionalitat, és a dir, el nombre d’estacions en un any. Per exemple, per a dades trimestrals \(m=4\) i per a dades mensuals \(m=12\).

Hi ha dues variants d’aquest mètode que es diferencien en la naturalesa del component estacional. Es prefereix el mètode additiu quan les variacions estacionals són aproximadament constants a través de la sèrie, mentre que el mètode multiplicatiu es prefereix quan les variacions estacionals canvien proporcionalment al nivell de la sèrie.

Amb el mètode additiu*, la component estacional s’expressa en termes absoluts a l’escala de la sèrie observada, i en l’equació de nivells la sèrie s’ajusta estacionalment restant la component estacional. Dins de cada any, el component estacional sumarà aproximadament zero.
Amb el mètode multiplicatiu, el component estacional s’expressa en termes relatius (percentatges) i la sèrie s’ajusta estacionalment dividint-lo pel component estacional. Dins de cada any, el component estacional sumarà aproximadament \(m\).

Mètode additu de Holt-Winters

La forma dels components per al mètode additiu és:

\[ \begin{align*} \hat{y}_{t+h|t} &= \ell_{t} + hb_{t} + s_{t+h-m(k+1)} \\ \ell_{t} &= \alpha(y_{t} - s_{t-m}) + (1 - \alpha)(\ell_{t-1} + b_{t-1})\\ b_{t} &= \beta^*(\ell_{t} - \ell_{t-1}) + (1 - \beta^*)b_{t-1}\\ s_{t} &= \gamma (y_{t}-\ell_{t-1}-b_{t-1}) + (1-\gamma)s_{t-m}, \end{align*} \]

\(k\) és la part sencera de \((h−1)/m\), que assegura que les estimacions dels índexs estacionals utilitzats per a la predicció provenen de l’últim any de la mostra.
L’equació de nivell mostra una mitjana ponderada entre l’observació desestacionalitzada \((y_t−s_{t−m})\) i la previsió no estacional \((ℓ_{t−1}+b_{t−1})\) per al temps \(t\) .
L’equació de tendència és idèntica al mètode lineal de Holt.
L’equació estacional mostra una mitjana ponderada entre l’índex estacional actual, \((y_t−ℓ_{t−1}−b_{t−1})\) i l’índex estacional de la mateixa temporada de l’any passat (és a dir, fa m períodes de temps).

L’equació del component estacional sovint s’expressa com

\[ s_{t} = \gamma^* (y_{t}-\ell_{t})+ (1-\gamma^*)s_{t-m}. \]

Si substituïm \(\ell_t\) amb l’equació del nivell, obtenim

\[ s_{t} = \gamma^*(1-\alpha) (y_{t}-\ell_{t-1}-b_{t-1})+ [1-\gamma^*(1-\alpha)]s_{t-m}, \]

que és identica a l’equció de suavització per al component estacional que acabem de veure, amb \(γ=γ∗(1−α)\). El paràmetre de restricció és \(0≤γ∗≤1\), qe dona lloc a \(0≤γ≤1−α\).

Mètode multiplicatiu de Holt-Winters

La forma dels components del mètode és

\[ \begin{align*} \hat{y}_{t+h|t} &= (\ell_{t} + hb_{t})s_{t+h-m(k+1)} \\ \ell_{t} &= \alpha \frac{y_{t}}{s_{t-m}} + (1 - \alpha)(\ell_{t-1} + b_{t-1})\\ b_{t} &= \beta^*(\ell_{t}-\ell_{t-1}) + (1 - \beta^*)b_{t-1} \\ s_{t} &= \gamma \frac{y_{t}}{(\ell_{t-1} + b_{t-1})} + (1 - \gamma)s_{t-m}. \end{align*} \]

Exemple: Viatges nacionals amb pernoctació a Austràlia

Apliquem el mètode de Holt-Winters amb estacionalitat additiva i multiplicativa17 per predir trimestralment les nits de visita a Austràlia que els turistes nacionals passen. La figura 8.7 mostra les dades del període 1998-2017 i les previsions del període 2018-2020. Les dades mostren un patró estacional evident, amb pics observats al trimestre de març de cada any, corresponents a l’estiu australià.

## opt_crit="mse"

aus_holidays <- tourism |>
  filter(Purpose == "Holiday") |>
  summarise(Trips = sum(Trips)/1e3)
fit <- aus_holidays |>
  model(
    additive = ETS(Trips ~ error("A") + trend("A") +
                                                season("A")),
    multiplicative = ETS(Trips ~ error("M") + trend("A") +
                                                season("M"))
  )
fc <- fit |> forecast(h = "3 years")
fc |>
  autoplot(aus_holidays, level = NULL) +
  labs(title="Australian domestic tourism (fable package)",
       y="Overnight trips (millions)") +
  guides(colour = guide_legend(title = "Forecast"))

# Conversió a format ts
tsaus_holidays <- ts(aus_holidays[,2], start = c(1998,1), frequency = 4)

# Model additiu i model multiplicatiu
addmod <- forecast::ets(tsaus_holidays, model = "AAA")
multmod <- forecast::ets(tsaus_holidays, model = "MAM")

# Previsió
fcadmod <-  addmod %>% forecast::forecast(h=12)
fcmultmod <-  multmod |> forecast::forecast(h=12)

# Gràfic
plot(tsaus_holidays, xlim = c(1997, 2021), ylim= c(7.5,14),
     lwd=1.8, col="dodgerblue", panel.first=grid(), xlab="",ylab="", 
     main= "Viatges nacionals amb pernoctació a Austràlia (ts() i forecast::)")
lines(fcadmod $mean, col= "blue",lwd= 2.5)
lines(fcmultmod$mean, col= "red", lwd=2.5)

legend("topleft", title = "Models de predicció",
       legend=c("Dades brutes", "Mèt. Addit.: 'AAA'", "Mèt. Multipl.: 'MAM'"), 
       fill = c("dodgerblue","blue","red"),
       bg=NULL,cex = 0.7, box.lty = NULL)

library(patchwork)


add <- ets(tsaus_holidays, model = "AAA") %>%  
  autoplot() +
  theme_bw()

mul <-ets(tsaus_holidays, model = "MAM") %>%  
  autoplot() +
  theme_bw()

add +mul

Com que ambdós mètodes tenen exactament el mateix nombre de paràmetres per estimar, podem comparar l’RMSE d’entrenament d’ambdós models. En aquest cas, el mètode amb estacionalitat multiplicativa s’ajusta lleugerament millor a les dades.

eval <- rbind(forecast::accuracy(multmod),
forecast::accuracy(addmod))
as.data.frame(eval) %>%  mutate(Mètode = c("Mulitplicatiu", "Additiu")) %>% select(8,2,3,5,6) %>% 
  flextable::flextable()

Mètode	RMSE	MAE	MAPE	MASE
Mulitplicatiu	0.4121669	0.3095759	3.248095	0.7467908
Additiu	0.4168762	0.3209553	3.377814	0.7742412

Mètode Holt-Winters amb amortiment

L’amortiment és possible tant amb mètodes additius com multiplicatius de Holt-Winters. Un mètode que sovint proporciona previsions precises i robustes per a dades estacionals és el mètode Holt-Winters amb una tendència amortida i estacionalitat multiplicativa:

\[ \begin{align*} \hat{y}_{t+h|t} &= \left[\ell_{t} + (\phi+\phi^2 + \dots + \phi^{h})b_{t}\right]s_{t+h-m(k+1)} \\ \ell_{t} &= \alpha(y_{t} / s_{t-m}) + (1 - \alpha)(\ell_{t-1} + \phi b_{t-1})\\ b_{t} &= \beta^*(\ell_{t} - \ell_{t-1}) + (1 - \beta^*)\phi b_{t-1} \\ s_{t} &= \gamma \frac{y_{t}}{(\ell_{t-1} + \phi b_{t-1})} + (1 - \gamma)s_{t-m}. \end{align*} \]

El mètode Holt-Winters també es pot utilitzar per a dades de tipus diari, on el període estacional és m=7, i la unitat de temps adequada per a h són els dies.

hwd <- forecast(ets(window(hyndsight, end= c(48,4)), model = "MAM", damped = T), h= 35)

autoplot(hyndsight, series = "Dades originals") +
  forecast::autolayer(hwd$mean, series = "Holt-Winters amb amortiment")+
  ggtitle("Sèrie original diària i ajustada am Mètode Holt-Winters \namb amortiment (forecast + autoplot)") +
  theme_bw()

plot(hyndsight, #xlim = c(0,53, 2021), ylim= c(7.5,14),
     lwd=1.8, col="dodgerblue3", panel.first=grid(), xlab="",ylab="", 
     main= "Sèrie original diària i ajustada am Mètode Holt-Winters \namb amortiment (forecast.plot)")
lines(hwd$mean, col= "red",lwd= 1.8)

Com que tots dos mètodes tenen exactament el mateix nombre de paràmetres per estimar, podem comparar l’RMSE d’entrenament dels dos models. En aquest cas, el mètode amb estacionalitat multiplicativa s’ajusta lleugerament millor a les dades.

El mètode Holt-Winters també es pot utilitzar per a dades de tipus diari, on el període estacional és m=7, i la unitat de temps adequada per a h és en dies.

4. Innovations state space models per a suavització exponencial

A la resta d’aquest capítol, estudiem els models estadístics que subjauen als mètodes de suavització exponencial que hem considerat fins ara. Són algoritmes que generen previsions puntuals. Els models estadístics d’aquesta secció generen les mateixes previsions puntuals, però també poden generar intervals de predicció (o previsió). Un model estadístic és un procés de generació de dades estocàstiques (o aleatòries) que pot produir una distribució de previsió completa. El marc estadístic general que introduirem també proporciona una manera d’utilitzar els criteris de selecció de models introduïts al permetent així que l’elecció del model es faci de manera objectiva.

Cada model consta d’una equació de mesura que descriu les dades observades i algunes equacions de transició que descriuen com els components o estats no observats (nivell, tendència, estacional) canvien amb el temps. Per tant, aquests es denominen “models d’espai d’estat”.

Per a cada mètode existeixen dos models: un amb errors additius i un amb errors multiplicatius. Les previsions puntuals produïdes pels models són idèntiques si utilitzen els mateixos valors dels paràmetres de suavització. Tanmateix, generaran intervals de predicció diferents.

Per distingir entre un model amb errors additius i un amb errors multiplicatius (i també per distingir els models dels mètodes), afegim una tercera lletra a la classificació de la Taula 7.6. Etiquetem cada model d’espai d’estat com a ETS(⋅,⋅,⋅) per a** (Error, Tendència, Estacional). Aquesta etiqueta també es pot considerar com a Suavització Exponencial. Utilitzant la mateixa notació que a la Taula 7.6, les possibilitats per a cada component són: \(Error =\{ A,M \}\), \(Tendència =\{ N,A,Ad \}\) i \(Estacional =\{ N,A,M \}\).

ETS(A,N,N): suavització exponencial simple amb errors additius

\[ \begin{align*} \text{Forecast equation} && \hat{y}_{t+1|t} & = \ell_{t}\\ \text{Smoothing equation} && \ell_{t} & = \alpha y_{t} + (1 - \alpha)\ell_{t-1}, \end{align*} \]

Si reorganitzem l’equació de suavització per al nivell, obtenim la forma de correcció d’errors:

\[ \begin{align*} \ell_{t} &= \alpha y_{t}+\ell_{t-1}-\alpha\ell_{t-1}\\ &= \ell_{t-1}+\alpha( y_{t}-\ell_{t-1})\\ &= \ell_{t-1}+\alpha e_{t} \end{align*} \]

On \(e_{t}=y_{t}-\ell_{t-1}=y_{t}-\hat{y}_{t|t-1}\) per a \(t=1,...,T\). És a dir, et és l’error d’un pas en el temps t calculat sobre les dades d’entrenament.

Els errors de les dades d’entrenament condueixen a l’ajust del nivell estimat durant tot el procés de suavització per a \(t=1,…,T\). Per exemple, si l’error en el temps \(t\) és negatiu, aleshores \(\hat y_{t|t−1}>y_t\) i, per tant, el nivell en el temps \(t−1\) s’ha sobreestimat. El nou nivell \(ℓ_t\) és aleshores el nivell anterior \(ℓ_t−1\) ajustat a la baixa. Com més a prop estigui \(α\) d’1, més “aproximada” serà l’estimació del nivell (es produiran grans ajustos). Com més petit sigui l’\(α\), més “suau” serà el nivell (es produeixen petits ajustaments).

També podem escriure \(y_t=ℓ_{t−1}+e_t\), de manera que cada observació sigui igual al nivell anterior més un error. Per convertir això en un model d’espai d’estats d’innovacions, tot el que hem de fer és especificar la distribució de probabilitat per a \(e_t\). Per a un model amb errors additius, suposem que els errors de previsió d’un pas et són soroll blanc distribuït normalment amb mitjana 0 i variància \(σ^2\). Una notació abreujada per a això és \(e_t=ε_t∼\text{NID}(0,σ^2)\); NID significa “distribuït normalment i independentment”.

Aleshores, les equacions del model es poden escriure com \[ \begin{align} y_t &= \ell_{t-1} + \varepsilon_t \tag{7.3}\\ \ell_t&=\ell_{t-1}+\alpha \varepsilon_t. \tag{7.4} \end{align} \]

Ens referim a (7.3) com l’equació de mesura (o observació) i a (7.4) com l’equació d’estat (o transició). Aquestes dues equacions, juntament amb la distribució estadística dels errors, formen un model estadístic completament especificat. Concretament, constitueixen un model d’espai d’estat d’innovacions subjacent a un suavització exponencial simple.

El terme “innovacions” prové del fet que totes les equacions d’aquest tipus d’especificació utilitzen el mateix procés d’error aleatori, \(ε_t\). Per la mateixa raó, aquesta formulació també es coneix com a model de “font única d’error”, en contrast amb la formulació alternativa de fonts múltiples d’error, que no presentem aquí.

L’equació de mesura mostra la relació entre les observacions i els estats no observats. En aquest cas, l’observació \(y_t\) és una funció lineal del nivell \(ℓ_{t−1}\), la part predictible de \(y_t\), i l’error aleatori \(ε_t\), la part imprevisible de \(y_t\). Per a altres models d’espai d’estat d’innovacions, aquesta relació pot ser no lineal.

L’equació de transició mostra l’evolució de l’estat al llarg del temps. La influència del paràmetre de suavització \(α\) és la mateixa que per als mètodes comentats anteriorment. Per exemple, \(α\) regeix el grau de canvi en nivells successius. Com més alt sigui el valor d’\(α\), més ràpids seran els canvis en el nivell; com més baix sigui el valor d’\(α\), més suaus seran els canvis. A l’extrem més baix, on \(α = 0\), el nivell de la sèrie no canvia amb el temps. A l’altre extrem, on \(α = 1\), el model es redueix a un model de passeig aleatori, \(y_t = y_{t − 1} + ε_t\).

ETS(M,N,N): suavització exponencial simple amb errors mmultiplicatius

\[ \begin{align*} y_t&=\ell_{t-1}(1+\varepsilon_t)\\ \ell_t&=\ell_{t-1}(1+\alpha \varepsilon_t). \end{align*} \]

ETS(A,A,N): Mètode linial de Holt amb errors additius

\[ \begin{align*} y_t&=\ell_{t-1}+b_{t-1}+\varepsilon_t\\ \ell_t&=\ell_{t-1}+b_{t-1}+\alpha \varepsilon_t\\ b_t&=b_{t-1}+\beta \varepsilon_t, \end{align*} \]

Per simplificar, \(\beta\) equival a \(\alpha\beta^*\)

ETS(M,A,N): mètode linia de Holt amb errors multiplicatius

\[ \begin{align*} y_t&=(\ell_{t-1}+b_{t-1})(1+\varepsilon_t)\\ \ell_t&=(\ell_{t-1}+b_{t-1})(1+\alpha \varepsilon_t)\\ b_t&=b_{t-1}+\beta(\ell_{t-1}+b_{t-1}) \varepsilon_t \end{align*} \]

La possibilitat de combinacions és en tot cas gran. Per això veure la taula a Forecasting Principles and Practices ed 2ª.

5. Estimació i selecció de models

Estimació de models ETS

Una alternativa a estimar els paràmetres minimitzant la suma d’errors quadràtics és maximitzar la “verosimilitud”.

La la verosimilitud és la probabilitat que les dades provinguin del model especificat. Per tant, una gran probabilitat s’associa amb un bon model. Per a un model d’error additiu, maximitzar la verosimilitud dóna els mateixos resultats que minimitzar la suma d’errors quadràtics. Tanmateix, s’obtindran resultats diferents per a models d’error multiplicatius. En aquesta secció, estimarem els paràmetres de suavització \(α, β, γ, ϕ\), i els estats inicials \(ℓ_0, b_0, s_0, s_{−1},…, s−m+1\), maximitzant la probabilitat.

Els possibles valors que poden prendre els paràmetres de suavització són restringits. Tradicionalment, els paràmetres s’han restringit a estar entre 0 i 1 perquè les equacions es puguin interpretar com a mitjanes ponderades. És a dir, \(0<α, β^∗, γ^∗, ϕ<1\).

Per als models d’espai d’estats, hem establert \(β=αβ^∗\) i \(γ=(1−α)γ^∗\). Per tant, les restriccions tradicionals es tradueixen a \(0<α<1\), \(0<β<α\) i \(0<γ<1−α\).

A la pràctica, el paràmetre d’amortiment \(ϕ\) normalment es restringeix encara més per evitar dificultats numèriques en l’estimació del model. A R, es restringeix de manera que \(0.8<ϕ<0.98\).

Una altra manera de veure els paràmetres és mitjançant una consideració de les propietats matemàtiques dels models d’espai d’estats. Els paràmetres estan restringits per evitar que les observacions en el passat llunyà tinguin un efecte continu sobre les previsions actuals. Això condueix a algunes restriccions d’admissibilitat sobre els paràmetres, que solen ser (però no sempre) menys restrictives que la regió habitual. Per exemple, per al model ETS(A,N,N), la regió de paràmetres habitual és 0<α<1, però la regió admissible és 0<α<2. Per al model ETS(A,A,N), la regió de paràmetres habitual és 0<α<1 i 0<β<α, però la regió admissible és 0<α<2 i 0<β<4−2α.

Selecció de models

Un gran avantatge del marc estadístic ETS¹ és que es poden utilitzar criteris d’informació per a la selecció de models. L’\(\text{AIC}\), l’\(\text{AIC}\text{_c}\) i el \(\text{BIC}\), introduïts a la secció 5.5, es poden utilitzar aquí per determinar quin dels models ETS és el més adequat per a una sèrie temporal determinada. \[ \text{AIC} = -2\log(L) + 2k, \]

on \(L\) és la verosimilitud del model i \(k\) és el nombre total de paràmetres i estats inicials que s’han estimat (inclosa la variància residual).

L’\(\text{AIC}\) corregit per biaix de mostra petita (\(\text{AIC}\text{_c}\)) es defineix com \[ \text{AIC}_{\text{c}} = \text{AIC} + \frac{k(k+1)}{T-k-1}, \]

i el Criteri d’informació bayesana, \(\text{BIC}\)

\[ \text{BIC} = \text{AIC} + k[\log(T)-2]. \]

Tres de les combinacions de (Error, Tendència, Estacional) poden provocar dificultats numèriques. Concretament, els models que poden causar aquestes inestabilitats són ETS(A,N,M), ETS(A,A,M) i ETS(A,Ad,M). Normalment no tenim en compte aquestes combinacions particulars a l’hora de seleccionar un model.

Els models amb errors multiplicatius són útils quan les dades són estrictament positives, però no són numèricament estables quan les dades contenen zeros o valors negatius. Per tant, els models d’error multiplicatiu no es consideraran si la sèrie temporal no és estrictament positiva. En aquest cas, només s’aplicaran els sis models totalment additius.

La funció `ets()` a R

Els models es poden estimar amb R utilitzant la funció ets() del paquet forecast. A diferència de les funcions ses, holt i hw, la funció ets() no produeix previsions. En canvi, estima els paràmetres del model i retorna informació sobre el model ajustat.

El codi R següent mostra els arguments més importants que pot prendre aquesta funció i els seus valors per defecte. Si només s’especifica la sèrie temporal i tots els altres arguments es deixen als seus valors per defecte, es seleccionarà automàticament un model adequat. Expliquem els arguments a continuació. Vegeu el fitxer d’ajuda per obtenir una descripció més completa.

ets(y, model="ZZZ", damped=NULL, alpha=NULL, beta=NULL, gamma=NULL,
    phi=NULL, lambda=NULL, biasadj=FALSE, additive.only=FALSE,
    restrict=TRUE, allow.multiplicative.trend=FALSE)

`y`

La sèrie temporal que es preveu.

`model`

Un codi de tres lletres que indica el model que s’estimarà utilitzant la classificació i la notació ETS. Les possibles entrades són “N” per a cap, “A” per a additiu, “M” per a multiplicatiu o “Z” per a la selecció automàtica. Si alguna de les entrades es deixa com a “Z”, aquest component es selecciona segons el criteri d’informació escollit. El valor per defecte de ZZZ garanteix que tots els components es seleccionin mitjançant el criteri d’informació.

`damped`

Si damped=TRUE, s’utilitzarà una tendència amortida (A o M). Si damped=FALSE, s’utilitzarà una tendència no amortida. Si damped=NULL (el valor per defecte), es seleccionarà una tendència amortida o no amortida, depenent de quin model tingui el valor més petit per al criteri d’informació.

`alpha, beta, gamma, phi`

Els valors dels paràmetres de suavització es poden especificar mitjançant aquests arguments. Si es defineixen com a NULL (el valor per defecte per a cadascun d’ells), els paràmetres s’estimen.

lambda

Paràmetre de transformació de Box-Cox. S’ignorarà si lambda=NULL (el valor per defecte). En cas contrari, la sèrie temporal es transformarà abans que s’estimi el model. Quan lambda no és NULL, additive.only s’estableix com a TRUE.

`biasadj`

Si TRUE i lambda no és NULL, els valors ajustats i les previsions transformades enrere s’ajustaran al biaix.

`additive.only`

Només es consideraran els models amb components additius si additive.only=TRUE. En cas contrari, es consideraran tots els models.

`restrict`

Si restrict=TRUE (el valor per defecte), els models que causen dificultats numèriques no es consideren en la selecció del model.

`allow.multiplicative.trend`

També hi ha models de tendència multiplicativa disponibles. Establiu aquest argument a TRUE per permetre que aquests models es considerin.

Treball amb objectes ets

La funció ets() retornarà un objecte de la classe ets. Hi ha moltes funcions R dissenyades per facilitar el treball amb objectes ets. Algunes d’elles es descriuen a continuació.

`coef()`

retorna tots els paràmetres ajustats.

`accuracy()`

retorna les mesures de precisió calculades sobre les dades d’entrenament.

`summary()`

imprimeix informació resumida sobre el model ajustat.

`autoplot()` i `plot()`

produeixen diagrames de temps dels components.

`residuals()`

retorna els residuals del model estimat.

`fitted()`

retorna les previsions d’un pas per a les dades d’entrenament.

`simulate()`

simularà futures rutes de mostra del model ajustat.

`forecast()`

calcula prediccions puntuals i intervals de predicció, tal com es descriu a la secció següent.

Exemple: Nits de visitants turístics internacionals a Austràlia

Ara utilitzem el marc estadístic ETS per predir les nits de visitants turístics a Austràlia segons les arribades internacionals durant el període 2016–2019. Deixem que la funció ets() **seleccioni el model minimitzant l’_{}**.

aust <- window(fpp2::austourists, start=2005)

fit <- ets(aust)

summary(fit)

## ETS(M,A,M) 
## 
## Call:
## ets(y = aust)
## 
##   Smoothing parameters:
##     alpha = 0.1908 
##     beta  = 0.0392 
##     gamma = 2e-04 
## 
##   Initial states:
##     l = 32.3679 
##     b = 0.9281 
##     s = 1.0218 0.9628 0.7683 1.2471
## 
##   sigma:  0.0383
## 
##      AIC     AICc      BIC 
## 224.8628 230.1569 240.9205 
## 
## Training set error measures:
##                      ME     RMSE     MAE        MPE     MAPE     MASE      ACF1
## Training set 0.04836907 1.670893 1.24954 -0.1845609 2.692849 0.409454 0.2005962

El model seleccionat és ETS(M,A,M) \[ \begin{align*} y_{t} &= (\ell_{t-1} + b_{t-1})s_{t-m}(1 + \varepsilon_t)\\ \ell_t &= (\ell_{t-1} + b_{t-1})(1 + \alpha \varepsilon_t)\\ b_t &=b_{t-1} + \beta(\ell_{t-1} + b_{t_1})\varepsilon_t\\ s_t &= s_{t-m}(1+ \gamma \varepsilon_t). \end{align*} \]

Les estimacions dels paràmetres són \(α=0.1908\), \(β=0.03919\) i \(γ=0.0001917\). L’ouput també retorna les estimacions per als estats inicials \(ℓ_0, b_0, s_0, s_{−1}, s_{−2}, s_{−3}\).

Compareu-les amb els valors obtinguts per al mètode Holt-Winters amb estacionalitat multiplicativa presentats a la Taula 7.5. Tot i que aquest model dóna previsions puntuals equivalents a un mètode Holt-Winters multiplicatiu, els paràmetres s’han estimat de manera diferent. Amb la funció ets(), el mètode d’estimació per defecte és la màxima versemblança en lloc de la suma mínima de quadrats.

autoplot(fit) +
  theme_bw()

Com que aquest model té errors multiplicatius, els residuals no són equivalents als errors de previsió d’un pas.² Els residuals es donen per \(\hat ε_t\), mentre que els errors de previsió d’un pas es defineixen com a \(y_t−\hat y_{t−1}\). Podem obtenir tots dos utilitzant la funció de residuals.

cbind('Residuals' = residuals(fit),
      'Errors de previsió' = residuals(fit, type='response')) %>%
  autoplot() + xlab("Any") + ylab("") +
  guides(colour= guide_legend(title = ""))+
  theme_bw()

6. Pronòstic amb models ETS

Les previsions de punts dels models s’obtenen per iteració de les equacions per a \(t=T+1,...,T+h\). Per exemple, per al model ETS(M,A,N), será \[ y_{T+1} = (\ell_T + b_T )(1+ \varepsilon_{T+1}). \] Així, com que és un punt, \(\hat{y}_{T+1|T}=\ell_{T}+b_{T}.\)

Análogament, fent les substitucions pertinents \[ \begin{align*} y_{T+2} &= (\ell_{T+1} + b_{T+1})(1 + \varepsilon_{T+1})\\ &= \left[ (\ell_T + b_T) (1+ \alpha\varepsilon_{T+1}) + b_T + \beta (\ell_T + b_T)\varepsilon_{T+1} \right] ( 1 + \varepsilon_{T+1}). \end{align*} \] i anul·lant \(\epsilon\), serà \(\hat{y}_{T+2|T}= \ell_{T}+2b_{T},\) i aíxí seguint.

Aquestes previsions són idèntiques a les previsions del mètode lineal de Holt, i també a les del model ETS(A,A,N). Per tant, les previsions puntuals obtingudes del mètode i dels dos models que hi fonamenten són idèntiques (suposant que s’utilitzen els mateixos valors dels paràmetres).

Les previsions puntuals d’ETS són iguals a les medianes de les distribucions de previsió. Per als models amb només components additius, les distribucions de previsió són normals, de manera que les medianes i les mitjanes són iguals. Per als models ETS amb errors multiplicatius o amb estacionalitat multiplicativa, les previsions puntuals no seran iguals a les mitjanes de les distribucions de previsió.

Per obtenir previsions d’un model ETS, utilitzem la funció forecast. El codi R següent mostra els possibles arguments que pren aquesta funció quan s’aplica a un model ETS. Expliquem cadascun dels arguments a continuació.

forecast(object, h=ifelse(object$m>1, 2*object$m, 10),
level=c(80,95), fan=FALSE, simulate=FALSE, bootstrap=FALSE, npaths=5000,
PI=TRUE, lambda=object$lambda, biasadj=NULL, ...)

`object`

L’objecte retornat per la funció ets().

`h`

L’horitzó de previsió: el nombre de períodes que es preveuen.

`level`

El nivell de confiança per als intervals de predicció.

`fan`

Si fan=TRUE, level=seq(50,99,by=1). Això és adequat per a diagrames de ventall.

`simulate`

Si simulate=TRUE, els intervals de predicció es produeixen mitjançant simulació en lloc d’utilitzar fórmules algebraiques. També s’utilitzarà la simulació (fins i tot si simulate=FALSE) on no hi hagi fórmules algebraiques disponibles per al model en particular.

`bootstrap`

Si bootstrap=TRUE i simulate=TRUE, els intervals de predicció simulats utilitzen errors remuestrejats en lloc d’errors distribuïts normalment.

`npaths`

El nombre de camins de mostrals en el càlcul d’intervals de predicció simulats.

`PI`

Si PI=TRUE, es produeixen intervals de predicció; en cas contrari, només es calculen prediccions puntuals.

`lambda`

El paràmetre de transformació de Box-Cox. Això s’ignora si lambda=NULL. En cas contrari, les prediccions es transformen cap enrere mitjançant una transformació inversa de Box-Cox. biasadj

Si lambda no és NULL, les prediccions (i els intervals de predicció) transformades cap enrere s’ajusten al biaix.

fit %>% forecast(h=8) %>%
  autoplot() +
  labs(title = "Previsió amb ETS(A,A,M)", 
       x= "Temps",
       y="Nits de visitants turístics a Austràlia (milions)",
       caption = "Data from fpp2 package"
       )+
  theme_bw()

Intervals de predicció

Un gran avantatge dels models és que també es poden generar intervals de predicció, cosa que no es pot fer amb els mètodes. Els intervals de predicció variaran entre els models amb mètodes additius i multiplicatius.

Per a la majoria dels models ETS, un interval de predicció es pot escriure com

\[ \hat{y}_{T+h|T} \pm k \sigma_h \]

\(k\) depèn de la probabilitat de cobertura: En aquest llibre normalment calculem intervals del 80% i intervals del 95%, tot i que es pot utilitzar qualsevol percentatge. La taula següent dóna el valor de k per a un rang de probabilitats de cobertura assumint errors de previsió distribuïts normalment.

Multiplicadors per a intervals de predicció.

Percentatge	Multiplicador
50	0.67
55	0.76
60	0.84
65	0.93
70	1.04
75	1.15
80	1.28
85	1.44
90	1.64
95	1.96
96	2.05
97	2.17
98	2.33
99	2.58

\(σ_h\) és la variància de la previsió. Els valors de \(k\) es donen a la Taula següent. Per als models ETS, les fórmules per a \(σ_h\) poden ser complicades; els detalls es donen a Rob J Hyndman et al. (2008), Capítol 6. A la taula següent, donem les fórmules per als models ETS additius, que són els més simples.

Expressions de variància de pronòstic per a cada model d’espai d’estats additiu, on \(σ_2\) és la variància residual, \(m\) és el període estacional, \(h_m = \lfloor(h-1)/m\rfloor\) i \(⌊u\)⌋ denota la part entera de \(u\).

\(\tiny\text{ETS(A,N,N) } \sigma_h = \sigma^2\big[1 + \alpha^2(h-1)\big]\)
\(\tiny\text{ETS(A,A,N)} \sigma_h = \sigma^2\Big[1 + (h-1)\big\{\alpha^2 + \alpha\beta h + \frac16\beta^2h(2h-1)\big\}\Big]\)
\(\tiny\text{ETS(A,A\damped,N)} \sigma_h = \sigma^2\biggl[1 + \alpha^2(h-1) + \frac{\beta\phi h}{(1-\phi)^2} \left\{2\alpha(1-\phi) +\beta\phi\right\} - \frac{\beta\phi(1-\phi^h)}{(1-\phi)^2(1-\phi^2)} \left\{ 2\alpha(1-\phi^2)+ \beta\phi(1+2\phi-\phi^h)\right\}\biggr]\)
\(\tiny\text{ETS(A,N,A) } \sigma_h = \sigma^2\Big[1 + \alpha^2(h-1) + \gamma h_m(2\alpha+\gamma)\Big]\)
\(\tiny\text{ETS(A,A,A)} \sigma_h = \sigma^2\Big[1 + (h-1)\big\{\alpha^2 + \alpha\beta h + \frac16\beta^2h(2h-1)\big\} + \gamma h_m \big\{2\alpha+ \gamma + \beta m (h_m+1)\big\} \Big]\)
\(\tiny\text{ETS(A,A\damped,A)} \sigma_h = \sigma^2\biggl[1 + \alpha^2(h-1) +\frac{\beta\phi h}{(1-\phi)^2} \left\{2\alpha(1-\phi) + \beta\phi \right\} - \frac{\beta\phi(1-\phi^h)}{(1-\phi)^2(1-\phi^2)} \left\{ 2\alpha(1-\phi^2)+ \beta\phi(1+2\phi-\phi^h)\right\}\biggl]\)

Per a alguns models ETS, no hi ha fórmules conegudes per als intervals de predicció. En aquests casos, la funció forecast.ets utilitza camins mostrals futurs simulats i calcula els intervals de predicció a partir dels percentils d’aquests camins futurs simulats.

Pronòstic ETS en blocs (bagged)

Una manera útil de millorar la precisió del pronòstic és generar múltiples versions de la sèrie temporal amb lleugeres variacions. A continuació, fer el pronòstic a partir de cadascuna d’aquestes sèries temporals addicionals i fer la mitjana dels prnòstics resultants. Això s’anomena “bagging”, que significa “bootstrap aggregating”. Les sèries temporals addicionals són versions bootstrap de la sèrie original i agreguem els resultats per obtenir una millor previsió.

La mitjana d’aquestes previsions dóna les previsions encapsulades de les dades originals. Tot el procediment es pot gestionar amb la funció baggedETS:

etsfc <- aust %>% 
  ets %>% forecast(h=12)
baggedfc <- aust %>%
  baggedETS %>% forecast(h=12)
autoplot(aust) +
  geom_line(lwd=1)+
  forecast::autolayer(baggedfc$mean, series="BaggedETS",lwd=1) +
  forecast::autolayer(etsfc$mean, series="ETS",lwd=1) +
  guides(colour=guide_legend(title="Pronòstics")) +
  theme_bw()

Per defecte, s’utilitzen 100 sèries bootstrapped i la longitud dels blocs utilitzats per obtenir residus bootstrapped s’estableix en 8 per a dades trimestrals.

En aquest cas, hi ha molt poca diferència. Bergmeir, Hyndman i Benıtez (2016) mostren que, de mitjana, dóna millors previsions que simplement aplicar ets directament. Naturalment, és més lent perquè calen molts més càlculs.

El marc estadístic ETS (Error, Tendència, Estacionalitat) és el que estem tractant. És un enfocament de modelització d’espai d’estats per a la previsió de sèries temporals que descompon les dades en components no observats.↩︎
Recordem que, en els models ETS (Error, Tendència, Estacionalitat), els residuals són les diferències dins la mostra entre els valors observats i els valors ajustats (\(y_t-\hat y_t\)), que actuen com a estimacions del terme d’error no observable. Els errors de previsió d’un pas són les diferències entre els valors reals i les seves previsions d’un pas endavant (\(y_t-\hat y_{t|t-1}\)), sovint utilitzades per a l’avaluació del model.↩︎

Suavització exponencial

Apunts sobre el llibre Forecasting Principles and Practice, Cap. 8

Stefano Sanfilippo

Abril 2026