Producción de Zinc (TMR)
Hernan Perci Nuñez Palomino
Última edición 11 Diciembre 2021
Total nacional de producción de Zinc en toneladas métricas de contenido recuperado (TMR). Fuente secundaria del INEI.
1 Librerias
library(readxl) # excel
library(knitr) # tablas
library(tidyverse) # varias funciones## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.6 v dplyr 1.0.7
## v tidyr 1.1.4 v stringr 1.4.0
## v readr 2.1.0 v forcats 0.5.1## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(ggthemes) # temas de graficos
library(tseries) ## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(forecast)
library(urca)
library(uroot) # prueba hegy
library(Metrics) # metricas rmse##
## Attaching package: 'Metrics'## The following object is masked from 'package:forecast':
##
## accuracylibrary(lmtest) #para los coeficientes## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric2 Cargar los datos
Zinc <- read_excel("Zinc.xlsx", sheet = "Data")
kable(Zinc,
col.names = c("Fecha", "Zinc (TMR)"),
align = "cc")| Fecha | Zinc (TMR) |
|---|---|
| 2007-01-01 | 98471.46 |
| 2007-02-01 | 96350.45 |
| 2007-03-01 | 107589.89 |
| 2007-04-01 | 93117.68 |
| 2007-05-01 | 109090.20 |
| 2007-06-01 | 118277.49 |
| 2007-07-01 | 116272.12 |
| 2007-08-01 | 101161.38 |
| 2007-09-01 | 98626.96 |
| 2007-10-01 | 107174.43 |
| 2007-11-01 | 87056.73 |
| 2007-12-01 | 98928.10 |
| 2008-01-01 | 110691.37 |
| 2008-02-01 | 102079.42 |
| 2008-03-01 | 112278.81 |
| 2008-04-01 | 106553.66 |
| 2008-05-01 | 114257.44 |
| 2008-06-01 | 121328.51 |
| 2008-07-01 | 114810.94 |
| 2008-08-01 | 126973.22 |
| 2008-09-01 | 110370.43 |
| 2008-10-01 | 118702.77 |
| 2008-11-01 | 105057.49 |
| 2008-12-01 | 123996.17 |
| 2009-01-01 | 116661.85 |
| 2009-02-01 | 97106.22 |
| 2009-03-01 | 101387.08 |
| 2009-04-01 | 104535.89 |
| 2009-05-01 | 104492.31 |
| 2009-06-01 | 102100.39 |
| 2009-07-01 | 108021.25 |
| 2009-08-01 | 115596.97 |
| 2009-09-01 | 93810.62 |
| 2009-10-01 | 118446.28 |
| 2009-11-01 | 118136.39 |
| 2009-12-01 | 110315.17 |
| 2010-01-01 | 106682.04 |
| 2010-02-01 | 101755.78 |
| 2010-03-01 | 98610.25 |
| 2010-04-01 | 107367.79 |
| 2010-05-01 | 111126.43 |
| 2010-06-01 | 115124.86 |
| 2010-07-01 | 114202.43 |
| 2010-08-01 | 103961.94 |
| 2010-09-01 | 102617.44 |
| 2010-10-01 | 106167.56 |
| 2010-11-01 | 93698.29 |
| 2010-12-01 | 93056.77 |
| 2011-01-01 | 103460.83 |
| 2011-02-01 | 89810.78 |
| 2011-03-01 | 90330.54 |
| 2011-04-01 | 86224.61 |
| 2011-05-01 | 103435.60 |
| 2011-06-01 | 89154.63 |
| 2011-07-01 | 88290.98 |
| 2011-08-01 | 82079.98 |
| 2011-09-01 | 80236.59 |
| 2011-10-01 | 91172.81 |
| 2011-11-01 | 80320.35 |
| 2011-12-01 | 87243.26 |
| 2012-01-01 | 87272.51 |
| 2012-02-01 | 90462.31 |
| 2012-03-01 | 93036.66 |
| 2012-04-01 | 91670.56 |
| 2012-05-01 | 89280.48 |
| 2012-06-01 | 98203.37 |
| 2012-07-01 | 92091.42 |
| 2012-08-01 | 96995.21 |
| 2012-09-01 | 92846.71 |
| 2012-10-01 | 86320.40 |
| 2012-11-01 | 85063.25 |
| 2012-12-01 | 89758.97 |
| 2013-01-01 | 95140.50 |
| 2013-02-01 | 87823.37 |
| 2013-03-01 | 99451.11 |
| 2013-04-01 | 99740.27 |
| 2013-05-01 | 102759.20 |
| 2013-06-01 | 109868.61 |
| 2013-07-01 | 96633.06 |
| 2013-08-01 | 92470.17 |
| 2013-09-01 | 83518.21 |
| 2013-10-01 | 95298.71 |
| 2013-11-01 | 91124.15 |
| 2013-12-01 | 98688.85 |
| 2014-01-01 | 86102.59 |
| 2014-02-01 | 81416.85 |
| 2014-03-01 | 80848.11 |
| 2014-04-01 | 82855.53 |
| 2014-05-01 | 96952.44 |
| 2014-06-01 | 86394.46 |
| 2014-07-01 | 100880.03 |
| 2014-08-01 | 115317.46 |
| 2014-09-01 | 92229.41 |
| 2014-10-01 | 95792.90 |
| 2014-11-01 | 101713.62 |
| 2014-12-01 | 101666.14 |
| 2015-01-01 | 96769.51 |
| 2015-02-01 | 96523.71 |
| 2015-03-01 | 100332.08 |
| 2015-04-01 | 97523.76 |
| 2015-05-01 | 93502.72 |
| 2015-06-01 | 99687.39 |
| 2015-07-01 | 108252.74 |
| 2015-08-01 | 104623.40 |
| 2015-09-01 | 111129.52 |
| 2015-10-01 | 105666.28 |
| 2015-11-01 | 98744.94 |
| 2015-12-01 | 99618.20 |
| 2016-01-01 | 87311.20 |
| 2016-02-01 | 91023.70 |
| 2016-03-01 | 94316.83 |
| 2016-04-01 | 82313.73 |
| 2016-05-01 | 86670.36 |
| 2016-06-01 | 94678.90 |
| 2016-07-01 | 90902.41 |
| 2016-08-01 | 99063.37 |
| 2016-09-01 | 99423.61 |
| 2016-10-01 | 104029.50 |
| 2016-11-01 | 108087.66 |
| 2016-12-01 | 102780.14 |
| 2017-01-01 | 97209.33 |
| 2017-02-01 | 92771.18 |
| 2017-03-01 | 93727.62 |
| 2017-04-01 | 104915.20 |
| 2017-05-01 | 107881.67 |
| 2017-06-01 | 107551.01 |
| 2017-07-01 | 98009.87 |
| 2017-08-01 | 106019.27 |
| 2017-09-01 | 115592.49 |
| 2017-10-01 | 107990.96 |
| 2017-11-01 | 118263.24 |
| 2017-12-01 | 106677.35 |
| 2018-01-01 | 93947.11 |
| 2018-02-01 | 100735.05 |
| 2018-03-01 | 101110.05 |
| 2018-04-01 | 115390.07 |
| 2018-05-01 | 117309.20 |
| 2018-06-01 | 105720.31 |
| 2018-07-01 | 106133.15 |
| 2018-08-01 | 116454.12 |
| 2018-09-01 | 102437.04 |
| 2018-10-01 | 99510.41 |
| 2018-11-01 | 96380.07 |
| 2018-12-01 | 102600.42 |
| 2019-01-01 | 86673.72 |
| 2019-02-01 | 91933.06 |
| 2019-03-01 | 100666.87 |
| 2019-04-01 | 99477.14 |
| 2019-05-01 | 101178.04 |
| 2019-06-01 | 98945.58 |
| 2019-07-01 | 91213.62 |
| 2019-08-01 | 104264.47 |
| 2019-09-01 | 101619.52 |
| 2019-10-01 | 112341.76 |
| 2019-11-01 | 96893.64 |
| 2019-12-01 | 112804.51 |
| 2020-01-01 | 110384.74 |
| 2020-02-01 | 98630.36 |
| 2020-03-01 | 90724.88 |
| 2020-04-01 | 14088.39 |
| 2020-05-01 | 25016.02 |
| 2020-06-01 | 102299.51 |
| 2020-07-01 | 96350.21 |
| 2020-08-01 | 114445.77 |
| 2020-09-01 | 114172.34 |
| 2020-10-01 | 121765.33 |
| 2020-11-01 | 117479.15 |
| 2020-12-01 | 133102.68 |
| 2021-01-01 | 104635.18 |
| 2021-02-01 | 114387.74 |
| 2021-03-01 | 114309.63 |
| 2021-04-01 | 110256.92 |
| 2021-05-01 | 121602.47 |
| 2021-06-01 | 111064.56 |
| 2021-07-01 | 102432.93 |
| 2021-08-01 | 110476.86 |
| 2021-09-01 | 110491.38 |
3 Visualizar la serie de tiempo
Zinc_ts <- ts(Zinc$Zinc_TMR, start=c(2007,1), frequency =12)
Zinc_ts## Jan Feb Mar Apr May Jun Jul
## 2007 98471.46 96350.45 107589.89 93117.68 109090.20 118277.49 116272.12
## 2008 110691.37 102079.42 112278.81 106553.66 114257.44 121328.51 114810.94
## 2009 116661.85 97106.22 101387.08 104535.89 104492.31 102100.39 108021.25
## 2010 106682.04 101755.78 98610.25 107367.79 111126.43 115124.86 114202.43
## 2011 103460.83 89810.78 90330.54 86224.61 103435.60 89154.63 88290.98
## 2012 87272.51 90462.31 93036.66 91670.56 89280.48 98203.37 92091.42
## 2013 95140.50 87823.37 99451.11 99740.27 102759.20 109868.61 96633.06
## 2014 86102.59 81416.85 80848.11 82855.53 96952.44 86394.46 100880.03
## 2015 96769.51 96523.71 100332.08 97523.76 93502.72 99687.39 108252.74
## 2016 87311.20 91023.70 94316.83 82313.73 86670.36 94678.90 90902.41
## 2017 97209.33 92771.18 93727.62 104915.20 107881.67 107551.01 98009.87
## 2018 93947.11 100735.05 101110.05 115390.07 117309.20 105720.31 106133.15
## 2019 86673.72 91933.06 100666.87 99477.14 101178.04 98945.58 91213.62
## 2020 110384.74 98630.36 90724.88 14088.39 25016.02 102299.51 96350.21
## 2021 104635.18 114387.74 114309.63 110256.92 121602.47 111064.56 102432.93
## Aug Sep Oct Nov Dec
## 2007 101161.38 98626.96 107174.43 87056.73 98928.10
## 2008 126973.22 110370.43 118702.77 105057.49 123996.17
## 2009 115596.97 93810.62 118446.28 118136.39 110315.17
## 2010 103961.94 102617.44 106167.56 93698.29 93056.77
## 2011 82079.98 80236.59 91172.81 80320.35 87243.26
## 2012 96995.21 92846.71 86320.40 85063.25 89758.97
## 2013 92470.17 83518.21 95298.71 91124.15 98688.85
## 2014 115317.46 92229.41 95792.90 101713.62 101666.14
## 2015 104623.40 111129.52 105666.28 98744.94 99618.20
## 2016 99063.37 99423.61 104029.50 108087.66 102780.14
## 2017 106019.27 115592.49 107990.96 118263.24 106677.35
## 2018 116454.12 102437.04 99510.41 96380.07 102600.42
## 2019 104264.47 101619.52 112341.76 96893.64 112804.51
## 2020 114445.77 114172.34 121765.33 117479.15 133102.68
## 2021 110476.86 110491.38ggplot(Zinc_ts, aes(x= time(Zinc_ts), y= Zinc_ts))+
geom_line(size=0.6)+
scale_x_continuous(breaks = seq(2007, 2021, 1))+
xlab("Fecha")+
ylab("Zinc (TMR)")+
theme_economist()+
labs(title = "Total nacional de producción de Zinc",
subtitle = "Toneladas métricas de contenido recuperado (TMR)",
caption = "Elaboración propia")## Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.
4 Descomposición temporal
plot(decompose(Zinc_ts)$random,
xlab="Fecha", ylab="Aleatorio")
plot(decompose(Zinc_ts)$seasonal,
xlab="Fecha", ylab="Estacional")
5 Identificación de estacionariedad
5.1 Correlogramas
Acf(Zinc_ts, lag.max = 24, plot=T, main="FAS del Zinc")
Pacf(Zinc_ts, lag.max = 24, plot=T, main="FAP del Zinc")
5.2 Pruebas ADF, PP y KPSS
5.2.1 ADF
adf.test(Zinc_ts)## Warning in adf.test(Zinc_ts): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: Zinc_ts
## Dickey-Fuller = -4.1848, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryADF<-ur.df(y=Zinc_ts,
type="trend",
lags=5,
selectlags = c("AIC"))
summary(ADF)##
## ###############################################
## # 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
## -80515 -5603 -322 5690 46878
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.074e+04 7.332e+03 5.557 1.06e-07 ***
## z.lag.1 -4.055e-01 6.967e-02 -5.820 2.94e-08 ***
## tt -8.047e-01 1.729e+01 -0.047 0.963
## z.diff.lag -6.045e-03 7.723e-02 -0.078 0.938
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11130 on 167 degrees of freedom
## Multiple R-squared: 0.2056, Adjusted R-squared: 0.1913
## F-statistic: 14.41 on 3 and 167 DF, p-value: 2.157e-08
##
##
## Value of test-statistic is: -5.8195 11.3379 17.0055
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -3.99 -3.43 -3.13
## phi2 6.22 4.75 4.07
## phi3 8.43 6.49 5.475.2.2 Phillips-Perron
pp.test(Zinc_ts)## Warning in pp.test(Zinc_ts): p-value smaller than printed p-value##
## Phillips-Perron Unit Root Test
##
## data: Zinc_ts
## Dickey-Fuller Z(alpha) = -72.406, Truncation lag parameter = 4, p-value
## = 0.01
## alternative hypothesis: stationaryPP<-ur.pp(Zinc_ts,
type="Z-tau",
model="trend",
lags="short")
summary(PP)##
## ##################################
## # Phillips-Perron Unit Root Test #
## ##################################
##
## Test regression with intercept and trend
##
##
## Call:
## lm(formula = y ~ y.l1 + trend)
##
## Residuals:
## Min 1Q Median 3Q Max
## -80370 -5734 -229 6055 46496
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.144e+04 6.249e+03 6.631 4.1e-10 ***
## y.l1 5.881e-01 6.166e-02 9.538 < 2e-16 ***
## trend -4.743e+00 1.647e+01 -0.288 0.774
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11080 on 173 degrees of freedom
## Multiple R-squared: 0.3464, Adjusted R-squared: 0.3388
## F-statistic: 45.84 on 2 and 173 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic, type: Z-tau is: -6.6769
##
## aux. Z statistics
## Z-tau-mu 6.6430
## Z-tau-beta -0.2884
##
## Critical values for Z statistics:
## 1pct 5pct 10pct
## critical values -4.012792 -3.436124 -3.1418835.2.3 KPSS
kpss.test(Zinc_ts)## Warning in kpss.test(Zinc_ts): p-value greater than printed p-value##
## KPSS Test for Level Stationarity
##
## data: Zinc_ts
## KPSS Level = 0.24888, Truncation lag parameter = 4, p-value = 0.1KPSS<-ur.kpss(Zinc_ts,
type="tau",
lags="short")
summary(KPSS)##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 4 lags.
##
## Value of test-statistic is: 0.2353
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.2165.3 Primera diferencia
Zinc_d <- diff(Zinc_ts)
ggplot(Zinc_d, aes(x= time(Zinc_d), y= Zinc_d))+
geom_line(size=0.6)+
scale_x_continuous(breaks = seq(2007, 2021, 1))+
xlab("Fecha")+
ylab("Diferencia")+
theme_economist()+
labs(title = "Primera diferencia del Zinc",
caption = "Elaboración propia")## Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.
5.4 Pruebas ADF, PP y KPSS para Zinc_d
5.4.1 ADF
adf.test(Zinc_d)## Warning in adf.test(Zinc_d): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: Zinc_d
## Dickey-Fuller = -7.7237, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryADF<-ur.df(y=Zinc_d,
type="trend",
lags=5,
selectlags = c("AIC"))
summary(ADF)##
## ###############################################
## # 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
## -79557 -5969 -489 6452 55193
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -846.78484 1888.98745 -0.448 0.654547
## z.lag.1 -1.89403 0.19886 -9.525 < 2e-16 ***
## tt 9.05770 18.35331 0.494 0.622306
## z.diff.lag1 0.61229 0.16515 3.708 0.000286 ***
## z.diff.lag2 0.33732 0.12236 2.757 0.006500 **
## z.diff.lag3 0.21492 0.07586 2.833 0.005189 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11730 on 164 degrees of freedom
## Multiple R-squared: 0.6414, Adjusted R-squared: 0.6305
## F-statistic: 58.68 on 5 and 164 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -9.5245 30.2471 45.367
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -3.99 -3.43 -3.13
## phi2 6.22 4.75 4.07
## phi3 8.43 6.49 5.475.4.2 Phillips-Perron
pp.test(Zinc_d)## Warning in pp.test(Zinc_d): p-value smaller than printed p-value##
## Phillips-Perron Unit Root Test
##
## data: Zinc_d
## Dickey-Fuller Z(alpha) = -178.21, Truncation lag parameter = 4, p-value
## = 0.01
## alternative hypothesis: stationaryPP<-ur.pp(Zinc_d,
type="Z-tau",
model="trend",
lags="short")
summary(PP)##
## ##################################
## # Phillips-Perron Unit Root Test #
## ##################################
##
## Test regression with intercept and trend
##
##
## Call:
## lm(formula = y ~ y.l1 + trend)
##
## Residuals:
## Min 1Q Median 3Q Max
## -78533 -5764 157 6404 79426
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 94.78090 920.24409 0.103 0.91809
## y.l1 -0.21459 0.07447 -2.882 0.00446 **
## trend 1.49128 18.21563 0.082 0.93485
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 12170 on 172 degrees of freedom
## Multiple R-squared: 0.04608, Adjusted R-squared: 0.03498
## F-statistic: 4.154 on 2 and 172 DF, p-value: 0.0173
##
##
## Value of test-statistic, type: Z-tau is: -17.7959
##
## aux. Z statistics
## Z-tau-mu 0.1569
## Z-tau-beta 0.0902
##
## Critical values for Z statistics:
## 1pct 5pct 10pct
## critical values -4.01308 -3.436262 -3.1419655.4.3 KPSS
kpss.test(Zinc_d)## Warning in kpss.test(Zinc_d): p-value greater than printed p-value##
## KPSS Test for Level Stationarity
##
## data: Zinc_d
## KPSS Level = 0.018696, Truncation lag parameter = 4, p-value = 0.1KPSS<-ur.kpss(Zinc_d,
type="tau",
lags="short")
summary(KPSS)##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 4 lags.
##
## Value of test-statistic is: 0.0163
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.2165.5 Correlogramas de Zinc_d
Acf(Zinc_d, lag.max = 24, plot=T, main="FAS de d1 del Zinc")
Pacf(Zinc_d, lag.max = 24, plot=T, main = "FAP de d1 del Zinc")
6 Arima
# Arima(1,1,1)
# Arima(1,1,2)
# Arima(1,1,4)
# Arima(2,1,1)
# Arima(2,1,2)
# Arima(2,1,4)
# Arima(4,1,1)
# Arima(4,1,2)
# Arima(4,1,4)7 Identificación estacional
ggseasonplot(Zinc_ts, year.labels = TRUE, year.labels.left = TRUE)+
ggtitle("Zinc 2007-2021")
8 Prueba de raíz unitaria estacional
8.1 DHF
8.1.1 Diferencia estacional
Zincs<-diff(Zinc_ts,12)
Zincs<-na.omit(Zincs)8.1.2 Rezagando en 12 periodos
Zinc12<-stats::lag(Zinc_ts,-12)8.1.3 Regresión
DHF <- lm(Zincs ~ Zinc12[13:177])
summary(DHF)##
## Call:
## lm(formula = Zincs ~ Zinc12[13:177])
##
## Residuals:
## Min 1Q Median 3Q Max
## -25938 -7128 -591 6755 86405
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8.930e+04 7.605e+03 -11.74 <2e-16 ***
## Zinc12[13:177] 8.985e-01 7.509e-02 11.97 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13330 on 163 degrees of freedom
## Multiple R-squared: 0.4676, Adjusted R-squared: 0.4644
## F-statistic: 143.2 on 1 and 163 DF, p-value: < 2.2e-168.2 Prueba HEGY
hegy<-hegy.test(Zinc_ts)
summary(hegy)##
## HEGY test for unit roots
##
## data: Zinc_ts
##
## Fitted model
## ------------
##
## Call:
## lm(formula = dx ~ 0 + ypi + xreg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -77880 -5362 -771 6164 43839
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## ypiYpi1 -2.531e-02 1.033e-02 -2.450 0.015438 *
## ypiYpi2 -1.668e-01 4.342e-02 -3.841 0.000180 ***
## ypiYpi3 -7.110e-02 1.972e-02 -3.605 0.000423 ***
## ypiYpi4 -4.083e-02 1.796e-02 -2.273 0.024404 *
## ypiYpi5 -1.263e-01 2.949e-02 -4.282 3.27e-05 ***
## ypiYpi6 -7.712e-02 2.812e-02 -2.742 0.006833 **
## ypiYpi7 -1.770e-01 3.923e-02 -4.512 1.28e-05 ***
## ypiYpi8 -1.315e-01 3.905e-02 -3.367 0.000964 ***
## ypiYpi9 -1.810e-01 3.995e-02 -4.529 1.19e-05 ***
## ypiYpi10 6.309e-02 4.116e-02 1.533 0.127392
## ypiYpi11 -2.072e-01 4.474e-02 -4.631 7.77e-06 ***
## ypiYpi12 9.972e-02 4.881e-02 2.043 0.042775 *
## xreg 3.058e+04 1.242e+04 2.462 0.014931 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11320 on 152 degrees of freedom
## Multiple R-squared: 0.6427, Adjusted R-squared: 0.6121
## F-statistic: 21.03 on 13 and 152 DF, p-value: < 2.2e-16
##
## Test statistic
## --------------
## statistic p-value
## t_1 -2.4496 0.1022
## t_2 -3.8405 1e-04 ***
## F_3:4 9.4647 1e-04 ***
## F_5:6 13.877 0 ***
## F_7:8 17.4612 0 ***
## F_9:10 11.7354 0 ***
## F_11:12 13.3985 0 ***
## F_2:12 24.6763 0 ***
## F_1:12 22.7063 0 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Deterministic terms: constant
## Lag selection criterion and order: fixed, 0
## P-values: based on response surface regressions9 Diferencia estacional
Zincs## Jan Feb Mar Apr May Jun
## 2008 12219.9024 5728.9739 4688.9161 13435.9815 5167.2424 3051.0162
## 2009 5970.4866 -4973.2001 -10891.7292 -2017.7707 -9765.1331 -19228.1166
## 2010 -9979.8154 4649.5509 -2776.8319 2831.9007 6634.1240 13024.4668
## 2011 -3221.2087 -11944.9962 -8279.7059 -21143.1782 -7690.8298 -25970.2222
## 2012 -16188.3184 651.5288 2706.1154 5445.9528 -14155.1196 9048.7333
## 2013 7867.9932 -2638.9363 6414.4578 8069.7047 13478.7214 11665.2426
## 2014 -9037.9137 -6406.5200 -18603.0037 -16884.7381 -5806.7586 -23474.1537
## 2015 10666.9194 15106.8570 19483.9669 14668.2303 -3449.7194 13292.9372
## 2016 -9458.3109 -5500.0082 -6015.2486 -15210.0265 -6832.3602 -5008.4966
## 2017 9898.1275 1747.4795 -589.2143 22601.4685 21211.3077 12872.1124
## 2018 -3262.2205 7963.8700 7382.4356 10474.8643 9427.5273 -1830.7025
## 2019 -7273.3834 -8801.9933 -443.1798 -15912.9290 -16131.1642 -6774.7267
## 2020 23711.0162 6697.2999 -9941.9929 -85388.7440 -76162.0127 3353.9336
## 2021 -5749.5536 15757.3840 23584.7547 96168.5236 96586.4461 8765.0495
## Jul Aug Sep Oct Nov Dec
## 2008 -1461.1779 25811.8377 11743.4686 11528.3374 18000.7641 25068.0693
## 2009 -6789.6937 -11376.2533 -16559.8059 -256.4812 13078.8945 -13680.9990
## 2010 6181.1835 -11635.0231 8806.8159 -12278.7255 -24438.1017 -17258.3975
## 2011 -25911.4511 -21881.9660 -22380.8441 -14994.7447 -13377.9349 -5813.5105
## 2012 3800.4384 14915.2347 12610.1132 -4852.4112 4742.8978 2515.7113
## 2013 4541.6391 -4525.0407 -9328.5017 8978.3035 6060.9003 8929.8785
## 2014 4246.9747 22847.2903 8711.2025 494.1979 10589.4669 2977.2892
## 2015 7372.7057 -10694.0581 18900.1140 9873.3803 -2968.6787 -2047.9422
## 2016 -17350.3316 -5560.0278 -11705.9171 -1636.7886 9342.7258 3161.9380
## 2017 7107.4596 6955.8949 16168.8851 3961.4592 10175.5784 3897.2118
## 2018 8123.2795 10434.8478 -13155.4474 -8480.5447 -21883.1747 -4076.9283
## 2019 -14919.5308 -12189.6429 -817.5220 12831.3490 513.5752 10204.0875
## 2020 5136.5948 10181.2910 12552.8135 9423.5692 20585.5138 20298.1740
## 2021 6082.7180 -3968.9038 -3680.9598ggplot(Zincs, aes(x=time(Zincs), y=Zincs))+
geom_line()+
scale_x_continuous(breaks = seq(2007, 2021,1))+
xlab("Fecha")+
ylab("Diferencia estacional")+
theme_economist()+
labs(title = "Diferencia estacional del Zinc",
caption = "Elaboración propia")## Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.
Acf(Zincs, lag.max = 24, plot=T, main="FAS del PIB en d(12)")
Pacf(Zincs, lag.max = 24, plot=T, main="FAP del PIB en d(12)")
10 Pruebas de estacionariedad para Zincs
adf.test(Zincs)## Warning in adf.test(Zincs): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: Zincs
## Dickey-Fuller = -4.2172, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationarypp.test(Zincs)## Warning in pp.test(Zincs): p-value smaller than printed p-value##
## Phillips-Perron Unit Root Test
##
## data: Zincs
## Dickey-Fuller Z(alpha) = -66.07, Truncation lag parameter = 4, p-value
## = 0.01
## alternative hypothesis: stationarykpss.test(Zincs)## Warning in kpss.test(Zincs): p-value greater than printed p-value##
## KPSS Test for Level Stationarity
##
## data: Zincs
## KPSS Level = 0.16433, Truncation lag parameter = 4, p-value = 0.111 identificación de los ordenes estacionales
# Q= 1,2; D=0; P=112 modelos
# SArima(1,1,1)(1,0,1)[12]
mod1<-stats::arima(Zinc_ts, order = c(1,1,1),
seasonal =list(order=c(1,0,1), period=12))
# SArima(1,1,2)(1,0,1)[12]
mod2<-stats::arima(Zinc_ts, order = c(1,1,2),
seasonal =list(order=c(1,0,1), period=12))
# SArima(1,1,4)(1,0,1)[12]
mod3<-stats::arima(Zinc_ts, order = c(1,1,4),
seasonal =list(order=c(1,0,1), period=12))
# SArima(2,1,1)(1,0,1)[12]
mod4<-stats::arima(Zinc_ts, order = c(2,1,1),
seasonal =list(order=c(1,0,1), period=12))
# SArima(2,1,2)(1,0,1)[12]
mod5<-stats::arima(Zinc_ts, order = c(2,1,2),
seasonal =list(order=c(1,0,1), period=12))
# SArima(2,1,4)(1,0,1)[12]
mod6<-stats::arima(Zinc_ts, order = c(2,1,4),
seasonal =list(order=c(1,0,1), period=12))
# SArima(4,1,1)(1,0,1)[12]
mod7<-stats::arima(Zinc_ts, order = c(4,1,1),
seasonal =list(order=c(1,0,1), period=12))
# SArima(4,1,2)(1,0,1)[12]
mod8<-stats::arima(Zinc_ts, order = c(4,1,2),
seasonal =list(order=c(1,0,1), period=12))
# SArima(4,1,4)(1,0,1)[12]
mod9<-stats::arima(Zinc_ts, order = c(4,1,4),
seasonal =list(order=c(1,0,1), period=12))## Warning in log(s2): Se han producido NaNs# SArima(1,1,1)(1,0,2)[12]
mod10<-stats::arima(Zinc_ts, order = c(1,1,1),
seasonal =list(order=c(1,0,2), period=12))
# SArima(1,1,2)(1,0,2)[12]
mod11<-stats::arima(Zinc_ts, order = c(1,1,2),
seasonal =list(order=c(1,0,2), period=12))
# SArima(1,1,4)(1,0,2)[12]
mod12<-stats::arima(Zinc_ts, order = c(1,1,4),
seasonal =list(order=c(1,0,2), period=12))
# SArima(2,1,1)(1,0,2)[12]
mod13<-stats::arima(Zinc_ts, order = c(2,1,1),
seasonal =list(order=c(1,0,2), period=12))
# SArima(2,1,2)(1,0,2)[12]
mod14<-stats::arima(Zinc_ts, order = c(2,1,2),
seasonal =list(order=c(1,0,2), period=12))
# SArima(2,1,4)(1,0,2)[12]
mod15<-stats::arima(Zinc_ts, order = c(2,1,4),
seasonal =list(order=c(1,0,2), period=12))
# SArima(4,1,1)(1,0,2)[12]
mod16<-stats::arima(Zinc_ts, order = c(4,1,1),
seasonal =list(order=c(1,0,2), period=12))
# SArima(4,1,2)(1,0,2)[12]
mod17<-stats::arima(Zinc_ts, order = c(4,1,2),
seasonal =list(order=c(1,0,2), period=12))
# SArima(4,1,4)(1,0,2)[12]
mod18<-stats::arima(Zinc_ts, order = c(4,1,4),
seasonal =list(order=c(1,0,2), period=12))
#auto
auto<-auto.arima(Zinc_ts, seasonal = TRUE)12.1 Bondad de ajuste
AIC<-c(mod1$aic,mod2$aic,
mod3$aic,mod4$aic,
mod5$aic,mod6$aic,
mod7$aic,mod8$aic,
mod9$aic,mod10$aic,
mod11$aic,mod12$aic,
mod13$aic,mod14$aic,
mod15$aic,mod16$aic,
mod17$aic,mod18$aic,auto$aic)
modelos<-c("mod1","mod2",
"mod3","mod4",
"mod5","mod6",
"mod7","mod8",
"mod9","mod10",
"mod11","mod12",
"mod13","mod14",
"mod15","mod16",
"mod17","mod18","auto")
RMSE1<-rmse(Zinc_ts, fitted.values(mod1))
RMSE2<-rmse(Zinc_ts, fitted.values(mod2))
RMSE3<-rmse(Zinc_ts, fitted.values(mod3))
RMSE4<-rmse(Zinc_ts, fitted.values(mod4))
RMSE5<-rmse(Zinc_ts, fitted.values(mod5))
RMSE6<-rmse(Zinc_ts, fitted.values(mod6))
RMSE7<-rmse(Zinc_ts, fitted.values(mod7))
RMSE8<-rmse(Zinc_ts, fitted.values(mod8))
RMSE9<-rmse(Zinc_ts, fitted.values(mod9))
RMSE10<-rmse(Zinc_ts, fitted.values(mod10))
RMSE11<-rmse(Zinc_ts, fitted.values(mod11))
RMSE12<-rmse(Zinc_ts, fitted.values(mod12))
RMSE13<-rmse(Zinc_ts, fitted.values(mod13))
RMSE14<-rmse(Zinc_ts, fitted.values(mod14))
RMSE15<-rmse(Zinc_ts, fitted.values(mod15))
RMSE16<-rmse(Zinc_ts, fitted.values(mod16))
RMSE17<-rmse(Zinc_ts, fitted.values(mod17))
RMSE18<-rmse(Zinc_ts, fitted.values(mod18))
RMSEa<-rmse(Zinc_ts, fitted.values(auto))
RMSE<-c(RMSE1,RMSE2,
RMSE3,RMSE4,
RMSE5,RMSE6,
RMSE7,RMSE8,
RMSE9,RMSE10,
RMSE11,RMSE12,
RMSE13,RMSE14,
RMSE15,RMSE16,
RMSE17,RMSE18,RMSEa)
TablaBA<-cbind(modelos, AIC, RMSE)
TablaBA## modelos AIC RMSE
## [1,] "mod1" "3784.19482028127" "10772.8373023539"
## [2,] "mod2" "3786.11767301223" "10762.802261832"
## [3,] "mod3" "3787.1129709152" "10699.5459244005"
## [4,] "mod4" "3786.14877456933" "10774.1477828802"
## [5,] "mod5" "3785.84135526016" "10669.5276811543"
## [6,] "mod6" "3791.93602529138" "10883.348936175"
## [7,] "mod7" "3787.88338186222" "10611.3926414245"
## [8,] "mod8" "3791.73897253402" "10875.8724687957"
## [9,] "mod9" "3789.71314473983" "10484.0238676865"
## [10,] "mod10" "3786.01894557604" "10742.8579815702"
## [11,] "mod11" "3787.91479088499" "10737.0079121452"
## [12,] "mod12" "3789.08813312976" "10688.9983823294"
## [13,] "mod13" "3787.96751303484" "10758.4266161549"
## [14,] "mod14" "3787.80614006213" "10719.7558126449"
## [15,] "mod15" "3793.88300773591" "10881.4741184415"
## [16,] "mod16" "3789.7540841053" "10574.208882615"
## [17,] "mod17" "3793.68836026902" "10874.9104979863"
## [18,] "mod18" "3795.0536358407" "10701.5435579717"
## [19,] "auto" "3801.59178098687" "10958.6493279383"En este ejemplo el modelo 9 tiene menor RMSE
13 Evaluación
13.1 Coeficientes
coeftest(mod9) ##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.787221 0.208326 -3.7788 0.0001576 ***
## ar2 -0.687345 0.132917 -5.1713 2.325e-07 ***
## ar3 0.157078 0.157129 0.9997 0.3174697
## ar4 0.283806 0.161105 1.7616 0.0781337 .
## ma1 0.458860 0.185711 2.4708 0.0134801 *
## ma2 0.124441 0.098667 1.2612 0.2072263
## ma3 -0.717016 0.083999 -8.5360 < 2.2e-16 ***
## ma4 -0.680528 0.191045 -3.5621 0.0003679 ***
## sar1 0.994803 0.043065 23.0999 < 2.2e-16 ***
## sma1 -0.968705 0.136597 -7.0917 1.325e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1autoplot(mod9)
13.2 Residuales
log(177)## [1] 5.17615res<-mod9$residuals
Box.test(res, lag=6, type="Ljung-Box")##
## Box-Ljung test
##
## data: res
## X-squared = 0.19956, df = 6, p-value = 0.9998Box.test(res, lag=6, type = "Box-Pierce")##
## Box-Pierce test
##
## data: res
## X-squared = 0.19398, df = 6, p-value = 0.9999ggtsdisplay(res)
13.3 Pruebas de estacionariedad de los residuos
adf.test(res) ## Warning in adf.test(res): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: res
## Dickey-Fuller = -5.381, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationarypp.test(res) ## Warning in pp.test(res): p-value smaller than printed p-value##
## Phillips-Perron Unit Root Test
##
## data: res
## Dickey-Fuller Z(alpha) = -180.11, Truncation lag parameter = 4, p-value
## = 0.01
## alternative hypothesis: stationarykpss.test(res) ## Warning in kpss.test(res): p-value greater than printed p-value##
## KPSS Test for Level Stationarity
##
## data: res
## KPSS Level = 0.18539, Truncation lag parameter = 4, p-value = 0.114 Alternativa para revisar la heterocedasticidad
res2<-res^2
ggtsdisplay(res2) 
Box.test(res2, lag=6, type="Ljung-Box")##
## Box-Ljung test
##
## data: res2
## X-squared = 15.945, df = 6, p-value = 0.0140515 Anexo: JB
jarque.bera.test(res)##
## Jarque Bera Test
##
## data: res
## X-squared = 1423.9, df = 2, p-value < 2.2e-16qqnorm(res);qqline(res)
16 Pronostico
fcast<-forecast(mod9, h=5, level=95)
autoplot(fcast, include = 177)
fcast## Point Forecast Lo 95 Hi 95
## Oct 2021 110488.13 89702.60 131273.7
## Nov 2021 103263.18 78222.20 128304.2
## Dec 2021 106614.17 80434.91 132793.4
## Jan 2022 104474.80 77687.55 131262.1
## Feb 2022 99722.78 72840.01 126605.6