Codigo Parcial 1 Series de Tiempo

Punto 1

# AutoCovarianza
rk <- function(k){
  if(k==0){r=31.2}
  if (k==1) {r=14.4}
  if (k==2) {r=20.4}
  if (k==3 | k == 4) {r=4.8}
  if(k>=5){r=0}
  return(r)
}

for (i in 0:7) {
  print(rk(i))
}

[1] 31.2
[1] 14.4
[1] 20.4
[1] 4.8
[1] 4.8
[1] 0
[1] 0
[1] 0

# AutoCorrelación
pk <- function(k){
  p=rk(k)/rk(0)
  return(p)
}


for (i in 0:7) {
  print(pk(i))
}

[1] 1
[1] 0.4615385
[1] 0.6538462
[1] 0.1538462
[1] 0.1538462
[1] 0
[1] 0
[1] 0

# Gráfica del ACF
ma.acf <- vector()

k =5

for(i in 0:k){
  ma.acf[i+1] <- pk(i)
}

# PACF
pacf_ma <- vector()
pacf_ma[1] <- ma.acf[1]
for (i in 2:(k+1)){
  deno <- toeplitz(c(1, ma.acf[1:(i-1)]))
  aux_1 <- deno
  aux_1[,i] <- ma.acf[1:i]
  nume <- aux_1
  pacf_ma[i] <- det(nume)/det(deno)
  if(i==2){
    pacf_ma[i] <- (pk(2)-(pk(1)^2))/(1-(pk(1)^2))
  }
}

par(mfrow=c(1,2))
barplot(ma.acf, main="ACF teórica de Xt", las=1, ylim=c(-1,1),
        names.arg = c(1,2,3,4,5,6))
abline(h=1.96/sqrt(200), lty=2)
abline(h=0)
abline(h=-1.96/sqrt(200), lty=2)
box()
barplot(pacf_ma, main="PACF teórica de Xt", las=1, ylim=c(-1,1),
        names.arg = c(1,2,3,4,5,6))
abline(h=1.96/sqrt(200), lty=2)
abline(h=0)
abline(h=-1.96/sqrt(200), lty=2)
box()

###### Simulación del Xt
w <- rnorm(205,mean = 0,sd = sqrt(4.8))

Xt <- c()

for(t in 3:203){
  Xt[t-3] <- w[t-2]+0.5*w[t-1]+2*w[t]+0.5*w[t+1]+w[t+2]
}

par(mfrow=c(1,1))
plot(1:200,Xt, type = "l", main="Serie de Tiempo Xt (Simulada)",
     ylab=expression(X[t] == w[t-2]+0.5*yw[t-1]+2*w[t]+0.5*w[t+1]+w[t+2]),
     xlab = "t")


library(TSstudio)
library(tseries)

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

library(magrittr)

# Simulado

par(mfrow=c(1,2))
Xt %>% acf(plot = F) %>% plot(main="")
Xt %>% pacf(plot = F) %>% plot(main="")
mtext("ACF & PACF Simulados para Xt",side = 3,line = - 2,outer = TRUE)

# Teórico vs Simulado

par(mfrow=c(1,2))
barplot(ma.acf, main="ACF teórica", las=1, ylim=c(-1,1),
        names.arg = c(1,2,3,4,5,6))
abline(h=1.96/sqrt(200), lty=2)
abline(h=0)
abline(h=-1.96/sqrt(200), lty=2)
box()
acf(Xt, lag.max = 6, plot = F) %>% plot(main="ACF Simulada")

par(mfrow=c(1,2))
barplot(pacf_ma[-1], main="PACF teórica de Xt", las=1, ylim=c(-1,1),
        names.arg = c(1,2,3,4,5))
abline(h=1.96/sqrt(200), lty=2)
abline(h=0)
abline(h=-1.96/sqrt(200), lty=2)
box()
pacf(Xt, lag.max = 6, plot = F) %>% plot(main="PACF Simulada")

Punto 2

library(magrittr)
# Creación del ACF Teórico:
par(mfrow=c(1,2))
ARMAacf(ar=c(0.9,-0.6), lag.max = 20) %>%
  barplot(.,las=1,ylim=c(-1,1),main="",ylab="ACF")
  abline(h=1.96/sqrt(180), lty=2)
  abline(h=0)
  abline(h=-1.96/sqrt(180), lty=2)
  box()

# Creación del PACF Teórico:
ARMAacf(ar=c(0.9,-0.6), lag.max = 20, pacf = T) %>%
  barplot(.,las=1,ylim=c(-1,1),main="",ylab="PACF")
  abline(h=1.96/sqrt(180), lty=2)
  abline(h=0)
  abline(h=-1.96/sqrt(180), lty=2)
  box()

mtext("ACF & PACF Teóricos para Xt",side = 3,line = - 2,outer = TRUE)

# Simulación del Xt:
Xt <- arima.sim(model=list(order=c(2,0,0), ar=c(0.9,-0.6)), n=180, sd=sqrt(6.2))+3.1

par(mfrow=c(1,1))
plot(1:180,Xt, type = "l", main="Serie de Tiempo Xt (Simulada)",
     ylab = expression(X[t]==3.1+0.9*X[t-1]-0.6*X[t-2]+w[t]),
     xlab = "t")

par(mfrow=c(1,2))
Xt %>% acf(plot = F) %>% plot(main="")
Xt %>% pacf(plot = F) %>% plot(main="")
mtext("ACF & PACF Simulados para Xt",side = 3,line = - 2,outer = TRUE)

Comparación de lo teórico y lo simulado:

# ACF (vs)
par(mfrow=c(1,2))
ARMAacf(ar=c(0.9,-0.6), lag.max = 20) %>% barplot(.,las=1,ylim=c(-1,1),main="ACF Teórico")+box()

numeric(0)

abline(h=1.96/sqrt(180), lty=2)
abline(h=0)
abline(h=-1.96/sqrt(180), lty=2)

Xt %>% acf(plot = F) %>% plot(main="ACF Simulado")

# PACF (vs)
par(mfrow=c(1,2))
ARMAacf(ar=c(0.9,-0.6), lag.max = 20, pacf = T) %>%
  barplot(.,las=1,ylim=c(-1,1),main="PACF Teórico",
          names.arg=1:20)
  box()
  abline(h=1.96/sqrt(180), lty=2)
  abline(h=0)
  abline(h=-1.96/sqrt(180), lty=2)

Xt %>% pacf(plot = F) %>% plot(main="PACF Simulado")

Punto 3:

library(readxl)
library(tidyverse)

── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
✔ ggplot2 3.3.6     ✔ purrr   0.3.4
✔ tibble  3.1.8     ✔ dplyr   1.0.9
✔ tidyr   1.2.0     ✔ stringr 1.4.0
✔ readr   2.1.2     ✔ forcats 0.5.1
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ tidyr::extract()   masks magrittr::extract()
✖ dplyr::filter()    masks stats::filter()
✖ dplyr::lag()       masks stats::lag()
✖ purrr::set_names() masks magrittr::set_names()

library(magrittr)
library(janitor)

Attaching package: 'janitor'

The following objects are masked from 'package:stats':

    chisq.test, fisher.test

library(lubridate)

Attaching package: 'lubridate'

The following objects are masked from 'package:base':

    date, intersect, setdiff, union

library(naniar)
library(radiant)

Loading required package: radiant.data

Attaching package: 'radiant.data'

The following objects are masked from 'package:lubridate':

    month, wday

The following object is masked from 'package:forcats':

    as_factor

The following objects are masked from 'package:purrr':

    is_double, is_empty, is_numeric

The following object is masked from 'package:ggplot2':

    diamonds

The following object is masked from 'package:magrittr':

    set_attr

The following object is masked from 'package:base':

    date

Loading required package: radiant.design
Loading required package: mvtnorm
Loading required package: radiant.basics
Loading required package: radiant.model
Loading required package: radiant.multivariate

Contexto de la base de datos:

Según la pagina de datos abiertos del metro de Medellín ellos definen afluencia como: “Cantidad de viajeros que acuden a un lugar o sitio determinado para acceder a un medio de transporte.”

Con esto en mano podemos ver nuestra base datos:

Datos de Afluencia por linea, hora y día, con total de pasajeros en ese día:

m19 <- read_excel("Afluencia_2019.xlsx")
m20 <- read_excel("Afluencia_2020.xlsx")
m21 <- read_excel("Afluencia_2021.xlsx")

Exploración y determinación para organizar la base de datos:

# Limpiar variables:
m19 %<>% clean_names()
m20 %<>% clean_names()
m21 %<>% clean_names()

####
m19 %>% head(5)

m19 %>% str()

tibble [3,563 × 23] (S3: tbl_df/tbl/data.frame)
 $ dia                              : POSIXct[1:3563], format: "2019-01-01" "2019-01-02" ...
 $ linea_de_servicio                : chr [1:3563] "Línea 1" "Línea 1" "Línea 1" "Línea 1" ...
 $ x4am                             : num [1:3563] 12 1716 1860 1845 1550 ...
 $ x5am                             : num [1:3563] 826 4498 4869 4989 3852 ...
 $ x6am                             : num [1:3563] 1134 6629 7587 7723 4898 ...
 $ x7am                             : num [1:3563] 987 6284 7088 7281 5221 ...
 $ x8am                             : num [1:3563] 884 4145 4483 4857 4016 ...
 $ x9am                             : num [1:3563] 975 3784 4486 4583 3815 ...
 $ x10am                            : num [1:3563] 1259 3746 4126 4269 3704 ...
 $ x11am                            : num [1:3563] 1348 3854 4136 4217 4048 ...
 $ x12am                            : num [1:3563] 1678 3787 4210 4209 4800 ...
 $ x1pm                             : num [1:3563] 1861 3760 4007 4261 4893 ...
 $ x2pm                             : num [1:3563] 1891 4628 4916 5002 5116 ...
 $ x3pm                             : num [1:3563] 2085 4709 5061 5280 4696 ...
 $ x4pm                             : num [1:3563] 2230 5750 6093 6056 4528 ...
 $ x5pm                             : num [1:3563] 2636 7517 8023 7992 4106 ...
 $ x6pm                             : num [1:3563] 2884 6661 7034 7081 4377 ...
 $ x7pm                             : num [1:3563] 2603 4307 4481 4805 3842 ...
 $ x8pm                             : num [1:3563] 2024 2975 3244 3466 2973 ...
 $ x9pm                             : num [1:3563] 1257 2709 2905 3075 2866 ...
 $ x10pm                            : num [1:3563] 271 1890 1739 1900 1983 ...
 $ x11pm                            : num [1:3563] 5 135 124 176 167 2 2 79 90 90 ...
 $ total_general_numero_de_pasajeros: num [1:3563] 28850 83484 90472 93067 75451 ...

m20 %>% head(5)

m20 %>% str()

tibble [3,764 × 23] (S3: tbl_df/tbl/data.frame)
 $ dia                              : POSIXct[1:3764], format: "2020-01-01" "2020-01-01" ...
 $ linea_de_servicio                : chr [1:3764] "Línea 1" "Línea 2" "Línea A" "Línea B" ...
 $ x4am                             : num [1:3764] NA NA 1 NA NA NA NA NA NA NA ...
 $ x5am                             : num [1:3764] 902 63 4901 750 NA ...
 $ x6am                             : num [1:3764] 1126 94 6168 863 NA ...
 $ x7am                             : num [1:3764] 881 66 4822 674 NA ...
 $ x8am                             : num [1:3764] 743 68 4158 693 1 ...
 $ x9am                             : num [1:3764] 924 61 5152 891 39 ...
 $ x10am                            : num [1:3764] 1137 79 6467 1161 60 ...
 $ x11am                            : num [1:3764] 1299 127 8488 1600 65 ...
 $ x12am                            : num [1:3764] 1614 112 10554 2023 51 ...
 $ x1pm                             : num [1:3764] 1804 134 11800 1994 78 ...
 $ x2pm                             : num [1:3764] 1868 160 12717 2281 113 ...
 $ x3pm                             : num [1:3764] 2056 169 12785 2140 66 ...
 $ x4pm                             : num [1:3764] 2374 188 15911 2775 130 ...
 $ x5pm                             : num [1:3764] 2727 287 19522 2942 105 ...
 $ x6pm                             : num [1:3764] 2908 302 20582 3083 136 ...
 $ x7pm                             : num [1:3764] 2978 236 20358 2319 121 ...
 $ x8pm                             : num [1:3764] 2457 188 16450 1738 49 ...
 $ x9pm                             : num [1:3764] 1369 85 10701 725 14 ...
 $ x10pm                            : num [1:3764] 174 16 1177 39 2 ...
 $ x11pm                            : num [1:3764] NA NA NA NA NA NA NA NA NA NA ...
 $ total_general_numero_de_pasajeros: num [1:3764] 29341 2435 192714 28691 1030 ...

m21 %>% head(5)

m21 %>% str()

tibble [4,122 × 23] (S3: tbl_df/tbl/data.frame)
 $ dia                              : POSIXct[1:4122], format: "2021-01-01" "2021-01-01" ...
 $ linea_de_servicio                : chr [1:4122] "LÍNEA A" "LÍNEA B" "LÍNEA K" "LÍNEA J" ...
 $ x4am                             : num [1:4122] NA NA NA NA NA NA NA NA NA NA ...
 $ x5am                             : num [1:4122] 3487 532 NA NA 614 ...
 $ x6am                             : num [1:4122] 3405 438 NA NA 590 ...
 $ x7am                             : num [1:4122] 2011 267 NA NA 355 ...
 $ x8am                             : num [1:4122] 1148 164 40 1 240 ...
 $ x9am                             : num [1:4122] 974 145 131 77 217 16 104 31 8 46 ...
 $ x10am                            : num [1:4122] 888 135 104 74 198 16 127 37 15 17 ...
 $ x11am                            : num [1:4122] 1011 142 156 92 231 ...
 $ x12am                            : num [1:4122] 1201 149 160 112 212 ...
 $ x1pm                             : num [1:4122] 1285 165 115 76 302 ...
 $ x2pm                             : num [1:4122] 1571 155 131 118 327 ...
 $ x3pm                             : num [1:4122] 1833 212 167 146 394 ...
 $ x4pm                             : num [1:4122] 2900 321 258 179 579 51 274 52 23 57 ...
 $ x5pm                             : num [1:4122] 3269 421 189 132 624 ...
 $ x6pm                             : num [1:4122] 3314 432 122 135 574 ...
 $ x7pm                             : num [1:4122] 2273 215 103 69 399 ...
 $ x8pm                             : num [1:4122] 1329 157 38 42 231 ...
 $ x9pm                             : num [1:4122] 12 NA NA 3 17 NA 15 3 NA 1 ...
 $ x10pm                            : num [1:4122] NA NA NA NA NA NA NA NA NA NA ...
 $ x11pm                            : num [1:4122] NA NA NA NA NA NA NA NA NA NA ...
 $ total_general_numero_de_pasajeros: num [1:4122] 31911 4050 1714 1256 6104 ...

Convertir a factor las lineas y estandarizarlas.

(2019) Valores Faltantes:

# Porcentajes de valores faltantes según el número de observaciones
vis_miss(m19)

Warning: `gather_()` was deprecated in tidyr 1.2.0.
Please use `gather()` instead.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

vis_miss(m20)

vis_miss(m21)

Los missing values solo representan el 6.7% de los datos totales.

# Número de valores faltantes por variable
gg_miss_var(m19)

Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
use `guide = "none"` instead.

gg_miss_var(m20)

Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
use `guide = "none"` instead.

gg_miss_var(m21)

Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
use `guide = "none"` instead.

Lo hora que más registra valores NA es a las 11PM

¿Que hacer con los valores NA?

1. Horarios del metro y comparar. (Festivos, Domingos, días que no funciona x o y linea, horarios por linea)

2. Suponer circunstancias dentro de los horarios del metro (dejar NA) y para los valores fuera de los horarios de ese año poner 0

# Quitar acentos y poner todo en minúscula
m19$linea_de_servicio <- chartr("áéíóú", "aeiou",
                                tolower(m19$linea_de_servicio))

# Factor la variable "Linea"
m19$linea_de_servicio %<>% gsub(" ","_",.) %>% as.factor()

# Nuevo nombre de variables
names(m19)[2] <- "linea"
names(m19)[23] <- "total_n_pasajeros"

m20$linea_de_servicio <- chartr("áéíóú", "aeiou",
                                tolower(m20$linea_de_servicio))

# Factor la variable "Linea"
m20$linea_de_servicio %<>% gsub(" ","_",.) %>% as.factor()

# Nuevo nombre de variables
names(m20)[2] <- "linea"
names(m20)[23] <- "total_n_pasajeros"

m21$linea_de_servicio <- chartr("áéíóú", "aeiou",
                                tolower(m21$linea_de_servicio))

# Factor la variable "Linea"
m21$linea_de_servicio %<>% gsub(" ","_",.) %>% as.factor()

# Nuevo nombre de variables
names(m21)[2] <- "linea"
names(m21)[23] <- "total_n_pasajeros"

Literal A:

Dimensiones de las bases de datos

m19 %>% dim()

[1] 3563   23

m20 %>% dim()

[1] 3764   23

m21 %>% dim()

[1] 4122   23

Literal B:

Unión de las bases de datos

datos_juntos <- merge(merge(m19,m20, all=T), m21, all=T)

datos_juntos %>% dim()

[1] 11449    23

Literal C:

Nuevo data frame

# Data frame temporal
temp <- gather(datos_juntos,"hora", "n_pasajerosXhora",
               c(-dia,-linea,-total_n_pasajeros))

# Orden de las horas
orden.hora <- paste0("x",c(4:12,1:11),c(rep("am",9),rep("pm",11)))

or.dsem <- c("lunes", "martes", "miércoles",
             "jueves", "viernes", "sábado", "domingo")

meses <- c("Enero","Febrero","Marzo","Abril","Mayo","Junio","Julio",
           "Agosto","Septiembre","Octubre","Noviembre","Diciembre")

# Conversión y creación de nuevas variables:
temp$fecha <- temp$dia
temp$hora %<>% factor(.,levels = orden.hora)
temp$dia <- day(temp$fecha) %>% as.factor()
temp$dsem <- temp$fecha %>% weekdays() %>% factor(., levels = or.dsem)
temp$sem <- temp$fecha %>% week()
temp$sem %<>% factor(.,)
temp$mes <- temp$fecha %>% month()
temp$mes %<>% factor(., labels = meses)
temp$year <- temp$fecha %>% year()
temp$year %<>% as.factor()
temp %>% head(5)

temp %>% str()

'data.frame':   228980 obs. of  10 variables:
 $ dia              : Factor w/ 31 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 2 ...
 $ linea            : Factor w/ 12 levels "linea_1","linea_2",..: 1 2 3 4 5 6 7 8 11 1 ...
 $ total_n_pasajeros: num  28850 2281 183664 25852 1080 ...
 $ hora             : Factor w/ 20 levels "x4am","x5am",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ n_pasajerosXhora : num  12 NA 14 NA NA ...
 $ fecha            : POSIXct, format: "2019-01-01" "2019-01-01" ...
 $ dsem             : Factor w/ 7 levels "lunes","martes",..: 2 2 2 2 2 2 2 2 2 3 ...
 $ sem              : Factor w/ 53 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ mes              : Factor w/ 12 levels "Enero","Febrero",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ year             : Factor w/ 3 levels "2019","2020",..: 1 1 1 1 1 1 1 1 1 1 ...

# Datos 
dc <- temp %>% dplyr::select(c("fecha", "hora", "dia", "dsem",
                         "sem", "mes", "year","n_pasajerosXhora",
                         "total_n_pasajeros", "linea"))

dc %>% dim()

[1] 228980     10

Literal D:

2 nuevos data frames por Linea A y B

dat_lin_A <- dc %>% filter(.,linea=="linea_a") %>% dplyr::select(c(-10))
dat_lin_B <- dc %>% filter(.,linea=="linea_b") %>% dplyr::select(c(-10))

dat_lin_A %>% dim()

[1] 21860     9

dat_lin_B %>% dim()

[1] 21860     9

# Los Missing values no son significativos:
gg_miss_var(dat_lin_A)

Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
use `guide = "none"` instead.

vis_miss(dat_lin_A)

Literal E:

Transformación de un data frame

dc$Pandemia <- ifelse(dc$fecha < "2020-03-23", "antes", "después")
dc$Pandemia <- factor(dc$Pandemia, levels = c("antes", "después"))

tabla de resumen:

result <- pivotr(
  dc, 
  cvars = c("hora", "dsem", "Pandemia", "linea"), 
  nvar = "n_pasajerosXhora", 
  data_filter = "linea %in% c('linea_a', 'linea_b')", 
  nr = Inf
)
# summary()
#dtab(result, dec = 2) %>% render()
r_ajs_total <- result$tab
register("r_ajs_total")


# Filtro

r_ajs_total %>% filter(.,Pandemia=="antes") %>%
  dplyr::select(c(-2))

Resúmenes Estadísticos:

# Promedio por día de la semana y la hora en la linea A
r1.promA <- dat_lin_A %>% filter(.,fecha<"2020-03-23") %>%
  aggregate(n_pasajerosXhora~dsem+hora,mean)

r2.promA <- dat_lin_A %>% filter(.,fecha>="2020-03-23") %>%
  aggregate(n_pasajerosXhora~dsem+hora,mean)

Se puede mirar que no registran promedios los días miércoles y sábados

# Promedio por día de la semana y la hora en la linea B
r1.promB <- dat_lin_B %>% filter(.,fecha<"2020-03-23") %>%
  aggregate(n_pasajerosXhora~dsem+hora,mean)

r2.promB <- dat_lin_B %>% filter(.,fecha>="2020-03-23") %>%
  aggregate(n_pasajerosXhora~dsem+hora,mean)

Se puede mirar que no registran promedios los días miércoles y sábados

r1PA_ajs

result <- pivotr(
  r1.promA, 
  cvars = c("dsem", "hora"), 
  nvar = "n_pasajerosXhora", 
  nr = Inf
)
# summary()
#dtab(result) %>% render()
r1PA_ajs <- result$tab
print(r1PA_ajs)

    hora      lunes     martes miércoles     jueves    viernes     sábado
1   x4am 12604.7812 14031.4219 13019.938 13386.7778 13622.4603 10916.0952
2   x5am 40700.4062 47065.6094 44013.828 44527.8095 45018.5238 29231.5397
3   x6am 55943.3750 64057.0156 60102.391 61765.0794 62941.3810 37216.2857
4   x7am 48081.0312 55371.7500 51198.688 51654.8750 52492.7500 37450.2381
5   x8am 29665.7969 33658.2812 30925.344 33269.2344 32650.8095 28640.4688
6   x9am 26509.5781 30619.7619 27873.129 29003.5000 29342.2812 25232.5312
7  x10am 25159.2951 27009.7500 26542.194 27665.3438 29201.7969 25211.5469
8  x11am 26659.3810 29697.8906 31597.953 27499.0469 30350.9062 29104.8254
9  x12am 29071.6667 32012.5469 30337.969 30850.3750 31639.7344 43610.0476
10  x1pm 30919.6094 31605.3281 30033.969 30641.5781 31755.0938 46326.6667
11  x2pm 34401.0625 36172.6719 35248.094 34431.2812 38454.0938 46777.0317
12  x3pm 36483.7500 39162.6406 37229.250 37210.2258 41562.5312 43248.9062
13  x4pm 49518.9062 55519.2812 52362.375 55365.1905 58713.0312 38730.7656
14  x5pm 71095.6406 79922.9375 74257.672 81557.3281 76402.8438 35820.6250
15  x6pm 57980.3750 64876.2969 60206.781 62729.6094 60449.8906 33434.5312
16  x7pm 34127.5156 37457.8281 35539.047 37714.4219 38199.3125 28215.9844
17  x8pm 22440.5156 24316.5469 27196.875 26240.2812 25596.9531 20914.7344
18  x9pm 18385.4219 20359.8594 19652.891 20610.9688 22117.1875 18833.0312
19 x10pm  8937.6562 10552.5079 10374.016 11390.4844 13736.3281 12391.8438
20 x11pm   131.6316   131.8852  2889.607   551.5469   505.0476   800.1111
21 Total 32940.8698 36680.0906 35030.100 35903.2479 36737.6478 29605.3905
     domingo      Total
1   3129.984 11530.2083
2  10516.266 37296.2832
3  13115.203 50734.3901
4  12281.375 44075.8153
5  12395.453 28743.6268
6  14033.141 26087.7032
7  15952.938 25248.9805
8  15996.484 27272.3554
9  17785.766 30758.3007
10 20107.766 31627.1443
11 21689.328 35310.5090
12 21031.328 36561.2332
13 22474.562 47526.3018
14 23371.094 63204.0201
15 23912.297 51941.3973
16 20407.547 33094.5223
17 18420.797 23589.5290
18 14429.609 19198.4241
19  1488.625  9838.7802
20    30.000   719.9756
21 15128.478 31717.9750

register("r1PA_ajs")

r2PA_ajs

result <- pivotr(
  r2.promA, 
  cvars = c("dsem", "hora"), 
  nvar = "n_pasajerosXhora", 
  nr = Inf
)
# summary()
#dtab(result) %>% render()
r2PA_ajs <- result$tab
print(r2PA_ajs)

    hora       lunes      martes   miércoles      jueves     viernes     sábado
1   x4am 10029.45652 11639.64516 11598.73118 11662.64130 10996.39130  9097.8876
2   x5am 25103.59783 29101.09677 29311.49462 29509.05435 27639.62366 20406.1444
3   x6am 30749.33696 36101.65591 36278.38710 36267.10870 34046.43011 23488.0111
4   x7am 24248.00000 28161.19355 28243.73118 28182.01087 26325.69892 20013.1778
5   x8am 16004.90217 18532.07527 18355.60215 18285.63043 17282.60215 14414.2778
6   x9am 13613.44565 15444.45161 15286.73118 15210.53261 14475.50538 12723.7889
7  x10am 12595.41304 13882.05376 13860.64516 13822.40217 13237.27957 12719.2333
8  x11am 12729.78261 13913.55914 13962.64516 13949.82609 13490.77419 14609.1889
9  x12am 14571.18478 16153.51613 16153.33333 16189.16304 15602.47312 23017.2667
10  x1pm 14699.56522 15947.33333 16038.07527 16192.17391 15700.66667 26121.9444
11  x2pm 18807.65217 20115.25806 20480.34409 20335.03261 20447.79570 27319.2778
12  x3pm 20128.86957 22107.30108 22672.21505 22619.50000 23029.93548 22621.1000
13  x4pm 30530.34783 35087.77419 35041.98925 35219.15217 34350.26882 21806.3667
14  x5pm 38784.56522 44966.19355 44431.87097 44851.43478 41391.37634 20039.0889
15  x6pm 30286.57609 35066.23656 34652.38710 35317.31522 32623.00000 18632.6000
16  x7pm 17127.46739 19283.23656 19227.59140 19834.10870 19600.45161 14872.3820
17  x8pm 10766.92391 11751.61290 11788.94624 12081.34783 12255.09677 10736.2360
18  x9pm  8320.97826  9254.75269  9506.89130  9511.64130  9809.18280  9241.2697
19 x10pm  3381.35870  4136.50538  4262.44565  4415.16484  4845.49451  4566.1250
20 x11pm    31.30645    34.91429    36.10959    41.52174    69.56061   170.7692
21 Total 17625.53652 20034.01829 20059.50835 20174.83813 19360.98039 16330.8068
        domingo       Total
1   2629.703297  9664.92234
2   7381.217391 24064.60415
3   8406.695652 29333.94650
4   6813.630435 23141.06325
5   6430.195652 15615.04080
6   7438.934783 13456.19859
7   8074.554348 12598.79734
8   8343.858696 12999.94782
9   9171.532609 15836.92424
10 10207.260870 16415.28853
11 11769.119565 19896.35428
12 11109.782609 20612.67197
13 12045.434783 29154.47624
14 12337.130435 35257.38003
15 11862.195652 28348.61580
16  9904.858696 17121.44234
17  8574.065217 11136.31840
18  6358.902174  8857.65974
19   578.402174  3740.78518
20     2.153846    55.19082
21  7971.981444 17365.38142

register("r2PA_ajs")

r1PB_ajs

result <- pivotr(
  r1.promB, 
  cvars = c("dsem", "hora"), 
  nvar = "n_pasajerosXhora", 
  nr = Inf
)
# summary()
#dtab(result) %>% render()
r1PB_ajs <- result$tab
print(r1PB_ajs)

    hora     lunes      martes miércoles      jueves     viernes    sábado
1   x4am 2342.9062  2653.19048 2474.3651  2475.95238  2486.47619 2026.2063
2   x5am 6590.9375  7614.31250 7142.5000  7216.06349  7333.63492 4807.2063
3   x6am 8309.7969  9427.92188 8899.3438  9032.87302  9326.01587 5707.8254
4   x7am 7642.7656  8840.82812 8202.7656  8273.00000  8379.10938 6014.8413
5   x8am 4918.0000  5586.79688 5162.6719  5569.87500  5524.57143 4889.1094
6   x9am 4334.3750  5032.79365 4571.6613  4717.03125  4838.89062 4533.3651
7  x10am 4160.2787  4437.68750 4329.1290  4521.98438  4846.46875 4530.0625
8  x11am 4287.0000  4690.03125 5048.9375  4371.90625  4840.10938 4927.3333
9  x12am 4640.6984  5109.92188 4876.4062  4920.32812  5084.67188 6587.9841
10  x1pm 4708.8095  4707.34375 4544.4688  4600.84375  4776.64062 6714.5714
11  x2pm 4608.1875  4860.25397 4699.6719  4599.54688  5172.25000 6260.8413
12  x3pm 5129.5938  5439.65625 5213.3281  5242.62903  5592.93750 6233.7812
13  x4pm 6366.7969  7221.64062 6809.9844  7197.23810  7489.87500 5574.5000
14  x5pm 9358.8281 10618.87500 9826.9531 10923.03125 10043.54688 5264.7031
15  x6pm 7302.5625  8245.76562 7746.7812  8073.79688  8016.34375 4778.4375
16  x7pm 4306.2656  4874.95312 4657.0156  5044.10938  5002.65625 3993.0000
17  x8pm 3265.5781  3867.25000 4162.9531  4135.29688  3867.45312 3241.6406
18  x9pm 2423.6094  3066.98438 3010.0000  3122.78125  3064.93750 2769.2656
19 x10pm 1110.0469  1431.03175 1542.4375  1616.65625  1828.82812 2124.9688
20 x11pm   23.5283    43.59016  441.6066    85.69841    80.16129  182.6452
21 Total 4791.5282  5388.54144 5168.1490  5287.03210  5379.77892 4558.1144
     domingo     Total
1   576.4219 2147.9312
2  1700.1094 6057.8234
3  2107.9688 7544.5351
4  2100.5312 7064.8345
5  2302.4844 4850.5013
6  2665.0000 4384.7310
7  3023.0156 4264.0895
8  3009.4531 4453.5387
9  3273.5938 4927.6578
10 3614.7969 4809.6392
11 3605.7344 4829.4980
12 3293.6508 5163.6538
13 3356.7344 6288.1099
14 3448.0781 8497.7165
15 3630.2031 6827.6987
16 3355.9219 4461.9888
17 2622.3438 3594.6451
18 2090.6250 2792.6004
19  362.6562 1430.9465
20   58.7500  130.8543
21 2509.9036 4726.1497

register("r1PB_ajs")

r2PB_ajs

result <- pivotr(
  r2.promB, 
  cvars = c("dsem", "hora"), 
  nvar = "n_pasajerosXhora", 
  nr = Inf
)
# summary()
#dtab(result) %>% render()
r2PB_ajs <- result$tab
print(r2PB_ajs)

    hora       lunes      martes   miércoles      jueves    viernes     sábado
1   x4am 1870.630435 2155.387097 2171.419355 2170.173913 2028.76087 1684.75281
2   x5am 4260.086957 5009.333333 5030.688172 5074.478261 4773.87097 3515.16667
3   x6am 4770.456522 5581.107527 5615.440860 5623.478261 5280.15054 3671.53333
4   x7am 3836.891304 4541.311828 4545.053763 4542.804348 4169.98925 3217.87778
5   x8am 2604.217391 3043.698925 3013.075269 3030.000000 2789.55914 2410.70000
6   x9am 2239.228261 2516.354839 2488.591398 2502.945652 2316.43011 2162.20000
7  x10am 2074.793478 2258.236559 2247.903226 2253.695652 2136.68817 2167.00000
8  x11am 2066.728261 2231.258065 2224.268817 2222.902174 2134.83871 2327.14444
9  x12am 2246.869565 2465.505376 2465.075269 2468.521739 2373.79570 3298.23333
10  x1pm 2134.673913 2291.333333 2288.268817 2317.032609 2249.43011 3576.08889
11  x2pm 2447.750000 2570.956989 2580.967742 2593.597826 2572.96774 3427.10000
12  x3pm 2714.913043 2912.849462 2956.462366 2990.826087 2979.35484 3089.00000
13  x4pm 3831.326087 4299.677419 4329.118280 4306.119565 4180.04301 3043.93333
14  x5pm 5120.250000 5896.344086 5852.602151 5877.913043 5362.27957 2735.51111
15  x6pm 3785.750000 4355.892473 4332.774194 4406.543478 4095.73118 2493.25843
16  x7pm 2138.391304 2395.268817 2422.075269 2483.010870 2445.15054 1901.58427
17  x8pm 1499.532609 1672.462366 1679.731183 1713.282609 1698.68817 1511.41573
18  x9pm 1048.141304 1198.526882 1268.728261 1232.413043 1236.44565 1177.23596
19 x10pm  430.097826  545.376344  637.717391  665.922222  654.67033  738.69318
20 x11pm    7.269231    7.241935    7.109375    7.344262   12.83333   23.68421
21 Total 2556.399875 2897.406183 2907.853558 2924.150281 2774.58390 2408.60567
      domingo       Total
1   455.29670 1790.917312
2  1175.26087 4119.840747
3  1284.35870 4546.646534
4  1132.97826 3712.415219
5  1108.27174 2571.360352
6  1392.03261 2231.111838
7  1475.77174 2087.726975
8  1520.01087 2103.878763
9  1618.65217 2419.521879
10 1703.08696 2365.702089
11 1837.93478 2575.896440
12 1724.48913 2766.842133
13 1842.29348 3690.358739
14 1854.09783 4671.285398
15 1777.54348 3606.784748
16 1365.45652 2164.419656
17 1195.45652 1567.224170
18  821.84783 1140.476989
19   84.31522  536.684645
20    1.00000    9.497478
21 1268.50777 2533.929605

register("r2PB_ajs")

Exploración de los datos:

Número de pasajeros promedio por Hora en la linea A (Antes del 23-03-2020)

visualize(
  r1.promA, 
  xvar = "hora", 
  yvar = "n_pasajerosXhora", 
  type = "bar", 
  facet_col = "dsem", 
  axes = "flip", 
  theme = "theme_bw", 
  base_family = "sans", 
  labs = list(
    title = "# De pasajeros promedio\npor Hora en la linea A", 
    subtitle = "Antes del 23-03-2020", 
    caption = "Facetado por días de la semana", 
    x = "Hora", y = "# de Pasajeros promedio"
  ), 
  custom = FALSE
)

Número de pasajeros promedio por Hora en la linea A (después del 23-03-2020)

visualize(
  r2.promA, 
  xvar = "hora", 
  yvar = "n_pasajerosXhora", 
  type = "bar", 
  facet_col = "dsem", 
  axes = "flip", 
  theme = "theme_bw", 
  base_family = "sans", 
  labs = list(
    title = "# de Pasajeros promedios \npor hora en la linea A", 
    subtitle = "Después del 23-03-2020", 
    caption = "Facetado por días de la semana", 
    x = "Hora", y = "# de Pasajeros promedio"
  ), 
  custom = FALSE
)

Número de pasajeros promedio por Hora en la linea B (Antes del 23-03-2020)

visualize(
  r1.promB, 
  xvar = "hora", 
  yvar = "n_pasajerosXhora", 
  type = "bar", 
  facet_col = "dsem", 
  axes = "flip", 
  theme = "theme_bw", 
  base_family = "sans", 
  labs = list(
    title = "# de Pasajeros Promedios\npor hora en la linea B", 
    subtitle = "Antes del 23-03-2020", 
    caption = "Facetado por días de la semana", 
    x = "Hora", y = "# de Pasajeros Promedio"
  ), 
  custom = FALSE
)

Número de pasajeros promedio por Hora en la linea B (Después del 23-03-2020)

visualize(
  r2.promB, 
  xvar = "hora", 
  yvar = "n_pasajerosXhora", 
  type = "bar", 
  facet_col = "dsem", 
  axes = "flip", 
  theme = "theme_bw", 
  base_family = "sans", 
  labs = list(
    title = "# de Pasajeros Promedios\npor hora en la linea B", 
    subtitle = "Después del 23-03-2020", 
    caption = "Facetado por días de la semana", 
    x = "Hora", y = "# de Pasajeros Promedio"
  ), 
  custom = FALSE
)

Plot Line

r1.promA

visualize(
  r1.promA, 
  xvar = "hora", 
  yvar = "n_pasajerosXhora", 
  type = "line", 
  color = "dsem", 
  axes = "flip", 
  theme = "theme_bw", 
  base_family = "sans", 
  labs = list(
    title = "# de pasajeros Promedio\npor hora en la linea A", 
    subtitle = "Antes del 23-03-2020", 
    x = "Hora", y = "# de Pasajeros promedio"
  ), 
  custom = FALSE
)

r2.promA

visualize(
  r2.promA, 
  xvar = "hora", 
  yvar = "n_pasajerosXhora", 
  type = "line", 
  color = "dsem", 
  axes = "flip", 
  theme = "theme_bw", 
  base_family = "sans", 
  labs = list(
    title = "# de pasajeros Promedio\npor hora en la linea A", 
    subtitle = "Después del 23-03-2020", 
    x = "Hora", y = "# de Pasajeros promedio"
  ), 
  custom = FALSE
)

r1.promB

visualize(
  r1.promB, 
  xvar = "hora", 
  yvar = "n_pasajerosXhora", 
  type = "line", 
  color = "dsem", 
  axes = "flip", 
  theme = "theme_bw", 
  base_family = "sans", 
  labs = list(
    title = "# de pasajeros Promedio\npor hora en la linea B", 
    subtitle = "Antes del 23-03-2020", 
    x = "Hora", y = "# de Pasajeros promedio"
  ), 
  custom = FALSE
)

r2.promB

visualize(
  r2.promB, 
  xvar = "hora", 
  yvar = "n_pasajerosXhora", 
  type = "line", 
  color = "dsem", 
  axes = "flip", 
  theme = "theme_bw", 
  base_family = "sans", 
  labs = list(
    title = "# de pasajeros Promedio\npor hora en la linea B", 
    subtitle = "Después del 23-03-2020", 
    x = "Hora", y = "# de Pasajeros promedio"
  ), 
  custom = FALSE
)

Literal F:

Nuevo data. frame

# Data frame temporal
df <- datos_juntos %>% dplyr::select(c(1,2,23))

df$fecha <- df$dia
df$dia <- df$fecha %>% day() #%>% as.factor()
df$dsem <- df$fecha %>% weekdays() %>% factor(., levels = or.dsem)
df$sem <- df$fecha %>% week()
df$mes <- df$fecha %>% month() %>% factor(., labels = meses)
df$mesn <- df$fecha %>% month()
df$year <- df$fecha %>% year() %>% factor()

df %<>% dplyr::select(c("linea","fecha", "dia", "dsem", "sem",
                       "mes", "mesn" ,"year", "total_n_pasajeros"))


names(df)[9] <- "total_pasajerosXdia"


# Filtrado por las lineas A y B:

dfA <- df %>% filter(.,linea=="linea_a") %>% dplyr::select(c(-1))
dfB <- df %>% filter(.,linea=="linea_b") %>% dplyr::select(c(-1))

Literal G:

# Series de Tiempo de la Linea A y B:

tsA <- ts(dfA$total_pasajerosXdia,
          start = c(2019, as.numeric(format(dfA$fecha[1], "%j"))),
          frequency = 365)                                                                       

tsB <- ts(dfB$total_pasajerosXdia,
          start = c(2019, as.numeric(format(dfB$fecha[1], "%j"))),
          frequency = 365)

Graficas Series de Tiempo

plot(tsA, main="# de Pasajeros X día en la Linea A",
     xlab="t en días", ylab="# de Pasajeros totales X día")

Se puede ver un corte abrupto a mediados del 2020 por lo que podemos decir que es el cierre repentino debido a la pandemia lo que baja dramáticamente la afluencia llegando a valores aproximados a cero, dado que solo aquellos con permisos especiales estaban habilitados a viajar en el metro.

además de ver después del corte del 2020 un inicio con expansión multiplicativa, es decir al pasar los días se puede ver un incremento en los pasajeros, dado que poco a poco se iban estableciendo estas rutas para un crecimiento en la afluencia en el metro.

vemos una tendencia alcista por lo que se puede decir que poco a poco el metro iba ganando la cantidad de pasajeros anteriores al corte del 2020, recuperándose de manera significativa (menos de 2 años).

No vemos una serie estacional, bueno la primera parte si representa una vista clara de la serie con algo de estacionalidad, pero después del corte se ve una tendencia y crecimiento de una varianza multiplicativa, es decir, que al paso del tiempo iba aumentando la cantidad de pasajeros de forma de cono (aumentando la varianza)

par(mfrow=c(1,1))
plot(tsB, main="# de Pasajeros X día en la Linea B",
     xlab="t en días", ylab="# de Pasajeros totales X día")

Para la Linea B, es un comportamiento similar…

Solo se diferencian en que la linea B, al no ser una linea de mucho uso en comparación con la linea A, tocaron valores más bajos pero con comportamientos no estacionales.

La serie demuestra estacionalidad en la primera parte, pero para la segunda solo se ve una tendencia con menos pendiente que la linea A, esto representa que la recuperación en la afluencia del metro en la linea B es más demorada dado que es una linea de menor uso en comparación con la A.

Para probar lo anterior veamos sus ACF

par(mfrow=c(1,2))
acf(tsA, plot = F) %>% plot(.,main="ACF de la Linea A")
acf(tsB, plot = F) %>% plot(.,main="ACF de la Linea B")

Como podemos ver las ACF de cada serie de tiempo, demuestra que no son estacionales por lo que se recomendaría dividir los datos.

En este ultimo gráfico veremos los promedios por días de la semana, por lo cual decimos que en el año 2020 fue el más bajo en cuestión de número de pasajeros promedio por día, esto también nos dice que aun en el 2021 no se han recuperado de las cifras, comparadas con el 2019, por lo que esperaríamos que en el 2022 lleguen a un equilibrio.

library(gridExtra)

Attaching package: 'gridExtra'

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

    combine

library(forecast)
# Pasajeros promedio por día de la semana (Linea A)

p1 <- aggregate(dfA,total_pasajerosXdia~dsem+year, mean)[,3]  %>%
  ts(., start = c(2019,1), frequency = 7) %>% ggseasonplot(.,polar = T, size =2)+
  labs(title = "Linea A")+
  xlab("Día")+theme_bw()

# Pasajeros promedio por día de la semana (Linea B)

p2 <- aggregate(dfB,total_pasajerosXdia~dsem+year, mean)[,3]  %>%
  ts(., start = c(2019,1), frequency = 7) %>% ggseasonplot(.,polar = T, size=2)+
  labs(title = "Linea B")+
  xlab("Día")+theme_bw()

grid.arrange(p1,p2,nrow=1, top="# de Pasajeros promedio por días de la semana")

Los sábados el metro no recibe casi pasajeros en promedio con los demás días de la semana.

Literal H:

división de los datos.

en un antes y después (23 de Marzo del 2020)

# Antes y después del 23 de marzo (Linea A)
varH <- ifelse(dfA$fecha<="2020-03-23","before", "after")

dfA$varH <- varH %>% factor(.,levels = c("before", "after"))

# Antes y despues del 23 de marzo (Linea B)
varH <- ifelse(dfB$fecha<="2020-03-23","before", "after")

dfB$varH <- varH

################################### LINEA A ###################################

before <-  dfA %>% filter(varH=="before")

after <- dfA %>% filter(varH=="after")

win.graph()
par(mfrow=c(2,1))
plot(x=before$fecha, y=before$total_pasajerosXdia, type = "l",
     xlab = "t en días", ylab = "# total de Pasajeros",
     main= "Serie Xt antes del 2020-03-23")
abline(h=mean(before$total_pasajerosXdia), col="red", lty=2)
abline(reg=lm(total_pasajerosXdia~fecha, before), col ="blue")
legend("topleft",legend = c("Tendencia con lm()", "Media"),
       col =c("blue", "red"), lty =1:2, cex=0.8)

plot(x=after$fecha, y=after$total_pasajerosXdia, type = "l",
     xlab = "t en días", ylab = "# total de Pasajeros",
     main= "Serie Xt después del 2020-03-23")
abline(h=mean(after$total_pasajerosXdia), col="red", lty=2)
abline(reg=lm(total_pasajerosXdia~fecha, after), col="blue")

legend("topleft",legend = c("Tendencia con lm()", "Media"),
       col =c("blue", "red"), lty =1:2, cex=0.8)

par(mfrow=c(1,1))
p1 <- before %>% ggplot(aes(x=fecha, y=total_pasajerosXdia))+
  geom_line()+
  geom_smooth(method = 'lm',se =F,
                           aes(colour="Tendencia con lm()"))+
  scale_colour_manual(name="", values=c("blue"))+
  geom_hline(aes(yintercept=mean(before$total_pasajerosXdia),
                 linetype="media del proceso"), color = "red")+
   scale_linetype_manual(name = "Linea", values = c(2,1),
                         guide = guide_legend(override.aes=
                                                list(color =c("red"))))+
  labs(title = "Serie Xt antes del 2020-03-23")+
  xlab("t en días")+ylab("# total de pasajeros")+theme_bw()


p2 <- after %>% ggplot(aes(x=fecha, y=total_pasajerosXdia))+
  geom_line()+
  geom_smooth(method = 'lm',se =F,
                           aes(colour="Tendencia con lm()"))+
  scale_colour_manual(name="", values=c("blue"))+
  geom_hline(aes(yintercept=mean(after$total_pasajerosXdia),
                 linetype="media del proceso"), color = "red")+
   scale_linetype_manual(name = "Linea", values = c(2,1),
                         guide = guide_legend(override.aes=
                                                list(color =c("red"))))+
  labs(title = "Serie Xt después del 2020-03-23")+
  xlab("t en días")+ylab("# total de pasajeros")+
  theme_bw()


grid.arrange(p1,p2,nrow=2,top="Linea A")

Warning: Use of `before$total_pasajerosXdia` is discouraged. Use
`total_pasajerosXdia` instead.

`geom_smooth()` using formula 'y ~ x'

Warning: Use of `after$total_pasajerosXdia` is discouraged. Use
`total_pasajerosXdia` instead.

`geom_smooth()` using formula 'y ~ x'

################################### LINEA B ###################################

before <-  dfB %>% filter(varH=="before")

after <- dfB %>% filter(varH=="after")


par(mfrow=c(2,1))
plot(x=before$fecha, y=before$total_pasajerosXdia, type = "l",
     xlab = "t en días", ylab = "# total de Pasajeros",
     main= "Serie Xt antes del 2020-03-23")
plot(x=after$fecha, y=after$total_pasajerosXdia, type = "l",
     xlab = "t en días", ylab = "# total de Pasajeros",
     main= "Serie Xt después del 2020-03-23")

par(mfrow=c(1,1))
p1 <- before %>% ggplot(aes(x=fecha, y=total_pasajerosXdia))+
  geom_line()+
  geom_smooth(method = 'lm',se =F,
                           aes(colour="Tendencia con lm()"))+
  scale_colour_manual(name="", values=c("blue"))+
  geom_hline(aes(yintercept=mean(before$total_pasajerosXdia),
                 linetype="media del proceso"), color = "red")+
   scale_linetype_manual(name = "Linea", values = c(2,1),
                         guide = guide_legend(override.aes=
                                                list(color =c("red"))))+
  labs(title = "Serie Xt antes del 2020-03-23")+
  xlab("t en días")+ylab("# total de pasajeros")+
  theme_bw()


p2 <- after %>% ggplot(aes(x=fecha, y=total_pasajerosXdia))+
  geom_line()+
  geom_smooth(method = 'lm',se =F,
                           aes(colour="Tendencia con lm()"))+
  scale_colour_manual(name="", values=c("blue"))+
  geom_hline(aes(yintercept=mean(after$total_pasajerosXdia),
                 linetype="media del proceso"), color = "red")+
   scale_linetype_manual(name = "Linea", values = c(2,1),
                         guide = guide_legend(override.aes=
                                                list(color =c("red"))))+
  labs(title = "Serie Xt después del 2020-03-23")+
  xlab("t en días")+ylab("# total de pasajeros")+
  theme_bw()


grid.arrange(p1,p2,nrow=2,top="Linea B")

Warning: Use of `before$total_pasajerosXdia` is discouraged. Use
`total_pasajerosXdia` instead.

`geom_smooth()` using formula 'y ~ x'

Warning: Use of `after$total_pasajerosXdia` is discouraged. Use
`total_pasajerosXdia` instead.

`geom_smooth()` using formula 'y ~ x'

Literal I:

dfA$t <- 1:nrow(dfA)
dfB$t <- 1:nrow(dfA)

#modeloa.A <- dfA %>% filter(.,varH=="after")  %>%
#  lm(total_pasajerosXdia~t+dsem+mes+year,.)

#dfA %>% filter(.,varH=="after")  %>%
#lm(total_pasajerosXdia~t+dsem+mes,.) %>% summary()


#summary(modeloa.A)