Taller 2
Juan José Tobón
##Taller 2: En este trabajo se interpretarán los resultados de las tendencias, ciclos, componente estacional, descomposición de los modelos, comparación entre series, se harán pronosticos y se presentarán conclusiones de las siguientes variables: Precio del petróleo (brent), Precio de las acciones de Ecopetrol (Pecop), Exportaciones(Xv) e Importaciones(Xcol).
#Cargar base
library(readxl)
Datos_Taller2 <- read_excel("C:/Users/Usuario/Desktop/Prediccion y big data/Datos Taller2.xlsx")## New names:
## • `` -> `...1`library(readxl) Datos_Taller2 <- read_excel(“C:/Users/Usuario/Desktop/Prediccion y big data/Datos Taller2.xlsx”)
#Cargar paquetes
#install.packages("lubridate")
#install.packages("esquisse")
#library(lubridate)
#install.packages("tidyverse")
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.2 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.2 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(ggplot2)
library(esquisse)
#install.packages("forecast")
library(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo#install.packages("tseries")
library(tseries)
#install.packages("TSA")
library(TSA) ## Registered S3 methods overwritten by 'TSA':
## method from
## fitted.Arima forecast
## plot.Arima forecast
##
## Attaching package: 'TSA'
##
## The following object is masked from 'package:readr':
##
## spec
##
## The following objects are masked from 'package:stats':
##
## acf, arima
##
## The following object is masked from 'package:utils':
##
## tar#install.packages("urca")
library(urca)
library(ggplot2)
#install.packages("dplyr")
library(dplyr)
library(stats)
#install.packages("seasonal")
library(seasonal)##
## Attaching package: 'seasonal'
##
## The following object is masked from 'package:tibble':
##
## view#install.packages("stats")
#install.packages("decompose")
#install.packages("fpp2")
library(fpp2)## ── Attaching packages ────────────────────────────────────────────── fpp2 2.5 ──
## ✔ fma 2.5 ✔ expsmooth 2.3install.packages(“lubridate”) install.packages(“esquisse”)
library(lubridate) install.packages(“tidyverse”) library(tidyverse)
library(ggplot2) library(esquisse) install.packages(“forecast”)
library(forecast) install.packages(“tseries”) library(tseries)
install.packages(“TSA”) library(TSA) install.packages(“urca”)
library(urca)
library(ggplot2) install.packages(“dplyr”) library(dplyr)
library(stats)
install.packages(“seasonal”) library(seasonal) install.packages(“stats”)
install.packages(“decompose”) install.packages(“fpp2”) library(fpp2)
#Estadisticas descriptivas ESTADISTICAS DESCRIPTIVAS
head(Datos_Taller2)## # A tibble: 6 × 21
## ...1 `-` td ohotel licc dah vehc ipcc enerv xv
## <dttm> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <dbl>
## 1 2012-01-01 00:00:00 6.94e7 15.1 41.0 71805 4.07e12 NA 77.4 200. 1.60e8
## 2 2012-02-01 00:00:00 6.83e7 15.8 48.2 60566 4.08e12 NA 77.6 201. 1.80e8
## 3 2012-03-01 00:00:00 6.96e7 15.8 52.8 97722 3.94e12 NA 77.5 221 1.91e8
## 4 2012-04-01 00:00:00 6.70e7 15.5 46.2 95605 3.88e12 NA 77.5 203. 1.61e8
## 5 2012-05-01 00:00:00 6.84e7 14.9 51.7 54599 3.98e12 NA 78.1 221. 1.99e8
## 6 2012-06-01 00:00:00 6.87e7 14.6 50.5 54973 3.94e12 NA 78.0 217. 1.72e8
## # ℹ 11 more variables: can <dbl>, ipir <dbl>, vehv <dbl>, enercol <dbl>,
## # res <dbl>, ipc <dbl>, vehcol <chr>, xcol <dbl>, brent <dbl>, pecop <dbl>,
## # ise <dbl>summary(Datos_Taller2)## ...1 - td
## Min. :2012-01-01 00:00:00.00 Min. :49201417 Min. : 9.388
## 1st Qu.:2014-09-23 12:00:00.00 1st Qu.:68424506 1st Qu.:12.055
## Median :2017-06-16 00:00:00.00 Median :73022796 Median :13.272
## Mean :2017-06-16 08:43:38.18 Mean :72314074 Mean :14.023
## 3rd Qu.:2020-03-08 18:00:00.00 3rd Qu.:76688026 3rd Qu.:14.365
## Max. :2022-12-01 00:00:00.00 Max. :84518504 Max. :30.128
##
## ohotel licc dah vehc
## Min. : 3.326 Min. : 11855 Min. :3.878e+12 Length:132
## 1st Qu.:48.354 1st Qu.: 50811 1st Qu.:5.563e+12 Class :character
## Median :54.751 Median : 75687 Median :6.109e+12 Mode :character
## Mean :51.017 Mean : 86611 Mean :7.673e+12
## 3rd Qu.:59.304 3rd Qu.:109978 3rd Qu.:1.067e+13
## Max. :70.739 Max. :292808 Max. :1.347e+13
## NA's :1
## ipcc enerv xv can
## Min. : 77.45 Min. :149.2 Min. : 49767731 Min. : 138148
## 1st Qu.: 82.55 1st Qu.:222.0 1st Qu.:141128908 1st Qu.:1716445
## Median : 96.39 Median :229.7 Median :156797911 Median :2002946
## Mean : 95.76 Mean :228.8 Mean :157104611 Mean :1941597
## 3rd Qu.:105.94 3rd Qu.:236.9 3rd Qu.:171067180 3rd Qu.:2213750
## Max. :128.07 Max. :259.4 Max. :222055829 Max. :2493511
## NA's :3
## ipir vehv enercol res
## Min. : 37.79 Min. : 4 Min. :1371 Min. :32518
## 1st Qu.: 92.62 1st Qu.:2178 1st Qu.:1713 1st Qu.:46705
## Median : 98.88 Median :2411 Median :1771 Median :47402
## Mean : 99.16 Mean :2365 Mean :1796 Mean :48610
## 3rd Qu.:104.44 3rd Qu.:2604 3rd Qu.:1851 3rd Qu.:53723
## Max. :127.02 Max. :3816 Max. :2179 Max. :59144
##
## ipc vehcol xcol brent
## Min. : 76.75 Length:132 Min. :1894291 Min. : 23.34
## 1st Qu.: 82.11 Class :character 1st Qu.:2900069 1st Qu.: 53.80
## Median : 96.20 Mode :character Median :3279388 Median : 66.43
## Mean : 95.21 Mean :3513579 Mean : 74.39
## 3rd Qu.:104.96 3rd Qu.:4152421 3rd Qu.:103.05
## Max. :126.03 Max. :5559477 Max. :124.93
##
## pecop ise
## Min. : 6.44 Min. : 81.92
## 1st Qu.:10.44 1st Qu.: 95.55
## Median :16.20 Median :102.00
## Mean :22.37 Mean :103.03
## 3rd Qu.:28.02 3rd Qu.:109.03
## Max. :64.70 Max. :133.63
## #plot(Datos_Taller2
#data <- Datos_Taller2[1:180,]
#test <- Datos_Taller2[181:188,] head(Datos_Taller2) summary(Datos_Taller2) plot(Datos_Taller2)
data <- Datos_Taller2[1:180,] test <- Datos_Taller2[181:188,]
#Crecación de serie de tiempo historica data <- ts(Datos_Taller2$pecop, frequency=12,start=2012)
data <- ts(Datos_Taller2$pecop, frequency=12,start=2012)#creación de la serie de tiempo para validación test <- ts(test$pecop, frequency=12,start=2012)
#test <- ts(test$pecop, frequency=12,start=2012)#Grafico de datos 1 autoplot(data)+xlab(“Tiempo”)+ylab(“Dólares”)+ggtitle(“Precio de las acciones de Ecopetrol”)
autoplot(data)+xlab("Tiempo")+ylab("Dólares")+ggtitle("Precio de las acciones de Ecopetrol")
#Grafico de datos 2 ggseasonplot(data)+ggtitle(“Demanda agrupada por estaciones”) ggsubseriesplot(data)+ggtitle(“Demanda por estación”) autoplot(stl(data, s.window=7))
ggseasonplot(data)+ggtitle("Demanda agrupada por estaciones")
ggsubseriesplot(data)+ggtitle("Demanda por estación")
autoplot(stl(data, s.window=7))
#Cargar datos de variables pecop<-ts(Datos_Taller2\(pecop, frequency=12, start=c(2012,1)) brent<-ts(Datos_Taller2\)brent[1:96], frequency=12, start=c(2012,1)) xv<-ts(Datos_Taller2\(xv, frequency=12, start=c(2012,1)) xcol<-ts(Datos_Taller2\)xcol[1:96], frequency=12, start=c(2012,1))
pecop<-ts(Datos_Taller2$pecop, frequency=12, start=c(2012,1))
brent<-ts(Datos_Taller2$brent[1:96], frequency=12, start=c(2012,1))
xv<-ts(Datos_Taller2$xv, frequency=12, start=c(2012,1))
xcol<-ts(Datos_Taller2$xcol[1:96], frequency=12, start=c(2012,1))#Utilizamos la función decompose (del paquete cargado previamente “STATS”) pecop_decomp<- decompose(pecop) brent_decomp <- decompose(brent) xv_decomp <- decompose(xv) xcol_decomp <- decompose(xcol)
pecop_decomp<- decompose(pecop)
brent_decomp <- decompose(brent)
xv_decomp <- decompose(xv)
xcol_decomp <- decompose(xcol)Graficar los componentes
#Precio de Acciones de Ecopetrol Se descomponen la serie de tiempo de “Pecop” desde su serie original, luego marcando su tendencia, porteriormente la estacionalidad y por ultimo mostrando el factor irregular.
par(mfrow = c(2, 2)) #Se utiliza para dividir la ventana gráfica en una matriz de 2 filas y 2 columnas plot(pecop_decomp\(x, main = "Precio acciones Ecopetrol", col = "grey", ylab = "Serie de tiempo") plot(pecop_decomp\)trend, main = “Tendencia”, col = “orange”, ylab = “Valores”) plot(pecop_decomp\(seasonal, main = "Estacionalidad", col = "yellow", ylab = "Valores") plot(pecop_decomp\)random, main = “Irregularidad”, col = “green”, ylab = “Valores”)
par(mfrow = c(2, 2)) #Se utiliza para dividir la ventana gráfica en una matriz de 2 filas y 2 columnas
plot(pecop_decomp$x, main = "Precio acciones Ecopetrol", col = "grey", ylab = "Serie de tiempo")
plot(pecop_decomp$trend, main = "Tendencia", col = "orange", ylab = "Valores")
plot(pecop_decomp$seasonal, main = "Estacionalidad", col = "yellow", ylab = "Valores")
plot(pecop_decomp$random, main = "Irregularidad", col = "green", ylab = "Valores")
#Precio del Petróleo Brent Se descomponen la serie de tiempo de “Brent” desde su serie original, luego marcando su tendencia, porteriormente la estacionalidad y por ultimo mostrando el factor irregular.
par(mfrow = c(2, 2)) #Se utiliza para dividir la ventana gráfica en una matriz de 2 filas y 2 columnas plot(brent_decomp\(x, main = "Precio del Petróleo Brent-Original", col = "red", ylab = "Serie de tiempo") plot(brent_decomp\)trend, main = “Tendencia”, col = “yellow”, ylab = “Valores”) plot(brent_decomp\(seasonal, main = "Estacionalidad", col = "pink", ylab = "Valores") plot(brent_decomp\)random, main = “Irregularidad”, col = “blue”, ylab = “Valores”)
par(mfrow = c(2, 2)) #Se utiliza para dividir la ventana gráfica en una matriz de 2 filas y 2 columnas
plot(brent_decomp$x, main = "Precio del Petróleo Brent-Original", col = "red", ylab = "Serie de tiempo")
plot(brent_decomp$trend, main = "Tendencia", col = "yellow", ylab = "Valores")
plot(brent_decomp$seasonal, main = "Estacionalidad", col = "pink", ylab = "Valores")
plot(brent_decomp$random, main = "Irregularidad", col = "blue", ylab = "Valores")
#Exportaciones
Se descomponen la serie de tiempo de “Xv” desde su serie original, luego marcando su tendencia, porteriormente la estacionalidad y por ultimo mostrando el factor irregular.
par(mfrow = c(2, 2)) #Se utiliza para dividir la ventana gráfica en una matriz de 2 filas y 2 columnas plot(xv_decomp\(x, main = "Exportaciones-Original", col = "red", ylab = "Serie de tiempo") plot(xv_decomp\)trend, main = “Tendencia”, col = “yellow”, ylab = “Valores”) plot(xv_decomp\(seasonal, main = "Estacionalidad", col = "pink", ylab = "Valores") plot(xv_decomp\)random, main = “Irregularidad”, col = “blue”, ylab = “Valores”)
par(mfrow = c(2, 2)) #Se utiliza para dividir la ventana gráfica en una matriz de 2 filas y 2 columnas
plot(xv_decomp$x, main = "Exportaciones-Original", col = "red", ylab = "Serie de tiempo")
plot(xv_decomp$trend, main = "Tendencia", col = "yellow", ylab = "Valores")
plot(xv_decomp$seasonal, main = "Estacionalidad", col = "pink", ylab = "Valores")
plot(xv_decomp$random, main = "Irregularidad", col = "blue", ylab = "Valores")
#Importaciones Se descomponen la serie de tiempo de “Xcol” desde su serie original, luego marcando su tendencia, porteriormente la estacionalidad y por ultimo mostrando el factor irregular.
par(mfrow = c(2, 2)) #Se utiliza para dividir la ventana gráfica en una matriz de 2 filas y 2 columnas plot(xcol_decomp\(x, main = "Importaciones-Original", col = "red", ylab = "Serie de tiempo") plot(xcol_decomp\)trend, main = “Tendencia”, col = “yellow”, ylab = “Valores”) plot(xcol_decomp\(seasonal, main = "Estacionalidad", col = "pink", ylab = "Valores") plot(xcol_decomp\)random, main = “Irregularidad”, col = “blue”, ylab = “Valores”)
par(mfrow = c(2, 2)) #Se utiliza para dividir la ventana gráfica en una matriz de 2 filas y 2 columnas
plot(xcol_decomp$x, main = "Importaciones-Original", col = "red", ylab = "Serie de tiempo")
plot(xcol_decomp$trend, main = "Tendencia", col = "yellow", ylab = "Valores")
plot(xcol_decomp$seasonal, main = "Estacionalidad", col = "pink", ylab = "Valores")
plot(xcol_decomp$random, main = "Irregularidad", col = "blue", ylab = "Valores")
#Aplicación del modelo ARIMA con Brent
Se va a utilizar el modelo ARIMA para la predicción de la variable “Brent” por medio de los parametros p, d y q.
Realizar prueba de raíz unitaria
adf.test(brent)
adf.test(brent)##
## Augmented Dickey-Fuller Test
##
## data: brent
## Dickey-Fuller = -1.1392, Lag order = 4, p-value = 0.9122
## alternative hypothesis: stationary#Para hacer este proceso se usa la libreria seasonal instalada y cargada previamente. #Colocamos el nombre a la nueva variable: brent_SA (serie ajustada por estacionalidad) brent_SA <- seasadj(brent_decomp)
brent_SA <- seasadj(brent_decomp)Graficar la serie de tiempo original y desestacionalizada en un mismo gráfico
Se muestra como la grafica original y desestacionalizada están una encima de la otra denotando diferencias en algunos periodos especificos, sin embargo, mantienen una tendencia similar.
plot(brent, main = “Figura 3. BRENT original y desestacionalizada”) lines(brent_SA, col = “red”) legend(“topleft”, legend = c(“Serie original”, “Serie desestacionalizada”), col = c(“black”, “red”), lty = 1)
plot(brent, main = "Figura 3. BRENT original y desestacionalizada")
lines(brent_SA, col = "red")
legend("topleft", legend = c("Serie original", "Serie desestacionalizada"), col = c("black", "red"), lty = 1)
#Para observar el comportamiento de la variable Brent Por medio de la diferenciación en niveles podemos nostar que los datos presentan tendencias y varios cambios significativos durante la serie de tiempo, los cuales no se pueden apreciar simplemente por medio de la serie original.
brentSA_d1= diff(brent, differences = 1) brentSA_d2= diff(brent, differences = 2)
plot(brentSA_d1, main = “Figura 4. Brent- Diferenciación en (niveles)”) lines(brentSA_d2, col = “red”) legend(“topleft”, legend = c(” Una diferencia”, “Dos diferencias”), col = c(“black”, “red”), lty = 1)
brentSA_d1= diff(brent, differences = 1)
brentSA_d2= diff(brent, differences = 2)
plot(brentSA_d1, main = "Figura 4. Brent- Diferenciación en (niveles)")
lines(brentSA_d2, col = "red")
legend("topleft", legend = c(" Una diferencia", "Dos diferencias"), col = c("black", "red"), lty = 1)
#Hacemos el Test Dickey-Fuller para corroborar: Por medio de este Test se deretminará si la serie es estacionaria o no.
H0: No estacionaria H1: Estacionaria
#H0: No estacionaria H1: EstacionariaRealizar prueba de raíz unitaria
adf.test(brentSA_d1) adf.test(brentSA_d2)
adf.test(brentSA_d1)## Warning in adf.test(brentSA_d1): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: brentSA_d1
## Dickey-Fuller = -4.1706, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationaryadf.test(brentSA_d2)## Warning in adf.test(brentSA_d2): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: brentSA_d2
## Dickey-Fuller = -7.0527, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationaryEl resultado del Test da a entender que la variable Brent es estacionaria, gracias al lo que indican los valores p=0.01 lo que indica que se rechaza la hipotesis por un valor signifacante de 1%.
#Determinamos p y q Acf(brentSA_d1, main=‘Figura 5.Función de autocorrelación (ACF) -brent diferenciado’) Pacf(brentSA_d1, main=‘Figura 6.Función de autocorrelación parcial (PACF) -brent diferenciado’)
Acf(brentSA_d1, main='Figura 5.Función de autocorrelación (ACF) -brent diferenciado')
Pacf(brentSA_d1, main='Figura 6.Función de autocorrelación parcial (PACF) -brent diferenciado')
#ESTIMACIÓN DE LOS POSIBLES MODELOS #Corremos la función auto.arima: auto.arima(brent_SA)
auto.arima(brent_SA)## Series: brent_SA
## ARIMA(0,1,1)
##
## Coefficients:
## ma1
## 0.3721
## s.e. 0.0890
##
## sigma^2 = 23.55: log likelihood = -284.43
## AIC=572.86 AICc=572.99 BIC=577.97#VALIDACIÓN
Ahora, para la validación del resultado, se hará un Test Ljung-Box que sirve para probar la autocorrelación de la serie de tiempo, lo que sirve para identificar un pronóstico mucho más acertado.
Model2=Arima(brent_SA,order = c(0,1,2))
Box.test(Model2$residuals, lag = 20, type = “Ljung-Box”)
Model2=Arima(brent_SA,order = c(0,1,2))
Box.test(Model2$residuals, lag = 20, type = "Ljung-Box")##
## Box-Ljung test
##
## data: Model2$residuals
## X-squared = 10.156, df = 20, p-value = 0.9653El resultado del Test indica que no existe suficiente evidencia para una autocorrelación gracias al valor p=0.9653.
#PRONÓSTICO
#Se hace pronpostico a 3 y 06 meses. Se elije el parametro c(95) para que sea a un nivel de confianza del 95%
Pronostico3=forecast(Model2,level= c(95), h=3) plot(Pronostico3) Pronostico3 Pronostico6=forecast(Model2,level= c(95), h=6) plot(Pronostico6) Pronostico6
Pronostico3=forecast(Model2,level= c(95), h=3)
plot(Pronostico3)
Pronostico3## Point Forecast Lo 95 Hi 95
## Jan 2020 70.60169 61.04572 80.15767
## Feb 2020 70.69839 54.30256 87.09423
## Mar 2020 70.69839 49.31457 92.08222Pronostico6=forecast(Model2,level= c(95), h=6)
plot(Pronostico6)
Pronostico6## Point Forecast Lo 95 Hi 95
## Jan 2020 70.60169 61.04572 80.15767
## Feb 2020 70.69839 54.30256 87.09423
## Mar 2020 70.69839 49.31457 92.08222
## Apr 2020 70.69839 45.28752 96.10926
## May 2020 70.69839 41.81662 99.58017
## Jun 2020 70.69839 38.72026 102.67653Ambos pronosticos dan resultados estables durante el tiempo, lo que permite una toma de decisión en el tiempo que se nos muestra y con los datos que se manejan en estos periodos.En resumidas cuentas, se pronostica de 3 a 6 meses, el precio del Petróleo se mantenga estable.
#Aplicación del modelo ARIMA con Pecop A continuación, se aplicará el modelo ARIMA para la predicción a la variable Pecop con la misma metodología de la variable Brent.
Realizar prueba de raíz unitaria
adf.test(pecop)
adf.test(pecop)##
## Augmented Dickey-Fuller Test
##
## data: pecop
## Dickey-Fuller = -1.6506, Lag order = 5, p-value = 0.7218
## alternative hypothesis: stationaryEl resultado del Test Dickey-Fuller indica que la variable no es estacionaria, gracias al valor p=0.7218 lo que descarta la hipotesis.
#Para hacer este proceso se usa la libreria seasonal instalada y cargada previamente. #Colocamos el nombre a la nueva variable: pecop_SA (serie ajustada por estacionalidad) pecop_SA <- seasadj(pecop_decomp)
pecop_SA <- seasadj(pecop_decomp)Graficar la serie de tiempo original y desestacionalizada en un mismo gráfico
plot(pecop, main = “Figura 3. PECOP original y desestacionalizada”) lines(pecop_SA, col = “red”) legend(“topleft”, legend = c(“Serie original”, “Serie desestacionalizada”), col = c(“black”, “red”), lty = 1)
plot(pecop, main = "Figura 3. PECOP original y desestacionalizada")
lines(pecop_SA, col = "red")
legend("topleft", legend = c("Serie original", "Serie desestacionalizada"), col = c("black", "red"), lty = 1)
Se puede observar la serie Original sobre la desestacionalizada y demuestra que la serie desestacionalizada se superpone a la original en casi todo el periodo.
#Para observar el comportamiento de la variable Pecop
pecopSA_d1= diff(pecop, differences = 1) pecopSA_d2= diff(pecop, differences = 2)
pecopSA_d1= diff(pecop, differences = 1)
pecopSA_d2= diff(pecop, differences = 2)Por medio de la diferenciación se puede interpretar las diferentes tendencias y patrones que contienen los datos.
plot(pecopSA_d1, main = “Figura 4. Peco- Diferenciación en (niveles)”) lines(pecopSA_d2, col = “red”) legend(“topleft”, legend = c(” Una diferencia”, “Dos diferencias”), col = c(“black”, “red”), lty = 1)
plot(pecopSA_d1, main = "Figura 4. Peco- Diferenciación en (niveles)")
lines(pecopSA_d2, col = "red")
legend("topleft", legend = c(" Una diferencia", "Dos diferencias"), col = c("black", "red"), lty = 1)
#Hacemos el Test Dickey-Fuller para corroborar:
H0: No estacionaria H1: Estacionaria
#H0: No estacionaria H1: EstacionariaRealizar prueba de raíz unitaria
adf.test(pecopSA_d1) adf.test(pecopSA_d2)
adf.test(pecopSA_d1)## Warning in adf.test(pecopSA_d1): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: pecopSA_d1
## Dickey-Fuller = -4.5208, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryadf.test(pecopSA_d2)## Warning in adf.test(pecopSA_d2): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: pecopSA_d2
## Dickey-Fuller = -7.6788, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryDe nuevo se corroboran los datos por medio del Test Dickey-Fuller, lo que arroja un resultado interesante al comprobar por medio de la serie estacionaria que la variable en efecto contiene una mayor evidencia de ser estacionaria gracias a su valor p=0.01 menor que 0.05, lo que indica que la evidencia es menor que el valor impreso, dando a entender que en efecto, es estacionaria.
#Determinamos p y q Acf(pecopSA_d1, main=‘Figura 5.Función de autocorrelación (ACF) -brent diferenciado’) Pacf(pecopSA_d1, main=‘Figura 6.Función de autocorrelación parcial (PACF) -brent diferenciado’)
Acf(pecopSA_d1, main='Figura 5.Función de autocorrelación (ACF) -brent diferenciado')
Pacf(pecopSA_d1, main='Figura 6.Función de autocorrelación parcial (PACF) -brent diferenciado')
#ESTIMACIÓN DE LOS POSIBLES MODELOS #Corremos la función auto.arima: auto.arima(pecop_SA)
auto.arima(pecop_SA)## Series: pecop_SA
## ARIMA(1,1,0)(0,0,2)[12] with drift
##
## Coefficients:
## ar1 sma1 sma2 drift
## 0.2216 -0.3411 -0.3252 -0.2754
## s.e. 0.0876 0.1252 0.1357 0.1162
##
## sigma^2 = 4.973: log likelihood = -292.05
## AIC=594.1 AICc=594.58 BIC=608.47#VALIDACIÓN Model2=Arima(pecop_SA,order = c(0,1,2))
Model2=Arima(pecop_SA,order = c(0,1,2))El resultado del modelo ARIMA permite entender que es ajustado y adecuado para modelar la serie temporal.
Box.test(Model2$residuals, lag = 20, type = “Ljung-Box”)
Box.test(Model2$residuals, lag = 20, type = "Ljung-Box")##
## Box-Ljung test
##
## data: Model2$residuals
## X-squared = 24.653, df = 20, p-value = 0.215El resultado de la prueba Ljung-Box permite ver el ruido blanco en los datos, lo que da a entender por medio del mismo que el modelo parece ser el correcto gracias a que los datos residuales no son los suficientmente relevantes para realizar una predicción correcta.
#PRONÓSTICO
#Se hace pronpostico a 3 y 06 meses. Se elije el parametro c(95) para que sea a un nivel de confianza del 95%
Pronostico3=forecast(Model2,level= c(95), h=3) plot(Pronostico3) Pronostico3 Pronostico6=forecast(Model2,level= c(95), h=6) plot(Pronostico6) Pronostico6
Pronostico3=forecast(Model2,level= c(95), h=3)
plot(Pronostico3)
Pronostico3## Point Forecast Lo 95 Hi 95
## Jan 2023 11.39595 6.767594 16.02430
## Feb 2023 11.37828 4.121806 18.63475
## Mar 2023 11.37828 1.700243 21.05631Pronostico6=forecast(Model2,level= c(95), h=6)
plot(Pronostico6)
Pronostico6## Point Forecast Lo 95 Hi 95
## Jan 2023 11.39595 6.7675942 16.02430
## Feb 2023 11.37828 4.1218063 18.63475
## Mar 2023 11.37828 1.7002428 21.05631
## Apr 2023 11.37828 -0.2265636 22.98312
## May 2023 11.37828 -1.8761677 24.63272
## Jun 2023 11.37828 -3.3420586 26.09861Se pronostica un aumento gradual en el Precio de las Acciones de Ecopetrol durante el periodo evaluado, lo que tiene un intervalo de confianza del 95%, esto permite una toma de decisiones certera sobre las acciones que se pueden tomar respecto al Precio de la Acción de Ecopetrol.
#Aplicación del modelo ARIMA con Xcol
A continuación, se aplicará el modelo ARIMA para la predicción a la variable XCol con la misma metodología de la variable anterior.
Realizar prueba de raíz unitaria
adf.test(xcol)
adf.test(xcol)##
## Augmented Dickey-Fuller Test
##
## data: xcol
## Dickey-Fuller = -1.4401, Lag order = 4, p-value = 0.8078
## alternative hypothesis: stationaryEl resultado del Test Dickey-Fuller indica que la variable no contiene suficiente evidencia para ser clasificada como no estacionaria gracias a su valo p=0.8078.
#Para hacer este proceso se usa la libreria seasonal instalada y cargada previamente. #Colocamos el nombre a la nueva variable: xcol_SA (serie ajustada por estacionalidad) xcol_SA <- seasadj(xcol_decomp)
xcol_SA <- seasadj(xcol_decomp)Graficar la serie de tiempo original y desestacionalizada en un mismo gráfico
plot(xcol, main = “Figura 3. XCOL original y desestacionalizada”) lines(xcol_SA, col = “red”) legend(“topleft”, legend = c(“Serie original”, “Serie desestacionalizada”), col = c(“black”, “red”), lty = 1)
plot(xcol, main = "Figura 3. XCOL original y desestacionalizada")
lines(xcol_SA, col = "red")
legend("topleft", legend = c("Serie original", "Serie desestacionalizada"), col = c("black", "red"), lty = 1)
Se puede observar una superpoción de las series original y desestacionalizada en periodos especificos, sin embargo, no marcan una tendencia.
#Para observar el comportamiento de la variable Xcol
xcolSA_d1= diff(pecop, differences = 1) xcolSA_d2= diff(pecop, differences = 2)
plot(xcolSA_d1, main = “Figura 4. Xcol- Diferenciación en (niveles)”) lines(xcolSA_d2, col = “red”) legend(“topleft”, legend = c(” Una diferencia”, “Dos diferencias”), col = c(“black”, “red”), lty = 1)
xcolSA_d1= diff(pecop, differences = 1)
xcolSA_d2= diff(pecop, differences = 2)
plot(xcolSA_d1, main = "Figura 4. Xcol- Diferenciación en (niveles)")
lines(xcolSA_d2, col = "red")
legend("topleft", legend = c(" Una diferencia", "Dos diferencias"), col = c("black", "red"), lty = 1)
La diferenciación en niveles permite observar evidencias de tendencias en el periodo de tiempo.
#Hacemos el Test Dickey-Fuller para corroborar:
H0: No estacionaria H1: Estacionaria
#H0: No estacionaria H1: EstacionariaRealizar prueba de raíz unitaria
adf.test(xcolSA_d1) adf.test(xcolSA_d2)
adf.test(xcolSA_d1)## Warning in adf.test(xcolSA_d1): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: xcolSA_d1
## Dickey-Fuller = -4.5208, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryadf.test(xcolSA_d2)## Warning in adf.test(xcolSA_d2): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: xcolSA_d2
## Dickey-Fuller = -7.6788, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryLos resultados del Test Dickey-Fuller indican que la serie es estacionaria con un alto nivel de confianza, por lo que permite entender que no contiene componentes aleatorios en la misma, esto gracias al valor p=0.01 y al mensaje que indica que el valor p es incluso menor que el impreso, lo que indica un mayor nivel de confianza ante este resultado.
#Determinamos p y q Acf(xvSA_d1, main=‘Figura 5.Función de autocorrelación (ACF) -XV diferenciado’) Pacf(xvSA_d1, main=‘Figura 6.Función de autocorrelación parcial (PACF) -XV diferenciado’)
#Acf(xvSA_d1, main='Figura 5.Función de autocorrelación (ACF) -XV diferenciado')
#Pacf(xvSA_d1, main='Figura 6.Función de autocorrelación parcial (PACF) -XV diferenciado')#ESTIMACIÓN DE LOS POSIBLES MODELOS #Corremos la función auto.arima: auto.arima(xcol_SA)
auto.arima(xcol_SA)## Series: xcol_SA
## ARIMA(0,1,1)(0,0,1)[12] with drift
##
## Coefficients:
## ma1 sma1 drift
## -0.4423 -0.3786 -23648.44
## s.e. 0.0829 0.1925 10419.49
##
## sigma^2 = 7.353e+10: log likelihood = -1322.81
## AIC=2653.62 AICc=2654.06 BIC=2663.83#VALIDACIÓN Model2=Arima(xv_SA,order = c(0,1,2))
#Model2=Arima(xv_SA,order = c(0,1,2))Los criterios arrojados por el modelo ARIMA sugieren que los datos son estacionales y contienen un buen ajuste para la predicción.
Box.test(Model2$residuals, lag = 20, type = “Ljung-Box”)
Box.test(Model2$residuals, lag = 20, type = "Ljung-Box")##
## Box-Ljung test
##
## data: Model2$residuals
## X-squared = 24.653, df = 20, p-value = 0.215Por medio del Test Ljung-Box se puede determinar gracias al valor p=0.3235 mayor que 0.05 que los datos residuales no son significantes en el resultado lo que permite entender que es un modelo adecuado.
#PRONÓSTICO
#Se hace pronpostico a 3 y 06 meses. Se elije el parametro c(95) para que sea a un nivel de confianza del 95%
Pronostico3=forecast(Model2,level= c(95), h=3) plot(Pronostico3) Pronostico3 Pronostico6=forecast(Model2,level= c(95), h=6) plot(Pronostico6) Pronostico6
Pronostico3=forecast(Model2,level= c(95), h=3)
plot(Pronostico3)
Pronostico3## Point Forecast Lo 95 Hi 95
## Jan 2023 11.39595 6.767594 16.02430
## Feb 2023 11.37828 4.121806 18.63475
## Mar 2023 11.37828 1.700243 21.05631Pronostico6=forecast(Model2,level= c(95), h=6)
plot(Pronostico6)
Pronostico6## Point Forecast Lo 95 Hi 95
## Jan 2023 11.39595 6.7675942 16.02430
## Feb 2023 11.37828 4.1218063 18.63475
## Mar 2023 11.37828 1.7002428 21.05631
## Apr 2023 11.37828 -0.2265636 22.98312
## May 2023 11.37828 -1.8761677 24.63272
## Jun 2023 11.37828 -3.3420586 26.09861Los resultados del pronostico tienen un intervalo de confianza del 95%, lo que permite tomar decisiones respecto a los valores de importaciones en colombia o por lo menos entender de alguna manera el mercado y su posible comportamiento a futuro.
#Aplicación del modelo ARIMA con Xv Se utilizará el modelo ARIMA para la predicción de la variable Xv.
Realizar prueba de raíz unitaria
adf.test(xv)
adf.test(xv)## Warning in adf.test(xv): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: xv
## Dickey-Fuller = -5.5039, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryEl Test Dickey-Fuller sirve para determinar la estacionalidad de la serie.En este caso, permite comprobar que en efecto, la serie es estacionaria.
#Para hacer este proceso se usa la libreria seasonal instalada y cargada previamente. #Colocamos el nombre a la nueva variable: xv_SA (serie ajustada por estacionalidad) xv_SA <- seasadj(xv_decomp)
xv_SA <- seasadj(xv_decomp)Graficar la serie de tiempo original y desestacionalizada en un mismo gráfico
plot(xv, main = “Figura 3. PECOP original y desestacionalizada”) lines(xv_SA, col = “red”) legend(“topleft”, legend = c(“Serie original”, “Serie desestacionalizada”), col = c(“black”, “red”), lty = 1)
plot(xv, main = "Figura 3. PECOP original y desestacionalizada")
lines(xv_SA, col = "red")
legend("topleft", legend = c("Serie original", "Serie desestacionalizada"), col = c("black", "red"), lty = 1)
Por medio de la superposición de la serie original y la desestacionalizada se permite entender que existen algunos periodos diferentes a la tendencia.
#Para observar el comportamiento de la variable XV
xvSA_d1= diff(pecop, differences = 1) xvSA_d2= diff(pecop, differences = 2)
plot(xvSA_d1, main = “Figura 4. XV- Diferenciación en (niveles)”) lines(xvSA_d2, col = “red”) legend(“topleft”, legend = c(” Una diferencia”, “Dos diferencias”), col = c(“black”, “red”), lty = 1)
xvSA_d1= diff(pecop, differences = 1)
xvSA_d2= diff(pecop, differences = 2)
plot(xvSA_d1, main = "Figura 4. XV- Diferenciación en (niveles)")
lines(xvSA_d2, col = "red")
legend("topleft", legend = c(" Una diferencia", "Dos diferencias"), col = c("black", "red"), lty = 1)
La diferenciación en niveles permite entender los patrones y cambios que tiene la serie de datos.
#Hacemos el Test Dickey-Fuller para corroborar:
H0: No estacionaria H1: Estacionaria
#H0: No estacionaria H1: EstacionariaRealizar prueba de raíz unitaria
adf.test(xvSA_d1) adf.test(xvSA_d2)
adf.test(xvSA_d1)## Warning in adf.test(xvSA_d1): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: xvSA_d1
## Dickey-Fuller = -4.5208, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryadf.test(xvSA_d2)## Warning in adf.test(xvSA_d2): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: xvSA_d2
## Dickey-Fuller = -7.6788, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationarySe aplica Test de Dickey-Fuller a la serie estacional, lo que indica que en efecto la serie estacioranaria gracias a p=0.01 menor que el 5%.
#Determinamos p y q Acf(xvSA_d1, main=‘Figura 5.Función de autocorrelación (ACF) -XV diferenciado’) Pacf(xvSA_d1, main=‘Figura 6.Función de autocorrelación parcial (PACF) -XV diferenciado’)
Acf(xvSA_d1, main='Figura 5.Función de autocorrelación (ACF) -XV diferenciado')
Pacf(xvSA_d1, main='Figura 6.Función de autocorrelación parcial (PACF) -XV diferenciado')
#ESTIMACIÓN DE LOS POSIBLES MODELOS #Corremos la función auto.arima: auto.arima(xv_SA)
auto.arima(xv_SA)## Series: xv_SA
## ARIMA(0,1,2)
##
## Coefficients:
## ma1 ma2
## -0.6100 -0.2670
## s.e. 0.0848 0.0824
##
## sigma^2 = 3.263e+14: log likelihood = -2374.47
## AIC=4754.94 AICc=4755.13 BIC=4763.57El modelo ARIMA(0,1,2) contiene dos datos que dan a entender el comportamiento de la varianza de los datos.
#VALIDACIÓN Model2=Arima(xv_SA,order = c(0,1,2)) Box.test(Model2$residuals, lag = 20, type = “Ljung-Box”)
Model2=Arima(xv_SA,order = c(0,1,2))
Box.test(Model2$residuals, lag = 20, type = "Ljung-Box")##
## Box-Ljung test
##
## data: Model2$residuals
## X-squared = 22.318, df = 20, p-value = 0.3235La prueba Ljung-Box permite entender que los residuos no son relevantes ya que se comportan de manera aleatoria.
#PRONÓSTICO
#Se hace pronpostico a 3 y 06 meses. Se elije el parametro c(95) para que sea a un nivel de confianza del 95%
Pronostico3=forecast(Model2,level= c(95), h=3) plot(Pronostico3) Pronostico3 Pronostico6=forecast(Model2,level= c(95), h=6) plot(Pronostico6) Pronostico6
Pronostico3=forecast(Model2,level= c(95), h=3)
plot(Pronostico3)
Pronostico3## Point Forecast Lo 95 Hi 95
## Jan 2023 145214493 109809895 180619092
## Feb 2023 142066917 104064558 180069276
## Mar 2023 142066917 103815613 180318221Pronostico6=forecast(Model2,level= c(95), h=6)
plot(Pronostico6)
Pronostico6## Point Forecast Lo 95 Hi 95
## Jan 2023 145214493 109809895 180619092
## Feb 2023 142066917 104064558 180069276
## Mar 2023 142066917 103815613 180318221
## Apr 2023 142066917 103568277 180565557
## May 2023 142066917 103322520 180811313
## Jun 2023 142066917 103078313 181055521La variable mantendrá un valor constante a lo largo del pronóstico con un intervalo de confianza del 95%, lo que da a entender cómo se comportará el mercado de las exportaciones para la toma de decisiones.