library("forecast")
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library("fpp2")
## Loading required package: ggplot2
## Loading required package: fma
## Loading required package: expsmooth

Métodos de pronóstico naive

Un pronóstico naive es simplemente el valor observado mpas recientemente.

\[ F_{t+k} = y_{t} \] Incluso si no es un método excelente para realizar pronósticos, es un buen benchmark para empezar a hacer predicciones más acertadas basadas en otro modelos.

fc_goog <- naive(goog, 10)
summary(fc_goog)
## 
## Forecast method: Naive method
## 
## Model Information:
## Call: naive(y = goog, h = 10) 
## 
## Residual sd: 8.7285 
## 
## Error measures:
##                     ME     RMSE      MAE        MPE      MAPE MASE       ACF1
## Training set 0.4212612 8.734286 5.829407 0.06253998 0.9741428    1 0.03871446
## 
## Forecasts:
##      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 1001         813.67 802.4765 824.8634 796.5511 830.7889
## 1002         813.67 797.8401 829.4999 789.4602 837.8797
## 1003         813.67 794.2824 833.0576 784.0192 843.3208
## 1004         813.67 791.2831 836.0569 779.4322 847.9078
## 1005         813.67 788.6407 838.6993 775.3910 851.9490
## 1006         813.67 786.2518 841.0882 771.7374 855.6025
## 1007         813.67 784.0549 843.2850 768.3777 858.9623
## 1008         813.67 782.0102 845.3298 765.2505 862.0895
## 1009         813.67 780.0897 847.2503 762.3133 865.0266
## 1010         813.67 778.2732 849.0667 759.5353 867.8047

La función naive también arroja intervalos de confianza para la predicción. Por ejemplo, con un nivel de confianza del 95% la predicción estará en:

\[ \hat{y_{t}} \pm 1.96\hat{\sigma} \] Donde \(\hat{\sigma}\) es la desviación estandar de la distribución del pronóstico

Entre más lejos el pronóstico más incertidumbre asociada al mismo.En el gráfico se puede ver como el intervalo se hace más grande a medida que se hacen pronósticos más alejados en el tiempo.

fc_goog <- naive(goog, 25)

autoplot(fc_goog)

### DAtos estacionales

Una modificación del pronóstico naive es usar los mismos valores de una periodicidad pasada:

$$ F_{t+k} = y_{t-M+k}

$$

fc_beer <- snaive(ausbeer, 16)
summary(fc_beer)
## 
## Forecast method: Seasonal naive method
## 
## Model Information:
## Call: snaive(y = ausbeer, h = 16) 
## 
## Residual sd: 19.1207 
## 
## Error measures:
##                    ME     RMSE      MAE      MPE    MAPE MASE       ACF1
## Training set 3.098131 19.32591 15.50935 0.838741 3.69567    1 0.01093868
## 
## Forecasts:
##         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2010 Q3            419 394.2329 443.7671 381.1219 456.8781
## 2010 Q4            488 463.2329 512.7671 450.1219 525.8781
## 2011 Q1            414 389.2329 438.7671 376.1219 451.8781
## 2011 Q2            374 349.2329 398.7671 336.1219 411.8781
## 2011 Q3            419 383.9740 454.0260 365.4323 472.5677
## 2011 Q4            488 452.9740 523.0260 434.4323 541.5677
## 2012 Q1            414 378.9740 449.0260 360.4323 467.5677
## 2012 Q2            374 338.9740 409.0260 320.4323 427.5677
## 2012 Q3            419 376.1020 461.8980 353.3932 484.6068
## 2012 Q4            488 445.1020 530.8980 422.3932 553.6068
## 2013 Q1            414 371.1020 456.8980 348.3932 479.6068
## 2013 Q2            374 331.1020 416.8980 308.3932 439.6068
## 2013 Q3            419 369.4657 468.5343 343.2438 494.7562
## 2013 Q4            488 438.4657 537.5343 412.2438 563.7562
## 2014 Q1            414 364.4657 463.5343 338.2438 489.7562
## 2014 Q2            374 324.4657 423.5343 298.2438 449.7562
autoplot(fc_beer)

Valores pronósticados y residuos

Es importante revisar que los errores se asemejan a ruido blanco.

Los invervalos de confianza de lo pronósticos son calculados asumiendo que los errores:

La función checkresiduals genera un histograma y calcula una prueba Ljung-Box de los errores del modelo

checkresiduals(fc_goog)

## 
##  Ljung-Box test
## 
## data:  Residuals from Naive method
## Q* = 13.123, df = 10, p-value = 0.2169
## 
## Model df: 0.   Total lags used: 10

Train -Test en series de tiempo

A diferencia de datos de corte transversal donde las partición puese ser aleatoria, en series de tiempo se toma como train la parte mas vieja de los datos y como test la parte más reciente, en el mismo orden.

train <- subset(gold, end = 1000)
test <- subset(gold, start = 1001, end = length(gold))

Para no hacer particiones muy estrictas (Semanas meses años) que podrían hacer perder datos estacionales, las particiones se pueden hacer por ventanas

train2 <- window(ausbeer, end = c(1995, 4))

Esta funciones es bastante util cuando queremos una ventana de tiempo bien definida

train3 <- window(ausbeer, start = c(1970, 1), end = c(1995, 4))

Precisión de prónosticos

El error de pronóstico está definido como:

\[ e_{t} = y_{t} - \hat{y_{t}} \]

Y estas son algunas medidas de error de pronóstico:

Mean absolute error \[ MAE = \sum_{t=1}^T|e_{t}| \]

Mean Squared error \[ MSE = \sum_{t=1}^Te_{t}^2 \]

Mean Absolute percent error

\[ MAE = \sum_{t=1}^T|\frac{e_{t}}{y_{t}}| \] Mean absolute scaled error

\[ MASE = MAE/Q \]

Donde \(Q\) es una constante de reescalamiento.

Revisar ACF1 y U de theil

Ejemplo de evaluación de pronósticos con métodos no estacionales

train <- subset(gold, end= 1000)

# Pronóstico basado en el ultima observación
naive_fc <- naive(train, h = 108)

# Pronóstico basado en la media de todas las observaciones de entrenamiento
mean_fc <- meanf(train, h = 108)

# Medidas de precisión del pronóstico naive
accuracy(naive_fc, gold)
##                      ME      RMSE      MAE        MPE      MAPE     MASE
## Training set  0.1079897  6.358087  3.20366  0.0201449 0.8050646 1.014334
## Test set     -6.5383495 15.842361 13.63835 -1.7462269 3.4287888 4.318139
##                    ACF1 Theil's U
## Training set -0.3086638        NA
## Test set      0.9793153  5.335899
# Medidas de precisión sobre la media
accuracy(mean_fc, gold)
##                         ME     RMSE      MAE       MPE      MAPE      MASE
## Training set -4.239671e-15 59.17809 53.63397 -2.390227 14.230224 16.981449
## Test set      1.319363e+01 19.55255 15.66875  3.138577  3.783133  4.960998
##                   ACF1 Theil's U
## Training set 0.9907254        NA
## Test set     0.9793153  6.123788
train <- subset(gold, end = 1000)
test <- subset(gold, start = 1001, end = length(gold))
naive_fc <- naive(train, h = 108)
accuracy(naive_fc, test)
##                      ME      RMSE      MAE        MPE      MAPE     MASE
## Training set  0.1079897  6.358087  3.20366  0.0201449 0.8050646 1.014334
## Test set     -6.5383495 15.842361 13.63835 -1.7462269 3.4287888 4.318139
##                    ACF1 Theil's U
## Training set -0.3086638        NA
## Test set      0.9793153  5.335899

Ejemplo de evaluación de pronóstico con métodos estacionales

train2 <- window(ausbeer, end = c(1995, 4))
test <- window(ausbeer, start = c(1996, 1), end = c(2004, 4))


snaive_fc <- snaive(train2, h = length(test))
mean_fc <- meanf(train2, h = length(test))


accuracy(snaive_fc, test)
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set  4.730769 20.60589 16.51282  1.284697 3.957622 1.0000000
## Test set     -8.527778 18.06854 13.08333 -2.011727 3.038430 0.7923137
##                     ACF1 Theil's U
## Training set  0.01783674        NA
## Test set     -0.39101873 0.2890478
accuracy(mean_fc, test)
##                         ME     RMSE      MAE       MPE      MAPE     MASE
## Training set -2.279427e-14 96.37835 78.64437 -6.646984 21.996740 4.762625
## Test set      2.393472e+01 49.38745 34.07500  4.640789  7.243654 2.063548
##                     ACF1 Theil's U
## Training set  0.72972747        NA
## Test set     -0.09479634 0.8164006

Cross-Validation en series de tiempo

el Procedimiento en general es similar al de datos de corte transversal pero el conjunto de entrenamiento siempre tiene que estar antes del de testeo en CV para series de tiempo. Esto implica:

Asumiento que vamos a hacer un pronóstico de solo un periodo hacia adelante, el CV sería:

  1. Seleccionar las observaciones hasta un periodo \(k\)
  2. Seleccionar la observación \(k+1\) como dato de test
  3. Se descartan las observaciones \(k+2, k+3,...,k+n\)
  4. Calculamos los errores de pronóstico para \(k+1\)
  5. Repetimos los pasos 1-4 para \(i=1,2,...,T-k\) donde \(T\) es el número total de observaciones
  6. Calcular las medidas de error de pronóstico basado en los errores obtenidos

Los pronósticos empiezan desde \(k+i-1\), que es el origen

Sin embargo la evaluación de solo un periodo adelante podría no ser util.

Supongamos que queremos un modelo que produzca un pronóstico \(h\) periodos hacia adelante en vez de 1. Lo que habría que hacer es ajustar el algoritmo anterior de manera que seleccionemos la observación en el momento \(k+h+i-1\) en el set de testeo y volver a aplicar cross validation

Ejemplo de Cross-validation con h=1

errors <- tsCV(goog, forecastfunction = naive, h = 1)
mean(errors^2, na.rm = TRUE)
## [1] 76.28775

Podemos cambiar los horizontes de comparación para ver la exactitud de los pronósticos cuando se aleja más el dato a pronosticar

MSE <- vector("numeric", 10)
for(h in 1:10) {
  errors <- tsCV(goog, forecastfunction = naive, h = h)
  MSE[h] <- mean(errors^2, na.rm = TRUE)
}

MSE
##  [1]  76.28775 117.60458 158.87772 197.78048 234.10338 268.50152 302.10622
##  [8] 335.04660 366.69461 397.21987