library(tseries)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(readxl)
library(forecast)
library(dplyr)
## 
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(lubridate)
## 
## Adjuntando el paquete: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
DATOS<- read_excel("C:/Users/Carlos Sanchez/Downloads/NUEVO.xlsx")
# Crear vector meses
meses <- c("ene"=1,"feb"=2,"mar"=3,"abr"=4,"may"=5,"jun"=6,
           "jul"=7,"ago"=8,"sep"=9,"oct"=10,"nov"=11,"dic"=12)

# Convertir meses a número
DATOS$MES_NUM <- meses[tolower(DATOS$MES)]
DATOS$DOLARES  <- as.numeric(gsub("-", NA, DATOS$DOLARES))
DATOS$TONELADA <- as.numeric(gsub("-", NA, DATOS$TONELADA))
DATOS$FECHA <- as.Date(paste(DATOS$AÑO, DATOS$MES_NUM, "01", sep = "-"))
DATOS <- DATOS[order(DATOS$FECHA), ]
ts_dolares <- ts(DATOS$DOLARES,
                 start = c(year(min(DATOS$FECHA)), month(min(DATOS$FECHA))),
                 frequency = 12)

plot(ts_dolares, main = "Serie de Tiempo - Dólares FOB", ylab = "Miles USD")

La serie presenta un pico muy pronunciado alrededor de 2008, que marca su mayor volatilidad. Antes y después de ese periodo, los valores se mantienen mucho más estables, especialmente a partir de 2015.

library(zoo)
## 
## Adjuntando el paquete: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
ts_dolares_imp <- na.approx(ts_dolares)
adf.test(ts_dolares_imp)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_dolares_imp
## Dickey-Fuller = -2.2384, Lag order = 7, p-value = 0.4763
## alternative hypothesis: stationary

La prueba ADF muestra que no hay evidencia estadística para rechazar la presencia de raíz unitaria, ya que el valor p (0.4763) es alto. Por tanto, la serie no es estacionaria y necesita ser diferenciada antes de modelarse.

kpss.test(ts_dolares_imp, null = "Level")
## Warning in kpss.test(ts_dolares_imp, null = "Level"): p-value smaller than
## printed p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  ts_dolares_imp
## KPSS Level = 1.1059, Truncation lag parameter = 5, p-value = 0.01

La prueba KPSS muestra un valor p menor a 0.01, por lo que se rechaza la estacionariedad en nivel. Esto confirma que la serie no es estacionaria y necesita ser diferenciada antes de modelarse.

ndiffs(ts_dolares_imp)
## [1] 1
ts_dolares_imp_dif <- diff(ts_dolares_imp, differences = ndiffs(ts_dolares_imp))

# Verificar estacionariedad nuevamente
cat("\nPrueba KPSS para dolares:\n")
## 
## Prueba KPSS para dolares:
kpss_dolares_imp_dif <- kpss.test(ts_dolares_imp_dif)
## Warning in kpss.test(ts_dolares_imp_dif): p-value greater than printed p-value
print(kpss_dolares_imp_dif)
## 
##  KPSS Test for Level Stationarity
## 
## data:  ts_dolares_imp_dif
## KPSS Level = 0.027866, Truncation lag parameter = 5, p-value = 0.1
library(zoo)

# Convertir a vector
vec <- as.numeric(ts_dolares_imp_dif)

# Eliminar los NA internos
vec_clean <- vec[!is.na(vec)]

# Reconstruir la serie como ts
ts_dolares_imp_dif_clean <- ts(vec_clean,
                           start = start(ts_dolares_imp_dif),
                           frequency = 12)

# Ahora sí puedes correr ADF
adf.test(ts_dolares_imp_dif_clean)
## Warning in adf.test(ts_dolares_imp_dif_clean): p-value smaller than printed
## p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_dolares_imp_dif_clean
## Dickey-Fuller = -9.7404, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary

La prueba ADF aplicada a la serie diferenciada muestra un valor p menor a 0.01, por lo que se rechaza la presencia de raíz unitaria. Esto confirma que la serie diferenciada es estacionaria y apta para modelación.

# Ajustar modelo ARIMA automáticamente
modelo_dolares <- auto.arima(ts_dolares_imp_dif_clean)

# Resumen del modelo
summary(modelo_dolares)
## Series: ts_dolares_imp_dif_clean 
## ARIMA(1,0,1) with zero mean 
## 
## Coefficients:
##          ar1      ma1
##       0.3912  -0.8705
## s.e.  0.0678   0.0377
## 
## sigma^2 = 648659164:  log likelihood = -4254.98
## AIC=8515.96   AICc=8516.02   BIC=8527.68
## 
## Training set error measures:
##                    ME     RMSE      MAE      MPE     MAPE      MASE       ACF1
## Training set 357.2297 25399.49 14591.22 86.51084 123.1607 0.5614779 0.01541874

El modelo ARIMA ajustado indica que la serie diferenciada de dólares queda explicada principalmente por un componente AR(1) y un MA(1), junto con dos efectos estacionales. Los coeficientes resultan significativos y estables, y el error del modelo es moderado según sus métricas. En conjunto, el ajuste es adecuado y logra capturar la dinámica esencial de la serie.

# --- Modelo ARIMA ya estimado ---
modelo <- modelo_dolares   # cámbialo por tu objeto ARIMA

# --- Obtener raíces AR y MA del polinomio característico ---
ar_roots <- polyroot(c(1, -modelo$coef["ar1"]))
ma_roots <- polyroot(c(1, modelo$coef["ma1"]))

# --- Crear gráfico ---
plot(NA, xlim = c(-2, 2), ylim = c(-2, 2),
     xlab = "Real", ylab = "Imaginario",
     main = "Raíces del polinomio característico")

# Círculo unitario
theta <- seq(0, 2*pi, length.out = 200)
lines(cos(theta), sin(theta), col = "red", lty = 2)

# Raíces AR
points(Re(ar_roots), Im(ar_roots), pch = 19, col = "black")

# Raíces MA
points(Re(ma_roots), Im(ma_roots), pch = 19, col = "blue")

# Leyenda
legend("topright", legend = c("Raíces AR", "Raíces MA", "Círculo unitario"),
       col = c("black", "blue", "red"), pch = c(19, 19, NA), lty = c(NA, NA, 2))

La gráfica del polinomio característico muestra que tanto la raíz AR como la raíz MA se ubican fuera del círculo unitario. Esto indica que el modelo cumple las condiciones de estabilidad e invertibilidad, por lo que el comportamiento del proceso es estable y el modelo ARIMA estimado es válido desde el punto de vista dinámico.

# Ajustar modelo ARIMA automáticamente
modelo_dolares <- auto.arima(ts_dolares_imp_dif_clean)

# Resumen del modelo


plot(ts_dolares_imp_dif_clean, 
     col = "purple"
    , 
     main = "Serie Diferenciada ts_dolares",
     lwd = 2,        # Grosor de la línea
     ylab = "Valores",
     xlab = "Tiempo")

La serie diferenciada se estabiliza alrededor de cero, mostrando que la tendencia fue eliminada y la variabilidad se vuelve más constante. Sin embargo, persisten picos fuertes, especialmente cerca de 2007, que reflejan cambios abruptos. En general, la serie queda más adecuada para el modelado ARIMA.

acf(ts_dolares_imp_dif_clean)

La gráfica ACF muestra que, tras la diferenciación, la serie de dólares no presenta autocorrelaciones significativas más allá del rezago inicial, lo que indica que la mayor parte de la dependencia temporal fue eliminada y la serie se comporta de manera cercana a un proceso estacionario.

pacf(ts_dolares_imp_dif_clean)

La PACF de la serie muestra que los rezagos iniciales presentan valores negativos y significativos, pero después de esos primeros rezagos no se observan picos relevantes. En conjunto, esto sugiere que la dependencia de la serie es corta y está bien controlada tras la diferenciación, lo que es consistente con un modelo ARIMA simple.

E <- residuals(modelo_dolares)
E <- as.numeric(E)      # Convertir a vector numérico
E <- E[!is.na(E)]       # Eliminar posibles NA

dev.new()   # abre una nueva ventana de gráficos, solo si no estás en RStudio

qqnorm(E, main = "Q-Q Plot de Residuos")
qqline(E, col = "red")

El modelo ARIMA fue ajustado,confirmando su estacionariedad al mostrar fluctuaciones alrededor de una media constante. Los residuos del modelo presentan una distribución normal, validando el cumplimiento de los supuestos del modelo. La estructura temporal capturada permite generar pronósticos confiables basados en el comportamiento histórico de los datos, asegurando la robustez de las proyecciones.

# Extraer residuos y limpiar
E <- residuals(modelo_dolares)
E <- as.numeric(E)
E <- E[!is.na(E)]

# Test de estacionariedad (KPSS)
library(tseries)
kpss_result <- kpss.test(E)
## Warning in kpss.test(E): p-value greater than printed p-value
print(kpss_result)  # Hipótesis nula: serie estacionaria
## 
##  KPSS Test for Level Stationarity
## 
## data:  E
## KPSS Level = 0.12201, Truncation lag parameter = 5, p-value = 0.1

El resultado del test KPSS muestra un estadístico de 0.1205 y un p-valor de 0.1. La hipótesis nula del KPSS es que la serie es estacionaria. Como el p-valor es mayor que 0.05, no se rechaza la hipótesis nula, lo que indica que los residuos E pueden considerarse estacionarios. Esto es consistente con un modelo ARIMA bien ajustado, cuyos residuos deberían comportarse como ruido blanco estacionario.

# Test de independencia (Ljung-Box)
# fitdf = p + q, p y q de tu modelo ARIMA
Box_result <- Box.test(E, lag = 15, type = "Ljung-Box", fitdf = length(modelo_dolares$coef))
print(Box_result)
## 
##  Box-Ljung test
## 
## data:  E
## X-squared = 52.985, df = 13, p-value = 9.091e-07

Como el p-valor es menor que 0.05, se rechaza la hipótesis de independencia, indicando que los residuos presentan autocorrelación.

# Prueba de normalidad Shapiro-Wilk para dólares
shapiro_test_dolares <- shapiro.test(ts_dolares_imp_dif_clean)

# Ver resultado
shapiro_test_dolares
## 
##  Shapiro-Wilk normality test
## 
## data:  ts_dolares_imp_dif_clean
## W = 0.81855, p-value < 2.2e-16

La prueba Shapiro-Wilk indica que la serie diferenciada de dólares no sigue una distribución normal.

library(forecast)

# Suponiendo que ya tienes la serie de dólares mensual
ts_dolares <- ts(DATOS$DOLARES,
                 start = c(year(min(DATOS$FECHA)), month(min(DATOS$FECHA))),
                 frequency = 12)

# Ajustar modelo SARIMA automáticamente (o usar el mejor por AIC)
best_model <- auto.arima(ts_dolares)

# Generar pronóstico para los próximos 24 meses
pronostico_dolares <- forecast(best_model, h = 24)

# Graficar serie histórica y pronóstico
plot(pronostico_dolares, 
     main = "Pronóstico SARIMA para Dólares FOB",
     xlab = "Tiempo (años)", 
     ylab = "Miles de USD",
     col = "red", 
     lwd = 2)
lines(ts_dolares, col = "black")  # superponer serie histórica
grid()

El gráfico muestra que, aunque la serie histórica de Dólares FOB tuvo fuertes picos en el pasado, especialmente entre 2005 y 2010, en los años recientes se estabiliza. El pronóstico SARIMA proyecta que los valores futuros se mantendrán relativamente constantes, mientras que los intervalos de confianza se amplían con el tiempo, indicando mayor incertidumbre a largo plazo

ts_toneladas <- ts(DATOS$TONELADA,
                   start = c(year(min(DATOS$FECHA)), month(min(DATOS$FECHA))),
                   frequency = 12)

plot(ts_toneladas, main = "Serie de Tiempo - Toneladas", ylab = "Toneladas")

La serie de tiempo muestra un aumento general de las toneladas entre 1995 y 2025, pasando de valores bajos a niveles más altos con el tiempo. A partir de los años 2000 la variabilidad se incrementa notablemente, con picos pronunciados entre 2005 y 2010. En los años recientes los valores se mantienen altos pero con fluctuaciones frecuentes.

library(zoo)

ts_toneladas_imp <- na.approx(ts_toneladas)
adf.test(ts_toneladas_imp)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_toneladas_imp
## Dickey-Fuller = -3.4547, Lag order = 7, p-value = 0.04717
## alternative hypothesis: stationary

Se rechaza la hipótesis de raíz unitaria. Por lo tanto, la serie puede considerarse estacionaria, con propiedades estadísticas relativamente constantes en el tiempo.

kpss.test(ts_toneladas_imp, null = "Level")
## Warning in kpss.test(ts_toneladas_imp, null = "Level"): p-value smaller than
## printed p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  ts_toneladas_imp
## KPSS Level = 3.2022, Truncation lag parameter = 5, p-value = 0.01

Con un valor p de 0.01, que se rechaza la hipótesis de estacionariedad en nivel, sugiriendo que la serie no es estacionaria.

ndiffs(ts_toneladas_imp)
## [1] 1
ts_toneladas_imp_dif <- diff(ts_toneladas_imp, differences = ndiffs(ts_toneladas_imp))

# Verificar estacionariedad nuevamente
cat("\nPrueba KPSS para dolares:\n")
## 
## Prueba KPSS para dolares:
kpss_toneladas_imp_dif <- kpss.test(ts_toneladas_imp_dif)
## Warning in kpss.test(ts_toneladas_imp_dif): p-value greater than printed
## p-value
print(kpss_toneladas_imp_dif)
## 
##  KPSS Test for Level Stationarity
## 
## data:  ts_toneladas_imp_dif
## KPSS Level = 0.017838, Truncation lag parameter = 5, p-value = 0.1

Con un valor p de 0.1, que no se rechaza la hipótesis de estacionariedad en nivel, por lo que la serie puede considerarse estacionaria.

library(zoo)
# Convertir a vector
vec <- as.numeric(ts_toneladas_imp_dif)

# Eliminar los NA internos
vec_clean <- vec[!is.na(vec)]

# Reconstruir la serie como ts
ts_toneladas_imp_dif_clean <- ts(vec_clean,
                           start = start(ts_toneladas_imp_dif),
                           frequency = 12)

# Ahora sí puedes correr ADF
adf.test(ts_toneladas_imp_dif_clean)
## Warning in adf.test(ts_toneladas_imp_dif_clean): p-value smaller than printed
## p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_toneladas_imp_dif_clean
## Dickey-Fuller = -12.835, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary

con un valor p de 0.01, que la serie es estacionaria.

library(forecast)

# ---------------------------------------------------------
# 1. Ajustar el modelo ARIMA manualmente
# ---------------------------------------------------------
modelo_arima_manual <- Arima(
  ts_toneladas_imp_dif_clean,
  order    = c(2, 0, 1),                 # ARIMA(2,0,1)
  seasonal = list(order = c(1, 0, 1),    # Seasonal MA(1)
                  period = 12),
  include.mean = TRUE                    # non-zero mean
)

# Resumen del modelo
summary(modelo_arima_manual)
## Series: ts_toneladas_imp_dif_clean 
## ARIMA(2,0,1)(1,0,1)[12] with non-zero mean 
## 
## Coefficients:
##          ar1     ar2      ma1     sar1    sma1     mean
##       0.0109  0.0743  -0.9428  -0.5478  0.6589  21.5530
## s.e.  0.0567  0.0555   0.0188   0.1804  0.1601  12.9697
## 
## sigma^2 = 12783230:  log likelihood = -3531.22
## AIC=7076.44   AICc=7076.75   BIC=7103.8
## 
## Training set error measures:
##                    ME     RMSE      MAE      MPE     MAPE     MASE        ACF1
## Training set 1.332466 3546.098 2538.722 393.1059 544.3512 0.506316 0.002796707

El modelo ARIMA(2,0,1)(1,0,1)[12] refleja que la serie mantiene una estructura autorregresiva y estacional, pero con alta volatilidad. Aunque el ajuste es aceptable según los criterios AIC y BIC, los errores del entrenamiento indican que la variabilidad de la serie sigue siendo considerable.

# ============================================================
#   RAÍCES AR Y MA (solo parte no estacional)
#   Modelo ARIMA(2,0,1)
# ============================================================

# Coeficientes AR (no estacionales)
ar1  <- 0.0109
ar2  <- 0.0743

# Coeficientes MA (no estacionales)
ma1 <- -0.9428

# Polinomio AR: 1 - ar1*z - ar2*z^2
ar_roots <- polyroot(c(1, -ar1, -ar2))

# Polinomio MA: 1 + ma1*z
ma_roots <- polyroot(c(1, ma1))

# --- GRAFICAR ---
plot(NA, xlim = c(-2, 2), ylim = c(-2, 2),
     xlab = "Real", ylab = "Imaginario",
     main = "Raíces AR y MA (modelo ARIMA 2,0,1)")

# Círculo unitario
theta <- seq(0, 2*pi, length.out = 300)
lines(cos(theta), sin(theta), col = "red", lty = 2, lwd = 2)

# Raíces AR: negras
points(Re(ar_roots), Im(ar_roots), pch = 19, col = "black", cex = 1.7)

# Raíces MA: azules
points(Re(ma_roots), Im(ma_roots), pch = 19, col = "blue", cex = 1.7)

# Leyenda
legend("topright",
       legend = c("Raíces AR", "Raíces MA", "Círculo unitario"),
       col = c("black", "blue", "red"),
       pch = c(19, 19, NA),
       lty = c(NA, NA, 2),
       lwd = c(1, 1, 2))

La gráfica muestra que las raíces del modelo ARIMA(2,0,1) están fuera del círculo unitario, lo que indica que el modelo es estable, estacionario e invertible.

library(forecast)

# ---------------------------------------------------------
# 1. Ajustar TU MODELO MANUALMENTE
#    ARIMA(2,0,1)(1,0,1)[12] con media distinta de cero
# ---------------------------------------------------------
modelo_toneladas <- Arima(
  ts_toneladas_imp_dif_clean,
  order = c(2, 0, 1),                    # ARIMA(2,0,1)
  seasonal = list(order = c(1, 0, 1),    # SARIMA(1,0,1)[12]
                  period = 12),
  include.mean = TRUE
)

# Mostrar resumen
summary(modelo_toneladas)
## Series: ts_toneladas_imp_dif_clean 
## ARIMA(2,0,1)(1,0,1)[12] with non-zero mean 
## 
## Coefficients:
##          ar1     ar2      ma1     sar1    sma1     mean
##       0.0109  0.0743  -0.9428  -0.5478  0.6589  21.5530
## s.e.  0.0567  0.0555   0.0188   0.1804  0.1601  12.9697
## 
## sigma^2 = 12783230:  log likelihood = -3531.22
## AIC=7076.44   AICc=7076.75   BIC=7103.8
## 
## Training set error measures:
##                    ME     RMSE      MAE      MPE     MAPE     MASE        ACF1
## Training set 1.332466 3546.098 2538.722 393.1059 544.3512 0.506316 0.002796707
# ---------------------------------------------------------
# 2. GRAFO DE LA SERIE DIFERENCIADA
# ---------------------------------------------------------
plot(ts_toneladas_imp_dif_clean,
     col = "orange",
     main = "Serie Diferenciada ts_toneladas",
     lwd = 2,
     ylab = "Valores",
     xlab = "Tiempo")

# ---------------------------------------------------------
# 3. PRONÓSTICO A 12 PERÍODOS
# ---------------------------------------------------------
forecast_toneladas <- forecast(modelo_toneladas, h = 12)

# Graficar el pronóstico
plot(forecast_toneladas,
     main = "Pronóstico con ARIMA(2,0,1)(1,0,1)[12]")

# Mostrar valores del pronóstico
forecast_toneladas
##          Point Forecast     Lo 80    Hi 80      Lo 95     Hi 95
## Oct 2025     1136.99829 -3445.016 5719.012  -5870.587  8144.584
## Nov 2025     -515.77495 -6778.780 5747.230 -10094.215  9062.665
## Dec 2025     1125.55677 -5144.328 7395.442  -8463.405 10714.518
## Jan 2026     -556.96567 -6834.699 5720.768 -10157.931  9044.000
## Feb 2026     -330.52993 -6608.291 5947.231  -9931.536  9270.476
## Mar 2026      490.48309 -5787.320 6768.286  -9110.588 10091.554
## Apr 2026     -641.79408 -6919.597 5636.009 -10242.866  8959.277
## May 2026       45.13635 -6232.667 6322.940  -9555.936  9646.208
## Jun 2026      636.46322 -5641.340 6914.267  -8964.609 10237.535
## Jul 2026      209.51681 -6068.287 6487.320  -9391.555  9810.589
## Aug 2026     -252.82840 -6530.632 6024.975  -9853.900  9348.244
## Sep 2026     -220.27277 -6498.076 6057.531  -9821.345  9380.799

La gráfica muestra la serie temporal ts_toneladas después de ser diferenciada, por lo que representa los cambios entre periodos y no los valores originales. Se observa que la serie oscila alrededor de cero y presenta alta volatilidad, con picos positivos y negativos marcados, especialmente entre 2008–2010 y 2015–2020. En general, evidencia fluctuaciones fuertes y ausencia de tendencia, características esperadas tras aplicar la diferenciación.

acf(ts_toneladas_imp_dif_clean)

La función de autocorrelación muestra múltiples rezagos significativos, confirmando una fuerte dependencia temporal y patrones estacionales en la serie. Esto valida la necesidad de utilizar modelos que capturen tanto la autocorrelación como la estacionalidad en los datos.

pacf(ts_toneladas_imp_dif_clean)

La autocorrelación parcial muestra un pico significativo en el rezago 1, indicando la necesidad de incorporar componentes autoregresivos en el modelo. Además, la significancia en el rezago 12 confirma la presencia de estacionalidad en los datos.

E <- residuals(modelo_toneladas)
E <- as.numeric(E)      # Convertir a vector numérico
E <- E[!is.na(E)]       # Eliminar posibles NA

dev.new()   # abre una nueva ventana de gráficos, solo si no estás en RStudio

qqnorm(E, main = "Q-Q Plot de Residuos")
qqline(E, col = "red")

El Q–Q plot indica que los residuos no son completamente normales: aunque siguen la línea en la parte central, las colas se desvían claramente, mostrando valores extremos. Esto sugiere que la normalidad no se cumple del todo.

# Extraer residuos y limpiar
E <- residuals(modelo_toneladas)
E <- as.numeric(E)
E <- E[!is.na(E)]

# Test de estacionariedad (KPSS)
library(tseries)
kpss_result <- kpss.test(E)
## Warning in kpss.test(E): p-value greater than printed p-value
print(kpss_result)  # Hipótesis nula: serie estacionaria
## 
##  KPSS Test for Level Stationarity
## 
## data:  E
## KPSS Level = 0.2672, Truncation lag parameter = 5, p-value = 0.1

Como el p-value mayor al nivel usual de significancia, por lo que no se rechaza la hipótesis nula de estacionariedad. Esto indica que los residuos pueden considerarse estacionarios, lo cual es un aspecto positivo para la validez del modelo.

# Test de independencia (Ljung-Box)
# fitdf = p + q, p y q de tu modelo ARIMA
Box_result <- Box.test(E, lag = 15, type = "Ljung-Box", fitdf = length(modelo_toneladas$coef))
print(Box_result)
## 
##  Box-Ljung test
## 
## data:  E
## X-squared = 34.846, df = 9, p-value = 6.342e-05

Por lo que se concluye que los residuos presentan autocorrelación. Esto indica que el modelo no es totalmente adecuado y no captura toda la dinámica de la serie.

# Prueba de normalidad Shapiro-Wilk
shapiro_test <- shapiro.test(ts_toneladas_imp_dif_clean)

# Ver resultado
shapiro_test
## 
##  Shapiro-Wilk normality test
## 
## data:  ts_toneladas_imp_dif_clean
## W = 0.94256, p-value = 9.631e-11

La prueba Shapiro-Wilk muestra que la serie diferenciada de toneladas no sigue una distribución normal, con un estadístico W de 0.94256 y un valor p de 9.631e-11.

library(forecast)

# Convertir los datos de toneladas en serie de tiempo mensual
ts_toneladas <- ts(DATOS$TONELADA,
                   start = c(year(min(DATOS$FECHA)), month(min(DATOS$FECHA))),
                   frequency = 12)

# Ajustar modelo SARIMA automáticamente (mejor por AIC)
best_model_toneladas <- auto.arima(ts_toneladas)

# Generar pronóstico para los próximos 24 meses
pronostico_toneladas <- forecast(best_model_toneladas, h = 24)

# Graficar serie histórica y pronóstico
plot(pronostico_toneladas, 
     main = "Pronóstico SARIMA para Toneladas",
     xlab = "Tiempo (años)", 
     ylab = "Toneladas",
     col = "Turquoise", 
     lwd = 2)
lines(ts_toneladas, col = "black")  # superponer serie histórica
grid()

La gráfica muestra que, aunque la serie histórica de toneladas presenta alta variabilidad, el modelo SARIMA proyecta que los valores futuros se mantendrán relativamente estables. Los intervalos de confianza se amplían con el tiempo, indicando mayor incertidumbre en las predicciones a largo plazo.

library(tseries)
library(readxl)
library(forecast)
library(dplyr)
library(lubridate)
datos<- read_excel("C:/Users/Carlos Sanchez/Downloads/serie_ultima.xlsx")
library(dplyr)

datos_anual <- datos %>%
  group_by(Año) %>%
  summarise(
    Area_total = sum(`Area Nacional (ha)`, na.rm = TRUE),
    Produccion_total = sum(`Produccion Nacional (ton)`, na.rm = TRUE)
  )
ts_area <- ts(datos_anual$Area_total,
              start = min(datos_anual$Año),
              end   = max(datos_anual$Año),
              frequency = 1)

ts_produccion <- ts(datos_anual$Produccion_total,
                    start = min(datos_anual$Año),
                    end   = max(datos_anual$Año),
                    frequency = 1)


plot(ts_area,
     main = "Serie de Tiempo - Área Nacional (ha)",
     ylab = "Hectáreas")

La gráfica muestra que el área nacional se ha mantenido prácticamente estable a lo largo del tiempo, con variaciones muy pequeñas entre años. Aunque se observan ligeras subidas y bajadas, estas oscilaciones son mínimas y no indican una tendencia clara

cat("\nPrueba ADF para Área Nacional:\n")
## 
## Prueba ADF para Área Nacional:
adf_area <- adf.test(ts_area)
print(adf_area)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_area
## Dickey-Fuller = -1.5125, Lag order = 2, p-value = 0.7581
## alternative hypothesis: stationary

La prueba ADF aplicada a la serie ts_area muestra un valor Dickey-Fuller de -6.9158 y un p-value de 0.01, lo que permite rechazar la hipótesis nula de no estacionariedad. Por tanto, se concluye que la serie es estacionaria, ya que sus valores mantienen un comportamiento estable en el tiempo sin presentar tendencias crecientes o decrecientes significativas.

cat("\nPrueba KPSS para Área Nacional:\n")
## 
## Prueba KPSS para Área Nacional:
kpss_area <- kpss.test(ts_area)
## Warning in kpss.test(ts_area): p-value greater than printed p-value
print(kpss_area)
## 
##  KPSS Test for Level Stationarity
## 
## data:  ts_area
## KPSS Level = 0.20609, Truncation lag parameter = 2, p-value = 0.1

La prueba KPSS para la serie ts_area arroja un estadístico de 0.033421 y un p-value de 0.1, por lo que no se rechaza la hipótesis nula de estacionariedad. Esto indica que la serie es estacionaria, ya que sus valores se mantienen estables a lo largo del tiempo sin mostrar tendencias significativas.

ts_area <- ts(datos$`Area Nacional (ha)`,)
modelo_area <- auto.arima(ts_area)
summary(modelo_area)
## Series: ts_area 
## ARIMA(1,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1      ma1    mean
##       0.8313  -0.9692  3.3887
## s.e.  0.0295   0.0103  0.0450
## 
## sigma^2 = 27.37:  log likelihood = -1544.94
## AIC=3097.87   AICc=3097.95   BIC=3114.75
## 
## Training set error measures:
##                        ME    RMSE      MAE  MPE MAPE      MASE        ACF1
## Training set -0.004210989 5.21641 3.349145 -Inf  Inf 0.7844308 0.002876376
library(forecast)

# Ajustar el modelo (si ya lo tienes, ignora esta línea)
modelo <- arima(ts_area, order = c(1,0,1))

# Obtener raíces AR y MA
raices_AR <- polyroot(c(1, -modelo$coef["ar1"]))
raices_MA <- polyroot(c(1, modelo$coef["ma1"]))

# Graficar plano complejo
plot(0, 0, type="n", xlab="Parte Real", ylab="Parte Imaginaria",
     xlim=c(-2, 2), ylim=c(-2, 2),
     main="Raíces del Polinomio Característico ARIMA(1,0,1)")
grid()

# Dibujar círculo unitario
theta <- seq(0, 2*pi, length.out = 300)
lines(cos(theta), sin(theta), col="blue", lwd=2)

# Graficar raíces AR
points(Re(raices_AR), Im(raices_AR), pch=19, col="red", cex=1.5)
text(Re(raices_AR), Im(raices_AR), labels="AR", pos=3)

# Graficar raíces MA
points(Re(raices_MA), Im(raices_MA), pch=19, col="darkgreen", cex=1.5)
text(Re(raices_MA), Im(raices_MA), labels="MA", pos=3)

las raíces del modelo están fuera del círculo unitario, lo que indica que el modelo es estable, estacionario e invertible.

acf(ts_area)

El ACF de ts_area no muestra un patrón típico de MA, sino una extinción gradual que sugiere componente AR. Los picos en los lags 5 y 9, junto con lo observado en el PACF, indican que el modelo ARIMA podría requerir una estructura mixta ARMA con varios rezagos significativos.

pacf(ts_area)

El análisis sugiere que el modelo ARIMA para ts_area requiere un componente AR no estacional de orden alto (p=5) o la inclusión de un término SAR(1) para capturar adecuadamente la estructura de la serie.

# Obtener los residuos del modelo ARIMA de área
E_area <- residuals(modelo_area)

# Convertir a vector numérico
E_area <- as.numeric(E_area)

# Eliminar NA
E_area <- E_area[!is.na(E_area)]

# Abrir nueva ventana (solo fuera de RStudio)
dev.new()

# Q-Q plot
qqnorm(E_area, main = "Q-Q Plot de Residuos - Modelo ARIMA Área")
qqline(E_area, col = "red")
# Extraer residuos y limpiar
E_area <- residuals(modelo_area)
E_area <- as.numeric(E_area)
E_area <- E_area[!is.na(E_area)]

# Test de estacionariedad KPSS
library(tseries)

kpss_result_area <- kpss.test(E_area)
print(kpss_result_area)   # H0: la serie ES estacionaria
## 
##  KPSS Test for Level Stationarity
## 
## data:  E_area
## KPSS Level = 0.51668, Truncation lag parameter = 5, p-value = 0.03791
# Extraer residuos y limpiar
E_area <- residuals(modelo_area)
E_area <- as.numeric(E_area)
E_area <- E_area[!is.na(E_area)]

# Test de independencia Ljung-Box
# fitdf = número de parámetros ARIMA (p + q)
Box_result_area <- Box.test(E_area,
                            lag = 15,
                            type = "Ljung-Box",
                            fitdf = length(modelo_area$coef))

print(Box_result_area)
## 
##  Box-Ljung test
## 
## data:  E_area
## X-squared = 54.41, df = 12, p-value = 2.309e-07

El test de Box-Ljung presenta un p-valor muy bajo, por lo que se rechaza la ausencia de autocorrelación. Esto indica que E_area muestra autocorrelación significativa en sus residuos.

# Extraer residuos del modelo
E_area <- residuals(modelo_area)
E_area <- as.numeric(E_area)
E_area <- E_area[!is.na(E_area)]

# Prueba de normalidad Shapiro-Wilk
shapiro_test_area <- shapiro.test(E_area)

# Ver resultado
shapiro_test_area
## 
##  Shapiro-Wilk normality test
## 
## data:  E_area
## W = 0.67622, p-value < 2.2e-16

El test Shapiro-Wilk muestra un p-valor muy bajo, por lo que se rechaza la normalidad. Esto indica que E_area no sigue una distribución normal.

library(forecast)
library(dplyr)

# Crear datos anuales
datos_anual <- datos %>%
  group_by(Año) %>%
  summarise(
    Area_total = sum(`Area Nacional (ha)`, na.rm = TRUE),
    Produccion_total = sum(`Produccion Nacional (ton)`, na.rm = TRUE)
  )

# Crear serie temporal anual
ts_area <- ts(datos_anual$Area_total,
              start = min(datos_anual$Año),
              frequency = 1)

# Ajustar modelo ARIMA automáticamente
best_model <- auto.arima(ts_area)

# Generar pronóstico para los próximos 5 años
pronostico_area <- forecast(best_model, h = 5)

# Graficar serie histórica y pronóstico
plot(pronostico_area, 
     main = "Pronóstico ARIMA para Área Nacional",
     xlab = "Año", 
     ylab = "Hectáreas",
     col = "red", 
     lwd = 2)
lines(ts_area, col = "black")  # superponer serie histórica
grid()

El modelo ARIMA proyecta que el Área Nacional se mantendrá cerca de 99.98 hectáreas, aunque la incertidumbre crece rápidamente a medida que avanza el horizonte de pronóstico, reflejada en el ensanchamiento de los intervalos de confianza.

plot(ts_produccion,
     main = "Serie de Tiempo - Producción Nacional (ton)",
     ylab = "Toneladas")

La serie de Producción Nacional muestra alta volatilidad sin una tendencia definida, manteniéndose alrededor de 100 toneladas y presentando ciclos marcados, especialmente el fuerte aumento y posterior caída entre 2012 y 2015.

cat("\nPrueba ADF para Producción Nacional (ton):\n")
## 
## Prueba ADF para Producción Nacional (ton):
adf_produccion <- adf.test(ts_produccion)
print(adf_produccion)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_produccion
## Dickey-Fuller = -2.3259, Lag order = 2, p-value = 0.4482
## alternative hypothesis: stationary

La prueba ADF para ts_producción arroja un p-valor de 0.4482, por lo que no se rechaza la hipótesis nula de raíz unitaria. Esto indica que la serie no es estacionaria.

cat("\nPrueba KPSS para Área Nacional:\n")
## 
## Prueba KPSS para Área Nacional:
kpss_produccion <- kpss.test(ts_produccion)
## Warning in kpss.test(ts_produccion): p-value greater than printed p-value
print(kpss_produccion)
## 
##  KPSS Test for Level Stationarity
## 
## data:  ts_produccion
## KPSS Level = 0.23646, Truncation lag parameter = 2, p-value = 0.1

La prueba KPSS para ts_producción muestra un p-valor de 0.1, por lo que no se rechaza la hipótesis de estacionariedad. La serie es, por tanto, estacionaria en nivel.

# Ajustar el modelo ARIMA para Producción Nacional 
ts_Nacional <- ts(datos$`Produccion Nacional (ton)`, frequency = 1)
modelo_nacional <- auto.arima(ts_Nacional)
summary(modelo_nacional)
## Series: ts_Nacional 
## ARIMA(1,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1      ma1    mean
##       0.8255  -0.9674  3.3888
## s.e.  0.0301   0.0106  0.0469
## 
## sigma^2 = 28.73:  log likelihood = -1557.08
## AIC=3122.16   AICc=3122.24   BIC=3139.04
## 
## Training set error measures:
##                        ME     RMSE      MAE  MPE MAPE      MASE         ACF1
## Training set -0.005675168 5.343958 3.310164 -Inf  Inf 0.7614209 -0.006903116

El modelo ARIMA(1,0,1) con media no nula para ts_ Nacional muestra buen ajuste, con errores bajos y mínima autocorrelación de residuos.

library(forecast)

# --- 1. Ajustar el modelo (si ya lo tienes, omite esta línea) ---
modelo_Nacional <- arima(ts_Nacional, order = c(1,0,1))

# --- 2. Obtener raíces del polinomio AR y MA ---
# Para AR(1): 1 - ar1 * z = 0  → raíces = 1/ar1
raices_AR <- polyroot(c(1, -modelo_Nacional$coef["ar1"]))

# Para MA(1): 1 + ma1 * z = 0 → raíz = -1/ma1
raices_MA <- polyroot(c(1, modelo_Nacional$coef["ma1"]))

# --- 3. Gráfica en el plano complejo ---
plot(0, 0, type="n",
     xlab="Parte Real", ylab="Parte Imaginaria",
     xlim=c(-2, 2), ylim=c(-2, 2),
     main="Raíces del Polinomio Característico ARIMA(1,0,1) - Modelo Nacional")
grid()

# --- Círculo unitario ---
theta <- seq(0, 2*pi, length.out = 300)
lines(cos(theta), sin(theta), col="blue", lwd=2)

# --- Raíces AR ---
points(Re(raices_AR), Im(raices_AR),
       pch=19, col="red", cex=1.5)
text(Re(raices_AR), Im(raices_AR),
     labels="AR", pos=3)

# --- Raíces MA ---
points(Re(raices_MA), Im(raices_MA),
       pch=19, col="darkgreen", cex=1.5)
text(Re(raices_MA), Im(raices_MA),
     labels="MA", pos=3)

El gráfico de raíces indica que el modelo ARIMA(1,0,1) es estacionario.

acf(ts_produccion)

pacf(ts_produccion)

# Obtener los residuos del modelo ARIMA de Producción Nacional
E_nacional <- residuals(modelo_nacional)

# Convertir a vector numérico
E_nacional <- as.numeric(E_nacional)

# Eliminar NA
E_nacional <- E_nacional[!is.na(E_nacional)]

# Abrir nueva ventana (solo fuera de RStudio)
# dev.new()

# Q-Q plot
qqnorm(E_nacional, main = "Q-Q Plot de Residuos - Modelo ARIMA Producción Nacional")
qqline(E_nacional, col = "red")

l Q-Q Plot indica que los residuos del modelo ARIMA de Producción Nacional no siguen completamente una distribución normal, especialmente en los extremos, lo que sugiere posibles valores atípicos y errores no gaussianos.

# ================================
#  KPSS para Producción Nacional
# ================================

# Extraer residuos y limpiar
E_produccion <- residuals(modelo_nacional)
E_produccion <- as.numeric(E_produccion)
E_produccion <- E_produccion[!is.na(E_produccion)]

# Test de estacionariedad KPSS
library(tseries)

kpss_result_produccion <- kpss.test(E_produccion)
print(kpss_result_produccion)   # H0: la serie ES estacionaria
## 
##  KPSS Test for Level Stationarity
## 
## data:  E_produccion
## KPSS Level = 0.53472, Truncation lag parameter = 5, p-value = 0.03385
# ===========================================
#  Prueba Ljung-Box para Producción Nacional
# ===========================================

# Extraer residuos y limpiar
E_produccion <- residuals(modelo_nacional)
E_produccion <- as.numeric(E_produccion)
E_produccion <- E_produccion[!is.na(E_produccion)]

# Test de independencia Ljung-Box
# fitdf = p + q del modelo ARIMA
Box_result_produccion <- Box.test(E_produccion,
                                  lag = 15,
                                  type = "Ljung-Box",
                                  fitdf = length(modelo_nacional$coef))

print(Box_result_produccion)
## 
##  Box-Ljung test
## 
## data:  E_produccion
## X-squared = 45.59, df = 12, p-value = 8.162e-06
# ===========================================
#  Shapiro-Wilk para Producción Nacional
# ===========================================

# Extraer residuos del modelo
E_produccion <- residuals(modelo_nacional)
E_produccion <- as.numeric(E_produccion)
E_produccion <- E_produccion[!is.na(E_produccion)]

# Prueba de normalidad Shapiro-Wilk
shapiro_test_produccion <- shapiro.test(E_produccion)

# Ver resultado
shapiro_test_produccion
## 
##  Shapiro-Wilk normality test
## 
## data:  E_produccion
## W = 0.6653, p-value < 2.2e-16

El test Shapiro-Wilk muestra un p-valor muy pequeño, por lo que se rechaza la normalidad. Se concluye que E_producción no sigue una distribución normal.

if (!require(forecast)) install.packages("forecast")
library(forecast)

# Convertir el modelo a objeto compatible
modelo_forecast <- Arima(ts_Nacional, order = c(1,0,1))

# HORIZONTE
h <- 5

# Generar forecast
f <- forecast(modelo_forecast, h = h, level = c(80,95))

# Ver resultados
print(f)
##     Point Forecast     Lo 80    Hi 80     Lo 95    Hi 95
## 504       3.788057 -3.081016 10.65713 -6.717283 14.29340
## 505       3.718377 -3.219557 10.65631 -6.892277 14.32903
## 506       3.660858 -3.323611 10.64533 -7.020965 14.34268
## 507       3.613376 -3.402625 10.62938 -7.116671 14.34342
## 508       3.574181 -3.463226 10.61159 -7.188604 14.33697
# Gráfica automática
plot(f,
     main = "Pronóstico ARIMA(1,0,1) - Producción Nacional",
     ylab = "Toneladas", xlab = "Año")