Caso 3 - Series de tiempo
Maria Camila Corredor Bernal, Nicol Daniela Peña Bermudez
31/5/2021
INDIA HOLIDAY BUREAU: PLANNING FOR THE TOURIST SEASON
Caso 3 - Procesos Estocasticos
Preguntas
1. Realizar un detallado análisis descriptivo de la serie de tiempo incluyendo los gráficos de la serie, los correlogramas, diagramas de dispersión y cualquier otro que considere de interés.
El turismo ha tenido una posición indispensable dentro de la economía india, formando parte en un 6,6% del producto interno bruto según el Consejo Mundial de Viajes y Turismo. Por ello, muchas empresas dedicadas a ofrecer este tipo de servicios se encuentran en la búsqueda de mejores métodos para poder ofrecer a sus clientes un servicio único, confiable y excelente; este es el caso de India Holiday Bureau, una agencia de viajes de la India que en 2013 reportó indicadores muy bajos de servicio debido a la poca planeación que se tenía en la empresa en cuanto a la demanda creciente de turistas.
A partir de lo anterior, a continuación se presenta la serie de tiempo que representa las llegadas mensuales de turistas extranjeros a la India en el periodo transcurrido entre 2001 y 2012, en donde podemos observar una clara tendencia ascendente a lo largo de estos años, además de un efecto en la amplitud de los picos, pues se puede apreciar como estos se van haciendo más longevos y pronunciados al pasar los años.
En consecuencia a su posicionamiento turístico dentro de esta industria, la India fue clasificada como el tercer destino turístico de más rápido crecimiento en el mundo. Además, según cifras oficiales, se espera que la industria turística mundial crezca a una tasa media anual del 8% durante la década de 2015 a 2024, situación a la que debe anteponerse la agencia India Holiday Bureau para lograr brindarle a sus clientes el servicio que ellos están esperando y poder mejorar sus indicadores de desempeño. De igual manera, podemos observar un patrón repetitivo año a año, por lo que se puede afirmar que la serie no es estacionaria, pues además de la marcada tendencia, también se presenta un efecto estacional.
Cabe aclarar que los datos presentados hacen referencia a la ventana de entrenamiento, tomada desde enero 1 de 2001 hasta diciembre 31 de 2012. En cuanto a la ventana de prueba, se tomó un periodo de 2 años, con inicio en enero 1 de 2013 y fin el 31 de diciembre de 2014.
# Grafico de la serie (Train)
autoplot(train, col= "cadetblue")+theme_minimal()+ ylab("Numero de turistas")+ xlab("")+ ggtitle("Llegadas mensuales de turistas extranjeros a la india")+theme_bw()
En el siguiente gráfico se presentan las llegadas mes a mes de turistas extranjeros a la India, y en las diferentes curvas de nivel se presentan los años evaluados. Como se concluyó anteriormente, la serie tiene un evidente comportamiento estacional, en donde el periodo con más llegadas reportadas hace referencia a la época vacacional de fin e inicio de año, igualmente se presenta un aumento en el mes de julio.
ggseasonplot(train, year.labels = TRUE, year.labels.left = TRUE)+ ggtitle("Gráfico estacional de las llegadas mensuales de turistas extranjeros")+ xlab("Mes")+ ylab("Numero de turistas")+ theme_minimal()
Con el propósito de conocer el método de pronóstico que ayude en mayor medida a la empresa IHB, a continuación se presentan los correlogramas de la serie, no obstante, cabe aclarar que sobre ellos no es posible elegir tentativamente un modelo pues la serie que se está trabajando no es estacionaria.
Acf(train, main="Gráfico ACF", col="maroon4")
Pacf(train, main="Gráfico PACF", col="deeppink")
Teniendo en cuenta los 2 gráficos anteriores correspondientes a los correlogramas, se podría pensar en realizar un pronóstico ARMA, pues tanto el ACF como el PACF tienen un comportamiento sinusoidal, sin embargo, para confirmar esta información será necesario realizar correcciones sobre la serie para convertirla en estacionaria y poder tener un análisis más acertado.
2. De ser necesario, realizar modificaciones a la serie: tomar logaritmo, diferenciación y/o diferenciación estacional. Explicar cuál es la razón para realizarlas
Como se presentó anteriormente, la serie se caracteriza por tener tanto tendencia como estacionalidad, por lo tanto, a continuación, se realizarán los ajustes necesarios para convertir la serie a una estacionaria.
Al tener las 2 características ya mencionadas, se procede en primera medida a realizar una diferenciación estacional con el propósito de eliminar el marcado efecto estacional. Esta modificación se realiza restando al valor real del instante t, el valor real de 12 periodos atrás, es decir, 12 meses, como se muestra en la siguiente ecuación:
\(x_t:y_t-y_{t-k}\)
Serie_dif <- diff(train, lag=12, differences = 1)
Serie_dif## Jan Feb Mar Apr May Jun Jul Aug Sep
## 2002 -55600 -21173 -32126 -25549 -6527 -42150 -46201 -33923 763
## 2003 46065 21559 1634 1152 -3063 41758 47128 42346 28250
## 2004 63130 69005 74712 62943 43994 46798 47097 48361 35434
## 2005 48632 38147 58909 24532 39892 23848 35414 20555 30411
## 2006 73512 69246 38915 60792 29614 31400 29462 30531 40707
## 2007 76142 62602 81485 41342 22009 31994 62534 54059 4001
## 2008 -23850 109801 7271 10551 27344 31175 32067 24891 39801
## 2009 -30473 -121706 -37703 -13557 822 10814 967 -13630 -10986
## 2010 87411 62365 70090 24412 26904 32289 33815 52466 39114
## 2011 53994 75567 23461 74555 51352 20822 8829 6317 47657
## Oct Nov Dec
## 2002 31662 35976 41930
## 2003 47302 44922 22797
## 2004 46878 94655 98256
## 2005 40310 38599 61884
## 2006 43642 18576 62160
## 2007 53165 90015 54989
## 2008 5449 -745 -62656
## 2009 8836 9841 81871
## 2010 48244 66654 64229
## 2011 52548 61589 56839autoplot(Serie_dif)+theme_minimal()+ylab("")+ xlab("")+ ggtitle("Gráfico para diferenciación estacional")
Al aplicar la diferenciación estacional podemos ver que este efecto se logró corregir en gran proporción y asimismo, se puede apreciar que el marcado efecto de tendencia disminuyó, sin embargo la serie continúa teniendo un efecto de tendencia a trozos, por lo tanto se procede a tratar esta característica con una diferenciación a un paso, es decir, a un mes.
Serie_dif2 <- diff(Serie_dif, lag=1, differences = 1)
autoplot(Serie_dif2)+theme_minimal()+ylab("")+ xlab("")+ ggtitle("Gráfico serie diferenciada a un paso")
Al graficar la serie diferenciada estacionalmente y a un paso, podemos ver finalmente que tanto el efecto estacional como el efecto de tendencia fueron controlados y además se logró ubicar su centro en el valor de 0, por tanto, la serie ya es estacionaria.
ggseasonplot(Serie_dif2, year.labels = TRUE, year.labels.left = TRUE)+ ggtitle("Gráfico estacional de la serie diferenciada (llegadas mensuales de turistas extranjeros)")+ xlab("Mes")+ ylab("Numero de turistas")+ theme_minimal()
Ahora bien, como la serie de tiempo fue modificada y tiene comportamiento estacionario, a continuación se presentan los correlogramas, sobre los cuales se analizará una propuesta de modelo de pronóstico.
Acf(Serie_dif2, main="Gráfico ACF", col="maroon4")
Pacf(Serie_dif2, main="Gráfico PACF", col="deeppink")
Teniendo en cuenta los resultados obtenidos en los graficos anteriores, se logra evidenciar que las graficas PACF Y ACF, presentan un comportamiento sinusoidal y un leve comportamiento exponencial, así mismo en el ACF podemos ver que las correlaciones 1,11,12 y 22 son significativas, esto puede deberse al comportamiento estacional inicial de la serie de tiempo y, ademas el PACF tiene como correlaciones significativas los lag’s 1,5,6,12 y 24, en donde estos ultimos tambien pueden verse afectados por el comportamiento estacional. Con base en lo anterior, tentativamente se puede pensar que estamos frente a un modelo SARIMA.
3. A partir de los dos puntos anteriores definir el posible modelo que se ajuste al comportamiento de la serie: ingenuo, AR, MA, ARMA, ARIMA, SARIMA.
Teniendo en cuenta el analisis propuesto en el punto anterior, se propone analizar un modelo SARIMA ya que para volver la serie estacionaria se tuvo que hacer una corrección estacional por medio de una diferención a 12 meses y una corrección para la tendencia con una diferenciación a un paso. Ademas, teniendo en cuenta que en ambos correlogramas se tiene un mismo comportamiento sinusoidal con leves rasgos de decaimiento exponencial y que aparentemente solo las primeras correlaciones en ambos graficos son significativas y estas se vuelven a repetir en periodos alejados, lo que se puede deber al efecto estacional, se propone, un valor de P=1 y Q=1, lo que finalmente resulta en un modelo SARIMA(0,1,0)(1,1,1)[12]:
\(x_t=y_{t}-y_{t-12}\) -> Diferenciación estacional
\(m_t=y_{t}-y_{t-1}\) -> Diferenciación a un paso
\(x_t=cte+ Φy_{t-12}-θ Ɛ_{t-12}+ Ɛ_t\) -> Modelo SARIMA
Sin embargo, para validar el modelo tambien se propone analizar los modelos ingenuos, cuyas metricas de error seran las minimas admitidas para el modelo SARIMA propuesto. A continuación se presentan los modelos ingenuos de promedio, ingenuo, ingenuo estacional y Drift respectivamente y posteriormente se presenta el modelo SARIMA.
#Calculo de ingenuos
#Promedio
Mod_prom <- meanf(train, h = 24)
summary(Mod_prom)##
## Forecast method: Mean
##
## Model Information:
## $mu
## [1] 356726
##
## $mu.se
## [1] 11891.01
##
## $sd
## [1] 136617.3
##
## $bootstrap
## [1] FALSE
##
## $call
## meanf(y = train, h = 24)
##
## attr(,"class")
## [1] "meanf"
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 1.412128e-11 136098.8 112798.1 -17.02068 38.3004 2.742713
## ACF1
## Training set 0.8673606
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2012 356726 180091.2 533360.8 85442.68 628009.4
## Feb 2012 356726 180091.2 533360.8 85442.68 628009.4
## Mar 2012 356726 180091.2 533360.8 85442.68 628009.4
## Apr 2012 356726 180091.2 533360.8 85442.68 628009.4
## May 2012 356726 180091.2 533360.8 85442.68 628009.4
## Jun 2012 356726 180091.2 533360.8 85442.68 628009.4
## Jul 2012 356726 180091.2 533360.8 85442.68 628009.4
## Aug 2012 356726 180091.2 533360.8 85442.68 628009.4
## Sep 2012 356726 180091.2 533360.8 85442.68 628009.4
## Oct 2012 356726 180091.2 533360.8 85442.68 628009.4
## Nov 2012 356726 180091.2 533360.8 85442.68 628009.4
## Dec 2012 356726 180091.2 533360.8 85442.68 628009.4
## Jan 2013 356726 180091.2 533360.8 85442.68 628009.4
## Feb 2013 356726 180091.2 533360.8 85442.68 628009.4
## Mar 2013 356726 180091.2 533360.8 85442.68 628009.4
## Apr 2013 356726 180091.2 533360.8 85442.68 628009.4
## May 2013 356726 180091.2 533360.8 85442.68 628009.4
## Jun 2013 356726 180091.2 533360.8 85442.68 628009.4
## Jul 2013 356726 180091.2 533360.8 85442.68 628009.4
## Aug 2013 356726 180091.2 533360.8 85442.68 628009.4
## Sep 2013 356726 180091.2 533360.8 85442.68 628009.4
## Oct 2013 356726 180091.2 533360.8 85442.68 628009.4
## Nov 2013 356726 180091.2 533360.8 85442.68 628009.4
## Dec 2013 356726 180091.2 533360.8 85442.68 628009.4#Ingenuo
Mod_ingenuo <-naive(train, h = 24)
summary(Mod_ingenuo)##
## Forecast method: Naive method
##
## Model Information:
## Call: naive(y = train, h = 24)
##
## Residual sd: 61705.8844
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 3458.725 61705.88 51715.55 -0.7718934 15.02912 1.257476 0.290612
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2012 736843 657763.7 815922.3 615901.7 857784.3
## Feb 2012 736843 625008.0 848678.0 565806.2 907879.8
## Mar 2012 736843 599873.7 873812.3 527366.5 946319.5
## Apr 2012 736843 578684.5 895001.5 494960.4 978725.6
## May 2012 736843 560016.4 913669.6 466410.0 1007276.0
## Jun 2012 736843 543139.1 930546.9 440598.5 1033087.5
## Jul 2012 736843 527618.9 946067.1 416862.4 1056823.6
## Aug 2012 736843 513173.0 960513.0 394769.3 1078916.7
## Sep 2012 736843 499605.2 974080.8 374019.1 1099666.9
## Oct 2012 736843 486772.4 986913.6 354393.0 1119293.0
## Nov 2012 736843 474566.7 999119.3 335726.0 1137960.0
## Dec 2012 736843 462904.4 1010781.6 317890.0 1155796.0
## Jan 2013 736843 451718.6 1021967.4 300782.9 1172903.1
## Feb 2013 736843 440955.5 1032730.5 284322.1 1189363.9
## Mar 2013 736843 430570.3 1043115.7 268439.3 1205246.7
## Apr 2013 736843 420525.9 1053160.1 253077.8 1220608.2
## May 2013 736843 410790.8 1062895.2 238189.2 1235496.8
## Jun 2013 736843 401338.1 1072347.9 223732.5 1249953.5
## Jul 2013 736843 392144.4 1081541.6 209672.0 1264014.0
## Aug 2013 736843 383189.7 1090496.3 195977.0 1277709.0
## Sep 2013 736843 374456.2 1099229.8 182620.3 1291065.7
## Oct 2013 736843 365928.3 1107757.7 169578.0 1304108.0
## Nov 2013 736843 357592.1 1116093.9 156828.8 1316857.2
## Dec 2013 736843 349435.3 1124250.7 144354.0 1329332.0#Ingenuo estacional
Mod_inestacional <- snaive(train, h = 24)
summary(Mod_inestacional)##
## Forecast method: Seasonal naive method
##
## Model Information:
## Call: snaive(y = train, h = 24)
##
## Residual sd: 47925.6868
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 31432.83 47925.69 41126.47 8.086573 11.68486 1 0.6135927
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2012 622713 561293.8 684132.2 528780.4 716645.6
## Feb 2012 627719 566299.8 689138.2 533786.4 721651.6
## Mar 2012 535613 474193.8 597032.2 441680.4 629545.6
## Apr 2012 446511 385091.8 507930.2 352578.4 540443.6
## May 2012 383439 322019.8 444858.2 289506.4 477371.6
## Jun 2012 405464 344044.8 466883.2 311531.4 499396.6
## Jul 2012 475544 414124.8 536963.2 381611.4 569476.6
## Aug 2012 428490 367070.8 489909.2 334557.4 522422.6
## Sep 2012 417478 356058.8 478897.2 323545.4 511410.6
## Oct 2012 559641 498221.8 621060.2 465708.4 653573.6
## Nov 2012 669767 608347.8 731186.2 575834.4 763699.6
## Dec 2012 736843 675423.8 798262.2 642910.4 830775.6
## Jan 2013 622713 535853.1 709572.9 489872.2 755553.8
## Feb 2013 627719 540859.1 714578.9 494878.2 760559.8
## Mar 2013 535613 448753.1 622472.9 402772.2 668453.8
## Apr 2013 446511 359651.1 533370.9 313670.2 579351.8
## May 2013 383439 296579.1 470298.9 250598.2 516279.8
## Jun 2013 405464 318604.1 492323.9 272623.2 538304.8
## Jul 2013 475544 388684.1 562403.9 342703.2 608384.8
## Aug 2013 428490 341630.1 515349.9 295649.2 561330.8
## Sep 2013 417478 330618.1 504337.9 284637.2 550318.8
## Oct 2013 559641 472781.1 646500.9 426800.2 692481.8
## Nov 2013 669767 582907.1 756626.9 536926.2 802607.8
## Dec 2013 736843 649983.1 823702.9 604002.2 869683.8#Drift
Mod_drift <- rwf(train, h = 24)
summary(Mod_drift)##
## Forecast method: Random walk
##
## Model Information:
## Call: rwf(y = train, h = 24)
##
## Residual sd: 61705.8844
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 3458.725 61705.88 51715.55 -0.7718934 15.02912 1.257476 0.290612
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2012 736843 657763.7 815922.3 615901.7 857784.3
## Feb 2012 736843 625008.0 848678.0 565806.2 907879.8
## Mar 2012 736843 599873.7 873812.3 527366.5 946319.5
## Apr 2012 736843 578684.5 895001.5 494960.4 978725.6
## May 2012 736843 560016.4 913669.6 466410.0 1007276.0
## Jun 2012 736843 543139.1 930546.9 440598.5 1033087.5
## Jul 2012 736843 527618.9 946067.1 416862.4 1056823.6
## Aug 2012 736843 513173.0 960513.0 394769.3 1078916.7
## Sep 2012 736843 499605.2 974080.8 374019.1 1099666.9
## Oct 2012 736843 486772.4 986913.6 354393.0 1119293.0
## Nov 2012 736843 474566.7 999119.3 335726.0 1137960.0
## Dec 2012 736843 462904.4 1010781.6 317890.0 1155796.0
## Jan 2013 736843 451718.6 1021967.4 300782.9 1172903.1
## Feb 2013 736843 440955.5 1032730.5 284322.1 1189363.9
## Mar 2013 736843 430570.3 1043115.7 268439.3 1205246.7
## Apr 2013 736843 420525.9 1053160.1 253077.8 1220608.2
## May 2013 736843 410790.8 1062895.2 238189.2 1235496.8
## Jun 2013 736843 401338.1 1072347.9 223732.5 1249953.5
## Jul 2013 736843 392144.4 1081541.6 209672.0 1264014.0
## Aug 2013 736843 383189.7 1090496.3 195977.0 1277709.0
## Sep 2013 736843 374456.2 1099229.8 182620.3 1291065.7
## Oct 2013 736843 365928.3 1107757.7 169578.0 1304108.0
## Nov 2013 736843 357592.1 1116093.9 156828.8 1316857.2
## Dec 2013 736843 349435.3 1124250.7 144354.0 1329332.0#SARIMA PROPUESTO
Mod_SARIMA <- arima(train, order = c(0,1,0), seasonal = list(order=c(1,1,1)))
summary(Mod_SARIMA)##
## Call:
## arima(x = train, order = c(0, 1, 0), seasonal = list(order = c(1, 1, 1)))
##
## Coefficients:
## sar1 sma1
## -0.0684 -0.5067
## s.e. 0.1418 0.1173
##
## sigma^2 estimated as 699241680: log likelihood = -1382.82, aic = 2771.64
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 1444.077 25107.46 16838.21 0.3182695 4.775846 0.3255929 -0.2315868#Grafico
autoplot(train)+ theme_minimal()+ ylab("Numero de turistas")+ xlab("")+ggtitle("Pronósticos ingenuos")+autolayer(meanf(train, h = 24), series="Promedio", PI=FALSE)+ autolayer(naive(train, h = 24), series="Ingenuo", PI=FALSE)+ autolayer(snaive(train, h = 24), series="Ingenuo estacional", PI=FALSE)+ autolayer(rwf(train, drift = TRUE, h = 24), series="Drift", PI=FALSE)
Al analizar los errores de ajuste de los 5 modelos propuestos, tentativamente se podria decir que el modelo que se podria ajustar mejor al pronostico para la cantidad de turistas extranjeros que llegan a la india es el modelo SARIMA con un error porcentual promedio de 4.8%, seguido del modelo ingenuo estacional, con un error porcentual promedio de 11.7%
4. Usar el modelador automático de R para estimar un modelo ARIMA o SARIMA sobre la serie. Analizar el modelo obtenido.
A continuación, se utiliza la función auto.arima para estimar un modelo sobre los datos de la ventana de entrenamiento de la serie de tiempo:
Modelo_arima <- auto.arima(train)
Modelo_arima## Series: train
## ARIMA(1,0,0)(0,1,1)[12] with drift
##
## Coefficients:
## ar1 sma1 drift
## 0.7316 -0.5383 2728.7625
## s.e. 0.0643 0.0783 360.5037
##
## sigma^2 estimated as 6.28e+08: log likelihood=-1386.66
## AIC=2781.32 AICc=2781.66 BIC=2792.47El resultado que se obtuvo hace referencia a un modelo SARIMA(1,0,0)(0,1,1)[12], modelado con la siguiente ecuación:
\(x_t=y_{t}-y_{t-12}\) -> Diferenciación estacional
\(m_t=y_{t}-y_{t-1}\) -> Diferenciación a un paso
\(m_t=2728.7+ 0.7316m_{t-1}+ 0.5383Ɛ_{t-12}+ Ɛ_t\) -> Modelo SARIMA
Este modelo nos indica en primera instancia que la serie tuvo que ser diferenciada estacionalmente para eliminar este efecto, e igualmente hace referencia a que el pronostico realizado para la llegada de turistas extranjeros a la india depende del error de estimación de hace 12 meses y tambien de la cantidad de turistas que arribaron el mes anterior, siendo este ultimo el que tiene mayor influencia sobre el pronostico.
5. Verifique los supuestos sobre los residuales del modelo obtenido automáticamente.
A continuación se presentan los correlogramas, pruebas de hipotesis y un AUTOARIMA, para determinar si los resuduales del modelo cumplen con los supuestos:
Residuales_arima <- Modelo_arima$residuals
Acf(Residuales_arima, main = "Gráfico ACF residuales ARIMA")
Pacf(Residuales_arima, main = "Gráfico PACF residuales ARIMA")
Para validar el cumplimiento de los supuestos de la serie de tiempo, se graficaron los correlogramas sobre los residuales del AUTOARIMA propuesto por el software R studio, y en ellos se pueden obeservar algunas correlaciones cercanas a las limites de significancia y particularmente en el ACF, el lag 11 resulta ser significativo y en el PACF es el lag 24, por lo tanto se procede a realizar otras pruebas para validar los supuestos.
Con el fin de analizar si los errrores son estacionarios se aplica la prueba de Dickey-Fuller que contrasta las siguentes hipótesis:
\(H_0:\)Los errores no son estacionarios
\(H_1:\)Los errores son estacionarios
adf.test(Residuales_arima)## Warning in adf.test(Residuales_arima): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: Residuales_arima
## Dickey-Fuller = -4.8692, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryA un nivel de significancia del 5% y un p-value del 0.01, se rechaza la hipotesis nula y se concluye que los residuales del modelo SARIMA cumplen el supuesto de no correlacion
Otra de las pruebas más utilizadas para verificar este supuesto es la de Ljung-Box que contrasta las siguientes hipotesis:
\(H_0\):Los errores se distribuyen de forma independiente (las correlaciones son cero)
\(H_1\):Los errores no se distribuyen de forma independiente (hay correlación)
Box.test(Residuales_arima)##
## Box-Pierce test
##
## data: Residuales_arima
## X-squared = 1.5586, df = 1, p-value = 0.2119A un nivel de significancia del 5% y un p-value del 0.2119, se acepta la hipotesis nula y se concluye que los residuales del modelo SARIMA se distribuyen de forma independiente, es decir las correlaciones son cero.
Finalmente, se evalua un modelo AUTOARIMA sobre los residuales con el proposito de establecer si estos conforman un proceso de ruido blanco y por tanto, conservan un comportamiento aleatorio de la serie, es decir lo que se busca con el AUTOARIMA es un resultado de: ARIMA(0,0,0)
auto.arima(Residuales_arima)## Series: Residuales_arima
## ARIMA(0,0,0) with zero mean
##
## sigma^2 estimated as 556622114: log likelihood=-1516.37
## AIC=3034.74 AICc=3034.77 BIC=3037.62De acuerdo con el resultado obtenido, tanto en el AUTOARIMA como en las pruebas anteriores se puede comcluir que los residuales del AUTOARIMA cumplen con los supuestos de no correlación.
6. Comparar el modelo propuesto en el punto 3 con el obtenido en el punto 4.
Teniendo en cuenta que el modelo propuesto y el modelo autocalculado arrojaron los siguientes resultados:
\(SARIMA(0,1,0)(1,1,1)[12]\) -> Modelo propuesto Punto 3
\(SARIMA(1,0,0)(0,1,1)[12]\) -> Modelo autocalculado Punto 4
Tanto el modelo propuesto como el modelo autocalculado responden a un modelo SARIMA cuyo pronostico depende del error de estimación de hace 12 meses y ademas en ninguno se considera la incidencia de los errores de estimacion inmediatamente anteriores, sin embargo, la gran diferencia radica en el parametro p que hace referencia al numero de observaciones anteriores de los cuales depende el pronostico, en este sentido, el modelo propuesto no depende de la observación del mes anterior sino del año anterior. Por otra parte, el modelo autocalculado no realiza una diferenciación a un paso como si se hizo en el modelo propuesto para controlar el fenomeno de tendencia que no se alcanzo a corregir con la diferenciación estacional.
7. Pronostique la cantidad de visitantes para la ventana de prueba y analice los indicadores de error MAE, MAPE y RMSE.
A continuación se realizan los pronosticos para los años 2013 y 2014 que fueron tomados como la ventana de test:
#Pronostico
#Pronostico Ingenuos
#Ingenuo promedio
PronosticoProm<-forecast(object = Mod_prom, h=24)
PronosticoProm## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2012 356726 180091.2 533360.8 85442.68 628009.4
## Feb 2012 356726 180091.2 533360.8 85442.68 628009.4
## Mar 2012 356726 180091.2 533360.8 85442.68 628009.4
## Apr 2012 356726 180091.2 533360.8 85442.68 628009.4
## May 2012 356726 180091.2 533360.8 85442.68 628009.4
## Jun 2012 356726 180091.2 533360.8 85442.68 628009.4
## Jul 2012 356726 180091.2 533360.8 85442.68 628009.4
## Aug 2012 356726 180091.2 533360.8 85442.68 628009.4
## Sep 2012 356726 180091.2 533360.8 85442.68 628009.4
## Oct 2012 356726 180091.2 533360.8 85442.68 628009.4
## Nov 2012 356726 180091.2 533360.8 85442.68 628009.4
## Dec 2012 356726 180091.2 533360.8 85442.68 628009.4
## Jan 2013 356726 180091.2 533360.8 85442.68 628009.4
## Feb 2013 356726 180091.2 533360.8 85442.68 628009.4
## Mar 2013 356726 180091.2 533360.8 85442.68 628009.4
## Apr 2013 356726 180091.2 533360.8 85442.68 628009.4
## May 2013 356726 180091.2 533360.8 85442.68 628009.4
## Jun 2013 356726 180091.2 533360.8 85442.68 628009.4
## Jul 2013 356726 180091.2 533360.8 85442.68 628009.4
## Aug 2013 356726 180091.2 533360.8 85442.68 628009.4
## Sep 2013 356726 180091.2 533360.8 85442.68 628009.4
## Oct 2013 356726 180091.2 533360.8 85442.68 628009.4
## Nov 2013 356726 180091.2 533360.8 85442.68 628009.4
## Dec 2013 356726 180091.2 533360.8 85442.68 628009.4accuracy(PronosticoProm, test)## ME RMSE MAE MPE MAPE MASE
## Training set 1.412128e-11 136098.8 112798.1 -17.02068 38.30040 2.742713
## Test set 2.076634e+05 244891.6 207663.4 33.40625 33.40625 5.049386
## ACF1 Theil's U
## Training set 0.8673606 NA
## Test set 0.6741652 2.316693#Ingenuo
PronosticoI<-forecast(object = Mod_ingenuo, h=24)
PronosticoI## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2012 736843 657763.7 815922.3 615901.7 857784.3
## Feb 2012 736843 625008.0 848678.0 565806.2 907879.8
## Mar 2012 736843 599873.7 873812.3 527366.5 946319.5
## Apr 2012 736843 578684.5 895001.5 494960.4 978725.6
## May 2012 736843 560016.4 913669.6 466410.0 1007276.0
## Jun 2012 736843 543139.1 930546.9 440598.5 1033087.5
## Jul 2012 736843 527618.9 946067.1 416862.4 1056823.6
## Aug 2012 736843 513173.0 960513.0 394769.3 1078916.7
## Sep 2012 736843 499605.2 974080.8 374019.1 1099666.9
## Oct 2012 736843 486772.4 986913.6 354393.0 1119293.0
## Nov 2012 736843 474566.7 999119.3 335726.0 1137960.0
## Dec 2012 736843 462904.4 1010781.6 317890.0 1155796.0
## Jan 2013 736843 451718.6 1021967.4 300782.9 1172903.1
## Feb 2013 736843 440955.5 1032730.5 284322.1 1189363.9
## Mar 2013 736843 430570.3 1043115.7 268439.3 1205246.7
## Apr 2013 736843 420525.9 1053160.1 253077.8 1220608.2
## May 2013 736843 410790.8 1062895.2 238189.2 1235496.8
## Jun 2013 736843 401338.1 1072347.9 223732.5 1249953.5
## Jul 2013 736843 392144.4 1081541.6 209672.0 1264014.0
## Aug 2013 736843 383189.7 1090496.3 195977.0 1277709.0
## Sep 2013 736843 374456.2 1099229.8 182620.3 1291065.7
## Oct 2013 736843 365928.3 1107757.7 169578.0 1304108.0
## Nov 2013 736843 357592.1 1116093.9 156828.8 1316857.2
## Dec 2013 736843 349435.3 1124250.7 144354.0 1329332.0accuracy(PronosticoI, test)## ME RMSE MAE MPE MAPE MASE
## Training set 3458.725 61705.88 51715.55 -0.7718934 15.02912 1.257476
## Test set -172453.583 215842.65 180859.17 -37.5541452 38.59215 4.397634
## ACF1 Theil's U
## Training set 0.2906120 NA
## Test set 0.6741652 2.809457#Ingenuo Estacional
PronosticoE<-forecast(object = Mod_inestacional, h=24)
PronosticoE## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2012 622713 561293.8 684132.2 528780.4 716645.6
## Feb 2012 627719 566299.8 689138.2 533786.4 721651.6
## Mar 2012 535613 474193.8 597032.2 441680.4 629545.6
## Apr 2012 446511 385091.8 507930.2 352578.4 540443.6
## May 2012 383439 322019.8 444858.2 289506.4 477371.6
## Jun 2012 405464 344044.8 466883.2 311531.4 499396.6
## Jul 2012 475544 414124.8 536963.2 381611.4 569476.6
## Aug 2012 428490 367070.8 489909.2 334557.4 522422.6
## Sep 2012 417478 356058.8 478897.2 323545.4 511410.6
## Oct 2012 559641 498221.8 621060.2 465708.4 653573.6
## Nov 2012 669767 608347.8 731186.2 575834.4 763699.6
## Dec 2012 736843 675423.8 798262.2 642910.4 830775.6
## Jan 2013 622713 535853.1 709572.9 489872.2 755553.8
## Feb 2013 627719 540859.1 714578.9 494878.2 760559.8
## Mar 2013 535613 448753.1 622472.9 402772.2 668453.8
## Apr 2013 446511 359651.1 533370.9 313670.2 579351.8
## May 2013 383439 296579.1 470298.9 250598.2 516279.8
## Jun 2013 405464 318604.1 492323.9 272623.2 538304.8
## Jul 2013 475544 388684.1 562403.9 342703.2 608384.8
## Aug 2013 428490 341630.1 515349.9 295649.2 561330.8
## Sep 2013 417478 330618.1 504337.9 284637.2 550318.8
## Oct 2013 559641 472781.1 646500.9 426800.2 692481.8
## Nov 2013 669767 582907.1 756626.9 536926.2 802607.8
## Dec 2013 736843 649983.1 823702.9 604002.2 869683.8accuracy(PronosticoE, test)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 31432.83 47925.69 41126.47 8.086573 11.684863 1.0000000 0.6135927
## Test set 38620.92 49792.10 40123.58 6.342514 6.708974 0.9756146 0.2647599
## Theil's U
## Training set NA
## Test set 0.4767471#Drift
PronosticoD<-forecast(object = Mod_drift, h=24)
PronosticoD## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2012 736843 657763.7 815922.3 615901.7 857784.3
## Feb 2012 736843 625008.0 848678.0 565806.2 907879.8
## Mar 2012 736843 599873.7 873812.3 527366.5 946319.5
## Apr 2012 736843 578684.5 895001.5 494960.4 978725.6
## May 2012 736843 560016.4 913669.6 466410.0 1007276.0
## Jun 2012 736843 543139.1 930546.9 440598.5 1033087.5
## Jul 2012 736843 527618.9 946067.1 416862.4 1056823.6
## Aug 2012 736843 513173.0 960513.0 394769.3 1078916.7
## Sep 2012 736843 499605.2 974080.8 374019.1 1099666.9
## Oct 2012 736843 486772.4 986913.6 354393.0 1119293.0
## Nov 2012 736843 474566.7 999119.3 335726.0 1137960.0
## Dec 2012 736843 462904.4 1010781.6 317890.0 1155796.0
## Jan 2013 736843 451718.6 1021967.4 300782.9 1172903.1
## Feb 2013 736843 440955.5 1032730.5 284322.1 1189363.9
## Mar 2013 736843 430570.3 1043115.7 268439.3 1205246.7
## Apr 2013 736843 420525.9 1053160.1 253077.8 1220608.2
## May 2013 736843 410790.8 1062895.2 238189.2 1235496.8
## Jun 2013 736843 401338.1 1072347.9 223732.5 1249953.5
## Jul 2013 736843 392144.4 1081541.6 209672.0 1264014.0
## Aug 2013 736843 383189.7 1090496.3 195977.0 1277709.0
## Sep 2013 736843 374456.2 1099229.8 182620.3 1291065.7
## Oct 2013 736843 365928.3 1107757.7 169578.0 1304108.0
## Nov 2013 736843 357592.1 1116093.9 156828.8 1316857.2
## Dec 2013 736843 349435.3 1124250.7 144354.0 1329332.0accuracy(PronosticoD, test)## ME RMSE MAE MPE MAPE MASE
## Training set 3458.725 61705.88 51715.55 -0.7718934 15.02912 1.257476
## Test set -172453.583 215842.65 180859.17 -37.5541452 38.59215 4.397634
## ACF1 Theil's U
## Training set 0.2906120 NA
## Test set 0.6741652 2.809457#Pronostico Modelo propuesto
PronosticoP<-forecast(object = Mod_SARIMA, h=24)
PronosticoP## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2012 684498.6 650610.3 718386.9 632670.9 736326.3
## Feb 2012 688562.0 640636.7 736487.3 615266.6 761857.4
## Mar 2012 619130.3 560434.1 677826.6 529362.2 708898.5
## Apr 2012 511632.0 443855.4 579408.6 407976.6 615287.4
## May 2012 459064.0 383287.5 534840.6 343173.8 574954.3
## Jun 2012 494972.5 411963.5 577981.6 368021.1 621923.9
## Jul 2012 571649.3 481989.3 661309.4 434526.2 708772.5
## Aug 2012 523735.5 427884.9 619586.1 377144.7 670326.4
## Sep 2012 493005.3 391340.4 594670.2 337522.2 648488.4
## Oct 2012 627196.1 520031.9 734360.3 463302.5 791089.6
## Nov 2012 725945.4 613550.6 838340.2 554052.4 897838.4
## Dec 2012 790750.6 673358.1 908143.1 611214.2 970287.0
## Jan 2013 738068.0 611131.7 865004.4 543935.6 932200.4
## Feb 2013 742195.9 606384.8 878007.0 534490.7 949901.0
## Mar 2013 671213.9 527073.4 815354.4 450770.0 891657.7
## Apr 2013 564973.4 412959.2 716987.5 332487.7 797459.0
## May 2013 511687.2 352187.5 671186.8 267753.5 755620.8
## Jun 2013 546646.4 379997.2 713295.6 291778.4 801514.4
## Jul 2013 622872.2 449367.7 796376.6 357520.0 888224.3
## Aug 2013 575017.1 394918.3 755116.0 299579.7 850454.6
## Sep 2013 545635.1 359174.8 732095.4 260468.7 830801.5
## Oct 2013 680371.0 487759.3 872982.7 385796.8 974945.2
## Nov 2013 779898.2 581325.6 978470.8 476207.6 1083588.8
## Dec 2013 844858.7 640498.9 1049218.4 532317.4 1157400.0accuracy(PronosticoP, test)## ME RMSE MAE MPE MAPE MASE
## Training set 1444.077 25107.46 16838.21 0.3182695 4.775846 0.4094253
## Test set -61176.368 69824.58 61176.37 -12.4887693 12.488769 1.4875182
## ACF1 Theil's U
## Training set -0.2315868 NA
## Test set 0.5582129 0.8892055#Pronostico Modelo AutoCalculado
PronosticoA<-forecast(object = Modelo_arima,h=24)
PronosticoA## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2012 672447.0 640331.8 704562.2 623331.0 721563.0
## Feb 2012 668114.5 628321.7 707907.2 607256.7 728972.2
## Mar 2012 590515.3 547168.1 633862.4 524221.6 656809.0
## Apr 2012 480175.6 435040.7 525310.4 411147.8 549203.3
## May 2012 423305.9 377242.7 469369.2 352858.3 493753.6
## Jun 2012 455384.0 408831.4 501936.6 384187.9 526580.1
## Jul 2012 529589.9 482777.5 576402.3 457996.4 601183.3
## Aug 2012 480550.2 433599.2 527501.1 408744.9 552355.4
## Sep 2012 450476.3 403451.4 497501.2 378557.9 522394.6
## Oct 2012 583823.8 536759.3 630888.2 511844.9 655802.6
## Nov 2012 682229.5 635143.9 729315.1 610218.3 754240.7
## Dec 2012 746414.9 699318.0 793511.8 674386.4 818443.4
## Jan 2013 688238.0 638629.7 737846.3 612368.6 764107.4
## Feb 2013 688455.5 637553.8 739357.3 610608.0 766303.1
## Mar 2013 614185.3 562604.5 665766.1 535299.3 693071.3
## Apr 2013 506281.1 454340.5 558221.7 426844.9 585717.4
## May 2013 451193.4 399061.2 503325.5 371464.1 530922.6
## Jun 2013 484575.1 432340.7 536809.6 404689.5 564460.8
## Jul 2013 559734.8 507445.8 612023.9 479765.6 639704.1
## Aug 2013 511393.0 459074.7 563711.3 431379.0 591406.9
## Sep 2013 481829.6 429495.7 534163.5 401791.7 561867.4
## Oct 2013 615550.6 563208.3 667892.9 535500.0 695601.3
## Nov 2013 714229.6 661882.9 766576.4 634172.2 794287.1
## Dec 2013 778615.0 726265.8 830964.2 698553.8 858676.1accuracy(PronosticoA, test)## ME RMSE MAE MPE MAPE MASE
## Training set -407.8078 23592.84 16058.87 -0.8742041 4.619900 0.3904752
## Test set -12998.4104 31663.02 28272.36 -3.4813304 5.625108 0.6874493
## ACF1 Theil's U
## Training set -0.1086629 NA
## Test set 0.4358972 0.3840704En la siguiente grafica podemos apreciar que los pronosticos que mejor se ajustan al comportamiento de la serie de tiempo son los modelos SARIMA y el modelo ingenuo estacional, pues se asemejan mas al estilo de la ventana de train, siguiendo la tendencia y estacionalidad propias de esta serie.
autoplot(train)+ theme_minimal()+ ylab("Numero de turistas")+ xlab("")+ggtitle("Pronósticos ingenuos")+autolayer(forecast(object = Mod_prom, h=24), series="Ingenuo Promedio", PI=FALSE)+autolayer(forecast(object = Mod_drift, h=24), series="Ingenuo Drift", PI=FALSE)+autolayer(forecast(object = Mod_inestacional, h=24), series="Ingenuo Estacional", PI=FALSE)+autolayer(forecast(object = Mod_ingenuo, h=24), series="Ingenuo", PI=FALSE)+autolayer(forecast(object = Mod_SARIMA, h=24), series="SARIMA Propuesto", PI=FALSE)+autolayer(forecast(object = Modelo_arima, h=24), series="SARIMA Autocalculado", PI=FALSE)
Teniendo en cuenta los resultados obtenidos en los errores de pronostico, se toma como valor minimo para los modelos SARIMA un error porcentual promedio de 6.7% correspondiente al error del modelo ingenuo estacional. Al comparar el modelo propuesto se encuenta un error porcentual promedio de 12.5% que al ser mayor al valor utilizado como benchemark (Ingenuo estacional), se descarta su aplicacion ya que lo que se busca es minimizar esta metrica. Por su parte el modelo autocalculado obtuvo un error porcentual promedio del 5.6%, por lo tanto se propone basarse en este modelo para las proyecciones en la agencia de viajes.
8. Utilice los resultados para resolver el dilema propuesto en el caso.
Teniendo en cuenta la situación en la cual se encuentra actualmente el INDIA HOLIDAY BUREAU, se realizaron pruebas para identificar un modelo que se ajustara a la situacion turistica caracteristica de esta compañia, con lo que se llego a la conclusión de que el modelo de pronostico que minimiza las metricas de error es un modelo SARIMA (1,0,0)(0,1,1)[12], el cual tiene en cuenta el error de pronostico del año anterior y la cantidad de turistas que llegaron el mes anterior.
plot(PronosticoA)
Como se muestra en la grafica anterior, este modelo permite predecir la llegada de los turistas extranjeros siguiendo con el comportamiento de tendencia de la serie de tiempo y tambien permite visualizar el efecto estacional anual en donde los periodos con mayor flujo de turistas corresponden a las epocas vacacionales, por ende se recomienda su aplicación para tener una mejor estimacion de la fuerza laboral con la que la agencia deberia contar para dichos periodos y asi mismo afrecer un servicio de calidad a sus clientes y aumentar sus rubros e indicadores de satisfacción.