Introducción

En este taller se escogen cuatro variables locales, sobre las cuales se realizara un procedimiento de extraccion de señales y prediccion mediante el modelo ARIMA.

Las variables son:

  1. td: tasa de desempleo en Cali (%)
  2. dah: deposito de ahorros en Cali (pesos)
  3. ipc: indice de precios del consumidor en Colombia (indice)
  4. xcol: exportaciones en colombia (miles de dolares FOB)

Las primeras dos variables se escogieron para analizar la situacion economica de las personas en Cali. Por otro lado, las dos ultimas variables se escogen para analizar la situacion economica desde el punto de vista monetario en Colombia.

Desarrollo

Importación de datos

en la primer seccion de codigo se importan las librerias y se cargan los datos del archivo de excel.

library(tseries)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr   1.1.1     ✔ readr   2.1.4
## ✔ forcats 1.0.0     ✔ stringr 1.5.0
## ✔ ggplot2 3.4.1     ✔ tibble  3.2.1
## ✔ purrr   1.0.1     ✔ tidyr   1.3.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## The following object is masked from 'package:purrr':
## 
##     some
library(astsa)
library(foreign)
library(timsac)
library(vars)
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## 
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## Loading required package: strucchange
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## Loading required package: sandwich
## 
## Attaching package: 'strucchange'
## 
## The following object is masked from 'package:stringr':
## 
##     boundary
## 
## Loading required package: urca
## Loading required package: lmtest
library(lmtest)
library(mFilter)
library(dynlm)
library(nlme)
## 
## Attaching package: 'nlme'
## 
## The following object is masked from 'package:dplyr':
## 
##     collapse
library(lmtest)
library(broom)
library(kableExtra)
## 
## Attaching package: 'kableExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     group_rows
library(knitr)
library(MASS)
library(parallel)
library(car)
library(mlogit)
## Loading required package: dfidx
## 
## Attaching package: 'dfidx'
## 
## The following object is masked from 'package:MASS':
## 
##     select
## 
## The following object is masked from 'package:stats':
## 
##     filter
library(dplyr)
library(tidyr)
library(forecast)
## 
## Attaching package: 'forecast'
## 
## The following object is masked from 'package:nlme':
## 
##     getResponse
## 
## The following object is masked from 'package:astsa':
## 
##     gas
library(fpp2)
## ── Attaching packages ────────────────────────────────────────────── fpp2 2.5 ──
## ✔ fma       2.5     ✔ expsmooth 2.3
## ── Conflicts ───────────────────────────────────────────────── fpp2_conflicts ──
## ✖ forecast::getResponse() masks nlme::getResponse()
## ✖ car::some()             masks purrr::some()
## 
## Attaching package: 'fpp2'
## 
## The following object is masked from 'package:astsa':
## 
##     oil
library(stats)
library(quantmod)
## Loading required package: xts
## 
## ################################### WARNING ###################################
## # We noticed you have dplyr installed. The dplyr lag() function breaks how    #
## # base R's lag() function is supposed to work, which breaks lag(my_xts).      #
## #                                                                             #
## # Calls to lag(my_xts) that you enter or source() into this session won't     #
## # work correctly.                                                             #
## #                                                                             #
## # All package code is unaffected because it is protected by the R namespace   #
## # mechanism.                                                                  #
## #                                                                             #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
## #                                                                             #
## # You can use stats::lag() to make sure you're not using dplyr::lag(), or you #
## # can add conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop   #
## # dplyr from breaking base R's lag() function.                                #
## ################################### WARNING ###################################
## 
## Attaching package: 'xts'
## 
## The following objects are masked from 'package:dplyr':
## 
##     first, last
## 
## Loading required package: TTR
library(readxl) 
library(seasonal) 
## 
## Attaching package: 'seasonal'
## 
## The following objects are masked from 'package:astsa':
## 
##     trend, unemp
## 
## The following object is masked from 'package:tibble':
## 
##     view
datos <- read_excel("C:/Users/migue/Desktop/taller 2/Datos Taller2.xlsx")

attach(datos)
names(datos)
## [1] "td"   "dah"  "ipc"  "xcol"

a continuacion se generan los objetos de series de tiempo mediante la funcion ts() y se grafican

td.ts = ts(td, start = c(2012,1), frequency = 12)
dah.ts = ts(dah, start = c(2012,1), frequency = 12)
ipc.ts = ts(ipc, start = c(2012,1), frequency = 12)
xcol.ts = ts(xcol, start = c(2012,1), frequency = 12)

autoplot(td.ts)+xlab("Tiempo")+ylab("%")+ggtitle("Tasa de desempleo en Cali")

autoplot(dah.ts)+xlab("Tiempo")+ylab("pesos")+ggtitle("Depositos de ahorro en Cali")

autoplot(ipc.ts)+xlab("Tiempo")+ylab("IPC")+ggtitle("Indice de Precios del consumidor en Colombia")

autoplot(xcol.ts)+xlab("Tiempo")+ylab("Miles de dolares FOB")+ggtitle("Exportaciones de Colombia")

Eleccion de modelo de descomposición

Se hace la descomposicion suponiendo un modelo aditivo usando la funcion decompose()

tdModAd = decompose(td.ts)
dahModAd = decompose(dah.ts)
ipcModAd = decompose(ipc.ts)
xcolModAd = decompose(xcol.ts)

plot(tdModAd)

plot(dahModAd)

plot(ipcModAd)

plot(xcolModAd)

tambien se hace la descomposicion de las series de tiempo suponiendo un modelo multiplicativo

# descomposicion para modelo multiplicativo
tdModMult = decompose(td.ts, type = "mult")
dahModMult = decompose(dah.ts, type = "mult")
ipcModMult = decompose(ipc.ts, type = "mult")
xcolModMult = decompose(xcol.ts, type = "mult")

plot(tdModMult)

plot(dahModMult)

plot(ipcModMult)

plot(xcolModMult)

se muestra la componente de tendencia individualmente

tdTendencia = tdModMult$trend
dahTendencia = dahModMult$trend
ipcTendencia = ipcModMult$trend
xcolTendencia = xcolModMult$trend



autoplot(tdTendencia)+xlab("Tiempo")+ylab("%")+ggtitle("Tendencia td")

autoplot(dahTendencia)+xlab("Tiempo")+ylab("pesos")+ggtitle("Tendencia dah")

autoplot(ipcTendencia)+xlab("Tiempo")+ylab("IPC")+ggtitle("Tendencia ipc")

autoplot(xcolTendencia)+xlab("Tiempo")+ylab("Miles de dolares FOB")+ggtitle("Tendencia xcol")

Analisis de tendencia

De la grafica de tasa de desempleo destaca un pico muy pronunciado en 2020, el cual esta totalmente relacionado con la pandemia del covid-19, la cual inicio ese año.

La grafida de depositos de ahorro muestra una tendencia alcista, en casi la totalidad del tiempo de la muestra.

La grafica de indice de precios tambien muestra una tendencia alcista, la cual es una caracteristica comun en un pais con una economia inflacionaria como lo es Colombia.

La grafica de exportaciones muestra algunos ciclos, es decir que las exportaciones han aumentado y disminuido a lo largo de los años de la muestra, estos ciclos pueden estar muy relacionados con la posicion del peso colombiano respecto a otras divisas.

se muestra la componente estacionaria individualmente

tdEstacionalidad = tdModMult$seasonal
dahEstacionalidad = dahModMult$seasonal
ipcEstacionalidad = ipcModMult$seasonal
xcolEstacionalidad = xcolModMult$seasonal

autoplot(tdEstacionalidad)+xlab("Tiempo")+ylab("%")+ggtitle("Estacionalidad td")

autoplot(dahEstacionalidad)+xlab("Tiempo")+ylab("pesos")+ggtitle("Estacionalidad dah")

autoplot(ipcEstacionalidad)+xlab("Tiempo")+ylab("IPC")+ggtitle("Estacionalidad ipc")

autoplot(xcolEstacionalidad)+xlab("Tiempo")+ylab("Miles de dolares FOB")+ggtitle("Estacionalidad xcol")

ggseasonplot(td.ts)+ggtitle("Tasa de desempleo en Cali agrupada por estaciones")

ggseasonplot(dah.ts)+ggtitle("Depositos de ahorro en Cali agrupado por estaciones")

ggseasonplot(ipc.ts)+ggtitle("IPC en Colombia agrupado por estaciones")

ggseasonplot(xcol.ts)+ggtitle("Exportaciones en Colombia agrupado por estaciones")

Analisis de estacionalidad

La estacionalidad en la tasa de desempleo muestra un crecimiento entre enero y marzo aproximadamente, para luego tener una tendencia a la baja hasta fin de año, el trazo atipico que se ve en la parte de arriba corresponde al 2020 donde el desempleo fue bastante amplio

En los depositos de ahorro se observa una tendencia alcista muy poco pronunciada, por lo cual la estacionalidad no afecta mucho esta variable

En el IPC se ve una tendencia alcista mucho mas estable, ya que es una variable que cada mes tiene a aumentar

Las exportaciones muestran una estacionalidad muy diferente en cada año, tiene muchas variaciones pero se puede percibir una tendencia a la alza a lo largo de los meses

se muestra la componente irregular individualmente

tdAleatorio = tdModMult$random
dahAleatorio = dahModMult$random
ipcAleatorio = ipcModMult$random
xcolAleatorio = xcolModMult$random

autoplot(tdAleatorio)+xlab("Tiempo")+ylab("%")+ggtitle("Aleatorio td")

autoplot(dahAleatorio)+xlab("Tiempo")+ylab("pesos")+ggtitle("Aleatorio dah")

autoplot(ipcAleatorio)+xlab("Tiempo")+ylab("IPC")+ggtitle("Aleatorio ipc")

autoplot(xcolAleatorio)+xlab("Tiempo")+ylab("Miles de dolares FOB")+ggtitle("Aleatorio xcol")

Analisis de varianza, tendencia-estacional

se muestra la combinacion de los componentes (tendencia*estacional) sobrepuestos sobre la tendencia

tdTendenciaCiclo = tdTendencia*tdEstacionalidad
dahTendenciaCiclo = dahTendencia*tdEstacionalidad
ipcTendenciaCiclo = ipcTendencia*tdEstacionalidad
xcolTendenciaCiclo = xcolTendencia*tdEstacionalidad

ts.plot(cbind(tdTendencia,tdTendenciaCiclo), lty = 1:2)

ts.plot(cbind(dahTendencia,dahTendenciaCiclo), lty = 1:2)

ts.plot(cbind(ipcTendencia,ipcTendenciaCiclo), lty = 1:2)

ts.plot(cbind(xcolTendencia,xcolTendenciaCiclo), lty = 1:2)

gracias a esta grafica se pueden apreciar las oscilaciones de cada variable al rededor de la tendencia, al identificar que estas oscilaciones cambian segun el tiempo, se puede considerar que la varianza de las series temporales son cambiantes, y por ende, el mejor modelo sera el multiplicativo.

Obtencion de parametro “d” y analisis de serie ajustada por estacionalidad

se estima la cantidad de diferenciaciones que se necesitaran para que las series sean estacionarias

Se realiza el test de Dickey-Fuller para verificar si las series de tiempo son estacionarias

#se ven las diferenciaciones
ndiffs(td.ts)
## [1] 0
ndiffs(dah.ts)
## [1] 1
ndiffs(ipc.ts)
## [1] 2
ndiffs(xcol.ts)
## [1] 1
#test de Dickey-Fuller
adf.test(td)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  td
## Dickey-Fuller = -2.5786, Lag order = 5, p-value = 0.3359
## alternative hypothesis: stationary
adf.test(dah)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  dah
## Dickey-Fuller = -1.6627, Lag order = 5, p-value = 0.7168
## alternative hypothesis: stationary
adf.test(ipc)
## Warning in adf.test(ipc): p-value greater than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  ipc
## Dickey-Fuller = 0.48669, Lag order = 5, p-value = 0.99
## alternative hypothesis: stationary
adf.test(xcol)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  xcol
## Dickey-Fuller = -1.5567, Lag order = 5, p-value = 0.7609
## alternative hypothesis: stationary

Ninguna serie es estacionaria originalmente por lo que se debe diferenciar, pero antes se revisa si la estacionalidad afecta mucho a la variable, de esta manera se decide por trabajar con datos originales o ajustados

#td
tdSa <- seasadj(tdModMult)
plot(td.ts)
lines(tdSa, col = "red")

observando la grafica se determina que se puede trabajar con los datos ajustados por estacionalidad, se procede a sacar la primer diferencia para lograr estacionariedad

# primer diferencia
tdDif = diff(tdSa)
autoplot(tdDif)+xlab("Tiempo")+ylab("%")+ggtitle("primera diferencia de td")

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

la prueba de dickey fuller muestra un valor p < 0.05 con una diferenciacion, por lo que se considera estacionaria con d = 1

se hace la diferenciacion y el test de Dickey Fuller con las demas series de tiempo para conocer el valor d

#dah
dahSa <- seasadj(dahModMult)
plot(dah.ts)
lines(dahSa, col = "red")

se comprueba que no se necesita ajustar por estacionalidad

# primer diferencia
dahDif = diff(dah.ts)
autoplot(dahDif)+xlab("Tiempo")+ylab("pesos")+ggtitle("primera diferencia de dah")

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

se logra la estacionariedad con d = 1

#ipc
ipcSa <- seasadj(ipcModMult)
plot(ipc.ts)
lines(ipcSa, col = "red")

se comprueba que no se necesita ajustar por estacionalidad

# dos difernecias
ipcDif = diff(ipc.ts)
ipcDif2 = diff(ipcDif)
autoplot(ipcDif2)+xlab("Tiempo")+ylab("IPC")+ggtitle("segunda diferencia de ipc")

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

se logra la estacionariedad con d = 2

#xcol
xcolSa <- seasadj(xcolModMult)
plot(xcol.ts)
lines(xcolSa, col = "red")

se comprueba que no se necesita ajustar por estacionalidad

# primer diferencia
xcolDif = diff(xcol.ts)
autoplot(xcolDif)+xlab("Tiempo")+ylab("Miles de dolares FOB")+ggtitle("primera diferencia de xcol")

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

se logra la estacionariedad con d = 1

Analisis de resultados contextualizados

antes de proceder con los modelos predictivos, se debe dar un espacio de analisis para cada una de las variables y una interpretacion en conjunto. 

Como se vio anteriormente, se tienen dos variables para la ciudad de Cali, las cuales estan muy relacionadas. El momento mas critico del desempleo en cali fue en 2020, debido a la crisis sanitaria y la cuarentena muchas empresas decidieron reducir su personal o simplemente no pudieron seguir operando dejando a muchas personas desempleadas. Por otro lado, los depositos de ahorro de los caleños continuaron en aumento e incluso lo hicieron mas rapido que años anteriores. Esto se debe a que las personas en una situacion de crisis prefieren resguardar su dinero y evitar el consumo de todos los productos que esten por fuera de la canasta basica. Hay que tener en cuenta la realidad del pais y del mundo para dicho periodo, ya que de ello depende en gran medida el analisis. Teniendo en cuenta lo anterior, observando el comportamiento de las graficas y situandonos en 2023 donde la pandemia ya ha quedado en un segundo plano y el desempleo ha disminuido notablemente, se puede suponer un escenario en el que las personas vuelven a tener ingresos consistentes y pueden darse mas oportunidad de gastar en lugar de ahorrar, esto significa que los depositos de ahorro pueden disminuir, tal como lo vienen haciendo desde el año pasado.

desde el punto de vista Nacional, se tienen la variable del Indice de Precios del Consumidor (IPC) y las exportaciones, estas variables se relacionan con la situacion monetaria del pais, sin embargo se deben tener muchas mas variables y elementos en cuenta para entender la situacion monetaria, al ver el comportamiento de las graficas nos damos cuenta que el IPC continua con su tendencia alcista e incluso con mayor velocidad de crecimiento que años anteriores, mientras que las exportaciones parece que vienen recuperandose de un declive en 2020, el pronostico de las exportaciones puede estar muy ligado a la tasa de cambio COP/USD ya que la mayoria de exportadores recibe ingresos en dolares mientras que el IPC se relaciona mas con las politicas monetarias del banco central las cuales pueden afectar tanto los costos de los productores como la demanda de bienes de canasta basica.

Obtencion de parametros “p” y “q”

Continuando con los modelos predictivos, se debe escoger q y p para cada modelo, gracias a las funciones de autocorrelacion acf() y pacf(), se puede determinar graficamente cuales son los posibles valores para p y q respectivamente, se hace observando los retrasos y cuales de ellos sobrepasan el umbral para ser considerados significativos

acf(tdDif)

Pacf(tdDif)

p = {0, 1} q = {0, 1, 2, 3}

posibles modelos ARIMA(p,d,q) para la variable td

  1. (0,1,0)
  2. (0,1,1)
  3. (0,1,2)
  4. (0,1,3)
  5. (1,1,0)
  6. (1,1,1)
  7. (1,1,2)
  8. (1,1,3)
#dah
acf(dahDif)

Pacf(dahDif)

p = {0, 1} q = {0, 1, 2, 3}

posibles modelos ARIMA(p,d,q) para la variable dah

  1. (0,1,0)
  2. (0,1,1)
  3. (0,1,2)
  4. (0,1,3)
  5. (1,1,0)
  6. (1,1,1)
  7. (1,1,2)
  8. (1,1,3)
#ipc
acf(ipcDif)

Pacf(ipcDif)

p = {0, 1} q = {0, 1, 2, 3}

posibles modelos ARIMA(p,d,q) para la variable ipc

  1. (0,2,0)
  2. (0,2,1)
  3. (0,2,2)
  4. (0,2,3)
  5. (1,2,0)
  6. (1,2,1)
  7. (1,2,2)
  8. (1,2,3)
#xcol
acf(xcolDif)

Pacf(xcolDif)

p = {0, 1, 2} q = {0, 1, 2}

posibles modelos ARIMA(p,d,q) para la variable xcol

  1. (0,1,0)
  2. (0,1,1)
  3. (0,1,2)
  4. (1,1,0)
  5. (1,1,1)
  6. (1,1,2)
  7. (2,1,0)
  8. (2,1,1)
  9. (2,1,2)

Generacion de modelos ARIMA

se calcula el modelo ARIMA optimo usando la funcion auto.arima() de esta manera se obtienen los coeficientes de los modelos que generan un menor AIC

#arima(td)
auto.arima(tdSa)
## Series: tdSa 
## ARIMA(1,0,2) with non-zero mean 
## 
## Coefficients:
##          ar1     ma1     ma2     mean
##       0.8592  0.7500  0.6842  13.9371
## s.e.  0.0464  0.0933  0.0681   1.0341
## 
## sigma^2 = 0.5381:  log likelihood = -146.55
## AIC=303.09   AICc=303.57   BIC=317.5
auto.arima(dah)
## Series: dah 
## ARIMA(0,1,0) 
## 
## sigma^2 = 5.839e+22:  log likelihood = -3619.48
## AIC=7240.97   AICc=7241   BIC=7243.84
auto.arima(ipc)
## Series: ipc 
## ARIMA(2,2,2) 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2
##       1.5702  -0.7954  -1.7911  0.8925
## s.e.  0.0821   0.0884   0.0488  0.0612
## 
## sigma^2 = 0.06906:  log likelihood = -9.47
## AIC=28.94   AICc=29.42   BIC=43.28
auto.arima(xcol)
## Series: xcol 
## ARIMA(0,1,1) 
## 
## Coefficients:
##           ma1
##       -0.5064
## s.e.   0.0731
## 
## sigma^2 = 1.248e+11:  log likelihood = -1859.04
## AIC=3722.09   AICc=3722.18   BIC=3727.84

se generan los modelos arima segun los indices obtenidos de la funcion auto.arima() y tambien modelos que, segun los criterios de las garficas de autocorrelacion, podrian funcionar bien

tdARIMA1 = Arima(tdSa,order = c(1,0,2))
dahARIMA1 = Arima(dah,order = c(0,1,0))
ipcARIMA1 = Arima(ipc,order = c(2,2,2))
xcolARIMA1 = Arima(xcol,order = c(0,1,1))


tdARIMA2 = Arima(tdSa,order = c(1,1,3))
dahARIMA2 = Arima(dah,order = c(1,1,3))
ipcARIMA2 = Arima(ipc,order = c(1,2,3))
xcolARIMA2 = Arima(xcol,order = c(2,1,2))

tdARIMA2
## Series: tdSa 
## ARIMA(1,1,3) 
## 
## Coefficients:
##          ar1      ma1      ma2      ma3
##       0.8742  -0.2555  -0.0639  -0.6807
## s.e.  0.0484   0.0949   0.0645   0.0690
## 
## sigma^2 = 0.5424:  log likelihood = -145.54
## AIC=301.08   AICc=301.56   BIC=315.46
dahARIMA2
## Series: dah 
## ARIMA(1,1,3) 
## 
## Coefficients:
##          ar1      ma1     ma2     ma3
##       0.3079  -0.0045  0.2164  0.0342
## s.e.  0.3896   0.3888  0.1270  0.1458
## 
## sigma^2 = 4.983e+22:  log likelihood = -3607.19
## AIC=7224.37   AICc=7224.85   BIC=7238.75
ipcARIMA2
## Series: ipc 
## ARIMA(1,2,3) 
## 
## Coefficients:
##          ar1      ma1      ma2      ma3
##       0.3910  -0.5860  -0.1859  -0.0666
## s.e.  0.3077   0.3145   0.1292   0.1574
## 
## sigma^2 = 0.07558:  log likelihood = -14.79
## AIC=39.58   AICc=40.06   BIC=53.91
xcolARIMA2
## Series: xcol 
## ARIMA(2,1,2) 
## 
## Coefficients:
##           ar1      ar2     ma1      ma2
##       -1.0566  -0.1102  0.5639  -0.4360
## s.e.   0.1746   0.1744  0.1648   0.1623
## 
## sigma^2 = 1.239e+11:  log likelihood = -1857.84
## AIC=3725.69   AICc=3726.17   BIC=3740.06

hay que destacar el hecho de que, para las variables td, dah y xcol, el AIC y el BIC son menores para el modelo escogido mediante las autocorrelaciones (modelo 2), que para el modelo automatico generado por el software (modelo 1). Solo para la variable ipc, el modelo automatico tuvo menor AIC y BIC. Sin embargo, la significancia algunas variables explicativas es un poco baja para el modelo 2, comparandolas con el modelo 1.

Analisis de componentes residuales de los modelos

A continuacion se realiza la preuba de box jenkins para identificar si los residuales son efectivamente ruido blanco

Box.test(tdARIMA1$residuals, lag = 20, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  tdARIMA1$residuals
## X-squared = 27.348, df = 20, p-value = 0.1258
Box.test(dahARIMA1$residuals, lag = 20, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  dahARIMA1$residuals
## X-squared = 46.834, df = 20, p-value = 0.0006185
Box.test(ipcARIMA1$residuals, lag = 20, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  ipcARIMA1$residuals
## X-squared = 32.202, df = 20, p-value = 0.04119
Box.test(xcolARIMA1$residuals, lag = 20, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  xcolARIMA1$residuals
## X-squared = 27.531, df = 20, p-value = 0.121
Box.test(tdARIMA2$residuals, lag = 20, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  tdARIMA2$residuals
## X-squared = 27.796, df = 20, p-value = 0.1143
Box.test(dahARIMA2$residuals, lag = 20, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  dahARIMA2$residuals
## X-squared = 18.687, df = 20, p-value = 0.5423
Box.test(ipcARIMA2$residuals, lag = 20, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  ipcARIMA2$residuals
## X-squared = 43.724, df = 20, p-value = 0.001639
Box.test(xcolARIMA2$residuals, lag = 20, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  xcolARIMA2$residuals
## X-squared = 22.063, df = 20, p-value = 0.3371

para el modelo 1, las variables td y xcol pasan la prueba dando un valor p > 0.05, demostrando que no hay correlacion en los residuos. Para el modelo 2, la unica variable que no pasa la prueba es ipc, la demas muestras residuos que no se correlacionan

Ahora se hace la prueba de shapiro para identifcar si los residuos se distribuyen de forma normal

shapiro.test(tdARIMA1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  tdARIMA1$residuals
## W = 0.86521, p-value = 1.319e-09
shapiro.test(dahARIMA1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  dahARIMA1$residuals
## W = 0.89709, p-value = 4.519e-08
shapiro.test(ipcARIMA1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  ipcARIMA1$residuals
## W = 0.96317, p-value = 0.001219
shapiro.test(xcolARIMA1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  xcolARIMA1$residuals
## W = 0.99382, p-value = 0.8386
shapiro.test(tdARIMA2$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  tdARIMA2$residuals
## W = 0.86385, p-value = 1.148e-09
shapiro.test(dahARIMA2$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  dahARIMA2$residuals
## W = 0.90444, p-value = 1.12e-07
shapiro.test(ipcARIMA2$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  ipcARIMA2$residuals
## W = 0.95637, p-value = 0.0003217
shapiro.test(xcolARIMA2$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  xcolARIMA2$residuals
## W = 0.99312, p-value = 0.7725

Para ambos modelos, las variables td, dah e ipc muestran que los residuos no se distribuyen de manera normal, mientras que para la variable xcol, en ambos modelos dio un valor p > 0.05, lo cual indica que los residuos tienen un comportamiento normal

Pronostico

se generan los pronosticos usando el modelo 1 definido por el software. Se pronostican 6 datos adicionales, es decir, 6 meses y luego se grafica la serie con el pronostico

tdPronostico1=forecast(tdARIMA1,level= c(95), h=6)
dahPronostico1=forecast(dahARIMA1,level= c(95), h=6)
ipcPronostico1=forecast(ipcARIMA1,level= c(95), h=6)
xcolPronostico1=forecast(xcolARIMA1,level= c(95), h=6)

autoplot(tdPronostico1)+xlab("Tiempo")+ylab("%")+ggtitle("Pronostico td - Modelo 1")

autoplot(dahPronostico1)+xlab("Tiempo")+ylab("pesos")+ggtitle("Pronostico dah - Modelo 1")

autoplot(ipcPronostico1)+xlab("Tiempo")+ylab("IPC")+ggtitle("Pronostico ipc - Modelo 1")

autoplot(xcolPronostico1)+xlab("Tiempo")+ylab("miles de dolares FOB")+ggtitle("Pronostico xcol - Modelo 1")

tdPron1.m2 = predict(tdARIMA1, n.ahead=6) 
dahPron1.m2 = predict(dahARIMA1, n.ahead=6) 
ipcPron1.m2 = predict(ipcARIMA1, n.ahead=6) 
xcolPron1.m2 = predict(xcolARIMA1, n.ahead=6)

se usa el modelo 2 para pronosticar 6 datos adicionales correspondientes a 6 meses y se grafican

tdPronostico2=forecast(tdARIMA2,level= c(95), h=6)
dahPronostico2=forecast(dahARIMA2,level= c(95), h=6)
ipcPronostico2=forecast(ipcARIMA2,level= c(95), h=6)
xcolPronostico2=forecast(xcolARIMA2,level= c(95), h=6)

autoplot(tdPronostico2)+xlab("Tiempo")+ylab("%")+ggtitle("Pronostico td - Modelo 2")

autoplot(dahPronostico2)+xlab("Tiempo")+ylab("pesos")+ggtitle("Pronostico dah - Modelo 2")

autoplot(ipcPronostico2)+xlab("Tiempo")+ylab("IPC")+ggtitle("Pronostico ipc - Modelo 2")

autoplot(xcolPronostico2)+xlab("Tiempo")+ylab("miles de dolares FOB")+ggtitle("Pronostico xcol - Modelo 2")

tdPron2.m2 = predict(tdARIMA2, n.ahead=6) 
dahPron2.m2 = predict(dahARIMA2, n.ahead=6) 
ipcPron2.m2 = predict(ipcARIMA2, n.ahead=6) 
xcolPron2.m2 = predict(xcolARIMA2, n.ahead=6)

Medicion del error

para medir el error se usan las metricas del Mean Absolute Error (MAE) y la Root Mean Square Deviation (RMSD) las cuales miden que tan alejados estan los pronosticos de la realidad

tdPronostico1 = round(tdPron1.m2$pred,6)
dahPronostico1 = round(dahPron1.m2$pred,6)
ipcPronostico1 = round(ipcPron1.m2$pred,6)
xcolPronostico1 = round(xcolPron1.m2$pred,6)

tdReales1 = datos[127:132,1]
dahReales1 = datos[127:132,2]
ipcReales1 = datos[127:132,3]
xcolReales1 = datos[127:132,4]

print(tdPronostico1)
##           Jan      Feb      Mar      Apr      May      Jun
## 2023 10.77929 10.81217 11.25215 11.63018 11.95499 12.23406
print(dahPronostico1)
## Time Series:
## Start = 133 
## End = 138 
## Frequency = 1 
## [1] 1.165763e+13 1.165763e+13 1.165763e+13 1.165763e+13 1.165763e+13
## [6] 1.165763e+13
print(ipcPronostico1)
## Time Series:
## Start = 133 
## End = 138 
## Frequency = 1 
## [1] 127.6401 129.2196 130.7190 132.1171 133.4198 134.6534
print(xcolPronostico1)
## Time Series:
## Start = 133 
## End = 138 
## Frequency = 1 
## [1] 3433769 3433769 3433769 3433769 3433769 3433769
print(tdReales1)
## # A tibble: 6 × 1
##      td
##   <dbl>
## 1 11.9 
## 2 10.7 
## 3 11.0 
## 4 10.0 
## 5 10.3 
## 6  9.39
print(dahReales1)
## # A tibble: 6 × 1
##       dah
##     <dbl>
## 1 1.22e13
## 2 1.20e13
## 3 1.18e13
## 4 1.17e13
## 5 1.18e13
## 6 1.17e13
print(ipcReales1)
## # A tibble: 6 × 1
##     ipc
##   <dbl>
## 1  120.
## 2  122.
## 3  123.
## 4  124.
## 5  124.
## 6  126.
print(xcolReales1)
## # A tibble: 6 × 1
##       xcol
##      <dbl>
## 1 4359502.
## 2 3445402.
## 3 3587364.
## 4 3136241.
## 5 3379029.
## 6 3473998.
tdPronosticoN1 = c(10.779, 10.812, 11.252, 11.63018, 11.95499, 12.23406)
dahPronosticoN1 = c(1.165763e+13, 1.165763e+13, 1.165763e+13, 1.165763e+13, 1.165763e+13, 1.165763e+13)
ipcPronosticoN1 = c(127.640, 129.220, 130.719, 132.1171, 133.4198, 134.6534)
xcolPronosticoN1 = c(3433769, 3433769, 3433769, 3433769, 3433769, 3433769)


tdRealesN1 = c(10.011868, 10.255740, 9.388159, 10.011868, 10.255740, 9.388159   )
dahRealesN1 = c(1.174042e+13, 1.178261e+13, 1.165763e+13, 1.174042e+13, 1.178261e+13, 1.165763e+13)
ipcRealesN1 = c(123.51, 124.46, 126.03, 123.51, 124.46, 126.03)
xcolRealesN1 = c(3136241, 3379029, 3473998, 3136241, 3379029, 3473998   )

accuracy(tdARIMA1)
##                        ME      RMSE       MAE       MPE     MAPE     MASE
## Training set -0.009158904 0.7223522 0.4751656 -0.287231 3.296839 0.175502
##                     ACF1
## Training set -0.00168021
tdARIMA1
## Series: tdSa 
## ARIMA(1,0,2) with non-zero mean 
## 
## Coefficients:
##          ar1     ma1     ma2     mean
##       0.8592  0.7500  0.6842  13.9371
## s.e.  0.0464  0.0933  0.0681   1.0341
## 
## sigma^2 = 0.5381:  log likelihood = -146.55
## AIC=303.09   AICc=303.57   BIC=317.5
MAEtd1 <- mean(abs(tdPronosticoN1 - tdRealesN1))
MAEdah1 <- mean(abs(dahPronosticoN1 - dahRealesN1))
MAEipc1 <- mean(abs(ipcPronosticoN1 - ipcRealesN1))
MAExcol1 <- mean(abs(xcolPronosticoN1 - xcolRealesN1))


RMSEtd1 <- sqrt(mean((tdPronosticoN1 - tdRealesN1)^2))
RMSEdah1 <- sqrt(mean((dahPronosticoN1 - dahRealesN1)^2))
RMSEipc1 <- sqrt(mean((ipcPronosticoN1 - ipcRealesN1)^2))
RMSExcol1 <- sqrt(mean((xcolPronosticoN1 - xcolRealesN1)^2))

print(MAEtd1)
## [1] 1.558449
print(MAEdah1)
## [1] 69256666667
print(MAEipc1)
## [1] 6.628217
print(MAExcol1)
## [1] 130832.3
print(RMSEtd1)
## [1] 1.730961
print(RMSEdah1)
## [1] 86552844167
print(RMSEipc1)
## [1] 6.957304
print(RMSExcol1)
## [1] 176198.5

Se obtienen los resultados del MAE y el RMSE para el modelo 1 usando los valores observados y los pronosticados

Se procede a obtener las metricas del modelo 2 para ser comparadas con las del modelo 1

tdPronostico2 = round(tdPron2.m2$pred,6)
dahPronostico2 = round(dahPron2.m2$pred,6)
ipcPronostico2 = round(ipcPron2.m2$pred,6)
xcolPronostico2 = round(xcolPron2.m2$pred,6)

tdReales2 = datos[127:132,1]
dahReales2 = datos[127:132,2]
ipcReales2 = datos[127:132,3]
xcolReales2 = datos[127:132,4]

print(tdPronostico2)
##           Jan      Feb      Mar      Apr      May      Jun
## 2023 10.75454 10.75070 11.15050 11.49999 11.80550 12.07257
print(dahPronostico2)
## Time Series:
## Start = 133 
## End = 138 
## Frequency = 1 
## [1] 1.164340e+13 1.161081e+13 1.159577e+13 1.159114e+13 1.158972e+13
## [6] 1.158928e+13
print(ipcPronostico2)
## Time Series:
## Start = 133 
## End = 138 
## Frequency = 1 
## [1] 127.4729 128.7450 129.9091 131.0310 132.1365 133.2354
print(xcolPronostico2)
## Time Series:
## Start = 133 
## End = 138 
## Frequency = 1 
## [1] 3449835 3408004 3454865 3409960 3452244 3412514
print(tdReales2)
## # A tibble: 6 × 1
##      td
##   <dbl>
## 1 11.9 
## 2 10.7 
## 3 11.0 
## 4 10.0 
## 5 10.3 
## 6  9.39
print(dahReales2)
## # A tibble: 6 × 1
##       dah
##     <dbl>
## 1 1.22e13
## 2 1.20e13
## 3 1.18e13
## 4 1.17e13
## 5 1.18e13
## 6 1.17e13
print(ipcReales2)
## # A tibble: 6 × 1
##     ipc
##   <dbl>
## 1  120.
## 2  122.
## 3  123.
## 4  124.
## 5  124.
## 6  126.
print(xcolReales2)
## # A tibble: 6 × 1
##       xcol
##      <dbl>
## 1 4359502.
## 2 3445402.
## 3 3587364.
## 4 3136241.
## 5 3379029.
## 6 3473998.
tdPronosticoN2 = c(110.75454, 10.75070, 11.15050, 11.49999, 11.80550, 12.07257)
dahPronosticoN2 = c(1.164340e+13, 1.161081e+13, 1.159577e+13, 1.159114e+13, 1.158972e+13, 1.158928e+13)
ipcPronosticoN2 = c(127.4729, 128.7450, 129.9091, 131.0310, 132.1365, 133.2354)
xcolPronosticoN2 = c(3449835, 3408004, 3454865, 3409960, 3452244, 3412514)


tdRealesN2 = c(11.863820,               10.707232,              11.039690,              10.011868,              10.255740,              9.388159    )
dahRealesN2 = c(1.219626e+13,               1.202048e+13,               1.179142e+13,               1.174042e+13,               1.178261e+13, 1.165763e+13  )
ipcRealesN2 = c(120.27,             121.50,             122.63,             123.51,         124.46,             126.03  )
xcolRealesN2 = c(4359502,               3445402,                3587364,                3136241,                3379029,                3473998 )

accuracy(tdARIMA2)
##                       ME      RMSE      MAE        MPE     MAPE      MASE
## Training set -0.02456326 0.7223669 0.470676 -0.3558722 3.273949 0.1738437
##                      ACF1
## Training set -0.009969689
tdARIMA2
## Series: tdSa 
## ARIMA(1,1,3) 
## 
## Coefficients:
##          ar1      ma1      ma2      ma3
##       0.8742  -0.2555  -0.0639  -0.6807
## s.e.  0.0484   0.0949   0.0645   0.0690
## 
## sigma^2 = 0.5424:  log likelihood = -145.54
## AIC=301.08   AICc=301.56   BIC=315.46
MAEtd2 <- mean(abs(tdPronosticoN2 - tdRealesN2))
MAEdah2 <- mean(abs(dahPronosticoN2 - dahRealesN2))
MAEipc2 <- mean(abs(ipcPronosticoN2 - ipcRealesN2))
MAExcol2 <- mean(abs(xcolPronosticoN2 - xcolRealesN2))


RMSEtd2 <- sqrt(mean((tdPronosticoN2 - tdRealesN2)^2))
RMSEdah2 <- sqrt(mean((dahPronosticoN2 - dahRealesN2)^2))
RMSEipc2 <- sqrt(mean((ipcPronosticoN2 - ipcRealesN2)^2))
RMSExcol2 <- sqrt(mean((xcolPronosticoN2 - xcolRealesN2)^2))

print(MAEtd2)
## [1] 17.46122
print(MAEdah2)
## [1] 2.6145e+11
print(MAEipc2)
## [1] 7.354983
print(MAExcol2)
## [1] 247997
print(RMSEtd2)
## [1] 40.39639
print(RMSEdah2)
## [1] 309818603917
print(RMSEipc2)
## [1] 7.357184
print(RMSExcol2)
## [1] 393808.6

Conclusion

Se obtienen los resultados de las metricas para el modelo 2. Para finalizar, se comparan los resultados y se observa que para todas las cuatro variables, el modelo 1 presenta un menor MAE y RMSE, respecto al modelo 2, en algunos casos con diferencias considerables. Esto significa que el modelo 1 introduce menor error en el pronostico. Teniendo en cuenta que en en terminos del comportamiento de los residuos y los valores de AIC, ambos modelos tienen comportamientos similares, se escoge el modelo 1 para predecir todas las variables, dando la ventaja de un error minimo en los pronosticos, siendo mas cercano a la realidad.