# Poner el directorio de trabajo con la función setwd() # Carga y exploración de datos #### cereales <- read.table("Cereales.txt", header=T, dec=".") # Cargamos los datos str(cereales) # Chequeamos la estructura de los datos op <- par(mfrow=c(2,1), cex=0.95) with(cereales, plot(capital, cereal, pch=16, col="darkgray", xlab="Stock de capital (millones de USD$)", ylab="Producción de cereales (ton)", cex.axis=0.9)) with(cereales, plot(trabajo, cereal, pch=16, col="darkgray", xlab="PEA en el agro (miles de personas)", ylab="Producción de cereales (ton)", cex.axis=0.9)) par(op) # INFERENCIA MULTIMODELO SEGÚN FIGURA 2 DE LA AYUDA DIDÁCTICA # (1) Planteamos múltiples modelos ##### # Ver sección "Producción de cereales" en el texto principal de la ayuda didáctica # Mnulo = Modelo nulo # M0.si = Modelo lineal sin interacción # M1.ci = Modelo lineal con interacción # M2.CD = Modelo Cobb Douglass # M3.ESC = Modelo de Elasticidad de Sustitución Constante # (2) Estimamos los parámetros de cada modelo #### # Dado que evaluamos modelos no-lineales el ajuste lo llevamos adelante con la función nls # Para saber más acerca de esta función podemos escribir en r ?nls # a0 - Modelo Nulo Mnulo <- nls(cereal ~ B0, start=list(B0=1), data = cereales) summary(Mnulo) # a1 - Modelo lineal sin interacción M0.si <- nls(cereal ~ B0 + B1 * capital + B2 * trabajo, start=list(B0=1, B1=1, B2=1), data=cereales) summary(M0.si) # a2 - Modelo lineal con interacción M1.ci <- nls(cereal ~ B0 + B1 * capital + B2 * trabajo + B3 * capital * trabajo, start=list(B0=1, B1=1, B2=1, B3=1), data=cereales) summary(M1.ci) # a3 - Modelo de Cobb Douglass M2.CD <- nls(cereal ~ A * (capital^B1 + trabajo^B2), start=list(A=1000, B1=1, B2=1), data=cereales) summary(M2.CD) # a4 - Modelo de Elasticidad de Sustitución Constante M3.ESC <- nls(cereal ~ A * (B * capital^(-r) + (1 - B) * trabajo^(-r))^(-v/r), data = cereales, start = c(A=800, B=0.98, r=1, v=1)) summary(M3.ESC) # (3) Verificamos los supuestos #### resM0.si <- resid(M0.si) # residuos del modelo lineal sin interacción predM0.si <- fitted(M0.si) # predichos del modelo lineal sin interacción resM1.ci <- resid(M1.ci) # residuos del modelo lineal con interacción predM1.ci <- fitted(M1.ci) # predichos del modelo lineal con interacción resM2.CD <- resid(M2.CD) # residuos del modelo de Cobb Douglass predM2.CD <- fitted(M2.CD) # predichos del modelo de Cobb Douglass resM3.ESC <- resid(M3.ESC) # residuos del modelo de Elasticidad de Sustitución Constante predM3.ESC <- fitted(M3.ESC) # predichos del modelo de Elasticidad de Sustitución Constante # 3.1 - Evaluamos la homocedasticidad, linealidad e independecia # Gráfico de residuos vs. predichos op <- par(mfrow=c(2,2)) plot(predM0.si, resM0.si, ylab="Residuos", xlab="Predichos", main="M0.si") abline(a=0,b=0, lw=2) plot(predM1.ci, resM1.ci, ylab="Residuos", xlab="Predichos", main="M1.ci") abline(a=0,b=0, lw=2) plot(predM2.CD, resM2.CD, ylab="Residuos", xlab="Predichos", main="M2.CD") abline(a=0,b=0, lw=2) plot(predM3.ESC, resM3.ESC, ylab="Residuos", xlab="Predichos", main="M3.ESC") abline(a=0,b=0, lw=2) par(op) # Gráfico de residuos en función del capital op <- par(mfrow=c(2,2)) plot(cereales$capital, resM0.si, ylab="Residuos", xlab="Stock de capital (millones de USD$)", main="M0.si") abline(a=0,b=0, lw=2) plot(cereales$capital, resM1.ci, ylab="Residuos", xlab="Stock de capital (millones de USD$)", main="M1.ci") abline(a=0,b=0, lw=2) plot(cereales$capital, resM2.CD, ylab="Residuos", xlab="Stock de capital (millones de USD$)", main="M2.CD") abline(a=0,b=0, lw=2) plot(cereales$capital, resM3.ESC, ylab="Residuos", xlab="Stock de capital (millones de USD$)", main="M3.ESC") abline(a=0,b=0, lw=2) par(op) # Gráfico de residuos en función del trabajo op <- par(mfrow=c(2,2)) plot(cereales$trabajo, resM0.si, ylab="Residuos", xlab="PEA en el agro (miles de personas)", main="M0.si") abline(a=0,b=0, lw=2) plot(cereales$trabajo, resM1.ci, ylab="Residuos", xlab="PEA en el agro (miles de personas)", main="M1.ci") abline(a=0,b=0, lw=2) plot(cereales$trabajo, resM2.CD, ylab="Residuos", xlab="PEA en el agro (miles de personas)", main="M2.CD") abline(a=0,b=0, lw=2) plot(cereales$trabajo, resM3.ESC, ylab="Residuos", xlab="PEA en el agro (miles de personas)", main="M3.ESC") abline(a=0,b=0, lw=2) par(op) # 3.2 - Evaluamos la normalidad # Histogramas y diagramas de caja y bigotes op <- par(mfrow=c(4,2)) hist(resM0.si, col="grey", xlab="Residuos", ylab="Frecuencia", main="M0.si") hist(resM1.ci, col="grey", xlab="Residuos", ylab="Frecuencia", main="M1.ci") boxplot(resM0.si, bty="l", range=1.5, col="grey", horizontal=T, xlab="Residuos") boxplot(resM1.ci, bty="l", range=1.5, col="grey", horizontal=T, xlab="Residuos") hist(resM2.CD, col="grey", xlab="Residuos", ylab="Frecuencia", main="M2.CD") hist(resM3.ESC, col="grey", xlab="Residuos", ylab="Frecuencia", main="M3.ESC") boxplot(resM2.CD, bty="l", range=1.5, col="grey", horizontal=T, xlab="Residuos") boxplot(resM3.ESC, bty="l", range=1.5, col="grey", horizontal=T, xlab="Residuos") par(op) # Gráficos cuantil-cuantil op <- par(mfrow=c(2,2)) qqnorm(resM0.si, main="M0.si", ylab="Cuantiles observados", xlab="Cuantiles teóricos") qqline(resM0.si) qqnorm(resM1.ci, main="M1.ci", ylab="Cuantiles observados", xlab="Cuantiles teóricos") qqline(resM1.ci) qqnorm(resM2.CD, main="M2.CD", ylab="Cuantiles observados", xlab="Cuantiles teóricos") qqline(resM2.CD) qqnorm(resM3.ESC, main="M3.ESC", ylab="Cuantiles observados", xlab="Cuantiles teóricos") qqline(resM3.ESC) par(op) # Asimetría y curtósis en los residuos install.packages("moments") # Instalamos el paquete moments library(moments) # Cargamos el paquete moments skewness(resM0.si) # la asimetría de una distribución normal es cero kurtosis(resM0.si) # la curtósis de una distribución normal es tres skewness(resM1.ci) kurtosis(resM1.ci) skewness(resM2.CD) kurtosis(resM2.CD) skewness(resM3.ESC) kurtosis(resM3.ESC) # Evaluamos normalidad mediante test de hipótesis (H0 supone normalidad) # Test de Kolmogorov-Smirnov ks.test(resM0.si, "pnorm", mean(resM0.si), sd(resM0.si)) ks.test(resM1.ci, "pnorm", mean(resM1.ci), sd(resM1.ci)) ks.test(resM2.CD, "pnorm", mean(resM2.CD), sd(resM2.CD)) ks.test(resM3.ESC,"pnorm", mean(resM3.ESC), sd(resM3.ESC)) # Modificacion de Lillefors al test de Kolmogorov-Smirnov install.packages("nortest") # Instalamos el paquete nortest library(nortest) # Cargamos el paquete nortest lillie.test(resM0.si) lillie.test(resM1.ci) lillie.test(resM2.CD) lillie.test(resM3.ESC) # Test de Shapiro-Wilk shapiro.test(resM0.si) shapiro.test(resM1.ci) shapiro.test(resM2.CD) shapiro.test(resM3.ESC) # Finalmente, graficamos los valores observados vs valores predichos op <- par(mfrow=c(2,2)) plot(cereales$cereal, predM0.si, ylab="Observados", xlab="Predichos", main="M0.si") abline(a=0,b=1, lw=2) plot(cereales$cereal, predM1.ci, ylab="Observados", xlab="Predichos", main="M1.ci") abline(a=0,b=1, lw=2) plot(cereales$cereal, predM2.CD, ylab="Obserrvados", xlab="Predichos", main="M2.CD") abline(a=0,b=1, lw=2) plot(cereales$cereal, predM3.ESC, ylab="Residuos", xlab="PEA en el agro (miles de personas)", main="M3.ESC") abline(a=0,b=1, lw=2) par(op) # (4) Ordenamos los modelos según AIC #### AIC(Mnulo, M0.si, M1.ci, M2.CD, M3.ESC) # También podemos usar la función AICtab del paquete bblme que ordena # los modelos según AIC y devuelve valores interesantes como la verosimlitud # y las diferencias de AIC entre modelos. install.packages("bbmle") # instalamos el paquete bbmle library("bbmle") # cargamos el paquete bbmle AICtab(Mnulo, M0.si, M1.ci, M2.CD, M3.ESC, weights = T, delta = TRUE, base=T, sort = TRUE, logLik = TRUE) # (5) Presentamos los modelos #### # Gráficos de los modelos ajustados op <- par(mfrow=c(1, 2)) # Factor Capital with(cereales, plot(capital, cereal, col = "darkgrey", pch = 16, xlab="Stock de capital (millones de USD$)", ylab="Producción de cereales (ton)", main="Factor Capital", cex.main=0.95, cex.axis=0.8, bty="l")) abline(a=coef(Mnulo)[1], b=0, lw=1) # Nulo curve(coef(M0.si)[1] + coef(M0.si)[2] * x + coef(M0.si)[3] * mean(cereales$trabajo), lwd=1, lty=3, add=T) # Lineal sin interacción curve(coef(M1.ci)[1] + coef(M1.ci)[2] * x + coef(M1.ci)[3] * mean(cereales$trabajo) + coef(M1.ci)[4] * x * mean(cereales$trabajo) , lwd=1, lty=2, add=T) # Lineal con interacción curve(coef(M2.CD)[1] * (x^(coef(M2.CD)[2]) + mean(cereales$trabajo)^(coef(M2.CD)[3])), lwd=2, lty=1, add=T) # Cobb Douglass curve( coef(M3.ESC)[1] * (coef(M3.ESC)[2] * x^(-coef(M3.ESC)[3]) + (1 - coef(M3.ESC)[2]) *mean(cereales$trabajo)^(-coef(M3.ESC)[3]))^(-coef(M3.ESC)[4]/coef(M3.ESC)[3]), lwd=1, lty=5, add=T) # ESC legend("topleft",legend=c("Nulo", "Lineal sin interacción", "Lineal con interacción", "Cobb-Douglass", "ESC"), lty=c(1, 3, 2, 1, 5), lwd=c(1, 1, 1, 2, 1), bty = "n", x.intersp=1, y.intersp=1.15, cex = 0.8) # Factor Trabajo with(cereales, plot(trabajo, cereal, col = "darkgrey", pch = 16, xlab="PEA en el agro (miles de personas)", ylab="", main="Factor Trabajo", cex.main=0.95, cex.axis=0.8, bty="l")) abline(a=coef(Mnulo)[1], b=0, lw=1) # Nulo curve(coef(M0.si)[1] + coef(M0.si)[2] * mean(cereales$capital) + coef(M0.si)[3] * x, lwd=1, lty=3, add=T) # Lineal sin interacción curve(coef(M1.ci)[1] + coef(M1.ci)[2] * mean(cereales$capital) + coef(M1.ci)[3] * x + coef(M1.ci)[4] * mean(cereales$capital) * x, lwd=1, lty=2, add=T) # Lineal con interacción curve(coef(M2.CD)[1] * (mean(cereales$capital)^(coef(M2.CD)[2]) + x^(coef(M2.CD)[3])), lwd=2, lty=1, add=T) # Cobb Douglass curve(4720.4 * ((0.424 * mean(cereales$capital)^0.888) + (1-0.424)* x^0.888)^(0.836/0.888), lwd=1, lty=5, add=T) curve( coef(M3.ESC)[1] * (coef(M3.ESC)[2] * mean(cereales$capital)^(-coef(M3.ESC)[3]) + (1 - coef(M3.ESC)[2]) *x^(-coef(M3.ESC)[3]))^(-coef(M3.ESC)[4]/coef(M3.ESC)[3]), lwd=1, lty=5, add=T) # ESC legend("topleft",legend=c("Nulo", "Lineal sin interacción", "Lineal con interacción", "Cobb-Douglass", "ESC"), lty=c(1, 3, 2, 1, 5), lwd=c(1, 1, 1, 2, 1), bty = "n", x.intersp=1, y.intersp=1.15, cex = 0.8) par(op)