Trabajo Final 3er Parcial

library(TSA)

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

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

    acf, arima

The following object is masked from <U+393C><U+3E31>package:utils<U+393C><U+3E32>:

    tar

library(ggplot2)
library(forecast)
library(gridExtra)
library(readxl)
library(urca)
library(vars)

Loading required package: MASS
Loading required package: strucchange
Loading required package: sandwich

library(zoo)
library(fpp2)

Loading required package: fma

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

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

    cement, housing, petrol

Loading required package: expsmooth

library(expsmooth)
library(readxl)
Produccion_mensual_de_miel_en_Mexico <- read_excel("Produccion mensual de miel en Mexico.xls")
View(Produccion_mensual_de_miel_en_Mexico)

1.

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

Se eligio 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.

Fuente: SAGARPA. Servicio de Información Agroalimentaria y Pesquera.

2.

Graficar los datos, analizar patrones y observaciones atipicas. Apoye su analisis con una descomposicion clasica.

-La serie refleja un ciclo decreciente en el periodo de 2003 a 2005 y de 2005 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

Declaracion de la serie y comportamiento general:

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

3.

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

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

Al aplicarle diferencia a los logaritmos se puede observar que el comportamiento de la serie se mantiene estacionario, aunque tiene variaciones pronunciadas desde el 2011 al 2015, pero en general, se observar con claridad el comportamiento estacionario.

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.

4.

Modelar la serie de tiempo elegida como funcion de una tendencia, tendencua cuadratica y/o medias estacionales. Presente el diagnostico del modelo elegido.

Se modela la serie en funcion del tiempo

mes.<-season(st)
modelo1<-lm(st~mes.+time(st)+ I(time(st)^2))
summary(modelo1)

Call:
lm(formula = st ~ mes. + time(st) + I(time(st)^2))

Residuals:
     Min       1Q   Median       3Q      Max 
-2407.75  -565.59   -22.22   512.61  2882.35 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)    2.109e+07  2.243e+07   0.940 0.348653    
mes.February   7.604e+02  3.427e+02   2.219 0.028091 *  
mes.March      2.714e+03  3.427e+02   7.919 6.26e-13 ***
mes.April      5.122e+03  3.428e+02  14.942  < 2e-16 ***
mes.May        4.730e+03  3.428e+02  13.798  < 2e-16 ***
mes.June       1.697e+03  3.428e+02   4.950 2.07e-06 ***
mes.July      -1.228e+03  3.429e+02  -3.583 0.000466 ***
mes.August    -1.586e+03  3.429e+02  -4.624 8.37e-06 ***
mes.September -6.066e+02  3.430e+02  -1.769 0.079064 .  
mes.October    2.660e+03  3.430e+02   7.754 1.57e-12 ***
mes.November   6.075e+03  3.431e+02  17.707  < 2e-16 ***
mes.December   4.608e+03  3.432e+02  13.429  < 2e-16 ***
time(st)      -2.102e+04  2.232e+04  -0.942 0.347912    
I(time(st)^2)  5.240e+00  5.555e+00   0.943 0.347111    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 873.8 on 142 degrees of freedom
Multiple R-squared:  0.9029,    Adjusted R-squared:  0.8941 
F-statistic: 101.6 on 13 and 142 DF,  p-value: < 2.2e-16

La R cuadrada nos demuestra que el modelo se ajusta en un 89% y el intercepto explica que en enero en promedio se producen 2.109e+07 miles de toneladas de miel.

Se puede apreciar que a pesar de las transformaciones a las que se ha sometido la serie esta refleja una tendencia creciente, sin embargo es casi imperceptible.

Grafico de los residuos

Los errores parecen estar distribuidos de manera aleatoria.

Grafico Cuantil Cuantil

No existe normalidad ya que no todos los datos no se encuentran dentro de la línea de normalidad

Histograma

El histograma muestra una ligera concentración a la izquierda

Autocorrelación

Estas dos gráficas muestran el grado de la correlacion que guardan los coeficientes estimados en la regresion. Podemos observar que hay tres momentos en los que los errores se salen de la media 0.

El gráfico de Q-q normal ratifica la conclusión anterior, ya que los valores observados no se situan sobre la recta esperada bajo el supuesto de normalidad.

Prueba Shapiro

shapiro.test(rstandard(modelo1))

    Shapiro-Wilk normality test

data:  rstandard(modelo1)
W = 0.98439, p-value = 0.07605

p-value = 0.07605 lo que esta arriba de 0.05 por lo tanto se acepta la hipotesis nula y se rechaza la alterna, es decir, el componente estocastico esta normalmente distribuido.

5.

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

Ya que la serie presentaba mucha estacionalidad se aplicaron diferencias estacionales. Con esto se estabiliza la varianza, sobre todo de los picos que presentaban los ultimos años de la serie.

6.

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

Ya que la R cuadrada ajustada arroja un valor de 0.1716 podemos rechazar la hipotesis nula y aceptar la alterna (pues el valor es menor a 0.5)

La serie es estacionaria

7.

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

eacf(miel.d)

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 o x o o o o o  x  x  o 
3 x x o o x o x o o o o  x  x  o 
4 x x o x o o x o o o o  x  o  o 
5 x o o x x o x x o o o  x  x  o 
6 x o o o x o x o o o o  x  x  x 
7 x o o o x x o o o o o  x  x  o 

8.

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

Propuesta 1

propuesta1

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(propuesta1)

Propuesta 2

Propuesta2 <-Arima(st, order=c(2,0,1), seasonal=c(2,1,0))
Propuesta2 <-Arima(st, order=c(2,0,1), seasonal=c(2,1,0))
Propuesta2

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(Propuesta2)

Propuesta 3

checkresiduals(Propuesta3)

    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

9.

Usar la funcion autoarima y compare sus resultados con el inciso anterior.

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.

10.

Presente la ecuacion final y describa.

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

11.

Realizar el pronostico ARIMA para dos años.

pronostico

$`mean`
           Jan       Feb       Mar       Apr       May       Jun       Jul       Aug       Sep
2016 2282.9644 3443.1732 5024.6335 9411.8680 9506.2051 4390.9253 1208.9613  622.5645 1443.3198
2017 2263.6895 3443.1732 5024.6335 9411.8680 9506.2051 4390.9253 1208.9613  622.5645 1443.3198
           Oct       Nov       Dec
2016 4999.3887 9214.6627 8838.9906
2017 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.

12.

Utilizandometodos de pronosticos simples y/o suavizamiento exponencial generar un pronostico para dos años. Elegir el metodo que presente la raiz del error medio al cuadrado mas pequeno.

accuracy(mfit3,m4)

'start' value not changedError in window.default(x, ...) : 'start' cannot be after 'end'

14

Analizar un VAR, necesita incluir los puntos siguientes:

a) Buscar otra variable (incluso mas variables) y redactar la teoría que se

encuentra detras de la relación entrevariables que propone.

Variable: Balanza de exportación de Miel Unidad de medida: Miles de dolares Fuente: SAT, SE, BANXICO, INEGI. Balanza Comercial de Mercancías de México. SNIEG. Información de Interés Nacional.

México es el 4to exportador de miel a nivel mundial. Un cambio en los niveles de exportacion provocoarían el aumento consecutivo de la producción de miel en el país para poder cubrir la demanda que se generó.

b) Evaluar la estacionalidad de la nueva variable.

Podemos observar que tanto la serie de produccion y de exportacion de miel en México presentan estacionalidad , para esto se aplico el método de diferencias de las diferencia estacional porque además tienen tendencia. “No hubo necesidad de crear otra base, ya que, las variables tienen la misma longitud”

c) Proponga el orden de rezagos optimo para el VAR y presentar la estimacion.

VARselect(base2, lag.max = 12, type = "const")[["selection"]]

AIC(n)  HQ(n)  SC(n) FPE(n) 
    12      3      2     12 

VARselect(base2, lag.max = 12, type = "const")

$`selection`
AIC(n)  HQ(n)  SC(n) FPE(n) 
    12      3      2     12 

$criteria
                  1            2            3            4            5            6            7
AIC(n) 2.972646e+01 2.955575e+01 2.950091e+01 2.947375e+01 2.946595e+01 2.950698e+01 2.954855e+01
HQ(n)  2.977997e+01 2.964493e+01 2.962577e+01 2.963428e+01 2.966215e+01 2.973886e+01 2.981610e+01
SC(n)  2.985815e+01 2.977523e+01 2.980818e+01 2.986882e+01 2.994881e+01 3.007763e+01 3.020699e+01
FPE(n) 8.129159e+12 6.853788e+12 6.488903e+12 6.316487e+12 6.269644e+12 6.535608e+12 6.817847e+12
                  8            9           10           11           12
AIC(n) 2.955416e+01 2.947774e+01 2.952014e+01 2.955391e+01 2.940972e+01
HQ(n)  2.985739e+01 2.981664e+01 2.989472e+01 2.996416e+01 2.985565e+01
SC(n)  3.030039e+01 3.031177e+01 3.044196e+01 3.056353e+01 3.050712e+01
FPE(n) 6.862558e+12 6.365170e+12 6.650552e+12 6.891205e+12 5.978531e+12

El número de rezagos maximo es 12 asi que se crean las siguientes propuestas:

VAR 1

serial.test(var1, lags.pt=1, type="PT.asymptotic")

    Portmanteau Test (asymptotic)

data:  Residuals of VAR object var1
Chi-squared = 3.1332, df = 0, p-value < 2.2e-16

VAR 2

serial.test(var2, lags.pt=10, type="PT.asymptotic")

    Portmanteau Test (asymptotic)

data:  Residuals of VAR object var2
Chi-squared = 45.333, df = 32, p-value = 0.05939

VAR 3

serial.test(var3, lags.pt=10, type="PT.asymptotic")

    Portmanteau Test (asymptotic)

data:  Residuals of VAR object var3
Chi-squared = 36.643, df = 28, p-value = 0.1269

Var 12

serial.test(var12, lags.pt=10, type="PT.asymptotic")

NaNs producedNaNs produced

    Portmanteau Test (asymptotic)

data:  Residuals of VAR object var12
Chi-squared = 9.8881, df = -8, p-value = NA

Se tomará el Var 3 y se analizaran los rezagos

d) Probar la correlacion serial y normalidad en los residuos.

plot(var3.serial, names = "Exportacion")

Invalid residual name(s) supplied, using residuals of first variable.

No se puede rechazar la hipotesis nula de que no hay correlacion serial

e) Resultados de la prueba de causalidad de Granger

causality(var3, cause='P')$Granger

    Granger causality H0: P do not Granger-cause T

data:  VAR object var3
F-Test = 1.3371, df1 = 3, df2 = 266, p-value = 0.2627

Se rechaza la hipotesis nula de no Granger causalidad

causality(var3, cause='T')$Granger

    Granger causality H0: T do not Granger-cause P

data:  VAR object var3
F-Test = 0.72403, df1 = 3, df2 = 266, p-value = 0.5384

Claramente no es posible rechazar la hipótesis nula de no Granger causalidad.

f) Crear un pronostico para dos años usando el VAR especificado.

predictions <- predict(var3, n.ahead = 10, ci = 0.95)
class(predictions)

[1] "varprd"

plot(predictions, names = "Exportaciones")

Invalid variable name(s) supplied, using first variable.

g) Interpretar las funciones impulso-respuesta.

En el caso de las produccion, se muestra como la cantidad de exportaciones responde ante un impulso de la cantidad de produccion , por lo tanto, las toneladas producidas infieren en la cantidad de exportaiones.

En el caso de la produccion, no nos muestra un resultado similar, entonces, un impulso en las exportaciones no corresponde a una respuesta de la cantidad de las toneladas producidas.

15.

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

Prueba

library(tseries)
po.test(base2[,1:2])

p-value smaller than printed p-value

    Phillips-Ouliaris Cointegration Test

data:  base2[, 1:2]
Phillips-Ouliaris demeaned = -180.27, Truncation lag parameter = 1, p-value = 0.01

Con H0 : las variables no están cointegradas.

Se rechaza la hipotesis nula es decir las variables están cointegradas