Este documento contiene el análisis de las tres partes del proyecto. Los datos provienen de dos bases entregadas en clase (convertidas a CSV) y de la API del Banco Mundial.

Archivos requeridos en la misma carpeta que este documento: patentes.csv y Market sizes.csv.

# Correr esta celda UNA sola vez, manualmente (no se ejecuta al compilar).
options(repos = c(CRAN = "https://cloud.r-project.org"))
install.packages(c("plm", "MASS", "car", "sandwich", "lmtest",
                   "ggplot2", "dplyr", "tidyr", "scales", "WDI"))
library(plm)        # modelos de datos de panel
library(MASS)       # binomial negativa
library(car)        # VIF
library(sandwich)   # errores estándar robustos
library(lmtest)     # pruebas de coeficientes
library(ggplot2)
library(dplyr)
library(tidyr)
library(scales)
library(WDI)        # datos del Banco Mundial
rmse <- function(obs, pred) sqrt(mean((obs - pred)^2, na.rm = TRUE))
mae  <- function(obs, pred) mean(abs(obs - pred), na.rm = TRUE)

# tasa de crecimiento anual compuesta
cagr <- function(valores, anios) {
  (valores[which.max(anios)] / valores[which.min(anios)])^(1 / (max(anios) - min(anios))) - 1
}

normalizar <- function(x) (x - min(x)) / (max(x) - min(x))

1 Parte 1: Patentes

pat <- read.csv("patentes.csv")
table(pat$year)
## 
## 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 
##  226  226  226  226  226  226  226  226  226  226

1.1 Estructura del panel

length(unique(pat$cusip))
## [1] 226
table(table(pat$cusip))
## 
##  10 
## 226

226 firmas con 10 años de datos cada una: el panel está balanceado.

1.2 Diagnóstico: la variable patents está truncada

aggregate(cbind(patents, patentsg) ~ year, data = pat, sum)
##    year patents patentsg
## 1  2012    5931     7244
## 2  2013    5884     7095
## 3  2014    6036     7000
## 4  2015    6198     6515
## 5  2016    5800     6464
## 6  2017    5811     5916
## 7  2018    5518     5932
## 8  2019    5313     4345
## 9  2020    4153     5230
## 10 2021    1117     5585

Las solicitudes de patente (patents) caen de 5,313 en 2019 a 1,117 en 2021, mientras las patentes otorgadas (patentsg) se mantienen estables. Si la caída fuera un fenómeno real, las otorgadas deberían bajar con uno o dos años de rezago, no sostenerse.

El patrón corresponde a un corte de la base hecho antes de que terminaran de registrarse las solicitudes de los últimos años. Decisión: usar patentsg como variable dependiente.

1.3 Deflactar las variables monetarias

rnd y sales están en dólares corrientes. rndeflt es la misma serie de I+D expresada en dólares del año base, de modo que su razón entrega el deflactor implícito de la propia base.

deflactor <- tapply(pat$rnd / pat$rndeflt, pat$year, median, na.rm = TRUE)
deflactor
##  2012  2013  2014  2015  2016  2017  2018  2019  2020  2021 
## 1.000 1.064 1.170 1.285 1.361 1.459 1.573 1.718 1.870 2.010

El deflactor pasa de 1.00 a 2.01 en diez años, equivalente a 7% anual. Ese ritmo corresponde a Estados Unidos en los años setenta, no al periodo 2012-2021, y el diccionario de la base indica que los años van de “72 through 81”. Los años de la columna year parecen estar reetiquetados. [VERIFICAR con el profesor]; no afecta la estimación.

pat$sales_real <- pat$sales / deflactor[as.character(pat$year)]

# log1p = log(1 + x), necesario porque hay observaciones en cero
pat$lpatentsg <- log1p(pat$patentsg)
pat$lrndstck  <- log1p(pat$rndstck)
pat$lsales    <- log1p(pat$sales_real)
pat$lemploy   <- log1p(pat$employ)

1.4 Diagnóstico: sobredispersión

var(pat$patentsg) / mean(pat$patentsg)
## [1] 231.1284
mean(pat$patentsg == 0)
## [1] 0.1889381

La razón varianza/media supera 200, cuando el modelo Poisson supone que vale 1. Se requiere binomial negativa o errores estándar robustos.

1.5 Diagnóstico: multicolinealidad

m_vif <- lm(lpatentsg ~ lrndstck + lsales + lemploy + merger, data = pat)
vif(m_vif)
##  lrndstck    lsales   lemploy    merger 
##  4.804379 13.472763 15.292813  1.035471

sales y employ tienen VIF por encima de 10: ambas miden tamaño de la firma. Decisión: conservar sales como control de tamaño y excluir employ.

d <- pat %>% filter(!is.na(patentsg), !is.na(lrndstck), !is.na(lsales), !is.na(merger))
d$t <- d$year - min(d$year)
c(observaciones = nrow(d), firmas = length(unique(d$cusip)))
## observaciones        firmas 
##          2100           216

1.6 Modelo 1: MCO agrupado

m1 <- lm(lpatentsg ~ lrndstck + lsales + merger + factor(year), data = d)
summary(m1)
## 
## Call:
## lm(formula = lpatentsg ~ lrndstck + lsales + merger + factor(year), 
##     data = d)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.88635 -0.55369  0.06853  0.61265  2.30205 
## 
## Coefficients:
##                  Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)      -0.20602    0.08339  -2.471             0.013566 *  
## lrndstck          0.66501    0.02133  31.180 < 0.0000000000000002 ***
## lsales            0.10141    0.02089   4.854      0.0000012993072 ***
## merger            0.44545    0.14589   3.053             0.002292 ** 
## factor(year)2013 -0.02930    0.08322  -0.352             0.724810    
## factor(year)2014 -0.09871    0.08312  -1.187             0.235177    
## factor(year)2015 -0.19377    0.08325  -2.328             0.020026 *  
## factor(year)2016 -0.30453    0.08339  -3.652             0.000267 ***
## factor(year)2017 -0.48318    0.08381  -5.765      0.0000000093839 ***
## factor(year)2018 -0.54722    0.08398  -6.516      0.0000000000901 ***
## factor(year)2019 -0.85633    0.08407 -10.185 < 0.0000000000000002 ***
## factor(year)2020 -0.78802    0.08432  -9.346 < 0.0000000000000002 ***
## factor(year)2021 -0.88911    0.08458 -10.512 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8403 on 2087 degrees of freedom
## Multiple R-squared:  0.7202, Adjusted R-squared:  0.7186 
## F-statistic: 447.6 on 12 and 2087 DF,  p-value: < 0.00000000000000022

La elasticidad del acervo de I+D es 0.66. Sin embargo, este modelo compara firmas distintas entre sí, por lo que confunde el nivel característico de cada firma con el efecto de aumentar el gasto.

1.7 Modelo 2: efectos fijos de firma y año

pdatos <- pdata.frame(d, index = c("cusip", "year"))
m2_fe <- plm(lpatentsg ~ lrndstck + lsales + merger, data = pdatos,
             model = "within", effect = "twoways")
summary(m2_fe)
## Twoways effects Within Model
## 
## Call:
## plm(formula = lpatentsg ~ lrndstck + lsales + merger, data = pdatos, 
##     effect = "twoways", model = "within")
## 
## Unbalanced Panel: n = 216, T = 2-10, N = 2100
## 
## Residuals:
##    Min. 1st Qu.  Median 3rd Qu.    Max. 
## -2.0611 -0.2911  0.0039  0.2899  1.6132 
## 
## Coefficients:
##          Estimate Std. Error t-value  Pr(>|t|)    
## lrndstck 0.274817   0.067781  4.0545 0.0000523 ***
## lsales   0.141325   0.050346  2.8071  0.005051 ** 
## merger   0.096788   0.108298  0.8937  0.371585    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    455.35
## Residual Sum of Squares: 446.89
## R-Squared:      0.018572
## Adj. R-Squared: -0.10044
## F-statistic: 11.8079 on 3 and 1872 DF, p-value: 0.00000011592

La elasticidad cae de 0.66 a 0.27 al controlar por firma. La mayor parte de la asociación original provenía de comparar firmas entre sí, no de la variación dentro de cada una.

m2_re <- plm(lpatentsg ~ lrndstck + lsales + merger, data = pdatos,
             model = "random", effect = "twoways")
phtest(m2_fe, m2_re)
## 
##  Hausman Test
## 
## data:  lpatentsg ~ lrndstck + lsales + merger
## chisq = 5.2329, df = 3, p-value = 0.1555
## alternative hypothesis: one model is inconsistent

El test de Hausman no rechaza la hipótesis nula (p = 0.16). Aun así se opta por efectos fijos: las características no observadas de cada firma (estrategia de propiedad intelectual, sector, cultura de innovación) plausiblemente correlacionan con su gasto en I+D, que es precisamente lo que los efectos fijos controlan. Es la especificación conservadora.

1.8 Modelo 3: Poisson con efectos fijos

patentsg es una variable de conteo, por lo que un modelo Poisson resulta más apropiado que MCO sobre el logaritmo. Los efectos fijos se incorporan como variables dummy.

m3_poisson <- glm(patentsg ~ lrndstck + lsales + merger + factor(cusip) + factor(year),
                  family = poisson, data = d)
coef(summary(m3_poisson))[c("lrndstck", "lsales", "merger"), ]
##            Estimate Std. Error   z value               Pr(>|z|)
## lrndstck 0.36029900 0.04741299 7.5991629 0.00000000000002980523
## lsales   0.22809657 0.03124494 7.3002721 0.00000000000028718631
## merger   0.03593022 0.04130652 0.8698438 0.38438578383738941646

Dada la sobredispersión detectada, los errores estándar de Poisson subestiman la varianza. Se corrigen agrupando por firma:

coeftest(m3_poisson, vcov. = vcovCL(m3_poisson, cluster = d$cusip))[c("lrndstck", "lsales", "merger"), ]
##            Estimate Std. Error   z value  Pr(>|z|)
## lrndstck 0.36029900  0.2300194 1.5663853 0.1172584
## lsales   0.22809657  0.1620094 1.4079219 0.1591542
## merger   0.03593022  0.1232704 0.2914747 0.7706883

Con errores robustos, la elasticidad del gasto en I+D deja de ser estadísticamente significativa (p = 0.12). Controlando por firma y por la correlación de los errores dentro de cada una, los datos no sostienen que un mayor gasto en I+D genere más patentes en el corto plazo.

1.9 Modelo 4: binomial negativa

m4_nb <- glm.nb(patentsg ~ lrndstck + lsales + merger + factor(year), data = d)
coeftest(m4_nb, vcov. = vcovCL(m4_nb, cluster = d$cusip))
## 
## z test of coefficients:
## 
##                   Estimate Std. Error  z value              Pr(>|z|)    
## (Intercept)      -0.691348   0.187686  -3.6835             0.0002300 ***
## lrndstck          0.707050   0.066459  10.6390 < 0.00000000000000022 ***
## lsales            0.208021   0.065711   3.1657             0.0015472 ** 
## merger            0.377248   0.284567   1.3257             0.1849427    
## factor(year)2013 -0.096053   0.047602  -2.0178             0.0436104 *  
## factor(year)2014 -0.174142   0.051908  -3.3548             0.0007941 ***
## factor(year)2015 -0.238204   0.067600  -3.5237             0.0004255 ***
## factor(year)2016 -0.428506   0.058971  -7.2664    0.0000000000003693 ***
## factor(year)2017 -0.637120   0.067076  -9.4985 < 0.00000000000000022 ***
## factor(year)2018 -0.698260   0.074468  -9.3766 < 0.00000000000000022 ***
## factor(year)2019 -1.089774   0.076267 -14.2889 < 0.00000000000000022 ***
## factor(year)2020 -0.930017   0.083712 -11.1097 < 0.00000000000000022 ***
## factor(year)2021 -1.047488   0.092200 -11.3611 < 0.00000000000000022 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
c(AIC_poisson = AIC(glm(patentsg ~ lrndstck + lsales + merger + factor(year),
                        family = poisson, data = d)),
  AIC_binomial_negativa = AIC(m4_nb))
##           AIC_poisson AIC_binomial_negativa 
##              44290.29              13282.32

El AIC baja de aproximadamente 44,000 a 13,000 frente a Poisson. Es el modelo mejor especificado de los cuatro.

1.10 Validación fuera de muestra

La partición es temporal, no aleatoria. Con datos de panel, mezclar años permitiría al modelo aprovechar información futura a través de los efectos fijos, inflando artificialmente el desempeño.

train <- d %>% filter(year <= 2018)
test  <- d %>% filter(year >= 2019)

m_val <- glm.nb(patentsg ~ lrndstck + lsales + merger + t, data = train)
pred_modelo <- predict(m_val, newdata = test, type = "response")

pred_media   <- rep(mean(train$patentsg), nrow(test))
ultimo_valor <- train$patentsg[match(test$cusip, train$cusip[train$year == 2018])]

data.frame(
  modelo = c("media histórica", "valor del año anterior", "binomial negativa"),
  RMSE   = c(rmse(test$patentsg, pred_media),
             rmse(test$patentsg, ultimo_valor),
             rmse(test$patentsg, pred_modelo)),
  MAE    = c(mae(test$patentsg, pred_media),
             mae(test$patentsg, ultimo_valor),
             mae(test$patentsg, pred_modelo))
)
##                   modelo     RMSE      MAE
## 1        media histórica 69.56094 35.71745
## 2 valor del año anterior 88.40002 34.57257
## 3      binomial negativa 73.36828 19.58772
aggregate(cbind(patents, patentsg) ~ year, data = pat, sum) %>%
  pivot_longer(-year, names_to = "tipo", values_to = "total") %>%
  ggplot(aes(year, total, color = tipo)) +
  geom_line(linewidth = 1) +
  geom_point() +
  labs(title = "Solicitudes de patente contra patentes otorgadas",
       subtitle = "Las solicitudes caen en 2020-2021; las otorgadas se mantienen",
       x = NULL, y = "Total (226 firmas)", color = NULL) +
  theme_minimal()

2 Parte 2: Cuidado de la piel

crudo <- read.csv("Market sizes.csv", check.names = FALSE)

mercado <- crudo %>%
  select(-Geography, -`Data Type`, -Unit, -`Current Constant`) %>%
  pivot_longer(-Category, names_to = "anio", values_to = "valor") %>%
  mutate(anio = as.integer(anio)) %>%
  rename(categoria = Category) %>%
  filter(!is.na(valor))

unique(mercado$categoria)
##  [1] "Beauty and Personal Care"               
##  [2] "Bath and Shower"                        
##  [3] "Deodorants"                             
##  [4] "Depilatories"                           
##  [5] "Fragrances"                             
##  [6] "Hair Care"                              
##  [7] "Men's Grooming"                         
##  [8] "Skin Care"                              
##  [9] "Sun Care"                               
## [10] "Premium Beauty and Personal Care"       
## [11] "Prestige Beauty and Personal Care"      
## [12] "Mass Beauty and Personal Care"          
## [13] "Dermocosmetics Beauty and Personal Care"

2.1 Validación de la taxonomía

La base contiene 13 filas: el total del mercado, 8 categorías de producto y 4 segmentos de posicionamiento por precio. Antes de sumar cualquier cosa hay que verificar que las filas no se traslapen entre sí.

total <- mercado %>%
  filter(categoria == "Beauty and Personal Care") %>%
  select(anio, total = valor)

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

suma_8 <- mercado %>%
  filter(categoria %in% categorias_producto) %>%
  group_by(anio) %>%
  summarise(suma = sum(valor))

left_join(total, suma_8, by = "anio") %>% mutate(razon = suma / total)
## # A tibble: 15 × 4
##     anio   total    suma razon
##    <int>   <dbl>   <dbl> <dbl>
##  1  2011 126329. 107413. 0.850
##  2  2012 136020. 115957. 0.853
##  3  2013 142539. 121621. 0.853
##  4  2014 149165. 127034. 0.852
##  5  2015 157656. 134716. 0.854
##  6  2016 170488. 145430. 0.853
##  7  2017 183580. 156175. 0.851
##  8  2018 192727. 164808. 0.855
##  9  2019 200916. 172726. 0.860
## 10  2020 191490. 166336. 0.869
## 11  2021 211921. 186692. 0.881
## 12  2022 236406. 209792. 0.887
## 13  2023 268852. 240652. 0.895
## 14  2024 296120. 267536. 0.903
## 15  2025 313593. 283735. 0.905
suma_7 <- mercado %>%
  filter(categoria %in% setdiff(categorias_producto, "Men's Grooming")) %>%
  group_by(anio) %>%
  summarise(suma = sum(valor))

left_join(total, suma_7, by = "anio") %>% mutate(residual = total - suma)
## # A tibble: 15 × 4
##     anio   total    suma residual
##    <int>   <dbl>   <dbl>    <dbl>
##  1  2011 126329.  88741.   37588 
##  2  2012 136020.  95263.   40757.
##  3  2013 142539.  99487.   43052.
##  4  2014 149165. 103858.   45307.
##  5  2015 157656. 109942.   47715.
##  6  2016 170488. 117984.   52504.
##  7  2017 183580. 126678.   56902 
##  8  2018 192727. 133310    59417.
##  9  2019 200916. 140032.   60883.
## 10  2020 191490. 136666.   54824.
## 11  2021 211921. 153754.   58168.
## 12  2022 236406. 172090.   64316.
## 13  2023 268852. 197285.   71567.
## 14  2024 296120. 220202.   75918.
## 15  2025 313593. 234227.   79366

El export de Euromonitor no incluye Oral Care, Color Cosmetics ni Baby & Child-specific, que en México suman alrededor de 55 a 60 mil millones de pesos en 2021.

  • Excluyendo Men's Grooming, el residual sin explicar es de 58 mil millones, consistente con esas tres categorías ausentes.
  • Incluyéndola, el residual baja a 25 mil millones, insuficiente para contenerlas.

Conclusión: Men's Grooming es una categoría transversal que ya incluye los productos masculinos contados en Bath and Shower, Deodorants, Depilatories, Fragrances y Skin Care. Sumarla junto con las demás duplica parte del mercado.

mercado %>%
  filter(categoria %in% c("Premium Beauty and Personal Care",
                          "Prestige Beauty and Personal Care")) %>%
  pivot_wider(names_from = categoria, values_from = valor) %>%
  mutate(razon = `Prestige Beauty and Personal Care` / `Premium Beauty and Personal Care`)
## # A tibble: 15 × 4
##     anio `Premium Beauty and Personal Care` Prestige Beauty and Persona…¹  razon
##    <int>                              <dbl>                         <dbl>  <dbl>
##  1  2011                             13022.                           NA  NA    
##  2  2012                             14580                            NA  NA    
##  3  2013                             15698.                           NA  NA    
##  4  2014                             17023.                           NA  NA    
##  5  2015                             18878.                           NA  NA    
##  6  2016                             21783.                        19715   0.905
##  7  2017                             24478.                        21722.  0.887
##  8  2018                             26344                         23190.  0.880
##  9  2019                             27806.                        24245.  0.872
## 10  2020                             23213.                        19409.  0.836
## 11  2021                             30225.                        25615.  0.847
## 12  2022                             37375.                        31798.  0.851
## 13  2023                             46018.                        38291.  0.832
## 14  2024                             53798.                        44287.  0.823
## 15  2025                             58846.                        48538.  0.825
## # ℹ abbreviated name: ¹​`Prestige Beauty and Personal Care`

La razón Prestige/Premium se mantiene entre 0.82 y 0.91 en todos los años: Prestige está contenido dentro de Premium. Premium más Mass cubren 88-90% del total, lo que confirma que estos cuatro segmentos son un corte del mismo mercado por nivel de precio, no componentes adicionales.

Para calcular el tamaño total del mercado se usan las 7 categorías sumables. Para comparar tendencias de crecimiento se conservan las 8, porque Men's Grooming es una serie válida en sí misma; lo que no puede hacerse es sumarla con las demás.

categorias_comparar <- categorias_producto

2.2 Deflactar a pesos constantes

La base está en precios corrientes, como indica su columna Current Constant. Se descargan tres indicadores de México del Banco Mundial: índice de precios, consumo privado per cápita y población.

mexico_wb <- WDI(
  indicator = c(ipc = "FP.CPI.TOTL", consumo_pc = "NE.CON.PRVT.PC.KD",
                poblacion = "SP.POP.TOTL"),
  country = "MX", start = 2011, end = 2025
) %>%
  select(anio = year, ipc, consumo_pc, poblacion) %>%
  arrange(anio)

inpc <- mexico_wb %>% select(anio, ipc)
inpc$factor <- inpc$ipc[inpc$anio == 2025] / inpc$ipc

# inflación anual promedio 2011-2025
((inpc$ipc[inpc$anio == 2025] / inpc$ipc[inpc$anio == 2011])^(1/14)) - 1
## [1] 0.04497999

La inflación corrió a 4.5% anual entre 2011 y 2025. Comparar categorías en pesos corrientes confundiría crecimiento de mercado con aumento de precios.

datos_piel <- mercado %>%
  filter(categoria %in% categorias_comparar) %>%
  left_join(inpc, by = "anio") %>%
  mutate(valor_real = valor * factor)

2.3 Crecimiento nominal contra crecimiento real

datos_piel %>%
  filter(anio >= 2021) %>%
  group_by(categoria) %>%
  summarise(crecimiento_nominal = cagr(valor, anio),
            crecimiento_real    = cagr(valor_real, anio)) %>%
  arrange(desc(crecimiento_real))
## # A tibble: 8 × 3
##   categoria       crecimiento_nominal crecimiento_real
##   <chr>                         <dbl>            <dbl>
## 1 Sun Care                     0.172           0.111  
## 2 Fragrances                   0.156           0.0957 
## 3 Men's Grooming               0.107           0.0497 
## 4 Skin Care                    0.106           0.0487 
## 5 Deodorants                   0.0966          0.0397 
## 6 Hair Care                    0.0922          0.0355 
## 7 Bath and Shower              0.0838          0.0275 
## 8 Depilatories                 0.0462         -0.00817

Depilatories pasa de +4.6% nominal a -0.8% real: en términos reales el mercado se está contrayendo. Deflactar reordena el ranking de categorías, y es la diferencia entre identificar una categoría sana y una en declive.

2.4 Modelo de panel: tendencia y choque de 2020

datos_piel <- datos_piel %>%
  mutate(log_valor = log(valor_real),
         t = anio - min(anio),
         covid = as.integer(anio == 2020))

pdatos_piel <- pdata.frame(datos_piel, index = c("categoria", "anio"))
m_piel <- plm(log_valor ~ t + covid, data = pdatos_piel, model = "within")
summary(m_piel)
## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = log_valor ~ t + covid, data = pdatos_piel, model = "within")
## 
## Balanced Panel: n = 8, T = 15, N = 120
## 
## Residuals:
##     Min.  1st Qu.   Median  3rd Qu.     Max. 
## -0.19570 -0.04606 -0.00327  0.04541  0.21506 
## 
## Coefficients:
##         Estimate Std. Error t-value              Pr(>|t|)    
## t      0.0235850  0.0016272 14.4941 < 0.00000000000000022 ***
## covid -0.1036209  0.0281841 -3.6766             0.0003673 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    1.8904
## Residual Sum of Squares: 0.64244
## R-Squared:      0.66016
## Adj. R-Squared: 0.63236
## F-statistic: 106.841 on 2 and 110 DF, p-value: < 0.000000000000000222
c(tendencia_real_anual = exp(coef(m_piel)["t"]) - 1,
  efecto_2020          = exp(coef(m_piel)["covid"]) - 1)
## tendencia_real_anual.t      efecto_2020.covid 
##             0.02386537            -0.09843296

Un coeficiente común para todas las categorías oculta la heterogeneidad. Se estima la misma especificación categoría por categoría:

tendencias <- data.frame()
for (cat in categorias_comparar) {
  sub <- datos_piel %>% filter(categoria == cat)
  m <- lm(log_valor ~ t + covid, data = sub)
  tendencias <- rbind(tendencias, data.frame(
    categoria       = cat,
    tendencia_anual = exp(coef(m)["t"]) - 1,
    choque_2020     = exp(coef(m)["covid"]) - 1,
    r2              = summary(m)$r.squared
  ))
}
tendencias %>% arrange(desc(tendencia_anual))
##          categoria tendencia_anual choque_2020        r2
## t7        Sun Care      0.04648242 -0.20392729 0.8739771
## t3      Fragrances      0.03241716 -0.22592584 0.8334621
## t6       Skin Care      0.03085562 -0.06797606 0.9308877
## t5  Men's Grooming      0.02162387 -0.11782962 0.8398062
## t  Bath and Shower      0.02092592  0.01096959 0.9390379
## t2    Depilatories      0.01580745  0.03135244 0.5808480
## t1      Deodorants      0.01323894 -0.12939910 0.4886910
## t4       Hair Care      0.01006495 -0.05091508 0.7016514

Sun Care presenta la tendencia real más alta (4.6% anual) pero también el mayor choque en 2020 (-20%). Bath and Shower y Skin Care son las más estables: crecen menos, pero apenas resintieron la pandemia.

2.5 Elasticidad ingreso

datos_piel2 <- datos_piel %>%
  left_join(mexico_wb %>% select(anio, consumo_pc), by = "anio") %>%
  mutate(log_consumo = log(consumo_pc))

pdatos_piel2 <- pdata.frame(datos_piel2, index = c("categoria", "anio"))
m_ingreso_niveles <- plm(log_valor ~ log_consumo + covid, data = pdatos_piel2,
                         model = "within")
coef(m_ingreso_niveles)["log_consumo"]
## log_consumo 
##    2.952284

La elasticidad en niveles es cercana a 3, una magnitud implausible: implicaría que el gasto en belleza crece tres veces más rápido que el ingreso. Ambas series tienen tendencia creciente y el modelo atribuye a la relación lo que en realidad es tendencia compartida. Es una regresión espuria.

En primeras diferencias la tendencia común desaparece:

datos_piel2 <- datos_piel2 %>%
  arrange(categoria, anio) %>%
  group_by(categoria) %>%
  mutate(d_log_valor   = log_valor - lag(log_valor),
         d_log_consumo = log_consumo - lag(log_consumo)) %>%
  ungroup() %>%
  filter(!is.na(d_log_valor))

pdatos_dif <- pdata.frame(datos_piel2, index = c("categoria", "anio"))
m_ingreso_dif <- plm(d_log_valor ~ d_log_consumo + covid, data = pdatos_dif,
                     model = "within")
coef(m_ingreso_dif)["d_log_consumo"]
## d_log_consumo 
##     0.5976338

Esta es la estimación que debe reportarse.

2.6 Validación fuera de muestra

Se entrena con 2011-2021 y se prueba contra 2022-2025, años que ya conocemos.

entrena <- datos_piel %>% filter(anio <= 2021)
prueba  <- datos_piel %>% filter(anio >= 2022)

m_pred <- lm(log_valor ~ categoria * t + covid, data = entrena)
prueba$prediccion <- exp(predict(m_pred, newdata = prueba))

ultimo_2021 <- entrena %>% filter(anio == 2021) %>% select(categoria, base = valor_real)
prueba <- prueba %>% left_join(ultimo_2021, by = "categoria")

data.frame(
  enfoque = c("sin crecimiento", "modelo de panel"),
  RMSE    = c(rmse(prueba$valor_real, prueba$base),
              rmse(prueba$valor_real, prueba$prediccion))
)
##           enfoque     RMSE
## 1 sin crecimiento 6346.408
## 2 modelo de panel 4661.179

El modelo de panel mejora sustancialmente frente a la referencia ingenua.

2.7 Decisión de inversión

El criterio discutido en clase es que la categoría de mayor crecimiento no es necesariamente la mejor inversión, porque el crecimiento atrae competencia. El puntaje combina cuatro criterios: crecimiento, resiliencia, escala y competencia.

tabla <- tendencias %>%
  left_join(datos_piel %>% filter(anio >= 2021) %>% group_by(categoria) %>%
              summarise(crecimiento_reciente = cagr(valor_real, anio)),
            by = "categoria") %>%
  left_join(datos_piel %>% filter(anio == 2025) %>%
              select(categoria, mercado_2025 = valor_real),
            by = "categoria")

2.7.1 Especificación inicial y por qué se descarta

Un primer enfoque penaliza el nivel de crecimiento reciente como aproximación a la competencia esperada:

tabla$puntaje_inicial <- 0.35 * normalizar(tabla$tendencia_anual) +
                         0.25 * normalizar(tabla$choque_2020) +
                         0.20 * normalizar(log(tabla$mercado_2025)) +
                         0.20 * (1 - normalizar(tabla$crecimiento_reciente))

tabla %>% arrange(desc(puntaje_inicial)) %>%
  select(categoria, crecimiento_reciente, puntaje_inicial)
##         categoria crecimiento_reciente puntaje_inicial
## 1       Skin Care          0.048698068       0.6579173
## 2 Bath and Shower          0.027495686       0.6082023
## 3    Depilatories         -0.008166441       0.5051898
## 4  Men's Grooming          0.049745106       0.4985838
## 5       Hair Care          0.035495754       0.4834556
## 6      Fragrances          0.095739354       0.4276105
## 7        Sun Care          0.111073188       0.4172951
## 8      Deodorants          0.039667053       0.3787376

Esta especificación coloca a Depilatories en primer lugar, que es justamente la única categoría con crecimiento real negativo. El defecto es que el término 1 - crecimiento normalizado premia el estancamiento en lugar de penalizar la competencia: cuanto menos crece una categoría, más puntos recibe.

Se corrige de dos formas:

  1. Un filtro previo: el mercado debe crecer en términos reales per cápita. Si crece menos que la población, el consumo por persona está cayendo y la categoría queda fuera antes de puntuar.
  2. Penalizar el exceso de crecimiento sobre la mediana del mercado, no el nivel. Crecer al ritmo del mercado no debería castigarse.

2.7.2 Especificación corregida

crecimiento_pc <- datos_piel %>%
  left_join(mexico_wb %>% select(anio, poblacion), by = "anio") %>%
  mutate(valor_pc = valor_real / poblacion) %>%
  filter(anio >= 2021) %>%
  group_by(categoria) %>%
  summarise(crecimiento_pc = cagr(valor_pc, anio))

tabla <- tabla %>% left_join(crecimiento_pc, by = "categoria")
tabla$pasa_filtro <- tabla$crecimiento_pc > 0

mediana <- median(tabla$crecimiento_reciente)
exceso  <- pmax(0, tabla$crecimiento_reciente - mediana)

tabla$puntaje <- 0.35 * normalizar(tabla$tendencia_anual) +
                 0.25 * normalizar(tabla$choque_2020) +
                 0.20 * normalizar(log(tabla$mercado_2025)) +
                 0.20 * (1 - exceso / max(exceso))

tabla %>% arrange(desc(puntaje)) %>%
  select(categoria, tendencia_anual, crecimiento_pc, pasa_filtro, puntaje)
##         categoria tendencia_anual crecimiento_pc pasa_filtro   puntaje
## 1       Skin Care      0.03085562     0.04005019        TRUE 0.7397946
## 2 Bath and Shower      0.02092592     0.01902265        TRUE 0.6680182
## 3  Men's Grooming      0.02162387     0.04108860        TRUE 0.5790867
## 4       Hair Care      0.01006495     0.02695675        TRUE 0.5566900
## 5    Depilatories      0.01580745    -0.01634539       FALSE 0.5051898
## 6      Deodorants      0.01323894     0.03109365        TRUE 0.4589685
## 7      Fragrances      0.03241716     0.08670356        TRUE 0.4477386
## 8        Sun Care      0.04648242     0.10191095        TRUE 0.4172951

Depilatories queda descartada por el filtro: su gasto real per cápita se contrae 1.6% anual. Skin Care encabeza el tablero, por combinar una tendencia sólida (3.1% real), la mayor escala del mercado y el menor choque en 2020 entre las categorías grandes.

2.7.3 Sensibilidad a los pesos

Los pesos son un supuesto del analista, no un resultado del modelo. Conviene verificar cuánto depende la conclusión de ellos:

esquemas <- list(
  solo_crecimiento = c(1.00, 0.00, 0.00, 0.00),
  balanceado       = c(0.35, 0.25, 0.20, 0.20),
  averso_al_riesgo = c(0.20, 0.45, 0.20, 0.15)
)

for (nombre in names(esquemas)) {
  w <- esquemas[[nombre]]
  p <- w[1] * normalizar(tabla$tendencia_anual) +
       w[2] * normalizar(tabla$choque_2020) +
       w[3] * normalizar(log(tabla$mercado_2025)) +
       w[4] * (1 - exceso / max(exceso))
  elegibles <- tabla$pasa_filtro
  cat(nombre, "->", tabla$categoria[elegibles][which.max(p[elegibles])], "\n")
}
## solo_crecimiento -> Sun Care 
## balanceado -> Skin Care 
## averso_al_riesgo -> Bath and Shower

El ganador cambia según los pesos: Skin Care con pesos balanceados, Bath and Shower si se prioriza la resiliencia, y Sun Care únicamente si el crecimiento es el único criterio.

Esto no es una debilidad del análisis, es el resultado: la decisión depende de la tolerancia al riesgo del inversionista, no solo de los datos. La recomendación se sostiene en el esquema balanceado, y la sensibilidad debe reportarse de forma explícita.

ggplot(tabla, aes(tendencia_anual, choque_2020, size = mercado_2025)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_point(color = "steelblue", alpha = 0.7) +
  geom_text(aes(label = categoria), vjust = -1.3, size = 3) +
  scale_x_continuous(labels = percent) +
  labs(title = "Crecimiento real contra sensibilidad al ciclo",
       subtitle = "Tendencia real anual y efecto estimado de 2020",
       x = "Tendencia real anual", y = "Choque de 2020", size = "Mercado 2025") +
  theme_minimal()

3 Parte 3: Banco Mundial

Se construye un panel de países con patentes, gasto en I+D, PIB per cápita y población. El análisis se agrupa por décadas, según la consigna.

wb <- WDI(
  indicator = c(patentes = "IP.PAT.RESD", id_pib = "GB.XPD.RSDV.GD.ZS",
                pib_pc = "NY.GDP.PCAP.KD", poblacion = "SP.POP.TOTL"),
  country = "all", start = 1996, end = 2021, extra = TRUE
)

El panel arranca en 1996 porque el indicador de gasto en I+D no existe antes de ese año. La columna region permite excluir los agregados regionales (World, OECD members, etc.), que la API devuelve mezclados con los países.

wb <- wb %>%
  filter(region != "Aggregates",
         !is.na(patentes), !is.na(id_pib), !is.na(pib_pc), !is.na(poblacion))

# al menos 15 años por país, para no estimar efectos fijos con muy pocas
# observaciones
paises_ok <- names(which(table(wb$iso3c) >= 15))
wb <- wb %>% filter(iso3c %in% paises_ok)

c(paises = length(unique(wb$iso3c)), observaciones = nrow(wb))
##        paises observaciones 
##            69          1566
range(wb$year)
## [1] 1996 2021
wb$decada     <- paste0(floor(wb$year / 10) * 10, "s")
wb$lid        <- log(wb$id_pib)
wb$lpib       <- log(wb$pib_pc)
wb$lpob       <- log(wb$poblacion)
wb$pat_millon <- wb$patentes / (wb$poblacion / 1e6)

3.1 Agrupación por décadas

wb %>%
  group_by(decada) %>%
  summarise(paises              = n_distinct(iso3c),
            id_pib_promedio     = mean(id_pib),
            patentes_por_millon = median(pat_millon))
## # A tibble: 4 × 4
##   decada paises id_pib_promedio patentes_por_millon
##   <chr>   <int>           <dbl>               <dbl>
## 1 1990s      59            1.04                62.4
## 2 2000s      69            1.12                52.5
## 3 2010s      69            1.24                46.6
## 4 2020s      66            1.42                36.9

El gasto en I+D como porcentaje del PIB sube en cada década, mientras las patentes por millón de habitantes bajan. El mundo invierte más en investigación y obtiene menos patentes por habitante.

3.2 Escalera de modelos

var(wb$patentes) / mean(wb$patentes)
## [1] 496690.5

La sobredispersión es aún mayor que en la Parte 1: los conteos van de cero a más de 180 mil patentes anuales.

m1_wb <- lm(log(patentes + 1) ~ lid + lpib + lpob + factor(year), data = wb)
summary(m1_wb)$coefficients["lid", ]
##                                                                                                                                                   Estimate 
##  1.1759441012001290438604428345570340752601623535156250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 
##                                                                                                                                                 Std. Error 
##  0.0410682873554084892919746607731212861835956573486328125000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 
##                                                                                                                                                    t value 
## 28.6338724335740550941409310325980186462402343750000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 
##                                                                                                                                                   Pr(>|t|) 
##  0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000007115002
pdatos_wb <- pdata.frame(wb, index = c("iso3c", "year"))
m2_wb_fe <- plm(log(patentes + 1) ~ lid + lpib + lpob, data = pdatos_wb,
                model = "within", effect = "twoways")
m2_wb_re <- plm(log(patentes + 1) ~ lid + lpib + lpob, data = pdatos_wb,
                model = "random", effect = "twoways")
summary(m2_wb_fe)$coefficients
##       Estimate Std. Error   t-value
## lid  0.6650576 0.03997443 16.637075
## lpib 0.7670455 0.07699884  9.961779
## lpob 2.6575502 0.17681280 15.030304
##                                                               Pr(>|t|)
## lid  0.000000000000000000000000000000000000000000000000000000004473928
## lpib 0.000000000000000000000115097655518850853347577513462684563121358
## lpob 0.000000000000000000000000000000000000000000000013284751672421986
phtest(m2_wb_fe, m2_wb_re)
## 
##  Hausman Test
## 
## data:  log(patentes + 1) ~ lid + lpib + lpob
## chisq = 282.18, df = 3, p-value < 0.00000000000000022
## alternative hypothesis: one model is inconsistent

A diferencia de la Parte 1, aquí el test de Hausman rechaza con claridad (p cercano a cero). Se usan efectos fijos sin ambigüedad.

m3_wb_poisson <- glm(patentes ~ lid + lpib + lpob + factor(iso3c) + factor(year),
                     family = poisson, data = wb)
coeftest(m3_wb_poisson, vcov. = vcovCL(m3_wb_poisson, cluster = wb$iso3c))["lid", ]
##   Estimate Std. Error    z value   Pr(>|z|) 
##  0.2108873  0.3600017  0.5857953  0.5580131
m4_wb_nb <- glm.nb(patentes ~ lid + lpib + lpob + factor(iso3c) + factor(year),
                   data = wb, control = glm.control(maxit = 100))
coeftest(m4_wb_nb, vcov. = vcovCL(m4_wb_nb, cluster = wb$iso3c))["lid", ]
##       Estimate     Std. Error        z value       Pr(>|z|) 
## 0.677305279256 0.145376535423 4.658972490205 0.000003177917
c(AIC_poisson = AIC(m3_wb_poisson), AIC_binomial_negativa = AIC(m4_wb_nb))
##           AIC_poisson AIC_binomial_negativa 
##             634761.96              22252.33

Con 69 variables dummy de país, el estimador de binomial negativa puede presentar sesgo por parámetros incidentales. Se reporta con esa reserva; la conclusión principal de esta parte se apoya en las regresiones por década.

3.3 Elasticidad I+D-patentes por década

elasticidad_decada <- data.frame()
for (dc in sort(unique(wb$decada))) {
  sub <- wb %>% filter(decada == dc)
  if (n_distinct(sub$iso3c) < 10) next
  m <- plm(log(patentes + 1) ~ lid + lpib + lpob,
           data = pdata.frame(sub, index = c("iso3c", "year")), model = "within")
  elasticidad_decada <- rbind(elasticidad_decada, data.frame(
    decada      = dc,
    elasticidad = summary(m)$coefficients["lid", "Estimate"],
    p           = summary(m)$coefficients["lid", "Pr(>|t|)"]
  ))
}
elasticidad_decada
##   decada elasticidad                  p
## 1  1990s  0.43323070 0.0002854830940060
## 2  2000s  0.51927646 0.0000000001030732
## 3  2010s  0.22900362 0.0000289181280036
## 4  2020s  0.02330779 0.9427515410928404

La elasticidad baja de 0.43-0.52 en los noventa y dos mil, a 0.23 en los dos mil diez, y no se distingue de cero en los veinte. Este patrón solo es visible al agrupar por décadas, que es precisamente la razón por la que la consigna pide ese corte: año con año, la señal se pierde en el ruido.

ggplot(elasticidad_decada, aes(decada, elasticidad, group = 1)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
  geom_line() +
  geom_point(size = 3) +
  labs(title = "La elasticidad I+D-patentes se debilita por década",
       x = NULL, y = "Elasticidad estimada") +
  theme_minimal()

3.4 Validación fuera de muestra

entrena_wb <- wb %>% filter(year <= 2016) %>% mutate(t = year - min(year))
prueba_wb  <- wb %>%
  filter(year >= 2017, iso3c %in% unique(entrena_wb$iso3c)) %>%
  mutate(t = year - min(entrena_wb$year))

# se incluye una tendencia lineal además de los efectos fijos; sin ella el
# modelo no tiene con qué proyectar hacia adelante
m_pred_wb <- lm(log(patentes + 1) ~ lid + lpib + lpob + t + factor(iso3c),
                data = entrena_wb)
prueba_wb$prediccion <- exp(predict(m_pred_wb, newdata = prueba_wb)) - 1

valor_2016 <- entrena_wb %>% filter(year == 2016) %>% select(iso3c, base = patentes)
prueba_wb  <- prueba_wb %>% left_join(valor_2016, by = "iso3c")

data.frame(
  enfoque = c("último valor observado (2016)", "modelo con efectos fijos"),
  RMSE    = c(rmse(prueba_wb$patentes, prueba_wb$base),
              rmse(prueba_wb$patentes, prueba_wb$prediccion))
)
##                         enfoque      RMSE
## 1 último valor observado (2016)  19046.95
## 2      modelo con efectos fijos 129935.27

La referencia ingenua predice mejor que el modelo. El conteo de patentes de un país es muy persistente año con año, de modo que el último valor observado ya contiene casi toda la información predictiva disponible.

Los modelos de esta parte sirven para explicar la relación entre gasto en I+D y patentes, y cómo esa relación cambió en el tiempo, no para pronosticar niveles. Reportar únicamente el R² dentro de muestra daría la impresión contraria.

3.5 México en contexto

wb %>%
  filter(iso3c %in% c("MEX", "BRA", "CHL", "KOR", "ESP"),
         year %in% c(2000, 2010, 2021)) %>%
  select(country, year, patentes, pat_millon, id_pib) %>%
  arrange(country, year)
##        country year patentes  pat_millon  id_pib
## 1       Brazil 2000     3179   18.268196 1.04752
## 2       Brazil 2010     4228   21.827351 1.15992
## 3       Brazil 2021     4666   22.266731 1.13371
## 4        Chile 2010      328   19.090341 0.33165
## 5        Chile 2021      402   20.661652 0.36063
## 6  Korea, Rep. 2000    72831 1549.328370 2.04940
## 7  Korea, Rep. 2010   131805 2659.819633 3.17913
## 8  Korea, Rep. 2021   186245 3597.578877 4.59673
## 9       Mexico 2000      431    4.370064 0.29205
## 10      Mexico 2010      951    8.369718 0.47353
## 11      Mexico 2021     1117    8.750617 0.27305
## 12       Spain 2000     2710   66.801644 0.88315
## 13       Spain 2010     3566   76.561562 1.35436
## 14       Spain 2021     1308   27.569449 1.39617
mexico <- wb %>% filter(iso3c == "MEX") %>% arrange(year)
c(id_pib_2010 = mexico$id_pib[mexico$year == 2010],
  id_pib_2021 = mexico$id_pib[mexico$year == 2021])
## id_pib_2010 id_pib_2021 
##     0.47353     0.27305

México pasa de 0.47% del PIB en I+D en 2010 a 0.27% en 2021, una caída de 42%. Es la mayor reducción entre los cinco países comparados, y con 8.8 patentes por millón de habitantes está por debajo de Brasil y Chile.

ggplot(wb %>% filter(iso3c %in% c("MEX", "BRA", "CHL", "KOR", "ESP")),
       aes(year, id_pib, color = country)) +
  geom_line(linewidth = 0.9) +
  labs(title = "Gasto en investigación y desarrollo como porcentaje del PIB",
       x = NULL, y = "I+D / PIB", color = NULL) +
  theme_minimal()

4 Conclusiones

Parte 1 (Patentes). La variable de solicitudes está truncada en 2020-2021, de modo que el análisis usa patentes otorgadas. La elasticidad del gasto en I+D cae de 0.66 a 0.27 al pasar a efectos fijos y deja de ser significativa con errores robustos: la asociación entre I+D y patentes es fuerte al comparar firmas, pero no se sostiene dentro de una misma firma en el corto plazo. El mejor modelo es binomial negativa, no MCO ni Poisson.

Parte 2 (Cuidado de la piel). Las filas de la base no son aditivas: Men's Grooming es transversal y Prestige está contenido en Premium. El crecimiento debe medirse en términos reales, porque con 4.5% de inflación anual Depilatories pasa de aparentar crecimiento a contraerse. La recomendación de inversión es Skin Care, con la salvedad de que el ganador cambia según los pesos asignados a crecimiento, riesgo y escala.

Parte 3 (Banco Mundial). La relación entre gasto en I+D y patentes se debilita década tras década hasta desaparecer en los años veinte, patrón visible solo al agrupar por décadas. Ningún modelo supera al pronóstico ingenuo, lo que delimita el uso de estos modelos a la explicación y no a la predicción. México redujo su gasto en I+D 42% desde 2010, la mayor caída entre países comparables.