Datos

#install.packages("tsDyn")
library(tsDyn)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

library(readxl)
library(ggplot2)
library(tseries)
library(forecast)
#file.choose()

Data = read_excel("C:\\Users\\Camilo\\Downloads\\Proyecto seminario\\DATOS LIMPIOS.xlsx", sheet = "Datos para trabajar en R")

## New names:
## • `` -> `...6`
## • `` -> `...7`
## • `` -> `...8`
## • `` -> `...9`

set.seed(123)
Datos = Data[10:11]
# Eliminar los Na de los datos 
Datos = na.omit(Datos) 

plot(Datos$`Cantidad de empresas registratas por trimestre`,type = "l")

La gráfica revela una tendencia creciente en los datos, sugiriendo además una posible estacionariedad. Con el objetivo de obtener predicciones más robustas, dividiremos la base de datos asignando una porción mayor para entrenamiento y reservando una parte para evaluar el rendimiento predictivo de los modelos propuestos.”

Propuesta de modelos

Datos1 = ts(Datos[2:2]) # Me quede solo con la columna de las observaciones 
Datos1 <- window(Datos1, end = time(Datos1)[length(Datos1) - 6])
Datos_prueba = window(Datos1, start  = time(Datos1)[length(Datos1) - 5])
length(Datos1)

## [1] 74

length(Datos_prueba)

## [1] 6

ts.plot(Datos1,col = c("red", "blue"), lwd = 1.5,
        main = " número de empresas registradas por trimestre menos los ultimos registros",
        ylab = "Valor", xlab = "Tiempo")

Es necesario que apliquemos logaritmo a los datos para estabilizarlos en Varianza y hacer una diferenciacion para eliminar tendencias.

Con el ACF se evidensia que existe mucha autocorrelacion en los datos y además podemos ver que existe extacionaridad en los periodos 3 y 4.

Propuesta de modelos

“Propondremos diversos modelos predictivos. Algunos serán generados mediante funciones o algoritmos del software, mientras que otros se basarán en la observación empírica de las funciones de autocorrelación (ACF) y autocorrelación parcial (PACF).

set.seed(123)
modelo_propuestos <- auto.arima(
  Wt,
  seasonal = TRUE,       # Habilitar componente estacional
  stepwise = FALSE,      # Búsqueda exhaustiva de parámetros
  approximation = FALSE, # Usar máxima precisión
  trace = TRUE           # Mostrar progreso
)

## 
##  ARIMA(0,0,0) with zero mean     : 97.94308
##  ARIMA(0,0,0) with non-zero mean : 99.68753
##  ARIMA(0,0,1) with zero mean     : 58.30182
##  ARIMA(0,0,1) with non-zero mean : 48.28555
##  ARIMA(0,0,2) with zero mean     : 43.46524
##  ARIMA(0,0,2) with non-zero mean : 40.0208
##  ARIMA(0,0,3) with zero mean     : 43.44556
##  ARIMA(0,0,3) with non-zero mean : 41.71955
##  ARIMA(0,0,4) with zero mean     : 39.39784
##  ARIMA(0,0,4) with non-zero mean : 37.64341
##  ARIMA(0,0,5) with zero mean     : 19.80356
##  ARIMA(0,0,5) with non-zero mean : 13.38677
##  ARIMA(1,0,0) with zero mean     : 62.94624
##  ARIMA(1,0,0) with non-zero mean : 63.40917
##  ARIMA(1,0,1) with zero mean     : 53.93497
##  ARIMA(1,0,1) with non-zero mean : 44.33628
##  ARIMA(1,0,2) with zero mean     : 43.82816
##  ARIMA(1,0,2) with non-zero mean : 41.85887
##  ARIMA(1,0,3) with zero mean     : Inf
##  ARIMA(1,0,3) with non-zero mean : Inf
##  ARIMA(1,0,4) with zero mean     : 15.41788
##  ARIMA(1,0,4) with non-zero mean : 13.27869
##  ARIMA(2,0,0) with zero mean     : 63.29713
##  ARIMA(2,0,0) with non-zero mean : 63.13152
##  ARIMA(2,0,1) with zero mean     : 56.06103
##  ARIMA(2,0,1) with non-zero mean : 46.37603
##  ARIMA(2,0,2) with zero mean     : 43.23138
##  ARIMA(2,0,2) with non-zero mean : 41.34535
##  ARIMA(2,0,3) with zero mean     : Inf
##  ARIMA(2,0,3) with non-zero mean : Inf
##  ARIMA(3,0,0) with zero mean     : 7.89026
##  ARIMA(3,0,0) with non-zero mean : -13.44335
##  ARIMA(3,0,1) with zero mean     : -1.242213
##  ARIMA(3,0,1) with non-zero mean : -14.89054
##  ARIMA(3,0,2) with zero mean     : Inf
##  ARIMA(3,0,2) with non-zero mean : -12.47063
##  ARIMA(4,0,0) with zero mean     : -6.013763
##  ARIMA(4,0,0) with non-zero mean : -15.12667
##  ARIMA(4,0,1) with zero mean     : Inf
##  ARIMA(4,0,1) with non-zero mean : -12.91923
##  ARIMA(5,0,0) with zero mean     : -6.463208
##  ARIMA(5,0,0) with non-zero mean : -12.84088
## 
## 
## 
##  Best model: ARIMA(4,0,0) with non-zero mean

print(summary(modelo_propuestos))

## Series: Wt 
## ARIMA(4,0,0) with non-zero mean 
## 
## Coefficients:
##           ar1      ar2      ar3     ar4    mean
##       -0.7272  -0.6379  -0.6048  0.2386  0.0365
## s.e.   0.1141   0.1222   0.1246  0.1163  0.0084
## 
## sigma^2 = 0.04021:  log likelihood = 14.2
## AIC=-16.4   AICc=-15.13   BIC=-2.66
## 
## Training set error measures:
##                       ME      RMSE       MAE      MPE     MAPE     MASE
## Training set -0.00518775 0.1935406 0.1515188 60.31877 104.1615 0.209071
##                      ACF1
## Training set -0.008806644

# el algoritmo nos recomienda trabajar con un ARIMA(4,0,0) 
Arima_modelo1 = arima(Wt_center, order = c(4, 0, 0))

# Proponemos algunos modelos basandonos en los resagos de las funciones ACF y PACF

Arima_modelo2 = arima(Wt_center, order = c(3, 0, 3))
Arima_modelo3 = arima(Wt_center, order = c(0, 0, 3))
Arima_modelo4 = arima(Wt_center, order = c(4, 0, 1))

# hasta este punto es necesario plantear un modelo SARIMA por el comportamiento de los datos en el ACF
Arima_modelo5= arima(Wt_center, order = c(4, 0, 1), seasonal = c(3,0,3))

#------------------------------------------------------------------------------------------------------
# Ahora vamos a proponer  los modelos TAR 

library(tsDyn)


mejores_resultados <- data.frame(m = integer(),
                                 thDelay = integer(),
                                 nthresh = integer(),
                                 AIC = numeric())

# Rango de búsqueda de parámetros
for (m in 1:5) {
  for (thDelay in 1:m) {
    for (nthresh in 1:2) {
      modelo_try <- try(setar(Wt, m = m, thDelay = thDelay, nthresh = nthresh), silent = TRUE)
      if (!inherits(modelo_try, "try-error")) {
        aic_val <- AIC(modelo_try)
        mejores_resultados <- rbind(mejores_resultados,
                                    data.frame(m = m, thDelay = thDelay, nthresh = nthresh, AIC = aic_val))
      }
    }
  }
}

# Mostrar el mejor modelo según AIC
mejor_modelo_info <- mejores_resultados[which.min(mejores_resultados$AIC), ]
print(mejor_modelo_info)

##    m thDelay nthresh       AIC
## 12 4       3       2 -266.6946

# Ajustar el mejor modelo
modelo_tar_optimo <- setar(Wt,
                           m = mejor_modelo_info$m,
                           thDelay = mejor_modelo_info$thDelay,
                           nthresh = mejor_modelo_info$nthresh)

# El modelo dos es la propuesta por el algoritmo y el modelo uno es de manera empirica
modelo_tar1 <- setar(Wt, m=2, thDelay=1, nthresh=1) 
modelo_tar2 <- setar(Wt, m=4, thDelay=3, nthresh=2)

verificar el AIC

Una ves tenemos ya los modelos propuestos vamos hacer la eleccion los medelos que mejor se ajusten a los datos basandonos en el criterio del AIC.

##    Modelo        AIC
## 1 ARIMA 1  -16.39940
## 2 ARIMA 2  -12.04054
## 3 ARIMA 3   40.82403
## 4 ARIMA 4  -14.64231
## 5 ARIMA 5  -24.19496
## 6   TAR 1 -169.66958
## 7   TAR 2 -266.69456

## 
##  Jarque Bera Test
## 
## data:  res1
## X-squared = 1.4634, df = 2, p-value = 0.4811

## 
##  Jarque Bera Test
## 
## data:  res2
## X-squared = 0.20817, df = 2, p-value = 0.9011

Los modelo 1, sugerido por la función auto.arima, y el modelo empírico 5 destacaron por presentar los mejores valores de AIC. De manera similar, en los modelos TAR, el algoritmo identificó un modelo con un AIC superior. En consecuencia, seleccionaremos estos modelos específicos para llevar a cabo la comparación de sus predicciones. Es así, que vamos a elegir los modelos mencionados para hacer la comparacion de las predicciones.

Predicciones de prueba a los modelos seleccionados

## [1] 3415 2499 3791 3368 3885 3009

## [1] 3806 2151 3484 3598 3772 2433

## [1] 3635 2763 1758 1766 2346 2600

##   Modelo Prediccion1 Prediccion2 Prediccion3 Prediccion4 Prediccion5
## 1  ARIMA        3415        2499        3791        3368        3885
## 2 SARIMA        3806        2151        3484        3598        3772
## 3    TAR        3635        2763        1758        1766        2346
##   Prediccion6
## 1        3009
## 2        2433
## 3        2600

Verificación de los modelos

Para evaluar la eficiencia predictiva de los modelos seleccionados, los compararemos con los datos de prueba reservados. Utilizaremos la función \(Accuracy\) para cuantificar los errores a través de las métricas \(ME,RMSE.MAE,MAPE,Theils_U\).

##   Modelo        ME     RMSE      MAE     MAPE  Theils_U
## 1  ARIMA -441.1667 697.7195 592.5000 22.65131 0.8136142
## 2 SARIMA -320.6667 857.1814 763.3333 29.18912 0.9096555
## 3    TAR  408.6667 776.4073 763.3333 29.18912 0.9096555

Conclusiones de los modelos

El modelo ARIMA exhibe errores menores en \(RMSE, MAE , MAPE\), lo que señala una mayor proximidad de sus predicciones a los valores reales. Adicionalmente, su valor de Theils_U de \(0.81\) indica un rendimiento superior a una predicción ingenua y una mayor precisión en comparación con los otros dos modelos. Por otro lado, el modelo \(TAR\) presenta un ME positivo y de menor magnitud, sugiriendo una tendencia a la sobreestimación, a diferencia de los modelos ARIMA y SARIMA que tienden a subestimar. No obstante, el modelo ARIMA continúa demostrando el menor error, incluso en comparación con el modelo SARIMA.

Es importante destacar que el error medio absoluto porcentual \((MAPE)\) del modelo ARIMA, aunque el más bajo de los tres, se sitúa en un \(22.65%\). Este porcentaje aún podría considerarse elevado si se busca una mayor precisión en las predicciones.

Prediccion con los datos completos

Datos2 <- tail(Datos[[2]], 80)  
length(Datos2)

## [1] 80

Xt1=log(Datos2)
Wt1 = diff(ts(Xt1))
mu1= mean(Wt1)
Modelo = arima(Wt1-mu1,order = c(4,0,0))
summary(Modelo)

## 
## Call:
## arima(x = Wt1 - mu1, order = c(4, 0, 0))
## 
## Coefficients:
##           ar1      ar2      ar3     ar4  intercept
##       -0.6904  -0.6048  -0.5716  0.2765     0.0066
## s.e.   0.1083   0.1156   0.1178  0.1104     0.0084
## 
## sigma^2 estimated as 0.03623:  log likelihood = 16.82,  aic = -21.65
## 
## Training set error measures:
##                        ME     RMSE       MAE       MPE     MAPE      MASE
## Training set -0.004720159 0.190332 0.1481422 -721.6201 809.5509 0.2116445
##                     ACF1
## Training set -0.01249152

# Verificacion dle comportamiento de los residuales
tsdiag(Modelo)

# Prueba de normalidad de los residuales
library(tseries)
jarque.bera.test(Modelo$residuals)

## 
##  Jarque Bera Test
## 
## data:  Modelo$residuals
## X-squared = 1.3154, df = 2, p-value = 0.518

# el valos del P-Value >0.05 no rechaza la hipotesis nula de normalidad luego los residuales son normales
pred_Modelo <- predict(Modelo, n.ahead = 4)$pred

Yt <- exp(pred_Modelo + mu1)  # Usa +mu si centraste antes

# Vector para guardar los valores reconstruidos (6 pasos de predicción)
tforecas <- numeric(4)
tforecas[1] <- Datos2[80]* Yt[1]  # 2959 es el último dato real antes del pronóstico

# Reconstruir los siguientes valores
for (i in 2:length(Yt)) {
  tforecas[i] <- tforecas[i - 1] * Yt[i]
}

# Redondear los resultados
tforecas <- round(tforecas)
tforecas

## [1] 3690 3379 3610 2753

#------------- 

# GRAFICA 
library(ggplot2)
Datos2 <- as.numeric(Datos2)
tforecas <- as.numeric(tforecas)

serie_completa <- c(Datos2, tforecas)


tiempo <- seq_along(serie_completa)

n_reales <- length(Datos2)
par(mfrow = c(1, 1))
df_plot <- data.frame(
  Tiempo = 1:length(serie_completa),
  Valor = serie_completa,
  Tipo = c(rep("Real", n_reales), rep("Predicción", length(tforecas)))
)

ggplot(df_plot, aes(x = Tiempo, y = Valor, color = Tipo)) +
  geom_line(size = 0.6) +
  geom_point(size = 2) +
  scale_color_manual(values = c("Real" = "blue", "Predicción" = "red")) +
  geom_vline(xintercept = n_reales, linetype = "dashed", color = "green", size = 0.5) +
  labs(title = "Serie con predicción a 4 pasos", x = "Tiempo", y = "Número de empresas registradas") +
  theme_minimal()

Predicciones de la Camara de Comercio en la ciudad de Ibagué-Tolima

Camilo Yaraslov Pita C - Laura Sofia Viera

Datos

Propuesta de modelos

Propuesta de modelos

verificar el AIC

Predicciones de prueba a los modelos seleccionados

Verificación de los modelos

Conclusiones de los modelos

Prediccion con los datos completos