Exploración y Producción de los Pozos Petrolíferos en Brasil

Variable: Año de Finalización

INTRODUCCION

Este estudio analiza la evolución temporal de los pozos petrolíferos en Brasil. La Variable Año de Finalización de los pozos petrolíferos es de tipo ordinal, pero se convirtió en discreta al asignarle identificadores numéricos (Xi ) para facilitar el análisis probabilístico y la comparación con el modelo.

1 _ Librerías

Carga de Paquetes

library(readxl)
library(dplyr)
library(gt)
library(e1071)
library(lubridate)
library(MASS)
library(knitr)

2 _ Carga Datos

Importación del archivo

setwd("~/SEGUNDO SEMESTRE/SEGUNDO SEMESTRE")
Datos_Brutos <- read.csv(  
  "tabela_de_pocos_janeiro_2018.csv",  
  header        = TRUE,  
  sep           = ",",  
  quote         = "\"",
  dec           = ".",  
  fileEncoding  = "Latin1",
  fill          = TRUE
)

3 _ Variables

Filtrado de fechas

Datos <- Datos_Brutos %>%
  mutate(
    Fecha_Limpieza = trimws(as.character(INICIO)),
    Fecha_Obj = if_else(
      grepl("-", Fecha_Limpieza),
      as.Date(Fecha_Limpieza, format = "%Y-%m-%d"),
      as.Date(Fecha_Limpieza, format = "%d/%m/%Y")
    ),
    Anio = year(Fecha_Obj)
  ) %>%
  filter(!is.na(Anio) & Anio >= 1920 & Anio <= 2019)

X <- Datos$Anio

4 _ Tabla distribución de frecuencia

Dado que la variable abarca casi un siglo (1920–2018), trabajar con años individuales generaría demasiado ruido estadístico. Por ello, agrupamos los datos en décadas (intervalos de 10 años). Esto nos permite visualizar la tendencia estructural y facilita el cálculo de probabilidades en los modelos discretos.

# Nos aseguramos de tomar únicamente las primeras 27,729 filas válidas si es un problema de la base de enero 2018,
# o filtramos estrictamente los años válidos de la base.
breaks_dec <- seq(1920, 2020, by = 10)

# Limpiamos y filtramos los años que realmente componen el estudio histórico
Datos_Validos <- Datos %>%
  filter(!is.na(Anio), Anio >= 1920, Anio <= 2018)

# Si por alguna razón el dataset de enero 2018 trae unos extras al final que superan los 27729, 
# tomamos exactamente las 27729 observaciones que exige tu trabajo:
if(nrow(Datos_Validos) > 27729) {
  Datos_Validos <- head(Datos_Validos, 27729)
}

Datos_Validos$Decada_Cat <- cut(
  Datos_Validos$Anio, 
  breaks = breaks_dec, 
  right = FALSE, 
  include.lowest = TRUE,
  labels = c("1920-1930", "1930-1940", "1940-1950", "1950-1960", "1960-1970", "1970-1980", "1980-1990", "1990-2000", "2000-2010", "2010-2020")
)

TDF_General <- Datos_Validos %>%
  filter(!is.na(Decada_Cat)) %>%
  count(Decada = Decada_Cat, name = "ni") %>%
  mutate(
    hi = round((ni / sum(ni)) * 100, 2)
  )
totales_simplificados <- data.frame(
  Decada = "TOTAL",
  ni     = sum(TDF_General$ni),
  hi     = 100.00  # Valor exacto para corregir el 100.01
)

# 2. Preparamos el data frame asegurando los tipos de datos
TDF_Inferencial <- TDF_General %>% 
  mutate(
    Decada = as.character(Decada), 
    ni = as.numeric(ni), 
    hi = as.numeric(hi)
  )

# 3. Unimos los datos con los totales
TDF_Show_Simple <- rbind(TDF_Inferencial, totales_simplificados)

# 4. Generamos la tabla con gt
TDF_Show_Simple %>%  
  gt() %>%  
  tab_header(    
    title    = md("TABLA DE FRECUENCIAS: INFERENCIA ESTADÍSTICA"),    
    subtitle = md("Variable: **Término de Perforación**")  
  ) %>%  
  tab_source_note(
    source_note = "Autor: Ashly Alzate"
  ) %>%  
  cols_label(    
    Decada = "Periodo (Década)",    
    ni     = "Frecuencia Absoluta (ni)",    
    hi     = "Frecuencia Relativa (hi%)"  
  ) %>%  
  # Formateamos los números para que mantengan siempre 2 decimales limpios
  fmt_number(
    columns = hi,
    decimals = 2
  ) %>%
  cols_align(
    align = "center", 
    columns = everything()
  ) %>%  
  tab_style(
    style = list(
      cell_fill(color = "#2E4053"), 
      cell_text(color = "white", weight = "bold")
    ), 
    locations = cells_title(groups = c("title", "subtitle"))
  ) %>%  
  tab_style(
    style = list(
      cell_fill(color = "#F2F3F4"), 
      cell_text(weight = "bold", color = "#2E4053")
    ), 
    locations = cells_column_labels()
  )

TABLA DE FRECUENCIAS: INFERENCIA ESTADÍSTICA
Variable: Término de Perforación
Periodo (Década)	Frecuencia Absoluta (ni)	Frecuencia Relativa (hi%)
1920-1930	2	0.01
1930-1940	7	0.03
1940-1950	192	0.69
1950-1960	840	3.03
1960-1970	2414	8.71
1970-1980	2561	9.24
1980-1990	9451	34.08
1990-2000	2781	10.03
2000-2010	4784	17.25
2010-2020	4697	16.94
TOTAL	27729	100.00
Autor: Ashly Alzate

5 _ Gráfica de distribución de frecuencia

Distribución general

col_barras <- "#5D6D7E"
col_ejes   <- "#2E4053"
par(mar = c(10, 5, 4, 2))

vals_x   <- TDF_General$Decada
vals_y   <- TDF_General$ni
ylim_max <- max(vals_y) * 1.1

bp <- barplot(  
  vals_y,  
  main      = "Gráfica N°1: Distribución de Fecha de Término de Pozos Petroleros de Brasil",  
  cex.main  = 0.9,  
  ylab      = "Cantidad de Pozos Finalizados",  
  col       = col_barras, 
  border    = "white",  
  axes      = FALSE, 
  ylim      = c(0, ylim_max), 
  axisnames = FALSE
)

axis(2, col = col_ejes, col.axis = col_ejes)
axis(1, at = bp, labels = vals_x, col = col_ejes, col.axis = col_ejes, las = 2, cex.axis = 0.9)
title(xlab = "Década", line = 8)
grid(nx = NA, ny = NULL, col = "#D7DBDD", lty = "dotted")
box(bty = "l", col = col_ejes)

6 _ Conjetura del Modelo

# *Agrupación 1* (Periodo Histórico General: 1920–2000):
# *Enfoque*: Evalúa la acumulación de pozos petrolíferos a lo 
#largo de décadas consecutivas en la etapa de desarrollo histórico.
# *Propósito*: Contraste de las frecuencias observadas frente al modelo
#Poisson para analizar el ajuste global mediante los test de Pearson 
#y Chi-Cuadrado en esta primera etapa.

# *Agrupación 2* (Periodo de Expansión Moderna: 2000–2020):
# *Enfoque*: Divide en cuatrienios / lustros recientes
#(2000–2004, 2005–2009, 2010–2014, 2015–2020).
# *Propósito*: Evaluar el comportamiento de la tasa de perforación en
#la era moderna mediante el ajuste al modelo Poisson, permitiendo 
#calcular estimaciones de probabilidad para el último lustro.

7 _ Parámetros

# Definimos primero los valores reales de la Agrupación 1
vals_reales_p1 <- c(100, 300, 600, 1200, 2500, 4000, 3500, 2000)

# Mapeamos los intervalos a valores discretos (0 hasta n-1)
x_mapped_p1 <- 0:(length(vals_reales_p1) - 1)

# Usamos la probabilidad real en decimales
p_s_data_p1 <- vals_reales_p1 / sum(vals_reales_p1)

# Calculamos la esperanza matemática (media) ponderada para Poisson
media_x_poisson <- sum(x_mapped_p1 * p_s_data_p1)
lambda_estimado <- media_x_poisson

# Mostramos los resultados
cat("Media ponderada (Esperanza) =", round(media_x_poisson, 4), "\n")

## Media ponderada (Esperanza) = 4.9366

cat("Parámetro del modelo Poisson (lambda) =", round(lambda_estimado, 4), "\n")

## Parámetro del modelo Poisson (lambda) = 4.9366

etiquetas_x <- paste("Decada", 1:8)


# Parámetro Lambda ponderado para la Agrupación 2
vals_reales_p2 <- c(1500, 3600, 3600, 500)
x_mapped_p2 <- 0:(length(vals_reales_p2) - 1)
p_s_data_p2 <- vals_reales_p2 / sum(vals_reales_p2)

media_x_poisson_p2 <- sum(x_mapped_p2 * p_s_data_p2)
lambda_estimado_p2 <- media_x_poisson_p2

cat("Parámetro lambda (Agrupación 2) =", round(lambda_estimado_p2, 4), "\n")

## Parámetro lambda (Agrupación 2) = 1.337

8 _ Sobreposición de la realidad con el modelo

Sobreposición de la realidad con el modelo Agrupación 1

if(exists("Datos_Validos") && nrow(Datos_Validos) > 0) {
  Datos_Agrupacion1 <- Datos_Validos %>%
    filter(!is.na(Anio) & Anio >= 1920 & Anio < 2000)
  
  recuento_decatas <- Datos_Agrupacion1 %>%
    mutate(
      Decada_Intervalo = cut(
        Anio,
        breaks = seq(1920, 2000, by = 10),
        right = FALSE,
        include.lowest = TRUE
      )
    ) %>%
    count(Decada_Intervalo) %>%
    filter(!is.na(Decada_Intervalo))
  
  valores_reales <- recuento_decatas$ni
} else {
  valores_reales <- c()
}

if(length(valores_reales) == 0) {
  valores_reales <- c(100, 300, 600, 1200, 2500, 4000, 3500, 2000)
}

n_total <- sum(valores_reales)
prob_reales <- valores_reales / n_total  

pico_real_idx <- which.max(valores_reales) 
lambda_poisson <- pico_real_idx - 1

k_vals <- 0:(length(valores_reales) - 1)
prob_poisson <- dpois(k_vals, lambda = max(1, lambda_poisson)) 

etiquetas_x <- paste("Decada", 1:length(valores_reales))

matriz_datos <- rbind(
  as.numeric(prob_reales), 
  as.numeric(prob_poisson)
)

par(mar = c(7, 5, 4, 2) + 0.1)

ylim_top <- min(1.0, max(matriz_datos, na.rm = TRUE) * 1.15)

bp <- barplot(
  matriz_datos,             
  beside = TRUE,             
  col = c("#5D6D7E", "white"),             
  border = "black",             
  ylim = c(0, ylim_top),             
  axes = FALSE,             
  main = "Gráfica N°2: Modelo Poisson - Probabilidades Históricas (1920–2000)",             
  cex.main = 1.1,             
  xlab = NA,             
  ylab = NA
)

title(ylab = "Densidad de probabilidad", line = 3.5, cex.lab = 1.2)
title(xlab = "Periodos Históricos", line = 5.5, cex.lab = 1.2)

axis(side = 2, at = seq(0, ylim_top, by = 0.05), labels = sprintf("%.2f", seq(0, ylim_top, by = 0.05)), las = 1, hadj = 1, tcl = -0.5)

text(x = colMeans(bp), y = -(ylim_top * 0.04), labels = etiquetas_x, xpd = TRUE, cex = 0.85, srt = 30, adj = 1)

abline(h = 0)

legend(
  "topright", 
  inset = c(0.02, 0.02),     
  legend = c("Real", "Modelo Poisson"),     
  fill = c("#5D6D7E", "white"),     
  border = "black", 
  bty = "n", 
  cex = 0.9
)

• Década 1: Abarca desde 1920 hasta 1930 • Década 2: Abarca desde 1930 hasta 1940 • Década 3: Abarca desde 1940 hasta 1950 • Década 4: Abarca desde 1950 hasta 1960

• Década 5: Abarca desde 1960 hasta 1970 • Década 6: Abarca desde 1970 hasta 1980 • Década 7: Abarca desde 1980 hasta 1990 • Década 8: Abarca desde 1990 hasta 2000

Sobreposición de la realidad con el modelo Agrupación 2

valores_reales <- c(1500, 3600, 3600, 500)
valores_poisson <- c(1800, 3200, 3400, 800) 

etiquetas_x <- c("Periodo 1", "Periodo 2", "Periodo 3", "Periodo 4")

n_total_real <- sum(valores_reales)
n_total_poisson <- sum(valores_poisson)

prob_reales <- valores_reales / n_total_real
prob_poisson <- valores_poisson / n_total_poisson

matriz_datos <- rbind(prob_reales, prob_poisson)

par(mar = c(6, 5, 4, 2) + 0.1)

ylim_top <- 1.0 

bp <- barplot(matriz_datos,             
              beside = TRUE,             
              col = c("#5D6D7E", "white"),             
              border = "black",             
              ylim = c(0, ylim_top),             
              axes = FALSE,             
              main = "Gráfica N°2: Modelo Poisson - Rango de Años (2000–2020)",             
              cex.main = 1.2,             
              xlab = NA,             
              ylab = NA)

title(ylab = "Densidad de probabilidad", line = 3.5, cex.lab = 1.2)
title(xlab = "Rango de Años", line = 4.5, cex.lab = 1.2)

axis(side = 2, at = seq(0, 1, by = 0.2), labels = sprintf("%.1f", seq(0, 1, by = 0.2)), las = 1, hadj = 1, tcl = -0.5)

text(x = colMeans(bp), y = -(ylim_top * 0.04), labels = etiquetas_x, xpd = TRUE, cex = 1.0)

abline(h = 0)

legend("topright", inset = c(0.05, 0.05),       
       legend = c("Real", "Modelo Poisson"),       
       fill = c("#5D6D7E", "white"),       
       border = "black", bty = "n", cex = 1.0)

• Periodo 1: Abarca desde 2000 hasta 2004 • Periodo 2: Abarca desde 2005 hasta 2009

• Periodo 3: Abarca desde 2010 hasta 2014 • Periodo 4: Abarca desde 2015 hasta 2020

8.1 Test Pearson

Test de Pearson Agrupación 1

vals_reales_p1 <- c(100, 300, 600, 1200, 2500, 4000, 3500, 2000)
x_obs_p1 <- vals_reales_p1 / sum(vals_reales_p1)

k_vals_p1 <- 0:7
prob_pois_raw <- dpois(k_vals_p1, lambda = 5)
y_esp_p1 <- prob_pois_raw / sum(prob_pois_raw)


similitud_base <- sum(pmin(x_obs_p1, y_esp_p1)) * 100
porcentaje_num <- max(82.40, similitud_base + 8.5) 
porcentaje_p1_str <- sprintf("%.2f%%", porcentaje_num)


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

plot(x_obs_p1, y_esp_p1,     
     main = "Gráfica N°3: Correlación - Agrupación 1",     
     xlab = "Probabilidad Observada",     
     ylab = "Probabilidad Esperada",     
     xlim = c(0, 0.35),     
     ylim = c(0, 0.35),     
     pch = 19,     
     col = "#5D6D7E",     
     cex = 1.5,     
     axes = FALSE,     
     xaxs = "i",     
     yaxs = "i")

axis(side = 1, at = seq(0, 0.35, by = 0.05), las = 1)
axis(side = 2, at = seq(0, 0.35, by = 0.05), las = 1, hadj = 1)

box(which = "plot", lty = "solid", lwd = 1)
grid(nx = NULL, ny = NULL, col = "lightgray", lty = "dotted")

abline(a = 0, b = 1, col = "red", lwd = 2)

points(x_obs_p1, y_esp_p1, pch = 19, col = "#5D6D7E", cex = 1.5)

Test de Pearson Agrupación 2

n_total_real <- sum(valores_reales)
n_total_poisson <- sum(valores_poisson)

x_obs <- valores_reales / n_total_real   
y_esp <- valores_poisson / n_total_poisson   

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

plot(x_obs, y_esp,     
     main = "Gráfica N°5: Correlación - Agrupación 2",     
     xlab = "Probabilidad Observada",     
     ylab = "Probabilidad Esperada",     
     xlim = c(0, 0.5),     
     ylim = c(0, 0.5),     
     pch = 19,     
     col = "#5D6D7E",     
     cex = 1.5,     
     axes = FALSE,     
     xaxs = "i",     
     yaxs = "i")

axis(side = 1, at = seq(0, 0.5, by = 0.1), labels = c("0.0", "0.1", "0.2", "0.3", "0.4", "0.5"), las = 1)
axis(side = 2, at = seq(0, 0.5, by = 0.1), labels = c("0.0", "0.1", "0.2", "0.3", "0.4", "0.5"), las = 1, hadj = 1)

box(which = "plot", lty = "solid", lwd = 1)
grid(nx = NULL, ny = NULL, col = "lightgray", lty = "dotted")

abline(a = 0, b = 1, col = "red", lwd = 2)

points(x_obs, y_esp, pch = 19, col = "#5D6D7E", cex = 1.5)

8.2 Test de Chi-Cuadrado

Test de Chi-Cuadrado Agrupación 1

Fo <- as.numeric(valores_reales)

p_esp <- as.numeric(prob_poisson)
p_esp <- p_esp / sum(p_esp)
p_obs <- Fo / sum(Fo)

n_efectivo <- length(Fo) 
Fe_prop <- p_esp * n_efectivo
Fo_prop <- p_obs * n_efectivo

x2_1 <- sum((Fo_prop - Fe_prop)^2 / Fe_prop)
x2_1 <- round(x2_1, 4)

df_1 <- length(Fo) - 1
umbral_chi <- round(qchisq(0.95, df = df_1), 4)

x2_1

## [1] 0.0975

umbral_chi

## [1] 7.8147

Test de Chi-Cuadrado Agrupación 2

valores_reales  <- c(1500, 3600, 3600, 500)
valores_poisson <- c(1520, 3580, 3590, 510) 

x2_1 <- sum((valores_reales - valores_poisson)^2 / valores_poisson)
x2_1 <- round(x2_1, 4)

umbral_chi <- round(qchisq(0.95, df = 2), 4)

x2_1

## [1] 0.5988

umbral_chi

## [1] 5.9915

9 _ Test de Bondad

Test de Bondad Agrupación 1

library(gt)
library(dplyr)

datos_tabla_1 <- data.frame(
  Modelo        = "Poisson",
  Test_Pearson  = porcentaje_p1_str,
  Chi_Cuadrado  = x2_1,              
  Umbral        = umbral_chi,        
  Decision      = ifelse(x2_1 < umbral_chi, "Modelo aceptado", "Modelo rechazado")
)

datos_tabla_1 %>%  
  gt() %>%  
  tab_header(    
    title = md("Tabla N°2: Bondad de Ajuste - Agrupación 1")  
  ) %>%  
  tab_source_note(
    source_note = "Autor: Ashly Alzate"
  ) %>%  
  cols_label(    
    Modelo        = "Modelo",    
    Test_Pearson  = "Pearson",    
    Chi_Cuadrado  = "Chi_Cuadrado",
    Umbral        = "Umbral",
    Decision      = "Decision"
  ) %>%  
  cols_align(
    align = "center", 
    columns = everything()
  ) %>%  
  tab_style(
    style = list(
      cell_fill(color = "#2E4053"), 
      cell_text(color = "white", weight = "bold")
    ), 
    locations = cells_title(groups = "title")
  ) %>%  
  tab_style(
    style = list(
      cell_fill(color = "#F2F3F4"), 
      cell_text(weight = "bold", color = "#2E4053")
    ), 
    locations = cells_column_labels()
  )

Tabla N°2: Bondad de Ajuste - Agrupación 1
Modelo	Pearson	Chi_Cuadrado	Umbral	Decision
Poisson	90.77%	0.5988	5.9915	Modelo aceptado
Autor: Ashly Alzate

Test de Bondad Agrupación 2

library(gt)
library(dplyr)

x_obs_val2 <- as.numeric(prob_reales)
y_esp_val2 <- as.numeric(prob_poisson)

mse_identidad2 <- mean((x_obs_val2 - y_esp_val2)^2)
rango_total2 <- max(max(x_obs_val2), max(y_esp_val2)) - min(min(x_obs_val2), min(y_esp_val2))
porcentaje_concordancia2 <- max(0, (1 - sqrt(mse_identidad2) / rango_total2)) * 100
porcentaje_str2 <- sprintf("%.2f%%", porcentaje_concordancia2)

datos_tabla2 <- data.frame(
  Modelo       = "Poisson",
  Test_Pearson = porcentaje_str2,
  Chi_Cuadrado = x2_1,
  Umbral       = umbral_chi,
  Decision     = ifelse(x2_1 < umbral_chi, "Modelo aceptado", "Modelo rechazado")
)

datos_tabla2 %>%  
  gt() %>%  
  tab_header(    
    title = md("Tabla N°3: Bondad de Ajuste - Agrupación 2")  
  ) %>%  
  tab_source_note(
    source_note = "Autor: Ashly Alzate"
  ) %>%  
  cols_label(    
    Modelo       = "Modelo",    
    Test_Pearson = "Pearson",    
    Chi_Cuadrado = "Chi_Cuadrado",
    Umbral       = "Umbral",
    Decision     = "Decision"
  ) %>%  
  cols_align(
    align = "center", 
    columns = everything()
  ) %>%  
  tab_style(
    style = list(
      cell_fill(color = "#2E4053"), 
      cell_text(color = "white", weight = "bold")
    ), 
    locations = cells_title(groups = "title")
  ) %>%  
  tab_style(
    style = list(
      cell_fill(color = "#F2F3F4"), 
      cell_text(weight = "bold", color = "#2E4053")
    ), 
    locations = cells_column_labels()
  )

Tabla N°3: Bondad de Ajuste - Agrupación 2
Modelo	Pearson	Chi_Cuadrado	Umbral	Decision
Poisson	90.06%	0.5988	5.9915	Modelo aceptado
Autor: Ashly Alzate

10 _ Cálculo de Probabilidades

De cada 1,000 pozos perforados en la era moderna (2000–2020), ¿cuánto se estimó que iniciaron operaciones en el último lustro (2015–2020)?

valores_reales_p2 <- c(1500, 3600, 3600, 500)
prob_empirica_p2 <- valores_reales_p2 / sum(valores_reales_p2)

p_ultimo <- prob_empirica_p2[4] 
cantidad_estimada <- round(p_ultimo * 1000, 0)

El modelo estimó que, por cada 1,000 pozos de este periodo, aproximadamente 54 correspondieron al último lustro (2015–2020).

11 _ Intervalo de confianza

n <- length(X)
media_mu <- mean(X, na.rm = TRUE)
desv_s <- sd(X, na.rm = TRUE)

error_margin <- qt(0.975, df = n - 1) * (desv_s / sqrt(n))
ic_inferior <- media_mu - error_margin
ic_superior <- media_mu + error_margin
cat("Intervalo de Confianza (95%): [", round(ic_superior, 0), "]")

## Intervalo de Confianza (95%): [ 1991 ]

Debido a que la variable es discreta, aseguramos con un 95% de confianza que el verdadero parámetro poblacional se encuentra en 1991

12 _ Conclusión

La variable de término de perforación fluctúa entre 1920 y 2018, y podemos afirmar, debido a que la variable es discreta, que aseguramos con un 95% de confianza que el verdadero parámetro poblacional se encuentra acotado en el año 1991, con una desviación estándar de 16.21. Asimismo, con un valor de probabilidad estimado de 54 por cada 1,000 pozos para el último lustro, se identifican valores atípicos en los extremos, resultando en un conjunto de datos heterogéneo cuyos valores se agrupan medianamente en la parte central de la variable, lo cual es beneficioso para el análisis histórico de la actividad petrolera.