Ejercicio ARIMA con Tesla (TSLA)

install.packages('quantmod') # librería para hacer scraping
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.3'
## (as 'lib' is unspecified)
install.packages('tseries')
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.3'
## (as 'lib' is unspecified)
install.packages('timeSeries')
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.3'
## (as 'lib' is unspecified)
install.packages('forecast')
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.3'
## (as 'lib' is unspecified)
install.packages('xts')
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.3'
## (as 'lib' is unspecified)
install.packages('ggplot2')
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.3'
## (as 'lib' is unspecified)
library(quantmod)
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: TTR
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(tseries)
library(timeSeries)
## Loading required package: timeDate
## 
## Attaching package: 'timeSeries'
## The following object is masked from 'package:zoo':
## 
##     time<-
library(forecast)
library(xts)
library(ggplot2)

1. Identificación de datos

TSLA <- getSymbols('TSLA', src='yahoo', from = as.Date("2021-11-01"),to=as.Date("2023-07-27"), auto.assign = FALSE)

Graficando la serie

chartSeries(TSLA, name="TSLA", subset="last 6 months", theme=chartTheme("white"))

Datos de yahoo finance:

data_TSLA <- data.frame(TSLA, tiempo = as.Date(rownames(data.frame(TSLA))))
head(data_TSLA)
##            TSLA.Open TSLA.High TSLA.Low TSLA.Close TSLA.Volume TSLA.Adjusted
## 2021-11-01  381.6667  403.2500 372.8867   402.8633   168146100      402.8633
## 2021-11-02  386.4533  402.8633 382.0000   390.6667   128213400      390.6667
## 2021-11-03  392.4433  405.1300 384.2067   404.6200   103885500      404.6200
## 2021-11-04  411.4700  414.4967 405.6667   409.9700    76192200      409.9700
## 2021-11-05  409.3333  413.2900 402.6667   407.3633    64886400      407.3633
## 2021-11-08  383.2633  399.0000 377.6667   387.6467   100337100      387.6467
##                tiempo
## 2021-11-01 2021-11-01
## 2021-11-02 2021-11-02
## 2021-11-03 2021-11-03
## 2021-11-04 2021-11-04
## 2021-11-05 2021-11-05
## 2021-11-08 2021-11-08
attach(data_TSLA)

Separando la serie close:

base1 = data.frame(tiempo, TSLA.Close)
names (base1) = c("tiempo","TSLA")
base1 <- na.omit(base1) #eliminando datos ominitidos "NA"
#base1
head(base1, n = 10)
##        tiempo     TSLA
## 1  2021-11-01 402.8633
## 2  2021-11-02 390.6667
## 3  2021-11-03 404.6200
## 4  2021-11-04 409.9700
## 5  2021-11-05 407.3633
## 6  2021-11-08 387.6467
## 7  2021-11-09 341.1667
## 8  2021-11-10 355.9833
## 9  2021-11-11 354.5033
## 10 2021-11-12 344.4733

Graficando la serie

ggplot(base1, aes(x = tiempo, y = TSLA)) + geom_line() +  labs(title = "Cotizaciones de Tesla", x = "Fecha", y = "Precios por accion")

Suavizando la serie temporal

ggplot(base1, aes(x = tiempo, y = TSLA)) + geom_line() + geom_smooth(se = FALSE)+ labs(title = "Cotizaciones de Tesla", x = "Fecha", y = "Precios por accion")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

SEGUNDA ETAPA: Análisis y diferenciación de la serie temporal

Calculamos del componente estacional

TSLA_ma = ts(na.omit(base1$TSLA), frequency=30)
decomp = stl(TSLA_ma, s.window="periodic")
deseasonal_base1 <- seasadj(decomp)
plot(decomp)

adf.test(TSLA_ma, alternative = "stationary")
## 
##  Augmented Dickey-Fuller Test
## 
## data:  TSLA_ma
## Dickey-Fuller = -2.122, Lag order = 7, p-value = 0.5261
## alternative hypothesis: stationary

Prueba aumentada de Dickey-Fuller

p-value = 0.5261 > 0.05 NSRHo no existe estacionareidad

pp.test(TSLA_ma, alternative = "stationary")
## 
##  Phillips-Perron Unit Root Test
## 
## data:  TSLA_ma
## Dickey-Fuller Z(alpha) = -8.8064, Truncation lag parameter = 5, p-value
## = 0.6179
## alternative hypothesis: stationary

Prueba aumentada de Phillips-Perron

p-value = 0.6179 > 0.05 NSRHo no existe estacionareidad

Continuamos con las autocorrelaciones.

Las ACF proporcionan información sobre cómo una observación influye en las siguientes.

Acf(TSLA_ma, main='')

Las PACF proporcionan la relación directa existente entre observaciones separadas por k retardos.

Pacf(TSLA_ma, main='')

Para el modelo ARIMA, la serie temporal debe ser estacionaria. Para conseguir esta estacionariedad, la diferenciaremos.

TSLA_d1 = diff(deseasonal_base1, differences = 1)
plot(TSLA_d1)

Para comprobar que la serie es, efectivamente, estacionaria, hacemos de nuevo el test aumentado de Dickey-Fuller.

adf.test(TSLA_d1, alternative = "stationary")
## Warning in adf.test(TSLA_d1, alternative = "stationary"): p-value smaller than
## printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  TSLA_d1
## Dickey-Fuller = -7.2859, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary

El p-value = 0.01 < 0.05. RH0, la serie temporal es estacionaria.

Al trazar la serie diferenciada, vemos un patron oscilante al rededor de 0, sin una tendencia fuerte visible. Esto sugiere que la diferenciacion de los terminos del orden 1 es suficiente y debe incluirse en el modelo.

A continuación, los picos en rezagos particulares de la serie diferenciada pueden ayudar a informar la elección de p o q para nuestro modelo:

Acf(TSLA_d1, main='ACF para la serie diferenciada 1 vez')

Pacf(TSLA_d1, main='ACF para la serie diferenciada 1 vez')

TERCERA ETAPA: Ajuste de un modelo ARIMA.

Ajustamos el modelo.

auto.arima(deseasonal_base1, seasonal=FALSE)
## Series: deseasonal_base1 
## ARIMA(0,1,0) 
## 
## sigma^2 = 101.1:  log likelihood = -1617.55
## AIC=3237.09   AICc=3237.1   BIC=3241.17

CUARTA ETAPA: Predicción

modeloarima<-auto.arima(deseasonal_base1, seasonal=FALSE)
modeloarima
## Series: deseasonal_base1 
## ARIMA(0,1,0) 
## 
## sigma^2 = 101.1:  log likelihood = -1617.55
## AIC=3237.09   AICc=3237.1   BIC=3241.17
tsdisplay(residuals(modeloarima), lag.max=10, main='(0,1,0) Model Residuals')

La autocorrelación no está dentro de las bandas, por lo que no está “bien comportada”.

QUINTA ETAPA

Podemos especificar el horizonte de pronóstico h periodos por delante para que se realicen las predicciones, y usar el modelo ajustado para generar dichas predicciones:

prediccion <- forecast(modeloarima, h=30)
plot(prediccion)

Los datos proyectados para los siguientes 30 días son:

tail(prediccion$mean,30)
## Time Series:
## Start = c(15, 16) 
## End = c(16, 15) 
## Frequency = 30 
##  [1] 255.7919 255.7919 255.7919 255.7919 255.7919 255.7919 255.7919 255.7919
##  [9] 255.7919 255.7919 255.7919 255.7919 255.7919 255.7919 255.7919 255.7919
## [17] 255.7919 255.7919 255.7919 255.7919 255.7919 255.7919 255.7919 255.7919
## [25] 255.7919 255.7919 255.7919 255.7919 255.7919 255.7919