Descripción básica

# Lectura de los datos
data(tcm)

# Gráfico
plot(tcm10y,
     main = 'Rendimiento mensual de los bonos del Tesoro de EE.UU.',
     xlab = 'Año',
     ylab = 'Rendimiento [%]',
     lwd = 2)
grid()

De la gráfica de la serie es claro que esta no es estacionaria en sentido débil, puesto que la media no es constante, lo cual se refleja en una tendencia creciente en el rendimiento mensual de los bonos hasta los primeros años de la década de los ochentas, momento a partir del cual se comienza a dar un decrecimiento en el rendimiento mensual. En todo caso, para despejar dudas relacionadas con la existencia de tendencia, se puede ver la gráfica de autocorrelación muestral:

acf(tcm10y,
    main = 'Función de autocorrelación muestral',
    xlab = 'Rezago',
    ylab = 'ACF',
    lwd = 3)
grid()

Como se puede ver, la realización observada de este proceso estacionario no es ergódica, ya que la función de autocorrelación muestral no tiene a cero rápidamente, sino que aún para rezagos grandes, la ACF muestral es significativa estadísticamente. Además, este decrecimiento refleja precisamente la existencia de tendencia en esta serie de tiempo. De igual forma, apelando a la función de autocorrelación parcial muestral:

pacf(tcm10y,
     main = 'Función de autocorrelación parcial muestral',
     xlab = 'Rezago',
     ylab = 'PACF',
     lwd = 3)
grid()

Para este gráfico se evidencia que las series rezagadas uno, dos y tres periodos en el tiempo, ignorando los rezagos internos, son significativas, corroborando que la serie en cuestión no es estacionaria.

Modelo uno

Ahora bien, teniendo en cuenta que la ACF es tipo cola (no amortiguada) y que la PACF es tipo corte, con último corte significativo \(p = 3\), entonces se tiene que es razonable considerar un modelo \(AR(3)\), esto es, si \(X_t\) es el rendimiento mensual del bono público del tesoro estadounidense en un periodo \(t\), entonces el modelo a considerar es:

\[ X_t = \phi_0 + \phi_1 X_{t-1} + \phi_2 X_{t-2} + \phi_3 X_{t-3} + w_t, \ w_t \sim R.B. \]

modelo1 <- arima(tcm10y, order = c(3, 0, 0))
modelo1

## 
## Call:
## arima(x = tcm10y, order = c(3, 0, 0))
## 
## Coefficients:
##          ar1      ar2     ar3  intercept
##       1.4033  -0.6438  0.2351     5.9815
## s.e.  0.0411   0.0669  0.0411     1.6701
## 
## sigma^2 estimated as 0.06829:  log likelihood = -45.39,  aic = 100.79

Ahora bien, el modelo estimado está dado por:

\[ \hat{X}_t = 5.9815 + 1.4033 X_{t-1} - 0.6438 X_{t-2} + 0.2351 X_{t-3} \]

Y gráficamente, se tiene la siguiente serie estimada:

# Serie estimada
modelo1_ajuste <- tcm10y - modelo1$residuals

# Gráfico de la serie original
plot(tcm10y,
     main = 'Rendimiento mensual de los bonos del Tesoro de EE.UU.',
     xlab = 'Año',
     ylab = 'Rendimiento [%]',
     lwd = 3,
     col = 'gray')
grid()

# Gráfico de la serie estimada -- AR(3)

lines(modelo1_ajuste,
      col = 'red',
      lty = 2)

# Leyenda
legend('topleft',
       legend = c('Serie observada',
                  'Serie estimada -- AR(3)'),
       lwd = c(3, 0),
       lty = c(1, 2),
       col = c('gray', 'red'))

Como se puede observar, esta serie parece ajustar de forma casi idéntica la serie ajustada. Esto puede lucir muy bien, sin embargo, indica que existe un sobreajuste en el modelo. De todos modos, es importante chequear que la serie ajustada sea estacionaria, y eso se consigue verificando que las raíces del polinomio \(\phi(B) = 1 - 1.4033 B + 0.6438 B^2 - 0.2351 B^3\) estén por fuera del círculo unitario (en el plano \(\mathbb{R} \ - \ \mathbb{C}\)). Esto se chequea fácilmente con las funciones polyroot y abs de R.

c(1, -1.4033, 0.6438, -0.2351) %>% polyroot %>%  abs

## [1] 1.006574 2.055658 2.055658

Y en efecto, todas las raíces son mayores estrictamente que uno (de momento no se considerará si esto se puede considerar de forma significativa, en especial para la raíz más pequeña que resultó muy cercana a uno), lo cual indica que el proceso no es estacionario.

Modelo dos

Para solventar lo anterior, se va a considerar la primera diferencia regular del modelo, la cual elimina la tendencia, esto es, se va a modelar la serie \(Y_t = X_{t} - X_{t-1}\).

# Serie con primera diferencia regular
y <- diff(tcm10y)

# Gráfico
plot(y,
     main = 'Primera diferencia regular del rendimiento mensual
     de los bonos del Tesoro de EE.UU.',
     xlab = 'Año',
     ylab = 'Rendimiento [%]',
     lwd = 1)
grid()

Y esta serie, si bien aún tiene problemas relacionados con la estacionariedad puesto que la varianza no puede ser considerada constante, sí solventa uno de los problemas considerados inicialmente y es que la media es estacionaria y nula. Así, vale la pena comenzar revisando la ACF:

acf(y,
    main = 'Función de autocorrelación muestral',
    xlab = 'Rezago',
    ylab = 'ACF',
    lwd = 3)
grid()

Y como se evidencia, se tiene corte en \(p = 3\). Ahora bien, para la PACF:

pacf(y,
     main = 'Función de autocorrelación parcial muestral',
     xlab = 'Rezago',
     ylab = 'PACF',
     lwd = 3)
grid()

Se observa un patrón tipo cola sinusoidal amortiguado. Esto indica que es razonable considerar un modelo de medias móviles de orden tres, cuya ecuación viene dada por:

\[ X_t = \theta_0 + \theta_1 w_{t-1} + \theta_2 w_{t-2} + \theta_3 w_{t-3} + w_t, \ w_t \sim R.B. \]

modelo2 <- arima(y, order = c(0, 0, 3))
modelo2

## 
## Call:
## arima(x = y, order = c(0, 0, 3))
## 
## Coefficients:
##          ma1      ma2      ma3  intercept
##       0.4494  -0.1219  -0.0077     0.0055
## s.e.  0.0423   0.0467   0.0386     0.0145
## 
## sigma^2 estimated as 0.06698:  log likelihood = -37.61,  aic = 85.23

Lo que quiere decir que el modelo ajustado es:

\[ \hat{X}_t = 0.0055 + 0.4494 w_{t-1} - 0.1219 w_{t-2} - 0.0077 w_{t-3} \] Gráficamente:

# Serie estimada
modelo2_ajuste <- y - modelo2$residuals

# Gráfico de la serie original
plot(y,
     main = 'Primera diferencia regular del rendimiento
     mensual de los bonos del Tesoro de EE.UU.',
     xlab = 'Año',
     ylab = 'Rendimiento [%]',
     lwd = 3,
     col = 'gray')
grid()

# Gráfico de la serie estimada -- AR(3)

lines(modelo2_ajuste,
      col = 'red',
      lty = 2)

# Leyenda
legend('topleft',
       legend = c('Serie observada',
                  'Serie estimada -- AR(3)'),
       lwd = c(3, 0),
       lty = c(1, 2),
       col = c('gray', 'red'))

Ejemplo dos

Ahora se va a considerar la serie asociada a los incidentes viales en la ciudad de Medellín asociados a motos y que fueron registrados entre 2015 y 2021.

direccion <- './incidentes_viales_motos.csv'
bd_incidentes <- read_delim(direccion,
                            delim = ';')
bd_incidentes %<>% clean_names()
bd_incidentes %>% dim()

## [1] 223439      9

bd_incidentes %>% head(3)

## # A tibble: 3 × 9
##   nro_radicado ano_incidente fecha_incidente hora_incidente clase_incidente
##          <dbl>         <dbl> <chr>           <time>         <chr>          
## 1      1473523          2015 27/01/15        22:00          Choque         
## 2      1473525          2015 27/01/15        15:40          Choque         
## 3      1473526          2015 27/01/15        18:00          Choque         
## # … with 4 more variables: gravedad_incidente <chr>, direccion <chr>,
## #   zona <chr>, diseno_via <chr>

bd_incidentes %>% tail(n = 3)

## # A tibble: 3 × 9
##   nro_radicado ano_incidente fecha_incidente hora_incidente clase_incidente
##          <dbl>         <dbl> <chr>           <time>         <chr>          
## 1      1759665          2021 25/08/2021      20:15          Choque         
## 2      1759743          2021 25/08/2021      20:25          Choque         
## 3      1759627          2021 25/08/2021      21:00          Otro           
## # … with 4 more variables: gravedad_incidente <chr>, direccion <chr>,
## #   zona <chr>, diseno_via <chr>

A continuación se va a tratar el problema de las fechas:

aux1 <- as.Date(rep(NA, nrow(bd_incidentes)))
aux2 <- bd_incidentes$fecha_incidente

ind1 <- which(nchar(aux2) == 8)
ind2 <- which(nchar(aux2) == 10)
              
aux1[ind1] <- as.Date(aux2[ind1], format = '%d/ %m/ %y')
aux1[ind2] <- as.Date(aux2[ind2], format = '%d/ %m/ %Y')

bd_incidentes$fecha_incidente <- aux1

bd_incidentes$mes <- month(bd_incidentes$fecha_incidente)
bd_incidentes$ano <- year(bd_incidentes$fecha_incidente)

resumen1 <- bd_incidentes %>% 
  group_by(ano, mes) %>% 
  summarise(accidentes = n())

resumen1$fecha <- as.Date(paste(resumen1$ano, resumen1$mes, 1, sep = "-"),
                          format = '%Y- %m- %d')

resumen1 %>% ggplot(aes(x = fecha, y = accidentes)) +
  geom_line(lwd = 1) +
  ggtitle('Incidentes viales por mes en Medellín con motos') +
  xlab('Fecha') + ylab('Incidentes') +
  theme_light() -> fig1

fig1 %>% plotly::ggplotly()

acf(resumen1$accidentes,
    main = 'Función de autocorrelación muestral',
    xlab = 'Rezago',
    ylab = 'ACF',
    lwd = 3)
grid()

Y como se evidencia, tiene forma tipo cola. Luego, con respecto a la PACF:

pacf(resumen1$accidentes,
     main = 'Función de autocorrelación parcial muestral',
     xlab = 'Rezago',
     ylab = 'PACF',
     lwd = 3)
grid()

modelo3 <- arima(resumen1$accidentes,
                 order =  c(6, 0, 0))
modelo3

## 
## Call:
## arima(x = resumen1$accidentes, order = c(6, 0, 0))
## 
## Coefficients:
##          ar1     ar2      ar3      ar4     ar5      ar6  intercept
##       0.7024  0.1013  -0.0705  -0.0515  0.2258  -0.2969  2794.4827
## s.e.  0.1055  0.1307   0.1328   0.1328  0.1409   0.1117    83.4591
## 
## sigma^2 estimated as 83525:  log likelihood = -567.52,  aic = 1151.04

modelo3$var.coef

##                     ar1           ar2           ar3           ar4           ar5
## ar1        0.0111371854 -0.0079573683 -1.530962e-03  0.0006076713  0.0010697036
## ar2       -0.0079573683  0.0170784445 -7.253965e-03 -0.0013472436 -0.0007386468
## ar3       -0.0015309624 -0.0072539653  1.763696e-02 -0.0075239542 -0.0011924228
## ar4        0.0006076713 -0.0013472436 -7.523954e-03  0.0176447174 -0.0083707326
## ar5        0.0010697036 -0.0007386468 -1.192423e-03 -0.0083707326  0.0198445167
## ar6        0.0001320265  0.0012186402  1.168201e-05 -0.0005955539 -0.0098152564
## intercept -0.0519470861 -0.0231660516  8.702025e-02  0.1867075033 -0.0469766878
##                     ar6     intercept
## ar1        1.320265e-04   -0.05194709
## ar2        1.218640e-03   -0.02316605
## ar3        1.168201e-05    0.08702025
## ar4       -5.955539e-04    0.18670750
## ar5       -9.815256e-03   -0.04697669
## ar6        1.248158e-02   -0.05590311
## intercept -5.590311e-02 6965.41814376

error_estandar <- modelo3$var.coef %>% diag() %>% sqrt() %>% round(4)

modelo3$coef - 1.96 * error_estandar

##           ar1           ar2           ar3           ar4           ar5 
##    0.49562913   -0.15488017   -0.33079428   -0.31174290   -0.05035449 
##           ar6     intercept 
##   -0.51587748 2630.90282921

modelo3$coef + 1.96 * error_estandar

##           ar1           ar2           ar3           ar4           ar5 
##    0.90918913    0.35746383    0.18978172    0.20883310    0.50197351 
##           ar6     intercept 
##   -0.07801348 2958.06250121