segundo parcial

series de tiempo

produccion de miel

```

library(forecast)
library(TSA)
library(urca)
library(ggplot2)

rm(list=ls())
base <- read.csv("C:\Users\Personal\Desktop\produccion de plata.csv", header=T, sep=',')

Error: '\U' used without hex digits in character string starting ""C:\U"

attach(base)

# Declarar una base mensual
ts.miel_o<-ts(base$Toneladas,start = c(2003/01,1),frequency = 12)
ts.miel_d<-ts(base$Toneladas,start = c(2003/01,1),frequency = 12)
ts.miel_t<-ts(base$Toneladas,start = c(2003/01,1),frequency = 12)

########## ETAPA: ESPECIFICACIÃ<U+0093>N  ##########
# Graficar serie en niveles 
autoplot(ts.miel_o, series='Original') +
 autolayer(ts.miel_t, series='Tendencia-ciclo')

ts.miel <- diff(ts.miel_t) # se elimina la tendencia
autoplot(ts.miel)

BoxCox.ar(ts.miel_d) # no es posible encontrar un valor de alfa adecuado

log(ts.miel)

NaNs produced

          Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
2003          6.732211 6.898715 7.978996      NaN      NaN      NaN      NaN
2004      NaN 7.134094 5.645447 7.847763      NaN      NaN      NaN      NaN
2005      NaN 4.234107 6.632002 8.003697 7.044033      NaN      NaN 5.710427
2006      NaN 6.248043 7.935230 5.529429 4.248495      NaN      NaN      NaN
2007      NaN 6.845880 7.881937 7.290293      NaN      NaN      NaN      NaN
2008      NaN 6.410175 8.065894 7.495542      NaN      NaN      NaN      NaN
2009      NaN      NaN 7.480428 7.773174 7.109879      NaN      NaN 4.304065
2010      NaN 6.035481 7.820440 6.839476 6.408529      NaN      NaN 2.890372
2011      NaN 6.864848 8.309185 7.640604      NaN      NaN      NaN      NaN
2012      NaN 7.458763 7.547502 7.393878 6.113682      NaN      NaN      NaN
2013      NaN      NaN 7.215975 8.280964 5.043425      NaN      NaN 3.526361
2014      NaN 7.398174 7.544861 7.820038 6.320768      NaN      NaN      NaN
2015      NaN 7.164720 7.170888 8.680502      NaN      NaN      NaN      NaN
2016      NaN 6.679599 4.779123                                             
          Sep      Oct      Nov      Dec
2003 7.628031 8.178358 6.682109 6.598509
2004 7.311886 7.301822 8.335192      NaN
2005 6.688355 7.733246 8.344743      NaN
2006 7.315218 7.862497 8.308199      NaN
2007 6.289716 7.999343 7.931285      NaN
2008 6.837333 8.130059 8.165932      NaN
2009 6.769642 8.087025 8.126223      NaN
2010 6.561031 8.191740 7.821643      NaN
2011 6.620073 8.266678 7.858254      NaN
2012 7.090910 8.295798 8.187299      NaN
2013 5.257495 8.404920 8.251142      NaN
2014 6.645091 8.235626 8.475329      NaN
2015 6.854355 8.070594 8.348538      NaN
2016                                    

plot(log(ts.miel)

Error: Incomplete expression: plot(log(ts.miel)

########## ANÁLISIS DE LA SERIE DESESTACIONALIZADA ##########
# Prueba de raÃ?z unitaria (H0: serie no estacionaria vs. H1: serie estacionaria)
# Serie en niveles
summary(ur.df(ts.miel_t)) # no es estac.

############################################### 
# 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 
 -5830   -689   1156   2531   6243 

Coefficients:
           Estimate Std. Error t value Pr(>|t|)    
z.lag.1    -0.17017    0.03764  -4.522 1.21e-05 ***
z.diff.lag  0.35251    0.07512   4.693 5.89e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2485 on 155 degrees of freedom
Multiple R-squared:  0.1797,    Adjusted R-squared:  0.1692 
F-statistic: 16.98 on 2 and 155 DF,  p-value: 2.143e-07


Value of test-statistic is: -4.5215 

Critical values for test statistics: 
      1pct  5pct 10pct
tau1 -2.58 -1.95 -1.62

summary(ur.df(ts.miel)) # ya es estac.

############################################### 
# 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 
-5348.5 -1523.3   -85.7  1582.5  6028.0 

Coefficients:
           Estimate Std. Error t value Pr(>|t|)    
z.lag.1    -1.11156    0.08348 -13.316  < 2e-16 ***
z.diff.lag  0.51663    0.06897   7.491 4.98e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2271 on 154 degrees of freedom
Multiple R-squared:  0.5359,    Adjusted R-squared:  0.5299 
F-statistic: 88.91 on 2 and 154 DF,  p-value: < 2.2e-16


Value of test-statistic is: -13.3157 

Critical values for test statistics: 
      1pct  5pct 10pct
tau1 -2.58 -1.95 -1.62

# funciones de autocorrelación
ggAcf(ts.miel)

ggPacf(ts.miel)

eacf(ts.miel)

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

ggtsdisplay(ts.miel)

# Propuesta de modelo con base en BIC (criterio de información de Bayes)
res <- armasubsets(y=ts.miel, nar=12, nma=12, y.name='test', ar.method='ols')

1  linear dependencies found

Reordering variables and trying again:

plot(res)

########## ETAPA: ESTIMACIÃ<U+0093>N  ##########
# Estimación (metodo: Maxima verosimilitud 'ML')
propuesta1 <- arima(ts.miel, order=c(2,0,0),method='ML') # AR(2)

propuesta2 <- arima(ts.miel, order=c(0,0,4),method='ML') # MA(4)

propuesta3 <- arima(ts.miel, order=c(1,0,1),method='ML') # ARMA(1,1)

########## ETAPA: DIAGNÃ<U+0093>STICO  ##########
##### Residuos estandarizados ##### 
# Grafica de los residuos estandarizados
plot(rstandard(propuesta1),ylab ='Residuos estandarizados', type='o'); abline(h=0)

plot(rstandard(propuesta2),ylab ='Residuos estandarizados', type='o'); abline(h=0)

plot(rstandard(propuesta3),ylab ='Residuos estandarizados', type='o'); abline(h=0)

# Prueba Ljung-Box (independencia de residuos)
tsdiag(propuesta1,gof=8,omit.initial=F)

tsdiag(propuesta2,gof=8,omit.initial=F)

tsdiag(propuesta3,gof=8,omit.initial=F)

# revisión de residuales
checkresiduals(propuesta1)

    Ljung-Box test

data:  Residuals from ARIMA(2,0,0) with non-zero mean
Q* = 215.46, df = 21, p-value < 2.2e-16

Model df: 3.   Total lags used: 24

checkresiduals(propuesta2)

    Ljung-Box test

data:  Residuals from ARIMA(0,0,4) with non-zero mean
Q* = 257.84, df = 19, p-value < 2.2e-16

Model df: 5.   Total lags used: 24

checkresiduals(propuesta3)

    Ljung-Box test

data:  Residuals from ARIMA(1,0,1) with non-zero mean
Q* = 434.67, df = 21, p-value < 2.2e-16

Model df: 3.   Total lags used: 24

########## PROPUESTA DE auto.arima()  ##########
propuesta.auto <- auto.arima(ts.miel)
propuesta.auto

Series: ts.miel 
ARIMA(1,0,0)(1,1,0)[12] with drift 

Coefficients:
          ar1     sar1    drift
      -0.3558  -0.3455  -0.9301
s.e.   0.0778   0.0824   4.1859

sigma^2 estimated as 1194056:  log likelihood=-1227.96
AIC=2463.92   AICc=2464.2   BIC=2475.85

#Utilizar los residuos de auto.arima
plot(rstandard(propuesta.auto),ylab ='Residuos estandarizados', type='o'); abline(h=0)

qqnorm(residuals(propuesta.auto)); qqline(residuals(propuesta.auto))

ggAcf(residuals(propuesta.auto))

tsdiag(propuesta.auto, gof=12, omit.initial=F)

checkresiduals(propuesta.auto)

    Ljung-Box test

data:  Residuals from ARIMA(1,0,0)(1,1,0)[12] with drift
Q* = 42.421, df = 21, p-value = 0.003723

Model df: 3.   Total lags used: 24

##### Pronostico a partir de auto.arima (2 años) #####
pronostico <- plot(forecast(propuesta.auto,h=24))

pronostico

$`mean`
             Jan         Feb         Mar         Apr         May         Jun
2016                                      5186.07416   -40.34241 -5496.97838
2017 -7125.08098   952.71867   512.37384  5413.23808   -64.57107 -5821.58302
2018 -7059.56576   883.55288   361.44051                                    
             Jul         Aug         Sep         Oct         Nov         Dec
2016 -3161.46060  -711.92426   868.35262  3383.29766  4405.91838   -58.51641
2017 -3070.93734  -759.79156   880.85446  3304.60297  4328.04575   -55.73418
2018                                                                        

$lower
                80%        95%
Apr 2016  3785.6861  3044.3655
May 2016 -1526.7495 -2313.6059
Jun 2016 -6993.9250 -7786.3607
Jul 2016 -4659.7365 -5452.8758
Aug 2016 -2210.3684 -3003.5968
Sep 2016  -630.1128 -1423.3525
Oct 2016  1884.8295  1091.5884
Nov 2016  2907.4499  2114.2086
Dec 2016 -1556.9849 -2350.2262
Jan 2017 -8623.5495 -9416.7908
Feb 2017  -545.7498 -1338.9912
Mar 2017  -986.0947 -1779.3360
Apr 2017  3656.6917  2726.8322
May 2017 -1851.1390 -2796.8909
Jun 2017 -7611.9164 -8559.6617
Jul 2017 -4861.7470 -5809.7444
Aug 2017 -2550.6615 -3498.6908
Sep 2017  -910.0231 -1858.0564
Oct 2017  1513.7244   565.6906
Nov 2017  2537.1671  1589.1332
Dec 2017 -1846.6129 -2794.6468
Jan 2018 -8850.4444 -9798.4784
Feb 2018  -907.3258 -1855.3597
Mar 2018 -1429.4382 -2377.4721

$upper
                80%        95%
Apr 2016  6586.4622  7327.7828
May 2016  1446.0647  2232.9211
Jun 2016 -4000.0317 -3207.5961
Jul 2016 -1663.1847  -870.0454
Aug 2016   786.5199  1579.7483
Sep 2016  2366.8180  3160.0577
Oct 2016  4881.7658  5675.0069
Nov 2016  5904.3868  6697.6281
Dec 2016  1439.9521  2233.1934
Jan 2017 -5626.6125 -4833.3711
Feb 2017  2451.1872  3244.4285
Mar 2017  2010.8424  2804.0837
Apr 2017  7169.7845  8099.6439
May 2017  1721.9968  2667.7487
Jun 2017 -4031.2496 -3083.5043
Jul 2017 -1280.1277  -332.1303
Aug 2017  1031.0784  1979.1077
Sep 2017  2671.7320  3619.7654
Oct 2017  5095.4815  6043.5154
Nov 2017  6118.9244  7066.9583
Dec 2017  1735.1445  2683.1784
Jan 2018 -5268.6871 -4320.6532
Feb 2018  2674.4316  3622.4655
Mar 2018  2152.3192  3100.3531

pronostico <- plot(forecast(propuesta.auto,h=24))

pronostico

$`mean`
             Jan         Feb         Mar         Apr         May         Jun
2016                                      5186.07416   -40.34241 -5496.97838
2017 -7125.08098   952.71867   512.37384  5413.23808   -64.57107 -5821.58302
2018 -7059.56576   883.55288   361.44051                                    
             Jul         Aug         Sep         Oct         Nov         Dec
2016 -3161.46060  -711.92426   868.35262  3383.29766  4405.91838   -58.51641
2017 -3070.93734  -759.79156   880.85446  3304.60297  4328.04575   -55.73418
2018                                                                        

$lower
                80%        95%
Apr 2016  3785.6861  3044.3655
May 2016 -1526.7495 -2313.6059
Jun 2016 -6993.9250 -7786.3607
Jul 2016 -4659.7365 -5452.8758
Aug 2016 -2210.3684 -3003.5968
Sep 2016  -630.1128 -1423.3525
Oct 2016  1884.8295  1091.5884
Nov 2016  2907.4499  2114.2086
Dec 2016 -1556.9849 -2350.2262
Jan 2017 -8623.5495 -9416.7908
Feb 2017  -545.7498 -1338.9912
Mar 2017  -986.0947 -1779.3360
Apr 2017  3656.6917  2726.8322
May 2017 -1851.1390 -2796.8909
Jun 2017 -7611.9164 -8559.6617
Jul 2017 -4861.7470 -5809.7444
Aug 2017 -2550.6615 -3498.6908
Sep 2017  -910.0231 -1858.0564
Oct 2017  1513.7244   565.6906
Nov 2017  2537.1671  1589.1332
Dec 2017 -1846.6129 -2794.6468
Jan 2018 -8850.4444 -9798.4784
Feb 2018  -907.3258 -1855.3597
Mar 2018 -1429.4382 -2377.4721

$upper
                80%        95%
Apr 2016  6586.4622  7327.7828
May 2016  1446.0647  2232.9211
Jun 2016 -4000.0317 -3207.5961
Jul 2016 -1663.1847  -870.0454
Aug 2016   786.5199  1579.7483
Sep 2016  2366.8180  3160.0577
Oct 2016  4881.7658  5675.0069
Nov 2016  5904.3868  6697.6281
Dec 2016  1439.9521  2233.1934
Jan 2017 -5626.6125 -4833.3711
Feb 2017  2451.1872  3244.4285
Mar 2017  2010.8424  2804.0837
Apr 2017  7169.7845  8099.6439
May 2017  1721.9968  2667.7487
Jun 2017 -4031.2496 -3083.5043
Jul 2017 -1280.1277  -332.1303
Aug 2017  1031.0784  1979.1077
Sep 2017  2671.7320  3619.7654
Oct 2017  5095.4815  6043.5154
Nov 2017  6118.9244  7066.9583
Dec 2017  1735.1445  2683.1784
Jan 2018 -5268.6871 -4320.6532
Feb 2018  2674.4316  3622.4655
Mar 2018  2152.3192  3100.3531

##### Pronostico a partir de auto.arima (2 años) #####
pronostico <- plot(forecast(propuesta.auto,h=24))
pronostico

pronostico <- plot(forecast(propuesta.auto,h=24))
pronostico

plot(cars)

---

plot(cars)

---