library(forecast)
library(tseries)
library(TSA)
library(urca)
library(ggplot2)
library(dplyr)
library(stats)
library(seasonal)
library(readxl)

Analisis de extraccion de señales

Se seleccionaron como variables las importaciones y exportaciones totales a nivel nacional, debido a su relación con el comercio internacional

Preparacion de los datos

Se seleccionan los primeros 145 datos de un total de 150 para comparar los pronosticos

data <- read_excel("D:/Uni/PEBD/taller2Data.xlsx")
m <- ts(data$m[1:145], frequency=12, start=c(2012, 1))
x <- ts(data$x[1:145], frequency=12, start=c(2012, 1))

autoplot(m)

autoplot(x)

Ajuste Estacional

Para las importaciones:

#Importaciones
decompM <- decompose(m)
sasM <- seasadj(decompM)
sasMDiff1 <- diff(m, differences=1)
autoplot(sasM)

autoplot(sasMDiff1)

Para las exportaciones:

#Exportaciones
decompX <- decompose(x)
sasX <- seasadj(decompX)
sasXDiff1 <- diff(x, differences=1)
autoplot(sasX)

autoplot(sasXDiff1)

Los datos muestran que es posible usar tanto la serie original como la ajustada, en adelante se usara la serie ajustada para los pronosticos

Pronostico utilizando modelo ARIMA

Identificacion del modelo

Segun las graficas, podemos concluir que d=1 es suficiente para el modelo

Se evaluan los retrasos para determinar p y q:

#Importaciones
Pacf(sasMDiff1)

Acf(sasMDiff1)

#Exportaciones
Pacf(sasXDiff1)

Acf(sasXDiff1)

Estimacion de los modelos

De acuerdo con los graficos, Se podrian evaluar los siguientes modelos para m:

(1, 1, 1) (1, 1, 4) (2, 1, 1) (2, 1, 4) (4, 1, 1) (4, 1, 4)

Se usara el modelo (1, 1, 1) para comparar con el automatico.

Además, se podrian evaluar los siguientes modelos para q:

(1, 1, 1) (1, 1, 5) (1, 1, 6) (2, 1, 1) (2, 1, 5) (2, 1, 6) (6, 1, 1) (6, 1, 5) (6, 1, 6)

Tambien se usara el modelo (1, 1, 1) para comparar con el automatico.

#Importaciones
mAuto = auto.arima(sasM)
mAuto
## Series: sasM 
## ARIMA(0,1,1) 
## 
## Coefficients:
##           ma1
##       -0.5433
## s.e.   0.0647
## 
## sigma^2 = 1.106e+11:  log likelihood = -2034.91
## AIC=4073.82   AICc=4073.91   BIC=4079.76
#Exportaciones
xAuto = auto.arima(sasX)
xAuto
## Series: sasX 
## ARIMA(0,1,1) 
## 
## Coefficients:
##           ma1
##       -0.3747
## s.e.   0.0743
## 
## sigma^2 = 1.224e+11:  log likelihood = -2042.11
## AIC=4088.21   AICc=4088.3   BIC=4094.15

Validacion

#Importacions
mModel = Arima(sasM,order = c(1,1,1))
Box.test(mModel$residuals, lag = 20, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  mModel$residuals
## X-squared = 19.138, df = 20, p-value = 0.5129
xModel = Arima(sasX,order = c(1,1,1))
Box.test(xModel$residuals, lag = 20, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  xModel$residuals
## X-squared = 12.568, df = 20, p-value = 0.8952

Este test confirma que los residuos son ruido blanco

Pronostico

Se hace un pronostico de 2 meses para ambas variables usando el modelo automatico y el seleccionado (1, 1, 1)

mAutoForecast = forecast(mAuto,level= c(95), h=7)
mModelForecast = forecast(mModel,level= c(95), h=7)
plot(mAutoForecast)

plot(mModelForecast)

mAutoForecast
##          Point Forecast   Lo 95   Hi 95
## Feb 2024        5069531 4417671 5721390
## Mar 2024        5069531 4352903 5786158
## Apr 2024        5069531 4293522 5845539
## May 2024        5069531 4238373 5900689
## Jun 2024        5069531 4186662 5952400
## Jul 2024        5069531 4137816 6001246
## Aug 2024        5069531 4091407 6047655
mModelForecast
##          Point Forecast   Lo 95   Hi 95
## Feb 2024        5091177 4440463 5741891
## Mar 2024        5070594 4369106 5772082
## Apr 2024        5074336 4297956 5850715
## May 2024        5073655 4233975 5913336
## Jun 2024        5073779 4174413 5973145
## Jul 2024        5073757 4118569 6028944
## Aug 2024        5073761 4065814 6081707
mAutoPreds = c(5069531, 5069531, 5069531, 5069531, 5069531)
mModelPreds = c(5091177, 5070594, 5074336, 5073655, 5073779)
mOgVals = data$m[146:150]
xAutoForecast = forecast(xAuto,level= c(95), h=7)
xModelForecast = forecast(xModel,level= c(95), h=7)
plot(xAutoForecast)

plot(xModelForecast)

xAutoForecast
##          Point Forecast   Lo 95   Hi 95
## Feb 2024        4015680 3329945 4701414
## Mar 2024        4015680 3206931 4824428
## Apr 2024        4015680 3100302 4931057
## May 2024        4015680 3004859 5026500
## Jun 2024        4015680 2917681 5113678
## Jul 2024        4015680 2836933 5194426
## Aug 2024        4015680 2761373 5269986
xModelForecast
##          Point Forecast   Lo 95   Hi 95
## Feb 2024        4017568 3329462 4705674
## Mar 2024        4014990 3205037 4824942
## Apr 2024        4015047 3096583 4933512
## May 2024        4015046 2999655 5030437
## Jun 2024        4015046 2911206 5118886
## Jul 2024        4015046 2829337 5200755
## Aug 2024        4015046 2752767 5277326
xAutoPreds = c(4015680, 4015680, 4015680, 4015680, 4015680)
xModelPreds = c(4017568, 4014990, 4015047, 4015046, 4015046)
xOgVals = data$x[146:150]

MAE y RMSE

#Importaciones
mAutoMAE <- mean(abs(mAutoPreds - mOgVals))
mModelMAE <- mean(abs(mModelPreds - mOgVals))
mAutoRMSE <- sqrt(mean((mAutoPreds - mOgVals)^2))
mModelRMSE <- sqrt(mean((mModelPreds - mOgVals)^2))
mAutoMAE
## [1] 441278.2
mModelMAE
## [1] 444883.8
mAutoRMSE
## [1] 476843.2
mModelRMSE
## [1] 477723.9
#Exportaciones
xAutoMAE <- mean(abs(xAutoPreds - xOgVals))
xModelMAE <- mean(abs(xModelPreds - xOgVals))
xAutoRMSE <- sqrt(mean((xAutoPreds - xOgVals)^2))
xModelRMSE <- sqrt(mean((xModelPreds - xOgVals)^2))
xAutoMAE
## [1] 275264.9
xModelMAE
## [1] 275631.1
xAutoRMSE
## [1] 298595
xModelRMSE
## [1] 299055.3

Para ambas variables, los dos indicadores muestran que el modelo automatico tiene mejor capacidad de prediccion