NETFLIX INC.
En el presente trabajo se va a desarrollar el análisis del comportamiento que presento la cotización de la empresa dentro del mercado de capitales, para ello en esta primera sección se hará una breve descripción de la empresa y se obtendrán los datos que serán utilizado en lo que sigue del trabajo.
Descripción de la empresa.
Netflix, Inc. es el servicio de entretenimiento de transmisión por suscripción líder en el mundo con más de 167 millones de pagos de membresías en más de 190 países, disfrutando de series de televisión, documentales y largometrajes en una amplia variedad de géneros e idiomas. Los miembros pueden mira todo lo que quieran, en cualquier momento, en cualquier lugar, en cualquier pantalla conectada a Internet. Los miembros pueden reproducir, pausar y reanudar la visualización, todo sin comerciales. Además, más de dos millones de miembros en los Estados Unidos (“EE.UU.”) se suscriben al antiguo servicio de DVD por correo.
La empresa es pionera en la entrega de entretenimiento en streaming, lanzando el servicio de streaming en 2007. Desde este lanzamiento, se ha ido desarrollando un ecosistema para pantallas conectadas a Internet y han agregado cantidades cada vez mayores de contenido que permiten a los consumidores disfrutar del entretenimiento directamente en sus pantallas conectadas a Internet. Como resultado de estos esfuerzos, han experimentado una creciente aceptación e interés por parte de los consumidores en la entrega de entretenimiento en streaming.
La estrategia principal de la compañía es hacer el nuestro negocio de membresía de transmisión a nivel mundial dentro de los parámetros de su objetivo sobre el margen operativo. Están continuamente mejorando la experiencia de sus miembros al expandir el contenido de transmisión con un enfoque de programación de contenido que deleita a los miembros y atrae a su vez nuevos miembros. Además, continúan mejorando la interfaz de usuario y ampliando el servicio de transmisión a más pantallas conectadas a Internet.
Al 31 de diciembre de 2019, la empresa tenía aproximadamente 8,600 empleados de tiempo completo. Además, obtuvo una utilidad neta de 1,866,916 de dólares.
Obtención de datos.
A continuación, se obtendrá los datos de la empresa del periodo 2017-01-01 y 2020-01-01, estos datos tendrán una frecuencia diaria. Para lograr esto, se usará los siguientes códigos:
getSymbols.yahoo("NFLX", env = globalenv(), from = as.Date("2016-01-01"), to = as.Date("2020-01-01"), periodicity = 'monthly')## [1] "NFLX"Como se observa, se utiliza la función getSymbols.yahoo para poder obtener los datos de la empresa indicando en el parámetro from la fecha desde se iniciará la descarga de datos, en el parámetro to colocamos la fecha final del periodo que se analizara y en el parámetro periodicity colocaremos daily pues se trabajaran con datos diarios.
Para poder confirmar que se tuvieron los datos deseados, se utilizara la función head para mostrar los 5 primeros datos descargados y luego se ejecutara la función tail para exponer los últimos datos de la serie descargada, además se utilizara la función length para determinar el número de datos descargados, para ello se ejecutara el siguiente codigo:
## NFLX.Open NFLX.High NFLX.Low NFLX.Close NFLX.Volume NFLX.Adjusted
## 2016-01-01 109.00 122.18 90.11 91.84 488193200 91.84
## 2016-02-01 91.79 97.48 79.95 93.41 389268900 93.41
## 2016-03-01 94.58 104.91 93.61 102.23 311333700 102.23
## 2016-04-01 102.93 111.85 88.21 90.03 340174300 90.03
## 2016-05-01 90.41 104.00 85.74 102.57 264997900 102.57## NFLX.Open NFLX.High NFLX.Low NFLX.Close NFLX.Volume NFLX.Adjusted
## 2019-08-01 324.25 328.58 287.20 293.75 137076700 293.75
## 2019-09-01 290.82 301.55 252.28 267.62 175411300 267.62
## 2019-10-01 267.35 308.75 257.01 287.41 231556400 287.41
## 2019-11-01 288.70 316.82 281.14 314.66 113645900 314.66
## 2019-12-01 314.39 338.00 292.02 323.57 124723600 323.57## [1] 288Por lo que se ve que se la serie cuenta con 3018 datos. A continuación, se mostrará la gráfica de la serie de tiempo descargada.
Por lo que se puede observar, la cotización de la empresa se ha mantenido altamente comercializada y a partir del año 2019 se puede observar que el precio cotizado se mantuvo alrededor de 350 dólares por acción. Para el 31 de diciembre del 2019, la acción cotizaba en \(323.57\$\) dólares y el volumen de comercialización es de \(3,713,300\$\).
Por otro lado, se va separar los precios al cierre ajustados de toda la serie descargada, ya que estos precios no están afectados por el pago de dividendos. Para ello se ejecuta el siguiente código, donde se guardarán los datos en el objeto pre.net.
Luego de esto, se graficará los precios ajustados. Para ello se ejecuta el siguiente código:
Además, cabe resaltar que se analizará los rendimientos de la acción de Netflix y a partir de estos se proyectaran, para lograr esto se utilizara el siguiente codigo:
Luego se elimina el valor NA que se da como resultado de obtener la primera diferenciacion de logaritmo Y se grafica la serie obtenida:
Por lo tanto, se utilizaran los retornos de las acciones de Netflix entre el periodo 2016-02-01 hasta 2019-12-01 y a partir de estos se proyectaran los retornos futuros.
ANALISIS DE DATOS
En esta sección se aplicará la metodología de Box-Jenkis para identificar los posibles modelos ARIMA a partir de los datos obtenidos.
Por esa razón, primero se analizarán los componentes de la serie temporal, segundo se estudiará el orden de integración que presenta la serie, tercero se estima un modelo ARIMA o ARMA y se transformará a la serie estacionaria, finalmente se graficará y analizará las funciones de autocorrelación simple y parcial con el fin de elegir posibles modelos ARIMA.
Componentes de la serie temporal.
Primero se definirá a la serie de retorno como una variable temporal además indicando que es de frecuencia diaria, para ello se desarrolla el siguiente código:
Por otro lado, se apartan una porcion de la serie para poder efectuar en lo que sigue del trabajo los modeloss y las proyecciones de los retornos de las acciones de Netflix.
Luego de ello, se procede a descomponer la serie de retornos con la función decompose, tal como se muestra en el siguiente código
Al obtener los diferentes componentes que posee el retorno de Netflix.
Componente Tendencial
Como se observa en el gráfico, se podría decir está bien definido hacia una tendencia negativa, sobre todo a partir del año 2018, pero se podría concluir que de manera gráfica el componente tendencial de los retornos de Netflix es negativo a largo plazo.
Componente Estacional
Como se observa en el siguiente gráfico, el componente estacional está presente en gran parte de la serie, es decir, no hay un periodo amplio entre cada ciclo estacional, lo cual puede decir que los rendimiento que presenta las acciones de Netflix están muy influenciados por el componente estacional.
Integración de la serie
En esta sección se determinará el grado de integración que necesita la serie para volverse estacionario, en caso esto sea necesario.
Iniciaremos aplicando ciertos test de raíz unitaria, con el objetivo de determinar si los retornos de las acciones de Netflix son estacionara.
Test de Raíz Unitaria
Primero se va a determinar si la serie presenta raíz unitaria o no, es decir si es o no estacionaria. Para ello se aplicará el test de Dickey - Fuller Aumentado.
Pero se ira analizando tomando en cuenta sus componentes para llegar luego a una conclusión. Primero se utilizará el test tomando el componente tendencial, luego se utilizará el test solo tomando en cuenta la constante y finalmente no se tomará en cuenta ninguna de los dos, es decir ni la constante ni el componente tendencial para aplicar el test de raíz unitaria.
Test de Raíz Unitaria con componente Tendencial
Primero se aplicó el test de raíz unitaria tomando en cuenta el componente tendencial, para ello se ejecutó el siguiente código, donde se especifica en el parámetro type que se tomara en cuenta el componente tendencial. Adicionalmente se utiliza el comando summary para poder observar de manera ordenada los resultados obtenidos.
##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression trend
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.232501 -0.073020 0.000395 0.050890 0.307416
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.057732 0.036939 1.563 0.126
## z.lag.1 -1.181790 0.240756 -4.909 1.5e-05 ***
## tt -0.001135 0.001296 -0.875 0.386
## z.diff.lag 0.009622 0.157328 0.061 0.952
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1101 on 41 degrees of freedom
## Multiple R-squared: 0.5859, Adjusted R-squared: 0.5556
## F-statistic: 19.34 on 3 and 41 DF, p-value: 5.729e-08
##
##
## Value of test-statistic is: -4.9087 8.0494 12.0738
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -4.15 -3.50 -3.18
## phi2 7.02 5.13 4.31
## phi3 9.31 6.73 5.61Por lo que, como se observa dado que el componente tendencial no es significativo dentro del modelo, no se puede tomar en cuenta el resultado del test de raíz unitaria.
Test de Raíz Unitaria con constante.
Dado que en la ejecución anterior, el componente tendencial no tiene significancia por lo que ahora se apicara el test tomando en cuenta solo una constante, para ello se ejecutara el siguiente código colocando en el parámetro type la categoría drift, para indicar que solo se desea aplicar el test considerando una constante. Adicionalmente se utiliza el comando summary para poder observar de manera ordenada los resultados obtenidos.
##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression drift
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.242868 -0.059068 0.005724 0.056798 0.310310
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.02926 0.01747 1.675 0.101
## z.lag.1 -1.13459 0.23399 -4.849 1.74e-05 ***
## z.diff.lag -0.01523 0.15431 -0.099 0.922
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1098 on 42 degrees of freedom
## Multiple R-squared: 0.5781, Adjusted R-squared: 0.5581
## F-statistic: 28.78 on 2 and 42 DF, p-value: 1.344e-08
##
##
## Value of test-statistic is: -4.8489 11.7563
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.58 -2.93 -2.60
## phi1 7.06 4.86 3.94Igualmente, al rechazar el estimado del parámetro no se toma en cuenta los resultados del test estadístico.
Test de Raíz Unitaria sin constante y tendencia.
A continuación, se aplicará el test sin tomar en cuenta la constante y el componente tendencial. Para ello se coloca el none dentro del parámetro type. Adicionalmente, se utiliza la función summary para poder obtener los resultados del test.
##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression none
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.22060 -0.04714 0.02941 0.07782 0.34117
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 -0.99723 0.22370 -4.458 5.84e-05 ***
## z.diff.lag -0.08257 0.15208 -0.543 0.59
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.112 on 43 degrees of freedom
## Multiple R-squared: 0.55, Adjusted R-squared: 0.5291
## F-statistic: 26.28 on 2 and 43 DF, p-value: 3.501e-08
##
##
## Value of test-statistic is: -4.4579
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.62 -1.95 -1.61Dado que los parámetros estimados no son significativos estadísticamente, no podríamos tomar en cuenta el resultado de este test estadístico.
Test Raíz Unitaria
Dado que los test aplicados anteriormente no fueron efectivos, pues los resultados no podrían ser estadísticamente significativos, se utilizara la función adf.test el cual es otra forma de aplicar el test de Dickey - Fuller Aumentado.
##
## Augmented Dickey-Fuller Test
##
## data: pre.ts
## Dickey-Fuller = -3.3353, Lag order = 3, p-value = 0.07769
## alternative hypothesis: stationaryPor lo tanto, la serie log de los precios de la acción de Netflix no es estacionaria.
Grado de Integración
Dado que los retornos de las acciones de Netflix no es estacionaria, será necesario determinar su grado de integración con el objetivo de volver a la serie estacionaria. Para ello se obtendrá la primera diferencia de los retornos y se grafica el resultado, usando el siguiente código:
Obtenida la primera diferenciación, se aplicará nuevamente el test de Dickey - Fuller Aumentado, para ello utilizamos la función adf.test como se muestra en el siguiente código:
## Warning in adf.test(dif1.ren): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: dif1.ren
## Dickey-Fuller = -5.5104, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationaryComo se observa, el test estadístico nos sugiere rechazar la hipótesis nula pudiendo concluir que, al obtener las primeras diferencias de los retornos de las acciones de Netflix, esta se vuelve estacionaria.
En suma, la serie de los retornos de las acciones de Netflix no es estacionaria y necesita ser diferenciada para poder ser estacionaria, en otras palabras el orden de integración es 1 (\(I_{(1)}\)).
Por lo dicho anteriormente el orden de integración es 1 (\(I_{(1)}\)), entonces tenemos que estimar un modelo ARIMA, pero para saber cómo se puede estimar dicho modelo usamos el código de Autoarima que busca el “mejor” modelo ARIMA para dicha serie de tiempo.
m1 = auto.arima(pre.ts1 , d=1, D=0, max.p = 12, max.q = 12, max.P = 1, max.Q = 1,
max.d = 0, max.D = 0, stationary = FALSE, seasonal = TRUE, ic = c("bic"))
summary(m1)## Series: pre.ts1
## ARIMA(2,1,0)
##
## Coefficients:
## ar1 ar2
## -0.8582 -0.4667
## s.e. 0.1463 0.1482
##
## sigma^2 estimated as 0.01603: log likelihood=24.54
## AIC=-43.07 AICc=-42.34 BIC=-38.24
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.002433068 0.1215186 0.08740656 72.60342 195.5708 0.6276713
## ACF1
## Training set -0.1259215Por lo que se puede observar que el mejor modelo ARIMA es (2,1,0) donde calcula un AR(2), un orden de integración 1 y MA nulo. Donde los coeficiente del AR(1) es de -0.8582 y del AR(2) es de -0.4667; por último el AIC del modelo ARIMA(2,1,0) es de -43.07.
A continuación, se graficará y analizará la gráfica de autocorrelación simple y autocorrelación parcial con los siguientes códigos:
Como sabemos por la teoría se analiza las líneas que salen de las bandas y de acuerdo con eso se realiza un determinado número de (p) y (q) para un AR y MA respectivamente, añadiendo que la autocorrelación parcial es para el AR y la autocorrelación es para el MA.
A raíz del análisis se puede proponer dos modelos ARIMA que serían: ARIMA(0,1,1)(0,1,0) y ARIMA(1,1,1)(0,1,0).
Se explica la propuesta de estos dos modelos porque el AR en la auctocorrelación parcial no cuenta con rezagos que salgan de su banda de confianza, en cambio en el MA por ende en la autocorrelación el primer rezago sale de su banda. Por esto, a raíz de lo explicado anteriormente se propone como primer modelo tentativo un ARIMA(0,1,1) y como segundo modelo tentativo un ARIMA(1,1,1).
Arima Elegidos
En este apartado se va estimar los tres modelos ARIMA tentativos del anterior apartado, los tres ARIMAS planteados fueron:
- AUTOARIMA
- ARIMA(0,1,1)(0,1,0)
- ARIMA(1,1,1)(0,1,0)
Auto.arima
m1 = auto.arima(pre.ts1 , d=1, D=0, max.p = 12, max.q = 12, max.P = 1, max.Q = 1,
max.d = 0, max.D = 0, stationary = FALSE, seasonal = TRUE, ic = c("bic"))Modelo 2 según criterio del investigador ARIMA(0,1,1)(0,1,0)
Modelo 3 según criterio del investigador ARIMA(1,1,1)(0,1,0)
Resumen para ver los criterios de información de los modelos propuestos.
## Series: pre.ts1
## ARIMA(2,1,0)
##
## Coefficients:
## ar1 ar2
## -0.8582 -0.4667
## s.e. 0.1463 0.1482
##
## sigma^2 estimated as 0.01603: log likelihood=24.54
## AIC=-43.07 AICc=-42.34 BIC=-38.24
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.002433068 0.1215186 0.08740656 72.60342 195.5708 0.6276713
## ACF1
## Training set -0.1259215##
## Call:
## arima(x = pre.ts1, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 0), period = 12))
##
## Coefficients:
## ma1
## -0.9493
## s.e. 0.1819
##
## sigma^2 estimated as 0.01846: log likelihood = 13.31, aic = -22.62
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.01254894 0.1101885 0.06421174 4.565272 123.8264 0.461108
## ACF1
## Training set -0.3274939##
## Call:
## arima(x = pre.ts1, order = c(1, 1, 1), seasonal = list(order = c(0, 1, 0), period = 12))
##
## Coefficients:
## ar1 ma1
## -0.3464 -0.8369
## s.e. 0.1999 0.1471
##
## sigma^2 estimated as 0.01688: log likelihood = 14.62, aic = -23.25
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.01149771 0.1053942 0.05992389 5.83656 112.3321 0.4303167
## ACF1
## Training set -0.09239702Se puede observar en el resumen de los tres modelos tentativos, el criterio de información de AIC(AKAIKE), es una medida de la calidad relativa de un modelo estadístico, para un conjunto dado de datos. Como tal, el AIC proporciona un medio para la selección del modelo. AIC maneja un trade-off entre la bondad de ajuste del modelo y la complejidad del modelo. Se basa en la entropía de información: se ofrece una estimación relativa de la información perdida cuando se utiliza un modelo determinado para representar el proceso que genera los datos. AIC no proporciona una prueba de un modelo en el sentido de probar una hipótesis nula, es decir AIC no puede decir nada acerca de la calidad del modelo en un sentido absoluto. Si todos los modelos candidatos encajan mal, AIC no dará ningún aviso de ello.
Para el primer modelo se tiene un AIC de -43.07, considerando que este modelo es un ARIMA(2,1,0)(0,1,0), este fue obtenido por el autoarima. El segundo modelo de ARIMA(0,1,1)(0,1,0), tiene un AIC de -22.62 este es mucho mayor al primer modelo que nos dio al ejecutar el autoarima, el tercer modelo que es un ARIMA(1,1,1)(0,1,0), tiene un AIC de -23.25, en este caso dado a los criterios de información se seleccionarían el modelo ARIMA(2,1,0) (0,1,0) y el tercer modelo ARIMA(1,1,1)(0,1,1), se va preferir siempre el que tenga menor criterio de información, en nuestro caso el que sea más negativo, en este caso el primer modelo, pero para fines de seguir el análisis se va elegir estos dos modelos.
Graficando
dif2.pre1 = m1$fitted
dif2.pre1<-ts(dif2.pre1,start = c(2016,2),frequency=12)
dif2.pre2 = fitted(m2) #se usa el comando fitted porque a diferencia del auto.arima, no es automático
dif2.pre2<-ts(dif2.pre2,start = c(2016,2),frequency=12)
dif2.pre3 = fitted(m3)
dif2.pre3<-ts(dif2.pre3,start = c(2016,2),frequency=12)
Model.Series1 = cbind(pre.ts1, dif2.pre1)
Model.Series2 = cbind(pre.ts1, dif2.pre2)
Model.Series3 = cbind(pre.ts1, dif2.pre3)
names(Model.Series1) = c('Observado', 'Estimado m1')
names(Model.Series2) = c('Observado', 'Estimado m2')
names(Model.Series3) = c('Observado', 'Estimado m3')
par(mfrow=c(3,1))
chart.TimeSeries(Model.Series1, main = "auto.arima:ARIMA(2,1,0)(0,1,0)", legend.loc = "bottomleft")
chart.TimeSeries(Model.Series2, main = "arima(0,1,1)(0,1,0)", legend.loc = "bottomleft")
chart.TimeSeries(Model.Series3, main = "arima(1,1,1)(0,1,0)", legend.loc = "bottomleft")
Observando los gráficos de los tres modelos en este caso, se puede observar que ninguno sigue el comportamiento esperado, ni los modelos que habíamos seleccionado que tenían un menor criterio de información. Para seguir con el análisis en los siguientes apartados se seguirá con la selección del modelo 1 y el modelo 3, ya que estos presentaron un menor criterio de información de AKAIKE. Por ello se seguirá con los demás análisis de estos modelos seleccionados.
PRUEBAS DE DIAGNOSTICO.
En esta sección se analizarán los residuos de los modelos 1 y 3, con el objetivo de conocer si estos poseen los parámetros estadísticos de un ruido blanco, es decir que tengan media cero y varianza constante, y además no presentan autocorrelación. De confirmar que uno o ambos de los modelos presentan residuos un comportamiento como el de ruido blanco, se podrá concluir que los modelos son eficientes y se podrán generar proyecciones a partir de estos.
Por esa razón esta sección estará dividida en dos partes, en la primera se estudiará la autocorrelación que presentan los residuos de cada modelo y en la segunda sección se evaluara si los ruidos siguen una distribución normal para finalmente poder concluir si los modelos son válidos estadísticamente.
Diagnostico de autocorrelación de los residuos (Ljung-Box test)
Para poder analizar la autocorrelación serial de los residuos, se usara la función checkresiduals() donde se podrá observar la función de autocorrelación, como también la distribución de los residuos, para ello se ejecutara los siguientes comandos:
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,1,0)
## Q* = 6.3701, df = 6, p-value = 0.383
##
## Model df: 2. Total lags used: 8
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,1)(0,1,0)[12]
## Q* = 8.3923, df = 6, p-value = 0.2107
##
## Model df: 2. Total lags used: 8Observando los resultados, tomando como análisis el grafico de autocorrelación simple de los residuos, se puede decir que el Modelo 3 presenta menor grado de autocorrelación serial que los presentado por el Modelos 1. Por otro lado, al analizar la gráfica de distribución de los residuos que se compara con la distribución de una normal, se puede notar que la distribución de los residuos del Modelo 1 sigue la forma de una distribución normal y para el Modelo 3 la distribución de sus residuos, al menos desde el análisis gráfico, parece no seguir una distribución normal.
Normalidad (Jarque-Bera)
A continuación, se analizará la distribución de los residuos con el objetivo de determinar si estos siguen una distribución normal o no. Par a ello se ejecutara la función jb.norm.test para aplicar el test estadístico de Jarque-Bera:
##
## Jarque-Bera test for normality
##
## data: m1$residuals
## JB = 7.4563, p-value = 0.024##
## Jarque-Bera test for normality
##
## data: m3$residuals
## JB = 17.778, p-value = 0.0025Cabe recordar que la hipótesis nula del test estadístico aplicado es que los datos pertenecen a una distribución normal, es decir analiza si el exceso de la asimetría y curtosis que presente la distribución sea cero al compararlo con los parámetros de una distribución normal.
Por lo tanto, se puede concluir que estadísticamente el Modelo 1 tiene residuos que no siguen una distribución normal, lo mismo se conluye con los residuos del Modelo 3, los cuales si siguen una distribución normal.
En suma, se puede decir que el mejor modelo que se puede utilizar es el Modelo 3 ya que presenta niveles bajos autocorrelacion serial y además los residuos siguen una distribución normal.
PREDICCIÓN DEL MODELO
Para poder predecir con el modelo empleado, inicialmente se extrajo 10 datos para poder analizar la capacidad predictiva, ejecutamos el siguiente código:
M1_F = forecast(m1, h = 4, level = c(30,60,90))
M3_F = forecast(m3, h = 4, level = c(30,60,90))
par(mfrow=c(1,2))
plot(M1_F, main = "Forecast arima(2,1,0)(0,1,0)")
plot(M3_F, main = "Forecast arima(0,1,1)(0,1,0)")
Un criterio de elección de acuerdo a la capacidad de predicción es el del Error Cuadrático Medio (MSE) y el modelo que tenga el menor MSE será el seleccionado y se proyectará, se ejecutara el siguiente código:
pre.ts.2 = pre.ts[38:47]
pre.ts.2 <-ts(pre.ts.2, start = c(2019,3),frequency=12)
M1_F$mean<-ts(M1_F$mean,start = c(2019,3),frequency=12)
M3_F$mean<-ts(M3_F$mean,start = c(2019,3),frequency=12)
#MCE
MSE1.F = MSE(M1_F$mean, pre.ts.2)
MSE3.F = MSE(M3_F$mean, pre.ts.2)
MSE.F = cbind(MSE1.F, MSE3.F )
MSE.F## MSE1.F MSE3.F
## [1,] 0.009011287 0.02208866Como se puede observar el modelo 3 es el que tiene menor MSE, por lo tanto, es el seleccionado.
Para hacer proyecciones finales se necesita toda la data, por ello se agregará de nuevo los datos extraídos anteriormente con el fin predictivo y se graficará la predicción, ejecutaremos el siguiente código:
pre.ts.s = pre.ts[1:47]
pre.ts.f = forecast(m3, h = 4, level = c(30,60,90))
plot(pre.ts.f, main = "Forecast arima(3,1,1)(0,1,0)")
Se pueden observar con mas detalles los valores con el siguiente código:
##
## Forecast method: ARIMA(1,1,1)(0,1,0)[12]
##
## Model Information:
##
## Call:
## arima(x = pre.ts1, order = c(1, 1, 1), seasonal = list(order = c(0, 1, 0), period = 12))
##
## Coefficients:
## ar1 ma1
## -0.3464 -0.8369
## s.e. 0.1999 0.1471
##
## sigma^2 estimated as 0.01688: log likelihood = 14.62, aic = -23.25
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.01149771 0.1053942 0.05992389 5.83656 112.3321 0.4303167
## ACF1
## Training set -0.09239702
##
## Forecasts:
## Point Forecast Lo 30 Hi 30 Lo 60 Hi 60 Lo 90
## 39 0.001946121 -0.048122784 0.05201503 -0.10741494 0.11130718 -0.2117877
## 40 0.076303146 0.025400006 0.12720629 -0.03488006 0.18748635 -0.1409919
## 41 0.061213726 0.009061292 0.11336616 -0.05269820 0.17512565 -0.1614143
## 42 -0.192969178 -0.245293195 -0.14064516 -0.30725588 -0.07868248 -0.4163296
## Hi 90
## 39 0.21567995
## 40 0.29359816
## 41 0.28384172
## 42 0.03039126