Modelos ARIMA - Ejemplo: Oro
Pam
2023-08-04
Introducción
El motivo de la introducción de los modelos ARIMA nace del hecho de que no se puede trabajar con una serie temporal no estacionaria. Se dice que una serie es estacionaria cuando su media, varianza y autocovarianza son invariantes en el tiempo. La mayoría de series temporales económicas no son estacionarias pero diferenciándolas un número determinado de veces la serie original se transforma en estacionaria, con lo cual ya se podría aplicar la metodología de los modelos ARIMA.
Concepto ARIMA
En general, se dice que una serie temporal Yt admite una representación autoregresiva integrada y de medias móviles de órdenes p , d y q respectivamente, y se denota por ARIMA( p , d , q ).
\(Y_{t} = c + \phi_{1} y_{dt-1}+ \phi_{2} y_{dt-p} +...+ \phi_{p} y_{dt-p} +\theta_{1} e_{t-1}+\theta_{2} e_{t-2}+...+\theta_{q} e_{t-q}\)
Donde p denota el número de términos autoregresivos, d el número de veces que la serie debe ser diferenciada para hacerla estacionaria y q el número de términos de la media móvil invertible.
La construcción de los modelos se lleva de manera iterativa mediante un proceso en el que se puede distinguir cuatro etapas:
Identificación: Utilizando los datos ordenados cronológicamente se intentara sugerir un modelo que merezca la pena ser investigada. El objetivo es determinar los valores que sean apropiados para reproducir la serie de tiempo. Comprende los pasos 1, 2 y 3 del ejemplo próximo.
Análisis y diferenciación de la serie temporal. Consiste en examinar la estacionariedad, los diagramas de autocorrelación, también conocidos como ACF y PACF, y elección del orden del modelo. Comprende los pasos 4 y 5.
Ajuste de un modelo ARIMA. Obtención de los coeficientes de determinación. Paso 6.
Predicción. Una vez seleccionado el mejor modelo candidato se pueden hacer pronósticos en términos probabilísticos de los valores futuros.
El objetivo del siguiente análisis es: Trazar, examinar y preparar series para modelar. Extrae el componente de estacionalidad de la serie temporal. Probar la estacionalidad y aplicar las transformaciones apropiadas. Elija el orden de un modelo ARIMA. Pronostica la serie.
La base de datos con la que hemos realizado el análisis está extraída de Yahoo Finance (a través de la técnica scraping) y trata sobre el tipo de cambio del peso mexicano con el dolar americano.
Para realizar el ejemplo, comenzamos instalando y ejecutando las librerías necesarias:
Instalación de librerias
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)Ejecutando las librerías
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 zoolibrary(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)PRIMERA ETAPA: Identificación
Lectura y preparación de los datos.
Quantmod es un paquete de R que proporciona un conjunto de herramientas para el modelado y análisis financiero cuantitativo. Permite a los usuarios acceder y manipular datos financieros de varias fuentes, incluido Yahoo Finance.
# Este comando importa todos los datos disponibles en yahoo finance
#ORO <- getSymbols('GC=F', src='yahoo', auto.assign=FALSE)
# Este comando importa los datos de acuerdo con un periodo de tiempo determinado
ORO <- getSymbols('GC=F', src='yahoo', from = as.Date("2023-01-01"),to=as.Date("2023-08-04"), auto.assign = FALSE)Graficando la serie
chartSeries(ORO, name="GC=F", subset="last 6 months", theme=chartTheme("white"))
Datos de yahoo finance:
data_ORO <- data.frame(ORO, tiempo = as.Date(rownames(data.frame(ORO))))
head(data_ORO)## GC.F.Open GC.F.High GC.F.Low GC.F.Close GC.F.Volume GC.F.Adjusted
## 2023-01-03 1836.2 1839.7 1836.2 1839.7 29 1839.7
## 2023-01-04 1845.6 1859.1 1845.6 1852.8 25 1852.8
## 2023-01-05 1855.2 1855.2 1834.8 1834.8 24 1834.8
## 2023-01-06 1838.4 1868.2 1835.3 1864.2 26 1864.2
## 2023-01-09 1867.0 1880.0 1867.0 1872.7 62 1872.7
## 2023-01-10 1877.8 1878.1 1871.6 1871.6 101 1871.6
## tiempo
## 2023-01-03 2023-01-03
## 2023-01-04 2023-01-04
## 2023-01-05 2023-01-05
## 2023-01-06 2023-01-06
## 2023-01-09 2023-01-09
## 2023-01-10 2023-01-10attach(data_ORO)Separando la serie close:
base1 = data.frame(tiempo, GC.F.Close)
names (base1) = c("tiempo","ORO")
base1 <- na.omit(base1) #eliminando datos ominitidos "NA"
#base1
head(base1, n = 10)## tiempo ORO
## 1 2023-01-03 1839.7
## 2 2023-01-04 1852.8
## 3 2023-01-05 1834.8
## 4 2023-01-06 1864.2
## 5 2023-01-09 1872.7
## 6 2023-01-10 1871.6
## 7 2023-01-11 1874.6
## 8 2023-01-12 1895.5
## 9 2023-01-13 1918.4
## 10 2023-01-17 1907.2Graficando la serie
ggplot(base1, aes(x = tiempo, y = ORO)) + geom_line() + labs(title = "Tiempo", x = "Fecha", y = "peso / quilates")
Con la función geom_smooth en ggplot2, puedes suavizar una serie temporal para analizar tendencias o patrones. Por defecto, geom_smooth también mostrará un intervalo de confianza alrededor de la suavización. Si deseas eliminar el intervalo de confianza, puedes utilizar el argumento *se = FALSE. Aquí tienes un ejemplo de cómo usarlo:
ggplot(base1, aes(x = tiempo, y = ORO)) + geom_line() + geom_smooth(se = FALSE)+ labs(title = "Tipo de cambio", x = "Fecha", y = "Pesos / quilates")## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Primero, calculamos el componente estacional del uso de datos stl(). A continuación,hallamos el componente estacional de la serie mediante suavizado y ajusta la serie original restando la estacionalidad.
ORO_ma = ts(na.omit(base1$ORO), frequency=30)
decomp = stl(ORO_ma, s.window="periodic")
deseasonal_base1 <- seasadj(decomp)
plot(decomp)
SEGUNDA ETAPA: Análisis y diferenciación de la serie temporal
Análisis de estacionariedad
La instalación de un modelo ARIMA requiere que la serie sea estacionaria . Se dice que una serie es estacionaria cuando su media, varianza y autocovarianza son invariantes en el tiempo.
Esta suposición tiene un sentido intuitivo: dado que ARIMA usa retardos previos de series para modelar su comportamiento
adf.test(ORO_ma, alternative = "stationary")##
## Augmented Dickey-Fuller Test
##
## data: ORO_ma
## Dickey-Fuller = -2.0366, Lag order = 5, p-value = 0.5612
## alternative hypothesis: stationaryContraste de hipótesis:
- H0: No estacionaria (tiene raíz unitaria)
- H1: Estacionaria (No tiene raíz unitaria)
Tras realizar la prueba aumentada de Phillips-prron (ADF), obtenemos un p-valor = 0.5612. Como el p-valor > 0.05, no rechazamos H0. Podemos concluir que nuestra serie temporal no es estacionaria.
pp.test(ORO_ma, alternative = "stationary")##
## Phillips-Perron Unit Root Test
##
## data: ORO_ma
## Dickey-Fuller Z(alpha) = -6.7634, Truncation lag parameter = 4, p-value
## = 0.7263
## alternative hypothesis: stationaryContraste de hipótesis:
- H0: No estacionaria (tiene raíz unitaria)
- H1: Estacionaria (No tiene raíz unitaria)
Tras realizar la prueba aumentada de Dickey-Fuller (ADF), obtenemos un p-valor =0.7263. Como el p-valor > 0.05, no rechazamos H0. Podemos concluir que nuestra serie temporal no es estacionaria.
Continuamos con las autocorrelaciones.
Las ACF proporcionan información sobre cómo una observación influye en las siguientes.
Acf(ORO_ma, main='')
Las PACF proporcionan la relación directa existente entre observaciones separadas por k retardos.
Pacf(ORO_ma, main='')
Para realizar un modelo ARIMA, la serie temporal debe ser estacionaria. Para conseguir esta estacionariedad, la diferenciaremos.
ORO_d1 = diff(deseasonal_base1, differences = 1)
plot(ORO_d1)
Para comprobar que la serie es, efectivamente, estacionaria, hacemos de nuevo el test aumentado de Dickey-Fuller.
adf.test(ORO_d1, alternative = "stationary")## Warning in adf.test(ORO_d1, alternative = "stationary"): p-value smaller than
## printed p-value##
## Augmented Dickey-Fuller Test
##
## data: ORO_d1
## Dickey-Fuller = -4.0911, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryEl p-valor obtenido es de 0.01 < 0.05. Rechazamos H0, 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(ORO_d1, main='ACF para la serie diferenciada 1 vez')
Pacf(ORO_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 = 255.3: log likelihood = -611.75
## AIC=1225.5 AICc=1225.53 BIC=1228.49Usando la notación ARIMA presentada anteriormente, el modelo ajustado se puede escribir como:
\(Y_{d_{t}} = -0.0082 -0.014 Y_{t-1}+E\)
donde E es un error y la serie original se diferencia con la orden 1.
CUARTA ETAPA: Predicción
Si los parámetros de orden del modelo y la estructura se especifican correctamente, no esperaríamos que existieran autocorrelaciones significativas.
Hemos ajustado un modelo que puede producir un pronóstico, pero comenzamos examinando los gráficos ACF y PACF para los residuos del modelo. Si los parámetros de orden del modelo y la estructura se especifican correctamente, no esperaríamos que existieran autocorrelaciones significativas.
modeloarima<-auto.arima(deseasonal_base1, seasonal=FALSE)
modeloarima## Series: deseasonal_base1
## ARIMA(0,1,0)
##
## sigma^2 = 255.3: log likelihood = -611.75
## AIC=1225.5 AICc=1225.53 BIC=1228.49tsdisplay(residuals(modeloarima), lag.max=10, main='(1,1,0) Model Residuals')
La autocorrelación está dentro de las bandas, por lo que 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(5, 28)
## End = c(6, 27)
## Frequency = 30
## [1] 1924.84 1924.84 1924.84 1924.84 1924.84 1924.84 1924.84 1924.84 1924.84
## [10] 1924.84 1924.84 1924.84 1924.84 1924.84 1924.84 1924.84 1924.84 1924.84
## [19] 1924.84 1924.84 1924.84 1924.84 1924.84 1924.84 1924.84 1924.84 1924.84
## [28] 1924.84 1924.84 1924.84