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.

Plotear un corazón 3D en R

En esta ocasión, después de no publicar en un largo tiempo, me doy a la tarea escribir algo con menor carga académica, pero que quizás en determinado momento, ya sea por curiosidad o por alguna necesidad tengas el interés en replicar e inclusive mejorar.

Plotear un corazón (de San Valentín) en dos o tres dimensiones en algún Software no es algo novedoso, sin duda que no lo es, existen numerosos códigos en Matlab, Mathematica, C, etc., regados por Internet. Pero para los que utilizamos mayormente R, aquí esta la alternativa.

El único paquete que necesitamos es misc3d, y plotear la siguiente función:

En donde, tanto x, y como z, son vectores que van desde -3 a 3, con intervalos k, una vez explicado esto, lo único que debemos hacer escribir las siguientes lineas de código:

# install.packages('misc3d')
library(misc3d)

k <- 0.030
x <- seq(-3, 3, by = k);x
y <- seq(-3, 3, by = k);y
z <- seq(-3, 3, by = k);z

f  <- function(x ,y, z)(-(x^2)*(z^3)-(9/80)*(y^2)*(z^3)) + ((x^2)+(9/4)*(y^2)+(z^2)-1)^3
contour3d(f,0,x,y,z, color='red', smooth=TRUE, fill= TRUE)

Existen algunos argumentos de la función que bien valdría la pena explorar, todo acorde a tus necesidades. También es valioso comentar que el intervalo k, entre más pequeño sea, mejor apariencia tendrá nuestro gráfico, pero el tiempo que tomará R para plotearlo será mayor.