BASE DE DATOS

library(rmarkdown)
library(foreign)
base = read.csv(file.choose())
class(base)
## [1] "data.frame"

Base de datos del precio del petróleo crudo, fijado por la Organización de países exportadores de petróleo (OPEP).

Periodicidad

Se tienen registros a partir del mes de Agosto de 1993, en una frecuencia mensual.

Unidad de medida

Dólares americanos

Fuente

Organización de países exportadores de petróleo (OPEP) https://www.indexmundi.com/es/precios-de-mercado/?mercancia=petroleo-crudo&meses=300

DECLARAR COMO SERIE DE TIEMPO

PPDP <- ts(base$Precios, frequency = 12, start = c(1993,8))
class(PPDP)
## [1] "ts"

GRÁFICA DE LOS PRECIOS DEL PETRÓLEO

library(forecast)
plot(PPDP, col = "red", main = "Precios del Barril Petróleo", ylab = "Dólares", xlab ="Tiempo", lwd=2, type="l", pch=10)

DESCOMPOSICIÓN DE LA SERIE

fit<-decompose(PPDP)
library(ggplot2)
plot (fit, col = "blue", xlab ="Tiempo", lwd=1, type="l", pch=10)

En la pesente gráfica se puede observar que la serie de los precios del barril de petróleo presenta todos los componentes de una serie de tiempo.

plot(fit$trend, col= "pink", main = "Tendencia de los Precios del Barril Petróleo", ylab = "Dólares", xlab ="Tiempo", lwd=1, type="o", pch=2)

Analizando especificamente la tendencia de la serie, se observa una tendencia positiva hasta el año 2008 donde la tendencia cambia y tiene un comportamiento negativo durante el año 2009 y una recuperación en el 2010 y se estabiliza hasta el 2014 donde se presenta una pendiente negativa muy pronunciada hasta el año 2016 donde se presenta el ultimo punto de inflexión de la tendencia y se presenta de nuevo una pendiente negativa.

plot(fit$random, col= "turquoise", xlab ="Tiempo", lwd=2, type="l", pch=10)

Se puede observar que se tiene un componente irregular aparentemente estacionario ya que a simple vista se observa una media cero, pero tambien se puede observar que apartir del año 2009 se presenta una varinza creciente, lo que sugiere que el componente irregular de la serie no es estacionario.

plot(fit$seasonal, col= "orange",main = "Estacionalidad de los Precios del Barril Petróleo", xlab ="Tiempo", lwd=2, type="o", pch=5 )

En la gráfica de estacionalidad se observa que la serie presenta estacionlidad, ya que este gráfico no es muy confiable, se procede a realizar otros tipos de gráficas

library(forecast)
ggseasonplot(PPDP, year.labels=TRUE, year.labels.left=TRUE)

En este gráfico se puede observar mejor la estacionalidad

#Estacionalidad
ggseasonplot(PPDP, polar="TRUE")

En estas últimas gráficas nos muestran mejor la estacionalidad y se puede observar claramente que la serie no presenta estacionalidad

Por lo tanto se muestra que no hay factortes temporales dentro de un año que influyen en los precios del petróleo.

ESTIMACIÓN DEL MODELO DE REGRESIÓN

library(TSA)
## 
## Attaching package: 'TSA'
## The following objects are masked from 'package:stats':
## 
##     acf, arima
## The following object is masked from 'package:utils':
## 
##     tar
mes.<-season(PPDP)
MOD1<-lm(PPDP~mes.)
summary(MOD1)
## 
## Call:
## lm(formula = PPDP ~ mes.)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -40.824 -28.814  -6.085  21.624  79.306 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    48.2132     6.5882   7.318 2.55e-12 ***
## mes.February    1.0504     9.3171   0.113    0.910    
## mes.March       2.7736     9.3171   0.298    0.766    
## mes.April       4.3984     9.3171   0.472    0.637    
## mes.May         5.1600     9.3171   0.554    0.580    
## mes.June        4.7035     9.4137   0.500    0.618    
## mes.July        5.3106     9.4137   0.564    0.573    
## mes.August      3.6192     9.3171   0.388    0.698    
## mes.September   3.0296     9.3171   0.325    0.745    
## mes.October     2.0276     9.3171   0.218    0.828    
## mes.November    0.9236     9.3171   0.099    0.921    
## mes.December   -0.2304     9.3171  -0.025    0.980    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 32.94 on 286 degrees of freedom
## Multiple R-squared:  0.003412,   Adjusted R-squared:  -0.03492 
## F-statistic: 0.089 on 11 and 286 DF,  p-value: 1

La regresión muestra que los meses no muestran una representación estadistica significativa lo que corrobora lo que se habia visto en las gráficas de estacionalidad. La R cuadrada nos muestra que el modelo no reprenta la relidad, por lo que se concluye que los meses no tienen importancia en el movimiento de los precios del petróleo.

Box Cox

BoxCox.ar(PPDP)
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

Se puede ver que los datos de la serie parecen normales con sesgos a la derecha, lo que puede significar que tiene tendencia que ya se habia observado en las gráficas anteriores, por lo que se procede a realizar otras pruebas.

Como la concentración es al cero se acepta la Hipótesis nula por lo cual se concluye que no hay estacionariedad.

Usar prueba Dickey-Fuller para evaluar el orden de integracion de la serie.

library(urca)
summary(ur.df(PPDP))
## 
## ############################################### 
## # 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 
## -20.4247  -1.7172   0.5162   2.5842  13.5305 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -0.002638   0.004055   -0.65    0.516    
## z.diff.lag  0.421359   0.053135    7.93 4.61e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.202 on 294 degrees of freedom
## Multiple R-squared:  0.1762, Adjusted R-squared:  0.1706 
## F-statistic: 31.45 on 2 and 294 DF,  p-value: 4.192e-13
## 
## 
## Value of test-statistic is: -0.6505 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

Tras realizar la prueba Dickey-Fulller a la serie oridinal se puede observar que la serie no es estacionaria lo que significa que tiene tendencia, por lo tanto se prosede a transformar la serie con logaritmos para reducir la varinza, aplicar diferencia para estabilizar la media.

D1PPDP<-diff(log(PPDP))
autoplot(diff(log(PPDP)))

Con la aplicación de los logaritmos y de la difencia se puede observar que la serie a simple vista ya luce estacionaria, lo cual se comprueba a través de la prueba Dickey-Fuller.

library(urca)
summary(ur.df(D1PPDP))
## 
## ############################################### 
## # 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 
## -0.27510 -0.04863  0.01616  0.05359  0.22353 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -0.71184    0.07103 -10.022   <2e-16 ***
## z.diff.lag -0.03348    0.05844  -0.573    0.567    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08033 on 293 degrees of freedom
## Multiple R-squared:  0.3687, Adjusted R-squared:  0.3643 
## F-statistic: 85.55 on 2 and 293 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -10.022 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

Tras hacer la prueba a la seriere diferenciada se puede ver que la serie ya es estacionaria.

Usar herramientas ACF, PACF, EACF y criterios de Akaike/ Bayes para construir propuestas de modelos.

Analizar los residuos.

ggAcf(D1PPDP)

Tas hacer la autocorrelación parcial se puede obeservar que se sale el resago uno significativamente. Con lo que con base en esta función se realiza la propuesta de un MA(1).

ggPacf(D1PPDP)

Tras realizar la función de autocorrelación parcial se observa que tres residuos salen de los limites. Por lo que se propone un AR(1) y un AR(2).

eacf(D1PPDP)
## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x o o o o o o o o o x  o  o  o 
## 1 x o o o o o o o o o o  o  o  o 
## 2 x o o o o o o o o o o  o  o  o 
## 3 x x o o o o o o o o o  o  o  o 
## 4 x x o o o o o o o o o  o  x  o 
## 5 x x x x o o o o o o o  o  x  o 
## 6 x x x x o o o o o o o  o  x  o 
## 7 x x o x o x o o o o o  o  x  o

No se puede realizar un ARMA

ggtsdisplay(D1PPDP)

Tras ver las gráficas de autocorrelación se mantienen los modelos ya propuestos anteriormente.

bayes <- armasubsets(y=D1PPDP,nar=11,nma=10,y.name='test',ar.method='ols')
## Warning in leaps.setup(x, y, wt = wt, nbest = nbest, nvmax = nvmax,
## force.in = force.in, : 10 linear dependencies found
plot(bayes)

Con la pueba de BAYES se propone un AR(1) que ya se habia propuesto en la función de autocorrelación parcial.

Propuestas.

MA(1)

propuesta1<-arima(PPDP,order = c(0,1,1)) 
propuesta1
## 
## Call:
## arima(x = PPDP, order = c(0, 1, 1))
## 
## Coefficients:
##          ma1
##       0.3297
## s.e.  0.0459
## 
## sigma^2 estimated as 18.32:  log likelihood = -853.3,  aic = 1708.6
checkresiduals(propuesta1)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)
## Q* = 46.694, df = 23, p-value = 0.00245
## 
## Model df: 1.   Total lags used: 24

AR(1)

propuesta2<-arima(PPDP,order = c(1,1,0)) 
propuesta2
## 
## Call:
## arima(x = PPDP, order = c(1, 1, 0))
## 
## Coefficients:
##          ar1
##       0.4177
## s.e.  0.0526
## 
## sigma^2 estimated as 17.51:  log likelihood = -846.62,  aic = 1695.25
checkresiduals(propuesta2)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,0)
## Q* = 33.729, df = 23, p-value = 0.06918
## 
## Model df: 1.   Total lags used: 24

AR(2)

propuesta3<-arima(PPDP,order = c(2,1,0)) 
propuesta3
## 
## Call:
## arima(x = PPDP, order = c(2, 1, 0))
## 
## Coefficients:
##          ar1     ar2
##       0.3971  0.0494
## s.e.  0.0579  0.0579
## 
## sigma^2 estimated as 17.46:  log likelihood = -846.26,  aic = 1696.52
checkresiduals(propuesta3)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,0)
## Q* = 33.585, df = 22, p-value = 0.05411
## 
## Model df: 2.   Total lags used: 24

ARMA(1, 1)

propuesta4<-arima(PPDP,order = c(1,1,1)) 
propuesta4
## 
## Call:
## arima(x = PPDP, order = c(1, 1, 1))
## 
## Coefficients:
##          ar1      ma1
##       0.4807  -0.0756
## s.e.  0.1006   0.1089
## 
## sigma^2 estimated as 17.48:  log likelihood = -846.39,  aic = 1696.78
checkresiduals(propuesta4)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,1)
## Q* = 33.582, df = 22, p-value = 0.05414
## 
## Model df: 2.   Total lags used: 24

Tras el análisis de los residuos y de los pronósticos asi como de las pruebas Ljung-Box se llega a la conclusión que el mejor modelo es el AR(1) no estacional con una diferencia no estacional ya que es la que mejor se comporta y la que mejor criterio de Akaike tiene.

Usar la función auto.arima() y compare sus resultados con el inciso anterior.

propuesta.auto<-auto.arima(PPDP, stepwise=FALSE, approximation=FALSE)
propuesta.auto
## Series: PPDP 
## ARIMA(0,1,2)(2,0,0)[12] 
## 
## Coefficients:
##          ma1     ma2    sar1     sar2
##       0.4115  0.2528  0.0761  -0.1484
## s.e.  0.0562  0.0552  0.0573   0.0571
## 
## sigma^2 estimated as 17.04:  log likelihood=-840.89
## AIC=1691.77   AICc=1691.98   BIC=1710.24
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(0,1,2)(2,0,0)[12]
## Q* = 22.14, df = 20, p-value = 0.3329
## 
## Model df: 4.   Total lags used: 24

Realizando el auto.arima nos arroja una propuesta que no contemplamos la cual se comporta mejor en los residuos y el criterio de Akaike es un poco menor por lo que se tomara la propuesta del auto auto.arima

Propuesta5

propuesta5<-Arima(PPDP,order = c(0,1,2),seasonal = c(2,0,0))
propuesta5
## Series: PPDP 
## ARIMA(0,1,2)(2,0,0)[12] 
## 
## Coefficients:
##          ma1     ma2    sar1     sar2
##       0.4115  0.2528  0.0761  -0.1484
## s.e.  0.0562  0.0552  0.0573   0.0571
## 
## sigma^2 estimated as 17.04:  log likelihood=-840.89
## AIC=1691.77   AICc=1691.98   BIC=1710.24

Presente la ecuacion final y describa.

\[ Precio=\theta e_{t-1}+ \theta e_{t-2}+ \Phi Y_{t-1}+ \Phi Y_{t-2}+2e_t \] \[ Precio= 0.4115e_{t-1} + 0.2528e_{t-2} + 0.0761Y_{t-1} -0.1484Y_{t-2}+2e_t \]

Crear un pronóstico a dos años.

autoplot(forecast(propuesta5))

El pronóstico nos muestra que la serie se estabiliza en su tendencia, lo que significa que los precios del petróleo crudo se mantendran en los próximos dos años.