Modelos SARIMA(Arima Estacionales).¿Qué son y cómo usarlos para Predecir?

CODIGO R DEL VIDEO

## PROCESOS ARIMA Y SARIMA

#Nuestro objetivo será estudiar la evolución de una serie de tiempo
#buscando ser capaz de modelizar y obtener predicciones

#ARIMA <- modelo estadístico que utiliza variaciones y
#regresiones de datos estadísticos con el fin de encontrar
#patrones para una predicción hacia el futuro.
#Se trata de un modelo dinámico de series temporales, es decir,
#las estimaciones futuras vienen explicadas por los datos del pasado
#y no por variables independientes.

# ARIMA(p,d,q) Compuesto por 3 parámetros (Parte o Retardos Regulares)

# SARIMA(p,d.q),(P,D,Q), periodo= «x» (Parte estacional)

## Los modelos SARIMA tienen los momentos y correlogramas iguales
# que los ARIMA pero en vez de considerar los retardos regulares
## nos fijaremos en los retardos estacionales

### Vamos a obtener el mejor Arima Estacional

### Diremos que una serie es estacional cuando su media no es constante en
### en el tiempo pero varía de forma periódica, ciclicamente E(Xt) = E(Xt+s)

# En serie diarias con estacionalidad Semanal <- S = 7
# En series Mensuales con estacionalidad Anual <- S = 12
# En series Horarias con estacionalidad diaria <- S= 24
# En series Bimensuales la estacionalidad anual es <- S= 6 y
# analogamente con las Cuatrimestrales <- S = 3 y las trimestrales S=4

### Dos formas de obtener la Estacionalidad en Arima

## 1. Periodrama + funciones de autocorrelacion para el modelo

## #Estacionalidad
#La realizamos llevando a cabo un análisis espectral, toda serie temporal aleatoria se puede representar mediante una combinación de ondas seno(x) o coseno (x) siendo x la serie temporal. Dichas ondas difieren en el tiempo que tardan en completarse.
#No nos vamos a extender en el cálculo debido a que se usa un cálculo complejo (Transformada de Fourier).

## 2. Fuerza bruta + AIC

#install.packages(«quantmod»)
library(quantmod)
#install.packages(«tseries»)
library(tseries)
#install.packages(«forecast»)
library(forecast)
#install.packages(«TSA»)
library(TSA)
#install.packages(«Quandl»)
library(Quandl)

getSymbols(«CPIAUCNS»,src=»FRED») #IPC consumidores urbanos

#getSymbols(«UNRATE»,src=»FRED») #tasa de desmepleo EEUU
#getSymbols(«GS10″,src=»FRED») #letras del tesoro a 10 añs
#getSymbols(«FEDFUNDS»,src=»FRED») #tasa interes fondos federales
#getSymbols(«EXUSEU»,src=»FRED») #EURSUD mensual
#getSymbols(«CCRETT01USM661N»,src=»FRED») #Tipo de cambio efectivo IPC
#getSymbols(«M2NS»,src=»FRED») #Stock de dinero M2
#getSymbols(«PAYEMS»,src=»FRED») #trabajadores no rurales
#getSymbols(«BOPGSTB»,src=»FRED») #Balance cuenta, Diferencia entre exp e imp

IPC=CPIAUCNS[«2001-01::2019-12»]

#desempleo=UNRATE[«2001-01::2019-12»]
#bono10=GS10[«2001-01::2019-12»]
#Fedfunds=FEDFUNDS[«2001-01::2019-12»]
#EUR=EXUSEU[«2001-01::2019-12»]
#credit=CCRETT01USM661N[«2001-01::2019-12»]
#pay=PAYEMS[«2001-01::2019-12»]
#money=M2NS[«2001-01::2019-12»]
#trade=BOPGSTB[«2001-01::2019-12»]

#IPC

y= IPC$CPIAUCNS

plot(y,type=»l»)

adf.test(y,alternative=»stationary»,k=0)

#vemos que el p-value es mayor a 0.05

#Por el grafico vemos como hay una tendencia. Para estar seguros
#probamos el test de autocorrelacion

acf(y)

# -2.003

dy=na.omit(diff(y)) # Diferenciamos para sacar una serie estacionaria
adf.test(dy,alternative=»stationary»,k=0)

plot(dy,type=»l»)

## Usar periodograma

p=periodogram(dy)
max(p$spec)
p$freq[match(max(p$spec),p$spec)]
1/p$freq[match(max(p$spec),p$spec)]
# 12

## Eleccion del modelo con autocorrelacion

#ahora calcularemos los test de correlacion y correlacion parcial para ver que
#que modelos utilizaremos para intentar predecir la serie

acf(dy)
pacf(dy)

##FUNCION AUTOARIMA + PRUEBA DE MODELOS

# NO he podido hacer autoarima estacional

Modelo1=auto.arima(y)
Modelo1 #arima 2,1,4

Ms010= Arima(y,order=c(2,1,4),seasonal=list(order=c(0,1,0),period=12))
Ms010
#537.85

ms111= Arima(y,order=c(2,1,4),seasonal=list(order=c(1,1,1),period=12))
ms111
#409.01

ms212= Arima(y,order=c(2,1,4),seasonal=list(order=c(2,1,2),period=12))
ms212
#410.07

# Modelo final

#Parte1 con autoarima
#parte 2 probando con varios modelos y eligiendo el menor AIC
#parte 3 — hecho al principio con la estacionalidad

MSAR= Arima(y,order=c(2,1,4),seasonal=list(order=c(1,1,1),period=12))
MSAR
#Arima()

pred2 <- forecast(MSAR, h=10)
plot(pred2,xlab=»IPC»,col=»red»)
Prediccion2<- pred2$mean
Prediccion2
PredYERROR_2=data.frame(predict(MSAR,n.ahead=3))
PredYERROR_2

## ALGORITMO FUERZA BRUTA

FuerzaUniversal = function( datos, valores ){

#Parallel
require(doParallel)
no_cores <- (detectCores() – 1)
cl<- makeCluster(no_cores)
registerDoParallel(cl)

y= datos


primeraVuelta = FALSE

primeraVuelta = FALSE
for(order_x in 0:2){
for(order_y in 0:2){
for(order_z in 0:2){
for(periodo in 12:12){

for(seasonal_x in 0:2){
for(seasonal_y in 0:2){
for(seasonal_z in 0:2){
actual = tryCatch(arima(datos,order=c(order_x,order_y,order_z),seasonal=list(order=c(seasonal_x,seasonal_y,seasonal_z),period=periodo)),
error =function(err)FALSE)

if( !is.logical(actual)){
print(c(order_x,order_y,order_z,seasonal_x,seasonal_y,seasonal_z,periodo, actual$aic))
if(primeraVuelta == FALSE){
primeraVuelta = TRUE
menor = actual$aic

aOK = order_x
bOK = order_y
cOK = order_z

xOk = seasonal_x
yOk = seasonal_y
zOk = seasonal_z
period = periodo
}
if (actual$aic < menor){
menor = actual$aic
aOK = order_x
bOK = order_y
cOK = order_z


xOk = seasonal_x
yOk = seasonal_y
zOk = seasonal_z
period = periodo
}
}
}
}
}
}
}
}
}
stopCluster(cl)
coordenadas = c(aOK,bOK,cOK,xOk,yOk,zOk,periodo, menor)
}

## POnemos la serie a la que se lo queremos aplicar
y= y

goku = FuerzaUniversal(y)

goku

#ModeloFinalSARIMA <- (0,1,1)(1,1,2) (12)