R Notebook

Series con Tendencia Aleatoria- ARIMA

Consultar un activo financiero desde Yahoo finance

library(quantmod) # api yahoo finance
library(ggplot2) # graficos
library(dplyr) # pandas de r
## serie de tiempo
library(pdfetch)
library(urca)
library(TSA)
library(fBasics)
library(TSstudio)
library(dygraphs)
library(car)
library(forecast)
library(corrplot)
library(fpp2)
library(MASS)
library(tseries)
library(lmtest)

Concexion con Api

datos=as.data.frame(pdfetch_YAHOO(c("GOOG"),from = as.Date("2010-01-01"),to=as.Date("2022-11-18")))
tail(datos)

Dimensión

dim(datos)

[1] 3243    6

names(datos)=c( "open", "high","low","close","adjclose","volume"  )
datos$date=as.Date(row.names(datos))
head(datos)

Separamos en Train y Test

min(datos$date)

[1] "2010-01-04"

max(datos$date)

[1] "2022-11-17"

df_train= datos[datos$date <= "2022-10-31",]
df_test= datos[datos$date>"2022-10-31",]

Análisis Gráfico

ggplot(aes(x=date,y=close),data = df_train)+geom_line(color="darkblue")

Evolución del test

ggplot(aes(x=date,y=close),data = df_test)+geom_line(color="darkblue")

ts_plot(df_test[,c("date","close")],title = "Evolución Precio Google día a día")

¿Es Necesario Estabilizar la Varianza?

• Logaritmo Natural
• Potencia2
• Raiza cuadrada
• Invertir la serie

library(MASS)
box_cox=boxcox(close~ date,
               data=df_train,
               lambda=c(0,0.5,1))

lambda=box_cox$x[which.max(box_cox$y)]
lambda

[1] 0

df_train$LogClose=log(df_train$close)

ggplot(aes(x=date,y=LogClose),data = df_train)+geom_line(color="darkblue")

Análisis de las funciones de autocorrelación simple y parcial

datos.train.ts=ts(df_train$LogClose) # formato de series de tiempo
plot.ts(datos.train.ts,tipe="0", col="darkgreen")

Función de autocorrelación simple

acf(datos.train.ts,lag.max = 550, main="Funcion Autocorrelacion Simple")

pacf(datos.train.ts,lag.max = 50, main="Funcion Autocorrelacion Parcial")

• Que no es un modelo estacionario, esta serie tiene tendencia

• por lo tanto la serie tiene tendencia

¿Tendencia es Aleatoria o Deterministica? : Test de Dickey - Fuller

adf.test(datos.train.ts,alternative="stationary")

    Augmented Dickey-Fuller Test

data:  datos.train.ts
Dickey-Fuller = -3.2499, Lag order = 14, p-value = 0.07928
alternative hypothesis: stationary

• P valor es menor a 0.05; la serie es estacionaria en tendencia (modelo matematico para explicar la tendencia)
• P valor es mayor a 0.05, la serie tiene al menos una raiz unitaria por lo tanto es una seria con tendecia aleatoria

debemos aplicar un modelo ARIMA

¿Cuntas veces deberiamos diferenciar la serie?

plot.ts(diff(datos.train.ts))

Gráficas de Autocorrelación Simple y Parcial de la ST en diferencia

par(mfrow=c(2,1))
acf(diff(datos.train.ts))
pacf(diff(datos.train.ts))

• Modelo es ARIMA con p=0, y q=0 ARIMA(0,1,0)

plot.ts(diff(diff(datos.train.ts)))

• Por ahora la Seria es ARIMA, con d=1, cual es p y q

Selección Automatica

auto.arima(datos.train.ts)

Series: datos.train.ts 
ARIMA(3,1,3) 

Coefficients:
          ar1     ar2     ar3     ma1      ma2      ma3
      -0.7998  0.6785  0.8747  0.7577  -0.6899  -0.8580
s.e.   0.0818  0.1391  0.0769  0.0776   0.1267   0.0717

sigma^2 = 0.0002821:  log likelihood = 8617.13
AIC=-17220.27   AICc=-17220.23   BIC=-17177.71

• el Autorima nos sugiere un modelo ARIMA (1,1,3)

• y(1-b)=b0+b1yt-e1+e2+e3

Autorima con un BIC

auto.arima(datos.train.ts,d=1,max.p=5, max.q=5, ic=c("bic"))

Series: datos.train.ts 
ARIMA(0,1,0) 

sigma^2 = 0.0002845:  log likelihood = 8600.42
AIC=-17198.84   AICc=-17198.84   BIC=-17192.76

Autorima con un AIC

auto.arima(datos.train.ts,d=1,max.p=5, max.q=5, ic=c("aic"))

Series: datos.train.ts 
ARIMA(3,1,3) 

Coefficients:
          ar1     ar2     ar3     ma1      ma2      ma3
      -0.7998  0.6785  0.8747  0.7577  -0.6899  -0.8580
s.e.   0.0818  0.1391  0.0769  0.0776   0.1267   0.0717

sigma^2 = 0.0002821:  log likelihood = 8617.13
AIC=-17220.27   AICc=-17220.23   BIC=-17177.71

• Tenemos dos posibles modelos: ARIMA(0,1,0) y un Modelo ARIMA(1,1,3)

Etapa de Estimación

Modelo1 = Arima(datos.train.ts,order = c(1,1,3),include.drift = TRUE,method = "CSS-ML")
summary(Modelo1)

Series: datos.train.ts 
ARIMA(1,1,3) with drift 

Coefficients:

Warning: NaNs produced

          ar1      ma1     ma2      ma3  drift
      -0.0242  -0.0258  0.0077  -0.0181  6e-04
s.e.      NaN      NaN     NaN   0.0109  3e-04

sigma^2 = 0.0002838:  log likelihood = 8606.85
AIC=-17201.71   AICc=-17201.68   BIC=-17165.23

Training set error measures:
                       ME       RMSE        MAE          MPE      MAPE      MASE          ACF1
Training set 9.304636e-07 0.01683029 0.01142656 -0.001693577 0.3251243 0.9985824 -0.0001070563

• y(1-0.000006)=0.96yt-1–1.012e1+0.06e2 …..

coeftest(Modelo1)

Warning: NaNs produced

z test of coefficients:

         Estimate  Std. Error z value Pr(>|z|)  
ar1   -0.02418651         NaN     NaN      NaN  
ma1   -0.02584458         NaN     NaN      NaN  
ma2    0.00771704         NaN     NaN      NaN  
ma3   -0.01806271  0.01088007 -1.6602  0.09688 .
drift  0.00055851  0.00027932  1.9995  0.04555 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

• Verificamos si el proceso es invertible

autoplot(Modelo1)

• Normalidad de los residuales

residuales= Modelo1$residuals
acf(residuales)

pacf(residuales)

checkresiduals(residuales,lag=15)

Test Normalidad

shapiro.test(residuales)

    Shapiro-Wilk normality test

data:  residuales
W = 0.91669, p-value < 2.2e-16

• p valor menor a 0.05; no hay normalidad

jarque.bera.test(residuales)

    Jarque Bera Test

data:  residuales
X-squared = 9174.4, df = 2, p-value < 2.2e-16

• p valor menor a 0.05, no hay normalidad

• La normalidad puede estar afectada por valores atipicos; hay que retomar

Evaluar el poder predictivo

ajust=Modelo1$fitted # predicciones de train  
ts.plot(datos.train.ts,ajust)
lines(ajust, col="red", type="o", cex=.5)

Evaluacion Test

df_test$LogClose=log(df_test$close)
datos.test.ts=ts(df_test$LogClose) # formato de series de tiempo
ts.plot(datos.test.ts)

Preciccion 24 dias

Tiempo=13
z_pred=forecast(Modelo1,h=Tiempo,level = c(80,95),biasadj = TRUE)
z_pred=as.data.frame(z_pred)
rownames(z_pred) <- 1:nrow(z_pred)
z_pred

ts.plot(ts(z_pred$`Point Forecast`),datos.test.ts,type="o")
lines(ts(z_pred$`Point Forecast`),col="red")

library(MLmetrics)

MAPE(exp(df_test$LogClose),exp(z_pred$`Point Forecast`))

[1] 0.05294389

• 6%

RMSE

RMSE(exp(df_test$LogClose),exp(z_pred$`Point Forecast`))

[1] 5.923233

ts.plot(exp(ts(z_pred$`Point Forecast`)),exp(datos.test.ts),type="o")
lines(exp(ts(z_pred$`Point Forecast`)),col="red")

cbind(exp(df_test$LogClose),exp(z_pred$`Point Forecast`))

       [,1]     [,2]
 [1,] 90.50 94.88985
 [2,] 87.07 94.85301
 [3,] 83.49 94.94275
 [4,] 86.70 94.99490
 [5,] 88.65 95.04799
 [6,] 88.91 95.10109
 [7,] 87.40 95.15422
 [8,] 94.17 95.20738
 [9,] 96.73 95.26057
[10,] 96.03 95.31379
[11,] 98.72 95.36704
[12,] 98.99 95.42032
[13,] 98.50 95.47363

Componentes Estacionales Complejas

• ARIMA para la parte stacional

SARIMA

par(mfrow=c(2,1)) 
Acf(diff(datos.train.ts),208,ci=0)
Pacf(diff(datos.train.ts),208,ci=0)

par(mfrow=c(2,1)) 
Acf(diff(diff(datos.train.ts)),208,ci=0)
Pacf(diff(diff(datos.train.ts)),208,ci=0)

• Estacional ARIMA(0,1,1) para la estacionalidad

Modelo1 = Arima(datos.train.ts,order = c(1,1,2),seasonal = list(order=c(1,1,2),period=12),include.drift = TRUE,method = "CSS-ML")

Warning: No drift term fitted as the order of difference is 2 or more.

coeftest(Modelo1)

z test of coefficients:

      Estimate Std. Error  z value  Pr(>|z|)    
ar1   0.956077   0.034123  28.0182 < 2.2e-16 ***
ma1  -1.006639   0.038397 -26.2164 < 2.2e-16 ***
ma2   0.040317   0.018165   2.2195   0.02645 *  
sar1 -0.732106   0.151346  -4.8373 1.316e-06 ***
sma1 -0.249149   0.146974  -1.6952   0.09004 .  
sma2 -0.750849   0.146786  -5.1153 3.133e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residuales

checkresiduals(Modelo1$residuals,lag=13)

shapiro.test(Modelo1$residuals)

    Shapiro-Wilk normality test

data:  Modelo1$residuals
W = 0.92081, p-value < 2.2e-16

Pronostico

Tiempo=13
z_pred=forecast(Modelo1,h=Tiempo,level = c(80,95))
z_pred=as.data.frame(z_pred)
rownames(z_pred) <- 1:nrow(z_pred)
z_pred

ts.plot(ts(z_pred$`Point Forecast`),datos.test.ts,type="o")
lines(ts(z_pred$`Point Forecast`),col="red")

Datos Atipicos