Cargamos las principales bibliotecas que vamos a emplear en nuestro análisis:
library(reshape2)
library(dplyr)
library(car) # regression analysis
library(plotly)
library(ggplot2)
library(lmtest) # lm Breusch-Godfrey test for serial correlation
library(ecm) # to calculate Durbin's h-statistic
library(orcutt) # cochrane-orcutt procedure in Case of First Order Autocorrelation
library(dynlm) # dynamic models
library(AER) # ivreg
library(stargazer) # regression tableExplicar brevemente la funcionalidad de esto
knitr::opts_chunk$set(echo = TRUE, warning = FALSE)Considere las siguientes variables para especificar un modelo econométrico para estimar un sistema de ecuaciones simultáneas:
\(Y_t\) = Nivel de producción (variable GDP = rgdpna).
\(X_{2t}\) = Factor de producción trabajo (variable L = avh*emp).
\(X_{3t}\) = Factor de producción capital físico (variable K = rnna).
\(X_{4t}\) = Patentes (Solicitudes de patentes a la EPO por año de prioridad).
\(X_{5t}\) = Capital Humano (Índice de capital humano por persona, basado en años de escolaridad y retornos a la educación).
data = read.csv("~/Desktop/Scripts de R/Datos/FebPwtExport2232023 19.41.35.csv")
head(data, 10)data_wide <- dcast(data, RegionCode + YearCode ~ VariableCode, value.var = "AggValue")
head(data_wide, 10)Renombramos las variables y creamos una nueva:
data_wide <- data_wide %>%
dplyr::mutate(l=avh*emp) %>%
dplyr::rename(Year = YearCode,
pib = rgdpna,
k = rnna)
names(data_wide)## [1] "RegionCode" "Year" "avh" "emp" "pib"
## [6] "k" "l"head(data_wide, 10)Importamos los datos correspondientes al capital humano:
data_kh <- read.csv("~/Desktop/Scripts de R/Datos/información_CH.csv")
head(data_kh, 10)data_kh <- dcast(data_kh, RegionCode + YearCode ~ VariableCode, value.var = "AggValue") data_kh <- data_kh %>%
dplyr::rename(kh = hc,
Year = YearCode) %>%
dplyr::select(-RegionCode)
head(data_kh, 10)data_pat <- read.csv("~/Desktop/Scripts de R/Datos/información_patentes.csv")
head(data_pat, 10)data_pat = data_pat %>%
dplyr::select(TIME_PERIOD, OBS_VALUE) %>%
dplyr::rename(Patentes = OBS_VALUE,
Year = TIME_PERIOD)
head(data_pat, 10)data_wide <- full_join(data_wide, data_kh, by = c("Year"))
head(data_wide, 10)data_wide <- full_join(data_wide, data_pat, by = c("Year"))
head(data_wide, 10)\[\log(pib_t) = \beta_1 + \beta_2 \log(capital_{t}) + \beta_3 \log(trabajo_{t}) + u_{1t} + \beta_7 \log(patentes_t) + \beta_8 \log(kh_t)\] \[\log(patentes_t) = \beta_4 + \beta_5 log(pib_t) + \beta_6 \log(kh_t) + u_{2t}\]
El MES (modelo de ecuaciones simultáneas) que hemos presentado está escrito en su forma estructural, por lo que podríamos escribir el MES en su forma reducida con el objetivo de poder estimar sus ecuaciones, ahora bien, antes de esto hemos de aplicar a nuestras ecuaciones las condiciones de Orden y de Rango, con el objetivo de ver la identificación de las mismas.
\[\log(pib_t) - \beta_7 \log(patentes_t) = \beta_1 + \beta_2 \log(capital_{t}) + \beta_3 \log(trabajo_{t}) + \beta_8 \log(kh_t) + u_{1t}\] \[ - \beta_5 log(pib_t) + \log(patentes_t) = \beta_4 + \beta_6\log(kh_t) + u_{2t}\]
\[\begin{equation} \begin{bmatrix} 1 & -\beta_7\\ - \beta_5 & 1 \end{bmatrix} \begin{bmatrix} \log(pib_t) \\ \log(patentes_t) \end{bmatrix} = \begin{bmatrix} \beta_1 & \beta_2 & \beta_3 & \beta_8\\ \beta_4 & 0 & 0 & \beta_6 \end{bmatrix}\begin{bmatrix} 1 \\ \log(capital_{t})\\ \log(trabajo_{t})\\ \log(kh_t) \end{bmatrix} + \begin{bmatrix} u_{1t}\\ u_{2t} \end{bmatrix} \end{equation}\]
HACER, EXPLICAR LOS MÉTODOS Y EL QUE VAMOS A USAR (MC2E)
Recordemos cual era la ecuación 2 (eq.2):
\[\log(patentes_t) = \beta_4 + \beta_5 log(pib_t) + \beta_6 \log(kh_t) + u_{2t}\]
Estimamos la segunda ecuación por MC2E (conjuntamente, utilizando la función ivreg()):
eq.2_hat <- ivreg(formula = log(Patentes) ~ log(pib) + log(kh) | log(kh) + log(k) + log(l),
data = data_wide)
summary(eq.2_hat, diagnostics = TRUE)##
## Call:
## ivreg(formula = log(Patentes) ~ log(pib) + log(kh) | log(kh) +
## log(k) + log(l), data = data_wide)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.166823 -0.032855 -0.005174 0.040154 0.178989
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.5604 3.2671 4.763 3.09e-05 ***
## log(pib) -2.0236 0.3549 -5.702 1.74e-06 ***
## log(kh) 16.1206 1.2238 13.173 2.46e-15 ***
##
## Diagnostic tests:
## df1 df2 statistic p-value
## Weak instruments 2 35 197.757 <2e-16 ***
## Wu-Hausman 1 35 3.576 0.0669 .
## Sargan 1 NA 2.592 0.1074
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07568 on 36 degrees of freedom
## Multiple R-Squared: 0.9863, Adjusted R-squared: 0.9856
## Wald test: 1301 on 2 and 36 DF, p-value: < 2.2e-16Podemos ecribir el siguiente chunk con la función stargazer de manera que los resultados queden recogidos de manera más amigable:
summ.fit <- summary(eq.2_hat , diagnostics = TRUE)
stargazer(eq.2_hat , type ="text",
add.lines = list(c(rownames(summ.fit$diagnostics)[1],
round(summ.fit$diagnostics[1, "statistic"], 3)),
c("p-value", round(summ.fit$diagnostics[1, "p-value"], 3)),
c(rownames(summ.fit$diagnostics)[2],
round(summ.fit$diagnostics[2, "statistic"], 3)),
c("p-value", round(summ.fit$diagnostics[1, "p-value"], 3)),
c(rownames(summ.fit$diagnostics)[3],
round(summ.fit$diagnostics[3, "statistic"], 3)),
c("p-value", round(summ.fit$diagnostics[1, "p-value"], 3))
))##
## ===============================================
## Dependent variable:
## ---------------------------
## log(Patentes)
## -----------------------------------------------
## log(pib) -2.024***
## (0.355)
##
## log(kh) 16.121***
## (1.224)
##
## Constant 15.560***
## (3.267)
##
## -----------------------------------------------
## Weak instruments 197.757
## p-value 0
## Wu-Hausman 3.576
## p-value 0
## Sargan 2.592
## p-value 0
## Observations 39
## R2 0.986
## Adjusted R2 0.986
## Residual Std. Error 0.076 (df = 36)
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01Podemos obtener esa misma estimación realizando las dos etapas nosotros mismos como tradicionalmente lo hacemos, ahora bien, es importante considerar que si buscamos arrojar los mismos resultados que hemos obtenido previaemnte, debemos eliminar las observaciones incompletas que pudiera haber en los datos, en nuestro caso, sabemos que esas observaciones incompletas se encuentran dentro de la variable “Patentes”:
eq.2_hat.1 <- lm(log(pib) ~ log(kh) + log(k) + log(l), data = na.omit(data_wide))
stargazer(eq.2_hat.1, type = "text")##
## ===============================================
## Dependent variable:
## ---------------------------
## log(pib)
## -----------------------------------------------
## log(kh) 1.947***
## (0.205)
##
## log(k) 0.264***
## (0.062)
##
## log(l) 0.482***
## (0.075)
##
## Constant 2.700***
## (0.331)
##
## -----------------------------------------------
## Observations 39
## R2 0.998
## Adjusted R2 0.998
## Residual Std. Error 0.011 (df = 35)
## F Statistic 5,966.992*** (df = 3; 35)
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01Debemos eliminar las observaciones incompletas para poder trabajar correctamente con los datos:
obs_completas = na.omit(data_wide)
nrow(obs_completas) ## [1] 39eq.2_hat.2 <- lm(log(Patentes) ~ fitted(eq.2_hat.1) + log(kh), data = na.omit(data_wide))
stargazer(eq.2_hat.2, type = "text")##
## ===============================================
## Dependent variable:
## ---------------------------
## log(Patentes)
## -----------------------------------------------
## fitted(eq.2_hat.1) -2.024***
## (0.337)
##
## log(kh) 16.121***
## (1.162)
##
## Constant 15.560***
## (3.101)
##
## -----------------------------------------------
## Observations 39
## R2 0.988
## Adjusted R2 0.987
## Residual Std. Error 0.072 (df = 36)
## F Statistic 1,443.658*** (df = 2; 36)
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01Recuperamos el código escrito previamente para poder obtener conclusiones viendo los resultados del test de endogeneidad de DW-Hausman, después veremos el test de validez de Sargan (para el que tendremos que escribir código adicional):
summ.fit <- summary(eq.2_hat , diagnostics = TRUE)
stargazer(eq.2_hat , type ="text",
add.lines = list(c(rownames(summ.fit$diagnostics)[1],
round(summ.fit$diagnostics[1, "statistic"], 3)),
c("p-value", round(summ.fit$diagnostics[1, "p-value"], 3)),
c(rownames(summ.fit$diagnostics)[2],
round(summ.fit$diagnostics[2, "statistic"], 3)),
c("p-value", round(summ.fit$diagnostics[1, "p-value"], 3)),
c(rownames(summ.fit$diagnostics)[3],
round(summ.fit$diagnostics[3, "statistic"], 3)),
c("p-value", round(summ.fit$diagnostics[1, "p-value"], 3))
))##
## ===============================================
## Dependent variable:
## ---------------------------
## log(Patentes)
## -----------------------------------------------
## log(pib) -2.024***
## (0.355)
##
## log(kh) 16.121***
## (1.224)
##
## Constant 15.560***
## (3.267)
##
## -----------------------------------------------
## Weak instruments 197.757
## p-value 0
## Wu-Hausman 3.576
## p-value 0
## Sargan 2.592
## p-value 0
## Observations 39
## R2 0.986
## Adjusted R2 0.986
## Residual Std. Error 0.076 (df = 36)
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01EXPLICAR FUNCIONAMIENTO TEST DE HAUSMAMN E INTERPRETAR RESULTADOS OBTENIDOS
Para realizar el test de Sargan evaluamos los instrumentos que hemos empleado, por lo que el código quedaría de la siguiente manera:
model_sargan_test <- lm(residuals(eq.2_hat) ~ log(kh) + log(k) + log(l), data = obs_completas)
stargazer(model_sargan_test, type ="text")##
## ===============================================
## Dependent variable:
## ---------------------------
## residuals(eq.2_hat)
## -----------------------------------------------
## log(kh) 2.069
## (1.417)
##
## log(k) -0.661
## (0.428)
##
## log(l) 0.777
## (0.520)
##
## Constant 0.522
## (2.289)
##
## -----------------------------------------------
## Observations 39
## R2 0.066
## Adjusted R2 -0.014
## Residual Std. Error 0.074 (df = 35)
## F Statistic 0.831 (df = 3; 35)
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01EXPLICAR FUNCIONAMIENTO TEST DE SARGAN (COMO MARKO) E INTEPRETAR RESULTADOS OBTENIDOS
\[\log(pib_t) = \beta_1 + \beta_2 \log(capital_{t}) + \beta_3 \log(trabajo_{t}) + u_{1t} + \beta_7 \log(patentes_t)\] \[\log(patentes_t) = \beta_4 + \beta_5 log(pib_t) + \beta_6 \log(kh_t) + u_{2t}\]
HACER, USAR MATRICES
HACER, EXPLICAR LOS MÉTODOS Y EL QUE VAMOS A USAR (MC3E)
Para estimar el sistema completo por MC3E tenemos que estimar las dos ecuaciones conjuntamente, a la vez debemos matizar los instrumentos que vamos a emplear; utilizaremos la función systemfit(), que es una función de la que disponemos en R para poder estimar sistemas de ecuaciones simultáneas, asismismo especificaremos el método de estimación que estamos usando, en nuestro caso, como ya hemos mencionado, es el de los mínimos cuadrados en tres etapas (MC3E).
library(systemfit)
eq.1.question2 <- log(pib) ~ log(k) + log(l) + log(Patentes)
eq.2.question2 <- log(Patentes) ~ log(pib) + log(kh)
inst <- ~log(k) + log(l) + log(kh)
system <- list(eq.1.question2=eq.1.question2,
eq.2.question2=eq.2.question2)
MC3E_system <- systemfit(system, "3SLS", inst = inst, data = obs_completas)
summary(MC3E_system)##
## systemfit results
## method: 3SLS
##
## N DF SSR detRCov OLS-R2 McElroy-R2
## system 78 71 0.213743 1e-06 0.987537 0.996541
##
## N DF SSR MSE RMSE R2 Adj R2
## eq.1.question2 39 35 0.007536 0.000215 0.014674 0.996346 0.996032
## eq.2.question2 39 36 0.206207 0.005728 0.075683 0.986333 0.985574
##
## The covariance matrix of the residuals used for estimation
## eq.1.question2 eq.2.question2
## eq.1.question2 0.000199241 -0.000413382
## eq.2.question2 -0.000413382 0.005727962
##
## The covariance matrix of the residuals
## eq.1.question2 eq.2.question2
## eq.1.question2 0.000215325 -0.000492639
## eq.2.question2 -0.000492639 0.005727962
##
## The correlations of the residuals
## eq.1.question2 eq.2.question2
## eq.1.question2 1.00000 -0.44359
## eq.2.question2 -0.44359 1.00000
##
##
## 3SLS estimates for 'eq.1.question2' (equation 1)
## Model Formula: log(pib) ~ log(k) + log(l) + log(Patentes)
## Instruments: ~log(k) + log(l) + log(kh)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.1758179 0.3799320 3.09481 0.0038601 **
## log(k) 0.3917931 0.0550435 7.11789 2.6932e-08 ***
## log(l) 0.5667272 0.0946036 5.99054 7.9456e-07 ***
## log(Patentes) 0.1471380 0.0176789 8.32280 8.1910e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.014674 on 35 degrees of freedom
## Number of observations: 39 Degrees of Freedom: 35
## SSR: 0.007536 MSE: 0.000215 Root MSE: 0.014674
## Multiple R-Squared: 0.996346 Adjusted R-Squared: 0.996032
##
##
## 3SLS estimates for 'eq.2.question2' (equation 2)
## Model Formula: log(Patentes) ~ log(pib) + log(kh)
## Instruments: ~log(k) + log(l) + log(kh)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.560383 3.267061 4.76281 3.0945e-05 ***
## log(pib) -2.023633 0.354899 -5.70199 1.7390e-06 ***
## log(kh) 16.120571 1.223776 13.17281 2.4425e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.075683 on 36 degrees of freedom
## Number of observations: 39 Degrees of Freedom: 36
## SSR: 0.206207 MSE: 0.005728 Root MSE: 0.075683
## Multiple R-Squared: 0.986333 Adjusted R-Squared: 0.985574MC3E_system##
## systemfit results
## method: 3SLS
##
## Coefficients:
## eq.1.question2_(Intercept) eq.1.question2_log(k)
## 1.175818 0.391793
## eq.1.question2_log(l) eq.1.question2_log(Patentes)
## 0.566727 0.147138
## eq.2.question2_(Intercept) eq.2.question2_log(pib)
## 15.560383 -2.023633
## eq.2.question2_log(kh)
## 16.120571Vamos a comprobar si el coeficiente de la variable log(Patentes) en la primera ecuación es igual al coeficiente de la variable log(kh) en la segunda ecuación, la hipótesis nula es que ambos coeficientes son iguales, la alternativa es que son diferentes. Para realizar este test usamos la función linearHypothesis():
linearHypothesis(MC3E_system, c("eq.1.question2_log(Patentes)=eq.2.question2_log(kh)"),test="F")Podemos atender al p-value del estadístico de contraste F y, a la luz de su valor (muy cercano a 0), tenemos evidencia suficiente para rechazar la hipótesis nula en favor de la hipótesis alternativa, es decir, podemos concluir que la elasticidad del PIB con respecto a las patentes es significativamente diferente a la elasticidad de las patentes con respecto al capital humano.