ACT 1 - A01721951

#gif

Codigo sesiones 2, 3, y 4

#Instalar paquetes y llamar librerías
library(WDI)
library(wbstats)
library(ggplot2)
library(gplots)

## 
## Attaching package: 'gplots'

## The following object is masked from 'package:stats':
## 
##     lowess

library(plm)

#Obtener la informacioó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)

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

#Obtener la informacioó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")

#Obtener la informacioó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

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

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

#Interpretación:
# Alta Heterogenidad: Si losp untos (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.

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 = 1.8241e-12, 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

Link a Aplicación de Shiny

[Act 1 - Indicadores Banco Mundial] (https://jujqel-jenaro-martinez.shinyapps.io/UltimaApp/) [Act 2 - Indicadores Banco Mundial Economy and Growth] (https://jujqel-jenaro-martinez.shinyapps.io/econ_grow/)

Actividad de Patentes

#Instalar paquetes y llamar librerías
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.0     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ lubridate 1.9.4     ✔ tibble    3.2.1
## ✔ purrr     1.0.1     ✔ tidyr     1.3.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::between() masks plm::between()
## ✖ dplyr::filter()  masks stats::filter()
## ✖ dplyr::lag()     masks plm::lag(), stats::lag()
## ✖ dplyr::lead()    masks plm::lead()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

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

## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(pglm)

## Loading required package: maxLik

## Warning: package 'maxLik' was built under R version 4.2.3

## Loading required package: miscTools
## 
## Please cite the 'maxLik' package as:
## Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
## 
## If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
## https://r-forge.r-project.org/projects/maxlik/

#file.choose()
patentes <- read_excel("/Users/jenaromtzg/Desktop/Generación de Escenarios Futuros/R/PATENT 3.xls")

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)

## Construcción del modelo de datos en panel
panel_patentes <- pdata.frame(patentes1, index = c("cusip", "year"))

# Modelo de efectos Fijos y Aleatorios

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

#Modelo de efectos aleatorios
modelo_efectos_aleatorios_patentes <- plm(patents ~ merger + employ + return + patentsg + stckpr + rnd + rndstck + 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 = 1104.9, df = 9, p-value < 2.2e-16
## alternative hypothesis: one model is inconsistent

# 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 = 617.25, df = 10, 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 = 617.25, df = 10, 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 = 42.281, df1 = 1, df2 = 1866, p-value = 1.012e-10
## alternative hypothesis: serial correlation

# como el pvalue es < 0.05, hay autocorrelacion 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 = 1.0069
## alternative hypothesis: serial correlation in idiosyncratic errors

# como el valor es de < 1.5, hay autocorrelacion positiva significativa

#Correcion 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]

#Generar Pronosticos 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$sales))

## Warning in solo_coeficientes * c(datos_de_prueba$merger,
## datos_de_prueba$employ, : longer object length is not a multiple of shorter
## object length

prediccion

## [1] 69.54499