1. Librerías

library(readxl)
library(dplyr)
library(gt)
library(plotly)
library(DT)

2. Carga de datos

datos <- read_excel("dataset_mundial_petro.xlsx")
cat("Número de registros (filas crudas):", nrow(datos), "\n")
## Número de registros (filas crudas): 49212
cat("Número de variables disponibles:", ncol(datos), "\n")
## Número de variables disponibles: 32

El archivo trae múltiples filas por yacimiento (una fila distinta por cada registro de reservas/producción asociado a un mismo campo), por eso más adelante se agrupa por Unit ID antes de construir la tabla del modelo.

3. Selección de variables

Y = PRODUCTION_START_YEAR (efecto: año en que el yacimiento empieza a producir)

Se seleccionan tres variables independientes con lógica causal, no solo por su valor de correlación:

  • X1 = DISCOVERY_YEAR (Discovery year): año en que se descubrió el yacimiento.
  • X2 = FID_YEAR (FID Year): año de la Decisión Final de Inversión.
  • X3 = STATUS_YEAR (Status year): año del último estado operativo registrado del yacimiento (activo, cerrado, en construcción, etc.).

Un yacimiento se descubre (X1), luego recibe aprobación de inversión (X2), y a lo largo de su vida operativa acumula un estado registrado (X3) antes o junto con el arranque de producción (Y). Por eso X1, X2 y X3 se tratan como variables explicativas y Y como el efecto.

FID Year y Production start year llegan como texto porque incluyen valores como “2025 (expected)”; se extrae el año numérico de esos casos con una función auxiliar.

extraer_anio <- function(v) as.numeric(sub("^(\\d{4}).*", "\\1", as.character(v)))

datos <- datos %>%
  mutate(
    X1_raw = as.numeric(`Discovery year`),
    X2_raw = extraer_anio(`FID Year`),
    X3_raw = as.numeric(`Status year`),
    Y_raw  = extraer_anio(`Production start year`)
  )

4. Dominios de las variables

Antes de ajustar cualquier modelo, se define el dominio de cada variable, ya que el modelo solo se interpreta dentro de estos rangos, nunca fuera de ellos.

resumen_dominios <- data.frame(
  Variable  = c("X1 (DISCOVERY_YEAR)", "X2 (FID_YEAR)", "X3 (STATUS_YEAR)",
                "Y (PRODUCTION_START_YEAR)"),
  Dominio_teorico = c(
    "Enteros positivos (Z+): un año calendario no puede ser negativo, fraccionario ni anterior al origen de la industria petrolera",
    "Enteros positivos (Z+): puede incluir años próximos porque hay proyectos con FID ya aprobado a futuro",
    "Enteros positivos (Z+): fecha de la última actualización de estado del yacimiento",
    "Enteros positivos (Z+): año calendario de arranque de producción; no existen registros de años negativos"
  ),
  Rango_observado = c(
    paste0("[", min(datos$X1_raw, na.rm = TRUE), ", ", max(datos$X1_raw, na.rm = TRUE), "]"),
    paste0("[", min(datos$X2_raw, na.rm = TRUE), ", ", max(datos$X2_raw, na.rm = TRUE), "]"),
    paste0("[", min(datos$X3_raw, na.rm = TRUE), ", ", max(datos$X3_raw, na.rm = TRUE), "]"),
    paste0("[", min(datos$Y_raw, na.rm = TRUE), ", ", max(datos$Y_raw, na.rm = TRUE), "]")
  )
) %>%
  gt() %>%
  tab_header(title = md("**Dominio teórico y rango observado de cada variable**"))

resumen_dominios
Dominio teórico y rango observado de cada variable
Variable Dominio_teorico Rango_observado
X1 (DISCOVERY_YEAR) Enteros positivos (Z+): un año calendario no puede ser negativo, fraccionario ni anterior al origen de la industria petrolera [1869, 2023]
X2 (FID_YEAR) Enteros positivos (Z+): puede incluir años próximos porque hay proyectos con FID ya aprobado a futuro [1975, 2029]
X3 (STATUS_YEAR) Enteros positivos (Z+): fecha de la última actualización de estado del yacimiento [1989, 2024]
Y (PRODUCTION_START_YEAR) Enteros positivos (Z+): año calendario de arranque de producción; no existen registros de años negativos [1896, 2040]

Sobre Y en particular: como Y representa años, su dominio se restringe a los enteros positivos dentro del rango observado. El modelo solo se usa dentro de ese rango, y ahí Y siempre se mantiene positivo.

5. Tabla de valores (agrupada y balanceada)

En la versión anterior, algunas variables tenían muchos más datos disponibles que otras (por ejemplo, FID Year solo tenía unos pocos cientos de valores frente a miles del resto), lo que generaba tablas desbalanceadas. Esto ocurre porque el archivo trae varias filas por yacimiento (una por cada registro de reservas). Para corregirlo:

  1. Se agrupa por Unit ID (identificador único de yacimiento) y se toma un solo valor representativo de cada variable por yacimiento.
  2. Solo después se filtran los yacimientos que tienen las tres variables completas a la vez (X1, X2, X3, Y), de modo que todas las filas de la tabla final tengan exactamente la misma cantidad de datos disponibles.
agregado <- datos %>%
  group_by(`Unit ID`) %>%
  summarise(
    DISCOVERY_YEAR        = first(X1_raw),
    FID_YEAR               = first(X2_raw),
    STATUS_YEAR             = first(X3_raw),
    PRODUCTION_START_YEAR   = first(Y_raw),
    .groups = "drop"
  )

tabla_original <- agregado %>% select(-`Unit ID`)
cat("Yacimientos únicos disponibles:", nrow(tabla_original), "\n")
## Yacimientos únicos disponibles: 8334

5.1 Depuración: solo yacimientos con las 3 variables y Y completas

tabla <- tabla_original %>%
  filter(!is.na(DISCOVERY_YEAR), !is.na(FID_YEAR), !is.na(STATUS_YEAR),
         !is.na(PRODUCTION_START_YEAR)) %>%
  filter(DISCOVERY_YEAR > 0, FID_YEAR > 0, STATUS_YEAR > 0,
         PRODUCTION_START_YEAR > 0)

cat("Yacimientos con las 4 variables completas y válidas:", nrow(tabla), "\n")
## Yacimientos con las 4 variables completas y válidas: 259
# Verificación de balance: todas las columnas deben tener el mismo número de datos
sapply(tabla, function(col) sum(!is.na(col)))
##        DISCOVERY_YEAR              FID_YEAR           STATUS_YEAR 
##                   259                   259                   259 
## PRODUCTION_START_YEAR 
##                   259
datatable(tabla,
          options = list(pageLength = 10),
          colnames = c("Año de Descubrimiento (X1)", "Año de FID (X2)",
                       "Año del Último Estado (X3)",
                       "Año de Inicio de Producción (Y)"))

6. Gráfica de dispersión

Primero se muestra cada X por separado contra Y (más fácil de leer una por una), y luego un 3D con las tres variables juntas.

g1 <- plot_ly(tabla, x = ~DISCOVERY_YEAR, y = ~PRODUCTION_START_YEAR,
              type = "scatter", mode = "markers",
              marker = list(size = 5, color = "#2E86C1", opacity = 0.6)) %>%
  layout(xaxis = list(title = "X1 Discovery"), yaxis = list(title = "Y"))

g2 <- plot_ly(tabla, x = ~FID_YEAR, y = ~PRODUCTION_START_YEAR,
              type = "scatter", mode = "markers",
              marker = list(size = 5, color = "#2E86C1", opacity = 0.6)) %>%
  layout(xaxis = list(title = "X2 FID"), yaxis = list(title = "Y"))

g3 <- plot_ly(tabla, x = ~STATUS_YEAR, y = ~PRODUCTION_START_YEAR,
              type = "scatter", mode = "markers",
              marker = list(size = 5, color = "#2E86C1", opacity = 0.6)) %>%
  layout(xaxis = list(title = "X3 Status"), yaxis = list(title = "Y"))

subplot(g1, g2, g3, nrows = 1, shareX = FALSE, shareY = TRUE,
        titleX = TRUE, titleY = TRUE, margin = 0.05) %>%
  layout(showlegend = FALSE)
plot_ly(tabla, x = ~DISCOVERY_YEAR, y = ~FID_YEAR, z = ~PRODUCTION_START_YEAR,
        type = "scatter3d", mode = "markers",
        marker = list(size = 4, color = ~STATUS_YEAR, colorscale = "Viridis", showscale = TRUE)) %>%
  layout(scene = list(
    xaxis = list(title = "X1 DISCOVERY_YEAR"),
    yaxis = list(title = "X2 FID_YEAR"),
    zaxis = list(title = "Y PRODUCTION_START_YEAR")
  ))

El color de los puntos en el 3D representa a X3 (STATUS_YEAR), la tercera variable que no cabe en los ejes.

7. Conjetura

Los puntos crecen de forma conjunta con X2 principalmente, y X1 y X3 aportan variación adicional sin romper la tendencia general. Se prueba con regresión lineal múltiple de 3 variables:

\[y = b_0 + b_1 x_1 + b_2 x_2 + b_3 x_3\]

8. Cálculo de parámetros

modelo <- lm(PRODUCTION_START_YEAR ~ DISCOVERY_YEAR + FID_YEAR + STATUS_YEAR,
             data = tabla)
summary(modelo)
## 
## Call:
## lm(formula = PRODUCTION_START_YEAR ~ DISCOVERY_YEAR + FID_YEAR + 
##     STATUS_YEAR, data = tabla)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -25.474  -1.264   0.004   1.307   7.709 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -709.74418  286.59299  -2.476   0.0139 *  
## DISCOVERY_YEAR   -0.01276    0.01293  -0.986   0.3248    
## FID_YEAR          0.99566    0.01774  56.126   <2e-16 ***
## STATUS_YEAR       0.36935    0.14620   2.526   0.0121 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.855 on 255 degrees of freedom
## Multiple R-squared:  0.9584, Adjusted R-squared:  0.9579 
## F-statistic:  1960 on 3 and 255 DF,  p-value: < 2.2e-16
b0 <- coef(modelo)[1]
b1 <- coef(modelo)[2]
b2 <- coef(modelo)[3]
b3 <- coef(modelo)[4]

data.frame(
  Parámetro = c("b0 (intercepto)", "b1 (pendiente X1)", "b2 (pendiente X2)",
                "b3 (pendiente X3)"),
  Valor     = c(round(b0, 4), round(b1, 4), round(b2, 4), round(b3, 4))
) %>%
  gt() %>%
  tab_header(title = md("**Parámetros del Modelo**"))
Parámetros del Modelo
Parámetro Valor
b0 (intercepto) -709.7442
b1 (pendiente X1) -0.0128
b2 (pendiente X2) 0.9957
b3 (pendiente X3) 0.3693

Sobre el signo de b1: X1 (DISCOVERY_YEAR) sale con pendiente negativa; el peso explicativo principal lo lleva X2. Esto no es un problema porque, como se definió en la sección 4, el modelo solo se usa dentro del rango de X1 (1938 a 2022), y en ese rango Y se mantiene siempre positivo.

9. Realidad y modelo

x1_seq <- seq(min(tabla$DISCOVERY_YEAR), max(tabla$DISCOVERY_YEAR), length.out = 25)
x2_seq <- seq(min(tabla$FID_YEAR), max(tabla$FID_YEAR), length.out = 25)
x3_medio <- mean(tabla$STATUS_YEAR)
malla  <- outer(x1_seq, x2_seq, function(a, b) b0 + b1 * a + b2 * b + b3 * x3_medio)

plot_ly() %>%
  add_markers(data = tabla, x = ~DISCOVERY_YEAR, y = ~FID_YEAR, z = ~PRODUCTION_START_YEAR,
              marker = list(size = 4, color = "#2E86C1"), name = "Datos reales") %>%
  add_surface(x = x1_seq, y = x2_seq, z = t(malla), opacity = 0.6,
              showscale = FALSE, name = "Modelo (X3 en su media)") %>%
  layout(scene = list(
    xaxis = list(title = "X1 DISCOVERY_YEAR"),
    yaxis = list(title = "X2 FID_YEAR"),
    zaxis = list(title = "Y PRODUCTION_START_YEAR")
  ))

Como X3 también forma parte del modelo, la superficie de arriba se calcula dejando X3 fija en su promedio para poder graficarla en 3D. Para ver el efecto de las 3 variables a la vez, se compara además lo predicho contra lo observado:

tabla$Y_predicho <- predict(modelo)

plot_ly(tabla, x = ~PRODUCTION_START_YEAR, y = ~Y_predicho,
        type = "scatter", mode = "markers",
        marker = list(size = 6, color = "#2E86C1")) %>%
  add_trace(x = ~PRODUCTION_START_YEAR, y = ~PRODUCTION_START_YEAR,
             mode = "lines", line = list(color = "red", dash = "dash"),
             name = "Ajuste perfecto") %>%
  layout(title = "Y observado vs. Y predicho",
         xaxis = list(title = "Y observado"),
         yaxis = list(title = "Y predicho"))

Mientras más cerca estén los puntos de la línea roja, mejor es el ajuste del modelo.

10. Test

pearson_x1 <- cor(tabla$DISCOVERY_YEAR, tabla$PRODUCTION_START_YEAR)
pearson_x2 <- cor(tabla$FID_YEAR, tabla$PRODUCTION_START_YEAR)
pearson_x3 <- cor(tabla$STATUS_YEAR, tabla$PRODUCTION_START_YEAR)
f_pvalue   <- pf(summary(modelo)$fstatistic[1],
                 summary(modelo)$fstatistic[2],
                 summary(modelo)$fstatistic[3], lower.tail = FALSE)
r2         <- summary(modelo)$r.squared

data.frame(
  Test      = c("Pearson X1 (DISCOVERY_YEAR)", "Pearson X2 (FID_YEAR)",
                "Pearson X3 (STATUS_YEAR)", "F-test modelo"),
  Valor     = c(round(pearson_x1, 4), round(pearson_x2, 4), round(pearson_x3, 4),
                format(f_pvalue, scientific = TRUE, digits = 4)),
  Resultado = c(ifelse(abs(pearson_x1) > 0.7, "Supera 0.7", "No supera 0.7"),
                ifelse(abs(pearson_x2) > 0.7, "Supera 0.7", "No supera 0.7"),
                ifelse(abs(pearson_x3) > 0.7, "Supera 0.7", "No supera 0.7"),
                ifelse(f_pvalue < 0.05, "Significativo", "No significativo"))
) %>%
  gt() %>%
  tab_header(title = md("**Resumen del Test**"))
Resumen del Test
Test Valor Resultado
Pearson X1 (DISCOVERY_YEAR) 0.595 No supera 0.7
Pearson X2 (FID_YEAR) 0.9783 Supera 0.7
Pearson X3 (STATUS_YEAR) 0.399 No supera 0.7
F-test modelo 9.661e-176 Significativo
 
cat("R² del modelo:", round(r2, 4), "\n")
## R² del modelo: 0.9584

Nota honesta: de las 3 variables, solo X2 (FID_YEAR) supera 0.7 de correlación con Y por sí sola; X1 y X3 no llegan a ese nivel individualmente. Aun así, el modelo conjunto explica 95.8% de la variabilidad de Y (R² = 0.9584), y el F-test confirma que el modelo en conjunto es significativo.

11. Restricciones

  • El modelo es válido dentro del dominio muestral definido en la sección 4:
    • X1 (DISCOVERY_YEAR): 1938 – 2022
    • X2 (FID_YEAR): 1975 – 2029
    • X3 (STATUS_YEAR): 2018 – 2024
    • Y (PRODUCTION_START_YEAR): 1975 – 2034 (siempre entero positivo, como corresponde a un año)
  • No se debe extrapolar el modelo fuera de estos rangos: al ser Y un año calendario, su dominio real son los enteros positivos, y el modelo solo tiene sentido dentro del rango observado.
  • La relación podría dejar de comportarse de forma lineal fuera de este rango (por ejemplo, por cambios regulatorios o tecnológicos entre décadas).
  • Se usaron únicamente yacimientos con las 4 variables completas y reales (sin relleno artificial); no sirve para predecir yacimientos con datos faltantes en X1, X2 o X3.

12. Estimación

escenarios <- data.frame(
  Escenario = c("Típico (medianas)", "Descubrimiento antiguo", "FID reciente", "Estado más antiguo"),
  DISCOVERY_YEAR = c(median(tabla$DISCOVERY_YEAR), min(tabla$DISCOVERY_YEAR),
                      median(tabla$DISCOVERY_YEAR), median(tabla$DISCOVERY_YEAR)),
  FID_YEAR        = c(median(tabla$FID_YEAR), median(tabla$FID_YEAR),
                      max(tabla$FID_YEAR), median(tabla$FID_YEAR)),
  STATUS_YEAR     = c(median(tabla$STATUS_YEAR), median(tabla$STATUS_YEAR),
                      median(tabla$STATUS_YEAR), min(tabla$STATUS_YEAR))
)

escenarios$PRODUCTION_START_YEAR_ESTIMADA <- predict(modelo, newdata = escenarios)

escenarios %>%
  mutate(PRODUCTION_START_YEAR_ESTIMADA = round(PRODUCTION_START_YEAR_ESTIMADA, 2)) %>%
  gt() %>%
  tab_header(title = md("**Estimaciones del Modelo**"))
Estimaciones del Modelo
Escenario DISCOVERY_YEAR FID_YEAR STATUS_YEAR PRODUCTION_START_YEAR_ESTIMADA
Típico (medianas) 2005 2019 2022 2021.74
Descubrimiento antiguo 1938 2019 2022 2022.60
FID reciente 2005 2029 2022 2031.70
Estado más antiguo 2005 2019 2018 2020.26

13. Conclusión

Con DISCOVERY_YEAR (X1), FID_YEAR (X2) y STATUS_YEAR (X3) se obtiene el modelo de regresión lineal múltiple:

\[y = -709.74 + -0.0128x_1 + 0.9957x_2 + 0.3693x_3\]

Interpretación de cada pendiente (manteniendo las demás variables constantes):

  • Por cada año adicional en X1 (descubrimiento), Y cambia en -0.0128.
  • Por cada año adicional en X2 (FID), Y cambia en 0.9957.
  • Por cada año adicional en X3 (último estado), Y cambia en 0.3693.

El modelo explica 95.8% de la variabilidad de Y (R² = 0.9584). El F-test resulta significativo (p = 9.66e-176), lo que confirma que, en conjunto, las 3 variables sí aportan información útil para explicar Y — aunque, tomadas de forma individual, X1 y X3 no superen 0.7 de correlación con Y.

El dominio de Y se definió como los enteros positivos dentro del rango observado (1975–2034), y el modelo se interpreta únicamente dentro de ese dominio.