ACTIVIDAD 01 SERIES DE TIEMPO

Maestría en Investigación Operativa y Estadística

1. Motivación de la investigación

La transición hacia fuentes de energía más sostenibles representa un desafío crucial tanto para las naciones desarrolladas como para aquellas en vías de desarrollo. En este contexto, es imperativo que los países reconozcan los beneficios asociados con la adopción de fuentes de energía no convencionales (FNCER). El objetivo de este estudio es explorar y analizar la relación existente entre el capital (expresado en millones de dólares) de un estado en los Estados Unidos de América y su propensión al uso de fuentes de energía alternativas y sostenibles.

1.1 Lectura de datos

# Cargue y limpieza de datos.

EMISIONES <- read.xlsx("CAPITAL-FNCER.xlsx") %>% clean_names()

1.2 Función Flextable, que permite mejorar la presentación de la base de datos

# Funciones

## Función para crear flextable
ftable <- function(x) {
  x %>% 
    flextable() %>% 
    theme_vanilla() %>%
    color(part = "footer", color = "#666666") %>%
    color( part = "header", color = "#FFFFFF") %>%
    bg( part = "header", bg = "#2c7fb8") %>%
    fontsize(size = 11) %>%
    font(fontname = 'Calibri') %>%
    # Ajustes de ancho y tipo de alineación de las columnas
    set_table_properties(layout = "autofit") %>% 
    # width(j=1, width = 3) %>%
    align(i = NULL, j = c(2:ncol(x)), align = "right", part = "all")
}

EMISIONES%>% 
  head(5) 

##       estado estado_simbolo capacidad_instalada_mw
## 1     Alaska             AK                 3577.1
## 2    Alabama             AL                39620.5
## 3  Arkansas              AR                19856.8
## 4    Arizona             AZ                41605.9
## 5 California             CA               109745.3
##   consumo_de_energia_primaria_mm_btu generacion_electrica_total_m_wh
## 1                           56817777                         6068243
## 2                         1295527860                       142679433
## 3                          599795550                        63857739
## 4                         1074295701                       113368628
## 5                         1707003763                       201747828
##   emisiones_n_ox_t emisiones_pm2_5_t emisiones_so2_t emisiones_co2e_t
## 1         20241.13          1222.723        2167.876          2958055
## 2         21325.02          2151.446        9889.930         56042549
## 3         18346.78          2493.187       41351.871         36013942
## 4         29340.65          4038.691       12169.717         49503782
## 5         35485.85          2156.422        3400.458         39056658
##   ic_t_co2e_m_wh generacion_electrica_total_no_renovable_m_wh
## 1      0.4874649                                      4264738
## 2      0.3927865                                    127531820
## 3      0.5639715                                     58129774
## 4      0.4366621                                    101099802
## 5      0.1935915                                    104363268
##   generacion_electrica_total_renovable_m_wh      fncer
## 1                                   1803504 0.01012684
## 2                                  15147612 0.08505523
## 3                                   5727965 0.03216305
## 4                                  12268827 0.06889058
## 5                                  97384559 0.54682322
##   generacion_electrica_total_renovable_no_hidraulica_m_wh
## 1                                                180138.1
## 2                                               3742396.2
## 3                                               1542545.1
## 4                                               6061914.5
## 5                                              59060405.7
##   participacion_de_generacion_electrica_renovable_percent capital_musd
## 1                                              0.29720369     139635.7
## 2                                              0.10616535     232079.9
## 3                                              0.08969884     147335.4
## 4                                              0.10822065     139635.7
## 5                                              0.48270437    3347858.5
##      capital gdp_real_musd hdi_subnacional gdp_ghg_usd_t hdi_ghg_mt_co2e_1
## 1 0.02747137       53433.8           0.939     18063.827        0.31743828
## 2 0.04565849      203432.7           0.889      3629.969        0.01586295
## 3 0.02898618      117126.2           0.889      3252.246        0.02468488
## 4 0.02747137      325395.3           0.916      6573.140        0.01850364
## 5 0.65864448     2729225.8           0.939     69878.630        0.02404200
##   mano_de_obra_productiva poblacion_estado_2019 poblacion_estado_2020
## 1                  454785                731545                733391
## 2                 2712163               4903185               5024279
## 3                 1664214               3017804               3011524
## 4                 3929985               7278717               7151502
## 5                24227775              39512223              39538223
##   mano_de_obra_productiva_percent porcentaje_de_desempleo_percent
## 1                       0.6216774                             5.6
## 2                       0.5531431                             3.2
## 3                       0.5514652                             3.5
## 4                       0.5399283                             4.8
## 5                       0.6131717                             4.1
##   pobreza_percent
## 1           0.101
## 2           0.155
## 3           0.162
## 4           0.135
## 5           0.118

Con flextable

EMISIONES %>% 
  head(10) %>% 
  ftable()

estado	estado_simbolo	capacidad_instalada_mw	consumo_de_energia_primaria_mm_btu	generacion_electrica_total_m_wh	emisiones_n_ox_t	emisiones_pm2_5_t	emisiones_so2_t	emisiones_co2e_t	ic_t_co2e_m_wh	generacion_electrica_total_no_renovable_m_wh	generacion_electrica_total_renovable_m_wh	fncer	generacion_electrica_total_renovable_no_hidraulica_m_wh	participacion_de_generacion_electrica_renovable_percent	capital_musd	capital	gdp_real_musd	hdi_subnacional	gdp_ghg_usd_t	hdi_ghg_mt_co2e_1	mano_de_obra_productiva	poblacion_estado_2019	poblacion_estado_2020	mano_de_obra_productiva_percent	porcentaje_de_desempleo_percent	pobreza_percent
Alaska	AK	3,577.1	56,817,777	6,068,243	20,241.128	1,222.72281	2,167.876	2,958,055	0.4874649	4,264,738	1,803,504	0.0101268404	180,138.1	0.29720369	139,635.66	0.02747137	53,433.8	0.939	18,063.827	0.317438285	454,785	731,545	733,391	0.6216774	5.6	0.101
Alabama	AL	39,620.5	1,295,527,860	142,679,433	21,325.024	2,151.44580	9,889.930	56,042,549	0.3927865	127,531,820	15,147,612	0.0850552292	3,742,396.2	0.10616535	232,079.91	0.04565849	203,432.7	0.889	3,629.969	0.015862947	2,712,163	4,903,185	5,024,279	0.5531431	3.2	0.155
Arkansas	AR	19,856.8	599,795,550	63,857,739	18,346.780	2,493.18744	41,351.871	36,013,942	0.5639715	58,129,774	5,727,965	0.0321630486	1,542,545.1	0.08969884	147,335.39	0.02898618	117,126.2	0.889	3,252.246	0.024684885	1,664,214	3,017,804	3,011,524	0.5514652	3.5	0.162
Arizona	AZ	41,605.9	1,074,295,701	113,368,628	29,340.647	4,038.69073	12,169.717	49,503,782	0.4366621	101,099,802	12,268,827	0.0688905841	6,061,914.5	0.10822065	139,635.66	0.02747137	325,395.3	0.916	6,573.140	0.018503637	3,929,985	7,278,717	7,151,502	0.5399283	4.8	0.135
California	CA	109,745.3	1,707,003,763	201,747,828	35,485.852	2,156.42151	3,400.458	39,056,658	0.1935915	104,363,268	97,384,559	0.5468232173	59,060,405.7	0.48270437	3,347,858.50	0.65864448	2,729,225.8	0.939	69,878.630	0.024041995	24,227,775	39,512,223	39,538,223	0.6131717	4.1	0.118
Colorado	CO	22,303.6	545,610,682	56,336,873	20,009.472	1,252.51227	11,414.025	37,494,245	0.6655365	42,485,115	13,851,758	0.0777788880	12,232,638.8	0.24587374	428,078.38	0.08421845	358,438.5	0.951	9,559.827	0.025363893	3,895,377	5,758,736	5,773,714	0.6764292	2.7	0.093
Connecticut	CT	12,814.8	346,263,247	40,050,038	4,690.635	228.78021	426.104	9,560,385	0.2387110	38,141,741	1,908,297	0.0107152647	1,478,225.1	0.04764783	277,801.19	0.05465351	251,568.2	0.958	26,313.606	0.100205172	2,298,065	3,565,287	3,605,944	0.6445666	3.6	0.100
Delaware	DE	4,066.9	32,888,543	5,258,538	449.625	86.39776	305.777	1,870,668	0.3557392	5,138,072	120,466	0.0006764277	120,466.0	0.02290865	57,946.93	0.01140025	64,143.5	0.939	34,289.085	0.501959684	602,812	973,764	989,948	0.6190535	3.6	0.113
Florida	FL	108,054.7	2,001,477,421	244,949,921	40,368.821	7,273.67854	21,423.446	107,444,751	0.4386397	234,827,314	10,122,607	0.0568393662	9,912,725.2	0.04132521	978,101.06	0.19242775	965,672.5	0.919	8,987.619	0.008553233	12,761,585	21,477,737	21,538,187	0.5941774	3.2	0.127
Georgia	GA	53,593.9	1,174,724,543	128,595,117	24,658.740	2,200.40623	18,952.850	56,642,732	0.4404734	117,701,120	10,893,996	0.0611707860	7,536,461.2	0.08471547	646,428.88	0.12717587	557,364.4	0.913	9,839.999	0.016118573	6,399,395	10,617,423	10,711,908	0.6027258	3.6	0.133

1.3 Base de datos filtrada por variables de interés

A <- data.frame(EMISIONES$estado, EMISIONES$estado_s, EMISIONES$fncer, EMISIONES$capital)


 A%>% 
  head(10) %>% 
  ftable()

EMISIONES.estado	EMISIONES.estado_s	EMISIONES.fncer	EMISIONES.capital
Alaska	AK	0.0101268404	0.02747137
Alabama	AL	0.0850552292	0.04565849
Arkansas	AR	0.0321630486	0.02898618
Arizona	AZ	0.0688905841	0.02747137
California	CA	0.5468232173	0.65864448
Colorado	CO	0.0777788880	0.08421845
Connecticut	CT	0.0107152647	0.05465351
Delaware	DE	0.0006764277	0.01140025
Florida	FL	0.0568393662	0.19242775
Georgia	GA	0.0611707860	0.12717587

1.4 Diagrama de dispersión.

EMISIONES %>% 
  ggplot(aes(x= capital, y= fncer )) +
  geom_point()+ theme_light()

Según la información visualizada en el gráfico de dispersión, parece haber una tendencia en los datos. Sin embargo, la presencia de valores atípicos u outliers dificulta la interpretación y comprensión completa de la relación entre las variables.

Se llevará a cabo un proceso de limpieza de datos mediante el siguiente procedimiento:

Los valores que se encuentren por encima de Q3 + 1.5IQR (el tercer cuartil más 1.5 veces el rango intercuartil) y por debajo de Q1 - 1.5IQR (el primer cuartil menos 1.5 veces el rango intercuartil) serán identificados como outliers y eliminados de la muestra. Esta acción tiene como objetivo mejorar la calidad de los datos y permitir una exploración más precisa de la relación entre las variables.

Una vez que se hayan eliminado estos valores atípicos, se podrá realizar un análisis más sólido y revelador de la tendencia subyacente en los datos, lo que nos ayudará a obtener una comprensión más profunda de la relación en cuestión.

1.5 Se eliminan Outliers y se realiza nuevamente el gráfico de dispersión

#OUTLIERS capital    datos_fncer

ggplot(EMISIONES, aes(x = capital, y = fncer, fill = capital,fncer)) +
  geom_boxplot(outlier.colour = "red", outlier.shape = 19, position = position_dodge(width = 0.0001)) +
  geom_jitter(width = 0.0001, shape = 20, size = 3, alpha = 0.8) +  # Agregar los puntos con jitter
  labs(title = "Distribucuón capital-fncer - Boxplot (Outliers) ",
       x = "Capital",
       y = "Fncer") +
  theme_minimal()

# Calcular Q1 (cuantil 1)
q1 <- quantile(EMISIONES$capital, probs = 0.25)

# Calcular Q3 (cuantil 3)
q3 <- quantile(EMISIONES$capital, probs = 0.75)

#Calculo del IQR 

IQR <- q3 -q1

Valores_debajo <- q1 - 1.5*IQR
Valores_arriba <- q3 + 1.5*IQR


datos_filtrados_capital <- subset(EMISIONES$capital, EMISIONES$capital > Valores_debajo & EMISIONES$capital < Valores_arriba)



#OUTLIERS fncer 

datos_fncer <- EMISIONES$fncer

# Calcular Q1 (cuantil 1)
q1 <- quantile(EMISIONES$fncer, probs = 0.25)

# Calcular Q3 (cuantil 3)
q3 <- quantile(EMISIONES$fncer, probs = 0.75)

#Calculo del IQR 

IQR <- q3 -q1

Valores_debajo <- q1 - 1.5*IQR
Valores_arriba <- q3 + 1.5*IQR


datos_filtrados_fncer <- subset(EMISIONES$fncer , EMISIONES$fncer  > Valores_debajo & EMISIONES$fncer< Valores_arriba)

#Se crea un data frame con los datos sin outliers

#Filtrar outliers en datos_capital y datos_fncer

datos_filtrados_capital <- subset(EMISIONES$capital, EMISIONES$capital > Valores_debajo & EMISIONES$capital < Valores_arriba)

datos_filtrados_fncer <- subset(EMISIONES$fncer, EMISIONES$fncer > Valores_debajo & EMISIONES$fncer < Valores_arriba)


data_combinada <- data.frame(capital = EMISIONES$capital, fncer = EMISIONES$fncer)

# Elimina outliers en el dataframe combinado

data_filtrados <- subset(data_combinada, capital > Valores_debajo & capital < Valores_arriba & fncer > Valores_debajo & fncer < Valores_arriba)

e <- length(datos_filtrados_capital)
d <- length(datos_filtrados_fncer)

# Crear el gráfico de dispersión
ggplot(data = data_filtrados, aes(x = capital, y = fncer)) +
  geom_point() +
  labs(title= "Gráfico de dispersión sin outliers", x = "Capital", y = "FNCER") +
  theme_light()

ggplot(data_filtrados, aes(x = capital, y = fncer, fill = capital,fncer)) +
  geom_boxplot(outlier.colour = "red", outlier.shape = 19, position = position_dodge(width = 0.0001)) +
  geom_jitter(width = 0.0001, shape = 20, size = 3, alpha = 0.8) +  
  labs(title = "Distribucuón fncer - Boxplot ( Sin Outliers) ",
       x = "Capital",
       y = "Fncer") +
  theme_minimal()

1.5 Generando un modelo súper simplificado (modelo ingenuo)

ggplot(data = data_filtrados, aes(x = capital, y = fncer)) +
  geom_point() + 
  geom_hline(yintercept = mean(data_filtrados$fncer), color = "red") +
  annotate("text", x = max(data_filtrados$capital), y = mean(data_filtrados$fncer) + 0.001,
           label = paste("Participación en fuentes no convencionales:", round(mean(data_filtrados$fncer), 4)),
           hjust = 1, vjust = -0.2, color = "red") +
  labs(x = "Capital", y = "FNCER") +
  theme_light()

## 1.6 Calidad de ajuste del modelo ingenuo

Se aplica el siguiente método que permita hallar el desajute del modelo ingenuo.

\[MSE_{\text{ingenuo}} = \color{red}{\frac{1}{n}}\sum^{n}_{i=1}(y_i - \bar{y})²\]

Es decir, nuestra medida de desajuste será el promedio de la suma de las diferencias cuadradas de cada dato \(y_i\) respecto a su media \(\bar{y}\), dicho de otra manera la varianza de la variable \(Y\) a predecir!!!!

n <- length(datos_filtrados_fncer)
MSE_ingenuo  <- (1/n) * sum((datos_filtrados_fncer - mean(datos_filtrados_fncer))^2)
MSE_ingenuo

## [1] 0.001063608

A continuaciòn, se presenta el gràfico que ilustra las distancias que hay entre cada dato \(y_i\) observado y el promedio histórico (valor predicho \(\bar{y}\) por el modelo ingenuo)

data_filtrados %>%
  ggplot(aes(x = capital, y = fncer)) +
  geom_point() + 
  geom_hline(yintercept = mean(data_filtrados$fncer), color = "red") +
  geom_segment(aes(xend = capital, yend = mean(data_filtrados$fncer)),
               col = 'red', lty = 'dashed') +
  labs(title = "Distancia entre cada dato observado y el promedio histórico", x = "Capital", y = "FNCER") +
  theme_light()

1.7 Introducción a las medidas de ajuste para modelos no ingenuos

A continuación, se presenta un posible posible ajustado usando valores \(\beta_0=-0.00001\) y \(\beta_1 = 1.0005\). Se eligieron los valores por inspección visual y prueba y error:

b0 <- -0.00001
b1 <- 1.0005


# Predicciones
x <- data_filtrados$capital
y_pred <- b0 + b1 * x

# Gráfica

data_filtrados %>%
  ggplot(aes(x = capital, y = fncer)) +
  geom_point() +
  geom_line(aes(x = x, y = y_pred), col = "red") +
  labs(x = "Capital", y = "FNCER") +
  theme_light()

1.8 Una animación que ilustra la estimación de betas

# Extraemos las variables de interés
x = data_filtrados$capital
y = data_filtrados$fncer

# Calculamos medias muestrales para las variables X y Y

x_barra <- mean(x)
y_barra <- mean(y)

# Definimos posibles valores para b1, y por consiguiente para b0
b1_values <- seq(-1, 0.5, by = 0.08)
b0_values <- y_barra - b1_values * x_barra

# Fijamos límites en X y Y para cada gráfico
xmin = min(x) - 0
xmax = max(x) + 1
ymin = min(y) - 0
ymax = max(y) + 1

# Construimos una función para calcular el MSE en función de x_i, y_i, b1 y b0
MSE <- function(x, y, b0, b1) {
  errors <- y - (b0 + b1*x)
  squared_errors <- errors^2
  mse <- mean(squared_errors)
  return(mse)
}

# Crear un dataframe para almacenar los resultados
list_mse <- vector()

# Calcular el MSE para diferentes valores de B1
for (i in 1:length(b1_values)) {
  mse <- MSE(x, y, b0_values[i], b1_values[i])
  list_mse[i] <- mse
}

mse_df = data.frame(
  B0 = b0_values,
  B1 = b1_values,
  mse = list_mse,
  color_point = 'green'
)

# Capturando el renglón donde ocurre el mínimo mse
pos_min <- which(mse_df$mse==min(mse_df$mse))

# Impriendo las filas cercanas al mínimo
mse_df$color_point[(pos_min - 4):(pos_min + 4)] <- rep('red',9)
ftable(mse_df[(pos_min - 8):(pos_min + 8), ])

B0	B1	mse	color_point
0.07451374	-0.60	0.001643409	green
0.07048342	-0.52	0.001515246	green
0.06645311	-0.44	0.001404177	green
0.06242280	-0.36	0.001310204	green
0.05839248	-0.28	0.001233325	red
0.05436217	-0.20	0.001173541	red
0.05033185	-0.12	0.001130851	red
0.04630154	-0.04	0.001105257	red
0.04227123	0.04	0.001096757	red
0.03824091	0.12	0.001105352	red
0.03421060	0.20	0.001131042	red
0.03018028	0.28	0.001173826	red
0.02614997	0.36	0.001233705	red
0.02211966	0.44	0.001310679	green

Ahora graficaremos una curva que muestre esta relación:

if (!file.exists("beta_vs_mse.gif")){
  # Crear la animación
  animation <- ggplot(mse_df, aes(B1, mse)) +
    geom_line(aes( x= B1, y = mse, color = color_point), linetype = "dashed") +
    geom_point(aes(color = color_point), size = 3) +
    labs(title = "Valores de MSE para cada posible B1, respetando E(e) = 0") +
    geom_text(data = mse_df, 
              aes(label = paste("MSE: ", round(mse,0))), x= 0.5, y = 0.5, color='black'
    ) +
    geom_text(data = mse_df, 
              aes(label = paste("MSE mínimo: ", round(min(mse),0))),
              x= 1, y = min(mse_df$mse) - 1e5, color='red') +
    scale_color_manual(values = c("red" = "red", "green" = "green"),
                       labels = c("Subóptima", "Óptima"),
                       name = "Soluciones") +
    theme_minimal() +
    transition_reveal(B1)

  animate(animation, duration = 10, fps = 2, renderer = gifski_renderer())
  anim_save("beta_vs_mse.gif")
}

Img 1: Relación entre B1 y MSE

if (!file.exists("beta_estimation.gif")){
  
# Las pausas en la animación son fotogramas con información repetida correspondiente
# al renglón óptimo

fotogramas_pausa <- data.frame(
  B0 = rep(b0_values[pos_min], 2*nrow(mse_df)),
  B1 = rep(b1_values[pos_min], 2*nrow(mse_df)),
  mse = rep(list_mse[pos_min], 2*nrow(mse_df)),
  color_point = rep('red',2*nrow(mse_df))  
)

mse_df <- rbind(mse_df, fotogramas_pausa)
mse_df <- mse_df[order(mse_df$B1), ]
mse_df$estado <- seq(1,nrow(mse_df))

# Construimos el dataframe de predicciones para cada estado, con el cual poder
# graficar los residuales en cada fotograma

df_predicciones <- data.frame(B0_pred=numeric(0),
                              B1_pred=numeric(0),
                              x_from_pred=numeric(0),
                              y_from_pred=numeric(0),
                              y_pred_mod=numeric(0),
                              estado = numeric(0))
for (i in 1:nrow(mse_df)){
  df_nuevo <-data.frame(
    BO_pred = rep(mse_df[i, 'B0'], length(x)),
    B1_pred = rep(mse_df[i, 'B1'], length(x)),
    x_from_pred = x,
    y_from_pred = y,
    y_pred_mod = mse_df[i, 'B0'] + mse_df[i,'B1']*x,
    estado = rep(mse_df[i, 'estado'], length(x))
  ) 
  df_predicciones <- rbind(df_predicciones, df_nuevo)
}
w <- data.frame(df_predicciones)
# Creamos el gráfico base

base_plot <- ggplot(data.frame(x = x, y = y), aes(x = x, y = y)) +
  geom_point() +
  geom_point(x=mean(x), y=mean(y), color='green', size=2.5) +
  geom_abline(aes(slope = 0, intercept = y_barra), 
              linetype = "dashed", color = "blue") +
  labs(title = "Estimación de Betas por Minimización del MSE") +
  theme_minimal()

# Creamos la animación

animation <- base_plot +
  geom_abline(aes(slope = B1, intercept = B0, color = color_point), data=mse_df) +
  geom_text(aes(label = paste("MSE: ", round(mse,0))), x= 0.5, y = 0.5, data=mse_df,
            color = "red")+
  geom_segment(data = df_predicciones,
               aes(x=x_from_pred, y=y_from_pred, xend=x_from_pred, yend=y_pred_mod),
               col = "violet", lty='dashed') +
  transition_states(estado, transition_length = 2, state_length = 1) +
  scale_color_manual(values = c("red" = "red", "green" = "green"),
                     labels = c("Subóptima", "Óptima"),
                     name = "Rectas regresoras") +
  enter_fade() +
  exit_fade()

# Guardamos la animación en un archivo GIF
animate(animation, duration = 20, fps = 2, renderer = gifski_renderer())
anim_save("beta_estimation.gif", animation)
}

Img 1: Estimación de betas por minimización de MSE

# b1 = cov(x,y) / (sigma_x)²
b1 <- cov(x,y) / var(x)
print(b1)

## [1] 0.03977735

# bo = y_barra - b1*x_barra
b0 <- mean(y) - b1*mean(x)
print(b0)

## [1] 0.04228244

residuales <- y - (b0 + b1*x)
n <- length(x)

sigma_2 <- sum(residuales^2)/(n-2)
print(sigma_2)

## [1] 0.001153001

# Confirmamos usando la función lm, a continuación la documentación

# lm(formula, data, subset, weights, na.action,
#    method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE,
#    singular.ok = TRUE, contrasts = NULL, offset, ...)

mod1 <- data_filtrados %>% lm(fncer ~ capital,.)
print(mod1$coefficients)

## (Intercept)     capital 
##  0.04228244  0.03977735

resumen <-  summary(mod1)
resumen$sigma

## [1] 0.03395587

parametros <- data.frame(metodo = c("manual","lm"), 
                         beta_0 = c(b0,mod1$coefficients[1]),
                         beta_1 = c(b1,mod1$coefficients[2]),
                         sigma_2 = c(sigma_2, (resumen$sigma)^2))



parametros %>% 
  ftable()

metodo	beta_0	beta_1	sigma_2
manual	0.04228244	0.03977735	0.001153001
lm	0.04228244	0.03977735	0.001153001

De acuerdo con la tabla antes expuesta se comprueba que ambos métodos (manual y librería lm de R studio) permiten llegar a la misma respuesta B0=0.04228244; B1=0.03977735 y sigma_2= 0.001153001

El modelo de regresión lineal es de la forma: Y = b0 + b1X + ε, donde Y es la variable de respuesta (en este caso, generación de energía eléctrica empleando fuentes no convencionales de energía renovable), X es la variable predictora (capital de cada uno de los estados de Estados Unidos), b0 es la intersección (o intercepto), b1 es el coeficiente de la variable predictora, y ε es el error aleatorio.

Los resultados del modelo indican que la intersección (intercepto) es aproximadamente 0.04228244 y el coeficiente asociado a la variable predictora “capital” es aproximadamente 0.03977735. Esto significa que, la relación entre el capital de los estados de Estados Unidos y la generación de energía eléctrica de fuentes renovables es positiva. En otras palabras, a medida que el capital aumenta, la generación de energía eléctrica de fuentes renovables tiende a aumentar. Por cada unidad adicional de capital, se estima que la generación de energía eléctrica de fuentes renovables aumenta en aproximadamente 0.0398 unidades. Lo que indicaría que frente a los retos continuos de disminuir emisiones a la atmosfera producto de la generación de energía eléctrica, se puede replantear la asignación del presupuesto para las generación de energía eléctrica con fuentes no convencionales (FNCER).

2. COMPROBACIÓN DEL SUPUESTO DE NORMALIDAD.

# Paso 1: Ordenamos los datos de menor a mayor
n <- length(mod1$residuals)
df_residuals <- data.frame(residuals = mod1$residuals)
df_residuals <- df_residuals %>% 
                arrange(residuals)

# Paso 2: Estimamos los parámetros para la distribución normal teórica a partir de los datos

mean_residuals <- mean(df_residuals$residuals)
sd_residuals <- sd(df_residuals$residuals)

# Paso 3: Calculamos la función acumulada para cada valor de residual
cdf_teorica <- pnorm(df_residuals$residuals,
                      mean = mean_residuals,
                      sd = sd_residuals)

# Paso 4: Calculamos la máxima diferencia absoluta entre la densidad empírica y la teórica

# Método 1: Construir la función empírica usando fórmulas de R (ecdf) y luego tomar diferencias
#           absolutas entre valor empírico y teórico (enfoque KS)
cdf_empirica <- ecdf(df_residuals$residuals)
diferencias_absolutas <- abs(cdf_empirica(df_residuals$residuals) - cdf_teorica)
pos_estadistico_m1 <- which.max(diferencias_absolutas)
estadistico_m1 <- diferencias_absolutas[pos_estadistico_m1]

print(paste0("Este es el valor del estadístico: ", estadistico_m1," y esta es la posición en ",
             "la que ocurre la máxima diferencia ", pos_estadistico_m1))

## [1] "Este es el valor del estadístico: 0.168225649337558 y esta es la posición en la que ocurre la máxima diferencia 22"

# Método 2: Construir manualmente los cálculos que permiten hallar la máxima diferencia
#           entre la empírica y la teórica usando la definición de la empírica como
#           función a trozos 📈

e_plus <- seq(1:n)/n
e_minus <- seq(0:(n-1))/n

D_plus <- e_plus - cdf_teorica
D_minus <- cdf_teorica - e_minus

pos_max_d_plus <- which.max(D_plus)
pos_max_d_minus <- which.max(D_minus)

estadistico_m2 <- max(D_plus[pos_max_d_plus], D_minus[pos_max_d_minus])
pos_estadistico_m2 <- function() {
  if(max(D_plus)>=max(D_minus)) {
    return(pos_max_d_plus)
    }else{
      return(pos_max_d_minus)
    }
  }

print(paste0("Este es el valor del estadístico: ", estadistico_m2," y esta es la posición en ",
             "la que ocurre la máxima diferencia ", pos_estadistico_m2()))

## [1] "Este es el valor del estadístico: 0.168225649337558 y esta es la posición en la que ocurre la máxima diferencia 22"

# Comprobamos los resultados invocando la función de R que realiza estos cálculos automáticos

Lilliefors_test <- lillie.test(df_residuals$residuals)
print(Lilliefors_test)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  df_residuals$residuals
## D = 0.16823, p-value = 0.005079

estadistico_l <- Lilliefors_test$statistic
print(paste0("Este es el estadístico de Lilliefors usando paquetería de R: ", estadistico_l))

## [1] "Este es el estadístico de Lilliefors usando paquetería de R: 0.168225649337558"

Hasta acá hemos verificado que los cálculos realizados manualmente coinciden con el código implementado por la librería nortest de R, y más específicamente en el código fuente de la función lillie.test. Ahora ilustremos algo de estos pasos usando una gráfica:

La prueba de normalidad de Lilliefors (Kolmogorov-Smirnov) se utiliza para evaluar si los residuales del modelo de regresión se distribuyen de manera normal. Esto es importante porque muchos modelos de regresión asumen que los residuales siguen una distribución normal, sin embargo, el resultado de la prueba muestra que el valor de D (estadístico de prueba de Lilliefors) es aproximadamente 0.16823, y el valor p es aproximadamente 0.005079, como el valor p es menor que un nivel de significancia previamente establecido (como 0.05), se rechaza la hipótesis nula de que los residuales siguen una distribución normal. En este caso, el valor p es menor que 0.05, lo que sugiere que hay evidencia para afirmar que los residuales no siguen una distribución normal 🤦🤦 .

x_seg <- df_residuals$residuals[pos_estadistico_m1]
y_sup_seg <- cdf_teorica[pos_estadistico_m1]
y_inf_seg <- cdf_empirica(x_seg)

cat("x_seg: ", format(x_seg, digits=5), 
    " y_inf_seg: ", format(y_inf_seg, digits=5),
    " y_sup_seg: ", format(y_sup_seg, digits=5),
    sep='')

## x_seg: -0.011272 y_inf_seg: 0.53659 y_sup_seg: 0.36836

datos_para_grafico <- data.frame(
  residuales = df_residuals$residuals
  )

ggplot(datos_para_grafico, aes(x = residuales)) +
#  geom_density(color = "blue", fill = "lightblue", alpha = 0.5) +
  stat_function(fun = pnorm, args = list(mean = mean_residuals,
                                         sd = sd_residuals),
                color = "red", size = 1) +
  stat_ecdf(color = "green", size = 1) +
  geom_segment(aes(x = x_seg, y = y_inf_seg, xend = x_seg,
                   yend = y_sup_seg),
               color = "purple", size = 1, linetype = "solid") +
  geom_point(aes(x = x_seg, y = y_inf_seg),
             color = "purple",size = 3) +
  geom_point(aes(x = x_seg, y = y_sup_seg),
             color = "purple", size = 3) +
  annotate("text", x = x_seg , y = (y_inf_seg + y_sup_seg)/2, 
           label = paste0("Máxima Separación=",
                          format(estadistico_m1,
                                 digits=5)," ", sep=''),
           color = "purple", size = 3, hjust = -0.05, vjust = 0) +
  coord_cartesian(xlim = c(min(datos_para_grafico$residuales),
                           max(datos_para_grafico$residuales))) +
  xlab("Residuales del modelo") +
  ylab("Probabilidad acumulada") +
  labs(title = "Esquema de funcionamiento de la prueba de Lilliefors") +
  theme_minimal()

#Prueba Shapiro-Wilk

# Probando normalidad del error. 

shapiro.test(mod1$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  mod1$residuals
## W = 0.90137, p-value = 0.00182

ks.test(mod1$residuals, "pnorm")

## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  mod1$residuals
## D = 0.48323, p-value = 0.000000002462
## alternative hypothesis: two-sided

Pruebas de Normalidad:

La prueba de normalidad Shapiro-Wilk y la prueba de Kolmogorov-Smirnov se utilizan para evaluar si los residuales del modelo de regresión se distribuyen de manera normal. Esto es importante porque muchos modelos de regresión asumen que los residuales siguen una distribución normal.

Resultados de la prueba Shapiro-Wilk:

El valor de W (estadístico de prueba de Shapiro-Wilk) es aproximadamente 0.90137, y el valor p es aproximadamente 0.00182 y el valor p es menor que el nivel de significancia previamente establecido (como 0.05), entonces se rechaza la hipótesis nula de que los residuales siguen una distribución normal. En este caso, el valor p es menor que 0.05, lo que sugiere que hay evidencia para afirmar que los residuales no siguen una distribución normal.

Resultados de la prueba de Kolmogorov-Smirnov:

El valor de D (estadístico de prueba de Kolmogorov-Smirnov) es aproximadamente 0.48323, y el valor p es aproximadamente 0.000000002462. Al igual que con la prueba de Shapiro-Wilk, si el valor p es menor que un nivel de significancia previamente establecido, se podría rechazar la hipótesis nula de que los residuales siguen una distribución normal. En este caso, el valor p es muy pequeño, lo que también sugiere que los residuales no siguen una distribución normal. 🤦🤦

# Primero inspeccionamos las variables individuales para observar su
# comportamiento temporal en un mismo gráfico

# library(lubridate)
x <- data_filtrados$capital
y <- data_filtrados$fncer


data_filtrados$predicciones <- mod1$fitted.values
data_filtrados$residuales <-mod1$residuals

3. EL SUPUESTO DE INDEPENDENCIA DEL ERROR

# Gráfico de residuales vs valores ajustados

data_filtrados %>% 
  ggplot(aes(x= predicciones, y=residuales)) +
  geom_point()+ 
  theme_light()+
  ylab("Residuales") +
  xlab("Valor predicho")

# Gráfico de residuales con índice temporal
data_filtrados %>% 
  ggplot(aes(x= capital, y=fncer)) +
  geom_point()+ 
  geom_line(aes(y = residuales, color = "residuales"), size = 1)+
  theme_light()+
  xlab("capital")+
  ylab("Residuales")

# Prueba de Durbin-Watson tomando el mismo orden en el que venían los
# datos en el dataframe puesto que este orden se corresponde con el 
# orden temporal.

resultado_dw <- dwtest(mod1)
print(resultado_dw)

## 
##  Durbin-Watson test
## 
## data:  mod1
## DW = 1.818, p-value = 0.2977
## alternative hypothesis: true autocorrelation is greater than 0

# Prueba de Ljung-Box
# Como esperamos cierta estacionalidad intra año en los datos podría tener sentido esperarla también para los residuales. Por lo tanto eligiremos lag = 1,2,3,...,12

p_values <- vector()
for (i in 1:12) {
  resultado_ljung_box <- Box.test(mod1$residuals, lag = i,
                                  type = "Ljung-Box")
  p_values[i] <-resultado_ljung_box$p.value
}

resultados_lb <- data.frame(
  lag = seq(1,12),
  p_value = p_values
)

print(resultado_ljung_box)

## 
##  Box-Ljung test
## 
## data:  mod1$residuals
## X-squared = 9.8268, df = 12, p-value = 0.6311

print(resultados_lb)

##    lag   p_value
## 1    1 0.6100487
## 2    2 0.8071955
## 3    3 0.9255195
## 4    4 0.2816650
## 5    5 0.3386285
## 6    6 0.4341509
## 7    7 0.4273534
## 8    8 0.5115144
## 9    9 0.6084303
## 10  10 0.5597444
## 11  11 0.6060921
## 12  12 0.6311490

Durbin-Watson test:

Esta prueba se utiliza para verificar la autocorrelación de los residuos en un modelo de regresión. La estadística de Durbin-Watson (DW) evalúa si hay autocorrelación positiva (valores cercanos a 0) o autocorrelación negativa (valores cercanos a 4) en los residuos. En tu caso, la estadística DW es 1.818, lo que está cerca de 2. Esto sugiere que no hay una autocorrelación significativa en los residuos del modelo.

Box-Ljung test:

Esta prueba se utiliza para evaluar si hay autocorrelación en los residuos de un modelo de regresión en múltiples rezagos (lags). La estadística X-cuadrado se compara con una distribución chi-cuadrado y su correspondiente valor p para determinar si hay evidencia de autocorrelación en los residuos. En tu caso, la estadística X-cuadrado es 9.8268 con 12 grados de libertad y un valor p de 0.6311. Un valor p alto (como en este caso) sugiere que no hay evidencia significativa de autocorrelación en los residuos a lo largo de múltiples rezagos.

Conclusión: según los resultados de estas pruebas, no parece haber autocorrelación significativa en los residuos de tu modelo de regresión. Esto es una buena señal, ya que la autocorrelación en los residuos podría indicar que el modelo no está capturando adecuadamente la estructura de los datos.

# Residuales estandarizados vs valores predichos
data_filtrados %>% 
  ggplot(aes(x= predicciones, y=residuales/sd(residuales))) +
  geom_point()+ 
  theme_light()+
  ylab("Residuales estandarizados") +
  xlab("Valor predicho fncer")

# Residuales estandarizados vs fncer
data_filtrados %>% 
  ggplot(aes(x= capital, y=residuales/sd(residuales))) +
  geom_point()+ 
  theme_light()+
  ylab("Residuales estandarizados") +
  ylab("capital")

Como se puede apreciar, el gráfico de errores estandarizados vs valores predichos será muy parecido al gráfico de errores vs variable predictora en un modelo de regresión lineal simple, puesto que en este contexto las predicciones dependen de la variable de pronóstico en forma tal que la primera no es más que un rescalado y desplazamiento de la segunda.

# Pruebas estadísticas formales
resultado_breusch_pagan <- bptest(mod1)
print(resultado_breusch_pagan)

## 
##  studentized Breusch-Pagan test
## 
## data:  mod1
## BP = 0.3371, df = 1, p-value = 0.5615

Si el valor p es alto (en este caso, p = 0.5615), no hay evidencia significativa para rechazar la hipótesis nula de homocedasticidad. En otras palabras, los datos sugieren que la varianza de los residuos es constante y no varía de manera heterocedástica en función de las variables independientes. Conclusión: según el resultado del “studentized Breusch-Pagan test”, no hay evidencia de heterocedasticidad en los residuos de tu modelo de regresión. Esto es positivo, ya que la homocedasticidad es una suposición importante en la regresión lineal.

resultado_white <- bptest(mod1,
  varformula = ~ (capital + I(capital^2)),
  data = data_filtrados)  
print(resultado_white)

## 
##  studentized Breusch-Pagan test
## 
## data:  mod1
## BP = 0.535, df = 2, p-value = 0.7653

El resultado del test Breusch-Pagan (también conocido como el test de White) se utiliza para evaluar la presencia de heterocedasticidad en un modelo de regresión. La heterocedasticidad se refiere a la variabilidad no constante de los errores en un modelo de regresión, lo que significa que la varianza de los errores no es la misma en todos los niveles de las variables predictoras.

BP = 0.535 es el estadístico de prueba de Breusch-Pagan. df = 2 son los grados de libertad del estadístico de prueba. p-value = 0.7653 es el valor p asociado al estadístico de prueba.

Hipótesis nula (H0): No hay heterocedasticidad en el modelo (la varianza de los errores es constante). Hipótesis alternativa (H1): Hay evidencia de heterocedasticidad en el modelo (la varianza de los errores no es constante).

Dado que el valor p (0.7653) es mayor que un nivel de significancia = 0.05, no tenemos evidencia suficiente para rechazar la hipótesis nula. En otras palabras, no hay suficiente evidencia para concluir que existe heterocedasticidad en tu modelo.

Conclusión: según los resultados del test Breusch-Pagan, no hay indicios de heterocedasticidad en tu modelo de regresión, ya que el valor p es mayor que el nivel de significancia seleccionado. Esto sugiere que la suposición de homocedasticidad (varianza constante de los errores) es razonable.

Para ilustrar como trabajan estas pruebas por debajo, considere la siguiente implementación manual:

# Implementación manual del test de White

#1. Construir el dataframe que relaciona el cuadrado del error con
#los términos lineales de la predictora y el término cuadrático
#(por ser modelo lineal simple no hay productos cruzados)

# 1. Crear un dataframe con el cuadrado del error y las variables del modelo
df_squared_error <- data.frame(
  squared_error = mod1$residuals^2,
  capital = data_filtrados$capital,
  capital_2 = data_filtrados$capital^2
)

# 2. Ajustar un modelo lineal sobre el cuadrado del error regresado por el término lineal y cuadrático de la variable capital
modelo_white <- lm(data = df_squared_error,
                   squared_error ~ capital + capital_2)

# 3. Calcular el estadístico de prueba como n * R^2 (con n = número de datos del dataframe de errores)
R2 <- summary(modelo_white)$r.squared
n <- nrow(df_squared_error)
w_estadistico <- n * R2

# 4. Calcular el valor p utilizando la distribución chi-cuadrado con 2 grados de libertad
p_valor <- 1 - pchisq(w_estadistico, df = 2)

# Imprimir el estadístico de prueba y el valor p
paste("BP =", round(w_estadistico, 4), "y p_valor =", round(p_valor, 5))

## [1] "BP = 0.535 y p_valor = 0.76529"

En la salida impresa puede notar que los resultados son consistentes con lo generado por la prueba de heterocedasticidad arrojada con la función bptest.

# 1. Crear un dataframe con el cuadrado del error y las variables del modelo
df_squared_error <- data.frame(
  squared_error = mod1$residuals^2,
  capital = data_filtrados$capital,
  capital_2 = data_filtrados$capital^2
)

# 2. Ajustar un modelo lineal sobre el cuadrado del error regresado por el término lineal y cuadrático de la variable capital
modelo_white <- lm(data = df_squared_error,
                   squared_error ~ capital + capital_2)

# 3. Calcular el valor estimado del error cuadrado
estimated_squared_error <- modelo_white$coefficients[1] +
                          modelo_white$coefficients[2] * df_squared_error$capital +
                          modelo_white$coefficients[3] * df_squared_error$capital_2

# 4. Crear un gráfico
library(ggplot2)

p <- ggplot(data = df_squared_error, aes(x = capital, y = squared_error)) +
  geom_point() +
  geom_line(aes(y = estimated_squared_error, color = "red")) +
  labs(
    title = "Modelo white sobre cuadrados del error",
    subtitle = paste(
      "R² =", round(R2, 4),
      "W-Estadístico =", round(w_estadistico, 4),
      "p_valor =", round(p_valor, 5)
    )
  ) +
  scale_color_manual(
    values = c("red" = "blue"),
    labels = c("Error cuadrado estimado"),
    name = "Modelo de White"
  )

print(p)

conclusiones generales:

Modelo de Regresión Lineal Múltiple: El modelo de regresión lineal múltiple parece mostrar una relación estadísticamente significativa entre la variable dependiente (fncer, generación de energía eléctrica de fuentes renovables) y la variable independiente (capital, capital de los estados de Estados Unidos). Los coeficientes del modelo indican que, en promedio, un aumento en el capital de los estados se asocia con un aumento en la generación de energía eléctrica de fuentes renovables, aunque la magnitud de este efecto es relativamente pequeña.

Normalidad de los Residuos: Los resultados de las pruebas de normalidad de los residuos (Lilliefors y Shapiro-Wilk) sugieren que los residuos del modelo no siguen una distribución normal. Frente a esta violación de uno de los supuestos principales de la regresión lineal, se debe concluir que el modelo presentado es inservible para la serie de datos elegida.

Autocorrelación de los Residuos: El resultado del test de Durbin-Watson sugiere que los residuos del modelo no muestran autocorrelación positiva significativa.

Homocedasticidad de los Residuos: El resultado del test de Breusch-Pagan no muestra evidencia significativa de heteroscedasticidad en los residuos. Esto indica que la varianza de los errores parece ser constante en diferentes niveles de las variables independientes.

Modelo de White: El modelo de White sugiere que podría haber heteroscedasticidad condicional en el modelo, lo que podría indicar que la varianza de los errores está relacionada con la variable capital y capital^2. Sin embargo, el p-valor alto en la prueba de White sugiere que no hay suficiente evidencia para rechazar la hipótesis nula de homocedasticidad.

Análisis de Gráficos: Los gráficos generados muestran la relación entre las variables y cómo se ajusta el modelo. Los puntos dispersos y la línea de regresión en los gráficos ayudan a visualizar cómo se relacionan las variables.

En resumen, aunque el modelo de regresión lineal múltiple muestra una relación entre las variables capital y fncer, es importante tener en cuenta las violaciones a los supuestos fundamentales, como la falta de normalidad en los residuos, por lo que se descarta el presente modelo.