Fronteras estocásticas de producción

1. Introducción

2. Frontera estocástica de producción

3. Método del área común

4. Detección y eliminación de valores anómalos

5. Conclusiones

6. Librería stfa en R

7. Futuros pasos

El análisis de frontera estocástica es uno de los métodos más utilizados para analizar y medir la eficiencia.

Surgen en 1977, como una sofisticación de los modelos más tradicionales de frontera de producción.

Son ampliamente utilizados en los campos de la economía y la cienciometría.

El modelo de frontera de producción estocástica se puede expresar como: \[ y_{i} = f(\vec x_{i},\vec \beta) \cdot e^{v_{i}} \cdot E_{i} \] donde

• \( y_{i} \) y \( ~\vec x_{i} \) es el valor de los inputs y el output para el productor \( i \).

• La función de producción es una Cobb-Douglas cuya expresión \( f(\vec x_{i},\vec \beta)=\prod_{k=1}^{K} x_{ki}^{\beta_{k}} \).

• \( \vec \beta \) es el vector con los parámetros del modelo que deseamos estimar.

• \( v_i \) representa la componente aleatoria para el productor i-ésimo y \( ~E_{i} \) representa el término de la eficiencia técnica para el productor i-ésimo.

Sea \( ~E_i = e^{-u_i} \), teniendo en cuenta las expresiones anteriores y tomando logaritmos, podemos expresar el modelo con la siguiente expresión:

\[ \ln{y_i} = \overbrace{\beta_0 + \sum_n \beta_n \ln{x_{ni}}+v_i }^{\mbox{Frontera Estocástica}} -u_i \] \[ \ln{y_i} = \beta_0 + \sum_n \beta_n \ln{x_{ni}}+\underbrace {v_i -u_i}_{\epsilon_i} \]

donde \( v_i \equiv \mathfrak{N}(0,\sigma^2_v) \) y \( u_i \equiv | \mathfrak{N}(0,\sigma^2_u) | \)

La función de densidad de \( \epsilon \) a partir de los supuestos anteriores es:

\[ f(\epsilon)=\int_0^\infty f(\epsilon,u)=\frac{2}{\sigma} \cdot \phi\left(\frac{\epsilon}{\sigma}\right) \cdot \Phi\left(-\frac{\epsilon\lambda}{\sigma}\right) \]

donde

• \( \sigma = \sqrt{\sigma_{u}^{2}+\sigma_{v}^{2}} \)

• \( \lambda = \frac{\sigma_u}{\sigma_v} \), nos ofrece una estimación del peso de cada componente del error.

• \( \phi(\cdot) \) y \( \Phi(\cdot) \) son la función de densidad y distribución de una normal típica.

Bajo los supuestos del modelo, el coeficiente de asimetría de Fisher ha de ser negativo:

\[ \mathcal{A}=\frac{\sqrt{2}(\pi-4)\lambda^3}{(\pi+(\pi-2)\lambda^2)^{3/2}} < 0 \]

Este problema está presente cuando las muestras son pequeñas o \( \lambda \) es pequeño. A medida que va aumentando el tamaño de muestra se va corrigiendo lentamente el problema.

El logaritmo de la función de verosimilitud de \( \epsilon \) es:

\[ \ln{L}=-\frac{N}{2}\cdot \ln{\frac{\pi}{2}} -N\cdot \ln{\sigma}+\sum_{i=1}^{N}\ln{\Phi\left(-\frac{\epsilon_i\lambda}{\sigma}\right)}-\frac{1}{2\sigma^2}\cdot \sum_{i=1}^{N}\epsilon^2_i \]

Derivando respecto de los parámetros de interés obtendríamos donde se hace máxima la verosimilitud.

Para realizar las estimaciones existen distintas bibliotecas en R donde está ya implementado este método. En este caso se utilizará la biblioteca frontier.

Para simular los datos se propone un modelo de un único input con las siguientes especificaciones:

\[ y_i=100+5x_i+\epsilon_i ~~ \mbox{donde} ~~i=1,...,K~~~;~~~\sigma_v^2=1 \]

Utilizando el método MV, la idea es introducirle una modificación que ayude a tener en cuenta la rugosidad de la curva a estimar.

\[ \ln L(X_i|\psi)=\sum_{i}\log{f(X_i)}-\psi R(f) \]

donde \( \psi \) es un parámetro de suavización positivo y \( R(f) \) es una función que recoge el nivel de rugosidad de \( f \).

\( R(f) \) es la raíz cuadrada de la función de densidad resultando la siguiente función penalizadora. Sea \( \rho=\sqrt{f} \) donde \( f \) es una función de densidad:

\[ R(f)=\int{\rho^{'2}} \]

En este trabajo se ha explorado el comportamiento de este método con \( 0<\psi<1 \), sin embargo los mejores resultados se obtenian cuando \( \psi=0 \).

La estimación núcleo se define como: \[ \hat{f}(x)=\frac{1}{Nh}\sum_{i=1}^N \kappa\left(\frac{x-x_i}{h}\right) \]

donde

• \( \kappa \) es una función núcleo, generalmente una función de densidad unimodal y simétrica respecto a 0, y que cumple determinadas propiedades.
• \( h \) es una constante positiva que se conoce como ancho de ventana o parámetro de suavización.

En nuestro caso:

• \[ \kappa(x)=\frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}} \]
• \( h=0.9AN^{(-1/5)} \) donde \( A=\mbox{mín}(\sigma,IQ/1.349) \)

Podemos obtener una expresión de la función de densidad de los residuos en función de \( \lambda \) y la desviación típica de los residuos:

\[ f( \epsilon)=\frac{2\sqrt{\lambda^2(\pi-2)+\pi}}{\sigma_\epsilon \sqrt{\pi(\lambda^2+1)}}\cdot\phi\left(\frac{ \epsilon\sqrt{\lambda^2(\pi-2)+\pi}}{\sigma_\epsilon\sqrt{\pi(\lambda^2+1)}}\right)\cdot\Phi\left(\frac{- \epsilon\lambda\sqrt{\lambda^2(\pi-2)+\pi}}{\sigma_\epsilon\sqrt{\pi(\lambda^2+1)}}\right) \] donde \( \sigma_\epsilon=\sigma\sqrt{\frac{\lambda^2(\pi-2)+\pi}{\lambda^2 \pi+\pi}} \) es la desviación típica de los residuos.

En el caso de \( \sigma_\epsilon \) podemos utilizar su estimación máximo-verosímil, ya que este método nos permite obtener un estimación adecuada de este parámetro.

Utilizando los resultados anteriores podemos definir el área común como:

\[ \mathcal{AC}(\hat{\lambda})=\mbox{máx}_{\lambda} \int_{\mathbb{R}}\mbox{mín}\left\{\hat{f}(\vec \epsilon),f(\vec \epsilon~|~\lambda,\hat{\sigma}_\epsilon)\right\}d\vec \epsilon \]

Si recordamos, los residuos \( \epsilon \) presentaban una distribución asimétrica ya que el valor de la asimetría teórico era:

\[ \mathcal{A}=\frac{\sqrt{2}(\pi-4)\lambda^3}{(\pi+(\pi-2)\lambda^2)^{3/2}} \]

Esto es un problema ya que la funciones núcleo son por lo general simétricas, tal y como sucede en nuestro caso con la función núcleo gaussiana.

Sea \( X \) es una variable aleatoria cuya función de distribución, que denotamos por \( F \), es continua e invertible. Según el Teorema de inversión, la v.a. \( Y=F(X)\equiv \mathfrak{U}(0,1) \).

Con esta transformación tendríamos solucionado el problema, ya que la distribución \( \mathfrak{U}(0,1) \) es simétrica.

Esto es debido a que la distribución uniforme tiene una densidad discontinua, con saltos notables. Para corregir este problema, sobre la variable \( Y \) realizamos la siguiente transformación, \( Z=\Phi^{-1}(Y)~,~Z\equiv \mathfrak{N}(0,1) \). Esto nos permite obtener una transformación de los residuos que conduce a una densidad continua y derivable.

La metodología es es la siguiente:

1. Obtenemos una estimación de \( \sigma_\epsilon \) y de \( \vec \epsilon \) mediante el método de máxima verosimilitud.

2. Aplicamos la transformación \( \vec {\epsilon}_{T}=\Phi^{-1}(F(\vec{\epsilon})) \).

3. Obtenemos una estimación de \( \lambda \) mediante el método \( \mathcal{AC} \).

Para realizar el análisis utilizamos el mismo modelo que en la estimación MV para generar nuestras muestras:

\[ y=100+5x_i+\epsilon_i ~~ \mbox{donde} ~~i=1,...,K~~~;~~~\sigma_v^2=1 \]

Este tipo de funciones se conocen con el nombre de estimadores de núcleo adaptable y se pueden definir como:

\[ \hat{f}(x)=\frac{1}{N}\sum^N_{i=1}\frac{1}{ H_i}\kappa\left(\frac{x-x_i}{ H_i}\right) \] donde

• \( H_{i} \) es el ancho de ventana variable que puede ser expresado como \( H_{i}=h \tau_i \)
• \( h \) es el ancho de ventana constante
• \( \tau_i \) es el parámetro de ajuste.

1. Se realiza una estimación piloto \( \hat{f}^0 \) con una ancho de ventana \( h \) fijo.

2. Se calculan los valores de los anchos de ventana locales como

   \[ \tau_i = \left(\frac{\hat{f}^0(x_i)}{g}\right)^{-\alpha} \] donde \( g \) es la media geométrica de \( \{\hat{f}^0(x_i)\mid i=1,...,N\} \) y \( \alpha \) es un parámetro de sensibilidad, tal que \( ~0 \leq \alpha \leq 1~ \).

3. Se define la estimación mediante núcleos adaptables como

\[ \hat{f}(x)=\frac{1}{N}\sum^N_{i=1}\frac{1}{h \tau_i}\kappa\left(\frac{x-x_i}{h \tau_i}\right) \] donde \( \kappa \) será una función kernel, \( h \) será el ancho de ventana y será \( \alpha \) el parámetro que nos permita regular la sensibilidad (Véase la definición de \( \tau_i \)).

Para realizar la estimaciones se continúa utilizando el modelo planteado anteriormente. Además para calcular el área común tendremos en cuenta el valor del parámetro de ajuste \( \alpha \).

\[ \mathcal{AC_{(\lambda,\alpha)}}=\int_{\mathbb{R}} \mbox{mín}[\hat{f}(\vec \epsilon),f(\vec \epsilon~|~\lambda,\alpha,\sigma_\epsilon)]d \vec \epsilon \]

Problemas a la hora de evaluar la normalidad de los residuos:

• Inconsistencia de los test estadísticos tradicionales
• Medidas inadecuadas para tamaño muestral pequeño

Método basado en una medida conocida como neguentropía, utilizada en los campos de teoría de la señal y teoría de la información.

La entropía mide la cantidad incertidumbre de una señal.

Para una variable aleatoria discreta \( X \) se puede definir su entropía \( H(X) \) como:

\[ H(X)=-\sum_{i}{P(X=a_i)\log{P(X=a_i)}} \]

donde \( a_i \) son los posibles valores de \( X \).

El concepto de entropía se puede generalizar para variables aleatorias continuas. En este caso se conoce como entropía diferencial. Sea \( X \) una variable aleatoria continua con función de densidad \( f(x) \), se puede definir la entropía diferencial como :

\[ H(X)=-\int{f(x)\log{f(x)}dx} \]

De entre todas las distribuciones de probabilidad con soporte en \( (-\infty,\infty) \) y con una varianza \( \sigma^2 \) concreta, es la distribución normal la que tiene una mayor entropía. Es por ello que podemos utilizar la entropía diferencial como medida de no normalidad de una señal. A partir de esta expresión es posible obtener una medida de no normalidad que sea no negativa y que tome valor cero cuando la variable sea normal. Estamos hablando de la neguentropía.

Se puede entender como la diferencia entre la cantidad de entropía que tiene una señal y el máximo de entropía que puede llegar a tener. Es por tanto una medida en este caso de organización de un sistema o señal.

Sea \( X \) una variable aleatoria, se puede definir la neguentropía \( J(X) \) como :

\[ J(X)=H(X_{\mathfrak{N}})-H(X) \]

donde \( H(X_{\mathfrak{N}}) \) es la entropía de una variable aleatoria normal con la misma media y varianza que \( X \).

Dado que la estimación de la neguentropía es algo compleja existen distintas aproximaciones. La más sencilla se calcula en base a momentos de alto orden. Sea \( X\equiv \mathfrak{N}(0,1) \):

\[ J(X)\approx \frac{1}{12} E(X^3)^2 + \frac{1}{48} \mbox{curtosis}(X)^2 \]

Existen más opciones como : \[ J(X)\approx (E[G(X_{\mathfrak{N}})]-E[G(X)])^2 \]

donde \( G(z)=\frac{1}{a}\log{\cosh(az)} \) para \( 1\leq a\leq 2 \), \( a \) es una parámetro de ajuste y \( \cosh(z)=\frac{e^z + e^{-z}}{2} \)

Otra variante es usar \( G(z)=-e^{\frac{-u^2}{2}} \).

1. Obtenemos la estimación máximo verosímil de \( {\vec \epsilon} \) y de \( \sigma_{\vec \epsilon} \).
2. Aplicamos la transformación \( \vec{\epsilon_{T}}=\Phi^{-1}(F(\vec{\vec \epsilon})) \).
3. Obtenemos una estimación de \( \lambda \) mediante el método \( \mathcal{AC} \).
4. Con los valores de \( \vec{\epsilon_{T}} \) obtenidos a partir de \( \hat{\lambda} \) procedemos a calcular todos los valores \( J(\vec{\epsilon_{T}}) \) resultantes de eliminar una observación de la muestra.
5. Eliminamos la observación que haga mínimo el valor de la neguentropía.
6. Con los nuevos \( \overrightarrow{\epsilon_{T}} \) volvemos a estimar \( \hat{\lambda} \) mediante el método \( \mathcal{AC} \) y obtenemos los residuos transformados asociados a este nuevo \( \lambda \).
7. Realizamos de manera iterativa hasta excluir como máximo el 5% de las observaciones. Detendremos antes el procedimiento si la asimetría de la muestra pasa a ser negativa.
8. Nos quedamos con aquel \( \hat{\lambda} \) con una menor neguentropía en el caso de que se alcance el número máximo de iteraciones o con el ultimo valor en otro caso.

Estas técnicas ofrecen en términos generales mejores resultados, tanto para el problema de la asimetría positiva como para el resto de casos.

No obstante hay que ser consciente de que el grado de mejora, si bien es positivo, es moderado.

Esto es debido a que aunque tenemos cierta certeza de que cual es la causa del problema, resulta difícil identificar que observaciones generan las distorsiones y su posterior tratamiento.

En balance se ha comprobado que con los métodos que ofrecen un análisis más flexible se obtienen mejores resultados. Aquellos procedimientos que atenúan el efecto de las observaciones individuales para dar un mayor peso al comportamiento general de los datos resultan más adecuados para este problema.

Como producto de este TFM y de la tesis que ha servido de base para el desarrollo del trabajo se ha elaborado una librería del software estadístico R donde aparecen implementadas las técnicas anteriormente utilizadas.

Las funciones que recoge el paquete para distintos modelos de error compuesto son:

• area.comun
• ng.area.comun
• negentropy
• adkernel

El código puede consultarse y descargarse de https://github.com/javiergarcirod/stfa.

Para instalar y cargar el paquete cuando este disponible:

##library(devtools)
##url<-'https://github.com/javiergarcirod/stfa'
##install_git(url=url) 

## library(stfa)

set.seed(3)
b0<-100;b1<-5;landa<-1;n<-100;sigma2v<-1;
sigmav <- sqrt(sigma2v);sigma2u <- landa^2 * sigma2v;sigmau <- sqrt(sigma2u)
x <- runif(n, 1, 30);v <- rnorm(n, 0, sigmav);u <- abs(rnorm(n, 0, sigmau));e.real <- v - u
y <- b0 + b1 * x + e.real
l1<-stfa::area.comun(x,y,model="Normal-HalfNormal",method='Adaptative-Kernel')$landa_ac
l2<-stfa::area.comun(x,y,model="Normal-HalfNormal",method='Kernel')$landa_ac
cat("Para un valor del parámetro =",landa,"los resultados son:","\n","núcleo adaptable=",l1,'\n','núcleo=',l2)

Para un valor del parámetro = 1 los resultados son: 
 núcleo adaptable= 1 
 núcleo= 1.93204

• Método de MV con penalización utilizando otras funciones penalizadoras u otros valores del parámetro de sensibilidad.

• Nuevas técnicas de optimización cuántica aplicadas al campo de la estimación MV.

• Nuevos métodos de eliminación de valores anómalos basado en reglas heurísticas o metaheurísticas que tenga en cuenta la influencia de conjuntos de puntos.