Actividad 1

Importar paquetes y llamar librerias

#install.packages("readxl")
library(readxl)
#install.packages("ggplot2")
library(ggplot2)
#install.packages("dplyr")
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

#install.packages("tidyr")
library(tidyr)
#install.packages("plm")
library(plm)

## 
## Attaching package: 'plm'

## The following objects are masked from 'package:dplyr':
## 
##     between, lag, lead

#install.packages("gplots")
library(gplots)

## 
## ---------------------
## gplots 3.3.0 loaded:
##   * Use citation('gplots') for citation info.
##   * Homepage: https://talgalili.github.io/gplots/
##   * Report issues: https://github.com/talgalili/gplots/issues
##   * Ask questions: https://stackoverflow.com/questions/tagged/gplots
##   * Suppress this message with: suppressPackageStartupMessages(library(gplots))
## ---------------------

## 
## Attaching package: 'gplots'

## The following object is masked from 'package:stats':
## 
##     lowess

Parte 1. Patentes

Este archivo es un panel de datos empresa-año: cada fila es una combinación de una empresa (cusip) y un año (2012–2021). Es un tipo de dataset muy usado en econometría de innovación (gasto en I+D vs. patentes), popularizado por estudios como los de Hall, Griliches y Hausman sobre patentes y R&D. Diccionario de variables: cusip: identificador de la empresa (no es una variable predictiva, es un ID). merger: 1 si la empresa tuvo una fusión importante ese año, 0 si no. employ: empleados, en miles. return: retorno de la acción, en %. patents: patentes solicitadas en el año. patentsg: patentes concedidas en el año — esta es nuestra variable objetivo. stckpr: precio de la acción. rnd: gasto en I+D, en millones de dólares corrientes. rndeflt: gasto en I+D, en millones de dólares deflactados (precios constantes de 1972). rndstck: “stock” acumulado de I+D (inversión acumulada, no solo del año). sales: ventas, en millones de dólares corrientes. sic: código de industria a 4 dígitos. year: año de la observación.

# Cargamos la base de datos
df1 <- read_xls("C:\\Users\\me\\Desktop\\reto\\bds\\PATENT 3.xls")

# Relación con variables explicativas
vars_num <- df1 %>% select(employ, return, patents, stckpr, rnd, rndeflt, rndstck, sales, patentsg)
round(cor(vars_num, use = "complete.obs"), 2)

##          employ return patents stckpr  rnd rndeflt rndstck sales patentsg
## employ     1.00   0.04    0.71   0.38 0.86    0.90    0.85  0.81     0.74
## return     0.04   1.00    0.08   0.31 0.08    0.09    0.07  0.05     0.09
## patents    0.71   0.08    1.00   0.45 0.55    0.61    0.55  0.50     0.94
## stckpr     0.38   0.31    0.45   1.00 0.40    0.45    0.39  0.33     0.47
## rnd        0.86   0.08    0.55   0.40 1.00    0.98    0.99  0.79     0.59
## rndeflt    0.90   0.09    0.61   0.45 0.98    1.00    0.96  0.78     0.64
## rndstck    0.85   0.07    0.55   0.39 0.99    0.96    1.00  0.78     0.59
## sales      0.81   0.05    0.50   0.33 0.79    0.78    0.78  1.00     0.55
## patentsg   0.74   0.09    0.94   0.47 0.59    0.64    0.59  0.55     1.00

ggplot(df1, aes(x = rnd, y = patentsg)) +
  geom_point(alpha = 0.4, color = "#4C72B0") +
  labs(title = "Gasto en I+D vs. patentes concedidas",
       x = "Gasto en I+D (millones $ corrientes)", y = "Patentes concedidas") +
  theme_minimal()

# Convertimos nuestra bdd a datos panel
df_modelo <- df1 %>%
  select(cusip, year, patentsg, rnd, employ, sales, patents, merger)

pdf <- pdata.frame(df_modelo, index = c("cusip", "year"))

Variables explicativas elegidas: ** rnd (gasto en I+D), employ (empleados), sales (ventas), patents (patentes solicitadas ese año) y merger (si hubo fusión). La idea es explicar cuántas patentes le conceden a una empresa (patentsg) a partir de cuánto invierte, qué tan grande es, y cuánto solicitó. A continuación se presenta un modelo de regresión multiple para tener algo con qué comparar nuestros siguientes modelos.

modelo_referencia <- lm(patentsg ~ rnd + employ + sales + patents + merger, data = pdf)
summary(modelo_referencia)

## 
## Call:
## lm(formula = patentsg ~ rnd + employ + sales + patents + merger, 
##     data = pdf)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -266.57   -3.55   -1.49    0.79  644.42 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.3832054  0.5906000   2.342  0.01927 *  
## rnd          0.0265024  0.0090625   2.924  0.00349 ** 
## employ       0.1359658  0.0282759   4.809 1.62e-06 ***
## sales        0.0003127  0.0002654   1.178  0.23884    
## patents      0.9653663  0.0112721  85.642  < 2e-16 ***
## merger      -2.4621552  4.0822979  -0.603  0.54648    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 25.5 on 2230 degrees of freedom
##   (24 observations deleted due to missingness)
## Multiple R-squared:  0.8976, Adjusted R-squared:  0.8973 
## F-statistic:  3908 on 5 and 2230 DF,  p-value: < 2.2e-16

Este modelo muestra una R cuadrada de .89, lo que significa que las variables escogidas explican el 89 del comportamiento de la variable dependiente. Pasamos a los modelos Pooling, Within y Random

form <- patentsg ~ rnd + employ + sales + patents + merger

pooling <- plm(form, data = pdf, model = "pooling")
within  <- plm(form, data = pdf, model = "within")
random  <- plm(form, data = pdf, model = "random")

Resumen de Pooling

summary(pooling)

## Pooling Model
## 
## Call:
## plm(formula = form, data = pdf, model = "pooling")
## 
## Unbalanced Panel: n = 225, T = 8-10, N = 2236
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -266.56790   -3.54922   -1.49325    0.78669  644.42008 
## 
## Coefficients:
##                Estimate  Std. Error t-value  Pr(>|t|)    
## (Intercept)  1.38320543  0.59059998  2.3420  0.019267 *  
## rnd          0.02650238  0.00906254  2.9244  0.003486 ** 
## employ       0.13596581  0.02827591  4.8085 1.622e-06 ***
## sales        0.00031272  0.00026542  1.1782  0.238838    
## patents      0.96536633  0.01127212 85.6420 < 2.2e-16 ***
## merger      -2.46215521  4.08229788 -0.6031  0.546484    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    14152000
## Residual Sum of Squares: 1449500
## R-Squared:      0.89758
## Adj. R-Squared: 0.89735
## F-statistic: 3908.43 on 5 and 2230 DF, p-value: < 2.22e-16

Resumen de Within / Efectos Fijos

summary(within)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = form, data = pdf, model = "within")
## 
## Unbalanced Panel: n = 225, T = 8-10, N = 2236
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -223.29533   -1.88554   -0.31819    1.51373  262.80141 
## 
## Coefficients:
##           Estimate Std. Error  t-value  Pr(>|t|)    
## rnd     -0.1309985  0.0125503 -10.4379 < 2.2e-16 ***
## employ  -0.0663792  0.0607919  -1.0919 0.2750038    
## sales   -0.0013453  0.0003577  -3.7609 0.0001742 ***
## patents  0.0670273  0.0188908   3.5481 0.0003969 ***
## merger   1.8352766  3.5413607   0.5182 0.6043476    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    715600
## Residual Sum of Squares: 586560
## R-Squared:      0.18033
## Adj. R-Squared: 0.086759
## F-statistic: 88.2654 on 5 and 2006 DF, p-value: < 2.22e-16

Resumen de Random / Efectos Aleatorios

summary(random)

## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = form, data = pdf, model = "random")
## 
## Unbalanced Panel: n = 225, T = 8-10, N = 2236
## 
## Effects:
##                 var std.dev share
## idiosyncratic 292.4    17.1     1
## individual      0.0     0.0     0
## theta:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       0       0       0       0       0 
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -266.56790   -3.54922   -1.49325    0.78669  644.42008 
## 
## Coefficients:
##                Estimate  Std. Error z-value  Pr(>|z|)    
## (Intercept)  1.38320543  0.59059998  2.3420  0.019179 *  
## rnd          0.02650238  0.00906254  2.9244  0.003451 ** 
## employ       0.13596581  0.02827591  4.8085  1.52e-06 ***
## sales        0.00031272  0.00026542  1.1782  0.238712    
## patents      0.96536633  0.01127212 85.6420 < 2.2e-16 ***
## merger      -2.46215521  4.08229788 -0.6031  0.546422    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    14152000
## Residual Sum of Squares: 1449500
## R-Squared:      0.89758
## Adj. R-Squared: 0.89735
## Chisq: 19542.1 on 5 DF, p-value: < 2.22e-16

Pruebas de especificación entre los diferentes modelos

pFtest(within, pooling)

## 
##  F test for individual effects
## 
## data:  form
## F = 13.175, df1 = 224, df2 = 2006, p-value < 2.2e-16
## alternative hypothesis: significant effects

efectos_fijos <- fixef(within)
head(efectos_fijos)

##        800       4626       4671       7500       7603      20753 
## 41.2082685  5.9277947  1.4333872  0.8330017  0.5298820  2.0559662

Pasamos a la prueba de pronóstico

escenario_2022 <- df_modelo %>%
  filter(year == 2021, complete.cases(across(c(rnd, employ, sales, patents, merger))))

fe_df <- data.frame(cusip = names(efectos_fijos), efecto_fijo = as.numeric(efectos_fijos))
# pdata.frame guarda cusip como texto; igualamos el tipo para poder unir
fe_df$cusip <- as.numeric(as.character(fe_df$cusip))

coefs <- coef(within)

pronostico_2022 <- escenario_2022 %>%
  inner_join(fe_df, by = "cusip") %>%
  mutate(
    patentsg_pronostico_2022 =
      efecto_fijo +
      coefs["rnd"]     * rnd +
      coefs["employ"]  * employ +
      coefs["sales"]   * sales +
      coefs["patents"] * patents +
      coefs["merger"]  * merger
  ) %>%
  select(cusip, patentsg_2021 = patentsg, patentsg_pronostico_2022)

cat("Empresas con pronóstico calculado:", nrow(pronostico_2022), "de", n_distinct(df_modelo$cusip), "\n")

## Empresas con pronóstico calculado: 223 de 226

knitr::kable(head(pronostico_2022, 15), digits = 1)

cusip	patentsg_2021	patentsg_pronostico_2022
800	70	38.4
4626	7	4.5
4671	1	1.2
7500	1	0.7
7603	0	0.2
20753	0	1.6
21367	0	1.0
23519	6	14.5
29069	0	0.8
38213	2	1.9
54303	2	0.9
67131	0	2.7
67383	9	3.5
74077	18	19.5
77491	0	0.3

summary(pronostico_2022$patentsg_pronostico_2022)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -13.6016   0.8866   3.2913  22.7981  19.5489 772.9755

cat("Pronósticos negativos (sin sentido para un conteo):", sum(pronostico_2022$patentsg_pronostico_2022 < 0), "\n")

## Pronósticos negativos (sin sentido para un conteo): 8

ggplot(pronostico_2022, aes(x = patentsg_2021, y = patentsg_pronostico_2022)) +
  geom_point(alpha = 0.5, color = "#4C72B0") +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey40") +
  labs(title = "Patentes concedidas: 2021 real vs. 2022 pronosticado",
       subtitle = "La línea punteada marca 'sin cambio' (pronóstico = valor de 2021)",
       x = "Patentes concedidas 2021 (real)", y = "Patentes concedidas 2022 (pronóstico)") +
  theme_minimal()

Conclusiones finales para el ejercicio 1: Las pruebas de especificación (PFtest, phtest) coinciden en que efectos fijos es el modelo correcto para estos datos ya que sí hay características propias de cada empresa que importan, y están correlacionadas con las variables explicativas. El pronóstico para 2022 se hizo bajo el supuesto de “todo sigue igual que 2021” # Parte 2. Cuidado de la Piel Si tuvieras que invertir en alguna sub-categoría, ¿en cuál lo harías? Justifica ampliamente tu respuesta. 1. Arreglos visuales

COLOR_HIST   <- "#2a78d6"  # azul - serie 1 (histórico / lo observado)
COLOR_FCST   <- "#eb6834"  # naranja - serie 2 (pronóstico)
COLOR_POS    <- "#2a78d6"  # azul - variación positiva
COLOR_NEG    <- "#e34948"  # rojo - variación negativa
COLOR_NEUTRAL<- "#898781"

2. Cargamos la Base de Datos

df2 <- read.csv("C:\\Users\\me\\Desktop\\reto\\bds\\market_sizes_raw.csv")

3. Limpieza de Datos

# col_types = "text" fuerza a que TODAS las columnas se lean como texto desde
# el inicio (en vez de dejar que readxl adivine el tipo columna por columna).
# Esto evita el error "Can't combine `2011` <character> and `2016` <double>",
# que pasa porque algunas columnas de año traen "-" en otras filas del archivo
# (en Prestige/Dermocosmetics) y eso hace que readxl las adivine como texto,
# mientras que columnas sin "-" las adivina como número — con tipos distintos,
# pivot_longer() no las puede apilar en una sola columna.
raw <- read_excel("C:\\Users\\me\\Desktop\\reto\\bds\\Market_sizes.xlsx", sheet = "Statistics Data", skip = 5, col_types = "text")
names(raw)[1] <- "Geography"
raw <- raw %>% filter(!is.na(Geography), Geography == "Mexico")

categorias_producto <- c("Bath and Shower", "Deodorants", "Depilatories", "Fragrances",
                          "Hair Care", "Men's Grooming", "Skin Care", "Sun Care")

raw <- raw %>% filter(Category %in% categorias_producto)

datos <- raw %>%
  select(Category, `2011`:`2025`) %>%
  pivot_longer(-Category, names_to = "year", values_to = "value") %>%
  mutate(year = as.integer(year), value = suppressWarnings(as.numeric(value)))

knitr::kable(datos %>% filter(year %in% c(2011, 2015, 2019, 2020, 2021, 2025)) %>%
               pivot_wider(names_from = year, values_from = value),
             digits = 0, caption = "Valores de mercado (MXN millones), años seleccionados")

Valores de mercado (MXN millones), años seleccionados

Category	2011	2015	2019	2020	2021	2025
Bath and Shower	8410	10815	12960	14634	15469	21342
Deodorants	9151	12379	16745	13835	15091	21824
Depilatories	802	1106	1516	1491	1564	1873
Fragrances	18728	22568	28483	25516	31801	56742
Hair Care	24854	30416	37247	36951	39555	56290
Men’s Grooming	18672	24774	32693	29669	32939	49508
Skin Care	25650	30964	40638	42196	47981	71830
Sun Care	1145	1693	2444	2043	2294	4326

4. Exploración: ¿cómo se ha movido cada categoría?

ggplot(datos, aes(x = year, y = value)) +
  geom_line(color = COLOR_HIST, linewidth = 0.9) +
  geom_point(color = COLOR_HIST, size = 1.6) +
  facet_wrap(~ Category, scales = "free_y", ncol = 4) +
  labs(title = "Tamaño de mercado histórico por subcategoría (2011-2025)",
       x = "Año", y = "MXN millones") +
  theme_minimal(base_size = 10) +
  theme(strip.text = element_text(face = "bold"))

# 5.Crecimiento acumulado 2011 → 2025 (CAGR)

cagr <- datos %>%
  group_by(Category) %>%
  summarise(
    valor_2011 = value[year == 2011],
    valor_2025 = value[year == 2025],
    cagr = (valor_2025 / valor_2011)^(1 / (2025 - 2011)) - 1
  ) %>%
  arrange(desc(cagr))

knitr::kable(cagr, digits = c(0, 0, 0, 3),
             col.names = c("Categoría", "Valor 2011", "Valor 2025", "CAGR 2011-2025"))

Categoría	Valor 2011	Valor 2025	CAGR 2011-2025
Sun Care	1145	4326	0.100
Fragrances	18728	56742	0.082
Skin Care	25650	71830	0.076
Men’s Grooming	18672	49508	0.072
Bath and Shower	8410	21342	0.069
Deodorants	9151	21824	0.064
Depilatories	802	1873	0.062
Hair Care	24854	56290	0.060

cagr %>%
  mutate(Category = factor(Category, levels = Category)) %>%
  ggplot(aes(x = Category, y = cagr)) +
  geom_col(fill = COLOR_HIST, width = 0.65) +
  geom_text(aes(label = scales::percent(cagr, accuracy = 0.1)), hjust = -0.15, size = 3.3) +
  coord_flip(clip = "off") +
  scale_y_continuous(labels = scales::percent, expand = expansion(mult = c(0, 0.18))) +
  labs(title = "Crecimiento anual compuesto (CAGR) 2011-2025 por subcategoría",
       x = NULL, y = "CAGR") +
  theme_minimal(base_size = 11) +
  theme(panel.grid.minor = element_blank())

# 6. Corremos varios modelos para encontrar el que mejor se acomode a nuestros datos

comparar_modelos <- function(d) {
  m_lin <- lm(value ~ year, data = d)
  m_exp <- lm(log(value) ~ year, data = d)

  pred_lin <- predict(m_lin)
  pred_exp <- exp(predict(m_exp))

  rmse_lin <- sqrt(mean((d$value - pred_lin)^2))
  rmse_exp <- sqrt(mean((d$value - pred_exp)^2))

  data.frame(
    r2_lineal = summary(m_lin)$r.squared,
    r2_exponencial = summary(m_exp)$r.squared,
    rmse_lineal = rmse_lin,
    rmse_exponencial = rmse_exp,
    modelo_elegido = ifelse(rmse_exp < rmse_lin, "Exponencial", "Lineal"),
    crecimiento_anual_modelo = ifelse(rmse_exp < rmse_lin,
                                       exp(coef(m_exp)["year"]) - 1,
                                       coef(m_lin)["year"] / mean(d$value))
  )
}

seleccion <- datos %>%
  group_by(Category) %>%
  group_modify(~ comparar_modelos(.x)) %>%
  ungroup()

knitr::kable(seleccion, digits = 4, caption = "Comparación lineal vs. exponencial y modelo elegido por categoría")

Comparación lineal vs. exponencial y modelo elegido por categoría

Category	r2_lineal	r2_exponencial	rmse_lineal	rmse_exponencial	modelo_elegido	crecimiento_anual_modelo
Bath and Shower	0.9484	0.9866	904.9124	501.8314	Exponencial	0.0675
Deodorants	0.9141	0.9231	1098.0006	1065.2359	Exponencial	0.0583
Depilatories	0.9849	0.9689	42.2082	69.2755	Lineal	0.0576
Fragrances	0.8194	0.9033	4901.5038	3852.3353	Exponencial	0.0775
Hair Care	0.9279	0.9728	2447.2723	1773.3294	Exponencial	0.0557
Men’s Grooming	0.9240	0.9616	2531.3553	1948.1793	Exponencial	0.0672
Skin Care	0.9095	0.9703	4392.0379	2717.0337	Exponencial	0.0773
Sun Care	0.8627	0.9349	361.2479	277.1663	Exponencial	0.0924

7. Comparación de modelos

La tabla de la sección anterior ya es esa comparación: en 7 de las 8 categorías el modelo exponencial tiene menor RMSE (mejor ajuste) que el lineal; la excepción es Depilatories, donde el lineal se ajusta mejor. Por eso no se vuelve a correr el mismo cálculo aquí (la sección 6 original y esta tenían el mismo código pegado dos veces, así que se dejó solo una vez).

8. Pronósticos

futuro <- data.frame(year = 2026:2028)

pronosticar <- function(d) {
  m_lin <- lm(value ~ year, data = d)
  m_exp <- lm(log(value) ~ year, data = d)
  rmse_lin <- sqrt(mean((d$value - predict(m_lin))^2))
  rmse_exp <- sqrt(mean((d$value - exp(predict(m_exp)))^2))

  if (rmse_exp < rmse_lin) {
    valor <- exp(predict(m_exp, newdata = futuro))
    modelo <- "Exponencial"
  } else {
    valor <- predict(m_lin, newdata = futuro)
    modelo <- "Lineal"
  }
  data.frame(year = futuro$year, value = valor, modelo = modelo, tipo = "Pronóstico")
}

pron <- datos %>%
  group_by(Category) %>%
  group_modify(~ pronosticar(.x)) %>%
  ungroup()

datos_hist <- datos %>% mutate(tipo = "Histórico")

serie_completa <- bind_rows(
  datos_hist %>% select(Category, year, value, tipo),
  pron %>% select(Category, year, value, tipo)
)

knitr::kable(pron %>% select(Category, year, value, modelo) %>%
               pivot_wider(names_from = year, values_from = value),
             digits = 0, caption = "Pronóstico de valor de mercado 2026-2028 (MXN millones)")

Pronóstico de valor de mercado 2026-2028 (MXN millones)

Category	modelo	2026	2027	2028
Bath and Shower	Exponencial	22102	23595	25188
Deodorants	Exponencial	22763	24092	25497
Depilatories	Lineal	2001	2080	2159
Fragrances	Exponencial	52078	56113	60460
Hair Care	Exponencial	55318	58398	61650
Men’s Grooming	Exponencial	50702	54109	57746
Skin Care	Exponencial	72529	78134	84172
Sun Care	Exponencial	4413	4820	5265

ggplot(serie_completa, aes(x = year, y = value, color = tipo, linetype = tipo)) +
  geom_line(linewidth = 0.9) +
  geom_point(size = 1.4) +
  facet_wrap(~ Category, scales = "free_y", ncol = 4) +
  scale_color_manual(values = c("Histórico" = COLOR_HIST, "Pronóstico" = COLOR_FCST)) +
  scale_linetype_manual(values = c("Histórico" = "solid", "Pronóstico" = "22")) +
  labs(title = "Histórico (2011-2025) y pronóstico (2026-2028) por subcategoría",
       x = "Año", y = "MXN millones", color = NULL, linetype = NULL) +
  theme_minimal(base_size = 10) +
  theme(strip.text = element_text(face = "bold"), legend.position = "top")

9. ¿En qué subcategoría invertiría?

Skin Care. La justificación combina tres criterios, no solo el crecimiento histórico:

1. Tamaño y trayectoria: ya es la categoría más grande (~MXN 71,830 millones en 2025, proyectada a ~MXN 84,200 millones en 2028), así que invertir ahí no depende de que un nicho pequeño se vuelva grande — el mercado ya es grande y sigue creciendo a un ritmo alto (~7.6%/año compuesto).
2. Calidad y confiabilidad de la tendencia: el modelo exponencial ajusta con R² = 0.97 sobre 15 años de historia, uno de los mejores ajustes de toda la categoría. Eso da más confianza en el pronóstico que, por ejemplo, el de Fragrances (R² = 0.90) o Sun Care (R² = 0.93), cuyas series son más erráticas.
3. Resiliencia ante crisis: Skin Care fue una de las dos únicas categorías que no cayó en 2020 (+3.8%) y tuvo la mejor recuperación relativa hacia 2021 (+18.1% vs. 2019). Eso importa mucho para una decisión de inversión real: un mercado que crece rápido pero se desploma en la siguiente crisis (como Sun Care, -16.4% en 2020 y todavía -6.2% por debajo de 2019 un año después) es una apuesta más arriesgada.

Parte 3. Banco mundial

url_gdp   <- "https://raw.githubusercontent.com/datasets/gdp/master/data/gdp.csv"
url_infl  <- "https://raw.githubusercontent.com/datasets/inflation/master/data/inflation-consumer.csv"
url_pop   <- "https://raw.githubusercontent.com/datasets/population/master/data/population.csv"

gdp  <- read.csv(url_gdp)
infl <- read.csv(url_infl)
pop  <- read.csv(url_pop)

paises <- c("CHN", "JPN", "KOR")
nombres_paises <- c(CHN = "China", JPN = "Japón", KOR = "Corea del Sur")

gdp  <- gdp  %>% filter(Country.Code %in% paises) %>% select(Country.Code, Year, gdp_usd = Value)
pop  <- pop  %>% filter(Country.Code %in% paises) %>% select(Country.Code, Year, poblacion = Value)
infl <- infl %>% filter(Country.Code %in% paises) %>% select(Country.Code, Year, inflacion_pct = Inflation)

datos_bm <- gdp %>%
  inner_join(pop, by = c("Country.Code", "Year")) %>%
  inner_join(infl, by = c("Country.Code", "Year")) %>%
  filter(Year >= 1990, Year <= 2023) %>%
  mutate(
    pais = nombres_paises[Country.Code],
    gdp_percapita_usd = gdp_usd / poblacion
  ) %>%
  rename(anio = Year) %>%
  select(pais, Country.Code, anio, gdp_usd, poblacion, gdp_percapita_usd, inflacion_pct) %>%
  arrange(pais, anio)

cat("Filas:", nrow(datos_bm), " | Países:", n_distinct(datos_bm$pais), " | Años:", min(datos_bm$anio), "-", max(datos_bm$anio), "\n")

## Filas: 102  | Países: 3  | Años: 1990 - 2023

knitr::kable(datos_bm %>% filter(anio %in% c(1990, 2000, 2010, 2020, 2023)), digits = 1,
             caption = "Muestra de los datos importados")

Muestra de los datos importados

pais	Country.Code	anio	gdp_usd	poblacion	gdp_percapita_usd	inflacion_pct
China	CHN	1990	3.608579e+11	1135185000	317.9	5.7
China	CHN	2000	1.211332e+12	1262645000	959.4	2.1
China	CHN	2010	6.087192e+12	1337705000	4550.5	6.9
China	CHN	2020	1.468774e+13	1411100000	10408.7	0.5
China	CHN	2023	1.779478e+13	1410710000	12614.1	-0.6
Corea del Sur	KOR	1990	2.833658e+11	42869283	6610.0	10.1
Corea del Sur	KOR	2000	5.761794e+11	47008111	12257.0	1.0
Corea del Sur	KOR	2010	1.143672e+12	49554112	23079.3	2.7
Corea del Sur	KOR	2020	1.644313e+12	51836239	31721.3	1.6
Corea del Sur	KOR	2023	1.712793e+12	51712619	33121.4	2.1
Japón	JPN	1990	3.185905e+12	123478000	25801.4	2.6
Japón	JPN	2000	4.968359e+12	126843000	39169.4	-1.3
Japón	JPN	2010	5.759072e+12	128070000	44968.2	-1.9
Japón	JPN	2020	5.055587e+12	126261000	40040.8	0.9
Japón	JPN	2023	4.212945e+12	124516650	33834.4	3.8

Panel balanceado: 3 países × 34 años = 102 observaciones, sin datos faltantes.

Exploración

ggplot(datos_bm, aes(x = anio, y = gdp_percapita_usd, color = pais)) +
  geom_line(linewidth = 1) +
  scale_color_manual(values = c("China" = "#2a78d6", "Japón" = "#eb6834", "Corea del Sur" = "#1baf7a")) +
  labs(title = "PIB per cápita (US$ corrientes), 1990-2023",
       x = "Año", y = "PIB per cápita (US$)", color = NULL) +
  theme_minimal(base_size = 11)

Las diferencias de nivel son enormes: Japón y Corea del Sur ya eran economías de ingreso alto en 1990, mientras que China partió de un PIB per cápita muy bajo y ha ido cerrando la brecha con un crecimiento sostenido. Esta clase de diferencias estructurales entre países es exactamente lo que un modelo pooled (que ignora la identidad del país) no puede capturar bien.

ggplot(datos_bm, aes(x = anio, y = inflacion_pct, color = pais)) +
  geom_line(linewidth = 1) +
  geom_hline(yintercept = 0, color = "grey50", linetype = "dashed") +
  scale_color_manual(values = c("China" = "#2a78d6", "Japón" = "#eb6834", "Corea del Sur" = "#1baf7a")) +
  labs(title = "Inflación anual (%), 1990-2023",
       x = "Año", y = "Inflación (%)", color = NULL) +
  theme_minimal(base_size = 11)

Japón ha tenido inflación cercana a cero (o negativa, deflación) buena parte del periodo; China y Corea tuvieron picos de inflación más altos en los años 90 y se han moderado desde entonces.

Paso 1: estructura de panel y modelo de referencia

pdf_bm <- pdata.frame(datos_bm, index = c("pais", "anio"))

modelo_referencia_bm <- lm(gdp_percapita_usd ~ inflacion_pct, data = pdf_bm)
summary(modelo_referencia_bm)

## 
## Call:
## lm(formula = gdp_percapita_usd ~ inflacion_pct, data = pdf_bm)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -29089  -8392   2748   9522  23319 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    26913.7     1455.2  18.495  < 2e-16 ***
## inflacion_pct  -2413.5      322.9  -7.474 3.02e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12190 on 100 degrees of freedom
## Multiple R-squared:  0.3584, Adjusted R-squared:  0.352 
## F-statistic: 55.86 on 1 and 100 DF,  p-value: 3.02e-11

Paso 2: los tres modelos (pooling, within, random)

form_bm <- gdp_percapita_usd ~ inflacion_pct

pooling_bm <- plm(form_bm, data = pdf_bm, model = "pooling")
within_bm  <- plm(form_bm, data = pdf_bm, model = "within")
random_bm  <- plm(form_bm, data = pdf_bm, model = "random")

Pooling

summary(pooling_bm)

## Pooling Model
## 
## Call:
## plm(formula = form_bm, data = pdf_bm, model = "pooling")
## 
## Balanced Panel: n = 3, T = 34, N = 102
## 
## Residuals:
##     Min.  1st Qu.   Median  3rd Qu.     Max. 
## -29088.7  -8391.5   2748.3   9521.7  23318.5 
## 
## Coefficients:
##               Estimate Std. Error t-value  Pr(>|t|)    
## (Intercept)   26913.65    1455.17 18.4952 < 2.2e-16 ***
## inflacion_pct -2413.48     322.93 -7.4738  3.02e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    2.3176e+10
## Residual Sum of Squares: 1.487e+10
## R-Squared:      0.35839
## Adj. R-Squared: 0.35197
## F-statistic: 55.8571 on 1 and 100 DF, p-value: 3.0202e-11

Within / Efectos fijos

summary(within_bm)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = form_bm, data = pdf_bm, model = "within")
## 
## Balanced Panel: n = 3, T = 34, N = 102
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -13310.366  -4410.235    -26.126   3237.947  14535.611 
## 
## Coefficients:
##               Estimate Std. Error t-value  Pr(>|t|)    
## inflacion_pct  -849.76     179.68 -4.7294 7.559e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    4181900000
## Residual Sum of Squares: 3404800000
## R-Squared:      0.18583
## Adj. R-Squared: 0.1609
## F-statistic: 22.3675 on 1 and 98 DF, p-value: 7.559e-06

Random / Efectos aleatorios

summary(random_bm)

## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = form_bm, data = pdf_bm, model = "random")
## 
## Balanced Panel: n = 3, T = 34, N = 102
## 
## Effects:
##                    var  std.dev share
## idiosyncratic 34742880     5894 0.627
## individual    20630819     4542 0.373
## theta: 0.7828
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -13724.627  -5046.723    -19.408   4771.156  14505.859 
## 
## Coefficients:
##               Estimate Std. Error z-value  Pr(>|z|)    
## (Intercept)   23222.90    2958.95  7.8484 4.215e-15 ***
## inflacion_pct  -946.06     193.71 -4.8839 1.040e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    5078300000
## Residual Sum of Squares: 4100300000
## R-Squared:      0.19258
## Adj. R-Squared: 0.18451
## Chisq: 23.852 on 1 DF, p-value: 1.0403e-06

En los tres modelos el coeficiente de inflación es negativo y significativo: más inflación se asocia con menor PIB per cápita (o con países/años de menor PIB per cápita), consistente con la idea económica de que la inestabilidad de precios va de la mano de menor desarrollo relativo en esta muestra.

Paso 3: ¿cuál modelo es el correcto?

Pooling vs. efectos fijos

pFtest(within_bm, pooling_bm)

## 
##  F test for individual effects
## 
## data:  form_bm
## F = 165, df1 = 2, df2 = 98, p-value < 2.2e-16
## alternative hypothesis: significant effects

El p-value es prácticamente 0: se rechaza pooling. Hay efectos propios de cada país (su nivel de desarrollo, historia económica, etc.) que pooling está ignorando. Efectos fijos le gana a pooling.

Efectos fijos vs. efectos aleatorios (prueba de Hausman)

phtest(random_bm, within_bm)

## 
##  Hausman Test
## 
## data:  form_bm
## chisq = 1.7694, df = 1, p-value = 0.1835
## alternative hypothesis: one model is inconsistent

Aquí el resultado es distinto al del ejercicio de patentes: el p-value es 0.18, no significativo. Eso quiere decir que no se rechaza que random effects sea consistente, así que —a diferencia del ejercicio anterior— aquí sí conviene usar random effects, porque es más eficiente que efectos fijos (usa tanto la variación entre países como en el tiempo) y la prueba de Hausman no encontró evidencia de que esté sesgado.

Conclusión del Paso 3: el modelo elegido es efectos aleatorios (random).

Paso 4: predicciones 2024-2026

Random effects tiene un solo intercepto y una sola pendiente para los tres países (a diferencia de efectos fijos, que le da un intercepto propio a cada uno), así que pronosticar es tan simple como en tu ejemplo original: intercepto + pendiente * inflación. El único supuesto que necesitamos es un valor de inflación futura por país; usamos el último dato observado (2023) como escenario “si la inflación se mantiene igual”.

intercepto_bm <- coef(random_bm)["(Intercept)"]
pendiente_bm  <- coef(random_bm)["inflacion_pct"]

inflacion_2023 <- datos_bm %>% filter(anio == 2023) %>% select(pais, inflacion_pct)

pronostico_bm <- expand.grid(pais = unique(datos_bm$pais), anio = 2024:2026) %>%
  left_join(inflacion_2023, by = "pais") %>%
  mutate(gdp_percapita_pronostico = intercepto_bm + pendiente_bm * inflacion_pct)

knitr::kable(pronostico_bm %>% pivot_wider(names_from = anio, values_from = gdp_percapita_pronostico,
                                          id_cols = pais),
             digits = 0, caption = "PIB per cápita pronosticado 2024-2026 (US$, supone inflación = la de 2023)")

PIB per cápita pronosticado 2024-2026 (US$, supone inflación = la de 2023)

pais	2024	2025	2026
China	23774	23774	23774
Corea del Sur	21270	21270	21270
Japón	19634	19634	19634

serie_completa_bm <- bind_rows(
  datos_bm %>% select(pais, anio, valor = gdp_percapita_usd) %>% mutate(tipo = "Histórico"),
  pronostico_bm %>% select(pais, anio, valor = gdp_percapita_pronostico) %>% mutate(tipo = "Pronóstico")
)

ggplot(serie_completa_bm, aes(x = anio, y = valor, color = pais, linetype = tipo)) +
  geom_line(linewidth = 1) +
  scale_color_manual(values = c("China" = "#2a78d6", "Japón" = "#eb6834", "Corea del Sur" = "#1baf7a")) +
  scale_linetype_manual(values = c("Histórico" = "solid", "Pronóstico" = "22")) +
  labs(title = "PIB per cápita: histórico (1990-2023) y pronóstico (2024-2026)",
       x = "Año", y = "PIB per cápita (US$)", color = NULL, linetype = NULL) +
  theme_minimal(base_size = 11)

Por qué el pronóstico sale casi plano: como el modelo de random effects solo usa la inflación para explicar el PIB per cápita (y le pega el mismo intercepto/pendiente a los tres países), el pronóstico 2024-2026 no mueve mucho el nivel de cada país frente a su valor de 2023 — el modelo no tiene forma de “saber” que China ha venido creciendo de forma sostenida si esa tendencia no pasa por la inflación. Esto es una limitación real y vale la pena decirla en las conclusiones, no esconderla.

Conclusiones

• El PIB per cápita de China, Japón y Corea del Sur tiene diferencias de nivel tan grandes que un modelo pooled (que ignora la identidad del país) queda claramente rechazado por la prueba F (pFtest, p < 0.001).
• A diferencia del ejercicio de patentes, aquí la prueba de Hausman (phtest, p = 0.18) no rechaza efectos aleatorios, así que el modelo elegido es random effects — más eficiente que efectos fijos porque aprovecha también la variación entre países, no solo dentro de cada país en el tiempo.
• La relación estimada entre inflación y PIB per cápita es negativa y significativa en los tres modelos: cada punto porcentual adicional de inflación se asocia con, aproximadamente, US$946 menos de PIB per cápita (según el modelo de random effects).
• El ajuste del modelo es modesto (R² ≈ 0.19): la inflación por sí sola explica una parte relativamente pequeña de las diferencias de PIB per cápita entre estos tres países — tiene sentido, porque el desarrollo económico depende de muchos factores más (capital humano, instituciones, comercio, tecnología) que no están en este modelo de una sola variable. El pronóstico 2024-2026 hereda esa misma limitación: como el “shock” de inflación no cambia mucho de un año a otro, el pronóstico tampoco se mueve mucho — esto no es un error del código, es una consecuencia esperable de haber usado solo una variable explicativa.
• Nota de asociación vs. causalidad: esta correlación negativa entre inflación y PIB per cápita no debe leerse como que “bajar la inflación garantiza” mayor PIB per cápita — con 3 países y una sola variable, es evidencia sugestiva, no una prueba causal.