title: “Taller 3_Series de Tiempo_Inés Esperanza Rodríguez Rueda” output: html_document date: “2024-09-10” Taller # 3 - Series de Tiempo
Establecer si la serie es estacionaria en sus datos originales gráficamente y validando con Dickey-Fuller.
library(readxl)## Warning: package 'readxl' was built under R version 4.4.1Datos <- read_excel("C:/Users/Administrador/Downloads/Series de Tiempo/Taller 3/Libro3.xlsx")
View(Datos)
attach(Datos)
summary(Datos)## Valores PIB
## Min. : 78.88
## 1st Qu.: 86.58
## Median : 99.17
## Mean : 99.85
## 3rd Qu.:107.19
## Max. :134.81Convertir objeto de serie de tiempo en R |
Datos.ts=ts(Datos,start = c (1,1), end = c(9,12), frequency = 12)
Datos.ts## Jan Feb Mar Apr May Jun Jul
## 1 78.88267 79.35263 79.45924 79.61915 79.86456 80.08365 80.03355
## 2 80.82213 81.30427 81.53248 81.78099 81.99070 82.20187 82.29903
## 3 83.37303 84.14541 84.49016 84.89067 85.22691 85.37703 85.59892
## 4 88.44274 89.46566 90.09688 90.45238 90.97727 91.32911 91.59756
## 5 93.74717 94.75693 95.16495 95.70566 95.92626 96.04420 96.03690
## 6 97.38363 98.20417 98.44586 98.82558 99.03045 99.15752 99.18659
## 7 100.52529 101.18411 101.63472 102.08028 102.29658 102.47209 102.59699
## 8 104.16016 104.83761 105.32962 105.38722 105.04157 104.69569 104.86702
## 9 106.01345 106.65862 107.08805 107.48509 107.97680 108.20861 108.64259
## Aug Sep Oct Nov Dec
## 1 80.13879 80.35038 80.39552 80.30780 80.51214
## 2 82.44722 82.54750 82.70804 82.86633 83.08325
## 3 85.88003 86.28199 86.67747 87.21648 87.72510
## 4 91.89515 92.23799 92.41316 92.58247 92.96870
## 5 96.13943 96.37798 96.43003 96.59962 96.87858
## 6 99.33403 99.45646 99.56643 99.72406 100.00000
## 7 102.80122 103.09592 103.24693 103.41952 103.77564
## 8 105.01001 105.49416 105.54625 105.44685 105.72996
## 9 109.11453 109.48749 109.43891 109.84225 110.40086plot(Datos.ts)
View(Datos.ts)start(Datos.ts)## [1] 1 1end(Datos.ts)## [1] 9 12plot(Datos.ts, main= "Serie de tiempo",ylab= "PIB" , col="red")
acf(Datos.ts)
pacf(Datos.ts)
acf(ts(Datos.ts, frequency = 1))
pacf(ts(Datos.ts, frequency = 1))
serielog <- log(Datos.ts)
serielog## Jan Feb Mar Apr May Jun Jul Aug
## 1 4.367962 4.373902 4.375244 4.377255 4.380332 4.383072 4.382446 4.383760
## 2 4.392251 4.398199 4.401001 4.404045 4.406606 4.409178 4.410359 4.412158
## 3 4.423325 4.432546 4.436635 4.441364 4.445317 4.447077 4.449673 4.452951
## 4 4.482355 4.493855 4.500886 4.504824 4.510610 4.514470 4.517405 4.520648
## 5 4.540601 4.551315 4.555612 4.561277 4.563580 4.564809 4.564733 4.565800
## 6 4.578658 4.587049 4.589507 4.593357 4.595427 4.596710 4.597003 4.598488
## 7 4.610409 4.616942 4.621385 4.625760 4.627876 4.629590 4.630809 4.632797
## 8 4.645930 4.652413 4.657095 4.657641 4.654356 4.651058 4.652693 4.654056
## 9 4.663566 4.669633 4.673651 4.677352 4.681916 4.684061 4.688064 4.692398
## Sep Oct Nov Dec
## 1 4.386397 4.386958 4.385867 4.388408
## 2 4.413374 4.415317 4.417229 4.419843
## 3 4.457621 4.462194 4.468393 4.474208
## 4 4.524372 4.526269 4.528100 4.532263
## 5 4.568278 4.568818 4.570575 4.573458
## 6 4.599720 4.600825 4.602407 4.605170
## 7 4.635660 4.637124 4.638794 4.642231
## 8 4.658656 4.659149 4.658207 4.660888
## 9 4.695810 4.695367 4.699045 4.704118plot(serielog)
##Estacionareidad
library(forecast)## Warning: package 'forecast' was built under R version 4.4.1## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoondiffs(Datos.ts)## [1] 1#Prueba DickeyFuller
library(tseries)## Warning: package 'tseries' was built under R version 4.4.1df1 <- adf.test(Datos.ts, k=0)
df1##
## Augmented Dickey-Fuller Test
##
## data: Datos.ts
## Dickey-Fuller = -1.307, Lag order = 0, p-value = 0.8641
## alternative hypothesis: stationary##Modelo 1 sin intercepto
modelo1 <- arima(Datos.ts,order = c(2,1,1))
summary(modelo1)##
## Call:
## arima(x = Datos.ts, order = c(2, 1, 1))
##
## Coefficients:
## ar1 ar2 ma1
## 1.5662 -0.5662 -0.9931
## s.e. 0.0825 0.0825 0.0124
##
## sigma^2 estimated as 0.03788: log likelihood = 21.5, aic = -34.99
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.01598044 0.1938687 0.1495031 0.01853999 0.1575178 0.4776965
## ACF1
## Training set 0.03939188tsdiag(modelo1)
Box.test(residuals(modelo1), type="Ljung-Box")##
## Box-Ljung test
##
## data: residuals(modelo1)
## X-squared = 0.17228, df = 1, p-value = 0.6781##Modelo 2 con intercepto
modelo2 <- arima(Datos.ts,order = c(0,1,0))
summary(modelo2)##
## Call:
## arima(x = Datos.ts, order = c(0, 1, 0))
##
##
## sigma^2 estimated as 0.1403: log likelihood = -46.76, aic = 95.53
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.2925655 0.3729432 0.3107992 0.3113973 0.3294991 0.9930745
## ACF1
## Training set 0.5356554tsdiag(modelo2)
Box.test(residuals(modelo2), type="Ljung-Box")##
## Box-Ljung test
##
## data: residuals(modelo2)
## X-squared = 31.857, df = 1, p-value = 1.66e-08##Modelo 3 con intercepto
modelo3 <- arima(Datos.ts,order = c(1,0,0))
summary(modelo3)##
## Call:
## arima(x = Datos.ts, order = c(1, 0, 0))
##
## Coefficients:
## ar1 intercept
## 0.9997 94.5365
## s.e. 0.0004 15.4835
##
## sigma^2 estimated as 0.1403: log likelihood = -50.93, aic = 107.86
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.2882731 0.3745962 0.3133589 0.3056643 0.3325341 1.001253
## ACF1
## Training set 0.4887699tsdiag(modelo3)
Box.test(residuals(modelo3), type="Ljung-Box")##
## Box-Ljung test
##
## data: residuals(modelo3)
## X-squared = 26.524, df = 1, p-value = 2.603e-07###Modelo 4 (2,0,0) =AR(2)
modelo4 <- arima(Datos.ts,order = c(2,0,0))
summary(modelo4)##
## Call:
## arima(x = Datos.ts, order = c(2, 0, 0))
##
## Coefficients:## Warning in sqrt(diag(x$var.coef)): Se han producido NaNs## ar1 ar2 intercept
## 0 1 94.5842
## s.e. NaN NaN NaN
##
## sigma^2 estimated as 0.499: log likelihood = -116.72, aic = 241.44
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.5741692 0.7064288 0.5972177 0.6095275 0.6317467 1.908247
## ACF1
## Training set 0.7410155tsdiag(modelo4)
Box.test(residuals(modelo4), type="Ljung-Box")##
## Box-Ljung test
##
## data: residuals(modelo4)
## X-squared = 60.966, df = 1, p-value = 5.773e-15###Modelo 5 (0,0,1) =MA(1)
modelo5 <- arima(Datos.ts,order = c(0,0,1))
summary(modelo5)##
## Call:
## arima(x = Datos.ts, order = c(0, 0, 1))
##
## Coefficients:
## ma1 intercept
## 1.0000 94.3967
## s.e. 0.0309 0.9481
##
## sigma^2 estimated as 24.5: log likelihood = -328.31, aic = 662.62
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.0493165 4.949277 4.344516 -0.5015835 4.703141 13.88172 0.9411687tsdiag(modelo5)
Box.test(residuals(modelo5), type="Ljung-Box")##
## Box-Ljung test
##
## data: residuals(modelo5)
## X-squared = 98.348, df = 1, p-value < 2.2e-16###Modelo 6 (1,0,1) =ARMA(1.1)
modelo6 <- arima(Datos.ts,order = c(1,0,1))
summary(modelo6)##
## Call:
## arima(x = Datos.ts, order = c(1, 0, 1))
##
## Coefficients:
## ar1 ma1 intercept
## 0.9997 0.6307 94.4463
## s.e. 0.0005 0.0658 17.0135
##
## sigma^2 estimated as 0.07398: log likelihood = -17.03, aic = 42.05
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.1772834 0.2719857 0.2185324 0.1880787 0.2318456 0.698261
## ACF1
## Training set -0.09403431tsdiag(modelo6)
Box.test(residuals(modelo6), type="Ljung-Box")##
## Box-Ljung test
##
## data: residuals(modelo6)
## X-squared = 0.98176, df = 1, p-value = 0.3218Box.test(residuals(modelo1), type="Ljung-Box") #pasa##
## Box-Ljung test
##
## data: residuals(modelo1)
## X-squared = 0.17228, df = 1, p-value = 0.6781Box.test(residuals(modelo2), type="Ljung-Box") #pasa##
## Box-Ljung test
##
## data: residuals(modelo2)
## X-squared = 31.857, df = 1, p-value = 1.66e-08Box.test(residuals(modelo3), type="Ljung-Box") #no pasa##
## Box-Ljung test
##
## data: residuals(modelo3)
## X-squared = 26.524, df = 1, p-value = 2.603e-07Box.test(residuals(modelo4), type="Ljung-Box") #no pasa##
## Box-Ljung test
##
## data: residuals(modelo4)
## X-squared = 60.966, df = 1, p-value = 5.773e-15Box.test(residuals(modelo5), type="Ljung-Box") #no pasa##
## Box-Ljung test
##
## data: residuals(modelo5)
## X-squared = 98.348, df = 1, p-value < 2.2e-16Box.test(residuals(modelo6), type="Ljung-Box") #pasa##
## Box-Ljung test
##
## data: residuals(modelo6)
## X-squared = 0.98176, df = 1, p-value = 0.3218Pronosticos arima
pronostico1 <- forecast::forecast(modelo2, h=12)
pronostico1## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 10 110.4009 109.9208 110.8809 109.6667 111.1351
## Feb 10 110.4009 109.7219 111.0798 109.3625 111.4392
## Mar 10 110.4009 109.5694 111.2324 109.1292 111.6726
## Apr 10 110.4009 109.4407 111.3610 108.9324 111.8693
## May 10 110.4009 109.3274 111.4743 108.7591 112.0426
## Jun 10 110.4009 109.2249 111.5768 108.6024 112.1993
## Jul 10 110.4009 109.1307 111.6710 108.4583 112.3434
## Aug 10 110.4009 109.0430 111.7587 108.3242 112.4775
## Sep 10 110.4009 108.9606 111.8411 108.1982 112.6035
## Oct 10 110.4009 108.8827 111.9190 108.0791 112.7226
## Nov 10 110.4009 108.8086 111.9931 107.9658 112.8360
## Dec 10 110.4009 108.7378 112.0639 107.8575 112.9442plot(pronostico1)
Modelo 1 y modelo 2 con estacionalidad
library(forecast)
auto.arima(Datos.ts, seasonal=F, ic= "aic")## Series: Datos.ts
## ARIMA(1,1,0) with drift
##
## Coefficients:
## ar1 drift
## 0.5513 0.2995
## s.e. 0.0806 0.0412
##
## sigma^2 = 0.03803: log likelihood = 23.99
## AIC=-41.98 AICc=-41.75 BIC=-33.96modelo8<- auto.arima(Datos.ts, seasonal=F, ic= "aic")
summary(modelo8)## Series: Datos.ts
## ARIMA(1,1,0) with drift
##
## Coefficients:
## ar1 drift
## 0.5513 0.2995
## s.e. 0.0806 0.0412
##
## sigma^2 = 0.03803: log likelihood = 23.99
## AIC=-41.98 AICc=-41.75 BIC=-33.96
##
## Training set error measures:
## ME RMSE MAE MPE MAPE
## Training set -0.0004037835 0.1922967 0.1467918 -0.001079885 0.1544726
## MASE ACF1
## Training set 0.04128232 0.06161327auto.arima(Datos.ts, seasonal=F, ic= "aic", trace= T)##
## ARIMA(2,1,2) with drift : -38.38836
## ARIMA(0,1,0) with drift : -5.528668
## ARIMA(1,1,0) with drift : -41.98022
## ARIMA(0,1,1) with drift : -35.93621
## ARIMA(0,1,0) : 95.52941
## ARIMA(2,1,0) with drift : -41.20685
## ARIMA(1,1,1) with drift : -40.86237
## ARIMA(2,1,1) with drift : -40.40483
## ARIMA(1,1,0) : -24.71434
##
## Best model: ARIMA(1,1,0) with drift## Series: Datos.ts
## ARIMA(1,1,0) with drift
##
## Coefficients:
## ar1 drift
## 0.5513 0.2995
## s.e. 0.0806 0.0412
##
## sigma^2 = 0.03803: log likelihood = 23.99
## AIC=-41.98 AICc=-41.75 BIC=-33.96modelo9 <- auto.arima(Datos.ts, seasonal=F, ic= "aic", trace= T)##
## ARIMA(2,1,2) with drift : -38.38836
## ARIMA(0,1,0) with drift : -5.528668
## ARIMA(1,1,0) with drift : -41.98022
## ARIMA(0,1,1) with drift : -35.93621
## ARIMA(0,1,0) : 95.52941
## ARIMA(2,1,0) with drift : -41.20685
## ARIMA(1,1,1) with drift : -40.86237
## ARIMA(2,1,1) with drift : -40.40483
## ARIMA(1,1,0) : -24.71434
##
## Best model: ARIMA(1,1,0) with drift summary(modelo9)## Series: Datos.ts
## ARIMA(1,1,0) with drift
##
## Coefficients:
## ar1 drift
## 0.5513 0.2995
## s.e. 0.0806 0.0412
##
## sigma^2 = 0.03803: log likelihood = 23.99
## AIC=-41.98 AICc=-41.75 BIC=-33.96
##
## Training set error measures:
## ME RMSE MAE MPE MAPE
## Training set -0.0004037835 0.1922967 0.1467918 -0.001079885 0.1544726
## MASE ACF1
## Training set 0.04128232 0.06161327auto.arima(Datos.ts, seasonal=T, ic= "aic", trace= T)##
## ARIMA(2,1,2)(1,1,1)[12] : -70.09927
## ARIMA(0,1,0)(0,1,0)[12] : -16.44564
## ARIMA(1,1,0)(1,1,0)[12] : -64.56279
## ARIMA(0,1,1)(0,1,1)[12] : -62.59924
## ARIMA(2,1,2)(0,1,1)[12] : -71.31988
## ARIMA(2,1,2)(0,1,0)[12] : -49.12081
## ARIMA(2,1,2)(0,1,2)[12] : -70.2677
## ARIMA(2,1,2)(1,1,0)[12] : -61.68952
## ARIMA(2,1,2)(1,1,2)[12] : -68.51692
## ARIMA(1,1,2)(0,1,1)[12] : -73.18225
## ARIMA(1,1,2)(0,1,0)[12] : -51.08889
## ARIMA(1,1,2)(1,1,1)[12] : -71.97781
## ARIMA(1,1,2)(0,1,2)[12] : -72.14196
## ARIMA(1,1,2)(1,1,0)[12] : -63.66549
## ARIMA(1,1,2)(1,1,2)[12] : -70.36718
## ARIMA(0,1,2)(0,1,1)[12] : -66.71436
## ARIMA(1,1,1)(0,1,1)[12] : -73.27803
## ARIMA(1,1,1)(0,1,0)[12] : -52.18385
## ARIMA(1,1,1)(1,1,1)[12] : -72.13425
## ARIMA(1,1,1)(0,1,2)[12] : -72.37374
## ARIMA(1,1,1)(1,1,0)[12] : -63.9505
## ARIMA(1,1,1)(1,1,2)[12] : -70.61674
## ARIMA(1,1,0)(0,1,1)[12] : -74.1525
## ARIMA(1,1,0)(0,1,0)[12] : -54.05119
## ARIMA(1,1,0)(1,1,1)[12] : Inf
## ARIMA(1,1,0)(0,1,2)[12] : Inf
## ARIMA(1,1,0)(1,1,2)[12] : -72.17038
## ARIMA(0,1,0)(0,1,1)[12] : Inf
## ARIMA(2,1,0)(0,1,1)[12] : -72.88837
## ARIMA(2,1,1)(0,1,1)[12] : -73.15594
##
## Best model: ARIMA(1,1,0)(0,1,1)[12]## Series: Datos.ts
## ARIMA(1,1,0)(0,1,1)[12]
##
## Coefficients:
## ar1 sma1
## 0.6003 -0.6988
## s.e. 0.0834 0.1317
##
## sigma^2 = 0.02438: log likelihood = 40.08
## AIC=-74.15 AICc=-73.89 BIC=-66.49modelo10 <- auto.arima(Datos.ts, seasonal=T, ic= "aic", trace= T)##
## ARIMA(2,1,2)(1,1,1)[12] : -70.09927
## ARIMA(0,1,0)(0,1,0)[12] : -16.44564
## ARIMA(1,1,0)(1,1,0)[12] : -64.56279
## ARIMA(0,1,1)(0,1,1)[12] : -62.59924
## ARIMA(2,1,2)(0,1,1)[12] : -71.31988
## ARIMA(2,1,2)(0,1,0)[12] : -49.12081
## ARIMA(2,1,2)(0,1,2)[12] : -70.2677
## ARIMA(2,1,2)(1,1,0)[12] : -61.68952
## ARIMA(2,1,2)(1,1,2)[12] : -68.51692
## ARIMA(1,1,2)(0,1,1)[12] : -73.18225
## ARIMA(1,1,2)(0,1,0)[12] : -51.08889
## ARIMA(1,1,2)(1,1,1)[12] : -71.97781
## ARIMA(1,1,2)(0,1,2)[12] : -72.14196
## ARIMA(1,1,2)(1,1,0)[12] : -63.66549
## ARIMA(1,1,2)(1,1,2)[12] : -70.36718
## ARIMA(0,1,2)(0,1,1)[12] : -66.71436
## ARIMA(1,1,1)(0,1,1)[12] : -73.27803
## ARIMA(1,1,1)(0,1,0)[12] : -52.18385
## ARIMA(1,1,1)(1,1,1)[12] : -72.13425
## ARIMA(1,1,1)(0,1,2)[12] : -72.37374
## ARIMA(1,1,1)(1,1,0)[12] : -63.9505
## ARIMA(1,1,1)(1,1,2)[12] : -70.61674
## ARIMA(1,1,0)(0,1,1)[12] : -74.1525
## ARIMA(1,1,0)(0,1,0)[12] : -54.05119
## ARIMA(1,1,0)(1,1,1)[12] : Inf
## ARIMA(1,1,0)(0,1,2)[12] : Inf
## ARIMA(1,1,0)(1,1,2)[12] : -72.17038
## ARIMA(0,1,0)(0,1,1)[12] : Inf
## ARIMA(2,1,0)(0,1,1)[12] : -72.88837
## ARIMA(2,1,1)(0,1,1)[12] : -73.15594
##
## Best model: ARIMA(1,1,0)(0,1,1)[12]summary(modelo10)## Series: Datos.ts
## ARIMA(1,1,0)(0,1,1)[12]
##
## Coefficients:
## ar1 sma1
## 0.6003 -0.6988
## s.e. 0.0834 0.1317
##
## sigma^2 = 0.02438: log likelihood = 40.08
## AIC=-74.15 AICc=-73.89 BIC=-66.49
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.01410728 0.1448901 0.1031469 0.01616982 0.1071299 0.02900803
## ACF1
## Training set -0.0473832Pronostico mejor modelo
pronostico2 <- forecast::forecast(modelo9, h=12, level=95)
pronostico2## Point Forecast Lo 95 Hi 95
## Jan 10 110.8432 110.4610 111.2255
## Feb 10 111.2215 110.5160 111.9270
## Mar 10 111.5644 110.5641 112.5647
## Apr 10 111.8878 110.6235 113.1522
## May 10 112.2005 110.6998 113.7013
## Jun 10 112.5073 110.7933 114.2214
## Jul 10 112.8108 110.9025 114.7191
## Aug 10 113.1125 111.0257 115.1994
## Sep 10 113.4132 111.1607 115.6658
## Oct 10 113.7134 111.3059 116.1209
## Nov 10 114.0133 111.4599 116.5667
## Dec 10 114.3130 111.6214 117.0045autoplot(pronostico2)
Pruebas de normalidad correspondientes
library(moments)
library(nortest)
qqnorm(Datos.ts)
qqline(Datos.ts)
hist(Datos.ts)
plot((density(Datos.ts)), ylab='Densidad', col= 'blue3', xlab = 'x25',
las=1, lwd=4)
shapiro.test(Datos.ts)##
## Shapiro-Wilk normality test
##
## data: Datos.ts
## W = 0.92576, p-value = 1.442e-05