library (TSA)
## 
## Attaching package: 'TSA'
## The following objects are masked from 'package:stats':
## 
##     acf, arima
## The following object is masked from 'package:utils':
## 
##     tar
library (forecast)
library(readxl)
library(gridExtra)
library(fpp2)
## Loading required package: ggplot2
## Loading required package: fma
## Loading required package: expsmooth
library(urca)
miel <- read_excel("C:/Users/garpe/Downloads/miel.xls")
View(miel)

1.

Breve descripcion de la serie elegida, así como su frecuencia, unidad de medida, fuente, etc.

Se eligió la producción mensual de miel en México por toneladas en un periodo de doce años (2003-2015), la cual resulta ser una serie estacional y estacionaria porque, como podemos verlo en la siguiente gráfica, la media se mantiene constante.
st<-ts(miel$Dato,frequency=12,start=c(2003,1))
autoplot(st,ylab="Producción en toneladas",xlab="Años",main="Producción mensual de miel en México (2003-2015)")

##### Fuente: SAGARPA. Servicio de Información Agroalimentaria y Pesquera. ###### (http://www.inegi.org.mx/sistemas/bie/)

2.

Graficar los datos, analizar patrones y observaciones atípicas. Apoye su analisis con una descomposición clásica

-La serie refleja un ciclo decreciente en el periodo de 2002 a 2004 y de 2004 a 2008 mantiente un ciclo creciente que se invierte nuevamente de 2008 al 2011 aproximadamente, retomando un ciclo con pendiente positiva de 2013 a 2015.
-Tiene una tendencia positiva
-Refleja mucha estacionalidad
fit<-decompose(st,type = "additive")
plot(fit,col="darkgreen",ylab="eje y",xlab="Años",lwd=.5,type="l",pch=5)

Como podemos observar en la gráfica de estacionalidad, hay dos épocas del año en la que más miel se produce. Durante el mes de abril y finales de octubre/ principios de noviembre
ggseasonplot(st,year.labels = TRUE, year.labels.left = TRUE,ylab="Toneladas",xlab="",main="Producción Mensual de Miel (2013-2015)")

ggseasonplot(st,polar = TRUE,ylab="",xlab="",main="Producción Mensual de Miel (2013-2015)")

3.

Si es necesario, utilizar transformacion Box-Cox/ logaritmos para estabilizar la varianza.

p1<-autoplot(st,lab="Años", ylab="Toneladas")
p2<-autoplot(log(st),xlab="Años",ylab="Toneladas(Log)")
grid.arrange(p1,p2)

De acuerdo al grafico se observa, que la aplicacion de logaritmos no fue suficiente para estabilizar la varianza de la serie. En este caso se tiene que aplicar la diferencia de logaritmos para poder estabilizarla.
plot(st,type="l")

plot(log(st), type="l")

tasa<-na.omit((st-zlag(st))/zlag(st))
plot(diff(log(st)), type = "l")

cor(diff(log(st))[-1],tasa[-1])
## [1] 0.9204683
BoxCox.ar(st)

BoxCox.ar(st,lambda = seq(-5,2,5))

El primer grafico indica el punto maximo que alcanza el 95% de nivel de confianza y el segundo grafico se comprueba al incluir los parametros.
plot(diff(log(st))^2,type="l")

Se elevo la serie al cuadrado segun la prueba Box.Cox para estabilizarla, y podemos observar que en la mayoria de los años esta cercanos a cero, es decir se estabiliza la varianza.

4.

En caso de estacionalidad, aplicar diferencias estacionales.

ggtsdisplay(diff(st),xlab="Años",ylab="Toneladas",main="Diferencias Estacionales")

Al aplicar las diferencias estacionales podemos observar que en la gráfica del proceso autoregresivo, estan marcadaslas correlaciones en cada trimestre.

5.

Usar prueba Dickey-Fuller para evaluar el orden de integracion de la serie. Diferenciar hasta que la serie sea estacionaria.

summary(ur.df(st))
## 
## ############################################### 
## # 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 
##  -5890   -674   1106   2482   6188 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -0.16384    0.03766  -4.351 2.48e-05 ***
## z.diff.lag  0.37085    0.07601   4.879 2.67e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2454 on 152 degrees of freedom
## Multiple R-squared:  0.1823, Adjusted R-squared:  0.1716 
## F-statistic: 16.95 on 2 and 152 DF,  p-value: 2.266e-07
## 
## 
## Value of test-statistic is: -4.3507 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62
Somo el valor del test se concentra en 1pct, sechaza la hipotesis nula, lo cual indica que no existe raiz unitaria en la serie.

6.

Usar herramientas ACF, PACF, EACF y criterios de Akaike/ Bayes para construir propuestas de modelos.

ggAcf(st)

ggPacf(st)

eacf(st)
## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x x x x x x x x x x  x  x  x 
## 1 x x x x x x x x x x x  x  x  x 
## 2 x o o o x o x o o o o  x  o  o 
## 3 o o o o x o x x o o o  x  x  o 
## 4 o x o o x o x o o o o  x  x  x 
## 5 o x o o o o x o o o o  x  x  o 
## 6 o x o x x o x o o o o  x  o  o 
## 7 x o x o x o o o o o o  x  x  o
alm <- armasubsets(y=diff(log(st)),nar=12,nma=12,y.name='test',ar.method='ols')
## Warning in leaps.setup(x, y, wt = wt, nbest = nbest, nvmax = nvmax,
## force.in = force.in, : 1 linear dependencies found
## Reordering variables and trying again:
plot(alm)

ggtsdisplay(st)

7.

Como parte del diagnostico del modelo, analizar los residuos y aplicar la prueba Ljung-Box.

Propuesta 1

fit1<-Arima(st, order=c(2,1,0), seasonal=c(1,1,0))
fit1
## Series: st 
## ARIMA(2,1,0)(1,1,0)[12] 
## 
## Coefficients:
##           ar1      ar2     sar1
##       -0.4843  -0.3442  -0.3224
## s.e.   0.0789   0.0782   0.0845
## 
## sigma^2 estimated as 1055531:  log likelihood=-1193.92
## AIC=2395.84   AICc=2396.13   BIC=2407.69
checkresiduals(fit1)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,0)(1,1,0)[12]
## Q* = 31.135, df = 21, p-value = 0.07143
## 
## Model df: 3.   Total lags used: 24

Propuesta 2

fit2<-Arima(st, order=c(2,0,1), seasonal=c(2,1,0))
fit2
## Series: st 
## ARIMA(2,0,1)(2,1,0)[12] 
## 
## Coefficients:
##          ar1      ar2      ma1     sar1     sar2
##       0.2560  -0.1326  -0.0901  -0.4107  -0.2403
## s.e.  0.6102   0.1185   0.6152   0.0912   0.0942
## 
## sigma^2 estimated as 753225:  log likelihood=-1177.53
## AIC=2367.06   AICc=2367.68   BIC=2384.88
checkresiduals(fit2)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,0,1)(2,1,0)[12]
## Q* = 8.4561, df = 19, p-value = 0.9815
## 
## Model df: 5.   Total lags used: 24

Propuesta 3

fit3<-Arima(st, order=c(2,1,1), seasonal=c(0,1,1))
fit2
## Series: st 
## ARIMA(2,0,1)(2,1,0)[12] 
## 
## Coefficients:
##          ar1      ar2      ma1     sar1     sar2
##       0.2560  -0.1326  -0.0901  -0.4107  -0.2403
## s.e.  0.6102   0.1185   0.6152   0.0912   0.0942
## 
## sigma^2 estimated as 753225:  log likelihood=-1177.53
## AIC=2367.06   AICc=2367.68   BIC=2384.88
checkresiduals(fit3)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,1)(0,1,1)[12]
## Q* = 8.3147, df = 20, p-value = 0.9896
## 
## Model df: 4.   Total lags used: 24

8.

Usar la funcion auto.arima() y compare sus resultados con el inciso anterior.

autoarima<- auto.arima(st, stepwise=FALSE, approximation=FALSE)
autoarima
## Series: st 
## ARIMA(0,0,1)(0,1,1)[12] 
## 
## Coefficients:
##          ma1     sma1
##       0.1681  -0.4677
## s.e.  0.0869   0.0916
## 
## sigma^2 estimated as 741460:  log likelihood=-1177.99
## AIC=2361.99   AICc=2362.16   BIC=2370.9
checkresiduals(autoarima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,0,1)(0,1,1)[12]
## Q* = 8.6851, df = 22, p-value = 0.9948
## 
## Model df: 2.   Total lags used: 24
Ninguna de las propuestas realizadas se acercaron a la propuesta realizada por el modelo autoarima. Se rechaza la hipotesis nula por lo que se concluye que no hay autocorrelación serial.

9.

Presente la ecuacion final y describa.

\[ Y_(t)= \theta (0.1767 ) + \Theta (-0.4718)_{t-12} \]

10.

Crear un pronósitco a dos años

pronostico <- plot(forecast(autoarima,h=24))

pronostico
## $mean
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 2016 2282.9644 3443.1732 5024.6335 9411.8680 9506.2051 4390.9253 1208.9613
## 2017 2263.6895 3443.1732 5024.6335 9411.8680 9506.2051 4390.9253 1208.9613
##            Aug       Sep       Oct       Nov       Dec
## 2016  622.5645 1443.3198 4999.3887 9214.6627 8838.9906
## 2017  622.5645 1443.3198 4999.3887 9214.6627 8838.9906
## 
## $lower
##                 80%        95%
## Jan 2016 1179.44539   595.2777
## Feb 2016 2324.16387  1731.7961
## Mar 2016 3905.62414  3313.2564
## Apr 2016 8292.85869  7700.4909
## May 2016 8387.19578  7794.8280
## Jun 2016 3271.91598  2679.5482
## Jul 2016   89.95198  -502.4158
## Aug 2016 -496.44478 -1088.8125
## Sep 2016  324.31045  -268.0573
## Oct 2016 3880.37937  3288.0116
## Nov 2016 8095.65342  7503.2857
## Dec 2016 7719.98133  7127.6136
## Jan 2017  999.87139   330.8465
## Feb 2017 2175.50157  1504.4368
## Mar 2017 3756.96184  3085.8970
## Apr 2017 8144.19639  7473.1316
## May 2017 8238.53348  7567.4687
## Jun 2017 3123.25368  2452.1889
## Jul 2017  -58.71032  -729.7751
## Aug 2017 -645.10708 -1316.1719
## Sep 2017  175.64815  -495.4167
## Oct 2017 3731.71707  3060.6523
## Nov 2017 7946.99112  7275.9263
## Dec 2017 7571.31903  6900.2542
## 
## $upper
##                80%       95%
## Jan 2016  3386.484  3970.651
## Feb 2016  4562.182  5154.550
## Mar 2016  6143.643  6736.011
## Apr 2016 10530.877 11123.245
## May 2016 10625.214 11217.582
## Jun 2016  5509.935  6102.302
## Jul 2016  2327.971  2920.338
## Aug 2016  1741.574  2333.942
## Sep 2016  2562.329  3154.697
## Oct 2016  6118.398  6710.766
## Nov 2016 10333.672 10926.040
## Dec 2016  9958.000 10550.368
## Jan 2017  3527.508  4196.532
## Feb 2017  4710.845  5381.910
## Mar 2017  6292.305  6963.370
## Apr 2017 10679.540 11350.604
## May 2017 10773.877 11444.942
## Jun 2017  5658.597  6329.662
## Jul 2017  2476.633  3147.698
## Aug 2017  1890.236  2561.301
## Sep 2017  2710.991  3382.056
## Oct 2017  6267.060  6938.125
## Nov 2017 10482.334 11153.399
## Dec 2017 10106.662 10777.727

El pronostico a dos años muestra la serie se mantendra estable, sin ningún ciclo o alguna variación en su tendencia , presentando aún su fuerte estacionalidad por lo que es muy probable (95%) que la producción anual de miel se mantenga constante en el futuro.