1 Configuración y Carga de Datos

##### UNIVERSIDAD CENTRAL DEL ECUADOR #####
#### AUTOR: ESTALIN ####
### CARRERA: INGENIERÍA EN PETRÓLEOS ###

#### MODELO DE REGRESIÓN EXPONENCIAL ####
## DATASET: Global Oil and Gas Extraction Tracker (GOGET) ##

library(readxl)
library(dplyr)
library(ggplot2)
library(gt)

# Selección del archivo Excel
ruta_archivo <- file.choose()
# guess_max se sube porque las primeras filas del archivo son solo metadatos
# (sin reservas ni unidades); con el valor por defecto (1000), readxl
# adivina mal el tipo de esas columnas y termina leyendo todo como NA.
Datos <- read_excel(ruta_archivo, guess_max = 50000)

## Estructura de los datos
str(Datos)
## tibble [49,212 × 32] (S3: tbl_df/tbl/data.frame)
##  $ Unit ID                           : chr [1:49212] "OG0000001" "OG0000002" "OG0000006" "OG0000007" ...
##  $ Unit Name                         : chr [1:49212] "Matzen" "Abalone" "Aguilhada" "Agulha" ...
##  $ Unit name local script            : chr [1:49212] NA "Abalone" "Aguilhada" "Agulha" ...
##  $ Fuel type                         : chr [1:49212] "oil and gas" "oil and gas" "oil and gas" "oil and gas" ...
##  $ Unit type                         : chr [1:49212] "field" "field" "field" "field" ...
##  $ Country                           : chr [1:49212] "Austria" "Brazil" "Brazil" "Brazil" ...
##  $ Subnational unit (province, state): chr [1:49212] NA "Espírito Santo" "Sergipe" "Rio Grande do Norte" ...
##  $ Latitude                          : num [1:49212] 48.4 -21.4 -10.7 -4.9 -22.1 ...
##  $ Longitude                         : num [1:49212] 16.7 -39.6 -36.9 -36.3 -40 ...
##  $ Location accuracy                 : chr [1:49212] "approximate" "exact" "exact" "exact" ...
##  $ Status                            : chr [1:49212] "operating" "operating" "operating" "operating" ...
##  $ Status year                       : num [1:49212] 2023 2022 2022 2022 2022 ...
##  $ Discovery year                    : num [1:49212] 1949 2001 1966 1975 1984 ...
##  $ FID Year                          : chr [1:49212] NA NA NA NA ...
##  $ Production start year             : chr [1:49212] "1951" "2009" "1969" "1979" ...
##  $ Operator                          : chr [1:49212] "OMV" "Shell Brasil Petróleo Ltda." NA NA ...
##  $ Owner                             : chr [1:49212] "OMV (100%)" "Shell Brasil (50%);ONGC Campos (27%);Qatarenergy (23%)" "Petrobras (100%)" "Petrobras (100%)" ...
##  $ Parent                            : chr [1:49212] "OMV Aktiengesellschaft (100%)" "Shell plc (50%);Oil and Natural Gas Corporation (ONGC) (27%)" "Petróleo Brasileiro S.A. (100%)" "Petróleo Brasileiro S.A. (100%)" ...
##  $ Basin                             : chr [1:49212] NA NA NA NA ...
##  $ Concession / block                : chr [1:49212] NA NA NA NA ...
##  $ Project or complex                : chr [1:49212] "Matzen" NA NA NA ...
##  $ Government unit ID                : chr [1:49212] NA NA NA NA ...
##  $ Wiki URL                          : chr [1:49212] "https://www.gem.wiki/Matzen_Oil_and_Gas_Field_(Austria)" "https://www.gem.wiki/Abalone_Oil_and_Gas_Field_%28Esp%C3%ADrito_Santo%2C_Brazil%29" "https://www.gem.wiki/Aguilhada_Oil_and_Gas_Field_%28Sergipe%2C_Brazil%29" "https://www.gem.wiki/Agulha_Oil_and_Gas_Field_%28Rio_Grande_do_Norte%2C_Brazil%29" ...
##  $ Unit name2                        : chr [1:49212] NA NA NA NA ...
##  $ Production/reserves               : chr [1:49212] NA NA NA NA ...
##  $ Fuel description                  : chr [1:49212] NA NA NA NA ...
##  $ Reserves classification (original): chr [1:49212] NA NA NA NA ...
##  $ Quantity (original)               : num [1:49212] NA NA NA NA NA NA NA NA NA NA ...
##  $ Units (original)                  : chr [1:49212] NA NA NA NA ...
##  $ Data year                         : chr [1:49212] NA NA NA NA ...
##  $ Quantity (converted)              : num [1:49212] NA NA NA NA NA NA NA NA NA NA ...
##  $ Units (converted)                 : chr [1:49212] NA NA NA NA ...

2 Extracción y Depuración de Variables

Se aplica un modelo de crecimiento/decrecimiento exponencial sobre los datos agrupados mediante técnica de binning, para reducir el ruido estadístico y capturar la tendencia general del comportamiento de las reservas petroleras a lo largo del tiempo.

Se estableció el Año de Descubrimiento (Discovery year) como variable independiente (x), por representar el momento histórico en que cada yacimiento fue identificado. Las Reservas (Quantity converted, en millones de barriles) actúan como variable dependiente (y), ya que expresan el volumen de hidrocarburo estimado para ese yacimiento.

Esta relación busca modelar cómo ha evolucionado el tamaño de los yacimientos descubiertos a lo largo de la historia de la exploración petrolera: la hipótesis de trabajo es que los campos gigantes fueron descubiertos preferentemente en las primeras décadas de exploración, mientras que los descubrimientos más recientes tienden a corresponder a yacimientos de menor tamaño remanente, lo que generaría una tendencia de tipo exponencial decreciente.

Nota sobre datos faltantes: por indicación del docente, los registros con valores faltantes no se eliminan. En su lugar, tanto el año de descubrimiento como la reserva se completan con la mediana de los valores disponibles, preservando así el tamaño de la muestra.

  • Variable Independiente (X): Año de Descubrimiento (agrupado cada 10 años).
  • Variable Dependiente (Y): Reservas Promedio (millones de barriles).
# Año de descubrimiento por yacimiento (dato único a nivel de campo)
disc_por_unidad <- Datos %>%
  group_by(`Unit ID`) %>%
  summarise(Descubrimiento = suppressWarnings(max(as.numeric(`Discovery year`), na.rm = TRUE))) %>%
  mutate(Descubrimiento = ifelse(is.infinite(Descubrimiento), NA, Descubrimiento))

# Reservas máximas reportadas (millones de barriles) por yacimiento
reservas_por_unidad <- Datos %>%
  filter(`Production/reserves` == "reserves", `Units (converted)` == "million bbl") %>%
  group_by(`Unit ID`) %>%
  summarise(Reservas = suppressWarnings(max(as.numeric(`Quantity (converted)`), na.rm = TRUE))) %>%
  mutate(Reservas = ifelse(is.infinite(Reservas), NA, Reservas))

# Se parte de TODOS los yacimientos del dataset (no se elimina ninguno)
todos_yacimientos <- Datos %>% distinct(`Unit ID`)

datos_unidad <- todos_yacimientos %>%
  left_join(disc_por_unidad, by = "Unit ID") %>%
  left_join(reservas_por_unidad, by = "Unit ID")

# Imputación con la mediana (en lugar de eliminar registros con NA)
mediana_disc <- median(datos_unidad$Descubrimiento, na.rm = TRUE)
mediana_res  <- median(datos_unidad$Reservas, na.rm = TRUE)

datos_raw <- datos_unidad %>%
  mutate(
    x_raw = ifelse(is.na(Descubrimiento), mediana_disc, Descubrimiento),
    y_raw = ifelse(is.na(Reservas), mediana_res, Reservas)
  ) %>%
  filter(y_raw > 0)

cat("Yacimientos totales:", nrow(todos_yacimientos), "\n")
## Yacimientos totales: 8334
cat("Mediana usada para anio de descubrimiento:", mediana_disc, "\n")
## Mediana usada para anio de descubrimiento: 1976
cat("Mediana usada para reservas (millones bbl):", round(mediana_res, 2), "\n")
## Mediana usada para reservas (millones bbl): 74.51
cat("Registros finales tras imputacion:", nrow(datos_raw), "\n")
## Registros finales tras imputacion: 8290

3 Análisis Gráfico Exploratorio

datos_plot <- datos_raw %>% filter(y_raw < quantile(y_raw, 0.99))

par(mar = c(5, 5, 4, 2))
color_trans <- rgb(0.2, 0.6, 0.86, 0.4)

plot(datos_plot$x_raw, datos_plot$y_raw,
     main = "Gráfica N.1: Diagrama de Dispersión de las Reservas (millones bbl)\nen función del Año de Descubrimiento",
     xlab = "Año de Descubrimiento",
     ylab = "Reservas (millones bbl)",
     col = color_trans, pch = 16, cex = 0.6, cex.main = 0.9, frame.plot = FALSE)

grid(nx = NULL, ny = NULL, col = "#D7DBDD", lty = "dotted")

4 Aplicación de Binning

Debido a la alta variabilidad observada en la Gráfica N.1 y al efecto de la imputación por mediana, se implementa la técnica de binning con el objetivo de disminuir el ruido estadístico y visualizar con mayor claridad la tendencia general de los datos.

# Agrupación técnica (cada 10 años)
datos_agrupados <- datos_raw %>%
  mutate(disc_bin = floor(x_raw / 10) * 10) %>%
  group_by(disc_bin) %>%
  summarise(
    y = mean(y_raw, na.rm = TRUE),
    conteo = n()
  ) %>%
  rename(x = disc_bin) %>%
  filter(conteo >= 5)

x_val <- datos_agrupados$x
y_val <- datos_agrupados$y

datos_agrupados
## # A tibble: 15 × 3
##        x      y conteo
##    <dbl>  <dbl>  <int>
##  1  1860 2513.       6
##  2  1880   74.5      6
##  3  1900  157.      53
##  4  1910  162.      87
##  5  1920 2514.      81
##  6  1930  795.     154
##  7  1940 1013.     292
##  8  1950 1191.     609
##  9  1960 1009.     666
## 10  1970  266.    4091
## 11  1980  408.     619
## 12  1990  360.     490
## 13  2000  333.     588
## 14  2010  562.     430
## 15  2020  271.     114

5 Conjetura del Modelo de Regresión Exponencial

La ecuación es \(y = a \cdot e^{b \cdot x}\). Para linealizar, aplicamos logaritmo solo a la variable dependiente: \(\ln(y) = \ln(a) + b \cdot x\).

# Linealización
y_log <- log(y_val)
modelo_exponencial <- lm(y_log ~ x_val)

# Parámetros
log_a <- coef(modelo_exponencial)[1]
b_param <- coef(modelo_exponencial)[2]
a_param <- exp(log_a)

6 Gráfica del Modelo Exponencial

Se presenta el ajuste del modelo incluyendo la banda de incertidumbre estadística (Intervalo de Confianza del 95%).

par(mar = c(5, 5, 4, 2))

plot(x_val, y_val,
     main = "Gráfica N.2: Modelo Exponencial de las Reservas (millones bbl)\nen función del Año de Descubrimiento",
     xlab = "Año de Descubrimiento",
     ylab = "Reservas (millones bbl)",
     col = "#3498DB", pch = 16, cex = 1.0, cex.main = 0.9, frame.plot = FALSE)

grid(nx = NULL, ny = NULL, col = "#D7DBDD", lty = "dotted")

x_seq <- seq(min(x_val), max(x_val), length.out = 500)

pred_log <- predict(modelo_exponencial,
                     newdata = data.frame(x_val = x_seq),
                     interval = "confidence",
                     level = 0.95)

y_pred_fit <- exp(pred_log[, "fit"])
y_pred_lwr <- exp(pred_log[, "lwr"])
y_pred_upr <- exp(pred_log[, "upr"])

# Intervalo de Confianza
polygon(c(x_seq, rev(x_seq)),
        c(y_pred_lwr, rev(y_pred_upr)),
        col = rgb(0.5, 0.5, 0.5, 0.2), border = NA)

lines(x_seq, y_pred_fit, col = "#E74C3C", lwd = 3)

legend("topleft",
       legend = c("Datos promediados (binning)",
                  "Modelo Exponencial",
                  "I.C. 95%"),
       col = c("#3498DB", "#E74C3C", "gray"),
       pch = c(16, NA, 15),
       lwd = c(NA, 3, NA),
       pt.cex = c(1, NA, 2),
       bty = "n")

7 Test de Bondad del Modelo

7.1 Coeficiente de correlación del modelo linealizado

r <- cor(x_val, log(y_val))
cat("El coeficiente de correlacion es: ", round(r, 2))
## El coeficiente de correlacion es:  -0.08

7.2 Coeficiente de determinación

r2 <- summary(modelo_exponencial)$r.squared
cat(paste0("El coeficiente de determinacion (R2) es: ", round(r2, 2)))
## El coeficiente de determinacion (R2) es: 0.01

8 Ecuación del Modelo

# Coeficientes del modelo
log_a <- coef(modelo_exponencial)[1]
b_param <- coef(modelo_exponencial)[2]
a_param <- exp(log_a)

# Construcción elegante de ecuación
ecuacion <- paste0("y = ", round(a_param, 4),
                    " * e^(", round(b_param, 6), "x)")

cat("La ecuacion estimada del modelo es:\n\n", "y =", round(a_param, 4), "* e^(", round(b_param, 6), "x)")
## La ecuacion estimada del modelo es:
## 
##  y = 13996.99 * e^( -0.00172 x)

9 Tabla Resumen del Modelo

tabla_resumen <- data.frame(
  Variable = c("Año de Descubrimiento", "Reservas (millones bbl)"),
  Tipo = c("Independiente (X)", "Dependiente (Y)"),
  R = c("", round(r, 2)),
  R2 = c("", round(r2, 2)),
  Parametro_a = c("", round(a_param, 4)),
  Exponente_b = c("", round(b_param, 6)),
  Ecuacion = c("", ecuacion)
)

tabla_resumen %>%
  gt() %>%
  tab_header(title = md("**Tabla N.1 del Resumen del Modelo de Regresión Exponencial**")) %>%
  cols_label(Ecuacion = "Ecuación") %>%
  tab_source_note(source_note = "Autor: Estalin — UCE") %>%
  cols_align(align = "center", columns = everything())
Tabla N.1 del Resumen del Modelo de Regresión Exponencial
Variable Tipo R R2 Parametro_a Exponente_b Ecuación
Año de Descubrimiento Independiente (X)
Reservas (millones bbl) Dependiente (Y) -0.08 0.01 13996.9867 -0.00172 y = 13996.9867 * e^(-0.00172x)
Autor: Estalin — UCE

10 Cálculo de Estimaciones

¿Qué volumen de reservas se estima para un yacimiento descubierto en 1980?

anio_test <- 1980
res_est <- a_param * exp(b_param * anio_test)

cat("Para un anio de descubrimiento de", anio_test, ", las reservas estimadas son:", round(res_est, 4), "millones de bbl")
## Para un anio de descubrimiento de 1980 , las reservas estimadas son: 464.329 millones de bbl

11 Conclusiones

# Nota: se evitan tildes y simbolos especiales (n, 2 en vez de ny R2)
# dentro de este texto generado por codigo, porque cat()/paste0() puede
# corromper esos caracteres en sesiones de R que no usan locale UTF-8
# (el resto del documento, en Markdown normal, si usa tildes sin problema).

tendencia_txt <- if (b_param < 0) "decreciente" else "creciente"
interpretacion_txt <- if (b_param < 0) {
  "lo cual es consistente con la teoria de agotamiento exploratorio: los yacimientos de mayor tamanio tienden a descubrirse en las primeras etapas de exploracion de una cuenca, mientras que los descubrimientos mas recientes corresponden, en promedio, a campos remanentes de menor volumen"
} else {
  "lo cual sugiere que, en promedio, los descubrimientos mas recientes registran mayores volumenes de reservas, posiblemente asociado a mejoras tecnologicas de exploracion y estimacion"
}

cat(paste0(
  "Entre el **Anio de Descubrimiento** y las **Reservas (millones bbl)** existe una relacion de tipo exponencial ", tendencia_txt,
  ", explicada por un coeficiente de determinacion R2 de ", round(r2 * 100, 2), "%, lo que indica que el modelo logra explicar una proporcion ",
  if (r2 < 0.3) "modesta" else if (r2 < 0.6) "moderada" else "alta",
  " de la variabilidad de las reservas en funcion del anio de descubrimiento.\n\n",
  "El coeficiente de correlacion fue de ", round(r * 100, 2), "%, evidenciando una relacion ",
  if (abs(r) < 0.3) "debil" else if (abs(r) < 0.6) "moderada" else "fuerte",
  " entre las variables analizadas, ", interpretacion_txt, ".\n\n",
  "La ecuacion matematica estimada del modelo es: y = ", round(a_param, 4), " * e^(", round(b_param, 6), "x)\n\n",
  "Cabe senialar que, al haberse imputado con la mediana una parte considerable de los valores faltantes de anio de descubrimiento y de reservas (en lugar de eliminarlos, segun lo indicado por el docente), el ajuste del modelo puede verse atenuado respecto a un analisis realizado unicamente sobre los registros originales completos. Esto se evidencia en la concentracion de observaciones alrededor de los valores de la mediana utilizada para la imputacion."
))

Entre el Anio de Descubrimiento y las Reservas (millones bbl) existe una relacion de tipo exponencial decreciente, explicada por un coeficiente de determinacion R2 de 0.66%, lo que indica que el modelo logra explicar una proporcion modesta de la variabilidad de las reservas en funcion del anio de descubrimiento.

El coeficiente de correlacion fue de -8.14%, evidenciando una relacion debil entre las variables analizadas, lo cual es consistente con la teoria de agotamiento exploratorio: los yacimientos de mayor tamanio tienden a descubrirse en las primeras etapas de exploracion de una cuenca, mientras que los descubrimientos mas recientes corresponden, en promedio, a campos remanentes de menor volumen.

La ecuacion matematica estimada del modelo es: y = 13996.9867 * e^(-0.00172x)

Cabe senialar que, al haberse imputado con la mediana una parte considerable de los valores faltantes de anio de descubrimiento y de reservas (en lugar de eliminarlos, segun lo indicado por el docente), el ajuste del modelo puede verse atenuado respecto a un analisis realizado unicamente sobre los registros originales completos. Esto se evidencia en la concentracion de observaciones alrededor de los valores de la mediana utilizada para la imputacion.