library(readxl)
library(dplyr)

## 
## Attaching package: '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)

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

library(ggplot2)
library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(sandwich)
library(gvlma)
library(broom)

Importación y preparación ———————————————

datos <- read.csv("RSAFS.csv") %>%
  mutate(
    fecha = as.Date(observation_date),
    y = RSAFS,
    t = 1:n(),
    mes = factor(month(fecha), levels = 1:12,
                 labels = c("Ene","Feb","Mar","Abr","May","Jun",
                            "Jul","Ago","Sep","Oct","Nov","Dic"))
  )
head(datos)

1)Introducción a los pronósticos (3.1) ———————————————

Tendencia: patrón de crecimiento o decrecimiento subyacente en una serie de tiempo Estacional: patrón que se repite de manera anual Cíclica: se refiere a patrones que se pueden identificar dentro de un intervalo amplio de tiempo (2 años, 5 años, 10 años) Irregular: Refleja los aspectos aleatorios en la medición de la serie de tiempo

# Serie mensual de ventas
p_serie <- ggplot(datos, aes(x = fecha, y = y)) +
  geom_line(linewidth = 0.7) +
  labs(
    title = "Serie mensual de ventas RSAFS",
    subtitle = "Enero 1992 - mayo 2025",
    x = "Fecha", y = "Millones de dólares"
  ) +
  theme_minimal(base_size = 11)
print(p_serie)

# 2) Modelado de la tendencia por funciones polinomiales (3.2) ——————————————— Ajuste de modelos de regresión con tiempo (t) lineal, cuadrático y cúbico

mod_lin  <- lm(y ~ t, data = datos)
mod_cuad <- lm(y ~ t + I(t^2), data = datos)
mod_cub  <- lm(y ~ t + I(t^2) + I(t^3), data = datos)

comparacion_poly <- tibble(
  Modelo = c("Lineal", "Cuadrático", "Cúbico"),
  SSE = c(sum(resid(mod_lin)^2), sum(resid(mod_cuad)^2), sum(resid(mod_cub)^2)),
  AIC = c(AIC(mod_lin), AIC(mod_cuad), AIC(mod_cub)),
  R2  = c(summary(mod_lin)$r.squared,
          summary(mod_cuad)$r.squared,
          summary(mod_cub)$r.squared)
)
print(comparacion_poly)

## # A tibble: 3 × 4
##   Modelo               SSE   AIC    R2
##   <chr>              <dbl> <dbl> <dbl>
## 1 Lineal     617487392282. 9627. 0.930
## 2 Cuadrático 310191589367. 9353. 0.965
## 3 Cúbico     130204146136. 9007. 0.985

Selección del modelo óptimo por menor AIC / SSE

modelo_tendencia <- mod_cub
print(modelo_tendencia)

## 
## Call:
## lm(formula = y ~ t + I(t^2) + I(t^3), data = datos)
## 
## Coefficients:
## (Intercept)            t       I(t^2)       I(t^3)  
##   1.331e+05    1.993e+03   -8.174e+00    1.739e-02

Gráfico comparativo de ajustes

ajustes_poly <- datos %>%
  mutate(
    lineal = fitted(mod_lin),
    cuadratico = fitted(mod_cuad),
    cubico = fitted(mod_cub)
  )

p_poly <- ggplot(ajustes_poly, aes(x = fecha)) +
  geom_line(aes(y = y, color = "Observado"), linewidth = 0.6) +
  geom_line(aes(y = lineal, color = "Lineal"), linewidth = 0.7) +
  geom_line(aes(y = cuadratico, color = "Cuadrático"), linewidth = 0.7) +
  geom_line(aes(y = cubico, color = "Cúbico"), linewidth = 0.8) +
  labs(title = "Comparación de tendencias polinomiales",
       x = "Fecha", y = "Millones de dólares", color = "Serie") +
  theme_minimal(base_size = 11)
print(p_poly)

# 3) Estacionalidad con variables ficticias (3.3) ————————–

mod_dummy <- lm(y ~ t + I(t^2) + I(t^3) + mes, data = datos)
summary(mod_dummy)

## 
## Call:
## lm(formula = y ~ t + I(t^2) + I(t^3) + mes, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -145840   -9773   -2023    9589   48426 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.333e+05  4.734e+03  28.168   <2e-16 ***
## t            1.993e+03  7.958e+01  25.047   <2e-16 ***
## I(t^2)      -8.180e+00  4.597e-01 -17.796   <2e-16 ***
## I(t^3)       1.740e-02  7.517e-04  23.146   <2e-16 ***
## mesFeb      -8.673e+02  4.447e+03  -0.195    0.845    
## mesMar      -2.049e+02  4.448e+03  -0.046    0.963    
## mesAbr      -2.374e+03  4.448e+03  -0.534    0.594    
## mesMay      -6.806e+02  4.448e+03  -0.153    0.878    
## mesJun       1.251e+03  4.481e+03   0.279    0.780    
## mesJul       8.485e+02  4.481e+03   0.189    0.850    
## mesAgo       7.702e+02  4.481e+03   0.172    0.864    
## mesSep       3.326e+02  4.482e+03   0.074    0.941    
## mesOct       4.715e+02  4.482e+03   0.105    0.916    
## mesNov      -5.188e+02  4.482e+03  -0.116    0.908    
## mesDic      -1.542e+03  4.482e+03  -0.344    0.731    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 18340 on 386 degrees of freedom
## Multiple R-squared:  0.9852, Adjusted R-squared:  0.9847 
## F-statistic:  1839 on 14 and 386 DF,  p-value: < 2.2e-16

Tabla de significancia de meses (enero queda como categoría base)

coef_meses <- tidy(mod_dummy) %>%
  filter(grepl("^mes", term)) %>%
  mutate(significativo_05 = ifelse(p.value < 0.05, "Sí", "No"))
print(coef_meses)

## # A tibble: 11 × 6
##    term   estimate std.error statistic p.value significativo_05
##    <chr>     <dbl>     <dbl>     <dbl>   <dbl> <chr>           
##  1 mesFeb    -867.     4447.   -0.195    0.845 No              
##  2 mesMar    -205.     4448.   -0.0461   0.963 No              
##  3 mesAbr   -2374.     4448.   -0.534    0.594 No              
##  4 mesMay    -681.     4448.   -0.153    0.878 No              
##  5 mesJun    1251.     4481.    0.279    0.780 No              
##  6 mesJul     848.     4481.    0.189    0.850 No              
##  7 mesAgo     770.     4481.    0.172    0.864 No              
##  8 mesSep     333.     4482.    0.0742   0.941 No              
##  9 mesOct     472.     4482.    0.105    0.916 No              
## 10 mesNov    -519.     4482.   -0.116    0.908 No              
## 11 mesDic   -1542.     4482.   -0.344    0.731 No

4) Estacionalidad trigonométrica (3.4) ———————————–

datos <- datos %>%
  mutate(
    sen12 = sin(2 * pi * t / 12),
    cos12 = cos(2 * pi * t / 12)
  )

mod_trig <- lm(y ~ t + I(t^2) + I(t^3) + sen12 + cos12, data = datos)
summary(mod_trig)

## 
## Call:
## lm(formula = y ~ t + I(t^2) + I(t^3) + sen12 + cos12, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -147516   -9581   -2343    9659   46763 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.332e+05  3.659e+03  36.397   <2e-16 ***
## t            1.992e+03  7.872e+01  25.307   <2e-16 ***
## I(t^2)      -8.174e+00  4.547e-01 -17.977   <2e-16 ***
## I(t^3)       1.739e-02  7.436e-04  23.385   <2e-16 ***
## sen12       -8.214e+02  1.280e+03  -0.641    0.522    
## cos12       -4.524e+02  1.283e+03  -0.353    0.725    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 18140 on 395 degrees of freedom
## Multiple R-squared:  0.9852, Adjusted R-squared:  0.985 
## F-statistic:  5260 on 5 and 395 DF,  p-value: < 2.2e-16

comparacion_estacional <- tibble(
  Modelo = c("Cúbico + dummies", "Cúbico + sen/cos"),
  SSE = c(sum(resid(mod_dummy)^2), sum(resid(mod_trig)^2)),
  AIC = c(AIC(mod_dummy), AIC(mod_trig))
)
print(comparacion_estacional)

## # A tibble: 2 × 3
##   Modelo                     SSE   AIC
##   <chr>                    <dbl> <dbl>
## 1 Cúbico + dummies 129795624740. 9028.
## 2 Cúbico + sen/cos 130027755413. 9010.

Visualización de modelos estacionales

ajustes_est <- datos %>%
  mutate(
    dummy = fitted(mod_dummy),
    trig = fitted(mod_trig)
  )

p_est <- ggplot(ajustes_est, aes(x = fecha)) +
  geom_line(aes(y = y, color = "Observado"), linewidth = 0.6) +
  geom_line(aes(y = dummy, color = "Dummies"), linewidth = 0.8) +
  geom_line(aes(y = trig, color = "Trigonométrico"), linewidth = 0.8) +
  labs(title = "Comparación de estacionalidad",
       x = "Fecha", y = "Millones de dólares", color = "Modelo") +
  theme_minimal(base_size = 11)
print(p_est)

5) Curvas de crecimiento (3.5) ——————————————-

Exponencial: y = alpha * exp(beta * t)

mod_exp <- nls(y ~ alpha * exp(beta * t),
               data = datos,
               start = list(alpha = 170000, beta = 0.003))
summary(mod_exp)

## 
## Formula: y ~ alpha * exp(beta * t)
## 
## Parameters:
##        Estimate Std. Error t value Pr(>|t|)    
## alpha 1.738e+05  1.709e+03   101.7   <2e-16 ***
## beta  3.438e-03  3.282e-05   104.8   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 25830 on 399 degrees of freedom
## 
## Number of iterations to convergence: 4 
## Achieved convergence tolerance: 1.399e-06

Asintótica: y = Asym + (R0 - Asym) * exp(-k * t) Este modelo permite interpretar una asíntota de largo plazo (Asym).

mod_asym <- nls(y ~ Asym + (R0 - Asym) * exp(-k * t),
                data = datos,
                start = list(Asym = 900000, R0 = 150000, k = 0.003),
                control = nls.control(maxiter = 500, warnOnly = TRUE))

## Warning in nls(y ~ Asym + (R0 - Asym) * exp(-k * t), data = datos, start =
## list(Asym = 9e+05, : step factor 0.000488281 reduced below 'minFactor' of
## 0.000976562

summary(mod_asym)

## 
## Formula: y ~ Asym + (R0 - Asym) * exp(-k * t)
## 
## Parameters:
##       Estimate Std. Error t value Pr(>|t|)    
## Asym 4.287e+06  6.911e+06   0.620    0.535    
## R0   1.685e+05  7.785e+03  21.640   <2e-16 ***
## k    2.483e-04  4.379e-04   0.567    0.571    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 51130 on 398 degrees of freedom
## 
## Number of iterations till stop: 51 
## Achieved convergence tolerance: 1.505
## Reason stopped: step factor 0.000488281 reduced below 'minFactor' of 0.000976562

Ajustes y comparación de errores

pred_exp  <- predict(mod_exp)
pred_asym <- predict(mod_asym)

comparacion_crecimiento <- tibble(
  Modelo = c("Exponencial", "Asintótico"),
  SSE = c(sum((datos$y - pred_exp)^2), sum((datos$y - pred_asym)^2)),
  AIC_aprox = c(AIC(mod_exp), AIC(mod_asym))
)
print(comparacion_crecimiento)

## # A tibble: 2 × 3
##   Modelo          SSE AIC_aprox
##   <chr>         <dbl>     <dbl>
## 1 Exponencial 2.66e11     9290.
## 2 Asintótico  1.04e12     9838.

Gráfico de curvas de crecimiento

ajustes_crec <- datos %>%
  mutate(
    exp_fit = pred_exp,
    asym_fit = pred_asym
  )

p_crec <- ggplot(ajustes_crec, aes(x = fecha)) +
  geom_line(aes(y = y, color = "Observado"), linewidth = 0.6) +
  geom_line(aes(y = exp_fit, color = "Exponencial"), linewidth = 0.8) +
  geom_line(aes(y = asym_fit, color = "Asintótico"), linewidth = 0.8) +
  labs(title = "Curvas de crecimiento",
       x = "Fecha", y = "Millones de dólares", color = "Modelo") +
  theme_minimal(base_size = 11)
print(p_crec)

Conclusiones e interpretación

3.1- La serie muestra una tendencia claramente creciente de largo plazo 3.2- Entre los modelos polinomiales, el cúbico suele ser el mejor si minimiza AIC y SSE 3.3- Si los meses no son significativos, la estacionalidad residual es débil 3.4- Si el modelo trigonométrico rinde parecido a las ficticias, conviene por parsimonia 3.5- En curvas de crecimiento, beta en la exponencial representa la tasa de crecimiento 3.6- En el modelo asintótico, Asym representa la asíntota de largo plazo

Referencias