Econometría II: Media Movil

Departamento de Economía

Carlos A. Yanes G.

2024-08-12

Paquetes con que se trabaja la sesión

R libraries

Los paquetes que vamos a utilizar y que se deben instalar para ser usados en clase son:

install.packages(c("pacman","TSstudio", "urca", "forecast", "devtools", "zoo"))

# Si tiene problemas con TSstudio podemos mirar entonces en Github
devtools::install_github("RamiKrispin/TSstudio")

Ejecución

No olvide que debe desde luego tener presente en su documento de R Markdown cargar estos paquetes.

library(pacman)
p_load(dynlm, fpp2, tidyverse, TSstudio, urca, forecast, zoo)

Preambulo

No paramétricos

  • Las Medias Móviles, los modelos de Atenuación Exponencial (simples y dobles) y los estacional (Holt Winters) son conocidos como modelos no paramétricos.
  • Son fáciles de implementar, sirven para el pronóstico de corto plazo (Máximo 4 periodos).
  • Son útiles siempre y cuando las series no tengan demasiada varianza y cuando el tamaño muestral es \(n \leq 30\).

Media Movil simple

Este método usa el promedio aritmético de los \(n\) valores de datos en la serie para pronosticar los valores futuros. \[\textrm{Promedio Móvil}=\frac{\sum \limits_{i=1}^{n}\left ( X_{i} \right )}{k}\] Donde:

\[\left\{\begin{matrix} k=&\textrm{móvil definido} \\ X_{i}=& \textrm{valores de la serie} \end{matrix}\right.\]

Media Movil

Semana Ventas
1 150 mill
2 230 mill
3 345 mill
4 421 mill
5 434 mill
… …

Media Movil

Semana Ventas Movil K=2
1 150 mill
2 230 mill
3 345 mill 190 mill
4 421 mill 287 mill
5 434 mill 383 mill
… … 427 mill
… … …

La formula detras de los K Moviles

K movil 4

Si tiene un k=4 periodos \(\Rightarrow\)

\[y_{T|t+h}=\frac{y_{t-1}+y_{t}+y_{t+1}+y_{t+2}}{4}\]

K movil 6

Si tiene un k=6 periodos \(\Rightarrow\) \[y_{T|t+h}=\frac{y_{t-3}+y_{t-2}+y_{t-1}+y_{t}+y_{t+1}+y_{t+2}}{6}\]

Media Movil en R

Media Movil en R

Code 
prueba<-c(150,230,345,421,434)
# Media móvil simple
n <- 2
stats::filter(prueba, rep(1 / n, n), sides = 1)
Time Series:
Start = 1 
End = 5 
Frequency = 1 
[1]    NA 190.0 287.5 383.0 427.5
Code 
library(pacman)
p_load(fpp2, tidyverse, zoo)
prueba<-as.data.frame(prueba)
prueba|>
mutate(pr_ma01 = rollmean(prueba, k = 1, fill = NA),
       pr_ma02 = rollmean(prueba, k = 2, fill = NA),
       pr_ma03 = rollmean(prueba, k = 3, fill = NA),
       pr_ma04 = rollmean(prueba, k = 4, fill = NA))|>
  head()
  prueba pr_ma01 pr_ma02  pr_ma03 pr_ma04
1    150     150   190.0       NA      NA
2    230     230   287.5 241.6667   286.5
3    345     345   383.0 332.0000   357.5
4    421     421   427.5 400.0000      NA
5    434     434      NA       NA      NA
Code 
library(pacman)
p_load(fpp2, tidyverse, zoo)
prueba<-as.data.frame(prueba)
prueba|>
mutate(prma01 = ma(prueba, order=1, centre=FALSE),
       prma02 = ma(prueba, order=2, centre=FALSE),
       prma03 = ma(prueba, order=3, centre=FALSE),
       prma04 = ma(prueba, order=4, centre=FALSE))|>
  head()
  prueba prma01 prma02   prma03 prma04
1    150    150  190.0       NA     NA
2    230    230  287.5 241.6667  286.5
3    345    345  383.0 332.0000  357.5
4    421    421  427.5 400.0000     NA
5    434    434     NA       NA     NA

Media Movil Doble

Semanas Yt
1 150
2 230
3 345
4 421
5 434
6 423
7 425
8 430

Media Movil Doble

Semanas Yt Mt
1 150
2 230
3 345 190
4 421 287.5
5 434 383
6 423 427.5
7 425 428.5
8 430 424

Media Movil Doble

Semanas Yt Mt M’t
1 150
2 230
3 345 190
4 421 287.5 238.75
5 434 383 335.25
6 423 427.5 405.25
7 425 428.5 428
8 430 424 426.25

Media Movil Doble

Semanas Yt Mt M’t at
1 150
2 230
3 345 190
4 421 287.5 238.75 336.25
5 434 383 335.25 430.75
6 423 427.5 405.25 449.75
7 425 428.5 428 429
8 430 424 426.25 421.75

Media Movil Doble

Semanas Yt Mt M’t at bt
1 150
2 230
3 345 190
4 421 287.5 238.75 336.25 336.25
5 434 383 335.25 430.75 430.75
6 423 427.5 405.25 449.75 449.75
7 425 428.5 428 429 429
8 430 424 426.25 421.75 421.75

Media Movil Doble

Semanas Yt Mt M’t at bt Yt+p
1 150
2 230
3 345 190
4 421 287.5 238.75 336.25 336.25
5 434 383 335.25 430.75 430.75 672.5
6 423 427.5 405.25 449.75 449.75 861.5
7 425 428.5 428 429 429 899.5
8 430 424 426.25 421.75 421.75 858
9 843.5

Media Movil Doble

Su formula se consolida como

Media Movil Doble

Utiliza ya las medias móviles obtenidas y las promedia nuevamente.

\[\begin{gather*} Y_{t+p} =a_{t} +b_{t} p\\ Donde:\ \\ a_{t} =2M_{t} -M'_{t}\\ b_{t} =\frac{2}{k-1}( M_{t} -M'_{t})\\ p=Periodos\ a\ pronosticar \end{gather*}\]

Media Movil Doble implicita

Esta media viene a consolidarse como

Media movil doble implicita

El promedio del promedio de la media movil. P.e: \[y_t=\frac{1}{2}\left[\frac{y_{t-1}+y_{t}+y_{t+1}}{3}\right]+\frac{1}{2}\left[\frac{y_{t}+y_{t-1}+y_{t-2}}{3}\right]\]

Code 
prueba|>
mutate(pr2ma01 = ma(prueba, order=1, centre=TRUE),
       pr2ma02 = ma(prueba, order=2, centre=TRUE),
       pr2ma03 = ma(prueba, order=3, centre=TRUE),
       pr2ma04 = ma(prueba, order=4, centre=TRUE))|>
  head()
  prueba pr2ma01 pr2ma02  pr2ma03 pr2ma04
1    150     150      NA       NA      NA
2    230     230  238.75 241.6667      NA
3    345     345  335.25 332.0000     322
4    421     421  405.25 400.0000      NA
5    434     434      NA       NA      NA

Media Movil Ponderada

Movil Ponderado

Este método usa también el promedio aritmético de los \(n\) valores de datos en la serie solo que pondera cada una de las observaciones y le da un mayor peso al elemento mas reciente y menor al mas antiguo.

\[\textrm{Móvil ponderado}=\sum \limits_{i=1}^{n} X_{i} \left ( P_{i} \right ) \quad \left\{\begin{matrix} X_{i}= \textrm{valores de la serie}\\ P_{i}= \left \{ 0 \leq P_{i} \leq 1 \right \} \end{matrix}\right.\]

Media Movil Ponderada

Code 
prueba<-c(150,230,345,421,434)
# Media móvil ponderada
stats::filter(prueba, c(0.25, 0.5, 0.25), sides = 1)
Time Series:
Start = 1 
End = 5 
Frequency = 1 
[1]     NA     NA 238.75 335.25 405.25

Suavizado Exponencial

Suavizado Exponencial

Suavizado Exponencial

Este método también le llaman suavización exponencial. Suele ser un caso especial de los no paramétricos y se le da un peso especifico \(\alpha\) para la observación mas reciente.

\[ F_{t+1}=\alpha Y_{t} +(1-\alpha)F_{t}\] Donde:

\(F_{t+1}\) es el pronostico de la serie para el siguiente (t); \(Y_{t}\) es el valor de la serie. \(F_{t}\) es el pronostico para el periodo (t) y \(\alpha\) el valor o constante de suavización que \(\left ( 0 \leq \alpha \leq 1 \right)\)

La intuición del modelo AS

El modelo puede ser descrito como el pronostico para el periodo t=2:

\[\begin{align*} F_{2} &= \alpha Y_{1}+ \left ( 1- \alpha \right)F_{1} \\ &= \alpha Y_{1}+ \left ( 1- \alpha \right)Y_{1} \\ &= Y_{1} \end{align*}\]

Para un t=3 \(\Rightarrow\):

\[F_{3} = \alpha Y_{2}+ \left ( 1- \alpha \right)F_{2}\Rightarrow \alpha Y_{2}+ \left ( 1- \alpha \right)Y_{1}\]

La intuición del modelo AS

Para el siguiente periodo o pronostico del periodo \(t=4\), la ecuación queda: \[\begin{align*} F_{4} &= \alpha Y_{3}+ \left ( 1- \alpha \right)F_{3} \\ &= \alpha Y_{3}+ \left ( 1- \alpha \right) \left [ \alpha Y_{2} + \left ( 1-\alpha \right ) Y_{1} \right ] \\ &= \alpha Y_{3} + \alpha \left ( 1-\alpha \right ) Y_{2}+\left ( 1-\alpha \right )^{2}Y_{1} \end{align*}\]

Los pronósticos tienen en cuenta los pesos de todas las observaciones anteriores: determinar el peso de \(\alpha\) no debe ser tarea difícil.

Determinando un posible \(\alpha\)

El parámetro de suavizado \(\alpha\) se puede mejorar si se establece la siguiente aproximación:

\[\begin{align*} F_{t+1} &= \alpha Y_{t}+ \left ( 1- \alpha \right)F_{t} \\ &= \alpha Y_{t}+ F_{t}- \alpha F_{t} \\ &= \underbrace{F_{t}}_{\textrm{Pronostico en t}} + \alpha \left ( \underbrace{Y_{t} - F_{t}}_{\textrm{Error }} \right ) \end{align*}\]

Holt-Winters

Holt Winters

Holt Winters (no estacional)

Para esta parte se tiene que el método anterior se le puede incorporar una parte lineal que adjuntan dos componentes. El primero hace referencia al \(\alpha\) tradicional de cada observación. Y aparece uno nuevo que será conocido como \(\beta\), le da el mismo tratamiento pero a los efectos proporcionados por la Tendencia.

\[\begin{align*} F_{t}&= \alpha Y_{t}+ (1-\alpha)(F_{t-1}+ T_{t-1}) \\ T_{t}&= \beta (F_{t}-F_{t-1})+ (1-\beta) T_{t-1}\\ \widehat{Y}_{t+h}&= F_{t}+hT_{t} \end{align*}\]

Holt Winters

Holt Winters (estacional)

Muy similar al anterior, solo que se tiene en cuenta el efecto estacional y es configurado de la forma de los estimadores \((\beta) \; \text{y} \; (\alpha)\). Se le denomina gamma (\(\gamma\)) y se encuentra en el intervalo \(\left [ 0 \leq \gamma \leq 1 \right]\).

\[\begin{aligned} F_{t} &= \alpha \frac{Y_{t}}{S_{t-L}}+ (1-\alpha)(F_{t-1}+ T_{t-1}) \\ S_{t} &= \gamma \frac{Y_{t}}{F_{t}}+(1-\gamma)S_{t-L}\\ \widehat{Y}_{t+h}&= F_{t}+hT_{t}+S_{t-L+h} \end{aligned}\]

Holt Winters

Holt Winters (estacional multiplicativo)

En esencia es el mismo método, solo que ahora los componentes se toman todos por uno mismo y no separados como el anterior.

\[\begin{aligned} \widehat{Y}_{t+h}&= F_{t} \times hT_{t} \times S_{t-L+h} \end{aligned}\]

Holt Winters Tendencia ajustada

Holt Winters Damped

Es una versión ampliada de HW y trata de darle un mejor tratamiento a la Tendencia. Aparece un parámetro llamado phi \((\phi)\) quien se encuentra entre \(\left [ 0 \leq \phi \leq 1 \right]\) y es quien ayuda a que los pronósticos de largo plazo sean mas confiables. Amortigua los efectos tendenciales.

\[\begin{aligned} F_{t} &= \alpha Y_{t}+ (1-\alpha)(F_{t-1}+ \phi T_{t-1}) \\ T_{t} &= \beta (F_t-F_{t-1})+(1-\beta)\phi T_{t-1}\\ \widehat{Y}_{t+h}&= F_{t}+\left(\phi+\phi^2+\cdots+\phi^h\right)T_{t} \end{aligned}\]

Caso 3

Datos

Cartera comercial de bancos

Agradecimientos

Lo desarrollado en esta sesión fue construido por los profesores (Hyndman and Athanasopoulos 2018) y (Krispin 2019) en los respectivos paquetes de fpp2 y TSstudio respectivamente

Uso de TSstudio para la serie

Code 
library(pacman)
p_load(readxl, TSstudio, tidyverse, stats, urca, forecast, ggfortify, ggplot2)

The downloaded binary packages are in
    /var/folders/y7/xnq2gn1x2cs79dhd19wnj9wm0000gn/T//RtmpdRxclX/downloaded_packages
Code 
bd <- read_excel("cartera.xlsx")
bd <- bd |> select(Cartera)
cartera <- ts(bd, frequency=12, start=c(2016,1))
head(cartera)
          Jan      Feb      Mar      Apr      May      Jun
2016 367822.1 372017.9 371504.5 373540.2 378922.9 381252.9
Code 
# Para gráficar
ts_plot(cartera, 
        title = "Cartera comercial de bancos",
        Ytitle = "En Miles de Millones de Pesos")
Code 
ts_decompose(cartera)
Code 
ts_seasonal(cartera, type = "normal")
Code 
ts_surface(cartera)

Mapa de calor estacional

Code 
ts_heatmap(cartera, color = "Reds")

Hora de pronosticos

Veamos con las medias moviles

Code 
# Básico
b1<-ts_ma(cartera, n_left = 6, n = NULL) # ma 7 
b1_ma2<-ts_ma(cartera, n_left = 3, n_right = 3, n=NULL) # ma 7

# Listas
 b1_ma7 <- b1$unbalanced_ma_7
 b1_ma72 <- b1_ma2$unbalanced_ma_7
# Objeto base
ma <- cbind(cartera, b1_ma7, b1_ma72)
 p <- ts_plot(ma,
 Xgrid = TRUE,
 Ygrid = TRUE,
 type = "single",
 title = "MA 7 (un lado) vs. MA 7 (dos lados)")

# Grafico 
library(plotly)
 p <- p |> layout(legend = list(x = 0.05, y = 0.95),
 yaxis = list(title = "Miles de Millones"),
 xaxis = list(title = "Años"))
p
Code 
library(forecast)
carma07 = ma(cartera, order=7, centre=FALSE)
carma12 = ma(cartera, order=12, centre=FALSE) # adicional
carma06 = ma(cartera, order=6, centre=FALSE)
carmad12 = ma(cartera, order=12, centre=TRUE)

# Generar pronostico de 12 periodos
MA07 <- forecast(carma07, h = 12)
# Gráfico
plot_forecast(MA07)
Code 
# Generar pronostico de 12 periodos centrado
MAd212 <- forecast(carmad12, h = 12)
# Gráfico
plot_forecast(MAd212)

Números e intervalos

Code 
MAd212
         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
Aug 2023       644750.1 644127.3 645372.9 643797.6 645702.6
Sep 2023       645811.9 644419.0 647204.9 643681.6 647942.3
Oct 2023       646873.7 644542.0 649205.4 643307.6 650439.8
Nov 2023       647935.5 644520.9 651350.1 642713.3 653157.7
Dec 2023       648997.3 644372.0 653622.6 641923.5 656071.1
Jan 2024       650059.1 644107.2 656010.9 640956.5 659161.7
Feb 2024       651120.9 643735.5 658506.2 639825.9 662415.8
Mar 2024       652182.7 643264.0 661101.3 638542.7 665822.6
Apr 2024       653244.4 642698.5 663790.4 637115.8 669373.1
May 2024       654306.2 642043.9 666568.5 635552.7 673059.8
Jun 2024       655368.0 641304.5 669431.6 633859.7 676876.4
Jul 2024       656429.8 640483.7 672376.0 632042.3 680817.3
Code 
MA07
         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
Nov 2023       648056.1 646770.6 649341.6 646090.1 650022.1
Dec 2023       648387.2 645535.7 651238.7 644026.2 652748.2
Jan 2024       648711.7 643976.3 653447.1 641469.6 655953.8
Feb 2024       649029.7 642149.4 655910.0 638507.1 659552.3
Mar 2024       649341.3 640093.9 658588.8 635198.6 663484.1
Apr 2024       649646.8 637838.7 661454.8 631588.0 667705.6
May 2024       649946.1 635406.4 664485.7 627709.6 672182.6
Jun 2024       650239.4 632815.0 667663.7 623591.2 676887.6
Jul 2024       650526.8 630079.8 670973.8 619255.9 681797.8
Aug 2024       650808.5 627213.7 674403.4 614723.3 686893.7
Sep 2024       651084.6 624227.7 677941.5 610010.5 692158.7
Oct 2024       651355.1 621131.6 681578.6 605132.2 697578.0

Exponencial simple

  • El paquete forecast contiene la función de ses. Con esta podemos entonces trabajar el resto de modelos no paramétricos. Lo aprenderemos a usar a continuación.

  • El paquete puede incluso usted cambiar los parámetros como el \((\color{red}{\alpha})\) a su comodidad. Pero recuerde que eso no es recomendable. Debe tener en cuenta que pesa mas, si el pasado o presente de la serie.

Code 
expcar <- ses(cartera, h=12)
# Plot 
plot_forecast(expcar)
Code 
expcar$model
Simple exponential smoothing 

Call:
 ses(y = cartera, h = 12) 

  Smoothing parameters:
    alpha = 0.9999 

  Initial states:
    l = 379123.1825 

  sigma:  4855.33

     AIC     AICc      BIC 
2094.366 2094.624 2102.090
Code 
expcar
         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
Feb 2024       646722.3 640499.9 652944.6 637206.0 656238.5
Mar 2024       646722.3 637923.0 655521.6 633264.9 660179.6
Apr 2024       646722.3 635945.5 657499.0 630240.7 663203.8
May 2024       646722.3 634278.5 659166.0 627691.1 665753.4
Jun 2024       646722.3 632809.8 660634.8 625444.9 667999.6
Jul 2024       646722.3 631481.9 661962.6 623414.2 670030.3
Aug 2024       646722.3 630260.9 663183.7 621546.7 671897.8
Sep 2024       646722.3 629124.3 664320.2 619808.5 673636.0
Oct 2024       646722.3 628056.9 665387.7 618176.0 675268.5
Nov 2024       646722.3 627047.2 666397.3 616631.9 676812.6
Dec 2024       646722.3 626086.9 667357.6 615163.2 678281.3
Jan 2025       646722.3 625169.4 668275.2 613759.9 679684.6
Code 
expcar2 <- ses(cartera, alpha = 0.001, h=12)
expcar2$model
Simple exponential smoothing 

Call:
 ses(y = cartera, h = 12, alpha = 0.001) 

  Smoothing parameters:
    alpha = 0.001 

  Initial states:
    l = 488827.8393 

  sigma:  89841.57

     AIC     AICc      BIC 
2658.452 2658.580 2663.601
Code 
# Plot 
plot_forecast(expcar2)

Holt Winters

Code 
HWa <- hw(cartera,seasonal="additive", h=12)
HWm <- hw(cartera,seasonal="multiplicative", h=12)
HWd <- holt(cartera, damped=TRUE, phi = 0.9, h=12)
HWa$model
Holt-Winters' additive method 

Call:
 hw(y = cartera, h = 12, seasonal = "additive") 

  Smoothing parameters:
    alpha = 0.9999 
    beta  = 0.134 
    gamma = 1e-04 

  Initial states:
    l = 369770.9754 
    b = 1797.9297 
    s = 1090.855 1761.501 97.3342 119.7127 -280.0605 635.3594
           2048.243 984.5148 392.8223 -284.6739 -2310.449 -4255.159

  sigma:  3123.23

     AIC     AICc      BIC 
2021.306 2029.053 2065.077
Code 
plot_forecast(HWa) # Estacional Aditivo
Code 
plot_forecast(HWm) # Estacional Multiplicativo
Code 
plot_forecast(HWd) # Con tendencia amortiguada

Status

Code 
# Modelo Multiplicativo
print(round(accuracy(HWm),2))
                 ME    RMSE     MAE  MPE MAPE MASE ACF1
Training set 923.31 4589.24 3652.03 0.18 0.74  0.1 0.63
Code 
# Modelo Aditivo
print(round(accuracy(HWa),2))
               ME    RMSE     MAE  MPE MAPE MASE ACF1
Training set 21.8 2854.04 2069.87 0.02 0.41 0.06 0.47
Code 
# Modelo damped
print(round(accuracy(HWd),2))
                 ME    RMSE     MAE  MPE MAPE MASE  ACF1
Training set 506.96 3238.59 2329.25 0.11 0.48 0.06 -0.01

Gracias por su atención!!

Referencias

Hyndman, Rob J, and George Athanasopoulos. 2018. Forecasting: Principles and Practice. OTexts.
Krispin, Rami. 2019. Hands-on Time Series Analysis with r: Perform Time Series Analysis and Forecasting Using r. Packt Publishing Ltd.

Universidad del Norte