Indicaciones: - Analizar la serie de tiempo del paquete forecast: wineind: Australian total wine sales

En todos los casos responder a las preguntas:

  1. Ajustar el mejor modelo SARIMA(p, d, q) × (P, D, Q)s segun BIC. Recuerde que ademas de ajustar el mejor modelo debe probar que los residuales sean ruido blanco gaussiano #preferiblemente.

Escribir la ecuacion que describe el proceso generador de datos basado en esta familia de modelos.

  1. Proyectar h = 6 pasos adelante para la serie wineind y para las precipitaciones

ANALISIS DE LA SERIE DE TIEMPO DE WINEIND (Australian total wine sales).

#install.packages("forecats")
library(urca)
library(tseries)
library(ggplot2)
library(ggfortify)
library(dplyr)
library(plotly)
library(TSstudio)
library(readxl)
library(highcharter)
library(astsa)
#library(FitAR)
library(forecast)
data(wineind)

Los datos observados corresponden a las Ventas totales de vino australiano por parte de los productores de vino en botellas < = 1 litro entre Enero de 1980 a Agosto de 1994.

Aspectos descriptivos de la serie de tiempo wineind:

start(wineind) #Conocer el tiempo de comienzo o inicio de la serie
## [1] 1980    1
end(wineind) #Conocer el tiempo de fin de la serie
## [1] 1994    8
frequency(wineind) #Conocer la frecuencia de la serie
## [1] 12

La sere inicia en enero de 1980 y termina en diciembre de 1994.

Unidad de tiempo a la que pertenece cada observación de la serie:

cycle(wineind)
##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1980   1   2   3   4   5   6   7   8   9  10  11  12
## 1981   1   2   3   4   5   6   7   8   9  10  11  12
## 1982   1   2   3   4   5   6   7   8   9  10  11  12
## 1983   1   2   3   4   5   6   7   8   9  10  11  12
## 1984   1   2   3   4   5   6   7   8   9  10  11  12
## 1985   1   2   3   4   5   6   7   8   9  10  11  12
## 1986   1   2   3   4   5   6   7   8   9  10  11  12
## 1987   1   2   3   4   5   6   7   8   9  10  11  12
## 1988   1   2   3   4   5   6   7   8   9  10  11  12
## 1989   1   2   3   4   5   6   7   8   9  10  11  12
## 1990   1   2   3   4   5   6   7   8   9  10  11  12
## 1991   1   2   3   4   5   6   7   8   9  10  11  12
## 1992   1   2   3   4   5   6   7   8   9  10  11  12
## 1993   1   2   3   4   5   6   7   8   9  10  11  12
## 1994   1   2   3   4   5   6   7   8

La unidad de tiempo asignada es numerica por cada uno de los meses del año, es decir, mes 1 corresponde a enero y asi sucesivamente

Representación gráfica de la serie:

library(tseries)
library(ggplot2)
library(ggfortify)
library(dplyr)
library(plotly)
library(TSstudio)

Grafico sencillo:

g<-autoplot(wineind,ts.colour = "red")+
  ggtitle("Ventas totales de vino australiano")+
  xlab("Año")+
  ylab("Botellas de vino") 
g

Ahora con un grafico dinamico:

ggplotly(g)
summary(wineind)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   13652   22115   24669   25392   28461   40226

Con el grafico dinamico puedo ver el año y mes exacto con mayor y menor cantidad de ventas de vino, tal es el caso de diciembre de 1987 con 40226 botellas de vino vendidas frente al año 1994 con menor ventas de botellas de vino de 13652.

Descomposición aditiva de la serie de tiempo wineindYt= Tt+ St+ Et. En sus componentes de tendencia, estacional y componente aleatorio.

library(highcharter)
hchart(stl(wineind, s.window='periodic'))

En los graficos de la descomposicion aditiva de la seria se observan 4 graficas dentro de una misma.

El primer grafico presenta los datos reales, en el eje x para los años (1980 a 1994) y para el eje Y la cantidad de botellas de vino vendidas.

El segundo grafico presenta el componente estacional (seasonal), se observan ciclos, no se evidencia una variacion especifica para un cierto periodo.

El tercer grafico presenta el componente de tendencia (trend) mostrando una tendencia positiva, de crecimiento.

El cuarto y ultimo grafico representa el componente aleatorio, los residuos (remainder), es decir, lo que no absorbe el modelo y se espera se encuentren lo mas cercano a cero, sin embargo se observan punto alejados de este valor.

Boxplot por mes

boxplot(wineind ~ cycle(wineind))

En el diagrama de cajas se observa la cantidad de botellas de vino vendidas en el eje y frente a los meses de enero a diciembre del periodo comprendido entre 1980 a 1994 en el eje x con el fin de identificar una tendencia dentro de los meses con mayores o menores ventas segun sea el interes, para el caso de estudio se observa que en los meses de agosto y diciembre se da la mayor cantidad de ventas

Gráficos de estacionalidad

ts_seasonal(wineind, type = "all")

Otra visualizacion de la serie de tiempo son los graficos de estacionalidad, que presenta 3 graficas diferentes dentro de la misma.

El primer grafico representa en lineas las ventas de botellas de vino frente a los meses con el fin de observar tendencia de crecimiento o decreecimiento en meses especificos a traves de los años.

El segundo grafico similar al anterior, tambien representacion en lineas pero ahora las ventas de botellas de vino frente a los años con el fin de observar tendencias, es decir, el año con mayor o menos venta en Australia.

Por ultimo, encontramos un grafico de Box Plots y al lado de cada uno sus datos reales, en el eje y las ventas de botellas devino, en el eje x los meses del año; se observa el mes de dciembre con mayor numero de ventas.

Punto 1. Ajuste al modelo modelo SARIMA(p, d, q) × (P, D, Q)s segun BIC

z1=ts(wineind, start = c(1980,1), frequency = 12)
z1
##        Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
## 1980 15136 16733 20016 17708 18019 19227 22893 23739 21133 22591 26786 29740
## 1981 15028 17977 20008 21354 19498 22125 25817 28779 20960 22254 27392 29945
## 1982 16933 17892 20533 23569 22417 22084 26580 27454 24081 23451 28991 31386
## 1983 16896 20045 23471 21747 25621 23859 25500 30998 24475 23145 29701 34365
## 1984 17556 22077 25702 22214 26886 23191 27831 35406 23195 25110 30009 36242
## 1985 18450 21845 26488 22394 28057 25451 24872 33424 24052 28449 33533 37351
## 1986 19969 21701 26249 24493 24603 26485 30723 34569 26689 26157 32064 38870
## 1987 21337 19419 23166 28286 24570 24001 33151 24878 26804 28967 33311 40226
## 1988 20504 23060 23562 27562 23940 24584 34303 25517 23494 29095 32903 34379
## 1989 16991 21109 23740 25552 21752 20294 29009 25500 24166 26960 31222 38641
## 1990 14672 17543 25453 32683 22449 22316 27595 25451 25421 25288 32568 35110
## 1991 16052 22146 21198 19543 22084 23816 29961 26773 26635 26972 30207 38687
## 1992 16974 21697 24179 23757 25013 24019 30345 24488 25156 25650 30923 37240
## 1993 17466 19463 24352 26805 25236 24735 29356 31234 22724 28496 32857 37198
## 1994 13652 22784 23565 26323 23779 27549 29660 23356
par(mfrow=c(3,3))
plot(z1)
acf(z1, lag.max = 40)
pacf(z1, lag.max = 40)
plot(diff(z1), lag.max=40)
acf(diff(z1), lag.max = 40)
pacf(diff(z1), lag.max = 40)
plot(diff(z1,12), lag.max=40)
acf(diff(z1,12), lag.max = 40)
pacf(diff(z1, 12), lag.max = 40)

Observamos en el primer grafico los datos de la serie, en la segunda grafica hacia abajo se observan los datos con una diferencia ordinaria (d1) y al final, grafico con una diferencia estacional (D1).

Pruebas de dickey Fuller Para z1

adf.test(z1)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  z1
## Dickey-Fuller = -6.9449, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

Prueba de estacionariedad para los datos de la serie.

Para diff(z1) d=1

adf.test(diff(z1))
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(z1)
## Dickey-Fuller = -8.9901, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

Prueba de estacionariedad para los datos con una diferencia ordinaria.

df para diff(z1,12) D=1

adf.test(diff(z1,12))
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(z1, 12)
## Dickey-Fuller = -4.458, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

Prueba de estacionariedad para los datos con una diferencia Estacional.

ts_cor(diff(z1,12), lag.max = 50)

Con ayuda de los graficos ACF y PACF postulo los diferentes modelos para la serie.

PROBANDO DISTINTOS MODELOS SARIMA

MODELO_0: sugerido R

La función auto arima sugiere un SARIMA(1,1,2)(0,1,1)[12].

auto.arima(z1)
## Series: z1 
## ARIMA(1,1,2)(0,1,1)[12] 
## 
## Coefficients:
##          ar1      ma1     ma2     sma1
##       0.4299  -1.4673  0.5339  -0.6600
## s.e.  0.2984   0.2658  0.2340   0.0799
## 
## sigma^2 = 5399312:  log likelihood = -1497.05
## AIC=3004.1   AICc=3004.48   BIC=3019.57

MODELO_1: SARIMA(1,1,1)x(0,1,1)_12

modelo1<-stats::arima(z1,
                      order=c(1,1,1), 
                      seasonal=list(order=c(0,1,1),
                                    period=12), fixed=c(NA,NA,NA))

Se plantea el modelo, en este caso con fixed:3 y se verifica que coeficientes son significativos, es decir, menores a 0.05:

modelo1
## 
## Call:
## stats::arima(x = z1, order = c(1, 1, 1), seasonal = list(order = c(0, 1, 1), 
##     period = 12), fixed = c(NA, NA, NA))
## 
## Coefficients:
##           ar1      ma1     sma1
##       -0.1078  -0.8850  -0.6435
## s.e.   0.0853   0.0346   0.0831
## 
## sigma^2 estimated as 5363686:  log likelihood = -1498.37,  aic = 3004.75
tt <- modelo1$coef[which(modelo1$coef!=0)]/sqrt(diag(modelo1$var.coef))
1 - pt(abs(tt),(modelo1$nobs - length(modelo1$coef[which(modelo1$coef!=0)])))
##          ar1          ma1         sma1 
## 1.040881e-01 0.000000e+00 5.069278e-13

Despues de ver que coeficientes son significativos, se ajusta de nuevo el modelo planteado y tomamos el valor del BIC:

modelo1<-stats::arima(z1,
                      order=c(1,1,1), 
                      seasonal=list(order=c(0,1,1),
                                    period=12), fixed=c(0,NA,NA)) 
modelo1
## 
## Call:
## stats::arima(x = z1, order = c(1, 1, 1), seasonal = list(order = c(0, 1, 1), 
##     period = 12), fixed = c(0, NA, NA))
## 
## Coefficients:
##       ar1      ma1     sma1
##         0  -0.9009  -0.6242
## s.e.    0   0.0291   0.0864
## 
## sigma^2 estimated as 5431462:  log likelihood = -1499.16,  aic = 3004.32
BIC(modelo1)
## [1] 3013.599

DIAGNOSTICO: Pruebas sobre los residuales del modelo1:

et<-residuals(modelo1)
x1.fit <- fitted(modelo1)

par(mfrow=c(3,2))
plot(z1,type="l",lty=2)
#lines(x1.fit,type="l",col="red")
plot(scale(et),type="l",main="Residuales")
#abline(h=2*sqrt(var(scale(et))),col="red",lty=2)
#abline(h=-2*sqrt(var(scale(et))),col="red",lty=2)
acf(et)
pacf(et)
#qqPlot(scale(et))
acf(abs(et)) #Mide Estructura Heterocedastica

Test de Autocorrelacion de Ljung-Box

Ho: r1=r2=r3=…=rlag=0

Ha: Al menos una es dif de cero

tsdiag(modelo1, gof.lag = 30)

length(z1)
## [1] 176
log(length(z1))
## [1] 5.170484
Box.test(et,lag=7,type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  et
## X-squared = 8.4209, df = 7, p-value = 0.2969

Prueba normalidad: Test de Normalidad basado en Sesgo y Curtosis

Ho: Datos provienen de una Dist. Normal

Ha: Los datos no provienen de una Dist. Normal

jarque.bera.test(et)
## 
##  Jarque Bera Test
## 
## data:  et
## X-squared = 63.7, df = 2, p-value = 1.477e-14

Test de Aleatoriedad

Ho: Residuales exhiben un comport. de aleatoriedad

Ha: Residuales no exhiben estructura (Tendencia, o cualquier otro comportamiento predecible)

runs.test(as.factor(sign(et)), alternative="two.sided")
## 
##  Runs Test
## 
## data:  as.factor(sign(et))
## Standard Normal = -0.75432, p-value = 0.4507
## alternative hypothesis: two.sided

Para el modelo1 de la serie wineind cumplen las pruebas de Autocorrelacion y Aleatoridad para los residuos. Los residuos nos son normales, sin embargo por esta prueba no se descarta el modelo.

MODELO_2: SARIMA(1,1,1)X(1,1,1)_12

modelo2<-stats::arima(z1,
                      order=c(1,1,1), 
                      seasonal=list(order=c(1,1,1),
                                    period=12), fixed=c(NA,NA,NA,NA)) 
modelo2
## 
## Call:
## stats::arima(x = z1, order = c(1, 1, 1), seasonal = list(order = c(1, 1, 1), 
##     period = 12), fixed = c(NA, NA, NA, NA))
## 
## Coefficients:
##           ar1      ma1    sar1     sma1
##       -0.0831  -0.8888  0.1496  -0.7461
## s.e.   0.0890   0.0347  0.1333   0.1055
## 
## sigma^2 estimated as 5299571:  log likelihood = -1497.78,  aic = 3005.55
tt <- modelo2$coef[which(modelo2$coef!=0)]/sqrt(diag(modelo2$var.coef))
1 - pt(abs(tt),(modelo2$nobs - length(modelo2$coef[which(modelo2$coef!=0)])))
##          ar1          ma1         sar1         sma1 
## 1.758919e-01 0.000000e+00 1.317847e-01 2.275091e-11

Despues de ver que coeficientes son significativos:

modelo2<-stats::arima(z1,
                      order=c(1,1,1), 
                      seasonal=list(order=c(1,1,1),
                                    period=12), fixed=c(0,NA,0,NA)) 
BIC(modelo2)
## [1] 3013.599

DIAGNOSTICO: Pruebas sobre los residuales del modelo.

et2<-residuals(modelo2)
x1.fit <- fitted(modelo2)

par(mfrow=c(3,2))
plot(z1,type="l",lty=2)
#lines(x1.fit,type="l",col="red")
plot(scale(et2),type="l",main="Residuales")
abline(h=2*sqrt(var(scale(et2))),col="red",lty=2)
abline(h=-2*sqrt(var(scale(et2))),col="red",lty=2)
acf(et2)
pacf(et2)
#qqPlot(scale(et2))
acf(abs(et2)) #Mide Estructura Heteroced?stica

Test de Autocorrelacion de Ljung-Box

Ho: r1=r2=r3=…=rlag=0

Ha: Al menos una es dif de cero

tsdiag(modelo2, gof.lag = 20)

length(z1)
## [1] 176
log(length(z1))
## [1] 5.170484
Box.test(et2,lag=7,type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  et2
## X-squared = 8.4218, df = 7, p-value = 0.2969

Prueba normalidad: Test de Normalidad basado en Sesgo y Curtosis

Ho: Datos provienen de una Dist. Normal

Ha: Los datos no provienen de una Dist. Normal

jarque.bera.test(et2)
## 
##  Jarque Bera Test
## 
## data:  et2
## X-squared = 63.695, df = 2, p-value = 1.477e-14

Test de Aleatoriedad

Ho: Residuales exhiben un comport. de aleatoriedad

Ha: Residuales no exhiben estructura (Tendencia, o cualquier otro comportamiento predecible)

runs.test(as.factor(sign(et2)), alternative="two.sided")
## 
##  Runs Test
## 
## data:  as.factor(sign(et2))
## Standard Normal = -0.75432, p-value = 0.4507
## alternative hypothesis: two.sided

MODELO_3: SARIMA(0,1,11)X(0,1,1)_12

modelo3<-stats::arima(z1,
                      order=c(0,1,11), 
                      seasonal=list(order=c(0,1,1),
                                    period=12), fixed=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)) 
modelo3
## 
## Call:
## stats::arima(x = z1, order = c(0, 1, 11), seasonal = list(order = c(0, 1, 1), 
##     period = 12), fixed = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA))
## 
## Coefficients:
##           ma1      ma2     ma3      ma4     ma5      ma6     ma7      ma8
##       -1.0106  -0.0920  0.2292  -0.1589  0.1793  -0.1665  0.2289  -0.3790
## s.e.   0.1031   0.1674  0.1725   0.1759  0.1626   0.1155  0.1392   0.1713
##          ma9     ma10    ma11     sma1
##       0.3249  -0.1363  0.1154  -0.7323
## s.e.  0.1567   0.1242  0.0899   0.0751
## 
## sigma^2 estimated as 4727669:  log likelihood = -1491.95,  aic = 3009.9
tt <- modelo3$coef[which(modelo3$coef!=0)]/sqrt(diag(modelo3$var.coef))
1 - pt(abs(tt),(modelo3$nobs - length(modelo3$coef[which(modelo3$coef!=0)])))
##        ma1        ma2        ma3        ma4        ma5        ma6        ma7 
## 0.00000000 0.29161749 0.09293446 0.18391189 0.13596372 0.07579301 0.05107298 
##        ma8        ma9       ma10       ma11       sma1 
## 0.01419385 0.01993537 0.13717243 0.10064577 0.00000000

Despues de ver que coeficientes son significativos:

modelo3<-stats::arima(z1,
                      order=c(0,1,11), 
                      seasonal=list(order=c(0,1,1),
                                    period=12), fixed=c(NA,0,0,0,0,0,NA,NA,NA,0,NA,NA)) 
BIC(modelo3)
## [1] 3025.678

DIAGNOSTICO: Pruebas sobre los residuales del modelo.

et3<-residuals(modelo3)
x1.fit <- fitted(modelo3)

par(mfrow=c(3,2))
plot(z1,type="l",lty=2)
#lines(x1.fit,type="l",col="red")
plot(scale(et3),type="l",main="Residuales")
abline(h=2*sqrt(var(scale(et3))),col="red",lty=2)
abline(h=-2*sqrt(var(scale(et3))),col="red",lty=2)
acf(et3)
pacf(et3)
#qqPlot(scale(et3))
acf(abs(et3)) #Mide Estructura Heteroced?stica

Test de Autocorrelacion de Ljung-Box

Ho: r1=r2=r3=…=rlag=0

Ha: Al menos una es dif de cero

tsdiag(modelo3, gof.lag = 20)

length(z1)
## [1] 176
log(length(z1))
## [1] 5.170484
Box.test(et3,lag=7,type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  et3
## X-squared = 6.256, df = 7, p-value = 0.5102

Prueba normalidad: Test de Normalidad basado en Sesgo y Curtosis

Ho: Datos provienen de una Dist. Normal

Ha: Los datos no provienen de una Dist. Normal

jarque.bera.test(et3)
## 
##  Jarque Bera Test
## 
## data:  et3
## X-squared = 62.188, df = 2, p-value = 3.131e-14

Test de Aleatoriedad

Ho: Residuales exhiben un comport. de aleatoriedad

Ha: Residuales no exhiben estructura (Tendencia, o cualquier otro comportamiento predecible)

runs.test(as.factor(sign(et3)), alternative="two.sided")
## 
##  Runs Test
## 
## data:  as.factor(sign(et3))
## Standard Normal = -0.75594, p-value = 0.4497
## alternative hypothesis: two.sided

El modelo1 es el de mejor comportamiento y sobre este se realizaran los pronosticos.

Punto 2. Proyeccion h = 6 pasos adelante para la serie wineind

Utilizando el modelo con el BIC mas pequeño:

plot(forecast(modelo1,h=6, fan=T))
lines(fitted(modelo1), col="red")

Con base en el modelo1 se realiza el pronostico para la serie wineind obrteniendo un comportamiento de los fdatos proniosticados muy similar a los datos reales.

ANALISIS DE LA SERIE DE TIEMPO PRECIPITACIONES.

library(readxl)
datos <- read_excel("precipitaciones (3).xlsx", sheet = "Hoja1")

precipitaciones<-ts(datos, start=c(2000,1),end= c(2019,12), frequency = 12)

Representación gráfica de la serie:

Grafico dinamico:

ggplotly(autoplot(precipitaciones,ts.colour = "red")+
  ggtitle("Precipitaciones en la cuenca del Rio Negro")+
  xlab("Year")+
  ylab("precipitaciones"))
## Warning: Ignoring unknown parameters: ts.colour
summary(precipitaciones)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0   315.8   570.5   685.9   960.5  2721.0

Con el grafico dinamico se observa el año y mes exacto con mayor y menor cantidad de precipitaciones en la Cuenca del Rio Negro ubicado en Zipaquira, tal es el caso del año 2008 con 0 precipitaciones frente al año 2011 con 2721, la mayor cantidad de precipitaciones.

Descomposición aditiva de la serie de tiempo precipitaciones Yt= Tt+ St+ Et. En sus componentes de tendencia, estacional y componente aleatorio.

library(highcharter)
hchart(stl(precipitaciones, s.window='periodic'))

En los graficos de la descomposicion aditiva de la seria se observan 4 graficas dentro de una misma.

El primer grafico presenta los datos reales, en el eje x para los años (2000 a 2019) y para el eje cantidad de precipitaciones en la Cuenca del Rio Negro ubicado en Zipaquira.

El segundo grafico presenta el componente estacional (seasonal), se observan ciclos, no se evidencia una variacion especifica para un cierto periodo.

El tercer grafico presenta el componente de tendencia (trend) mostrando muy poca tendencia.

El cuarto y ultimo grafico representa el componente aleatorio, los residuos (remainder), es decir, lo que no absorbe el modelo y se espera se encuentren lo mas cercano a cero, se observan levemente alejados de este valor alejados de este valor.

Boxplot por mes

boxplot(precipitaciones ~ cycle(precipitaciones))

En el diagrama de cajas se observa la de cantidad de precipitaciones en la Cuenca del Rio Negro ubicado en Zipaquira en el eje y frente a los meses de enero a diciembre del periodo comprendido entre 2000 a 2019 en el eje x, con el fin de identificar una tendencia dentro de los meses con mayores o menores precipitaciones segun sea el interes, para el caso de estudio se observa que el mes de abril es el de mayor cantidad de precipitaciones.

Gráficos de estacionalidad

ts_seasonal(precipitaciones, type = "all")

Otra visualizacion de la serie de tiempo son los graficos de estacionalidad, que presenta 3 graficas diferentes dentro de la misma.

El primer grafico representa en lineas la cantidad de precipitaciones en la cuenca del Rio Negro frente a los meses con el fin de observar tendencia de crecimiento o decreecimiento en meses especificos a traves de los años, tal es el caso de Abril y Octubre con mayos cantidad de precipitaciones.

El segundo grafico similar al anterior, tambien representacion en lineas pero ahora la cantidad de precipitaciones frente a los años con el fin de observar tendencias, es decir, el año con mayor o menor la cantidad de precipitaciones.

Por ultimo, encontramos un grafico de Box Plots y al lado de cada uno sus datos reales, en el eje y la cantidad de precipitaciones, en el eje x los meses del año; se observa el mes de abril con mayor la cantidad de precipitaciones.

Punto 1. Ajuste al modelo modelo SARIMA (p, d, q) × (P, D, Q)s segun BIC

z1= precipitaciones

par(mfrow=c(3,3))
plot(z1)
acf(z1, lag.max = 60)
pacf(z1, lag.max = 60)
plot(diff(z1), lag.max=60)
acf(diff(z1), lag.max = 60)
pacf(diff(z1), lag.max = 60)
plot(diff(z1,12), lag.max=60)
acf(diff(z1,12), lag.max = 60)
pacf(diff(z1, 12), lag.max = 60)

Pruebas de dickey Fuller Para z1

adf.test(z1)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  z1
## Dickey-Fuller = -4.0974, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary

El resultado de la prueba DF es de 0.01 menos a 0.05, lo que indica que la serie es estacionaria.

Para diff(z1) d=1

adf.test(diff(z1))
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(z1)
## Dickey-Fuller = -8.2258, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary

El resultado de la prueba DF con una diferencia ordinaria (d1) es de 0.01 menos a 0.05, lo que indica que la serie es estacionaria.

df para diff(z1,12) D=1

adf.test(diff(z1,12))
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(z1, 12)
## Dickey-Fuller = -5.2199, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary

El resultado de la prueba DF con una diferencia estacionaria (D1) es de 0.01 menos a 0.05, lo que indica que la serie es estacionaria.

ts_cor(diff(z1,12), lag.max = 60)

PROBANDO DISTINTOS MODELOS SARIMA PARA LA SERIE PRECIPITACIONES

MODELO_0: sugerido R

La función auto arima sugiere un SARIMA(1,0,0)x(0,0,2)_12

modelo0<-auto.arima(z1)
modelo0
## Series: z1 
## ARIMA(1,0,0)(0,0,2)[12] with non-zero mean 
## 
## Coefficients:
##          ar1    sma1    sma2      mean
##       0.2880  0.1233  0.1415  686.8798
## s.e.  0.0623  0.0646  0.0633   54.2692
## 
## sigma^2 = 235749:  log likelihood = -1823.35
## AIC=3656.69   AICc=3656.95   BIC=3674.1
tt <- modelo0$coef[which(modelo0$coef!=0)]/sqrt(diag(modelo0$var.coef))
1 - pt(abs(tt),(modelo0$nobs - length(modelo0$coef[which(modelo0$coef!=0)])))
##          ar1         sma1         sma2    intercept 
## 3.135790e-06 2.872914e-02 1.316706e-02 0.000000e+00
BIC(modelo0)
## [1] 3674.098

DIAGNOSTICO: Pruebas sobre los residuales del modelo.

et<-residuals(modelo0)
x1.fit <- fitted(modelo0)

par(mfrow=c(3,2))
plot(z1,type="l",lty=2)
#lines(x1.fit,type="l",col="red")
plot(scale(et),type="l",main="Residuales")
#abline(h=2*sqrt(var(scale(et))),col="red",lty=2)
#abline(h=-2*sqrt(var(scale(et))),col="red",lty=2)
acf(et)
pacf(et)
#qqPlot(scale(et))
acf(abs(et)) #Mide Estructura Heterocedastica

Test de Autocorrelacion de Ljung-Box

Ho: r1=r2=r3=…=rlag=0

Ha: Al menos una es dif de cero

tsdiag(modelo0, gof.lag = 30)

length(z1)
## [1] 240
log(length(z1))
## [1] 5.480639
Box.test(et,lag=7,type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  et
## X-squared = 13.311, df = 7, p-value = 0.06489

Prueba normalidad: Test de Normalidad basado en Sesgo y Curtosis

Ho: Datos provienen de una Dist. Normal

Ha: Los datos no provienen de una Dist. Normal

jarque.bera.test(et)
## 
##  Jarque Bera Test
## 
## data:  et
## X-squared = 53.423, df = 2, p-value = 2.508e-12

Test de Aleatoriedad

Ho: Residuales exhiben un comport. de aleatoriedad

Ha: Residuales no exhiben estructura (Tendencia, o cualquier otro comportamiento predecible)

runs.test(as.factor(sign(et)), alternative="two.sided")
## 
##  Runs Test
## 
## data:  as.factor(sign(et))
## Standard Normal = 0.93009, p-value = 0.3523
## alternative hypothesis: two.sided

En resumen para el Modelo: Ljung-Box p-value = 0.06489 jarque.bera.test p-value = 2.508e-12 Runs Test p-value = 0.3523

MODELO_1: SARIMA(1,1,0)X(0,1,2)_12

modelo1<-stats::arima(z1,
                      order=c(1,1,0), 
                      seasonal=list(order=c(0,1,2),
                                    period=12), fixed=c(NA,NA,NA)) 
modelo1
## 
## Call:
## stats::arima(x = z1, order = c(1, 1, 0), seasonal = list(order = c(0, 1, 2), 
##     period = 12), fixed = c(NA, NA, NA))
## 
## Coefficients:
##           ar1     sma1    sma2
##       -0.4484  -1.0405  0.0405
## s.e.   0.0593   0.0923  0.0686
## 
## sigma^2 estimated as 252489:  log likelihood = -1752.48,  aic = 3512.96
tt <- modelo1$coef[which(modelo1$coef!=0)]/sqrt(diag(modelo1$var.coef))
1 - pt(abs(tt),(modelo1$nobs - length(modelo2$coef[which(modelo1$coef!=0)])))
##          ar1         sma1         sma2 
## 5.264678e-13 0.000000e+00 2.777091e-01

Despues de ver que coeficientes son significativos:

modelo1<-stats::arima(z1,
                      order=c(1,1,0), 
                      seasonal=list(order=c(0,1,2),
                                    period=12), fixed=c(NA,NA,0))  
BIC(modelo1)
## [1] 3521.588

DIAGNOSTICO: Pruebas sobre los residuales del modelo.

et2<-residuals(modelo1)
x1.fit <- fitted(modelo1)

par(mfrow=c(3,2))
plot(z1,type="l",lty=2)
#lines(x1.fit,type="l",col="red")
plot(scale(et2),type="l",main="Residuales")
abline(h=2*sqrt(var(scale(et2))),col="red",lty=2)
abline(h=-2*sqrt(var(scale(et2))),col="red",lty=2)
acf(et2)
pacf(et2)
#qqPlot(scale(et2))
acf(abs(et2)) #Mide Estructura Heteroced?stica

Test de Autocorrelacion de Ljung-Box

Ho: r1=r2=r3=…=rlag=0

Ha: Al menos una es dif de cero

tsdiag(modelo1, gof.lag = 20)

length(z1)
## [1] 240
log(length(z1))
## [1] 5.480639
Box.test(et,lag=7,type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  et
## X-squared = 13.311, df = 7, p-value = 0.06489

Prueba normalidad: Test de Normalidad basado en Sesgo y Curtosis Ho: Datos provienen de una Dist. Normal Ha: Los datos no provienen de una Dist. Normal

jarque.bera.test(et2)
## 
##  Jarque Bera Test
## 
## data:  et2
## X-squared = 12.744, df = 2, p-value = 0.001709

Test de Aleatoriedad Ho: Residuales exhiben un comport. de aleatoriedad Ha: Residuales no exhiben estructura (Tendencia, o cualquier otro comportamiento predecible)

runs.test(as.factor(sign(et2)), alternative="two.sided")
## 
##  Runs Test
## 
## data:  as.factor(sign(et2))
## Standard Normal = 0.66484, p-value = 0.5062
## alternative hypothesis: two.sided

El modelo1 es el elegido por presentar menor numero BIC y con este modelo se realizaran los pronosticos para la serie Precipitaciones.

Punto 2. Proyeccion h = 6 pasos adelante para la serie precipitaciones

Utilizando el modelo con el BIC mas pequeño, en este caso el modelo1:

plot(forecast(modelo1,h=6, fan=T))
lines(fitted(modelo1), col="red")

El comportamiento de los datos pronosticados es muy similar añl de los datos reales.

NOTA

Para la serie Precipitaciones, no se observa una tendencia marcada ni creciente ni decreciente, al analizar los ciclos sucedev algo similar, estos estan suavizados.

El pronostico presentado se realizo usando un modelo ajustado en tendencia (d1) y en ciclos (D1), a pesar de no tener comportamiento pronunciados para estos dos aspectos.

Se sugiere un nuevo modelo, pero ahora transformando los datos a traves de Box Cox, con el fin de mejorar el comportamiento de la varianza de los datos, la cual no es constante, requisito para una serie de tiempo.