Actividad 1

Importar paquetes y llamar librerías

#install.packages("readxl")
library(readxl)
#install.packages("plm")
library(plm)
#install.packages("gplots")
library(gplots)
library(tidyverse)
library(plm)
library(gplots)
#install.packages("WDI")
library(WDI)
#install.packages("wbstats")
library(wbstats)
library(lmtest)
library(sandwich)

Parte 1. Patentes

#¿Cómo se relaciona el aumento en la inversión en Research & Development con el número de patentes concedidas por las empresas, considerando que el efecto de la inversión puede presentarse hasta tres años después?

df1 <- read.csv("/Users/carlalievanoespinosa/Desktop/business analytics/8vo semestre/df1.csv")

df1 <- df1 %>%
  arrange(cusip, year)
df1_panel <- pdata.frame(
  df1,
  index = c("cusip", "year")
)

pdim(df1_panel)

## Balanced Panel: n = 226, T = 10, N = 2260

#Rezago de 3 años
df1_panel$rnd_lag1 <- stats::lag(
  df1_panel$rndeflt,
  k = 1
)

df1_panel$rnd_lag2 <- stats::lag(
  df1_panel$rndeflt,
  k = 2
)

df1_panel$rnd_lag3 <- stats::lag(
  df1_panel$rndeflt,
  k = 3
)

head(
  df1_panel[
    ,
    c(
      "cusip",
      "year",
      "rndeflt",
      "rnd_lag1",
      "rnd_lag2",
      "rnd_lag3"
    )
  ],
  15
)

##           cusip year rndeflt rnd_lag1 rnd_lag2 rnd_lag3
## 800-2012    800 2012       3       NA       NA       NA
## 800-2013    800 2013       3        3       NA       NA
## 800-2014    800 2014       3        3        3       NA
## 800-2015    800 2015       3        3        3        3
## 800-2016    800 2016       3        3        3        3
## 800-2017    800 2017       3        3        3        3
## 800-2018    800 2018       3        3        3        3
## 800-2019    800 2019       3        3        3        3
## 800-2020    800 2020       4        3        3        3
## 800-2021    800 2021       4        4        3        3
## 4626-2012  4626 2012       1       NA       NA       NA
## 4626-2013  4626 2013       1        1       NA       NA
## 4626-2014  4626 2014       1        1        1       NA
## 4626-2015  4626 2015       1        1        1        1
## 4626-2016  4626 2016       1        1        1        1

variables_modelo <- c(
  "patents",
  "patentsg",
  "rndeflt",
  "rnd_lag1",
  "rnd_lag2",
  "rnd_lag3"
)

df_model <- df1_panel[
  complete.cases(df1_panel[, variables_modelo]),
]

pdim(df_model)

## Balanced Panel: n = 226, T = 7, N = 1582

sort(unique(df_model$year))

## [1] 2015 2016 2017 2018 2019 2020 2021
## Levels: 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021

##PARTE A. R&D → PATENTES SOLICITADAS

# OPCIÓN 1 - MODELO POOLED

pooled_patents <- plm(
  patents ~ rndeflt +
    rnd_lag1 +
    rnd_lag2 +
    rnd_lag3,
  data = df_model,
  model = "pooling"
)

summary(pooled_patents)

## Pooling Model
## 
## Call:
## plm(formula = patents ~ rndeflt + rnd_lag1 + rnd_lag2 + rnd_lag3, 
##     data = df_model, model = "pooling")
## 
## Balanced Panel: n = 226, T = 7, N = 1582
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -407.1786  -11.5903  -10.3097   -5.0386  715.2093 
## 
## Coefficients:
##             Estimate Std. Error t-value  Pr(>|t|)    
## (Intercept) 11.59027    1.39586  8.3033 < 2.2e-16 ***
## rndeflt      0.34591    0.20470  1.6899  0.091250 .  
## rnd_lag1     0.18937    0.32489  0.5829  0.560063    
## rnd_lag2    -0.91870    0.30137 -3.0484  0.002339 ** 
## rnd_lag3     0.90251    0.20228  4.4616 8.711e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    7131900
## Residual Sum of Squares: 4577600
## R-Squared:      0.35815
## Adj. R-Squared: 0.35652
## F-statistic: 219.992 on 4 and 1577 DF, p-value: < 2.22e-16

# OPCIÓN 2 - EFECTOS FIJOS

within_patents <- plm(
  patents ~ rndeflt +
    rnd_lag1 +
    rnd_lag2 +
    rnd_lag3,
  data = df_model,
  model = "within",
  effect = "individual"
)

summary(within_patents)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = patents ~ rndeflt + rnd_lag1 + rnd_lag2 + rnd_lag3, 
##     data = df_model, effect = "individual", model = "within")
## 
## Balanced Panel: n = 226, T = 7, N = 1582
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -488.46648   -1.26762   -0.14286    1.60287  175.72053 
## 
## Coefficients:
##           Estimate Std. Error t-value  Pr(>|t|)    
## rndeflt  -0.424555   0.102322 -4.1492 3.544e-05 ***
## rnd_lag1 -0.054727   0.149961 -0.3649    0.7152    
## rnd_lag2 -0.205098   0.141414 -1.4503    0.1472    
## rnd_lag3 -0.572877   0.106157 -5.3965 8.014e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    908870
## Residual Sum of Squares: 774890
## R-Squared:      0.14742
## Adj. R-Squared: 0.00301
## F-statistic: 58.4433 on 4 and 1352 DF, p-value: < 2.22e-16

# PRUEBA F
# H0: Pooled es suficiente
# H1: existen efectos individuales por empresa

prueba_f_patents <- pFtest(
  within_patents,
  pooled_patents
)

prueba_f_patents

## 
##  F test for individual effects
## 
## data:  patents ~ rndeflt + rnd_lag1 + rnd_lag2 + rnd_lag3
## F = 29.488, df1 = 225, df2 = 1352, p-value < 2.2e-16
## alternative hypothesis: significant effects

# OPCIÓN 3 - EFECTOS ALEATORIOS

random_patents <- plm(
  patents ~ rndeflt +
    rnd_lag1 +
    rnd_lag2 +
    rnd_lag3,
  data = df_model,
  model = "random",
  effect = "individual"
)

summary(random_patents)

## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = patents ~ rndeflt + rnd_lag1 + rnd_lag2 + rnd_lag3, 
##     data = df_model, effect = "individual", model = "random")
## 
## Balanced Panel: n = 226, T = 7, N = 1582
## 
## Effects:
##                   var std.dev share
## idiosyncratic  573.14   23.94 0.362
## individual    1010.67   31.79 0.638
## theta: 0.7262
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -390.99395   -4.63077   -3.74683   -0.79097  326.02391 
## 
## Coefficients:
##              Estimate Std. Error z-value Pr(>|z|)    
## (Intercept) 14.943717   2.686516  5.5625 2.66e-08 ***
## rndeflt      0.207926   0.115444  1.8011  0.07169 .  
## rnd_lag1    -0.141619   0.177775 -0.7966  0.42567    
## rnd_lag2    -0.019263   0.166886 -0.1154  0.90811    
## rnd_lag3     0.294035   0.116054  2.5336  0.01129 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    1375200
## Residual Sum of Squares: 1280200
## R-Squared:      0.069071
## Adj. R-Squared: 0.066709
## Chisq: 117.006 on 4 DF, p-value: < 2.22e-16

# PRUEBA LM DE BREUSCH-PAGAN
# H0: no existen efectos individuales aleatorios
# H1: Random es preferible a Pooled

prueba_lm_patents <- plmtest(
  pooled_patents,
  effect = "individual",
  type = "bp"
)

prueba_lm_patents

## 
##  Lagrange Multiplier Test - (Breusch-Pagan)
## 
## data:  patents ~ rndeflt + rnd_lag1 + rnd_lag2 + rnd_lag3
## chisq = 2479.5, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effects

# PRUEBA DE HAUSMAN

prueba_hausman_patents <- phtest(
  within_patents,
  random_patents
)

prueba_hausman_patents

## 
##  Hausman Test
## 
## data:  patents ~ rndeflt + rnd_lag1 + rnd_lag2 + rnd_lag3
## chisq = 430.84, df = 4, p-value < 2.2e-16
## alternative hypothesis: one model is inconsistent

if(
  prueba_f_patents$p.value >= 0.05 &
  prueba_lm_patents$p.value >= 0.05
){
  
  modelo_patents_nombre <- "Pooled"
  modelo_patents <- pooled_patents
  
} else if(
  prueba_f_patents$p.value < 0.05 &
  prueba_lm_patents$p.value >= 0.05
){
  
  modelo_patents_nombre <- "Fixed Effects"
  modelo_patents <- within_patents
  
} else if(
  prueba_f_patents$p.value >= 0.05 &
  prueba_lm_patents$p.value < 0.05
){
  
  modelo_patents_nombre <- "Random Effects"
  modelo_patents <- random_patents
  
} else {
  
  if(prueba_hausman_patents$p.value < 0.05){
    
    modelo_patents_nombre <- "Fixed Effects"
    modelo_patents <- within_patents
    
  } else {
    
    modelo_patents_nombre <- "Random Effects"
    modelo_patents <- random_patents
    
  }
}

modelo_patents_nombre

## [1] "Fixed Effects"

summary(modelo_patents)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = patents ~ rndeflt + rnd_lag1 + rnd_lag2 + rnd_lag3, 
##     data = df_model, effect = "individual", model = "within")
## 
## Balanced Panel: n = 226, T = 7, N = 1582
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -488.46648   -1.26762   -0.14286    1.60287  175.72053 
## 
## Coefficients:
##           Estimate Std. Error t-value  Pr(>|t|)    
## rndeflt  -0.424555   0.102322 -4.1492 3.544e-05 ***
## rnd_lag1 -0.054727   0.149961 -0.3649    0.7152    
## rnd_lag2 -0.205098   0.141414 -1.4503    0.1472    
## rnd_lag3 -0.572877   0.106157 -5.3965 8.014e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    908870
## Residual Sum of Squares: 774890
## R-Squared:      0.14742
## Adj. R-Squared: 0.00301
## F-statistic: 58.4433 on 4 and 1352 DF, p-value: < 2.22e-16

##PARTE B. R&D → PATENTES CONCEDIDAS

# OPCIÓN 1 - POOLED

pooled_granted <- plm(
  patentsg ~ rndeflt +
    rnd_lag1 +
    rnd_lag2 +
    rnd_lag3,
  data = df_model,
  model = "pooling"
)

summary(pooled_granted)

## Pooling Model
## 
## Call:
## plm(formula = patentsg ~ rndeflt + rnd_lag1 + rnd_lag2 + rnd_lag3, 
##     data = df_model, model = "pooling")
## 
## Balanced Panel: n = 226, T = 7, N = 1582
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -450.6105  -13.8552  -11.8552   -4.6872  673.3196 
## 
## Coefficients:
##              Estimate Std. Error t-value  Pr(>|t|)    
## (Intercept) 13.855242   1.437257  9.6401 < 2.2e-16 ***
## rndeflt     -0.037879   0.210769 -0.1797 0.8573951    
## rnd_lag1     0.621334   0.334529  1.8573 0.0634493 .  
## rnd_lag2    -1.116372   0.310310 -3.5976 0.0003311 ***
## rnd_lag3     1.158833   0.208283  5.5638 3.097e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    8365600
## Residual Sum of Squares: 4853100
## R-Squared:      0.41987
## Adj. R-Squared: 0.4184
## F-statistic: 285.337 on 4 and 1577 DF, p-value: < 2.22e-16

# OPCIÓN 2 - EFECTOS FIJOS

within_granted <- plm(
  patentsg ~ rndeflt +
    rnd_lag1 +
    rnd_lag2 +
    rnd_lag3,
  data = df_model,
  model = "within",
  effect = "individual"
)

summary(within_granted)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = patentsg ~ rndeflt + rnd_lag1 + rnd_lag2 + rnd_lag3, 
##     data = df_model, effect = "individual", model = "within")
## 
## Balanced Panel: n = 226, T = 7, N = 1582
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -172.66694   -1.45774   -0.20473    1.46162  125.35634 
## 
## Coefficients:
##           Estimate Std. Error t-value  Pr(>|t|)    
## rndeflt  -0.526476   0.054287 -9.6980 < 2.2e-16 ***
## rnd_lag1  0.357361   0.079562  4.4916 7.671e-06 ***
## rnd_lag2 -0.204800   0.075028 -2.7297  0.006422 ** 
## rnd_lag3  0.184111   0.056322  3.2689  0.001107 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    241040
## Residual Sum of Squares: 218120
## R-Squared:      0.095086
## Adj. R-Squared: -0.058187
## F-statistic: 35.5164 on 4 and 1352 DF, p-value: < 2.22e-16

#Prueba F
prueba_f_granted <- pFtest(
  within_granted,
  pooled_granted
)

prueba_f_granted

## 
##  F test for individual effects
## 
## data:  patentsg ~ rndeflt + rnd_lag1 + rnd_lag2 + rnd_lag3
## F = 127.69, df1 = 225, df2 = 1352, p-value < 2.2e-16
## alternative hypothesis: significant effects

# OPCIÓN 3 - EFECTOS ALEATORIOS

random_granted <- plm(
  patentsg ~ rndeflt +
    rnd_lag1 +
    rnd_lag2 +
    rnd_lag3,
  data = df_model,
  model = "random",
  effect = "individual"
)

summary(random_granted)

## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = patentsg ~ rndeflt + rnd_lag1 + rnd_lag2 + rnd_lag3, 
##     data = df_model, effect = "individual", model = "random")
## 
## Balanced Panel: n = 226, T = 7, N = 1582
## 
## Effects:
##                   var std.dev share
## idiosyncratic  161.33   12.70 0.105
## individual    1375.90   37.09 0.895
## theta: 0.8716
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -105.63542   -3.34358   -2.47193   -0.17583  163.19022 
## 
## Coefficients:
##              Estimate Std. Error z-value  Pr(>|z|)    
## (Intercept) 19.258655   2.893187  6.6566 2.803e-11 ***
## rndeflt     -0.312896   0.059729 -5.2386 1.618e-07 ***
## rnd_lag1     0.323558   0.090629  3.5701 0.0003568 ***
## rnd_lag2    -0.133673   0.085224 -1.5685 0.1167662    
## rnd_lag3     0.469500   0.060632  7.7434 9.678e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    374890
## Residual Sum of Squares: 330960
## R-Squared:      0.11719
## Adj. R-Squared: 0.11495
## Chisq: 209.346 on 4 DF, p-value: < 2.22e-16

prueba_lm_granted <- plmtest(
  pooled_granted,
  effect = "individual",
  type = "bp"
)

prueba_lm_granted

## 
##  Lagrange Multiplier Test - (Breusch-Pagan)
## 
## data:  patentsg ~ rndeflt + rnd_lag1 + rnd_lag2 + rnd_lag3
## chisq = 4000.6, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effects

#Prueba Hausman
prueba_hausman_granted <- phtest(
  within_granted,
  random_granted
)

prueba_hausman_granted

## 
##  Hausman Test
## 
## data:  patentsg ~ rndeflt + rnd_lag1 + rnd_lag2 + rnd_lag3
## chisq = 265.5, df = 4, p-value < 2.2e-16
## alternative hypothesis: one model is inconsistent

if(
  prueba_f_granted$p.value >= 0.05 &
  prueba_lm_granted$p.value >= 0.05
){
  
  modelo_final_nombre <- "Pooled"
  modelo_final <- pooled_granted
  
} else if(
  prueba_f_granted$p.value < 0.05 &
  prueba_lm_granted$p.value >= 0.05
){
  
  modelo_final_nombre <- "Fixed Effects"
  modelo_final <- within_granted
  
} else if(
  prueba_f_granted$p.value >= 0.05 &
  prueba_lm_granted$p.value < 0.05
){
  
  modelo_final_nombre <- "Random Effects"
  modelo_final <- random_granted
  
} else {
  
  if(prueba_hausman_granted$p.value < 0.05){
    
    modelo_final_nombre <- "Fixed Effects"
    modelo_final <- within_granted
    
  } else {
    
    modelo_final_nombre <- "Random Effects"
    modelo_final <- random_granted
    
  }
}

modelo_final_nombre

## [1] "Fixed Effects"

summary(modelo_final)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = patentsg ~ rndeflt + rnd_lag1 + rnd_lag2 + rnd_lag3, 
##     data = df_model, effect = "individual", model = "within")
## 
## Balanced Panel: n = 226, T = 7, N = 1582
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -172.66694   -1.45774   -0.20473    1.46162  125.35634 
## 
## Coefficients:
##           Estimate Std. Error t-value  Pr(>|t|)    
## rndeflt  -0.526476   0.054287 -9.6980 < 2.2e-16 ***
## rnd_lag1  0.357361   0.079562  4.4916 7.671e-06 ***
## rnd_lag2 -0.204800   0.075028 -2.7297  0.006422 ** 
## rnd_lag3  0.184111   0.056322  3.2689  0.001107 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    241040
## Residual Sum of Squares: 218120
## R-Squared:      0.095086
## Adj. R-Squared: -0.058187
## F-statistic: 35.5164 on 4 and 1352 DF, p-value: < 2.22e-16

resultado_robusto <- coeftest(
  modelo_final,
  vcov = vcovHC(
    modelo_final,
    method = "arellano",
    type = "HC1",
    cluster = "group"
  )
)

resultado_robusto

## 
## t test of coefficients:
## 
##          Estimate Std. Error t value Pr(>|t|)   
## rndeflt  -0.52648    0.16145 -3.2609 0.001138 **
## rnd_lag1  0.35736    0.22267  1.6049 0.108755   
## rnd_lag2 -0.20480    0.20870 -0.9813 0.326613   
## rnd_lag3  0.18411    0.10777  1.7084 0.087799 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Efecto acumulado de R&D durante los cuatro años

coef_rnd <- coef(modelo_final)[
  c(
    "rndeflt",
    "rnd_lag1",
    "rnd_lag2",
    "rnd_lag3"
  )
]

coef_rnd

##    rndeflt   rnd_lag1   rnd_lag2   rnd_lag3 
## -0.5264760  0.3573614 -0.2048000  0.1841114

efecto_acumulado <- sum(coef_rnd)

efecto_acumulado

## [1] -0.1898031

cor(
  as.data.frame(
    df_model[, c(
      "rndeflt",
      "rnd_lag1",
      "rnd_lag2",
      "rnd_lag3"
    )]
  ),
  use = "complete.obs"
)

##            rndeflt  rnd_lag1  rnd_lag2  rnd_lag3
## rndeflt  1.0000000 0.9951703 0.9866716 0.9849683
## rnd_lag1 0.9951703 1.0000000 0.9947065 0.9883012
## rnd_lag2 0.9866716 0.9947065 1.0000000 0.9946290
## rnd_lag3 0.9849683 0.9883012 0.9946290 1.0000000

library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

## The following object is masked from 'package:purrr':
## 
##     some

linearHypothesis(
  modelo_final,
  "rndeflt + rnd_lag1 + rnd_lag2 + rnd_lag3 = 0",
  vcov. = vcovHC(
    modelo_final,
    method = "arellano",
    type = "HC1",
    cluster = "group"
  )
)

## 
## Linear hypothesis test:
## rndeflt  + rnd_lag1  + rnd_lag2  + rnd_lag3 = 0
## 
## Model 1: restricted model
## Model 2: patentsg ~ rndeflt + rnd_lag1 + rnd_lag2 + rnd_lag3
## 
## Note: Coefficient covariance matrix supplied.
## 
##   Res.Df Df  Chisq Pr(>Chisq)   
## 1   1353                        
## 2   1352  1 9.4683   0.002091 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# ==============================================================
# CORRECCIÓN DE MULTICOLINEALIDAD Y MODELO FINAL
# ==============================================================

# Crear una sola variable que resume el R&D
# de los tres años anteriores

df_model$rnd_prom_3 <- (
  df_model$rnd_lag1 +
  df_model$rnd_lag2 +
  df_model$rnd_lag3
) / 3


# --------------------------------------------------------------
# MODELO FINAL DE EFECTOS FIJOS
# Fixed Effects ya fue seleccionado previamente
# mediante las pruebas F, LM y Hausman
# --------------------------------------------------------------

modelo_final_fe <- plm(
  patentsg ~ rnd_prom_3,
  data = df_model,
  model = "within",
  effect = "individual"
)

summary(modelo_final_fe)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = patentsg ~ rnd_prom_3, data = df_model, effect = "individual", 
##     model = "within")
## 
## Balanced Panel: n = 226, T = 7, N = 1582
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -187.73029   -1.44924   -0.24304    1.38790  127.28454 
## 
## Coefficients:
##             Estimate Std. Error t-value Pr(>|t|)
## rnd_prom_3 -0.056013   0.040439 -1.3851   0.1662
## 
## Total Sum of Squares:    241040
## Residual Sum of Squares: 240700
## R-Squared:      0.0014139
## Adj. R-Squared: -0.16514
## F-statistic: 1.91861 on 1 and 1355 DF, p-value: 0.16624

# --------------------------------------------------------------
# ERRORES ESTÁNDAR ROBUSTOS
# --------------------------------------------------------------

resultado_robusto_final <- coeftest(
  modelo_final_fe,
  vcov = vcovHC(
    modelo_final_fe,
    method = "arellano",
    type = "HC1",
    cluster = "group"
  )
)

resultado_robusto_final

## 
## t test of coefficients:
## 
##             Estimate Std. Error t value Pr(>|t|)
## rnd_prom_3 -0.056013   0.093829  -0.597   0.5506

# ==============================================================
# NUEVA ESPECIFICACIÓN: STOCK DE R&D
# ==============================================================

# Modelo de efectos fijos
# Fixed Effects ya fue seleccionado con F, LM y Hausman

modelo_stock <- plm(
  patentsg ~ rndstck,
  data = df1_panel,
  model = "within",
  effect = "individual"
)

summary(modelo_stock)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = patentsg ~ rndstck, data = df1_panel, effect = "individual", 
##     model = "within")
## 
## Unbalanced Panel: n = 216, T = 2-10, N = 2103
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -228.16717   -2.01585   -0.32629    1.60338  268.63087 
## 
## Coefficients:
##           Estimate Std. Error t-value  Pr(>|t|)    
## rndstck -0.0317122  0.0016961 -18.698 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    713560
## Residual Sum of Squares: 601970
## R-Squared:      0.15638
## Adj. R-Squared: 0.059761
## F-statistic: 349.601 on 1 and 1886 DF, p-value: < 2.22e-16

# Errores estándar robustos

resultado_stock_robusto <- coeftest(
  modelo_stock,
  vcov = vcovHC(
    modelo_stock,
    method = "arellano",
    type = "HC1",
    cluster = "group"
  )
)

resultado_stock_robusto

## 
## t test of coefficients:
## 
##          Estimate Std. Error t value Pr(>|t|)   
## rndstck -0.031712   0.010369 -3.0584 0.002256 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# ============================================================
# VALIDACIÓN DEL MODELO CON STOCK DE R&D
# ============================================================

# Misma muestra para los tres modelos
df_stock <- df1_panel[
  complete.cases(df1_panel[, c("patentsg", "rndstck")]),
]


# 1. Pooled
pooled_stock <- plm(
  patentsg ~ rndstck,
  data = df_stock,
  model = "pooling"
)


# 2. Fixed Effects
fixed_stock <- plm(
  patentsg ~ rndstck,
  data = df_stock,
  model = "within",
  effect = "individual"
)


# 3. Random Effects
random_stock <- plm(
  patentsg ~ rndstck,
  data = df_stock,
  model = "random",
  effect = "individual"
)


# ------------------------------------------------------------
# PRUEBAS DE SELECCIÓN
# ------------------------------------------------------------

# Heterogeneidad: Pooled vs Fixed
pFtest(fixed_stock, pooled_stock)

## 
##  F test for individual effects
## 
## data:  patentsg ~ rndstck
## F = 123.64, df1 = 215, df2 = 1886, p-value < 2.2e-16
## alternative hypothesis: significant effects

# Pooled vs Random
plmtest(
  pooled_stock,
  effect = "individual",
  type = "bp"
)

## 
##  Lagrange Multiplier Test - (Breusch-Pagan)
## 
## data:  patentsg ~ rndstck
## chisq = 5795.9, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effects

# Fixed vs Random
phtest(
  fixed_stock,
  random_stock
)

## 
##  Hausman Test
## 
## data:  patentsg ~ rndstck
## chisq = 273.98, df = 1, p-value < 2.2e-16
## alternative hypothesis: one model is inconsistent

# ------------------------------------------------------------
# EFECTOS FIJOS DE EMPRESA Y AÑO
# ------------------------------------------------------------

fixed_stock_tw <- plm(
  patentsg ~ rndstck,
  data = df_stock,
  model = "within",
  effect = "twoways"
)

summary(fixed_stock_tw)

## Twoways effects Within Model
## 
## Call:
## plm(formula = patentsg ~ rndstck, data = df_stock, effect = "twoways", 
##     model = "within")
## 
## Unbalanced Panel: n = 216, T = 2-10, N = 2103
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -222.10922   -3.12579   -0.58684    2.78151  267.33331 
## 
## Coefficients:
##           Estimate Std. Error t-value  Pr(>|t|)    
## rndstck -0.0298525  0.0017259 -17.297 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    678060
## Residual Sum of Squares: 584840
## R-Squared:      0.13748
## Adj. R-Squared: 0.034092
## F-statistic: 299.191 on 1 and 1877 DF, p-value: < 2.22e-16

# Errores robustos
resultado_stock_tw <- coeftest(
  fixed_stock_tw,
  vcov = vcovHC(
    fixed_stock_tw,
    method = "arellano",
    type = "HC1",
    cluster = "group"
  )
)

resultado_stock_tw

## 
## t test of coefficients:
## 
##          Estimate Std. Error t value Pr(>|t|)   
## rndstck -0.029853   0.010193 -2.9286 0.003446 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# ==============================================================
# PREDICCIÓN DE PATENTES CONCEDIDAS PARA 2022
# ==============================================================

library(dplyr)
library(plm)

# --------------------------------------------------------------
# 1. Preparar datos para el modelo predictivo
# --------------------------------------------------------------

df_pred_model <- df1 %>%
  mutate(
    year = as.numeric(as.character(year))
  ) %>%
  filter(
    !is.na(patentsg),
    !is.na(rndstck)
  ) %>%
  arrange(cusip, year)

df_pred_panel <- pdata.frame(
  df_pred_model,
  index = c("cusip", "year")
)


# --------------------------------------------------------------
# 2. Modelo predictivo
# Efectos fijos por empresa + tendencia temporal
# --------------------------------------------------------------

modelo_pred <- plm(
  patentsg ~ rndstck + year,
  data = df_pred_panel,
  model = "within",
  effect = "individual"
)

summary(modelo_pred)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = patentsg ~ rndstck + year, data = df_pred_panel, 
##     effect = "individual", model = "within")
## 
## Unbalanced Panel: n = 216, T = 2-10, N = 2103
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -222.10922   -3.12579   -0.58684    2.78151  267.33331 
## 
## Coefficients:
##             Estimate  Std. Error  t-value  Pr(>|t|)    
## rndstck   -0.0298525   0.0017259 -17.2971 < 2.2e-16 ***
## year2013  -0.2894251   1.7532973  -0.1651  0.868903    
## year2014  -0.2517885   1.7522766  -0.1437  0.885759    
## year2015  -2.0185134   1.7544393  -1.1505  0.250077    
## year2016  -1.9026020   1.7580438  -1.0822  0.279291    
## year2017  -4.0596771   1.7612893  -2.3049  0.021278 *  
## year2018  -3.4622903   1.7698353  -1.9563  0.050580 .  
## year2019 -10.2173979   1.7744538  -5.7581 9.919e-09 ***
## year2020  -5.1674518   1.7811633  -2.9012  0.003761 ** 
## year2021  -2.5948339   1.7891830  -1.4503  0.147145    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    713560
## Residual Sum of Squares: 584840
## R-Squared:      0.18039
## Adj. R-Squared: 0.082142
## F-statistic: 41.3115 on 10 and 1877 DF, p-value: < 2.22e-16

# --------------------------------------------------------------
# 3. Proyectar rndstck de cada empresa para 2022
# usando su tendencia histórica
# --------------------------------------------------------------

rndstck_2022 <- df_pred_model %>%
  group_by(cusip) %>%
  filter(sum(!is.na(rndstck)) >= 2) %>%
  group_modify(~{
    
    modelo_rnd <- lm(
      rndstck ~ year,
      data = .x
    )
    
    tibble(
      year = 2022,
      rndstck_pred = predict(
        modelo_rnd,
        newdata = data.frame(year = 2022)
      )
    )
    
  }) %>%
  ungroup()


# --------------------------------------------------------------
# 4. Extraer coeficientes del modelo
# --------------------------------------------------------------

beta_rnd <- unname(
  coef(modelo_pred)["rndstck"]
)

beta_year <- unname(
  coef(modelo_pred)["year"]
)

efectos_empresa <- fixef(
  modelo_pred,
  type = "level"
)


# --------------------------------------------------------------
# 5. Predecir patentes concedidas para 2022
# --------------------------------------------------------------

prediccion_2022 <- rndstck_2022 %>%
  mutate(
    
    efecto_empresa = unname(
      efectos_empresa[as.character(cusip)]
    ),
    
    patentsg_pred =
      efecto_empresa +
      beta_rnd * rndstck_pred +
      beta_year * 2022,
    
    # Las patentes no pueden ser negativas
    patentsg_pred = pmax(
      0,
      patentsg_pred
    ),
    
    # Redondear a número entero de patentes
    patentsg_pred_round = round(
      patentsg_pred
    )
    
  ) %>%
  filter(
    !is.na(patentsg_pred)
  ) %>%
  select(
    cusip,
    year,
    rndstck_pred,
    patentsg_pred,
    patentsg_pred_round
  ) %>%
  arrange(
    desc(patentsg_pred)
  )


# --------------------------------------------------------------
# 6. Resultados de la predicción 2022
# --------------------------------------------------------------

prediccion_2022

## # A tibble: 0 × 5
## # ℹ 5 variables: cusip <int>, year <dbl>, rndstck_pred <dbl>,
## #   patentsg_pred <dbl>, patentsg_pred_round <dbl>

# Primeras 20 empresas con mayor predicción
head(prediccion_2022, 20)

## # A tibble: 0 × 5
## # ℹ 5 variables: cusip <int>, year <dbl>, rndstck_pred <dbl>,
## #   patentsg_pred <dbl>, patentsg_pred_round <dbl>

# Resumen de predicciones
summary(prediccion_2022$patentsg_pred)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 

# Total de patentes concedidas estimadas para 2022
sum(prediccion_2022$patentsg_pred_round)

## [1] 0

#CONCLUSIÓN: Los resultados muestran que existen diferencias significativas entre las empresas, por lo que el modelo de efectos fijos fue el más adecuado para analizar la relación entre la inversión en R&D y las patentes concedidas. Después de corregir el problema de multicolinealidad mediante el uso del stock acumulado de R&D (rndstck), se encontró una relación negativa y estadísticamente significativa con las patentes concedidas ((), (p=0.0034)), incluso al controlar por diferencias entre empresas y entre años. Esto indica que, dentro de una misma empresa, un aumento en el stock de R&D no se refleja necesariamente en un aumento inmediato de las patentes concedidas. Por lo tanto, los resultados no respaldan la hipótesis inicial de una relación positiva directa y sugieren que el proceso mediante el cual la inversión en investigación se convierte en patentes puede depender de otros factores y de periodos de tiempo más amplios.

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.

df2 <- read.csv("/Users/carlalievanoespinosa/Desktop/business analytics/8vo semestre/df2.csv")

df2 <- read.csv(
  "/Users/carlalievanoespinosa/Desktop/business analytics/8vo semestre/df2.csv",
  skip = 1,
  check.names = FALSE,
  na.strings = c("-", "")
  )

# Subcategorías
categorias <- c(
  "Bath and Shower",
  "Deodorants",
  "Depilatories",
  "Fragrances",
  "Hair Care",
  "Men's Grooming",
  "Skin Care",
  "Sun Care"
)

# Pasar de formato ancho a formato largo
df2_long <- df2 %>%
  select(Category, all_of(as.character(2011:2025))) %>%
  pivot_longer(
    cols = all_of(as.character(2011:2025)),
    names_to = "year",
    values_to = "market_size"
  ) %>%
  mutate(
    year = as.integer(year),
    market_size = parse_number(as.character(market_size))
  )


totales <- df2_long %>%
  filter(Category == "Beauty and Personal Care") %>%
  select(year, total_market = market_size)

# Calcular participación de mercado
df2_model <- df2_long %>%
  filter(Category %in% categorias) %>%
  left_join(totales, by = "year") %>%
  mutate(
    market_share = (market_size / total_market) * 100,
    trend = year - 2011,
    category = factor(Category)
  ) %>%
  arrange(category, year)

# Revisar
head(df2_model)

## # A tibble: 6 × 7
##   Category         year market_size total_market market_share trend category    
##   <chr>           <int>       <dbl>        <dbl>        <dbl> <dbl> <fct>       
## 1 Bath and Shower  2011       8410.      126329.         6.66     0 Bath and Sh…
## 2 Bath and Shower  2012       9085       136020.         6.68     1 Bath and Sh…
## 3 Bath and Shower  2013       9711.      142539.         6.81     2 Bath and Sh…
## 4 Bath and Shower  2014      10227.      149165.         6.86     3 Bath and Sh…
## 5 Bath and Shower  2015      10815.      157656.         6.86     4 Bath and Sh…
## 6 Bath and Shower  2016      11482.      170488.         6.73     5 Bath and Sh…

# Panel
df2_panel <- pdata.frame(
  as.data.frame(df2_model),
  index = c("category", "year")
)

# Revisar estructura del panel
pdim(df2_panel)

## Balanced Panel: n = 8, T = 15, N = 120

# Comprobar si es balanceado
is.pbalanced(df2_panel)

## [1] TRUE

# PRUEBA DE HETEROGENEIDAD

par(mar = c(8, 4, 3, 1))

plotmeans(
  market_share ~ category,
  data = df2_model,
  las = 2,
  xlab = "",
  ylab = "Participación de mercado (%)"
)

## Warning in arrows(x, li, x, pmax(y - gap, li), col = barcol, lwd = lwd, :
## zero-length arrow is of indeterminate angle and so skipped
## Warning in arrows(x, li, x, pmax(y - gap, li), col = barcol, lwd = lwd, :
## zero-length arrow is of indeterminate angle and so skipped
## Warning in arrows(x, li, x, pmax(y - gap, li), col = barcol, lwd = lwd, :
## zero-length arrow is of indeterminate angle and so skipped
## Warning in arrows(x, li, x, pmax(y - gap, li), col = barcol, lwd = lwd, :
## zero-length arrow is of indeterminate angle and so skipped
## Warning in arrows(x, li, x, pmax(y - gap, li), col = barcol, lwd = lwd, :
## zero-length arrow is of indeterminate angle and so skipped
## Warning in arrows(x, li, x, pmax(y - gap, li), col = barcol, lwd = lwd, :
## zero-length arrow is of indeterminate angle and so skipped
## Warning in arrows(x, li, x, pmax(y - gap, li), col = barcol, lwd = lwd, :
## zero-length arrow is of indeterminate angle and so skipped
## Warning in arrows(x, li, x, pmax(y - gap, li), col = barcol, lwd = lwd, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(x, ui, x, pmin(y + gap, ui), col = barcol, lwd = lwd, :
## zero-length arrow is of indeterminate angle and so skipped
## Warning in arrows(x, ui, x, pmin(y + gap, ui), col = barcol, lwd = lwd, :
## zero-length arrow is of indeterminate angle and so skipped
## Warning in arrows(x, ui, x, pmin(y + gap, ui), col = barcol, lwd = lwd, :
## zero-length arrow is of indeterminate angle and so skipped
## Warning in arrows(x, ui, x, pmin(y + gap, ui), col = barcol, lwd = lwd, :
## zero-length arrow is of indeterminate angle and so skipped
## Warning in arrows(x, ui, x, pmin(y + gap, ui), col = barcol, lwd = lwd, :
## zero-length arrow is of indeterminate angle and so skipped
## Warning in arrows(x, ui, x, pmin(y + gap, ui), col = barcol, lwd = lwd, :
## zero-length arrow is of indeterminate angle and so skipped
## Warning in arrows(x, ui, x, pmin(y + gap, ui), col = barcol, lwd = lwd, :
## zero-length arrow is of indeterminate angle and so skipped
## Warning in arrows(x, ui, x, pmin(y + gap, ui), col = barcol, lwd = lwd, :
## zero-length arrow is of indeterminate angle and so skipped

# OPCIÓN 1 - MODELO DE REGRESIÓN AGRUPADA (POOLED)

pooled <- plm(
  market_share ~ trend,
  data = df2_panel,
  model = "pooling"
)

summary(pooled)

## Pooling Model
## 
## Call:
## plm(formula = market_share ~ trend, data = df2_panel, model = "pooling")
## 
## Balanced Panel: n = 8, T = 15, N = 120
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -10.616051  -5.734911  -0.098855   6.435351  12.001541 
## 
## Coefficients:
##              Estimate Std. Error t-value Pr(>|t|)    
## (Intercept) 10.488539   1.292050  8.1177 5.16e-13 ***
## trend        0.051766   0.157070  0.3296   0.7423    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    6527
## Residual Sum of Squares: 6521
## R-Squared:      0.00091964
## Adj. R-Squared: -0.0075471
## F-statistic: 0.108618 on 1 and 118 DF, p-value: 0.74231

# OPCIÓN 2 - MODELO DE EFECTOS FIJOS (WITHIN)

within <- plm(
  market_share ~ trend,
  data = df2_panel,
  model = "within",
  effect = "individual"
)

summary(within)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = market_share ~ trend, data = df2_panel, effect = "individual", 
##     model = "within")
## 
## Balanced Panel: n = 8, T = 15, N = 120
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -1.8250909 -0.3761985 -0.0040966  0.2783801  2.6852699 
## 
## Coefficients:
##       Estimate Std. Error t-value Pr(>|t|)   
## trend 0.051766   0.016069  3.2215 0.001674 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    70.204
## Residual Sum of Squares: 64.202
## R-Squared:      0.085501
## Adj. R-Squared: 0.019591
## F-statistic: 10.3779 on 1 and 111 DF, p-value: 0.0016738

# PRUEBA F

prueba_f <- pFtest(within, pooled)

prueba_f

## 
##  F test for individual effects
## 
## data:  market_share ~ trend
## F = 1594.8, df1 = 7, df2 = 111, p-value < 2.2e-16
## alternative hypothesis: significant effects

# OPCIÓN 3 - MODELO DE EFECTOS ALEATORIOS (RANDOM)

random <- plm(
  market_share ~ trend,
  data = df2_panel,
  model = "random",
  effect = "individual"
)

summary(random)

## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = market_share ~ trend, data = df2_panel, effect = "individual", 
##     model = "random")
## 
## Balanced Panel: n = 8, T = 15, N = 120
## 
## Effects:
##                   var std.dev share
## idiosyncratic  0.5784  0.7605 0.009
## individual    61.4547  7.8393 0.991
## theta: 0.975
## 
## Residuals:
##     Min.  1st Qu.   Median  3rd Qu.     Max. 
## -1.72003 -0.46571 -0.12213  0.22992  2.79033 
## 
## Coefficients:
##              Estimate Std. Error z-value  Pr(>|z|)    
## (Intercept) 10.488539   2.774764  3.7800 0.0001568 ***
## trend        0.051766   0.016069  3.2215 0.0012753 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    74.253
## Residual Sum of Squares: 68.25
## R-Squared:      0.080839
## Adj. R-Squared: 0.073049
## Chisq: 10.3779 on 1 DF, p-value: 0.0012753

# PRUEBA DE HAUSMAN

prueba_hausman <- phtest(within, random)

prueba_hausman

## 
##  Hausman Test
## 
## data:  market_share ~ trend
## chisq = 1.1842e-15, df = 1, p-value = 1
## alternative hypothesis: one model is inconsistent

if(prueba_f$p.value >= 0.05){
  
  modelo_final <- "pooled"
  
} else {
  
  if(prueba_hausman$p.value < 0.05){
    
    modelo_final <- "within"
    
  } else {
    
    modelo_final <- "random"
    
  }
}

modelo_final

## [1] "random"

# CONCLUSIÓN:
# De acuerdo con la prueba de Hausman, el p-value es mayor a 0.05. Por lo tanto, el modelo de Efectos Aleatorios es más adecuado
# que el modelo de Efectos Fijos.

# Se utilizará el modelo Random para realizar la predicción de la participación de mercado de cada subcategoría para 2026.

# Coeficientes del modelo Random
beta0 <- coef(random)["(Intercept)"]
beta1 <- coef(random)["trend"]

beta0

## (Intercept) 
##    10.48854

beta1

##     trend 
## 0.0517657

# Calcular la parte común estimada por el modelo
df2_model <- df2_model %>%
  mutate(
    pred_parte_comun = beta0 + beta1 * trend,
    diferencia = market_share - pred_parte_comun
  )

# Extraer los efectos aleatorios estimados por el modelo
efectos_random <- ranef(random)

efectos_random

## Bath and Shower      Deodorants    Depilatories      Fragrances       Hair Care 
##       -3.988697       -3.273814      -10.150282        4.193032        7.908366 
##  Men's Grooming       Skin Care        Sun Care 
##        4.908848       10.099658       -9.697110

# Crear datos para 2026
nuevo_2026 <- data.frame(
  category = names(efectos_random),
  year = 2026,
  trend = 15
)

# Coeficientes
beta0 <- coef(random)["(Intercept)"]
beta1 <- coef(random)["trend"]

# Incorporar los efectos aleatorios
nuevo_2026$efecto_random <- as.numeric(efectos_random)

# Predicción 2026
nuevo_2026$market_share_2026 <-
  beta0 +
  beta1 * nuevo_2026$trend +
  nuevo_2026$efecto_random

nuevo_2026

##          category year trend efecto_random market_share_2026
## 1 Bath and Shower 2026    15     -3.988697          7.276328
## 2      Deodorants 2026    15     -3.273814          7.991210
## 3    Depilatories 2026    15    -10.150282          1.114743
## 4      Fragrances 2026    15      4.193032         15.458056
## 5       Hair Care 2026    15      7.908366         19.173390
## 6  Men's Grooming 2026    15      4.908848         16.173872
## 7       Skin Care 2026    15     10.099658         21.364682
## 8        Sun Care 2026    15     -9.697110          1.567914

actual_2025 <- df2_model %>%
  filter(year == 2025) %>%
  select(
    category,
    market_share_2025 = market_share
  )

resultado_final <- nuevo_2026 %>%
  select(
    category,
    market_share_2026
  ) %>%
  left_join(actual_2025, by = "category") %>%
  mutate(
    cambio_pp = market_share_2026 - market_share_2025
  ) %>%
  arrange(desc(market_share_2026))

resultado_final

##          category market_share_2026 market_share_2025  cambio_pp
## 1       Skin Care         21.364682        22.9054931 -1.5408111
## 2       Hair Care         19.173390        17.9499600  1.2234302
## 3  Men's Grooming         16.173872        15.7874748  0.3863973
## 4      Fragrances         15.458056        18.0941915 -2.6361351
## 5      Deodorants          7.991210         6.9593731  1.0318371
## 6 Bath and Shower          7.276328         6.8055431  0.4707847
## 7        Sun Care          1.567914         1.3796231  0.1882909
## 8    Depilatories          1.114743         0.5972074  0.5175352

#CONCLUSIÓN: Con base en los resultados del modelo de efectos aleatorios, Hair Care representa la alternativa de inversión más atractiva entre las subcategorías analizadas de Beauty & Personal Care para 2026. #Para 2026, Hair Care presenta una participación de mercado estimada de 19.17%, lo que la posiciona como la segunda subcategoría con mayor participación proyectada, únicamente por debajo de Skin Care, con 21.36%. Esto indica que Hair Care ya cuenta con una presencia importante dentro del mercado total de Beauty & Personal Care, por lo que una inversión en esta categoría no dependería del crecimiento de un segmento pequeño o poco consolidado, sino de una categoría que actualmente representa una proporción significativa del mercado. #Además, Hair Care presenta el mayor incremento absoluto esperado en participación de mercado entre las principales subcategorías, pasando de aproximadamente 17.95% en 2025 a 19.17% en 2026. Esto representa un crecimiento de alrededor de 1.22 puntos porcentuales.

Parte 3. Banco mundial

# obtener esperanza de vida de varios países
life_expectancy_data <- wb_data(country = c("MX", "US", "CA"),
                                indicator = "SP.DYN.LE00.IN",
                                start_date = 1950,end_date = 2025)

# Preparar datos
life_data <- life_expectancy_data %>%
  transmute(
    country = country,
    year = as.integer(date),
    life_expectancy = SP.DYN.LE00.IN
  ) %>%
  filter(!is.na(life_expectancy)) %>%
  mutate(
    trend = year - min(year)
  ) %>%
  arrange(country, year)

head(life_data)

## # A tibble: 6 × 4
##   country  year life_expectancy trend
##   <chr>   <int>           <dbl> <int>
## 1 Canada   1960            71.1     0
## 2 Canada   1961            71.3     1
## 3 Canada   1962            71.4     2
## 4 Canada   1963            71.4     3
## 5 Canada   1964            71.8     4
## 6 Canada   1965            71.9     5

life_panel <- pdata.frame(
  life_data,
  index = c("country", "year")
)

pdim(life_panel)

## Balanced Panel: n = 3, T = 65, N = 195

is.pbalanced(life_panel)

## [1] TRUE

plotmeans(
  life_expectancy ~ country,
  data = life_data,
  xlab = "País",
  ylab = "Esperanza de vida"
)

ggplot(
  life_data,
  aes(
    x = year,
    y = life_expectancy,
    color = country
  )
) +
  geom_line(linewidth = 1) +
  labs(
    title = "Evolución de la esperanza de vida",
    x = "Año",
    y = "Esperanza de vida",
    color = "País"
  ) +
  theme_minimal()

#Opción 1 - Efectos agrupados
pooled_life <- plm(
  life_expectancy ~ trend,
  data = life_panel,
  model = "pooling"
)

summary(pooled_life)

## Pooling Model
## 
## Call:
## plm(formula = life_expectancy ~ trend, data = life_panel, model = "pooling")
## 
## Balanced Panel: n = 3, T = 65, N = 195
## 
## Residuals:
##     Min.  1st Qu.   Median  3rd Qu.     Max. 
## -12.5473  -3.4926   1.8974   3.9455   5.0129 
## 
## Coefficients:
##              Estimate Std. Error t-value  Pr(>|t|)    
## (Intercept) 66.120307   0.656051 100.785 < 2.2e-16 ***
## trend        0.223581   0.017686  12.642 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    7575
## Residual Sum of Squares: 4143.7
## R-Squared:      0.45297
## Adj. R-Squared: 0.45013
## F-statistic: 159.814 on 1 and 193 DF, p-value: < 2.22e-16

#Opción 2 - Efectos fijos
fixed_life <- plm(
  life_expectancy ~ trend,
  data = life_panel,
  model = "within",
  effect = "individual"
)

summary(fixed_life)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = life_expectancy ~ trend, data = life_panel, effect = "individual", 
##     model = "within")
## 
## Balanced Panel: n = 3, T = 65, N = 195
## 
## Residuals:
##     Min.  1st Qu.   Median  3rd Qu.     Max. 
## -6.79967 -0.48439  0.29545  1.11777  3.37249 
## 
## Coefficients:
##        Estimate Std. Error t-value  Pr(>|t|)    
## trend 0.2235814  0.0076021   29.41 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    4188.9
## Residual Sum of Squares: 757.67
## R-Squared:      0.81912
## Adj. R-Squared: 0.81628
## F-statistic: 864.969 on 1 and 191 DF, p-value: < 2.22e-16

prueba_f_life <- pFtest(
  fixed_life,
  pooled_life
)

prueba_f_life

## 
##  F test for individual effects
## 
## data:  life_expectancy ~ trend
## F = 426.79, df1 = 2, df2 = 191, p-value < 2.2e-16
## alternative hypothesis: significant effects

#Opción 3 - Efectos aleatorios
random_life <- plm(
  life_expectancy ~ trend,
  data = life_panel,
  model = "random",
  effect = "individual"
)

summary(random_life)

## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = life_expectancy ~ trend, data = life_panel, effect = "individual", 
##     model = "random")
## 
## Balanced Panel: n = 3, T = 65, N = 195
## 
## Effects:
##                  var std.dev share
## idiosyncratic  3.967   1.992 0.132
## individual    25.986   5.098 0.868
## theta: 0.9516
## 
## Residuals:
##     Min.  1st Qu.   Median  3rd Qu.     Max. 
## -7.07789 -0.44781  0.37419  1.19318  3.09427 
## 
## Coefficients:
##               Estimate Std. Error z-value  Pr(>|z|)    
## (Intercept) 66.1203066  2.9565836  22.364 < 2.2e-16 ***
## trend        0.2235814  0.0076021  29.410 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    4196.8
## Residual Sum of Squares: 765.61
## R-Squared:      0.81757
## Adj. R-Squared: 0.81663
## Chisq: 864.969 on 1 DF, p-value: < 2.22e-16

prueba_hausman_life <- phtest(
  fixed_life,
  random_life
)

prueba_hausman_life

## 
##  Hausman Test
## 
## data:  life_expectancy ~ trend
## chisq = 8.1205e-15, df = 1, p-value = 1
## alternative hypothesis: one model is inconsistent

# Extraer efectos aleatorios por país
efectos_life <- ranef(random_life)

efectos_life

##        Canada        Mexico United States 
##      3.991447     -5.734168      1.742721

# Coeficientes del modelo Random
beta0 <- coef(random_life)["(Intercept)"]
beta1 <- coef(random_life)["trend"]

beta0

## (Intercept) 
##    66.12031

beta1

##     trend 
## 0.2235814

# Año inicial de la base
anio_inicial <- min(life_data$year)

# Crear base para predicción 2026
nuevo_2026_life <- data.frame(
  country = names(efectos_life),
  year = 2026,
  trend = 2026 - anio_inicial,
  efecto_random = as.numeric(efectos_life)
)

nuevo_2026_life

##         country year trend efecto_random
## 1        Canada 2026    66      3.991447
## 2        Mexico 2026    66     -5.734168
## 3 United States 2026    66      1.742721

# Predicción de esperanza de vida para 2026
nuevo_2026_life <- nuevo_2026_life %>%
  mutate(
    life_expectancy_2026 =
      beta0 +
      beta1 * trend +
      efecto_random
  )

nuevo_2026_life

##         country year trend efecto_random life_expectancy_2026
## 1        Canada 2026    66      3.991447             84.86813
## 2        Mexico 2026    66     -5.734168             75.14251
## 3 United States 2026    66      1.742721             82.61940

resultado_2026_life <- nuevo_2026_life %>%
  select(
    country,
    life_expectancy_2026
  ) %>%
  arrange(desc(life_expectancy_2026))

resultado_2026_life

##         country life_expectancy_2026
## 1        Canada             84.86813
## 2 United States             82.61940
## 3        Mexico             75.14251

max(life_data$year)

## [1] 2024

ultimo_anio <- max(life_data$year)

actual_life <- life_data %>%
  filter(year == ultimo_anio) %>%
  select(
    country,
    life_expectancy_actual = life_expectancy
  )

resultado_final_life <- resultado_2026_life %>%
  left_join(
    actual_life,
    by = "country"
  ) %>%
  mutate(
    cambio = life_expectancy_2026 -
      life_expectancy_actual
  ) %>%
  arrange(desc(life_expectancy_2026))

resultado_final_life

##         country life_expectancy_2026 life_expectancy_actual     cambio
## 1        Canada             84.86813               82.10805  2.7600790
## 2 United States             82.61940               78.89024  3.7291577
## 3        Mexico             75.14251               75.26400 -0.1214876

#COCLUSIÓN: De acuerdo con el modelo de efectos aleatorios, Canadá presenta la mayor esperanza de vida estimada para 2026, con aproximadamente 84.87 años, seguido por Estados Unidos con 82.62 años y México con 75.14 años. Los resultados muestran que las diferencias estructurales entre los países se mantendrían, con una brecha proyectada de aproximadamente 9.73 años entre Canadá y México y de 7.48 años entre Estados Unidos y México.

file <- "Actividad_1_patentes.Rmd"

x <- readLines(
  file,
  encoding = "UTF-8",
  warn = FALSE
)