Datos

Importación de los archivos voitures.xlsx y gasolina.xlsx

library(readxl)

voitures <- read_excel("voitures.xlsx")
gazole <- read_excel("gasolina.xlsx")


# Colores
primary <- "#7a7fff"
secondary <- "darkred"

gazole
# A tibble: 12 × 2
     año sp95_ttc
   <dbl>    <dbl>
 1  2011     1.5 
 2  2012     1.57
 3  2013     1.54
 4  2014     1.48
 5  2015     1.35
 6  2016     1.3 
 7  2017     1.38
 8  2018     1.5 
 9  2019     1.51
10  2020     1.36
11  2021     1.55
12  2022     1.81

Desagregación de la variable gasolina

library(tempdisagg)
gaz <- td(gazole$sp95_ttc ~ 1, to = 12, method = "denton-cholette", conversion="average")

gasolina_disagregated <- predict(gaz)
gasolina_disagregated
  [1] 1.479923 1.480765 1.482450 1.484977 1.488347 1.492559 1.497613 1.503510
  [9] 1.510249 1.517831 1.526255 1.535521 1.545630 1.554428 1.561917 1.568095
 [17] 1.572962 1.576520 1.578767 1.579704 1.579330 1.577647 1.574653 1.570349
 [25] 1.564734 1.559378 1.554279 1.549438 1.544855 1.540530 1.536462 1.532653
 [33] 1.529101 1.525807 1.522771 1.519993 1.517472 1.513961 1.509460 1.503967
 [41] 1.497485 1.490011 1.481547 1.472092 1.461646 1.450210 1.437783 1.424366
 [49] 1.409958 1.396359 1.383569 1.371589 1.360418 1.350056 1.340503 1.331760
 [57] 1.323826 1.316701 1.310385 1.304879 1.300182 1.296561 1.294016 1.292548
 [65] 1.292156 1.292841 1.294602 1.297439 1.301353 1.306342 1.312409 1.319551
 [73] 1.327770 1.336283 1.345092 1.354195 1.363593 1.373286 1.383274 1.393557
 [81] 1.404134 1.415006 1.426173 1.437635 1.449392 1.460559 1.471136 1.481124
 [89] 1.490522 1.499330 1.507549 1.515178 1.522217 1.528667 1.534527 1.539797
 [97] 1.544478 1.546632 1.546259 1.543359 1.537933 1.529979 1.519499 1.506493
[105] 1.490959 1.472899 1.452312 1.429198 1.403558 1.382007 1.364545 1.351173
[113] 1.341890 1.336697 1.335593 1.338579 1.345654 1.356818 1.372072 1.391416
[121] 1.414848 1.438544 1.462504 1.486726 1.511211 1.535959 1.560971 1.586246
[129] 1.611783 1.637584 1.663648 1.689975 1.716565 1.740940 1.763098 1.783041
[137] 1.800767 1.816278 1.829573 1.840653 1.849516 1.856163 1.860595 1.862811
voitures
# A tibble: 160 × 11
   mois    gazole essence hybride_gazole hybride_essence gazole_et_hbride
   <chr>    <dbl>   <dbl>          <dbl>           <dbl>            <dbl>
 1 2011_01 133299   50320             12             896           133311
 2 2011_02 144175   60340             25             993           144200
 3 2011_03 176451   78082             84            1322           176535
 4 2011_04 125996   43699             12            1241           126008
 5 2011_05 147561   48360             45             991           147606
 6 2011_06 152527   52262             67             732           152594
 7 2011_07 113933   44641             59             831           113992
 8 2011_08  75680   31809            166             787            75846
 9 2011_09 122287   46656             81            1273           122368
10 2011_10 134174   44120             70            1319           134244
# ℹ 150 more rows
# ℹ 5 more variables: essence_et_hybride <dbl>, hybride_rechargeable <dbl>,
#   electrique <dbl>, gaz <dbl>, Total <dbl>

Only select the first 144 rows of the voitures dataset

Dataframe

library(flextable)
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
voitures <- voitures[1:144,]

par(family = "Source Serif 4") 

set_flextable_defaults(
  font.size = 12, font.family = "Source Serif 4",
  font.color = "#333333",
  table.layout = "fixed",
  border.color = "gray",
  padding.top = 3, padding.bottom = 3,
  padding.left = 4, padding.right = 4,
  table.width= 0.3
  )

cantidad_autos = voitures$electrique
df <- data.frame(
  mes = voitures$mois, 
  cantidad_autos = (voitures$electrique + voitures$hybride_rechargeable), 
  gasolina= gasolina_disagregated
  )

ft <- flextable(df[1:10,]) %>%
  set_caption(caption = "Datos de autos eléctricos") %>%
  set_header_labels(ft, values = c("Mes", "Cantidad de autos eléctricos e híbridos recargables", "Gasolina Eur")) %>%
  color(j = 1, color = primary, part = "body") 

save_as_image(ft, path = "tabla_datos.png")
[1] "tabla_datos.png"
ft

Mes

Cantidad de autos eléctricos e híbridos recargables

Gasolina Eur

2011_01

103

1.479923

2011_02

140

1.480765

2011_03

251

1.482450

2011_04

192

1.484977

2011_05

145

1.488347

2011_06

143

1.492559

2011_07

67

1.497613

2011_08

109

1.503510

2011_09

300

1.510249

2011_10

208

1.517831

df$cantidad_autos
  [1]   103   140   251   192   145   143    67   109   300   208   491   544
 [13]   332   415   518   438   590   135   616  1059   613   869   416   337
 [25]   357   733  1418   975   719   948   577   441   739  1236   798   814
 [37]   501   483   942   955  1025  1223   872   452  2002  1247  1161  1936
 [49]   852  1409  1588  2195  1582  3042  1766  1031  2269  2389  2236  3001
 [61]  1899  2572  3782  2817  2394  3010  1600  1336  2694  2291  1829  3628
 [73]  2723  2870  3884  2148  2917  4170  2508  2012  3765  3019  3102  4396
 [85]  2367  3208  5480  2912  3186  4738  2963  2979  3896  3971  5052  5724
 [97]  4290  4233  6192  4611  4132  5944  4347  3623  5635  6537  5606  7266
[109] 14774 13456  7383  1665  7229 21194 17091 11035 18044 20482 18935 36383
[121] 14924 17745 29949 21094 24627 37618 18486 17641 29099 27516 28809 39361
[133] 18495 23488 32133 23360 26983 34365 20563 19519 34473 28332 33036 40053
tabla_resumen <- summarizor(df[,2:3]) %>% as_flextable()  %>%
  set_caption(caption = "Resumen de los datos") %>%
  set_header_labels(ft, values = c("Variable", NA, NA, "Estadístico N(144)")) %>%
  color(j = 1, color = primary, part = "body") 
  #set_table_properties(width = 0.8, layout = "autofit")

save_as_image(tabla_resumen, path = "tabla_resumen.png")
[1] "tabla_resumen.png"
tabla_resumen

Variable

Estadístico N(144)

cantidad_autos

Mean (SD)

7,364.9 (10,315.2)

Median (IQR)

2,770.0 (5,785.5)

Range

67.0 - 40,053.0

gasolina

Mean (SD)

1.5 (0.1)

Median (IQR)

1.5 (0.2)

Range

1.3 - 1.9

And save the plots as png images ## Plots

library(tseries)
Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 
cantidad<-ts(df$cantidad_autos, start=c(2011,1), frequency=12)
gasolina1<-ts(df$gasolina, start=c(2011,1), frequency=12)



# png("gasolina.png", width = 800, height = 600, bg = "transparent")  # Ajusta el ancho y la altura según tus preferencias
# Save the plot

par(family = "Source Serif 4") 
plot(gasolina1, main="Gasolina", xlab="Año", ylab="Gasolina Eur (log)", col=primary,col.main=primary, lwd=3,)

# Abre el dispositivo de salida PNG

# Código para generar el gráfico
# dev.off()  # Cierra el dispositivo de salida

# show plot
#png("noprop.png", width = 800, height = 600, bg = "transparent")  # Ajusta el ancho y la altura según tus preferencias
par(family = "Source Serif 4") 
plot(cantidad, main = "Cantidad de autos eléctricos e híbridos",col=primary, col.main=primary, lwd=3, xlab="Año",  ylab = "Cantidad de autos eléctricos")

# Código para generar el gráfico

#dev.off() 

Distribución

grafico de distribución gasolina y de cantidad, y que tenga la linea de distribución, y gardar graficas como imagenes png

#png("dist_gasolina.png", width = 800, height = 600, bg = "transparent")  # Ajusta el ancho y la altura según tus preferencias
par(family = "Source Serif 4")
hist(log(gasolina1), main="Distribución de la gasolina", xlab="Gasolina Eur (log)", col=primary, col.main=primary, lwd=3)
lines(density(log(gasolina1)), col=secondary, lwd=3)

#dev.off()  # Cierra el dispositivo de salida
#png("dist_cantidad.png", width = 800, height = 600, bg = "transparent") # Ajusta el ancho y la altura según tus preferencias
par(family = "Source Serif 4")
hist(log(cantidad), main="Distribución de la cantidad de autos eléctricos e híbridos vendidos al mes", xlab="Cantidad (log)", col=primary, col.main=primary, lwd=3)
lines(density(log(cantidad)), col=secondary, lwd=3)

#dev.off()  # Cierra el dispositivo de salida

Test ADF

library(stargazer)

Please cite as: 
 Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
 R package version 5.2.3. https://CRAN.R-project.org/package=stargazer 
library(flextable)
library(dplyr)
adf_gasolina <- adf.test(gasolina1)
adf_cantidad <- adf.test(cantidad)


adf_summary <- data.frame(
  Data = c("Gasolina", "Cantidad de autos eléctricos e híbridos"),
  Statistic = c(adf_gasolina$statistic, adf_cantidad$statistic),
  P.Value = c(adf_gasolina$p.value, adf_cantidad$p.value),
  Lags = c(adf_gasolina$parameter, adf_cantidad$parameter)
)



# Usar stargazer para mostrar los resultados
ft <- flextable(adf_summary) %>%
  set_caption(caption = "Resultados del Test de Dickey-Fuller Aumentado") %>%
  theme_vanilla() %>%
  set_table_properties(width = 0.8, layout = "autofit") %>%
  align(align = "center", part = "all") %>%
  color(j = 1, color = primary, part = "body") %>%
  color(j = 2, color = secondary, part = "body") %>%
  set_table_properties(width = 0.4, layout = "autofit") %>%
  set_header_labels(ft, values = c(NA, "Estadístico", "Valor p", "Lags"))

# save ft as png
save_as_image(ft, path ="test_Dickey_Fuller.png")
[1] "test_Dickey_Fuller.png"
ft

Estadístico

Valor p

Lags

Gasolina

-1.2372527

0.8940102

5

Cantidad de autos eléctricos e híbridos

0.1256636

0.9900000

5

Modelo

library(dynlm)
Loading required package: zoo

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

    as.Date, as.Date.numeric
?dynlm
starting httpd help server ...
 done
modelo<-dynlm(log10(cantidad)~L(log10(gasolina1))+ log10(gasolina1) )

tabla_regresion <- modelo %>% as_flextable() %>%
  color(j = 1, color = primary, part = "body") %>%
  set_table_properties(width = 0.4, layout = "autofit") %>%
  set_caption(caption = "Modelo de regresión dinámica")  %>%
  set_header_labels(values = c(NA, "Estadístico", "Error estándar", "Valor t", "Valor p", NA))

save_as_image(tabla_regresion, path = "modelo_regresion.png")
[1] "modelo_regresion.png"
tabla_regresion

Estadístico

Error estándar

Valor t

Valor p

(Intercept)

3.130

0.243

12.858

0.0000

***

L(log10(gasolina1))

-67.294

16.016

-4.202

0.0000

***

log10(gasolina1)

68.817

15.631

4.402

0.0000

***

Signif. codes: 0 <= '***' < 0.001 < '**' < 0.01 < '*' < 0.05

Residual standard error: 0.6144 on 140 degrees of freedom

Multiple R-squared: 0.1464, Adjusted R-squared: 0.1342

F-statistic: 12 on 140 and 2 DF, p-value: 0.0000

Test del modelo

library(lmtest)
goldfeld <- gqtest(modelo)

tabla_goldfeld <- goldfeld %>% as_flextable() %>%
  color(j = c(1,2), color = primary, part = "body") %>%
  set_table_properties(width = 0.8, layout = "autofit") %>%
  set_caption(caption = "Test de Goldfeld-Quandt") %>%
  set_header_labels(ft, values = c("df1", "df2", "Estadístico", "Valor p", "Metodología", "Alternativa"))
Multiple parameters; naming those columns df1, df2
save_as_image(tabla_goldfeld, path = "test_goldfeld.png")
[1] "test_goldfeld.png"
tabla_goldfeld

df1

df2

Estadístico

Valor p

Metodología

Alternativa

69.0

68.0

1.6

0.0302 *

Goldfeld-Quandt test

variance increases from segment 1 to 2

Signif. codes: 0 <= '***' < 0.001 < '**' < 0.01 < '*' < 0.05

Jarque Bera test

library(tseries)
library(flextable)
library(dplyr)
jarque <- jarque.bera.test(residuals(modelo)) %>% 
  as_flextable(jarque) %>%
  color(j = 1, color = primary, part = "body") %>%
  set_caption(caption = "Test de Jarque-Bera") %>%
  set_header_labels(values = c("Estadístico", "Valor p", "Parámetro", "Método")) 

save_as_image(jarque, path = "test_jarque.png")
[1] "test_jarque.png"
jarque

Estadístico

Valor p

Parámetro

Método

4.4

0.1129

2.0

Jarque Bera Test

Signif. codes: 0 <= '***' < 0.001 < '**' < 0.01 < '*' < 0.05

Test Durbin Watson

dw <- dwtest(modelo)

tabla_dw <- as_flextable(dw) %>%
  color(j = 1, color = primary, part = "body") %>%
  set_caption(caption = "Test de Durbin-Watson") %>%
  set_header_labels(values = c("Estadístico", "Valor p", "Metodología", "Alternativa")) %>%
  set_table_properties(width = 0.8, layout = "autofit") 
  
save_as_image(tabla_dw, path = "test_dw.png")
[1] "test_dw.png"
tabla_dw

Estadístico

Valor p

Metodología

Alternativa

0.1

0.0000***

Durbin-Watson test

true autocorrelation is greater than 0

Signif. codes: 0 <= '***' < 0.001 < '**' < 0.01 < '*' < 0.05

gr <- gen_grob(tabla_dw, fit = "fixed", just = "center")
dims <- dim(gr)

svglite::svglite(filename = "hello-grid-graphics.svg", width = dims$width + .1, height = dims$height + .1)
plot(gr)
dev.off()
png 
  2 

Residuos

residuos<-modelo$residuals
residuos<-as.numeric(residuos)

#png("residuos.png", width = 800, height = 600, bg = "transparent")  # Ajusta el ancho y la altura según tus preferencias
par(family = "Source Serif 4") 
plot(lag(residuos),residuos, main="Residuos vs. Lags de los residuos", xlab="Lag", ylab="Residuos", col=primary, lwd=2)
abline(lm(residuos~lag(residuos)), col=secondary, lwd=3)

#dev.off()  # Cierra el dispositivo de salida

Corrección de Cochrane-Orcutt

library(orcutt)

modelofinal<-cochrane.orcutt(modelo)

sumary <- summary.orcutt(modelofinal)

coefficients <- sumary$coefficients

# Extrae los errores estándar
std_errors <- modelo$residuals

sumary$stat
[1]   0.2170000 139.0000000   0.0031000  -0.0112000   0.2000000   2.0000000
[7]   0.8047944
# Extrae los valores de t
t_values <- coef(sumary)[, "t value"]

# Extrae los valores p
p_values <- coef(sumary)[, "Pr(>|t|)"]


# Crea un data frame con los resultados
result <- data.frame(
  Coefficient = rownames(coefficients),
  Estimate = coefficients
 # Std_Error = std_error
)


tabla_modelo_orcutt <- flextable(result) %>%
  set_caption(caption = "Modelo de regresión dinámica con corrección de Cochrane-Orcutt") %>%
  color(j = 1, color = primary, part = "body") %>%
  set_table_properties(width = 0.8, layout = "autofit") %>%
  set_caption(caption = "Modelo de regresión dinámica con corrección de Cochrane-Orcutt") %>%
  set_header_labels(ft, values = c("Coeficiente", "Estimación", "Error estándar", "Valor t", "Valor p"))
                    
save_as_image(tabla_modelo_orcutt, path = "modelo_orcutt.png")
[1] "modelo_orcutt.png"
## R cuadrado signif code
modelo$
tabla_modelo_orcutt
NULL
sumary
Call:
dynlm(formula = log10(cantidad) ~ L(log10(gasolina1)) + log10(gasolina1))

                    Estimate Std. Error t value  Pr(>|t|)    
(Intercept)          3.22336    0.77297   4.170 5.331e-05 ***
L(log10(gasolina1))  2.60289   34.40187   0.076    0.9398    
log10(gasolina1)     0.02853   34.17741   0.001    0.9993    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.217 on 139 degrees of freedom
Multiple R-squared:  0.0031 ,  Adjusted R-squared:  -0.0112
F-statistic: 0.2 on 2 and 139 DF,  p-value: < 8.048e-01

Durbin-Watson statistic 
(original):    0.13338 , p-value: 1.794e-30
(transformed): 2.49309 , p-value: 9.974e-01