Actividad 1

Importar paquetes y llamar librerias

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

Parte 1. Patentes

INSTRUCCIONES: Genera el mejor modelo de la Base de Datos “Patentes” y genera predicciones.

df1 <- read_excel("/Users/adriana/Downloads/PATENT 3.xls")
summary(df1)

##      cusip            merger           employ            return       
##  Min.   :   800   Min.   :0.0000   Min.   :  0.085   Min.   :-73.022  
##  1st Qu.:368514   1st Qu.:0.0000   1st Qu.:  1.227   1st Qu.:  5.128  
##  Median :501116   Median :0.0000   Median :  3.842   Median :  7.585  
##  Mean   :514536   Mean   :0.0177   Mean   : 18.826   Mean   :  8.003  
##  3rd Qu.:754688   3rd Qu.:0.0000   3rd Qu.: 15.442   3rd Qu.: 10.501  
##  Max.   :878555   Max.   :1.0000   Max.   :506.531   Max.   : 48.675  
##                                    NAs    :21        NAs    :8        
##     patents         patentsg           stckpr              rnd           
##  Min.   :  0.0   Min.   :   0.00   Min.   :  0.1875   Min.   :   0.0000  
##  1st Qu.:  1.0   1st Qu.:   1.00   1st Qu.:  7.6250   1st Qu.:   0.6847  
##  Median :  3.0   Median :   4.00   Median : 16.5000   Median :   2.1456  
##  Mean   : 22.9   Mean   :  27.14   Mean   : 22.6270   Mean   :  29.3398  
##  3rd Qu.: 15.0   3rd Qu.:  19.00   3rd Qu.: 29.2500   3rd Qu.:  11.9168  
##  Max.   :906.0   Max.   :1063.00   Max.   :402.0000   Max.   :1719.3535  
##                                    NAs    :2                             
##     rndeflt             rndstck              sales                sic      
##  Min.   :   0.0000   Min.   :   0.1253   Min.   :    1.222   Min.   :2000  
##  1st Qu.:   0.4788   1st Qu.:   5.1520   1st Qu.:   52.995   1st Qu.:2890  
##  Median :   1.4764   Median :  13.3532   Median :  174.065   Median :3531  
##  Mean   :  19.7238   Mean   : 163.8234   Mean   : 1219.601   Mean   :3333  
##  3rd Qu.:   8.7527   3rd Qu.:  74.5625   3rd Qu.:  728.964   3rd Qu.:3661  
##  Max.   :1000.7876   Max.   :9755.3516   Max.   :44224.000   Max.   :9997  
##                      NAs    :157         NAs    :3                         
##       year     
##  Min.   :2012  
##  1st Qu.:2014  
##  Median :2016  
##  Mean   :2016  
##  3rd Qu.:2019  
##  Max.   :2021  
## 

#patents y rnd están muy sesgados: la mitad de las empresas tiene 3 patentes o menos, pero hay una con 906; y hay 538 filas con 0 patentes. Con datos así, un modelo en niveles no es buena idea, por eso usamos logaritmos (log(x+1) para poder incluir los ceros).
df1$log_patents <- log(df1$patents + 1)
df1$log_rnd <- log(df1$rnd + 1)
df1 <- pdata.frame(df1, index = c("cusip", "year"))
modelo_regresion <- lm(log_patents ~ log_rnd, data = df1)
summary(modelo_regresion) #En logaritmos el modelo ajusta mucho mejor (R2 de 0.60 contra 0.31 en niveles) y los residuos son mucho más chicos y estables.

## 
## Call:
## lm(formula = log_patents ~ log_rnd, data = df1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.5640 -0.6233  0.0149  0.6954  2.8844 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.40127    0.03046   13.18   <2e-16 ***
## log_rnd      0.79002    0.01344   58.77   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9754 on 2258 degrees of freedom
## Multiple R-squared:  0.6047, Adjusted R-squared:  0.6046 
## F-statistic:  3454 on 1 and 2258 DF,  p-value: < 2.2e-16

#Prueba de Heterogeneidad
plotmeans(log_patents ~ cusip, data = df1)

#INTERPRETACIÓN: La línea no es horizontal: hay empresas con muchas más patentes que otras, así que sí hay heterogeneidad.

# Opcion 1 - Modelo de Regresión Agrupada (Pooled)
pooled <- plm(log_patents ~ log_rnd, data = df1, model = "pooling")
summary(pooled)

## Pooling Model
## 
## Call:
## plm(formula = log_patents ~ log_rnd, data = df1, model = "pooling")
## 
## Balanced Panel: n = 226, T = 10, N = 2260
## 
## Residuals:
##    Min. 1st Qu.  Median 3rd Qu.    Max. 
## -3.5640 -0.6233  0.0149  0.6954  2.8844 
## 
## Coefficients:
##             Estimate Std. Error t-value  Pr(>|t|)    
## (Intercept) 0.401272   0.030457  13.175 < 2.2e-16 ***
## log_rnd     0.790021   0.013441  58.775 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    5435.2
## Residual Sum of Squares: 2148.4
## R-Squared:      0.60473
## Adj. R-Squared: 0.60455
## F-statistic: 3454.49 on 1 and 2258 DF, p-value: < 2.22e-16

#Opcion 2 - Modelo de Efectos Fijos (within)
within <- plm(log_patents ~ log_rnd, data = df1, model = "within")
summary(within)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = log_patents ~ log_rnd, data = df1, model = "within")
## 
## Balanced Panel: n = 226, T = 10, N = 2260
## 
## Residuals:
##    Min. 1st Qu.  Median 3rd Qu.    Max. 
## -2.7726 -0.2953  0.0371  0.3467  1.7910 
## 
## Coefficients:
##          Estimate Std. Error t-value  Pr(>|t|)    
## log_rnd -0.522079   0.036694 -14.228 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    777.46
## Residual Sum of Squares: 707.06
## R-Squared:      0.090556
## Adj. R-Squared: -0.010543
## F-statistic: 202.432 on 1 and 2033 DF, p-value: < 2.22e-16

#Prueba F
#INTERPRETACIÓN: Si la p es menor a 0.05 no usar POOLED, si p es mayor a 0.05 usar pooled.
pFtest(within, pooled) 

## 
##  F test for individual effects
## 
## data:  log_patents ~ log_rnd
## F = 18.419, df1 = 225, df2 = 2033, p-value < 2.2e-16
## alternative hypothesis: significant effects

#Opcion 3 - Modelo de Efectos Aleatorios (Random)
random <- plm(log_patents ~ log_rnd, data = df1, model = "random")
summary(random)

## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = log_patents ~ log_rnd, data = df1, model = "random")
## 
## Balanced Panel: n = 226, T = 10, N = 2260
## 
## Effects:
##                  var std.dev share
## idiosyncratic 0.3478  0.5897 0.465
## individual    0.3999  0.6324 0.535
## theta: 0.7171
## 
## Residuals:
##    Min. 1st Qu.  Median 3rd Qu.    Max. 
## -3.0612 -0.4923  0.0738  0.5252  1.8901 
## 
## Coefficients:
##             Estimate Std. Error z-value  Pr(>|z|)    
## (Intercept) 1.193600   0.068539  17.415 < 2.2e-16 ***
## log_rnd     0.316866   0.026990  11.740 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    1150.1
## Residual Sum of Squares: 1084
## R-Squared:      0.057528
## Adj. R-Squared: 0.05711
## Chisq: 137.826 on 1 DF, p-value: < 2.22e-16

#Prueba de Hausman
#INTERPRETACIÓN: Si la p es menor a 0.05 usar efectos fijos, si sale mayor usar efectos aleatorios.
phtest(random, within) 

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

#Por lo tanto, el mejor modelo para este panel es el de EFECTOS FIJOS. La prueba F rechaza Pooled (sí hay heterogeneidad) y la prueba de Hausman rechaza Aleatorios (esa heterogeneidad está correlacionada con el gasto en I+D: las empresas que más patentan también son las que más invierten).

#Pronóstico: como el modelo está en logaritmos, hay que "devolver" la predicción a patentes con exp()-1 al final
df1_pronostico <- data.frame(
  cusip = c("149123", "369604", "459200"),
  year = c(2022, 2022, 2022),
  rnd = c(226.5315, 813.7559, 1611.4912) * 1.10
)
df1_pronostico$log_rnd <- log(df1_pronostico$rnd + 1)
intercepto <- fixef(within)
pendiente <- coef(within)["log_rnd"]
df1_pronostico$log_patentes <- intercepto[df1_pronostico$cusip] + pendiente * df1_pronostico$log_rnd
df1_pronostico$patentes <- exp(df1_pronostico$log_patentes) - 1
df1_pronostico[, c("cusip","patentes")]

##    cusip patentes
## 1 149123 141.5081
## 2 369604 498.0288
## 3 459200 281.9298

CONCLUSIÓN: El mejor modelo sigue siendo el de Efectos Fijos, pero usando logaritmos en vez de niveles, porque patents y rnd están muy sesgados (pocas empresas concentran la mayoría de las patentes y del gasto). En logaritmos el modelo ajusta mucho mejor y el coeficiente se interpreta como elasticidad: cuánto cambia el % de patentes por cada % de cambio en I+D. Revisando los datos empresa por empresa encontré algo importante: 2021 (el último año de la base) tiene una caída rarísima en patentes para casi todas las empresas grandes, después de años estables en cientos — es un problema típico de las bases de patentes, donde el último año sale bajo porque el trámite tarda y muchas patentes recientes todavía no se habían registrado cuando se armó la base. Por eso el pronóstico 2022 sale más alto que las patentes “oficiales” de 2021: el modelo se basa en el nivel histórico normal de cada empresa, no en ese último año atípico.

Parte 2. Cuidado de la piel

INSTRUCCIONES: ¿En cuál categoría invertirías y por qué? Si tuvieras que invertir en alguna subcategoría, ¿en cuál lo harías? Justifica de acuerdo al modelo.

df2_excel <- read_excel("/Users/adriana/Downloads/Market sizes.xls")
summary(df2_excel)

##      Geography       Category      Data Type         Unit     Current Constant
##  Length   :13   Length   :13   Length   :13   Length   :13   Length   :13     
##  N.unique : 1   N.unique :13   N.unique : 1   N.unique : 1   N.unique : 1     
##  N.blank  : 0   N.blank  : 0   N.blank  : 0   N.blank  : 0   N.blank  : 0     
##  Min.nchar: 6   Min.nchar: 8   Min.nchar:16   Min.nchar:11   Min.nchar:14     
##  Max.nchar: 6   Max.nchar:39   Max.nchar:16   Max.nchar:11   Max.nchar:14     
##                                                                               
##         2011           2012           2013           2014           2015   
##  Length   :13   Length   :13   Length   :13   Length   :13   Length   :13  
##  N.unique :12   N.unique :12   N.unique :12   N.unique :12   N.unique :12  
##  N.blank  : 0   N.blank  : 0   N.blank  : 0   N.blank  : 0   N.blank  : 0  
##  Min.nchar: 1   Min.nchar: 1   Min.nchar: 1   Min.nchar: 1   Min.nchar: 1  
##  Max.nchar:18   Max.nchar:18   Max.nchar:18   Max.nchar:18   Max.nchar:18  
##                                                                            
##       2016             2017             2018             2019       
##  Min.   :  1188   Min.   :  1320   Min.   :  1439   Min.   :  1516  
##  1st Qu.: 11482   1st Qu.: 12123   1st Qu.: 12490   1st Qu.: 12960  
##  Median : 21783   Median : 24478   Median : 26344   Median : 27806  
##  Mean   : 37710   Mean   : 40692   Mean   : 42891   Mean   : 44853  
##  3rd Qu.: 32384   3rd Qu.: 34260   3rd Qu.: 35597   3rd Qu.: 37247  
##  Max.   :170488   Max.   :183580   Max.   :192727   Max.   :200916  
##       2020             2021             2022             2023       
##  Min.   :  1491   Min.   :  1564   Min.   :  1689   Min.   :  1850  
##  1st Qu.: 13835   1st Qu.: 15091   1st Qu.: 16815   1st Qu.: 19001  
##  Median : 23213   Median : 30225   Median : 37173   Median : 43366  
##  Mean   : 42589   Mean   : 47822   Mean   : 53823   Mean   : 61738  
##  3rd Qu.: 36951   3rd Qu.: 39555   3rd Qu.: 42739   3rd Qu.: 46900  
##  Max.   :191490   Max.   :211921   Max.   :236406   Max.   :268852  
##       2024             2025       
##  Min.   :  1802   Min.   :  1873  
##  1st Qu.: 20235   1st Qu.: 21342  
##  Median : 47334   Median : 49508  
##  Mean   : 68537   Mean   : 72844  
##  3rd Qu.: 53798   3rd Qu.: 58846  
##  Max.   :296120   Max.   :313593

df2 <- data.frame(
  categoria = c(rep("Bath and Shower",15), rep("Deodorants",15), rep("Depilatories",15),
                rep("Fragrances",15), rep("Hair Care",15), rep("Mens Grooming",15),
                rep("Skin Care",15), rep("Sun Care",15)),
  anio = rep(2011:2025, 8),
  valor = c(
    8410.4, 9085, 9711.1, 10226.9, 10814.8, 11481.6, 12123.1, 12490.1, 12959.5, 14633.9, 15469.3, 16815.3, 19001, 20234.6, 21341.7,
    9151, 10148.3, 10697.7, 11469.7, 12379.2, 13723.4, 14891.6, 15913.6, 16745.1, 13834.8, 15091.1, 17872.4, 19454.9, 20787.8, 21824.1,
    802.5, 897.7, 969.9, 1036.3, 1106.4, 1187.9, 1319.5, 1439.1, 1516.1, 1491.2, 1563.5, 1688.6, 1850.2, 1802, 1872.8,
    18727.9, 19957.2, 20651.8, 20979.6, 22567.6, 24208.9, 26152.3, 27536, 28483.1, 25516.1, 31800.8, 37172.8, 45318.9, 51579.2, 56742.1,
    24853.9, 26565.6, 27672.5, 29154.6, 30416.2, 32383.6, 34260, 35597.2, 37247.4, 36950.7, 39554.7, 42738.8, 46900.3, 53024.9, 56289.8,
    18672.1, 20694.4, 22133.4, 23175.3, 24774, 27446, 29497.2, 31498.1, 32693.3, 29669.4, 32938.6, 37702.4, 43366.2, 47333.8, 49508.4,
    25650.3, 27306.2, 28366.9, 29438.3, 30964, 33054.3, 35787.8, 38067.5, 40637.5, 42196.5, 47980.6, 52939.6, 61023.6, 68590.4, 71830,
    1144.7, 1303.1, 1417.3, 1552.8, 1693.3, 1944.1, 2143.3, 2266.5, 2443.8, 2042.9, 2293.6, 2862, 3736.4, 4183.4, 4326.4
  )
)
df2$t <- df2$anio #copia numérica de anio; anio ya es el índice del panel, así que si lo usamos directo en el modelo R lo vuelve categórico (uno por año) en vez de número, y por eso salían los NA
df2 <- pdata.frame(df2, index = c("categoria","anio"))

modelo_regresion2 <- lm(valor ~ t, data = df2)
summary(modelo_regresion2)

## 
## Call:
## lm(formula = valor ~ t, data = df2)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -30062 -10733   -673  13521  39895 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2937049.2   670768.7  -4.379 2.60e-05 ***
## t               1466.2      332.4   4.411 2.29e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15730 on 118 degrees of freedom
## Multiple R-squared:  0.1415, Adjusted R-squared:  0.1343 
## F-statistic: 19.46 on 1 and 118 DF,  p-value: 2.287e-05

#Prueba de Heterogeneidad
plotmeans(valor ~ categoria, data = df2)

#INTERPRETACIÓN: Skin Care y Hair Care operan en un nivel de mercado mucho más grande que Depilatories o Sun Care: sí hay heterogeneidad.

# Opcion 1 - Modelo de Regresión Agrupada (Pooled)
pooled2 <- plm(valor ~ t, data = df2, model = "pooling")
summary(pooled2)

## Pooling Model
## 
## Call:
## plm(formula = valor ~ t, data = df2, model = "pooling")
## 
## Balanced Panel: n = 8, T = 15, N = 120
## 
## Residuals:
##    Min. 1st Qu.  Median 3rd Qu.    Max. 
##  -30062  -10733    -673   13521   39895 
## 
## Coefficients:
##                Estimate  Std. Error t-value  Pr(>|t|)    
## (Intercept) -2937049.22   670768.71 -4.3786 2.598e-05 ***
## t               1466.17      332.39  4.4110 2.287e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    3.4018e+10
## Residual Sum of Squares: 2.9203e+10
## R-Squared:      0.14155
## Adj. R-Squared: 0.13427
## F-statistic: 19.4565 on 1 and 118 DF, p-value: 2.2866e-05

#Opcion 2 - Modelo de Efectos Fijos (within)
within2 <- plm(valor ~ t, data = df2, model = "within")
summary(within2)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = valor ~ t, data = df2, model = "within")
## 
## Balanced Panel: n = 8, T = 15, N = 120
## 
## Residuals:
##    Min. 1st Qu.  Median 3rd Qu.    Max. 
##   -9760   -3015   -1523    2555   19311 
## 
## Coefficients:
##   Estimate Std. Error t-value  Pr(>|t|)    
## t  1466.17     116.12  12.626 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    8.168e+09
## Residual Sum of Squares: 3352800000
## R-Squared:      0.58952
## Adj. R-Squared: 0.55993
## F-statistic: 159.413 on 1 and 111 DF, p-value: < 2.22e-16

#Prueba F
pFtest(within2, pooled2) 

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

#Opcion 3 - Modelo de Efectos Aleatorios (Random)
random2 <- plm(valor ~ t, data = df2, model = "random")
summary(random2)

## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = valor ~ t, data = df2, model = "random")
## 
## Balanced Panel: n = 8, T = 15, N = 120
## 
## Effects:
##                     var   std.dev share
## idiosyncratic  30205839      5496  0.11
## individual    244180653     15626  0.89
## theta: 0.9096
## 
## Residuals:
##    Min. 1st Qu.  Median 3rd Qu.    Max. 
##  -11596   -3144    -709    1841   21173 
## 
## Coefficients:
##                Estimate  Std. Error z-value  Pr(>|z|)    
## (Intercept) -2937049.22   234403.59 -12.530 < 2.2e-16 ***
## t               1466.17      116.12  12.626 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    8379500000
## Residual Sum of Squares: 3564300000
## R-Squared:      0.57464
## Adj. R-Squared: 0.57104
## Chisq: 159.413 on 1 DF, p-value: < 2.22e-16

#Prueba de Hausman
#INTERPRETACIÓN: Si la p es menor a 0.05 usar efectos fijos, si sale mayor usar efectos aleatorios.
phtest(random2, within2) 

## 
##  Hausman Test
## 
## data:  valor ~ t
## chisq = 1.9543e-12, df = 1, p-value = 1
## alternative hypothesis: one model is inconsistent

#Por lo tanto, el mejor modelo para este panel es el de EFECTOS ALEATORIOS.

#Pronóstico 2026 para cada categoría de producto
df2_pronostico <- data.frame(
  categoria = c("Bath and Shower","Deodorants","Depilatories","Fragrances","Hair Care","Mens Grooming","Skin Care","Sun Care"),
  t = rep(2026, 8)
)
intercepto2 <- coef(random2)["(Intercept)"]
pendiente2 <- coef(random2)["t"]
efectos2 <- ranef(random2)
df2_pronostico$valor <- intercepto2 + efectos2[df2_pronostico$categoria] + pendiente2*df2_pronostico$t
df2_pronostico[order(-df2_pronostico$valor), c("categoria","valor")]

##         categoria    valor
## 7       Skin Care 53816.52
## 5       Hair Care 48512.05
## 6   Mens Grooming 43056.53
## 4      Fragrances 42150.12
## 2      Deodorants 26716.76
## 1 Bath and Shower 25448.13
## 8        Sun Care 14244.21
## 3    Depilatories 13264.96

df3 <- data.frame(
  categoria = c(rep("Premium",10), rep("Prestige",10), rep("Mass",10), rep("Dermocosmetics",10)),
  anio = rep(2016:2025, 4),
  valor = c(
    21782.8, 24478.4, 26344, 27806.3, 23212.9, 30224.8, 37374.7, 46017.8, 53798.1, 58845.9,
    19715, 21722.1, 23189.7, 24244.8, 19409.1, 25615.4, 31797.5, 38290.7, 44287.3, 48538.1,
    129002.7, 138197.6, 145032.6, 151178.5, 146414.9, 158941.4, 174599.9, 195814, 213682.8, 225083.1,
    3808.3, 4841.9, 5483.3, 6216.3, 6792.4, 8288.2, 9733.7, 12967.4, 15551.4, 17172.2
  )
)
df3$t <- df3$anio #misma razón: copia numérica para el modelo, anio se queda como índice del panel
df3 <- pdata.frame(df3, index = c("categoria","anio"))

#Prueba de Heterogeneidad
plotmeans(valor ~ categoria, data = df3)

#INTERPRETACIÓN: Mass es mucho más grandes que Dermocosmetics: sí hay heterogeneidad.

# Opcion 1 - Modelo de Regresión Agrupada (Pooled)
pooled3 <- plm(valor ~ t, data = df3, model = "pooling")
summary(pooled3)

## Pooling Model
## 
## Call:
## plm(formula = valor ~ t, data = df3, model = "pooling")
## 
## Balanced Panel: n = 4, T = 10, N = 40
## 
## Residuals:
##    Min. 1st Qu.  Median 3rd Qu.    Max. 
##  -64823  -37524  -27657    9359  143088 
## 
## Coefficients:
##               Estimate Std. Error t-value Pr(>|t|)
## (Intercept) -9641543.5  7277989.5 -1.3248   0.1932
## t               4801.7     3602.1  1.3331   0.1905
## 
## Total Sum of Squares:    1.7031e+11
## Residual Sum of Squares: 1.6271e+11
## R-Squared:      0.044675
## Adj. R-Squared: 0.019535
## F-statistic: 1.77703 on 1 and 38 DF, p-value: 0.19045

#Opcion 2 - Modelo de Efectos Fijos (within)
within3 <- plm(valor ~ t, data = df3, model = "within")
summary(within3)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = valor ~ t, data = df3, model = "within")
## 
## Balanced Panel: n = 4, T = 10, N = 40
## 
## Residuals:
##    Min. 1st Qu.  Median 3rd Qu.    Max. 
##  -18979   -7934   -1587    5709   35680 
## 
## Coefficients:
##   Estimate Std. Error t-value Pr(>|t|)    
## t  4801.75     667.31  7.1957 2.14e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    1.2752e+10
## Residual Sum of Squares: 5143300000
## R-Squared:      0.59667
## Adj. R-Squared: 0.55058
## F-statistic: 51.7777 on 1 and 35 DF, p-value: 2.1398e-08

#Prueba F
pFtest(within3, pooled3)

## 
##  F test for individual effects
## 
## data:  valor ~ t
## F = 357.4, df1 = 3, df2 = 35, p-value < 2.2e-16
## alternative hypothesis: significant effects

#Opcion 3 - Modelo de Efectos Aleatorios (Random)
random3 <- plm(valor ~ t, data = df3, model = "random")
summary(random3)

## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = valor ~ t, data = df3, model = "random")
## 
## Balanced Panel: n = 4, T = 10, N = 40
## 
## Effects:
##                     var   std.dev share
## idiosyncratic 1.470e+08 1.212e+04 0.027
## individual    5.237e+09 7.237e+04 0.973
## theta: 0.9471
## 
## Residuals:
##    Min. 1st Qu.  Median 3rd Qu.    Max. 
##  -16235   -7355   -3169    5035   41362 
## 
## Coefficients:
##                Estimate  Std. Error z-value  Pr(>|z|)    
## (Intercept) -9641543.53  1348787.76 -7.1483 8.786e-13 ***
## t               4801.75      667.31  7.1957 6.215e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    1.3193e+10
## Residual Sum of Squares: 5584100000
## R-Squared:      0.57673
## Adj. R-Squared: 0.56559
## Chisq: 51.7777 on 1 DF, p-value: 6.2154e-13

#Prueba de Hausman
phtest(random3, within3)

## 
##  Hausman Test
## 
## data:  valor ~ t
## chisq = 5.1971e-13, df = 1, p-value = 1
## alternative hypothesis: one model is inconsistent

#Por lo tanto, el mejor modelo para este panel también es el de EFECTOS ALEATORIOS, por la misma razón que en las categorías de producto.

#Pronóstico 2026 para cada subcategoría
categorias_tier <- c("Premium", "Prestige", "Mass", "Dermocosmetics")
df3_pronostico <- data.frame(categoria = categorias_tier, t = rep(2026, 4))
intercepto3 <- coef(random3)["(Intercept)"]
pendiente3 <- coef(random3)["t"]
efectos3 <- ranef(random3)
df3_pronostico$valor <- intercepto3 + efectos3[df3_pronostico$categoria] + pendiente3*df3_pronostico$t
df3_pronostico[order(-df3_pronostico$valor), c("categoria","valor")]

##        categoria     valor
## 3           Mass 193903.84
## 1        Premium  61469.25
## 2       Prestige  56176.50
## 4 Dermocosmetics  35638.66

CONCLUSIÓN: En los dos paneles el modelo que mejor funcionó fue el de Efectos Aleatorios, porque sí hay diferencias de tamaño entre categorías pero no interfieren con la tendencia en el tiempo. Invertiría en Skin Care porque ya es la categoría más grande del mercado y encima sigue proyectándose como la de mayor valor para 2026. Y a una subcategoría, invertiria en Dermocosmetics: todavía es chica, pero es la que más rápido está creciendo (casi el doble que las demás), así que es donde veo más potencial a futuro.

Parte 3. Banco mundial

PIB (crecimiento anual, %) y Desempleo (% de la fuerza laboral) — 5 países, 2000-2023

INSTRUCCIONES: Importa datos del Banco Mundial y genera el mejor modelo, incluyendo predicciones, conclusiones y gráficas.

wb <- tryCatch({
  library(WDI)
  d <- WDI(country = c("MX","BR","CA","US","AR"),
           indicator = c(pib_crecimiento = "NY.GDP.MKTP.KD.ZG",
                         desempleo = "SL.UEM.TOTL.ZS"),
           start = 2000, end = 2023)
  d <- na.omit(d[, c("country","year","pib_crecimiento","desempleo")])
  names(d) <- c("pais","anio","pib_crecimiento","desempleo")
  d
}, error = function(e) read.csv("panel_wb.csv"))

pdfWB <- pdata.frame(wb, index = c("pais", "anio"))

modelo_regresion4 <- lm(pib_crecimiento ~ desempleo, data = pdfWB)
summary(modelo_regresion4)

## 
## Call:
## lm(formula = pib_crecimiento ~ desempleo, data = pdfWB)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.1604  -1.4021   0.2238   1.3557   8.6690 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.4172     0.8283   4.126 6.92e-05 ***
## desempleo    -0.1882     0.1039  -1.811   0.0726 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.586 on 118 degrees of freedom
## Multiple R-squared:  0.02706,    Adjusted R-squared:  0.01881 
## F-statistic: 3.282 on 1 and 118 DF,  p-value: 0.07261

#Prueba de Heterogeneidad
plotmeans(pib_crecimiento ~ pais, data = wb, n.label = FALSE)

#INTERPRETACIÓN: Aquí la línea sale más plana el crecimiento promedio del PIB no varía tanto entre estos 5 países, así que podría no haber heterogeneidad fuerte.

# Opcion 1 - Modelo de Regresión Agrupada (Pooled)
pooled4 <- plm(pib_crecimiento ~ desempleo, data = pdfWB, model = "pooling")
summary(pooled4)

## Pooling Model
## 
## Call:
## plm(formula = pib_crecimiento ~ desempleo, data = pdfWB, model = "pooling")
## 
## Balanced Panel: n = 5, T = 24, N = 120
## 
## Residuals:
##    Min. 1st Qu.  Median 3rd Qu.    Max. 
## -11.160  -1.402   0.224   1.356   8.669 
## 
## Coefficients:
##             Estimate Std. Error t-value  Pr(>|t|)    
## (Intercept)  3.41718    0.82826  4.1257 6.915e-05 ***
## desempleo   -0.18823    0.10391 -1.8115   0.07261 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    1559.7
## Residual Sum of Squares: 1517.5
## R-Squared:      0.027057
## Adj. R-Squared: 0.018812
## F-statistic: 3.28151 on 1 and 118 DF, p-value: 0.072608

#Opcion 2 - Modelo de Efectos Fijos (within)
within4 <- plm(pib_crecimiento ~ desempleo, data = pdfWB, model = "within")
summary(within4)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = pib_crecimiento ~ desempleo, data = pdfWB, model = "within")
## 
## Balanced Panel: n = 5, T = 24, N = 120
## 
## Residuals:
##    Min. 1st Qu.  Median 3rd Qu.    Max. 
## -11.025  -0.982   0.106   1.504   9.647 
## 
## Coefficients:
##           Estimate Std. Error t-value Pr(>|t|)   
## desempleo -0.49609    0.15661 -3.1677 0.001972 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    1553
## Residual Sum of Squares: 1427.3
## R-Squared:      0.080901
## Adj. R-Squared: 0.040589
## F-statistic: 10.0345 on 1 and 114 DF, p-value: 0.001972

#Prueba F
#INTERPRETACIÓN: Si la p es menor a 0.05 no usar POOLED, si p es mayor a 0.05 usar pooled.
pFtest(within4, pooled4)

## 
##  F test for individual effects
## 
## data:  pib_crecimiento ~ desempleo
## F = 1.8007, df1 = 4, df2 = 114, p-value = 0.1335
## alternative hypothesis: significant effects

#Opcion 3 - Modelo de Efectos Aleatorios (Random)
random4 <- plm(pib_crecimiento ~ desempleo, data = pdfWB, model = "random")
summary(random4)

## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = pib_crecimiento ~ desempleo, data = pdfWB, model = "random")
## 
## Balanced Panel: n = 5, T = 24, N = 120
## 
## Effects:
##                  var std.dev share
## idiosyncratic 12.520   3.538     1
## individual     0.000   0.000     0
## theta: 0
## 
## Residuals:
##    Min. 1st Qu.  Median 3rd Qu.    Max. 
## -11.160  -1.402   0.224   1.356   8.669 
## 
## Coefficients:
##             Estimate Std. Error z-value  Pr(>|z|)    
## (Intercept)  3.41718    0.82826  4.1257 3.696e-05 ***
## desempleo   -0.18823    0.10391 -1.8115   0.07006 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    1559.7
## Residual Sum of Squares: 1517.5
## R-Squared:      0.027057
## Adj. R-Squared: 0.018812
## Chisq: 3.28151 on 1 DF, p-value: 0.070064

#Prueba de Hausman
phtest(random4, within4)

## 
##  Hausman Test
## 
## data:  pib_crecimiento ~ desempleo
## chisq = 6.9035, df = 1, p-value = 0.008603
## alternative hypothesis: one model is inconsistent

#Por lo tanto, el mejor modelo para este panel es el de POOLED. La prueba F da p > 0.05: con solo 5 países, no hay evidencia suficiente de heterogeneidad como para justificar Efectos Fijos o Aleatorios.

#Pronóstico 2024 para cada país, a partir de su desempleo observado en 2023
intercepto4 <- coef(pooled4)["(Intercept)"]
pendiente4 <- coef(pooled4)["desempleo"]
df4_pronostico <- wb[wb$anio == max(wb$anio), c("pais", "desempleo", "pib_crecimiento")]
names(df4_pronostico)[3] <- "pib_observado"
df4_pronostico$pib_2024_pred <- intercepto4 + pendiente4 * df4_pronostico$desempleo
df4_pronostico[order(-df4_pronostico$pib_2024_pred), ]

##              pais desempleo pib_observado pib_2024_pred
## 96         Mexico     2.765      3.106796      2.896729
## 120 United States     3.638      2.934363      2.732406
## 72         Canada     5.415      1.953086      2.397927
## 24      Argentina     6.139     -1.855788      2.261650
## 48         Brazil     7.947      3.241655      1.921335

#Gráfica: relación PIB vs. desempleo
plot(wb$desempleo, wb$pib_crecimiento, col = as.factor(wb$pais), pch = 19,
     xlab = "Desempleo (%)", ylab = "PIB, crecimiento anual (%)",
     main = "PIB vs. Desempleo, 2000-2023 (5 países)")
abline(a = intercepto4, b = pendiente4, col = "black", lwd = 2)
legend("topright", legend = levels(as.factor(wb$pais)), col = 1:5, pch = 19, cex = 0.8)

CONCLUSIÓN: El mejor modelo para el panel de PIB y desempleo es el Pooled, porque con solo 5 países la prueba F no encuentra evidencia suficiente de heterogeneidad entre ellos. Existe una relación negativa entre desempleo y crecimiento del PIB: por cada punto porcentual adicional de desempleo, el PIB crece en promedio 0.19 puntos porcentuales menos.