El modelo ARIMAX es similar a un modelo de regresión multivariante, pero permite aprovechar la autocorrelación que puede estar presente en los residuos de la regresión para mejorar la precisión de un pronóstico.
El conjunto de datos contiene las siguientes variables:
- Consumo de helado en los EE. UU. (en pintas, per cápita),
- Ingreso familiar promedio por semana (en USD),
- Precio del helado (por pinta), y
- Temperatura media (en grados Fahrenheit). El número de observaciones es 30. Corresponden a períodos de cuatro semanas en el lapso del 18 de marzo de 1951 al 11 de julio de 1953
Cargamos el conjunto de datos y graficamod las variables
cons(consumo de helado), temp(temperatura) y
income.
library(ggplot2)
library(gridExtra)
library(DT)getwd()## [1] "/cloud/project/Series/SARIMAX"df <- read.csv("/cloud/project/Series/SARIMAX/Icecream.csv")
df%>%DT::datatable()p1 <- ggplot(df, aes(x = X, y = cons)) +
ylab("Consumption") +
xlab("") +
geom_line() +
expand_limits(x = 0, y = 0)
p2 <- ggplot(df, aes(x = X, y = temp)) +
ylab("Temperature") +
xlab("") +
geom_line() +
expand_limits(x = 0, y = 0)
p3 <- ggplot(df, aes(x = X, y = income)) +
ylab("Income") +
xlab("Period") +
geom_line() +
expand_limits(x = 0, y = 0)
grid.arrange(p1, p2, p3, ncol=1, nrow=3)
Estimamos un modelo ARIMA para los datos sobre el consumo de helado
usando la función auto.arima. Luego, pase el modelo como
entrada a la forecastfunción para obtener un pronóstico para los
próximos 6 períodos (ambas funciones son del forecastpaquete).
library(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoofit_cons <- auto.arima(df$cons)
fcast_cons <- forecast(fit_cons, h = 6)Trace el pronóstico obtenido con la autoplot.forecastfunción del forecastpaquete
library(forecast)
#autoplot.forecast(fcast_cons)
plot(fcast_cons)
Utilice la función accuracy del paquete
forecast para encontrar el error escalado absoluto medio
(MASE) del modelo ARIMA ajustado.
accuracy(fit_cons)## ME RMSE MAE MPE MAPE MASE
## Training set 3.438782e-05 0.03389421 0.02475857 -0.8273305 6.629011 0.7542003
## ACF1
## Training set 0.03907652Estime un modelo ARIMA extendido para los datos de consumo con la
variable de temperatura como regresor adicional (usando la función
auto.arima). Luego haga un pronóstico para los próximos 6
períodos (tenga en cuenta que este pronóstico requiere una suposición
sobre la temperatura esperada; suponga que la temperatura para los
próximos 6 períodos estará representada por el siguiente vector:)
fcast_temp <- c(70.5, 66, 60.5, 45.5, 36, 28). Trace el
pronóstico obtenido.
fit_cons_temp <- auto.arima(df$cons, xreg = df$temp)
fcast_temp <- c(70.5, 66, 60.5, 45.5, 36, 28)
fcast_cons_temp <- forecast(fit_cons_temp, xreg = fcast_temp, h = 6)
plot(fcast_cons_temp)
Imprimir resumen del pronóstico obtenido. Encuentre el coeficiente de la variable de temperatura, su error estándar y el MASE del pronóstico. Compare el MASE con el de la previsión inicial.
summary(fcast_cons_temp)##
## Forecast method: Regression with ARIMA(0,1,0) errors
##
## Model Information:
## Series: df$cons
## Regression with ARIMA(0,1,0) errors
##
## Coefficients:
## xreg
## 0.0028
## s.e. 0.0007
##
## sigma^2 = 0.001108: log likelihood = 58.03
## AIC=-112.06 AICc=-111.6 BIC=-109.32
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.002563685 0.03216453 0.02414157 0.564013 6.478971 0.7354048
## ACF1
## Training set -0.1457977
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 31 0.5465774 0.5039101 0.5892446 0.4813234 0.6118313
## 32 0.5337735 0.4734329 0.5941142 0.4414905 0.6260566
## 33 0.5181244 0.4442225 0.5920263 0.4051012 0.6311476
## 34 0.4754450 0.3901105 0.5607796 0.3449371 0.6059529
## 35 0.4484147 0.3530078 0.5438217 0.3025024 0.5943270
## 36 0.4256524 0.3211393 0.5301654 0.2658135 0.5854913Compruebe la significación estadística del coeficiente de la variable de temperatura utilizando la coeftestfunción del lmtestpaquete. ¿Es el coeficiente estadísticamente significativo al nivel del 5%
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericcoeftest(fit_cons_temp)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## xreg 0.0028453 0.0007302 3.8966 9.756e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1La función que estima el modelo ARIMA puede ingresar más regresores adicionales, pero solo en forma de matriz. Cree una matriz con las siguientes columnas:
-valores de la variable de temperatura, -valores de la variable renta, -los valores de la variable de ingresos se retrasaron un período, -los valores de la variable renta se retrasaron dos períodos.
Imprime la matriz.
Nota: las últimas tres columnas se pueden crear
anteponiendo dos NA al vector de valores de la variable de
ingresos y utilizando el vector obtenido como entrada a la función
embed (con el dimensionparámetro igual al número de
columnas a crear).
temp_column <- matrix(df$temp, ncol = 1)
income <- c(NA, NA, df$income)
income_matrix <- embed(income, 3)
vars_matrix <- cbind(temp_column, income_matrix)
print(vars_matrix)## [,1] [,2] [,3] [,4]
## [1,] 41 78 NA NA
## [2,] 56 79 78 NA
## [3,] 63 81 79 78
## [4,] 68 80 81 79
## [5,] 69 76 80 81
## [6,] 65 78 76 80
## [7,] 61 82 78 76
## [8,] 47 79 82 78
## [9,] 32 76 79 82
## [10,] 24 79 76 79
## [11,] 28 82 79 76
## [12,] 26 85 82 79
## [13,] 32 86 85 82
## [14,] 40 83 86 85
## [15,] 55 84 83 86
## [16,] 63 82 84 83
## [17,] 72 80 82 84
## [18,] 72 78 80 82
## [19,] 67 84 78 80
## [20,] 60 86 84 78
## [21,] 44 85 86 84
## [22,] 40 87 85 86
## [23,] 32 94 87 85
## [24,] 27 92 94 87
## [25,] 28 95 92 94
## [26,] 33 96 95 92
## [27,] 41 94 96 95
## [28,] 52 96 94 96
## [29,] 64 91 96 94
## [30,] 71 90 91 96Use la matriz obtenida para ajustar tres modelos ARIMA extendidos que usan las siguientes variables como regresores adicionales:
- temperatura, ingresos,
- temperatura, ingresos en rezagos 0, 1,
- temperatura, ingresos en rezagos 0, 1, 2. Examine el resumen de cada modelo y busque el modelo con el valor más bajo del criterio de información de Akaike (AIC).
Tenga en cuenta que el AIC no se puede utilizar para comparar modelos ARIMA con diferentes órdenes de integración (expresados por los términos intermedios en las especificaciones del modelo) debido a una diferencia en el número de observaciones. Por ejemplo, un valor AIC de un modelo no diferenciado, ARIMA (p, 0, q), no se puede comparar con el valor correspondiente de un modelo diferenciado, ARIMA (p, 1, q).
fit_vars_0 <- auto.arima(df$cons, xreg = vars_matrix[, 1:2])
fit_vars_1 <- auto.arima(df$cons, xreg = vars_matrix[, 1:3])
fit_vars_2 <- auto.arima(df$cons, xreg = vars_matrix[, 1:4])
print(fit_vars_0$aic)## [1] -113.3357print(fit_vars_1$aic)## [1] -111.9228print(fit_vars_2$aic)## [1] -110.2497Utilice el modelo que se encontró en el ejercicio anterior para hacer
un pronóstico para los próximos 6 períodos y grafique el pronóstico. (El
pronóstico requiere una matriz de la temperatura y los ingresos
esperados para los próximos 6 períodos; cree la matriz utilizando la
variable fcast_temp y los siguientes valores para los
ingresos esperados:) 91, 91, 93, 96, 96, 96. Encuentre el error escalado
absoluto medio del modelo y compárelo con los de los dos primeros
modelos de este conjunto de ejercicios
expected_temp_income <- matrix(c(fcast_temp, 91, 91, 93, 96, 96, 96), ncol = 2, nrow = 6)
fcast_cons_temp_income <- forecast(fit_vars_0, xreg = expected_temp_income, h = 6)
#autoplot.forecast(fcast_cons_temp_income)
plot(fcast_cons_temp_income)
accuracy(fit_cons)[, "MASE"]## [1] 0.7542003## [1] 0.8200619
accuracy(fit_cons_temp)[, "MASE"]## [1] 0.7354048## [1] 0.7354048
accuracy(fit_vars_0)[, "MASE"]## [1] 0.7290753## [1] 0.7290753
# the model with two external regressors has the lowest
# mean absolute scaled error (0.7290753)