# Aplicación: Estimación de una Función de Producción


print(paste("Lo  siguiente es ejecutado el", date())) #estampa fecha y hora

# Lee datos de producción disponibles en R

library(Ecdat)
help(Metal)               # entrega información de la base de datos que se utilizará.
data(Metal)               # coloca a los datos en memora.
names(Metal)              # entrega los nombres de las columnas
attach(Metal)             # para evitar colocar $ cuando se llame a una variable
summary(Metal)            # estadística descrptiva de los datos. Observe
                          # que en Max existen datos que están muy lejos de la media
                          # podrían haber problemas a futuro
met=as.matrix(Metal)      # crea matriz de datos
Ly=log(met[,1])           # 1a columna de la matriz es logaritmo de y, o sea, Ly.
LL=log(met[,2])           
LK=log(met[,3])
datmtx=cbind(Ly,LK, LL)   #las une en una sola matriz, datmtx
fBasics::basicStats(datmtx) #entrega estadísticos
matplot(cbind(Ly, LK, LL), typ="b", pch=c("y","K","L"),
        main="Descripción de la data Metales", xlab="Número de la observación",
        ylab="Logaritmos del producto, capital, and trabajo") #grafica logaritmos

# Estimación Lineal, en logaritmos
met=as.matrix(datmtx) #met es la matriz de datos
reg1 = lm(Ly~LK+LL) #la regresión del log del producto en función del log del capital y del log del trabajo
summary(reg1) # todas significativas y positivas. R2 es 94%.
names(reg1) #variables de la regresión disponibles de trabajar.

# Estimación No lineal
#y=met[,1]; L=met[,2]; K=met[,3]
#regnon=nls(y~A*(K^alp)*(L^bet), # las variables no van en logaritmos
#           start=list(A=1.17064,alp=0.37571,bet=0.60300)) #A,alp y bet puse los de la estimacón lineal para comenzar a iterar
#summary(regnon) # todas significativas y positivas.

# Veamos estadísticos de regresión 1, reg1, regresión lineal
fitted.values(reg1)          # valores estimados
summary(fitted.values(reg1)) # descripción de estadísticos relevantes
hist(fitted.values(reg1))    # histograma de valores estimados
residuals(reg1)              # residuos
summary(residuals(reg1))     # estadísticos de los residuos
hist(residuals(reg1))        # histograma de residuos

# Gráfico de diagnóstico para la regresión 1, reg1, en logaritmos
# veremos estadísticos para outliers
influence.measures(reg1)  #observaciones 1,16 y 22, * es significancia estadística y posible outlier
plot(reg1,1) # 12, 22 y 26 
plot(reg1,2) # 12, 22 y 26
plot(reg1,3) # 12, 22 y 26
plot(reg1,4) # 12, 13 y 22
plot(reg1,5) # 12, 13 y 22 
plot(reg1,6) # 12, 13 y 22 # Cook expone sensibilidad de la línea de regresión a 
#todas las observaciones individuales
#el residuo estandarizado para la observación 22 esta fuera de los límites de la
#línea segmentada que marca número 0.5. El eje horizontal tiene un leverage de cada
#punto medido por el tamaño de su distancia de Cook.Si el hecho de insertar o 
#sacar el punto cambia mucho la regresión, tiene una gran importancia, ya que 
#modifica muchísimo la gráfica y se genera una especie de outlier.

#Primero vamos a ajustar el modelo y luego vamos a aplicar la prueba de Bonferroni para detectar outliers.
library(car)
outlierTest(reg1, cutoff=Inf, n.max=4) #obs 22 y 26 son outliers
influenceIndexPlot(reg1, vars="Bonf", las=2) #que sean outliers no implica
#necesariamente que deban ser eliminadas del modelo, quizás 
#no afecten la estimación de los parámetros.

#Segundo, veremos distancia de Cook. Es una medida de cómo influye la observación  
#i-ésima sobre la estimación de β al ser retirada del conjunto de datos
#Una distancia de Cook grande significa que una observación tiene un peso 
#grande en la estimación de β.
cooks.distance(reg1)
#Es mejor representar las distancias de Cook en forma gráfica
#para identificar los posible puntos influyentes así:
influenceIndexPlot(reg1, vars="Cook") # vea observaciones 12,13 y 22
#Ahora vamos a revisar los residuales del modelo
par(mfrow=c(2, 2))
plot(reg1, col='deepskyblue4', pch=19)  # vea observaciones 12,13 y 22
#tienen valores residuales grandes.
#Si hay una razón de peso para considerarlas como observaciones atípicas, 
#ellas deben salir del modelo. Si por el contrario, no hay nada raro con 
#las observaciones, ellas deben seguir en el modelo.
hv=hatvalues(reg1) #valores estimados de hv
hv
plot(hv, typ="l")
cd=cooks.distance(reg1)
cd
plot(cd, typ="l") #12, 13 y 22 se observan en datos y gráfica como outliers

#influyentes, según distancia de Cook, deben ser elminados.
Metal1=Metal[-c(12,13,22),] #elimina las observaciones
#con la data Metal1 se realiza la nueva regresión.


# Testeo de retornos constantes a escala. Trabajaremos con reg1
library(car)
linearHypothesis(reg1, "1*LK+1*LL=1")
# p-value 0.7366 > 0.05  acepta H0, hay retornos constantes a escala en base Metal.

#Gráficos 3D
library(scatterplot3d)
length(labor) #largo de la serie 27
mynumlet=c(as.character(1:9),letters[1:18]) #una serie de 9 números y 18 letras
print(mynumlet)
Met=data.frame(capital,labor,va) #va es valor agregado del producto, 
                                 #el input tiene que ser data.frame para scatterplot3d
m3d <- scatterplot3d(Met, type="h", highlight.3d=TRUE,
                     angle=55, scale.y=0.7, pch=mynumlet,
                     main="Scatterplot de la Producción de Metales (unidades originales)",
                     box=F,zlab="Valor agregado del producto")
my.lm <- lm(va ~ capital+ labor) 
summary(my.lm)  #Pero.....¡¡esta no es Cobb-Douglas!!, hay que linealizar, "usar" con logaritmos

#Uso de logaritmos (regresión linealizada)
Metlog=data.frame(LK, LL, Ly)
m3dL<- scatterplot3d(Metlog, type="h", highlight.3d=TRUE,
                     angle=55, scale.y=0.7, pch=mynumlet,
                     main="Metales: Función de Producción Cobb Douglass",
                     box=F, zlab="Log va (valor agregado del producto)", xlab="Logaritmo del capital",
                     ylab="Log del trabajo")
my.lm2 <- lm(Ly ~ LK+LL)
m3d$plane3d(my.lm2, lty.box = "solid")
summary(my.lm2) # signos correctos, individualmente significativas las variables y R2=94%

##################################################################################
# DE AQUÍ EN ADELANTE OCUPAREMOS OTRA BASE DE DATOS, ver archivo excel adjunto

# Leeremos datos y correremos 3 función de producción distintas
rm(list=ls()) #remueve objetos, limpia la memoria (Console y Environment)

url="https://faculty.fordham.edu/vinod/WECOData.csv" #nueva data de producción
weco=read.csv(url)
names(weco)            # lista los nombres, hay 6 variables
library(fBasics)       # llama al paquete de estadísticas descriptivas
basicStats(weco)       # calcula estadísticos
attach(weco)           # permite acceso a variables sin poner $
length(Y)              # largo de la serie, 59

#Tome el logaritmo natural de la variable del 
#producto cruzado de los tres factores.
lY=log(Y)
lK=log(K)
G=Engg
lG=log(G)
lL=log(L) #log de producto, capital, ingeniería y trabajo.
lKL=lK*lL #interacciones y forma cuadrática
lKG=lK*lG
lLG=lL*lG
lK2=lK^2# Los siguientes tres regresores son necesitados para la 
        # regresión de la función de producción translogaritmica.
lG2=lG^2
lL2=lL^2
reg=lm(lY~lK+lL+lG+lKL+lKG+lLG) #regresión sin términos al cuadrado
                                #homogenea no multiplicativa, es una
                                #versión restringida de Cobb-Douglass
names(reg)
summary(reg)#prints results
a0=reg$coefficients[1]#extrae intercepto
a1=reg$coefficients[2]#idem para log capital
a2=reg$coefficients[3]#idem para log trabajo
a3=reg$coefficients[4]#idem para ingeniería
a4=reg$coefficients[5]#idem para interacción (log K)(log L)
a5=reg$coefficients[6]#idem para interacción (log K)(log G)
a6=reg$coefficients[7]#idem para interaction (log G)(log L)

# los siguientes son vectores, no constantes.
MEK=a1+a4*lL+a5*lG  # elasticidad marginal del capital
MEL=a2+a4*lK+a6*lG  # elasticidad marginal del trabajo
MEG=a3+a5*lK+a6*lL  # elasticidad marginal de la ingeniería
SCE=(MEK+MEL+MEG)   # elasticidad de escala
inv.SCE=(MEK+MEL+MEG)^(-1) #inverso de la elasticidad de escala

#los siguientes comandos generan los gráficos
plot(MEK, type="l", xlab="mes",
     main="WECo: Elasticidad marginal del Capital")
plot(MEL, type="l", xlab="mes",
     main="WECo: Elasticidad marginal del Trabajo")
plot(MEG, type="l", xlab="mes",
     main="WECo: Elasticidad marginal de la Ingeniería")
plot(SCE/10, type="l", xlab="mes",
     main="WECo: Elasticidad de Escala =AC/MC",lwd=2)
plot(inv.SCE*10, type="l", xlab="mes", main="WECo: Cambio
en el costo total cuando el producto aumenta en 1%", ylab="MC/AC",lwd=2)
mrtsLK=(MEL/MEK)*(K/L)# tasa marginal de sustitución técnica entre trabajo L y capital K
cor(G, mrtsLK)# si la correlación es baja, entonces el input ingeniería es 
              # separable del capital y trabajo a la vez
              #corr=0.02236146. Por lo tanto, si es separable.

#########################################################################
# Función de R para definir contornos de la función de producción
# Aquí se estiman tres funciones: Cobb-Douglas, TransLogaritmica y Multiplicativa no homogenea
########### Primero correr la función que dibujará.######################

pfcontour=function(y,K, L, level=T, z=0, type=c("Cobb-Douglas",
                                                "MNH","TransLog"),n50=50) 
  {
  # ajusta la regresión y dibuja contornos usando paquete clines
  Ly=y; LK=K; LL=L                # vuelve a calcular los logaritmos
  if (level) Ly=log(y)            # esto lo podríamos haber sacado de este programa.
  T=length(y)
  if(level) LK=log(K);
  LK2=LK^2
  if (level) LL=log(L);
  LL2=LL^2; LKLL=LK*LL
  #print(c(T,length(z)))
  if(length(z)==T){               # z es el periodo
    reg1=switch(type,
                "Cobb-Douglas" = lm(Ly~LK+LL+z),
                "MNH" = lm(Ly~LK+LL+LKLL+z),
                "TransLog" = lm(Ly~LK+LL+LKLL+LK2+LL2+z))    
    print(reg1)
  }#end if length(z)
  if(length(z)==1){
    reg1=switch(type,
                "Cobb-Douglas" = lm(Ly~LK+LL),
                "MNH" = lm(Ly~LK+LL+LKLL),
                "TransLog" = lm(Ly~LK+LL+LKLL+LK2+LL2))
    print(reg1)
  }#end if length(z)
  
  Lymtx=matrix(NA,n50,n50)
  a=as.numeric(reg1$coefficients)
  x=rep(NA,n50)
  y=rep(NA,n50)
  rangeLK=(max(LK)-min(LK))/n50
  rangeLL=(max(LL)-min(LL))/ n50
  for (i in 1:n50){Lk=min(LK)+rangeLK*i;x[i]=Lk
  for (j in 1:n50){Ll=min(LL)+rangeLL*j;y[j]=Ll
  #
  if (length(z)==T){
    Lymtx[i,j]= switch(type,
                       "Cobb-Douglas" =a[1]+a[2]*Lk+a[3]*Ll+a[4]*z,
                       "MNH"=a[1]+a[2]*Lk+a[3]*Ll+a[4]*Lk*Ll+a[5]*z,
                       "TransLog" =a[1]+a[2]*Lk+a[3]*Ll+a[4]*Lk*Ll+a[5]*Lk^2+
                         a[6]*Ll^2+a[7]*z
    ) #end switch simple parenthesis
  }# end if length z
  #
  if (length(z)==1){
    Lymtx[i,j]= switch(type,
                       "Cobb-Douglas" =a[1]+a[2]*Lk+a[3]*Ll,
                       "MNH"=a[1]+a[2]*Lk+a[3]*Ll+a[4]*Lk*Ll,
                       "TransLog" =a[1]+a[2]*Lk+a[3]*Ll+a[4]*Lk*Ll+a[5]*Lk^2+a[6]*
                         Ll^2
    ) #end switch simple parenthesis
  }# end if length z
  }#end j loop
  }#end i loop
  #print(head(Lymtx))
  contour(exp(x),exp(y),exp(Lymtx),
          main=paste(c("Curvas de nivel de la función de producción",type),
                     sep=" "), xlab="capital", ylab="trabajo",lwd=2)
  #
  list(Lymtx=Lymtx, reg1=reg1)
  #estas listas son el producto de la función.
  #In Lymtx=Lymtx the Lymtx la igualdad de la izquierda es el nombre de salida publicado fuera de la función  y el mismo
  # Lymtx en la igualdad del lado derecho fuera de la función.
}  #FIN de la función

######################################################################################
# Aquí trabajaremos otra información y con la función de dibujo gráfico

#objects() # este comandos lista todos los objetos en memoria
#rm(list=ls()) #remueve todo en Consola y Enviroment

url="https://faculty.fordham.edu/vinod/belldata.csv"
bell=read.csv(url)
#bell=read.table(file="c:/data/belldata.csv", header=T,sep=",") #no lo usaremos

summary(bell)
attach(bell)
pfc= pfcontour(y,k,lab,level=T, type="MNH",n50=50,z=0) #pfcontour es la función anterior
pfc2= pfcontour(y,k,lab,level=T, type="Cobb-Douglas",n50=50,z=0)
pfc3= pfcontour(y,k,lab,level=T, type="TransLog",n50=50,z=0)
#arriba la línea crea un objeto llamado pfc que contiene el producto
#de la función pfcontour. De esta manera, R no imprimirá
#entero Lymtx y reg1 y desordenar la pantalla
#el menu de archivo de R permite grabar el gráfico en diversos formatos.

#Testeo de la estimación de la función de producción
#rm(list=ls()) #rm remueve objetos
#clean slate
mylink=("https://faculty.fordham.edu/vinod/belldata.csv")
#bell=read.table(file=mylink, header=T,sep=",")
summary(bell); attach(bell)
Ly=log(y); Lx1=log(k); Lx2=log(lab); Lx6=log(poiss6) 
Lx12=(Lx1*Lx2);Lx11=(Lx1^2); Lx22=(Lx2^2)           

regnh=lm(Ly~Lx1+Lx2+Lx12) #Función de producción multiplicativa no homogenea, no lleva al cuadrado ni Lx6
library(car)#hace los test de hipótesis
linearHypothesis(regnh, "1*Lx12 = 0") #p-val=8.585e-11 Reject. Interacción capital trabajo significativo

regb=lm(Ly~Lx1+Lx2+Lx12 +Lx11 +Lx22 +Lx6)  # Función de producción translogarítmica
summary(regb) # sólo sale la tecnología estadísticamente significativa

regbi=lm(Ly~Lx1+Lx2+Lx1:Lx2 +Lx6);summary(regbi) # los ":" significa que los 
#parámetros de las interacciones (multiplicaciones de variables) y elevaciones 
#al cuadrado son todas cero a la vez. #se rechaza H0. Son significativas.

bcross=regbi$co[5] #coeficientes de todos los términos de interacción
#test nulo, todos los coeficientes de interacción se igualan a cero  0.4729
coefs <- names(coef(regbi))
linearHypothesis(regbi, coefs[grep(":", coefs)],verbose=TRUE)
# the p-value for above test =0.0005175
#Por lo tanto, rechazamos la hipotesis nula que todas las interacciones son =0

# Cálculo de las elasticidades para Función de producción translogarítmica
a=regb$coe[2:7]
ME1=a[1]+a[3]* mean(Lx2) +2*a[4]* mean(Lx1) #elasticidad marginal para X1
ME2=a[2]+a[3]* mean(Lx1) +2*a[5]* mean(Lx2) #elasticidad marginal para x2
SCE=ME1+ME2;SCE #elasticidad de escala=1.341316 >>1, hay monopolio natural
EOS=SCE/(SCE+2*bcross);EOS #elasticidad de sustitución o EOS evaluado a la media=0.739

# Diagnóstico de colinealidad para  Bell System, translogarítmica
X=cbind(Lx1,Lx2,Lx12,Lx11,Lx22,Lx6)
# Cómputo de svd
svdx=svd(X)
cond.no=svdx$d[1]/svdx$d[ncol(X)]
cond.no # 115562.6 mucho mas grande que 10p=70, Por lo tanto, colinealidad

# Nueva función en R para profunda comprensión del rank de deficiencia
# Función obtiene m (rank deficiency) de k (ridge biasing parameter).  
## Comienza la nueva función en la siguiente línea
getmfromk=function(ei, maxk=2, showplot=T, n100=100){
  #Input       ei= eigenvalues
  #Input      maxk= mayores k permitidos
  p=length(ei)
  k=seq(0,maxk, maxk/n100)
  m=rep(0,length(k))
  del=rep(0,length(ei))
  for (i in 1:(n100+1)){
    for (j in 1: length(ei)){
      del[j]=ei[j]/(ei[j]+k[i])}
    sumd=sum(del)
    m[i]=p-sumd
    i=i+1}
  if (showplot){
    plot(k,m, type="l", main="Parámetro de polarización de cresta versus rango
deficiencia")}
  list(m=m,k=k)
}  #Fin de la función llamada getkfromm (get k de m)
# ahora correr
getmfromk(svdx$d^2,maxk=0.1)

# Regresión Ridgen del paquete MASS & RXshrink de R
#li=("https://faculty.fordham.edu/vinod/belldata.csv")
#bell=read.table(file=li, header=T,sep=",")
#attach(bell)
#Ly=log(y); Lx1=log(k); Lx2=log(lab); Lx6=log(poiss6)
#Lx12=(Lx1*Lx2);Lx11=(Lx1^2); Lx22=(Lx2^2)
formu=Ly~Lx1+Lx2+Lx12 +Lx11 +Lx22 +Lx6
library(MASS)#lm.ridge
select(lm.ridge(Ly~Lx1+Lx2+Lx12 +Lx11 +Lx22 +Lx6, lambda =
                  seq(0,0.01,0.0025)))
HKB=lm.ridge(formu,lambda=4.467141e-05) #obtenido comando anterior
LW=lm.ridge(formu,lambda=0.007425944)   #obtenido comando anterior
GCV=lm.ridge(formu,lambda=0.0025)       #obtenido comando anterior
library(RXshrink)                       # para dos parámetros ridge k y Q
mydf=data.frame(cbind(Ly,Lx1,Lx2,Lx12 ,Lx11 ,Lx22 ,Lx6))
#el comando de arriba es necesario para obtener el paquete a trabajar
qmr=qm.ridge(formu, mydf);qmr;plot(qmr) #no comentado por brevedad
mc=MLcalc(formu,data=mydf,rscale=2)
Ob=mc$beta[2,]#coeficientes relevantes a través de segunda fila
ME1=rep(NA,4); ME2=ME1# guarda elasticidades

for (j in 1:4){
  if(j==1) a=coef(HKB)[2:7]
  if(j==2) a=coef(LW)[2:7]
  if(j==3) a=coef(GCV)[2:7]
  if(j==4) a=Ob
  ME1[j]=a[1]+a[3]* mean(Lx2) +2*a[4]* mean(Lx1)
  ME2[j]=a[2]+a[3]* mean(Lx1) +2*a[5]* mean(Lx2)
}#end of j loop
mtx=cbind(ME1,ME2)#preparando a imprimir los resultados
rownames(mtx)=c("HKB","LW","GCV","RXshrink")
print(mtx)

#Función para estandarizar X'X
stdze = function(oldx){
  # estandariza oldx= Matriz Original
  bigt=nrow(oldx)
  sqb=sqrt(bigt-1)
  m=ncol(oldx)
  cmn=apply(oldx,2,mean)#column means c=column
  csd=apply(oldx,2,sd)#column standard deviations
  newx=oldx#---Initialize----#
  for (j in 1:m){
    newx[,j]=(oldx[,j]-cmn[j])/(sqb*csd[j])
  } #end loop for j
  return(newx)
}

#Function para desestandarizar datos y coeficientes

unstdze = function(y, oldx, b) {
  #  oldx=Original Matrix of y and regressors without column of
  #ones, INPUT: y= dependent variable data and oldx is X matrix
  #INPUT: b=regr coefficients based on standardized data
  bigt=nrow(oldx)
  sqb=sqrt(bigt-1)
  m=ncol(oldx)
  if (m != length(b) ) print("Error in unstdze function ncol
#not=length(b)")
  cmn=apply(oldx,2,mean)#column means c=column
  csd=apply(oldx,2,sd)#column standard deviations
  unstdb=b #---Initialize----#
  for (j in 1:m){
    unstdb[j]=(b[j])/(sqb*csd[j])
  } #end loop for j
  for (j in 1:m){
    unstdb[j]=(b[j])/(sqb*csd[j])
  } #end loop for j
  intercept=mean(y)-sum(cmn*unstdb)
  list(intercept=intercept, unstdb=unstdb) }
# fin de la función

# Función ridge regrresion con estandarización
ridge.std=function(y,x,k=c(0,.05,0.01), bestk=0, plot=T)
{  #estimate the ridge regression model regressing y on X
  #version standardizes the data
  ys=y-mean(y)
  xs=stdze(x)
  xx= t(xs) %*%xs
  mink=min(k)
  p=ncol(xs)
  if (length(k)>1) {
    mymtx=matrix(NA,length(k),p)
    for (i in 1:length(k)){
      ki=k[i]*diag(ncol(xs))
      xy= t(xs) %*%ys
      sinv=solve(xx+ki)
      bk=sinv %*%xy
      mymtx[i,]=bk
      if (k[i]==0) {lsfit=xs%*%bk; lscoef=bk}
      resid=ys-xs%*%bk
      s2=sum(resid^2)/(length(ys)-ncol(xs)-1)
      if(mink==0){
        if(k[i]==0){
          HKB <- (p - 2) * s2/sum(lscoef^2); print(c("HoerlKennard
#Baldwin
k=",HKB),q=F)
          LW <- (p - 2) * s2 * length(ys)/sum(lsfit^2);print(c("Lawless
Wang k=",
                                                               LW),q=F)}}
    }# end k loop
    if(plot){
      mynum=as.character(1:p)
      nam=1:p
      matplot(k,mymtx,main="Ridge Trace: Look for stable region",
              type="b",pch=mynum,
              xlab="Ridge biasing parameter k", ylab="Ridge regression
coefficients", lty=1:6)
    } }#end if length(k)>1 end if for plot=T
  print(c("value of ridge biasing parameter chosen is =",bestk),
        q=F)
  ki=bestk*diag(ncol(xs))
  sinv=solve(xx+ki)
  xy= t(xs) %*%ys
  bk=sinv %*%xy
  resid=ys-xs%*%bk
  s2=sum(resid^2)/(length(ys)-ncol(xs)-1)
  varbk=s2*(sinv %*%xx %*% sinv)
  se.bk=sqrt(diag(varbk))
  coef=as.numeric(bk)
  unst1=unstdze(y,x,coef)
  unst=as.numeric(unst1$unstdb)
  unst.int=as.numeric(unst1$intercept)
  tstat=coef/se.bk
  out=cbind(unst,coef,se.bk,tstat)
  print(out)
  list(bk=bk,varbk=varbk,se.bk=se.bk,unst=unst,
       intercept=unst.int)}
#example
ridge.std(y,X,k=0,bestk=0.1)

#Regresión Ridge estandarizada para Bell Data
X=cbind(Lx1,Lx2,Lx12 ,Lx11 ,Lx22 ,Lx6)
rr=ridge.std(Ly,X,k=seq(0,.1,.01),bestk=.025)
a=rr$unst  #coeficientes de regresión usar desestandarizado
ME1=a[1]+a[3]* mean(Lx2) +2*a[4]* mean(Lx1);ME1
#ME1=0.5966401 is ahora sensible
ME2=a[2]+a[3]* mean(Lx1) +2*a[5]* mean(Lx2);ME2
# ME2=0.8386455 es ahora sensible
print(c("Elasticidad marginal para datos de Bell en la media",ME1,ME2),q=F)
SCE=ME1+ME2
print(c("Elasticidad de escala en la media SCE =",SCE),q=F)
# SCE= 1.43529 >1 #hay economías de escala