Activdad 1. Análisis y aplicación de datos panel

Librerias

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

1. Código sesiones 2, 3 y 4

Analisis de Datos de Panel

Información de paises e indicadores

### A. Obtener la información de 1 país
PIB_MEX <- wb_data(country = "MX", indicator = "NY.GDP.PCAP.CD",
                   start_date = 1900, end_date = 2025)

summary(PIB_MEX)

##     iso2c              iso3c             country               date     
##  Length:64          Length:64          Length:64          Min.   :1960  
##  Class :character   Class :character   Class :character   1st Qu.:1976  
##  Mode  :character   Mode  :character   Mode  :character   Median :1992  
##                                                           Mean   :1992  
##                                                           3rd Qu.:2007  
##                                                           Max.   :2023  
##  NY.GDP.PCAP.CD        unit            obs_status          footnote        
##  Min.   :  355.1   Length:64          Length:64          Length:64         
##  1st Qu.: 1427.8   Class :character   Class :character   Class :character  
##  Median : 4006.5   Mode  :character   Mode  :character   Mode  :character  
##  Mean   : 5097.1                                                           
##  3rd Qu.: 8905.4                                                           
##  Max.   :13790.0                                                           
##   last_updated       
##  Min.   :2025-01-28  
##  1st Qu.:2025-01-28  
##  Median :2025-01-28  
##  Mean   :2025-01-28  
##  3rd Qu.:2025-01-28  
##  Max.   :2025-01-28

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")

### B. Obtener la información de varios países
PIB_PANEL <- wb_data(country = c("MX","US", "CA"), indicator = "NY.GDP.PCAP.CD",
                   start_date = 1900, end_date = 2025)
summary(PIB_PANEL)

##     iso2c              iso3c             country               date     
##  Length:192         Length:192         Length:192         Min.   :1960  
##  Class :character   Class :character   Class :character   1st Qu.:1976  
##  Mode  :character   Mode  :character   Mode  :character   Median :1992  
##                                                           Mean   :1992  
##                                                           3rd Qu.:2007  
##                                                           Max.   :2023  
##  NY.GDP.PCAP.CD        unit            obs_status          footnote        
##  Min.   :  355.1   Length:192         Length:192         Length:192        
##  1st Qu.: 4059.2   Class :character   Class :character   Class :character  
##  Median :10544.4   Mode  :character   Mode  :character   Mode  :character  
##  Mean   :19152.2                                                           
##  3rd Qu.:29010.1                                                           
##  Max.   :82769.4                                                           
##   last_updated       
##  Min.   :2025-01-28  
##  1st Qu.:2025-01-28  
##  Median :2025-01-28  
##  Mean   :2025-01-28  
##  3rd Qu.:2025-01-28  
##  Max.   :2025-01-28

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

### C. Obtener la información de varios indicadores en varios países

MEGAPIB <- wb_data(country = c("MX","US", "CA"), indicator = c("NY.GDP.PCAP.CD", "SP.DYN.LE00.IN"),
                   start_date = 1900, end_date = 2025)

summary(MEGAPIB)

##     iso2c              iso3c             country               date     
##  Length:192         Length:192         Length:192         Min.   :1960  
##  Class :character   Class :character   Class :character   1st Qu.:1976  
##  Mode  :character   Mode  :character   Mode  :character   Median :1992  
##                                                           Mean   :1992  
##                                                           3rd Qu.:2007  
##                                                           Max.   :2023  
##                                                                         
##  NY.GDP.PCAP.CD    SP.DYN.LE00.IN 
##  Min.   :  355.1   Min.   :55.02  
##  1st Qu.: 4059.2   1st Qu.:71.11  
##  Median :10544.4   Median :74.36  
##  Mean   :19152.2   Mean   :73.41  
##  3rd Qu.:29010.1   3rd Qu.:77.49  
##  Max.   :82769.4   Max.   :82.22  
##                    NA's   :3

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

Interpretación: - Alta Heterogenidad: Si los puntos (medias) están muy separados entre países. - Baja Heterogenidad: Si los puntos (medias) están muy cerca uno de otros. - En este caso, EUA y Canadá tienen un PIB per Cápita mayor que México, mostrando alta heterogenidad entre países.

Modelos y Prueba de Hausman

datos_panel <- pdata.frame(PIB_PANEL, index = c("country", "date"))

#Modelo de efectos fijos
modelo_efectos_fijos <- plm(NY.GDP.PCAP.CD ~ date, data=datos_panel, model = "within")

#Modelo de efectos aleatorios
modelo_efectos_aleatorios <- plm(NY.GDP.PCAP.CD ~ date, data=datos_panel, model ="random")

#Prueba de hausman
phtest(modelo_efectos_fijos, modelo_efectos_aleatorios)

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

#Como el p-value es mayor a 0.05, usamos el modelo de efectos aleatorios

2. Link Aplicación Shiny

[Indicadores Banco Mundial] (https://a00833172.shinyapps.io/PIBB/)

3. Código Sesión 5

Actividad Patentes

Importar Base de Datos

patentes <- read_excel("/Users/genarorodriguezalcantara/Desktop/Tec/Generacion de escenarios futuros con analítica (Gpo 101)/PIB/PATENT 3.xls")

Entender 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)) #contas NAs

## [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)

A. Construcción del modelo de datos en panel

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

B. Modelo de Efectos Fijos y Aleatorios

Modelo de efectos fijos

modelo_efectos_fijos_patentes <- plm(patents ~ merger + employ + return + patentsg + stckpr + rnd + rndeflt + sales + sic, data = panel_patentes, model = "within")

Modelo de efectos aleatorios

modelo_efectos_aleatorios_patentes <- plm(patents ~ merger + employ + return + patentsg + stckpr + rnd + rndeflt + sales + sic, data = panel_patentes, model ="random")

Prueba de hausman

phtest(modelo_efectos_fijos_patentes, modelo_efectos_aleatorios_patentes)

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

C. Pruebas de heterocedasticidad y autocorrelación serial

Prueba de Heterocedasticidad para el model de efectos fijos

bptest(modelo_efectos_aleatorios_patentes)

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_efectos_aleatorios_patentes
## BP = 616, df = 9, p-value < 2.2e-16

# como el pvalue es < 0.05, hay heterocedasticidad en los residuos

Prueba de Heterocedasticidad para el model de efectos aleatorios

bptest(modelo_efectos_aleatorios_patentes)

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_efectos_aleatorios_patentes
## BP = 616, df = 9, p-value < 2.2e-16

# como el pvalue es < 0.05, hay heterocedasticidad en los residuos

pwartest(modelo_efectos_fijos_patentes)

## 
##  Wooldridge's test for serial correlation in FE panels
## 
## data:  modelo_efectos_fijos_patentes
## F = 73.741, df1 = 1, df2 = 1866, p-value < 2.2e-16
## alternative hypothesis: serial correlation

# como el pvalue es < 0.05, hay autocorrelación serial en los errores

pbnftest(modelo_efectos_aleatorios_patentes)

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

# como el valor es de < 1.5, hay autocorrelación positiva significativa

D. Correción del modelo con errores estandar robustos

coeficientes_corregidos <- coeftest(modelo_efectos_fijos_patentes, vcov=vcovHC(modelo_efectos_fijos_patentes, type = "HC0"))

solo_coeficientes <- coeficientes_corregidos[,1]

E. Generar Pronósticos 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)

prediction <- 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$sales))

prediction

## [1] -5.620776