A continuación abordaremos la estimación de un modelo no lineal. En primer lugar, estudiaremos el caso de un modelo de regresión lineal múltiple (el caso simple ya ha sido estudiado en este blog) que se puede linealizar mediante transformaciones sencillas. Mientras que en segundo lugar, supondremos que no es factible la linealización del modelo, por lo que habrá que buscar una alternativa para la estimación.
Los datos con los que vamos a trabajar se cargan mediante las órdenes:
datos = read.table("ProduccionAgropecuaria.csv", header=T, sep=";")
names(datos) # vemos el nombre de cada una de las variables[1] "Anos" "P" "K" "W" # P = producción total
# K = capital
# W = trabajo
attach(datos) # incorporamos los datos a R para referenciar cada variable por el nombre anteriorEstos datos corresponden a la producción agropecuaria total (medida en toneladas), P, al capital (medido a partir del PIB en miles de pesos con año base 1993), K, y al trabajo (medido como el número de empleados, más concretamente, el promedio de personal ocupado remunerado), T, de un modelo de producción agropecuario para México a partir de la información para los años 1980-2007.
Considerando que la relación entre las variables corresponde a una función de producción de Cobb-Douglas del tipo \(P = a \cdot W^b \cdot K^c \cdot e^u\), se tiene que podemos estimar su versión linealizada \(log(P) = log(a) + b \cdot log(T) + c \cdot log(K) + u\). Para ello recurrimos al código lm:
regresion1 = lm(log(P)~log(W)+log(K))
summary(regresion1)
Call:
lm(formula = log(P) ~ log(W) + log(K))
Residuals:
Min 1Q Median 3Q Max
-0.078440 -0.025424 0.007202 0.025985 0.064056
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -10.38590 3.32728 -3.121 0.0045 **
log(W) 0.43576 0.25436 1.713 0.0991 .
log(K) 1.17094 0.08467 13.830 3.24e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.04163 on 25 degrees of freedom
Multiple R-squared: 0.9304, Adjusted R-squared: 0.9248
F-statistic: 167.1 on 2 and 25 DF, p-value: 3.402e-15Se puede observar que el coeficiente del término independiente y del capital son significativamente distintos de cero y que el modelo es válido en su conjunto. Puesto que los coeficientes se interpretan como elasticidades, se tiene que cuando el capital aumenta un 1% se tiene que la producción aumenta un 1.17%.
En la siguiente representación gráfica se tienen los valores observados y ajustados:
plot(P ~ Anos, xlab="Año", ylab="Producción agropecuaria", col="blue")
lines(Anos, fitted(regresion1), col="red", lwd=2)
Por otro lado, la estimación del modelo original sería \(P = 0.0000308 \cdot W^{0.43576} \cdot K^{1.17094}\).
Si consideramos que la perturbación aleatoria no se añade de forma multiplicativa sino aditiva se tendría que P = a * T^(b) * K^(c) + u y, en tal caso, no se puede linealizar el modelo facilmente teniendo que recurrir a la estimación de modelos no lineales. Para ello recurrimos al código nls (del paquete stats) que forma parte del paquete base de R:
regresion2 = nls(P ~ a*(W^(b))*(K^(c)), start=list(a=0.001, b=0.5, c=1))
summary(regresion2)
Formula: P ~ a * (W^(b)) * (K^(c))
Parameters:
Estimate Std. Error t value Pr(>|t|)
a 2.239e-05 7.108e-05 0.315 0.7554
b 4.606e-01 2.271e-01 2.028 0.0533 .
c 1.167e+00 7.256e-02 16.088 1.07e-14 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1888000 on 25 degrees of freedom
Number of iterations to convergence: 8
Achieved convergence tolerance: 1.217e-07Adviértase que como valores iniciales se han considerado una aproximación de las estimaciones obtenidas en la regresión anterior, obteniéndose unas estimaciones muy parecidas a las anteriores. Por defecto, se usa para la estimación el algoritmo de Gauss-Newton. Para más información ver help(nls).
Una vez más se podrían representar los valores observados y ajustados:
plot(P ~ Anos, xlab="Año", ylab="Producción agropecuaria", col="blue")
lines(Anos, fitted(regresion2), col="red", lwd=2)
Algunas direcciones de internet interesantes sobre este tema son: