Plot de la función de producción Cobb-Douglas en R

El propósito de esta entrada, es mostrarles principalmente una, de las seguramente muchas formas que existen en R, para generar un gráfico en tres dimensiones de la función de producción Cobb-Douglas (linealizada en este caso).

Cabe destacar que no se requiere de paquetería especial alguna, pero que el código necesario para producir el gráfico es ligeramente extenso y es posible que nos perdamos en algún momento. Mi recomendación es que lo copies, y que trates de "desmenuzarlo" poco a poco, con el objetivo de comprenderlo y posteriormente realizar las modificaciones que se ajusten a tus necesidades.

El primer paso es la base de datos, para los que se encuentran familiarizados con la función de producción Cobb-Douglas resultará sencillo entender la estructura que debe tener la base de datos, para los que no, recomiendo que lean un poco al respecto a fin de poder seguir todo el proceso. Los datos deben lucir de la siguiente manera:

En la imagen se muestran los primeros 6 datos de la base que utilizaremos en este ejemplo y que tiene el nombre de cd. La columna Y representa el valor de producción, K el valor del factor Capital, y finalmente L es el valor del factor Trabajo; todos estos valores son en realidad, el logaritmo neperiano del verdadero valor monetario ($), efectuado así para representar linealmente la función.

Lo siguiente que debemos hacer es estimar los coeficientes (elasticidades) de los factores de producción con el siguiente código:

cd.lm <-     lm(formula= I(Y) ~ I(K) + I(L),

             data  = cd, na.action=na.exclude );summary( cd.lm )

Call:

lm(formula = I(Y) ~ I(K) + I(L), data = cd, na.action = na.exclude)

Residuals:

    Min      1Q  Median      3Q     Max 

-0.2812 -0.0938 -0.0030  0.0904  0.3340 

Coefficients:

            Estimate Std. Error t value             Pr(>|t|)    

(Intercept)   2.7235     0.2644   10.30 < 0.0000000000000002 ***

I(K)          0.4566     0.0632    7.23       0.000000000067 ***

I(L)          0.8845     0.0577   15.32 < 0.0000000000000002 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.131 on 111 degrees of freedom

Multiple R-squared:  0.817,     Adjusted R-squared:  0.814 

F-statistic:  248 on 2 and 111 DF,  p-value: <0.0000000000000002

A partir de este punto, es que comenzamos a plotear la función de producción, si sigues paso a paso el siguiente código, podrás identificar que parte del gráfico es la que se esta desarrollando.

getOption("device") ()                    ## Abre ambiente gráfico
interval <- 1
x <-   seq(from = floor(min(cd$K)/interval)*interval,
           to   = ceiling(max(cd$K)/interval)*interval,
           by   = interval)

y <-   seq(from = floor(min(cd$L)/interval)*interval,
           to   = ceiling(max(cd$L)/interval)*interval,
           by   = interval)

f <-function(x,y){
   coef(cd.lm)[1] + x*coef(cd.lm)[2] + y*coef(cd.lm)[3]
   }

z <- outer(x,y,f)

zlim <- c(floor(min(z)), ceiling(max(z)))
#Alternativa zlim <- c(floor(min(z)), 1)

cd$Y.hat <- f(x = cd$K, y = cd$L)

## Comienza el plot 3D con argumentos específicos

cd.persp <- persp(
          x        = x,
          y        = y,
          z        = z,
          theta    = 310,
          phi      = 0,
          col      = 'cornflowerblue',
          ltheta   = -135,
          shade    = 0.01,
          ticktype = 'detailed',
          main     = 'Función de producción',
          xlab     = 'Capital',
          ylab     = 'Trabajo',
          zlab     = 'Valor producción',
          zlim     = zlim,
          scale    = TRUE,
          border   = NA,
          nticks   = 4 ,
          font.lab = 2)

# Se plotean los valores reales

cd.below <- subset(cd, cd$Y <= cd$Y.hat)

cd.below.trans3d <- trans3d(
       x    = cd.below$K,
       y    = cd.below$L,
       z    = cd.below$Y,
       pmat = cd.persp )

points( cd.below.trans3d, col = "red", pch = 16)

cd.below.hat.trans3d <- trans3d(
      x    = cd.below$K,
      y    = cd.below$L,
      z    = cd.below$Y.hat,
      pmat = cd.persp )

segments(
   x0  = cd.below.trans3d$x,
   y0  = cd.below.trans3d$y,
   x1  = cd.below.hat.trans3d$x,
   y1  = cd.below.hat.trans3d$y,
   lty = "solid",
   col = "red")

cd.above <- subset (cd, cd$Y >= cd$Y.hat)

cd.above.trans3d <- trans3d(
       x    = cd.above$K,
       y    = cd.above$L,
       z    = cd.above$Y,
       pmat = cd.persp )

points (cd.above.trans3d, col = "green", pch = 16)

cd.above.hat.trans3d <- trans3d(
        x    = cd.above$K,
        y    = cd.above$L,
        z    = cd.above$Y.hat,
        pmat = cd.persp )

segments(
     x0  = cd.above.trans3d$x,
     y0  = cd.above.trans3d$y,
     x1  = cd.above.hat.trans3d$x,
     y1  = cd.above.hat.trans3d$y,
     lty = "solid",
     col = "green")

# Agregamos la leyenda

legend("bottom",
       c("Ajustado","Real (superior)","Real (inferior)"),
       horiz=TRUE,
       pch=c(18,19,19),
       col=c("cornflowerblue","green","red"),
       inset=0.95,
       pt.cex=c(2,1.1,1.1),
       bty="n")

El resultado será (dependiendo de tus datos), un gráfico como el que se presenta al inicio. Sin embargo, tenemos otra vía que requiere muchas menos líneas de código; aunque cabe mencionar que existen algunas diferencias en cuanto a estética se refiere. Para elaborarlo requerimos de un paquete denominado scatterplot3d. Una vez instalada esta paquetería, procedemos con el siguiente código:

library(scatterplot3d)

cd.lm <-     lm(formula= I(Y) ~ I(K) + I(L),

             data  = cd, na.action=na.exclude );summary( cd.lm )

s3d <- scatterplot3d(cd$K,cd$L,cd$Y,
       pch=16,highlight.3d=TRUE,
       type="h",main="Función de producción")
s3d$plane3d(cd.lm)


Lo que produce un gráfico como el siguiente:

Claro está que existen más argumentos para esta función (scatterplot3d), lo que podría producir un gráfico de mejor calidad, aunque corresponde al usuario indagar para lograr el objetivo que busca.