Actividad 1

Actividad 1. Modelo Econométrico

Código Sesisón 2, 3, y 4. Análisis de Datos de Panel

library(WDI)
library(wbstats)
library(tidyverse)
library(ggplot2)
library(gplots)
library(plm)
library(lmtest)
library(pglm)

#Información de México
PIB_MEX <- wb_data(country = "MX", indicator = "NY.GDP.PCAP.CD", start_date = 1990, end_date = 2025)

ggplot(PIB_MEX, aes(x=date, y=NY.GDP.PCAP.CD)) + 
  geom_point() + 
  geom_line() + 
  labs(title = "PIB per Capita en México (Current USD$)", x = "Año", y = "Valor")

#Información de varios países
PIB_PANEL <- wb_data(country = c("MX", "US", "DK"), indicator = "NY.GDP.PCAP.CD", start_date = 1990, end_date = 2025)

ggplot(PIB_PANEL, aes(x=date, y=NY.GDP.PCAP.CD, color=iso3c)) + 
  geom_point() + 
  geom_line() + 
  labs(title = "PIB per Capita Panel (Current USD$)", x = "Año", y = "Valor")

#Información de varios países e indicadores
VI_PANEL <- wb_data(country = c("MX", "US", "DK"), indicator = c("NY.GDP.PCAP.CD", "SP.DYN.LE00.IN"), start_date = 1990, end_date = 2025)

ggplot(VI_PANEL, aes(x=date, y=NY.GDP.PCAP.CD, color=iso3c)) + 
  geom_point() + 
  geom_line() + 
  labs(title = "PIB per Capita Panel (Current USD$)", x = "Año", y = "Valor")

#Heterogeneidad: Variación entre individuos
plotmeans(NY.GDP.PCAP.CD ~ country, main = "Heterogeneidad entre países", xlab = "País", ylab = "PIB per Capita", data = VI_PANEL)

#Alta Heterogeneida: Si los puntos (medias) están muy separados entre sí. 


#MODELOS DE EFECTO FIJO Y ALEATORIOS

#Paso 1. Convertir la base de datos a formato de panel
datos_panel <- pdata.frame(PIB_PANEL, index = c("country", "date"))

#Modelo de Efectos Fijos
modelo_efecto_fijo <- plm(NY.GDP.PCAP.CD ~ date, data = datos_panel, model = "within")
summary(modelo_efecto_fijo)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = NY.GDP.PCAP.CD ~ date, data = datos_panel, model = "within")
## 
## Balanced Panel: n = 3, T = 34, N = 102
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -15362.18  -4826.15   -700.84   5571.86  15709.61 
## 
## Coefficients:
##          Estimate Std. Error t-value  Pr(>|t|)    
## date1991   375.05    6499.90  0.0577 0.9541617    
## date1992  1763.43    6499.90  0.2713 0.7870066    
## date1993  2029.53    6499.90  0.3122 0.7558421    
## date1994  3312.09    6499.90  0.5096 0.6120589    
## date1995  4755.69    6499.90  0.7317 0.4669682    
## date1996  5445.10    6499.90  0.8377 0.4052111    
## date1997  5292.53    6499.90  0.8142 0.4184310    
## date1998  6034.70    6499.90  0.9284 0.3565682    
## date1999  6837.61    6499.90  1.0520 0.2966562    
## date2000  6882.42    6499.90  1.0588 0.2935281    
## date2001  7308.70    6499.90  1.1244 0.2649026    
## date2002  8436.39    6499.90  1.2979 0.1988305    
## date2003 11171.57    6499.90  1.7187 0.0903514 .  
## date2004 14090.20    6499.90  2.1678 0.0337889 *  
## date2005 15930.73    6499.90  2.4509 0.0169060 *  
## date2006 18001.68    6499.90  2.7695 0.0072818 ** 
## date2007 20962.48    6499.90  3.2250 0.0019602 ** 
## date2008 23260.58    6499.90  3.5786 0.0006534 ***
## date2009 20034.17    6499.90  3.0822 0.0029968 ** 
## date2010 20851.75    6499.90  3.2080 0.0020632 ** 
## date2011 22870.48    6499.90  3.5186 0.0007910 ***
## date2012 22347.85    6499.90  3.4382 0.0010188 ** 
## date2013 24025.13    6499.90  3.6962 0.0004469 ***
## date2014 25095.79    6499.90  3.8610 0.0002597 ***
## date2015 22075.56    6499.90  3.3963 0.0011609 ** 
## date2016 22625.57    6499.90  3.4809 0.0008909 ***
## date2017 24521.44    6499.90  3.7726 0.0003480 ***
## date2018 26893.95    6499.90  4.1376 0.0001015 ***
## date2019 27149.91    6499.90  4.1770 8.855e-05 ***
## date2020 26769.79    6499.90  4.1185 0.0001084 ***
## date2021 32477.19    6499.90  4.9966 4.547e-06 ***
## date2022 34527.71    6499.90  5.3120 1.371e-06 ***
## date2023 37028.18    6499.90  5.6967 3.075e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    1.5466e+10
## Residual Sum of Squares: 4182600000
## R-Squared:      0.72956
## Adj. R-Squared: 0.58615
## F-statistic: 5.39549 on 33 and 66 DF, p-value: 3.151e-09

#Modelo de Efectos Aleatorios
modelo_efecto_aleatorio <- plm(NY.GDP.PCAP.CD ~ date, data = datos_panel, model = "random")
summary(modelo_efecto_aleatorio)

## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = NY.GDP.PCAP.CD ~ date, data = datos_panel, model = "random")
## 
## Balanced Panel: n = 3, T = 34, N = 102
## 
## Effects:
##                     var   std.dev share
## idiosyncratic  63373056      7961 0.112
## individual    500147856     22364 0.888
## theta: 0.9391
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -16937.44  -4810.28   -561.29   5934.15  16444.19 
## 
## Coefficients:
##             Estimate Std. Error z-value  Pr(>|z|)    
## (Intercept) 17976.26   13705.48  1.3116 0.1896516    
## date1991      375.05    6499.90  0.0577 0.9539874    
## date1992     1763.43    6499.90  0.2713 0.7861599    
## date1993     2029.53    6499.90  0.3122 0.7548577    
## date1994     3312.09    6499.90  0.5096 0.6103593    
## date1995     4755.69    6499.90  0.7317 0.4643782    
## date1996     5445.10    6499.90  0.8377 0.4021877    
## date1997     5292.53    6499.90  0.8142 0.4155028    
## date1998     6034.70    6499.90  0.9284 0.3531851    
## date1999     6837.61    6499.90  1.0520 0.2928199    
## date2000     6882.42    6499.90  1.0588 0.2896682    
## date2001     7308.70    6499.90  1.1244 0.2608297    
## date2002     8436.39    6499.90  1.2979 0.1943127    
## date2003    11171.57    6499.90  1.7187 0.0856637 .  
## date2004    14090.20    6499.90  2.1678 0.0301773 *  
## date2005    15930.73    6499.90  2.4509 0.0142492 *  
## date2006    18001.68    6499.90  2.7695 0.0056137 ** 
## date2007    20962.48    6499.90  3.2250 0.0012595 ** 
## date2008    23260.58    6499.90  3.5786 0.0003454 ***
## date2009    20034.17    6499.90  3.0822 0.0020546 ** 
## date2010    20851.75    6499.90  3.2080 0.0013366 ** 
## date2011    22870.48    6499.90  3.5186 0.0004338 ***
## date2012    22347.85    6499.90  3.4382 0.0005856 ***
## date2013    24025.13    6499.90  3.6962 0.0002188 ***
## date2014    25095.79    6499.90  3.8610 0.0001129 ***
## date2015    22075.56    6499.90  3.3963 0.0006831 ***
## date2016    22625.57    6499.90  3.4809 0.0004997 ***
## date2017    24521.44    6499.90  3.7726 0.0001616 ***
## date2018    26893.95    6499.90  4.1376 3.510e-05 ***
## date2019    27149.91    6499.90  4.1770 2.954e-05 ***
## date2020    26769.79    6499.90  4.1185 3.814e-05 ***
## date2021    32477.19    6499.90  4.9966 5.836e-07 ***
## date2022    34527.71    6499.90  5.3120 1.084e-07 ***
## date2023    37028.18    6499.90  5.6967 1.221e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    1.5593e+10
## Residual Sum of Squares: 4309400000
## R-Squared:      0.72363
## Adj. R-Squared: 0.58952
## Chisq: 178.051 on 33 DF, p-value: < 2.22e-16

#Prueba de Hausman
phtest(modelo_efecto_fijo, modelo_efecto_aleatorio)

## 
##  Hausman Test
## 
## data:  NY.GDP.PCAP.CD ~ date
## chisq = 1.3283e-13, df = 33, p-value = 1
## alternative hypothesis: one model is inconsistent

  #Como el p-value es mayor a 0.05, usamos el modelo de Efectos Aleatorios

Aplicación de Shiny (Ejemplo y ejercicio sesión 3 por mesas)

App de Shiny

Actividad de Patentes

Importar la base de datos

library(readxl)
patentes <- read_xls("/Users/lishdz/Desktop/8vo/periodo 1/modulo 1 - r/PATENT 3.xls")

Entender la base de datos

summary(patentes)

##      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  
##                                    NA's   :21        NA's   :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  
##                                    NA's   :2                             
##     rndeflt             rndstck             sales               sic      
##  Min.   :   0.0000   Min.   :   0.125   Min.   :    1.22   Min.   :2000  
##  1st Qu.:   0.4788   1st Qu.:   5.152   1st Qu.:   52.99   1st Qu.:2890  
##  Median :   1.4764   Median :  13.353   Median :  174.06   Median :3531  
##  Mean   :  19.7238   Mean   : 163.823   Mean   : 1219.60   Mean   :3333  
##  3rd Qu.:   8.7527   3rd Qu.:  74.563   3rd Qu.:  728.96   3rd Qu.:3661  
##  Max.   :1000.7876   Max.   :9755.352   Max.   :44224.00   Max.   :9997  
##                      NA's   :157        NA's   :3                        
##       year     
##  Min.   :2012  
##  1st Qu.:2014  
##  Median :2016  
##  Mean   :2016  
##  3rd Qu.:2019  
##  Max.   :2021  
## 

sum(is.na(patentes))

## [1] 191

sapply(patentes, function(x) sum(is.na(x)))

##    cusip   merger   employ   return  patents patentsg   stckpr      rnd 
##        0        0       21        8        0        0        2        0 
##  rndeflt  rndstck    sales      sic     year 
##        0      157        3        0        0

patentes1 <- na.omit(patentes)

Construir base de datos tipo panel

patentes_panel <- pdata.frame(patentes1, index = c("cusip", "year"))

Modelo de Efectos Fijos y Aleatorios

#Modelo de Efectos Fijos
modelo_efecto_fijo_patentes <- plm(patents ~ merger + employ + return + patentsg + stckpr + rnd + rndeflt + rndstck + sales + sic, data = patentes_panel, model = "within")
summary(modelo_efecto_fijo_patentes)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = patents ~ merger + employ + return + patentsg + 
##     stckpr + rnd + rndeflt + rndstck + sales + sic, data = patentes_panel, 
##     model = "within")
## 
## Unbalanced Panel: n = 215, T = 2-10, N = 2083
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -468.39577   -1.75634   -0.25666    1.85265  172.64513 
## 
## Coefficients:
##             Estimate  Std. Error  t-value  Pr(>|t|)    
## merger    6.02467998  4.30535335   1.3993    0.1619    
## employ   -0.09095534  0.08057733  -1.1288    0.2591    
## return   -0.01221444  0.12005904  -0.1017    0.9190    
## patentsg  0.03913907  0.02580379   1.5168    0.1295    
## stckpr   -0.03959771  0.03347713  -1.1828    0.2370    
## rnd      -2.04101003  0.15053766 -13.5581 < 2.2e-16 ***
## rndeflt   3.25369409  0.22523191  14.4460 < 2.2e-16 ***
## rndstck   0.19724166  0.01808942  10.9037 < 2.2e-16 ***
## sales    -0.00188938  0.00041715  -4.5293 6.294e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    1090400
## Residual Sum of Squares: 714450
## R-Squared:      0.34479
## Adj. R-Squared: 0.2662
## F-statistic: 108.696 on 9 and 1859 DF, p-value: < 2.22e-16

#Modelo de Efectos Aleatorios
modelo_efecto_aleatorio_patentes <- plm(patents ~ merger + employ + return + patentsg + stckpr + rnd + rndeflt + rndstck + sales + sic, data = patentes_panel, model = "random")
summary(modelo_efecto_aleatorio_patentes)

## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = patents ~ merger + employ + return + patentsg + 
##     stckpr + rnd + rndeflt + rndstck + sales + sic, data = patentes_panel, 
##     model = "random")
## 
## Unbalanced Panel: n = 215, T = 2-10, N = 2083
## 
## Effects:
##                 var std.dev share
## idiosyncratic 384.3    19.6     1
## individual      0.0     0.0     0
## theta:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       0       0       0       0       0 
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -525.42194   -2.59738   -0.31264    1.88763  277.92369 
## 
## Coefficients:
##                Estimate  Std. Error z-value  Pr(>|z|)    
## (Intercept)  1.19864916  2.94181986  0.4075   0.68368    
## merger       1.92231907  4.04770404  0.4749   0.63485    
## employ       0.12548448  0.03060149  4.1006 4.121e-05 ***
## return       0.06432167  0.10374558  0.6200   0.53526    
## patentsg     0.78696226  0.01016726 77.4016 < 2.2e-16 ***
## stckpr       0.00355791  0.02557045  0.1391   0.88934    
## rnd         -0.18291882  0.04480367 -4.0827 4.452e-05 ***
## rndeflt      0.26805014  0.03877619  6.9128 4.753e-12 ***
## rndstck     -0.00122890  0.00628664 -0.1955   0.84502    
## sales       -0.00054529  0.00025769 -2.1161   0.03434 *  
## sic         -0.00049485  0.00081918 -0.6041   0.54579    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    10910000
## Residual Sum of Squares: 1154800
## R-Squared:      0.89416
## Adj. R-Squared: 0.89365
## Chisq: 17504.4 on 10 DF, p-value: < 2.22e-16

#Prueba de Hausman
phtest(modelo_efecto_fijo_patentes, modelo_efecto_aleatorio_patentes)

## 
##  Hausman Test
## 
## data:  patents ~ merger + employ + return + patentsg + stckpr + rnd +  ...
## chisq = 1104.9, df = 9, p-value < 2.2e-16
## alternative hypothesis: one model is inconsistent

  #Como el p-value es mayor a 0.05, usamos el modelo de Efectos Aleatorios

Pueba de heterocedasticidad y autocorrelación

#HETEROCEDASTICIDAD
#Modelos Fijos
bptest(modelo_efecto_fijo_patentes)

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_efecto_fijo_patentes
## BP = 617.25, df = 10, p-value < 2.2e-16

  #Como el p-value es < 0.05, hay heterocedasticidad en los residuos.

#Modelo Aleatorio
bptest(modelo_efecto_aleatorio_patentes)

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_efecto_aleatorio_patentes
## BP = 617.25, df = 10, p-value < 2.2e-16

  #Como el p-value es < 0.05, hay heterocedasticidad en los residuos.

# AUTOCORRELACIÓN
#Modelos Efecto Fijo
pwartest(modelo_efecto_fijo_patentes)

## 
##  Wooldridge's test for serial correlation in FE panels
## 
## data:  modelo_efecto_fijo_patentes
## F = 42.281, df1 = 1, df2 = 1866, p-value = 1.012e-10
## alternative hypothesis: serial correlation

    #Como el p-value es < 0.05, hay autocorrelación en los residuos.

#Modelos Efecto Aleatorio
pbnftest(modelo_efecto_aleatorio_patentes)

## 
##  modified Bhargava/Franzini/Narendranathan Panel Durbin-Watson Test
## 
## data:  patents ~ merger + employ + return + patentsg + stckpr + rnd +  ...
## DW = 1.0069
## alternative hypothesis: serial correlation in idiosyncratic errors

  #Como el p-value es menos a 1.5 hay autocorrelación serial. 


#Corrección del Modelo con Errores Estándar Robustos
coeficientes_corregidos <- coeftest(modelo_efecto_fijo_patentes, vcov = vcovHC(modelo_efecto_fijo_patentes, type = "HC0"))

solo_coeficientes <- coeficientes_corregidos[, 1]
solo_coeficientes

##       merger       employ       return     patentsg       stckpr          rnd 
##  6.024679976 -0.090955344 -0.012214441  0.039139069 -0.039597715 -2.041010028 
##      rndeflt      rndstck        sales 
##  3.253694087  0.197241659 -0.001889381

4. Generar pronóstico y evaluar modelo

datos_de_prueba <- data.frame(merger=0, employ=10, return=6, patentsg = 24, stckpr = 48, rnd = 3, rndeflt = 3, rndstck = 16, sales = 344)

prediccion <- sum(solo_coeficientes*c( datos_de_prueba$merger, datos_de_prueba$employ, datos_de_prueba$return, datos_de_prueba$patentsg, datos_de_prueba$stckpr, datos_de_prueba$rnd, datos_de_prueba$rndeflt, datos_de_prueba$rndstck, datos_de_prueba$sales))

prediccion

## [1] 4.199779