Producción de Limones

library(ggfortify)
library(forecast)
library(fpp2)
library(gridExtra)
library(seasonal)
library(TSA)
library(urca)
library(ggplot2)

Breve referencia a los datos, frecuencia, unidad de medida, fuente, etc.

La serie representa la producción de limones con una frecuencia mensual, la unidad de medida es en toneladas y los datos recopilados fueron obtenidos de la página de INEGI del año 2010 al 2016 (última actualización).

http://www.inegi.org.mx/

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

base<-read.csv(file.choose("limones.csv"))
limones<-ts(base$limones,start = c(2010,1),frequency=12)
Preciolimon<-ts(base$Precio,start = c(2010,1),frequency=12)

View(base)

pl<-plot(limones,xlab="periodo",main="Producción de Limones")

Observamos que la serie presenta estacionalidad, y es completamente cíclica. Los comportamientos de la producción de limón son de manera atípica que se debe a diversos factores, el comportamiento en un año de la producción durante los primeros meses se observa en los cultivos ya que se preparan para que posteriormente los arboles den fruto, por lo que durante los primeros meses la producción, es la mas baja del año; otro de los factores del estiaje de la producción, es que el cítrico no se desarrolle con excelencia y los frutos no den jugo aunado con los problemas cada vez mas presentes de los cambios de temperatura, que afectan las cosechas y con ello directamente en el precio del producto y al sector alimentario.

Descomposición Clasica

limondesc<- decompose(limones, type='additive')
autoplot(limondesc)

Se realizó la descomposición clasica para observar los componentes de la serie (estacionalidad,tendencia y ciclo) por separado

Utilizar transformación Box-Cox / logaritmos para estabilizar la varianza.

Presentar grafica.

BoxCox.ar(limones)

En el caso de esta serie no es necesario aplicar logaritmos debido a que la serie muestra unicamente un comportamiento estacional.

Modelar la serie de tiempo elegida como función de una tendencia, tendencia cuadratica y/o medias estcionales.Presenta el Diagnostico del modelo elegido.

Usando la serie elegida como variable dependiente y el tiempo como independiente.

limones.lm<-lm(limones ~ time(limones))
summary(limones.lm)

Call:
lm(formula = limones ~ time(limones))

Residuals:
     Min       1Q   Median       3Q      Max 
-1060953  -612848   -49737   625882  1161898 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)
(Intercept)   -153408302   93041438  -1.649    0.104
time(limones)      76708      46221   1.660    0.101

Residual standard error: 679200 on 70 degrees of freedom
Multiple R-squared:  0.03786,   Adjusted R-squared:  0.02411 
F-statistic: 2.754 on 1 and 70 DF,  p-value: 0.1015

Graficamos

ggplot(limones, aes(y=limones, x=time(limones))) +
  geom_point() +
  geom_smooth(method = "lm")

Observamos que los residuos tienden a alejarse del componente de referencia lo que determina que el modelo no es el indicado, puede deberse este comportamiento de los residuos a que durante los primeros meses del año la producción es casi nula

Evaluación del Modelo

autoplot(limones.lm, which = 1:6, ncol = 2, label.size = 3)

La distancia de cook representa eventos no controlados, tales como la estacionalidad ya que cada año se da una ventana natural de menos producción de limón en los primeros meses del año, la afectación por las lluvias ya que los productores de los pricipales estados reportan lluvias, lo que tiene repercusiones en la floración y caida de los frutos.El retraso en la maduración ha provocado que se postergue la disponibilidad del producto y por ultimo las lluvias tardias generan entre otros efectos, la plaga conocida como antracnosis, la cual reduce temporalmente el rendimiento de la producción.

El grafico quantil-quantil nos da una referencia con respecto a la normalidad de nuestra serie; se observa una alta relación de los quantiles de los datos reales con los calculados. Por ultimo en la primer gráfica se observa que los residuos casi se ajustan al parametro.

Prueba Shapiro-Wilk

shapiro.test(rstudent(limones.lm))

    Shapiro-Wilk normality test

data:  rstudent(limones.lm)
W = 0.93868, p-value = 0.001626

En este caso el valo “p” < 0.5 por lo que podemos rechazar que ña distribuión sea de tipo normal.

Correlograma

autoplot(acf(rstudent(limones.lm),plot=F))

Se observa que las lineas horizontales representan los errores estandar las cuales se ven rebasadas notoriamente por los rezagos del modelo

En caso de estacionalidad de la serie de tiempo, aplicar diferencias estacionales.

Presentar gráfica.

ggtsdisplay(diff(limones,12),main="diferencia estacional")

ggtsdisplay(diff(diff(limones,12),main="Segunda Diferencia Estacional"))

difestacional<-(diff(diff(limones,12)))

Comprobamos que al aplicar una segunda diferencia estacional la serie se comporta de manera mas estable y los residuos permanecen dentro de la banda.

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

Ho: Existe raíz unitaria Ha: No existe raíz unitaria

summary(ur.df(limones,type="trend",selectlags="AIC"))

############################################### 
# 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 
-1708080  -160783    99130   313299   680866 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.612e+05  1.556e+05   2.322 0.023335 *  
z.lag.1     -4.271e-01  1.050e-01  -4.068 0.000129 ***
tt           2.398e+03  3.113e+03   0.770 0.443771    
z.diff.lag   1.998e-01  1.224e-01   1.633 0.107267    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 521300 on 66 degrees of freedom
Multiple R-squared:  0.2012,    Adjusted R-squared:  0.1648 
F-statistic:  5.54 on 3 and 66 DF,  p-value: 0.001875


Value of test-statistic is: -4.0676 5.6098 8.2969 

Critical values for test statistics: 
      1pct  5pct 10pct
tau3 -4.04 -3.45 -3.15
phi2  6.50  4.88  4.16
phi3  8.73  6.49  5.47

summary(ur.df(limones,type="none",selectlags="AIC"))

############################################### 
# 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 
-1914018   167553   274264   323151   674992 

Coefficients:
           Estimate Std. Error t value Pr(>|t|)
z.lag.1    -0.09093    0.05891  -1.543    0.127
z.diff.lag  0.03964    0.12436   0.319    0.751

Residual standard error: 565500 on 68 degrees of freedom
Multiple R-squared:  0.03406,   Adjusted R-squared:  0.005648 
F-statistic: 1.199 on 2 and 68 DF,  p-value: 0.3078


Value of test-statistic is: -1.5434 

Critical values for test statistics: 
     1pct  5pct 10pct
tau1 -2.6 -1.95 -1.61

adf.test(limones)

p-value smaller than printed p-value

    Augmented Dickey-Fuller Test

data:  limones
Dickey-Fuller = -5.3103, Lag order = 4, p-value = 0.01
alternative hypothesis: stationary

adf.test(difestacional)

p-value smaller than printed p-value

    Augmented Dickey-Fuller Test

data:  difestacional
Dickey-Fuller = -4.6259, Lag order = 3, p-value = 0.01
alternative hypothesis: stationary

El valor de P=0.01 en ambos casos, tanto en la serie original como en la diferenciada, por lo que se rechaza la hipotesis nula, por lo tanto no existe raíz unitaria en ambos modelos lo que indica que nuestra serie es completamente estacionaria.

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

acf(difestacional) 

pacf(difestacional)

eacf(difestacional)

AR/MA
  0 1 2 3 4 5 6 7 8 9 10 11 12 13
0 x o o o o o o o o o o  o  o  o 
1 x o o o o o o o o o o  o  o  o 
2 x x o o o o o o o o o  o  o  o 
3 x o o o o o o o o o o  o  o  o 
4 x o o o o o o o o o o  o  o  o 
5 x o o o o o o o o o o  o  o  o 
6 o x o o x o o o o o o  o  x  o 
7 x o o o x o o o o o o  o  x  o 

De tal modo que proponemos tres modelos Arima

Ajuste1<-arima(difestacional,order=c(1,0,0))
checkresiduals(Ajuste1)

    Ljung-Box test

data:  Residuals from ARIMA(1,0,0) with non-zero mean
Q* = 14.568, df = 22, p-value = 0.8801

Model df: 2.   Total lags used: 24

Ajuste2<-arima(difestacional,order=c(0,0,1))
checkresiduals(Ajuste2)

    Ljung-Box test

data:  Residuals from ARIMA(0,0,1) with non-zero mean
Q* = 16.233, df = 22, p-value = 0.8042

Model df: 2.   Total lags used: 24

Ajuste3 <-arima(difestacional,order=c(1,1,1))
checkresiduals(Ajuste3)

    Ljung-Box test

data:  Residuals from ARIMA(1,1,1)
Q* = 13.948, df = 22, p-value = 0.9033

Model df: 2.   Total lags used: 24

Usar la función auto.arima y compare sus resultados con el inciso anterior

auto<-auto.arima(difestacional,stepwise = F,approximation = F)
auto

Series: difestacional 
ARIMA(1,0,1) with zero mean 

Coefficients:
         ar1      ma1
      0.5528  -0.9621
s.e.  0.1257   0.0662

sigma^2 estimated as 4.191e+09:  log likelihood=-737.03
AIC=1480.05   AICc=1480.49   BIC=1486.28

checkresiduals(auto)

    Ljung-Box test

data:  Residuals from ARIMA(1,0,1) with zero mean
Q* = 11.794, df = 22, p-value = 0.9615

Model df: 2.   Total lags used: 24

Se observa un ligero pero significativo movimiento en el estadistico Q. Mientras que los valores de p solo muestran una pegueña variabilidad, pero se observa que existen valores mas pequeños que en el modelo auto.arima aunque mayores en el estadistico Q.

Presenta la ecuación final y describa

Ajuste1<-arima(difestacional,order=c(1,0,0))
Ajuste1

Call:
arima(x = difestacional, order = c(1, 0, 0))

Coefficients:
          ar1  intercept
      -0.2649   2723.892
s.e.   0.1253   7054.358

sigma^2 estimated as 4.664e+09:  log likelihood = -740.51,  aic = 1485.03

checkresiduals(Ajuste1)

    Ljung-Box test

data:  Residuals from ARIMA(1,0,0) with non-zero mean
Q* = 14.568, df = 22, p-value = 0.8801

Model df: 2.   Total lags used: 24

Ajuste2<-arima(difestacional,order=c(0,0,1))
Ajuste2

Call:
arima(x = difestacional, order = c(0, 0, 1))

Coefficients:
          ma1  intercept
      -0.3770   2344.605
s.e.   0.1658   5539.689

sigma^2 estimated as 4.562e+09:  log likelihood = -739.91,  aic = 1483.81

checkresiduals(Ajuste2)

    Ljung-Box test

data:  Residuals from ARIMA(0,0,1) with non-zero mean
Q* = 16.233, df = 22, p-value = 0.8042

Model df: 2.   Total lags used: 24

Presenta la ecuación final y describa

Ajuste4<-arima(difestacional,order=c(1,0,1))
Ajuste4

Call:
arima(x = difestacional, order = c(1, 0, 1))

Coefficients:
         ar1      ma1  intercept
      0.5246  -0.9998  1452.7242
s.e.  0.1148   0.0483   917.9811

sigma^2 estimated as 3.806e+09:  log likelihood = -735.96,  aic = 1477.91

checkresiduals(Ajuste4)

    Ljung-Box test

data:  Residuals from ARIMA(1,0,1) with non-zero mean
Q* = 12.944, df = 21, p-value = 0.9106

Model df: 3.   Total lags used: 24

Al aplicar la función de auto.arima no pronostica el mejor modelo para la serie, se observa que el modelo sugerido por autoarima resulta ser mejor que los propouestos(Ajuste1, Ajuste2, Ajuste3),dado que el ACI es aun mas bajo, lo cual es mejor, por lo tanto ningun resago se sale de las bandas de los errores estandar, asi como los remanentes de los rezagos se encuentran dentro del parámetro, la prueba lJung Box nos indica que no hya correlación serial.Por lo tanto el mejor modelo es el auto arima.

Realizar el Pronóstico ARIMA para dos años

fit4<-Arima(limones,order=c(1,0,1))
fit4

Series: limones 
ARIMA(1,0,1) with non-zero mean 

Coefficients:
         ar1     ma1       mean
      0.5622  0.1907  1008492.8
s.e.  0.1375  0.1456   160935.6

sigma^2 estimated as 2.734e+11:  log likelihood=-1048.98
AIC=2105.96   AICc=2106.55   BIC=2115.06

autoplot(forecast(fit4))

Utilizando métodos de pronósticos simples y/o suavizamiento exponencial generar un pronóstico para dos años. Elegir el método que presente la raíz del error medio al cuadrado más pequeño.

limonW<-window(limones,start=c(2010,1),end=c(2015,12))
autoplot(limonW)

Luego de haber realizado el recorte de nuestra serie, procedemos a realizar el pronostico para los dos años siguientes mediante el uso de los metodos ya anteriormente mencionados.

autoplot(limonW) + 
  forecast::autolayer(meanf(limonW, h=24), PI=FALSE, series="Mean") + 
  forecast::autolayer(naive(limonW, h=24), PI=FALSE, series="Naïve") + 
  forecast::autolayer(snaive(limonW, h=24), PI=FALSE, series="Seasonal naïve") + 
  ggtitle("Pronostico de Producción de Limones") + 
  xlab("Año") + ylab("Producción de Limones") + 
  guides(colour=guide_legend(title="Tipo de Pronostico"))

Dado que la serie de producción de limones es completamente estacional y comparando el pronostico de los dos años, se concluye que el metodo de naive estacional presenta un comportamiento mas logico y razonable para los siguientes años.

Analizar un VAR, necesita incluir los puntos siguientes

a) Buscar otra variable (incluso más variables) y redactar la teoría que se encuentra detrás de la relación entre variables que propone.

El limón es uno de los citricos mas importante para la economía de México ya que se usa en hogares e industrias de alimentos. Recordando la ley de la oferta la cual nos indica que la oferta es directamente proporcional al precio (cuanto mas alto sea el precio del producto mas unidades se ofreceran a la venta) y por el contrario la ley de la demanda indica que la demanda es ivnersamente proporcional al precio (cuanto mas alto sea el precio, menos demandaran los consumidores).

Por lo tanto se reconoce que el precio esta directamente relacionado con la producción, cabe mencionar que existen fatores que intervienen en el precio de los citricos como por ejemplo la estacionalidad ya que cada año se da una ventana natural de menos producción de limón en los primeros meses del año, la afectación por las lluvias ya que los productores de los pricipales estados reportan lluvias, lo que tiene repercusiones en la floración y caida de los frutos.El retraso en la maduración ha provocado que se postergue la disponibilidad del producto y por ultimo las lluvias tardias generan entre otros efectos, la plaga conocida como antracnosis, la cual reduce temporalmente el rendimiento de la producción y como consecuencia existe un aumento en los precios del producto.

b) Evaluar la estacionalidad de la nueva variable.

autoplot(Preciolimon)

logprecio<-(log(Preciolimon))
difprecio<-(diff(log(Preciolimon,12)))
diffprecio<- (diff(diff(log(Preciolimon,12))))
autoplot(difprecio)

autoplot(diffprecio)

Se observa que el indice de precios por si solos tiene tendencia y muestra quiza un pcoo de estacionalidad pero encontramos que la varianza no es estable, por lo tanto se aplico do diferencias para estabilizar la serie.

Proponga el orden de Rezagos optimos para el VAR y presentar la estimación.

library(zoo)

Attaching package: <U+393C><U+3E31>zoo<U+393C><U+3E32>

The following objects are masked from <U+393C><U+3E31>package:base<U+393C><U+3E32>:

    as.Date, as.Date.numeric

library(vars)

Loading required package: MASS

Attaching package: <U+393C><U+3E31>MASS<U+393C><U+3E32>

The following objects are masked from <U+393C><U+3E31>package:fma<U+393C><U+3E32>:

    cement, housing, petrol

Loading required package: strucchange
Loading required package: sandwich
Loading required package: lmtest

library(dynlm)
library(lmtest)

limones2<- cbind.zoo(limones = difestacional, Precio = diffprecio) 
View(limones2)

autoplot(limones2[,1:2],facet=TRUE)

VARselect(limones2, lag.max = 2, type="const")[["selection"]]

Error in VARselect(limones2, lag.max = 2, type = "const") : 
NAs in y.

Evaluar la cointengración usando la metodología de Engle-Granger y Johansen

par(mfrow=c(1,2))
plot(limones, main="Producción de Limones")
plot(Preciolimon,main="Precio General")

ur.df(limones)

############################################################### 
# Augmented Dickey-Fuller Test Unit Root / Cointegration Test # 
############################################################### 

The value of the test statistic is: -1.5434 

ur.df(Preciolimon)

############################################################### 
# Augmented Dickey-Fuller Test Unit Root / Cointegration Test # 
############################################################### 

The value of the test statistic is: 1.1175 

Ho:La serie no es estacionaria

Según la prueba de cointegración, indica que el término del error es mayor a los criticos (y1=-1, y2=1) por lo tanto no son estacionarias

lp.reg <- dynlm(Preciolimon ~ L(Preciolimon) + limones + L(limones))
summary(lp.reg)

Time series regression with "ts" data:
Start = 2010(2), End = 2015(12)

Call:
dynlm(formula = Preciolimon ~ L(Preciolimon) + limones + L(limones))

Residuals:
     Min       1Q   Median       3Q      Max 
-15.2849  -3.2170  -0.4302   2.1711  19.4172 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)     1.599e+01  6.406e+00   2.496  0.01501 *  
L(Preciolimon)  8.754e-01  5.702e-02  15.354  < 2e-16 ***
limones         4.151e-06  1.540e-06   2.695  0.00889 ** 
L(limones)     -5.031e-06  1.567e-06  -3.210  0.00204 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.628 on 67 degrees of freedom
Multiple R-squared:  0.7982,    Adjusted R-squared:  0.7892 
F-statistic: 88.36 on 3 and 67 DF,  p-value: < 2.2e-16

dwtest(lp.reg)

    Durbin-Watson test

data:  lp.reg
DW = 2.7365, p-value = 0.9987
alternative hypothesis: true autocorrelation is greater than 0

Se observa que el valor p resulta ser significativo p>0.5 por lo tanto se rechaza la Ho aceptando la alterna por lo que las variables tienen verdadera autocorrelación entre ellas.

error <- residuals(lp.reg)
plot(error)
abline(h=0)

ur.df(error)

############################################################### 
# Augmented Dickey-Fuller Test Unit Root / Cointegration Test # 
############################################################### 

The value of the test statistic is: -7.862 

Ho:La serie no es estacionaria

Según la prueba de cointegración, indica que el término del error es mayor a los criticos p=-7.86 por lo tanto existe una relacion de cointegración.

limones3 <- ts.intersect(limones,Preciolimon, error)
dyn <- dynlm(d(Preciolimon) ~ L(error) + L(d(limones)) + L(d(Preciolimon)), data=data)
summary(dyn)

Time series regression with "zooreg" data:
Start = mar. 2010, End = dic. 2015

Call:
dynlm(formula = d(Preciolimon) ~ L(error) + L(d(limones)) + L(d(Preciolimon)), 
    data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-19.4016  -3.5103  -0.4213   3.3281  13.5639 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)  
(Intercept)        4.782e-01  8.191e-01   0.584   0.5613  
L(error)          -9.723e-01  4.087e-01  -2.379   0.0203 *
L(d(limones))     -2.955e-06  2.282e-06  -1.295   0.1999  
L(d(Preciolimon))  4.366e-01  3.939e-01   1.108   0.2717  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.499 on 66 degrees of freedom
Multiple R-squared:  0.2484,    Adjusted R-squared:  0.2142 
F-statistic: 7.269 on 3 and 66 DF,  p-value: 0.0002759