INTRODUCCIÓN

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 \(Y_{t}\) 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_{d_{t-1}} + \phi_p Y_{d_{tp}} + \cdots + \theta_1 \varepsilon_1 + \theta_q \varepsilon_{tq} + \varepsilon_t \]

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.
  • 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.
  • Ajuste de un modelo ARIMA. Obtención de los coeficientes de determinación.
  • Predicción. Una vez seleccionado el mejor modelo candidato se pueden hacer pronósticos en términos probabilísticos de los valores futuros.

La base de datos con la que hemos realizado el análisis La data a usar del BCRP es la Producción agropecuaria (miles de toneladas) de Quinua desde los años 2007 al 2024.

Para realizar el ejemplo, comenzamos instalando y ejecutando las librerías necesarias:

library('forecast')
## Warning: package 'forecast' was built under R version 4.3.2
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library('ggplot2')
## Warning: package 'ggplot2' was built under R version 4.3.2
library('tseries')
## Warning: package 'tseries' was built under R version 4.3.2
library('seasonal')
## Warning: package 'seasonal' was built under R version 4.3.2
library('seasonalview')
## Warning: package 'seasonalview' was built under R version 4.3.2
## 
## Attaching package: 'seasonalview'
## The following object is masked from 'package:seasonal':
## 
##     view

PRIMERA ETAPA

Importamos la base de datos

library(readxl)
## Warning: package 'readxl' was built under R version 4.3.2
data <- read_excel("data.xlsx")
data
data.1=ts(data[,2], start = c(2007,1), frequency =12)

plot.ts(data.1, col="blue", xlab = "Mes/Año", ylab = "Miles de tonealdas",
        main = "Producción agropecuaria (miles de toneladas) - Agrícola - Quinua")

Podemos observar que la serie original cuenta con estacionalidad, por ende tenemos que removerlo.

Removemos la Estacionalidad

library(seasonal)
data.2=seas(data.1)
plot(data.2)

En la gráfica de rojo es la nueva serie que debemos obtener para ello haremos la descomposición de la serie original.

descomposicion=decompose(data.1, type = c("additive"), filter = NULL)
plot(descomposicion)

Podemos observar primero a la serie original y el trend que es la tendencia el cual se usara para remover esta estacionaliedad de la serie.

quinua=descomposicion$trend
plot(quinua, main = "Trend - Quinua")

Separamos la tendencia.

SEGUNDA ETAPA

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(na.omit(quinua))
## 
##  Augmented Dickey-Fuller Test
## 
## data:  na.omit(quinua)
## Dickey-Fuller = -2.9892, Lag order = 5, p-value = 0.1622
## alternative hypothesis: stationary

Contraste de hipótesis:

#Ho:Se asume que la serie temporal tiene una raíz unitaria. (No es estacionaria)

#H1: Se sugiere que serie temporal no tiene una raíz unitaria. (Es estacionaria)

Tras realizar la prueba aumentada de Dickey-Fuller (ADF), obtenemos un p-valor = 0.1713. 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(na.omit(quinua), lag.max=72,na.action=na.fail,col="blue",ci=0.95,ci.type="white",lwd=2,
    main = "Fac Serie temporal DIS - Producción agropecuaria (miles de toneladas) - Agrícola - Quinua al 95%")

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

pacf(na.omit(quinua), lag.max=72,na.action=na.pass,col="blue",ci=0.95,lwd=2,
     main = "FacP de Producción agropecuaria (miles de toneladas) - Agrícola - Quinua al 95%")

Para realizar un modelo ARIMA, la serie temporal debe ser estacionaria.

#Para conseguir esta estacionariedad, sacamos la primera diferencia a la tendencia de la serie.

pd_quinua = diff(quinua)
plot(pd_quinua)

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

adf.test(na.omit(pd_quinua))
## Warning in adf.test(na.omit(pd_quinua)): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  na.omit(pd_quinua)
## Dickey-Fuller = -4.1134, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

El 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(na.omit(pd_quinua), main='ACF for Differenced Series')

pacf(na.omit(pd_quinua), main='PACF for Differenced Series')

TERCERA ETAPA

auto.arima(pd_quinua, seasonal=FALSE)
## Series: pd_quinua 
## ARIMA(4,0,0) with zero mean 
## 
## Coefficients:
##          ar1      ar2     ar3      ar4
##       1.2488  -0.9594  0.6107  -0.1930
## s.e.  0.0710   0.1060  0.1056   0.0705
## 
## sigma^2 = 0.02241:  log likelihood = 92.38
## AIC=-174.77   AICc=-174.44   BIC=-158.53

Usando la notación ARIMA presentada anteriormente, el modelo ajustado se puede escribir como:

\[ \hat{Y_{d_t}} = 1.2488Y_{t-1} - 0.9594Y_{t-2} + 0.6107Y_{t-3} - 0.1930Y_{t-4} + \varepsilon \]

donde \(\varepsilon\) es un error y la serie original se diferencia con la orden 1.

CUARTA ETAPA

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(pd_quinua, seasonal=FALSE)
modeloarima
## Series: pd_quinua 
## ARIMA(4,0,0) with zero mean 
## 
## Coefficients:
##          ar1      ar2     ar3      ar4
##       1.2488  -0.9594  0.6107  -0.1930
## s.e.  0.0710   0.1060  0.1056   0.0705
## 
## sigma^2 = 0.02241:  log likelihood = 92.38
## AIC=-174.77   AICc=-174.44   BIC=-158.53
tsdisplay(residuals(modeloarima), lag.max=15, main='(4,0,0) Model Residuals')

Se observa que en ACF que se repite el rezago 4 sobrepasa el horizonte al igual que en el rezago 12 y eso debid a la serie diferenciada donde se tuvo un dato que estuvo por debajo de los demas(outliers). En el PACF se observa tambien que en los rezagos 4 y 12 estan fuera del límite.

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=24)
plot(prediccion)