Suavização Exponencial

ME607 - Séries Temporais

Prof. Carlos Trucíos
ctrucios@unicamp.br

Instituto de Matemática, Estatística e Computação Científica (IMECC),
Universidade Estadual de Campinas (UNICAMP).

Agenda

• Introdução
• Suavização Exponencial Simples
• Método de Holt
• Método de Holt-Winters

Introdução

Introdução

• Uma das limitações dos termos de tendência utilizados nos modelos de regressão é que raramente encontramos séries com tendência determinística (linear, quadrática, cúbica) na vida real.
• Uma alternativa a este problema é utilizar splines. Contudo, seu uso para fins de previsão é bastante limitado, pois não sabemos quantas observações passadas utilizar para estimar/predizer o valor futuro da série.
• Para superar estas limitações, surgiram os modelos de suavização exponencial. A ideia principal destes modelos é permitir com que as últimas observações da série tenham um maior peso nas previsões do que as observações mais antigas.
• Na aula de hoje veremos alguns dos métodos de suavização exponencial mais conhecidos.

Suavização Exponencial Simples

Suavização Exponencial Simples

Utilizado quando a série temporal não tem um padrão claro de tendência e/ou sazonalidade. Por exemplo:

Com este método, a previsão um passo à frente é obtida através de:

\[\hat{y}_{T+1|T} = \alpha y_{T} + (1-\alpha) \hat{y}_{T|T-1} = \hat{y}_{T|T-1} + \alpha (y_T - \hat{y}_{T|T-1}),\] em que \(0 \leq \alpha \leq 1\) é o parâmetro de suavização.

Suavização Exponencial Simples

\[\hat{y}_{T+1|T} = \alpha y_T + (1-\alpha) \underbrace{\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 + \alpha (1-\alpha)y_{T-1} + (1-\alpha)^2 \underbrace{\hat{y}_{T-1|T-2}}_{\alpha y_{T-2} + (1-\alpha) \hat{y}_{T-2|T-3}}\]

\[\hat{y}_{T+1|T} = \alpha y_{T} + \alpha (1-\alpha) y_{T-1} + \alpha (1-\alpha)^2 y_{T-2} + \ldots,\]

\[\hat{y}_{T+1|T} = \alpha [y_{T} + (1-\alpha) y_{T-1} + (1-\alpha)^2 y_{T-2} + \ldots],\]

Note que \((1 + (1 - \alpha) + (1-\alpha)^2+ \cdots) = 1/\alpha\). Ou seja, a previsão \(\hat{y}_{T+1|T}\) é uma média ponderada (com pesos que somam 1) das observações passadas, em que, as observações mais recentes recebem maior peso.

Suavização Exponencial Simples

O peso que recebem as observações mais antigas cai rapidamente.

Suavização Exponencial Simples

\[\hat{y}_{T+1|T} = \alpha y_{T} + \alpha (1-\alpha) y_{T-1} + \alpha (1-\alpha)^2 y_{T-2} + \ldots,\]

• Se \(\alpha = 1\), o método se reduz ao método NAIVE, \(\hat{y}_{T+1|T} = y_T\).
• A medida que \(\alpha \rightarrow 1\), as observações mais recentes recebem um peso maior.

Note que,

• \(\hat{y}_{T+1|T} = \alpha y_T + (1-\alpha) \hat{y}_{T|T-1},\)
• \(\hat{y}_{T|T-1} = \alpha y_{T-1} + (1-\alpha) \hat{y}_{T-1|T-2},\)
• \(\vdots\)
• \(\hat{y}_{3|2} = \alpha y_2 + (1-\alpha) \hat{y}_{2|1},\)
• \(\hat{y}_{2|1} = \alpha y_1 + (1-\alpha) l_0,\)

Assim como nos modelos de regressão, \(l_0\) e \(\alpha\) precisam ser estimados.

Suavização Exponencial Simples

Para obter os valores estimados \(\hat{\alpha}\) e \(\hat{l}_0\), buscamos os valores que minimizem a soma de quadrados dos resíduos, ou seja:

\[SQE = \displaystyle \sum_{t=1}^T (\underbrace{y_t - \hat{y}_{t|t-1}}_{e_t})^2 = \sum_{t=1}^T e_t^2\]

Infelizmente, a solução não tem forma fechada e precisaremos de métodos numéricos para encontrá-la. Na pratica, isto não será problema, utilizaremos a função ETS() do R para ajustar o modelo.

Não é dificil escrever um código para resolver o problema de minimização acima. Explore a função optim() do R (lhe sera muito útil no futuro) e tente fazer você mesmo.

Suavização Exponencial Simples

• Dados
• Gráfico
• Modelo
• Previsão
• Gráfico II

Code

library(fpp3)
algeria_economy <- global_economy |> 
  filter(Country == "Algeria") 
glimpse(algeria_economy)

Rows: 58
Columns: 9
Key: Country [1]
$ Country    <fct> "Algeria", "Algeria", "Algeria", "Algeria", "Algeria", "Alg…
$ Code       <fct> DZA, DZA, DZA, DZA, DZA, DZA, DZA, DZA, DZA, DZA, DZA, DZA,…
$ Year       <dbl> 1960, 1961, 1962, 1963, 1964, 1965, 1966, 1967, 1968, 1969,…
$ GDP        <dbl> 2723648552, 2434776646, 2001468868, 2703014867, 2909351793,…
$ Growth     <dbl> NA, -13.605441, -19.685042, 34.313729, 5.839413, 6.206898, …
$ CPI        <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 2.569025, 2.738580, 2.8…
$ Imports    <dbl> 67.14363, 67.50377, 20.81865, 36.82552, 29.43976, 25.83308,…
$ Exports    <dbl> 39.04317, 46.24456, 19.79387, 24.68468, 25.08406, 22.60394,…
$ Population <dbl> 11124888, 11404859, 11690153, 11985136, 12295970, 12626952,…

Code

algeria_economy %>% autoplot(Exports)

Sem tendência nem sazonalidade claras, canditato perfeito para usar suavização exponencial.

Code

modelo <- algeria_economy %>%
  model(ETS(Exports ~ error("A") )) # "A" para aditivo e "M" para multiplativo
report(modelo)

Series: Exports 
Model: ETS(A,N,N) 
  Smoothing parameters:
    alpha = 0.8399875 

  Initial states:
   l[0]
 39.539

  sigma^2:  35.6301

     AIC     AICc      BIC 
446.7154 447.1599 452.8968 

Code

fore_algeria <- modelo %>% forecast(h = 6)
fore_algeria

# A fable: 6 x 5 [1Y]
# Key:     Country, .model [1]
  Country .model                         Year    Exports .mean
  <fct>   <chr>                         <dbl>     <dist> <dbl>
1 Algeria "ETS(Exports ~ error(\"A\"))"  2018  N(22, 36)  22.4
2 Algeria "ETS(Exports ~ error(\"A\"))"  2019  N(22, 61)  22.4
3 Algeria "ETS(Exports ~ error(\"A\"))"  2020  N(22, 86)  22.4
4 Algeria "ETS(Exports ~ error(\"A\"))"  2021 N(22, 111)  22.4
5 Algeria "ETS(Exports ~ error(\"A\"))"  2022 N(22, 136)  22.4
6 Algeria "ETS(Exports ~ error(\"A\"))"  2023 N(22, 161)  22.4

Code

fore_algeria %>% autoplot(algeria_economy) +
  geom_line(aes(y = .fitted), color = "red",  data = augment(modelo))

Suavização Exponencial Simples

Seja o processo gerados de dados da forma \[y_t = \mu_t + e_t,\] em que \(\mu_t\) é a media (não constante) de \(y_t\).

• Seja \(\hat{\mu}_{T+1|T}\) a previsão um passo à frente do nível.
• Como a média muda ao longo do tempo, podemos obter \(\hat{\mu}_{T+1|T}\) como a solução de \[\displaystyle \sum_{t = 1}^T (y_t - \mu_{T+1})^2\omega_t.\]
• Derivando, temos que \(\hat{\mu}_{T+1|T} = \displaystyle \sum_{t = 1}^T y_t \omega_t.\)
• Vamos impor as restrições de que \(\omega_t > 0\), \(\sum \omega_t = 1\) e que observações recentes recebem maior peso.

Suavização Exponencial Simples

Uma forma de fazer com que as observações recentes recebam peso maior, é fazer \(\omega_t\) descrescer de forma geométrica.

• \(\omega_T = c\), \(\omega_{T-1} = c\phi\), \(\omega_{T-2} = c\phi^2\), \(\cdots\)
• \(1 = c(1 + \phi + \phi^2 + \phi^3+ \cdots) = \dfrac{c}{1-\phi}\)
• \(\hat{\mu}_{T+1|T} = (1-\phi)y_T + (1-\phi)\phi y_{T-1} + (1-\phi)\phi^2 y_{T-2} + \cdots\)

Se fizermos \(1-\phi = \alpha\), temos a mesma Eq. do slide 9.

Suavização Exponencial Simples

veja que o proceso da forma \(y_t = \mu_t + e_t\) é um processo em que nem a tendência nem a sazonalidade são claros.

Code

n <- 1000
e <- rnorm(n, 0, 5)
mu <- runif(n, -10, 10)
y <- mu + e
ts.plot(y)

Suavização Exponencial Simples tem previsão “flat”, i.e. \[\hat{y}_{T+h|T} = \hat{y}_{T+1|T} = \alpha y_T + (1 - \alpha) \hat{y}_{T|T-1}\]

Método de Holt

Método de Holt

Seja o processo gerador de dados da forma\[\begin{align} y_t = & \mu_t + e_t, \\ \mu_t = & \mu_{t-1} + \beta_{t-1}. \end{align}\]

Note que o processo acima representa uma série com tendência linear.

Code

n <- 1000
beta <- runif(n, 0, 0.1)
e <- rnorm(n, 0, 10)
mu <- rep(100, n)
y <- rep(100, n)
for (i in 2:n){
  mu[i] <- mu[i - 1] + beta[i - 1]
  y[i] <- mu[i] + e[i]
}
ts.plot(y)

Método de Holt

• Note que a diferença entre \(\mu_t - \mu_{t-1} = \beta_{t-1}\) (pendente no momento \(t-1\)).
• \(\hat{y}_{t+1|t} = \hat{\mu}_{t+1|t} = \hat{\mu}_{t|t} + \hat{\beta}_{t}\)
• Após observar \(y_{t+1}\), podemos calcular o erro de previsão (\(y_{t+1}-\hat{\mu}_{t+1|t}\)) e corrigir a previsão anterior por uma fração do erro cometido.
• Assim, para \(\theta<1\) \[\begin{align} \hat{\mu}_{t+1|t+1} = & \hat{\mu}_{t+1|t} + (1-\theta) (y_{t+1}-\hat{\mu}_{t+1|t}) \\ = & \hat{\mu}_{t|t} + \hat{\beta}_{t}+ (1-\theta) (y_{t+1}-\hat{\mu}_{t|t} - \hat{\beta}_{t})\\ = & (1 - \theta) y_{t+1} + \theta (\hat{\mu}_{t|t} + \hat{\beta}_t) \end{align}\]
• De forma semelhante, para \(\gamma <1\), \[\hat{\beta}_{t+1} = \hat{\beta}_t + (1-\gamma) (\hat{\mu}_{t+1|t+1} - \hat{\mu}_{t|t} - \hat{\beta}_t) = \gamma \hat{\beta}_t + (1-\gamma) (\hat{\mu}_{t+1|t+1} - \hat{\mu}_{t|t})\]

Método de Holt

• Utilizado quando a série apresenta tendência.
• A previsão envolve duas equações de suavização: uma para o nível e outra para a tendência:
  • \(\text{Previsão: } \quad \hat{y}_{T+h|T} = \hat{\mu}_{T|T} + h \hat{\beta}_T,\)
  • \(\text{Equação do Nível: } \hat{\mu}_{T|T} = \theta^{\ast} y_T + (1 - \theta^{\ast}) (\hat{\mu}_{T-1|T-1} + \hat{\beta}_{T-1}),\)
  • \(\text{Equação de Tendência: } \quad \hat{\beta}_T= \gamma^{\ast} \hat{\beta}_{T-1} + (1-\gamma^{\ast}) (\hat{\mu}_{T|T} - \hat{\mu}_{T-1|T-1}).\)

Em que \(0 \leq \theta^{\ast} \leq 1\) e \(0 \leq \gamma^{\ast} \leq 1\). Note que, no instante \(t=1\) precisaremos de \(\hat{\mu}_0\) e \(\hat{\beta}_0\)

Para estimar os parâmetros, minimizamos a soma de quadrados dos resíduos.

Método de Holt

• Dados
• Gráfico
• Modelo
• Previsão
• Gráfico II

Code

australia_economy <- global_economy |> 
  filter(Country == "Australia") |> 
  mutate(Pop = Population / 1e6) # População em milhoes 1e6 = 10^6
glimpse(australia_economy)

Rows: 58
Columns: 10
Key: Country [1]
$ Country    <fct> "Australia", "Australia", "Australia", "Australia", "Austra…
$ Code       <fct> AUS, AUS, AUS, AUS, AUS, AUS, AUS, AUS, AUS, AUS, AUS, AUS,…
$ Year       <dbl> 1960, 1961, 1962, 1963, 1964, 1965, 1966, 1967, 1968, 1969,…
$ GDP        <dbl> 18573188487, 19648336880, 19888005376, 21501847911, 2375853…
$ Growth     <dbl> NA, 2.4856050, 1.2964777, 6.2142784, 6.9787237, 5.9834500, …
$ CPI        <dbl> 7.960458, 8.142560, 8.116545, 8.168574, 8.402706, 8.688866,…
$ Imports    <dbl> 14.06175, 15.02508, 12.63093, 13.83405, 13.76450, 15.26734,…
$ Exports    <dbl> 12.99445, 12.40310, 13.94301, 13.00589, 14.93825, 13.22018,…
$ Population <dbl> 10276477, 10483000, 10742000, 10950000, 11167000, 11388000,…
$ Pop        <dbl> 10.27648, 10.48300, 10.74200, 10.95000, 11.16700, 11.38800,…

Code

australia_economy %>% autoplot(Pop)

Tendência clara.

Code

modelo <- australia_economy %>%
  model(ETS(Pop ~ error("A") + trend("A")))
report(modelo)

Series: Pop 
Model: ETS(A,A,N) 
  Smoothing parameters:
    alpha = 0.9999 
    beta  = 0.3266366 

  Initial states:
     l[0]      b[0]
 10.05414 0.2224818

  sigma^2:  0.0041

      AIC      AICc       BIC 
-76.98569 -75.83184 -66.68347 

Code

fore_aus <- modelo %>% forecast(h = 6)
fore_aus

# A fable: 6 x 5 [1Y]
# Key:     Country, .model [1]
  Country   .model                                    Year           Pop .mean
  <fct>     <chr>                                    <dbl>        <dist> <dbl>
1 Australia "ETS(Pop ~ error(\"A\") + trend(\"A\"))"  2018 N(25, 0.0041)  25.0
2 Australia "ETS(Pop ~ error(\"A\") + trend(\"A\"))"  2019  N(25, 0.011)  25.3
3 Australia "ETS(Pop ~ error(\"A\") + trend(\"A\"))"  2020  N(26, 0.023)  25.7
4 Australia "ETS(Pop ~ error(\"A\") + trend(\"A\"))"  2021  N(26, 0.039)  26.1
5 Australia "ETS(Pop ~ error(\"A\") + trend(\"A\"))"  2022  N(26, 0.061)  26.4
6 Australia "ETS(Pop ~ error(\"A\") + trend(\"A\"))"  2023   N(27, 0.09)  26.8

Code

fore_aus %>% autoplot(australia_economy) +
  geom_line(aes(y = .fitted), color = "red",  data = augment(modelo)) 

Método de Holt com amortização

• O método de Holt (1957) assume que, para fazer a previsão, a tendência é constante nos \(h\) passos à frente.
• Em 1985, é proposto um método que amortece a tendência. Com este método as equações são definidas como:

\[\begin{align} \hat{y}_{t+h|t} = & \hat{\mu}_{t|t} + (\phi + \phi^2 + \cdots + \phi^h) \hat{\beta}_t, \\ \hat{\mu}_{t|t} = & \alpha y_t + (1-\alpha) (\hat{\mu}_{t-1|t-1} + \phi \hat{\beta}_{t-1}), \\ \hat{\beta}_t = & \gamma (\hat{\mu}_{t|t} - \hat{\mu}_{t-1|t-1}) + (1 - \gamma) \phi \hat{\beta}_{t-1}. \end{align}\]

Note que se \(\phi = 1\) temos o método de Holt.

Método de Holt com tendência não constante

• Modelo
• Previsão
• Gráfico II

Code

modelo <- australia_economy |> 
  model(holt = ETS(Pop ~ error("A") + trend("A")),
        amortecimento = ETS(Pop ~ error("A") + trend("Ad")))
modelo |> select(amortecimento) |> report()

Series: Pop 
Model: ETS(A,Ad,N) 
  Smoothing parameters:
    alpha = 0.9986305 
    beta  = 0.4271838 
    phi   = 0.98 

  Initial states:
     l[0]      b[0]
 10.03648 0.2478533

  sigma^2:  0.0045

      AIC      AICc       BIC 
-71.01630 -69.36924 -58.65364 

Code

fore_aus <- modelo %>% forecast(h = 6)
fore_aus

# A fable: 12 x 5 [1Y]
# Key:     Country, .model [2]
   Country   .model         Year           Pop .mean
   <fct>     <chr>         <dbl>        <dist> <dbl>
 1 Australia holt           2018 N(25, 0.0041)  25.0
 2 Australia holt           2019  N(25, 0.011)  25.3
 3 Australia holt           2020  N(26, 0.023)  25.7
 4 Australia holt           2021  N(26, 0.039)  26.1
 5 Australia holt           2022  N(26, 0.061)  26.4
 6 Australia holt           2023   N(27, 0.09)  26.8
 7 Australia amortecimento  2018 N(25, 0.0045)  25.0
 8 Australia amortecimento  2019  N(25, 0.014)  25.3
 9 Australia amortecimento  2020  N(26, 0.029)  25.6
10 Australia amortecimento  2021  N(26, 0.051)  26.0
11 Australia amortecimento  2022  N(26, 0.082)  26.3
12 Australia amortecimento  2023   N(27, 0.12)  26.6

Code

fore_aus |> autoplot(australia_economy, level = NULL) +
  geom_line(aes(y = .fitted, color = .model),  data = augment(modelo)) 

Método de Holt com tendência não constante

Note que, quando \(h \rightarrow \infty\), \[y_{T+h|T} = \hat{\mu}_{T|T} + (\phi + \phi^2 + \cdots + \phi^h) \hat{\beta}_T \rightarrow l_t + \dfrac{\phi}{1-\phi} \hat{\beta}_T\]

Método de Holt-Winters

Método de Holt-Winters

Muitas vezes, a série temporal além de possuir tendência, possui também sazonalidade. Para lidar com estas características, temos o modelo de Holt-Winters, que é uma extensão do método de Holt. Seja \(m\) o período sazonal, então:

Holt-Winters Aditivo

\[\hat{y}_{t+h|t} = \hat{\mu}_{t|t} + h\hat{\beta}_t + s_{t+h-m(k+1)},\]

\[l_t = \alpha (y_t - s_{t-m}) + (1-\alpha) (\hat{\mu}_{t-1|t-1}+ \hat{\beta}_{t-1}),\]

\[b_t = \delta (\hat{\mu}_{t|t} - \hat{\mu}_{t-1|t-1}) + (1-\delta) \hat{\beta}_{t-1},\]

\[s_t = \gamma (y_t - \hat{\mu}_{t-1|t-1} - \hat{\beta}_{t-1}) + (1-\gamma) s_{t-m}.\] Em que \(k\) é a parte inteira de \((h-1)/m\).

Holt-Winters Multiplicativo

\[\hat{y}_{t+h|t} = (\hat{\mu}_{t|t} + h\hat{\beta}_{t})s_{t+h-m(k+1)},\]

\[l_t = \alpha \dfrac{y_t}{s_{t-m}} + (1-\alpha) (\hat{\mu}_{t-1|t-1} + \hat{\beta}_{t-1}),\]

\[b_t = \delta (l_t - \hat{\mu}_{t-1|t-1}) + (1-\delta) \hat{\beta}_{t-1},\]

\[s_t = \gamma \dfrac{y_t}{\hat{\mu}_{t-1|t-1} + \hat{\beta}_{t-1}} + (1-\gamma) s_{t-m}\]

Método de Holt-Winters

Se a variação sazonal é constante utilizamos o modelo aditivo, caso contrário o multiplicativo.

• Dados
• Gráfico
• Modelo
• Previsão
• Gráfico II

Code

aus_holidays <- tourism %>% 
  filter(Purpose == "Holiday") %>%
  summarise(Trips = sum(Trips)/1000) # Agrupamos por Quarter e os valores unidades de mil. 
glimpse(aus_holidays)

Rows: 80
Columns: 2
$ Quarter <qtr> 1998 Q1, 1998 Q2, 1998 Q3, 1998 Q4, 1999 Q1, 1999 Q2, 1999 Q3,…
$ Trips   <dbl> 11.806038, 9.275662, 8.642489, 9.299524, 11.172027, 9.607613, …

Code

aus_holidays %>% autoplot(Trips)

Code

modelo <- aus_holidays |> 
  model(HWA = ETS(Trips ~ error("A") + trend("A") + season("A")),
        HWM = ETS(Trips ~ error("M") + trend("A") + season("M")))
report(modelo)

# A tibble: 2 × 9
  .model  sigma2 log_lik   AIC  AICc   BIC   MSE  AMSE    MAE
  <chr>    <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
1 HWA    0.193     -105.  229.  231.  250. 0.174 0.184 0.321 
2 HWM    0.00212   -104.  227.  229.  248. 0.170 0.183 0.0328

Code

fore_holidays <- modelo %>% forecast(h = 6)
fore_holidays

# A fable: 12 x 4 [1Q]
# Key:     .model [2]
   .model Quarter       Trips .mean
   <chr>    <qtr>      <dist> <dbl>
 1 HWA    2018 Q1 N(13, 0.19)  12.9
 2 HWA    2018 Q2 N(11, 0.21)  11.2
 3 HWA    2018 Q3 N(11, 0.23)  11.0
 4 HWA    2018 Q4 N(11, 0.26)  11.2
 5 HWA    2019 Q1  N(13, 0.3)  13.4
 6 HWA    2019 Q2 N(12, 0.34)  11.7
 7 HWM    2018 Q1 N(13, 0.37)  13.3
 8 HWM    2018 Q2 N(11, 0.28)  11.2
 9 HWM    2018 Q3 N(11, 0.28)  10.8
10 HWM    2018 Q4 N(11, 0.32)  11.1
11 HWM    2019 Q1 N(14, 0.55)  13.8
12 HWM    2019 Q2 N(12, 0.43)  11.7

Code

fore_holidays |> autoplot(aus_holidays) +
  geom_line(aes(y = .fitted, color = .model), data = augment(modelo)) 

Método de Holt-Winters

Se quisermos incluir o amorcetimento na componente de tendência, basta trocar trend("A") por trend("Ad") nos códigos do slide anterior.

Comentários finais

• Além dos métodos mencionados, existem vários outros métodos de suavização.
• Esses métodos podem ser obtidos como as possíveis configurações que podem ser obtidas com:
  • error(): A (aditivo) ou M (multiplicativo)
  • season(): A, M ou N (nenhum).
  • trend(): N, A, M, Ad (aditivo amortecido), Md (multiplicativo amortecido).
• A escolha não sempre é facil. Embora não sejam infalíveis, existem métodos automáticos (teste fazendo, por exemplo, modelo <- aus_holidays |> model(auto = ETS(Trips)).
• Maiores detalhes a respeito do procedimento automático podem ser enconrados em Hydman et al. (2002).
• Para um maior aprofundamento nos metodos de suavização exponencial, ver Hyndman et al. (2008). Link com código e erratas, aqui.

Referências

• Hyndman, R. J., Koehler, A. B., Snyder, R. D., & Grose, S. (2002). A state space framework for automatic forecasting using exponential smoothing methods. International Journal of forecasting, 18(3), 439-454.
• Hyndman, R., Koehler, A. B., Ord, J. K., & Snyder, R. D. (2008). Forecasting with exponential smoothing: the state space approach. Springer Science & Business Media.
• Hyndman, R.J., & Athanasopoulos, G. (2021). Forecasting: principles and practice, 3rd edition, OTexts: Melbourne, Australia. OTexts.com/fpp3.. Chapter 8.
• Peña, D. (2005). Análisis de series temporales. Alianza. Capítulo 2.