El objetivo general consiste en estimar la distribución del ingreso en México utilizando la encuesta Nacional de Ingreso y Gasto de los Hogares (ENIGH), la cuál nos propociona el ingreso corriente de los hogares. Buscamos la compatibilidad de la ENIGH con dos fuentes de información adicionales: el Sistema de Cuentas Nacionales de México (SCNM) y el Servicio de Administración Tributaria (SAT).

El procedimiento se basa en realizar una modelación paramétrica a la variable ingreso corriente de la ENIGH suponiendo que es una variable aleatoria con una forma funcional intrínseca, denotada f(y|θ). Los parametros de esta función que denotamos por θ, resumirán la información clave de la muestra, permitiendonos reconstruir la distribución del ingreso.

La propuesta no es solo ajustar los datos que nos da la encuesta, sino en integrar y conciliar la información complementaria de dichas fuentes por medio de un ajuste con restrucciones.

Se explica, las fuentes de información utilizadas, en enfoque teorico que guia esta metodología y la implementación en R y los resultados de un ejercicio numérico.

DATOS

Las fuentes de información utilizadas y los detalles importantes que hay que saber de ellas

ENIGH

  • La encuesta es biaunal, por lo tanto, el proceso de ajuste de la distribución se ha implementado cada dos años a partir de 2004, abarcando tanto el ámbito nacional como el de entidad federativa.

  • Para la estimación a nivel entidad se toma en cuenta la fuente: Módulo de Condiciones Socioeconómicas (MSC), dada su capacitad de desagregación.

  • Estas encuestas presentan su información trimestralmente, por lo que consideramos los resultados de cada ejercicio que hacemos ya sea a nivel nacional o por entidad, también sean trimestralizados.

  • Se consideró excluir del análisis los hogares que registran un ingreso corriente trimestral igual o inferior a $2000 pesos. Esta decisión implica un ajuste en los factores de expansión para asegurar que la población total representada se mantenga constante.La reasignación del peso de los hogares eliminados se llevó a cabo de manera uniforme entre los hogares restantes, distribuyendo su representatividad proporcionalmente.

El siguiente código en R demuestra cómo descargar un archivo CSV directamente desde el sitio web del INEGI. Utilizaremos los datos de 2012 como ejemplo práctico.

url.zip<-"https://www.inegi.org.mx/contenidos/programas/enigh/tradicional/2012/microdatos/Tra_Concentrado_2012_concil_2010_csv.zip"
## para el caso de archivos más recientes por ejemplo, 2022 usar la siguiente ruta ## url.zip<-"https://www.inegi.org.mx/contenidos/programas/enigh/nc/2022/microdatos/enigh2022_ns_concentradohogar_csv.zip"

nombre_zip<-"Tra_Concentrado_2012_concil_2010_csv.zip" 

# para el caso 2022 el nombre del archivo es concentradohogar.csv

# crea un directorio "enigh_2012_datos" en la ruta actual de trabajo para que se #pueda guardar los datos. 
dir_unzip <- "enigh_2012_datos"
if (!dir.exists(dir_unzip)) {
  dir.create(dir_unzip)
}
 #descarga y descomprime  el archivo zip 
download.file(url.zip, destfile = nombre_zip, mode = "wb")
unzip(nombre_zip, exdir = dir_unzip)
#guarda la ruta donde se encuentra el archivo para despues leerlo
ruta_csv <- file.path(dir_unzip, "tra_concentrado_2012_concil_2010.csv") 
# caso 2022 concentradohogar.csv


#los datos con los que se trabajaran
data12<- read.csv(ruta_csv, encoding = "Latin-1", na.strings = c("", " "))
names(data12) # todas las variables que contiene el archivo 
file.remove(nombre_zip)
#los datos con los que se trabajaran
dir_unzip <- "enigh_2012_datos"
ruta_csv <- file.path(dir_unzip, "tra_concentrado_2012_concil_2010.csv")
data12<- read.csv(ruta_csv, encoding = "Latin-1", na.strings = c("", " "))
names(data12) # todas las variables que contiene el archivo 
##   [1] "folioviv"   "foliohog"   "ubica_geo"  "ageb"       "tam_loc"   
##   [6] "est_socio"  "est_dis"    "upm"        "factor_hog" "clase_hog" 
##  [11] "sexo_jefe"  "edad_jefe"  "educa_jefe" "tot_integ"  "hombres"   
##  [16] "mujeres"    "mayores"    "menores"    "p12_64"     "p65mas"    
##  [21] "ocupados"   "percep_ing" "perc_ocupa" "ing_total"  "ing_cor"   
##  [26] "ing_mon"    "trabajo"    "sueldos"    "horas_extr" "comisiones"
##  [31] "otra_rem"   "negocio"    "noagrop"    "industria"  "comercio"  
##  [36] "servicios"  "agrope"     "agricolas"  "pecuarios"  "reproducc" 
##  [41] "pesca"      "otros_trab" "rentas"     "utilidad"   "arrenda"   
##  [46] "transfer"   "jubilacion" "becas"      "donativos"  "remesas"   
##  [51] "bene_gob"   "otros_ing"  "gasto_nom"  "autoconsum" "remu_espec"
##  [56] "transf_esp" "transf_hog" "trans_inst" "estim_alqu" "percep_tot"
##  [61] "percep_mon" "retiro_inv" "prestamos"  "otras_perc" "erogac_nom"
##  [66] "gasto_tot"  "gasto_cor"  "gasto_mon"  "alimentos"  "ali_dentro"
##  [71] "cereales"   "carnes"     "pescado"    "leche"      "huevo"     
##  [76] "aceites"    "tuberculo"  "verduras"   "frutas"     "azucar"    
##  [81] "cafe"       "especias"   "otros_alim" "bebidas"    "ali_fuera" 
##  [86] "tabaco"     "vesti_calz" "vestido"    "calzado"    "vivienda"  
##  [91] "alquiler"   "pred_cons"  "agua"       "energia"    "limpieza"  
##  [96] "cuidados"   "utensilios" "enseres"    "salud"      "atenc_ambu"
## [101] "hospital"   "medicinas"  "transporte" "publico"    "foraneo"   
## [106] "adqui_vehi" "mantenim"   "refaccion"  "combus"     "comunica"  
## [111] "educa_espa" "educacion"  "esparci"    "paq_turist" "personales"
## [116] "cuida_pers" "acces_pers" "otros_gas"  "transf_gas" "erogac_tot"
## [121] "erogac_mon" "cuota_viv"  "mater_serv" "material"   "servicio"  
## [126] "deposito"   "prest_terc" "pago_tarje" "deudas"     "balance"   
## [131] "otras_erog" "smg"
#file.remove(nombre_zip)

El siguiente código es para obtener el vector de valores de ingresos y factores de expansión

 muestra<-data12$ing_cor
 fac<-data12$factor # aveces el factor de expansión tiene el nombre "factor_hog"
 total_ingreso <-sum(muestra*fac) # 
 tot_hog<-sum(fac) 
 print(paste("total del ingreso de la muestra: ",total_ingreso,sep=""))
## [1] "total del ingreso de la muestra: 1203202597623.71"
 print(paste("total de hogares en la población: ",tot_hog,sep=""))
## [1] "total de hogares en la población: 31559379"

Se quitan los valores menores a 2000 y se redistribuye la población

w<-c()
w<-which(muestra<2000)
if(length(w)==0)
{muestra<-muestra} else 
{muestra<-muestra[-w]
 fac<-fac[-w]
}
fquitados<-tot_hog-sum(fac)  ## re-austar los valores de expansión 
coci<-fquitados/length(fac)
nfac<-(fac+coci)
fac<-nfac  ## nuevos factores ajustados

SCNM

El sistema de Cuentas Nacionales aporta distintos tipos de ingresos anuales en hogares. Constantemente actualizan sus valores debido algunos cambios necesarios ya sea de carácter metodológico o de cambios de años base, por ejemplo el cambio de base de 2003 a 2008 y el reciente cambio de 2008 a 2013. El ingreso que se usa en este proyecto es el numerado por B.7n Ingreso Disponible Ajustado Neto (IDAN), ya que este ingreso no incluye la depreciación de activos y sí contempla las transferencias en especie, ingreso reportado por ENIGH.

  • Obtenemos el promedio trimestral nacional de hogares de ingreso según Cuentas Nacionales usando como denominador el total nacional de hogares (THogares) que reporta ENIGH:

\[mean_{SCNM} = \frac{\frac{IDAN}{4}}{Thogares}\]

  • Para el caso entidad, el promedio se calcula obteniendo una proporción: el ingreso corriente total de la entidad y el ingreso corriente nacional, ambos obtenidos de MCS. Este cociente lo denotamos por (PP). El promedio de ingreso de cualquier entidad es:

\[mean_{SCNM_{entidad}}= \frac{\frac{PP*IDAN}{4}}{Thogares_{entidad}}\]

Para obtener el promedio de cuentas nacionales

idan<-  11635413   # valor de IDAN del 2012 
tot_cn<-idan *1000000 # se multiplica por 1 millón  
ci<-c(tot_cn/4/tot_hog)  # trimestralizado 
print(paste("media de cuentas nacionales_2012= ",ci,sep=""))
## [1] "media de cuentas nacionales_2012= 92170.8012695687"
print(ci)
## [1] 92170.8

SAT

  • Los datos proporcionados por SAT son declaraciones por contribuyentes anualizados. Cada registro de los archivos están identificados de manera única por un número consecutivo RFC_ID. Se indican en la mayoría de los casos, su entidad de procedencia.

  • El archivo de datos usados para nuestros ajustes son las declaraciones didácticas donde se proporciona el total de ingresos acumulables denotado como i_dec_tiaonct1.

  • Nos interesa el ingreso declarado de los más ricos.

  • Como la información del contribuyente es anual, entonces el ingreso de los más ricos tendrá que ser dividido por 4 para obtener datos trimestralizados.

  • Para el caso del ajuste a nivel nacional para algún año de interés nos interesa: el promedio de ingresos declarado al SAT por el millonesimo percentil(.999999) que corresponde los 32 más ricos, y el ingreso del treintaidosavo más rico.

  • Para el caso de entidades la referencia es a partir de los 20 más ricos.

p<-0.999999
mean_Sat_anual<- 1624174001.00 # promedio del ingreso de los 32 más ricos anual                                     2012
mean_sat<-mean_Sat_anual/4 # trimestralizado 2012
qi<-177965502 # ingreso del 32 del top-32 más ricos ya trimestralizado
print(mean_sat)
## [1] 406043500

PROCESO DE AJUSTE

Se probaron varias funciones de densidad paramétricas para el ajuste de los datos de ingreso. Entre ellas Gamma Generalizada (GG) con 3 parámetros dada por:

\[f_{y}(y|\mu,\sigma, \nu)=\frac{|\nu|\theta^{\theta}z^{\theta}exp(-\theta z)}{\Gamma(\theta)y}\] Donde \[ z = (y / \mu)^{\nu}\]

y \[ \theta = 1 / \sigma^{2} \nu^{2}\].

La esperanza de Y es:

\[ E(Y,\mu,\sigma,\nu) = \mu\Gamma(\theta + \frac{1}{\nu})/[\theta^{1 / \nu} \Gamma(\theta)]\] El proceso se reduce a encontrar la estimación de los párametros de la función Gamma Generaliza por medio de log-verosimilitud sujeta a:

\[E(Y,\mu,\sigma,\nu)=mean_{SCNM}\]

\[ (\int_{q_{i}}^{Inf} yf(y)dy)/q = mean_{32másricos} \]

Donde \({\bf q}\) es la proporción de hogares que tienen ingresos igual o menor que \({\bf q_{i}}\).

Implementación de la Distribución gamma tres parámetros

LAs librerías en R necesarias para el ajuste

#install.packages("gamlss")
#install.packages("alabama")
library(gamlss)
library(alabama)

La siguiente función que llamamos fgamG es la de densidad Gama Generalizada (GG) con pesos que se ajusta en este análisis

fgamG<-function(x,m=muestra,w=fac) ## funcion de verosimilitud
{ 
  fi<-c()
  fi<- -1*w*dGG(m, x[1], x[2],x[3],log=TRUE)
  sumaf<-sum(fi)
  return(sumaf)
}

La función hin implementa las restricciones de desigualdad: \(\mu\),\(\sigma\) > 0

hin<-function(x)
{
  h <- rep(NA, 2)
  h[1]<-x[1]
  h[2]<-x[2]
  h
}

La función heq implementa las restricciones de igualdad:

  • La esperanza de la función de densidad sea exactamente el valor del promedio de SCN.
  • La esperanza truncada de la distribución a partir del quantil \(qi\) a \(Inf\) sea exactamente el promedio de los 32(20) más ricos.
# esta función cambia de acuerdo a la función de densidad de ajuste
heq<- function(x)  
{
  #print(x)
  p=p
  par<-x
  t<-1/(x[2]^2*x[3]^2)
  h <- rep(NA, 2)
  inte<-c()
  inte <-try(integrate(function(y) y*(dGG(y, mu=par[1],sigma=par[2],nu=par[3])), lower=qi, upper=1000000000000,stop.on.error = FALSE)$v,silent=TRUE)
  if (is.character(inte))
    inte<-1000
  h[1]<-inte-prom
  h[2]<-ci-(x[1]*gamma(t+(1/x[3]))/(t^(1/x[3])*gamma(t)))
  h
  # print(h)
}

Caso no restringido

El siguiente paso, es el de obtener una estimación preliminar de los parámetros de la función GG a través de la función gmag, posteriormente se utiliza la función constrOptim de la librería Alabama para que se finalice el proceso de estimación, el vector xGGsin es el que se obtiene como estimación final.

gammues<-gamlss(muestra~1, family=GG,weights=fac)
## GAMLSS-RS iteration 1: Global Deviance = 723264570 
## GAMLSS-RS iteration 2: Global Deviance = 720775085 
## GAMLSS-RS iteration 3: Global Deviance = 719896509 
## GAMLSS-RS iteration 4: Global Deviance = 719606898 
## GAMLSS-RS iteration 5: Global Deviance = 719511905 
## GAMLSS-RS iteration 6: Global Deviance = 719480407 
## GAMLSS-RS iteration 7: Global Deviance = 719469851 
## GAMLSS-RS iteration 8: Global Deviance = 719466285 
## GAMLSS-RS iteration 9: Global Deviance = 719465075 
## GAMLSS-RS iteration 10: Global Deviance = 719464663 
## GAMLSS-RS iteration 11: Global Deviance = 719464523 
## GAMLSS-RS iteration 12: Global Deviance = 719464475 
## GAMLSS-RS iteration 13: Global Deviance = 719464459 
## GAMLSS-RS iteration 14: Global Deviance = 719464453 
## GAMLSS-RS iteration 15: Global Deviance = 719464451 
## GAMLSS-RS iteration 16: Global Deviance = 719464450 
## GAMLSS-RS iteration 17: Global Deviance = 719464450 
## GAMLSS-RS iteration 18: Global Deviance = 719464450 
## GAMLSS-RS iteration 19: Global Deviance = 719464450 
## GAMLSS-RS iteration 20: Global Deviance = 719464450
## Warning in RS(): Algorithm RS has not yet converged
gmag<-refit(gammues) # repetir el proceso en caso de no convergencia
## GAMLSS-RS iteration 21: Global Deviance = 719464450 
## GAMLSS-RS iteration 22: Global Deviance = 719464450 
## GAMLSS-RS iteration 23: Global Deviance = 719464450
    mu<-fitted(gmag,"mu")[1]
    sigma<-fitted(gmag,"sigma")[1]
    nu<-fitted(gmag,"nu")[1]
    theta<-c(mu=mu,sigma=sigma,nu=nu) # estimación inicial 
# ajuste con alabama sin restricciones, notar que theta se usa como valor inicial # en el ajuste 
gammues<-c()
x<-c()
gammues<-constrOptim.nl(par=theta, fn=fgamG, hin=hin)
## Min(hin):  0.7990935 
## par:  24787.32 0.7990935 -0.2469241 
## fval:  359732225 
## Min(hin):  0.7990934 
## par:  24787.32 0.7990934 -0.2469241 
## fval:  359732225 
## Min(hin):  0.7990934 
## par:  24787.32 0.7990934 -0.2469241 
## fval:  359732225 
## Min(hin):  0.7990934 
## par:  24787.32 0.7990934 -0.2469241 
## fval:  359732225
xGGsin<-gammues$par
cat("vector estimado del parámetro de la función GG sin restricciones =\n")
## vector estimado del parámetro de la función GG sin restricciones =
print(xGGsin)
##          mu.1       sigma.1          nu.1 
## 24787.3206428     0.7990934    -0.2469241

Es importante tomar en cuenta que la calidad de las estimaciones depende significativamente de los valores iniciales que se proporcionen. Para asegurar que obtenemos la mejor solución posible, es necesario ejecutar el proceso de estimación con múltiples conjuntos de valores iniciales y, posteriormente, seleccionar aquella que maximice la verosimilitud.

Para obtener estadísticos de la función no restringida se tiene las siguientes instrucciones

x<-xGGsin
t<-1/(x[2]^2*x[3]^2)
ExGGsin<-(x[1]*gamma(t+(1/x[3])))/(t^(1/x[3])*gamma(t)) # Esperanza de GG no                                                                restringida
vGGsin<-gammues$value   # valor de la verosimilitud 
# quantiles 

# vector probabilidades para quantiles
probs = c(0.1, 0.5, 1, 2, 5, 10,20,30,40, 50,60,70,80,90,95,98,99,99.5,99.9)/100  
par<-xGGsin
qGGsin<-qGG(probs, mu=par[1],sigma=par[2], nu=par[3],lower.tail = TRUE, 
            log.p = FALSE)
cat("Quantiles para distintas probabilidades de GG sin restricciones =\n")
## Quantiles para distintas probabilidades de GG sin restricciones =
print(cbind(probs,quantiles=qGGsin))
##       probs  quantiles
##  [1,] 0.001   2744.726
##  [2,] 0.005   3878.791
##  [3,] 0.010   4605.042
##  [4,] 0.020   5571.655
##  [5,] 0.050   7460.109
##  [6,] 0.100   9730.087
##  [7,] 0.200  13533.767
##  [8,] 0.300  17269.693
##  [9,] 0.400  21354.214
## [10,] 0.500  26130.664
## [11,] 0.600  32084.474
## [12,] 0.700  40118.826
## [13,] 0.800  52383.905
## [14,] 0.900  76567.617
## [15,] 0.950 105677.408
## [16,] 0.980 153378.206
## [17,] 0.990 197800.189
## [18,] 0.995 250730.447
## [19,] 0.999 414404.476

Caso restringido

La función encargada de encontrar la solución al sistema no lineal con restricciones es constOptim.nl de la librería Alabama.

par<-xGGsin
q<-1-p  # porcentaje de hogares con ingreso menor a qi
prom<-mean_sat*q # proporción de hogares que tienen ingresos igual o menor que qi
ansgam<- constrOptim.nl(par=par, fn=fgamG,heq=heq,hin=hin)
## Min(hin):  0.7990934 Max(abs(heq)):  54297.3 
## Outer iteration:  1 
## Min(hin):  0.7990934 Max(abs(heq)):  54297.3 
## par:  24787.3 0.799093 -0.246924 
## fval =   359700000 
##  
## Outer iteration:  2 
## Min(hin):  0.8031024 Max(abs(heq)):  54297.3 
## par:  24887.4 0.803102 -0.258712 
## fval =   359700000 
##  
## Outer iteration:  3 
## Min(hin):  0.8608035 Max(abs(heq)):  10889.48 
## par:  20307.8 0.860804 -1.01438 
## fval =   361300000 
##  
## Outer iteration:  4 
## Min(hin):  0.8680892 Max(abs(heq)):  10889.48 
## par:  20278.2 0.868089 -1.03739 
## fval =   361500000 
##  
## Outer iteration:  5 
## Min(hin):  0.9059414 Max(abs(heq)):  10889.48 
## par:  21994.8 0.905941 -0.910707 
## fval =   361600000 
##  
## Outer iteration:  6 
## Min(hin):  0.985292 Max(abs(heq)):  1692.016 
## par:  25619.9 0.985292 -0.686764 
## fval =   362300000 
##  
## Outer iteration:  7 
## Min(hin):  1.013595 Max(abs(heq)):  1692.016 
## par:  27027.6 1.01359 -0.612439 
## fval =   362600000 
##  
## Outer iteration:  8 
## Min(hin):  1.043101 Max(abs(heq)):  274.0612 
## par:  28623.5 1.0431 -0.53475 
## fval =   3.63e+08 
##  
## Outer iteration:  9 
## Min(hin):  1.053588 Max(abs(heq)):  274.0612 
## par:  29218.7 1.05359 -0.507248 
## fval =   363100000 
##  
## Outer iteration:  10 
## Min(hin):  1.060648 Max(abs(heq)):  21.15123 
## par:  29635.2 1.06065 -0.488604 
## fval =   363200000 
##  
## Outer iteration:  11 
## Min(hin):  1.060647 Max(abs(heq)):  21.15123 
## par:  29635.2 1.06065 -0.488603 
## fval =   363200000 
##  
## Outer iteration:  12 
## Min(hin):  1.062819 Max(abs(heq)):  0.2946807 
## par:  29726.8 1.06282 -0.483538 
## fval =   363300000 
##  
## Outer iteration:  13 
## Min(hin):  1.062819 Max(abs(heq)):  0.2946807 
## par:  29726.8 1.06282 -0.483538 
## fval =   363300000 
##  
## Outer iteration:  14 
## Min(hin):  1.062819 Max(abs(heq)):  0.2946807 
## par:  29726.8 1.06282 -0.483538 
## fval =   363300000 
##  
## Outer iteration:  15 
## Min(hin):  1.062819 Max(abs(heq)):  0.2946807 
## par:  29726.8 1.06282 -0.483538 
## fval =   363300000 
##  
## Outer iteration:  16 
## Min(hin):  1.062586 Max(abs(heq)):  0.0005080144 
## par:  29735.5 1.06259 -0.483755 
## fval =   363300000 
##  
## Outer iteration:  17 
## Min(hin):  1.062586 Max(abs(heq)):  4.18819e-05 
## par:  29735.5 1.06259 -0.483755 
## fval =   363300000 
##  
## Outer iteration:  18 
## Min(hin):  1.062586 Max(abs(heq)):  4.18819e-05 
## par:  29735.5 1.06259 -0.483755 
## fval =   363300000 
##  
## Outer iteration:  19 
## Min(hin):  1.062586 Max(abs(heq)):  4.18819e-05 
## par:  29735.5 1.06259 -0.483755 
## fval =   363300000 
##  
## Outer iteration:  20 
## Min(hin):  1.062586 Max(abs(heq)):  4.18819e-05 
## par:  29735.5 1.06259 -0.483755 
## fval =   363300000 
##  
## Outer iteration:  21 
## Min(hin):  1.062586 Max(abs(heq)):  6.611343e-06 
## par:  29735.5 1.06259 -0.483755 
## fval =   363300000 
##  
## Outer iteration:  22 
## Min(hin):  1.062586 Max(abs(heq)):  6.611343e-06 
## par:  29735.5 1.06259 -0.483755 
## fval =   363300000 
##  
## Outer iteration:  23 
## Min(hin):  1.062586 Max(abs(heq)):  6.611343e-06 
## par:  29735.5 1.06259 -0.483755 
## fval =   363300000 
## 
xGGcon<-ansgam$par

Para asegurar la convergencia del ajuste, es importante verificar los resultados de la función constrOptim.nl. Por un lado, el valor de Max(abs(heq)) debe ser muy pequeño, cercano a cero. Esto garantiza que las restricciones impuestas por heq (h[1] y h[2]) se cumplan de manera efectiva, ie, que sus valores tiendan a cero. Por otro lado, el valor de fval debe ser monotónomente creciente en cada iteración, esto demuestra que el algoritmo está encontrando una mejor solución en cada paso.

cat("vector estimado del parámetro de la función GG con restricciones =\n")
## vector estimado del parámetro de la función GG con restricciones =
print(xGGcon)
##          mu.1       sigma.1          nu.1 
## 29735.4923450     1.0625863    -0.4837548

Para obtener estadisticos necesarios de esta función restringida, simplemente guardamos la salida xGGcon que contiene la estimación de los parámetros de la función GG bajo restricciones aplicadas.

x<-xGGcon
t<-1/(x[2]^2*x[3]^2)
ExGGcon<-x[1]*gamma(t+(1/x[3]))/(t^(1/x[3])*gamma(t))
cat("Esperanza de la función GG=\n")
## Esperanza de la función GG=
print(ExGGcon)
##    mu.1 
## 92170.8
vGGcon<-ansgam$value ## verosimilitud

par<-xGGcon
qGGcon<-qGG(probs, mu=par[1],sigma=par[2], nu=par[3],lower.tail = TRUE, log.p = FALSE)
cat("Quantiles para distintas probabilidades de GG con restricciones =\n")
## Quantiles para distintas probabilidades de GG con restricciones =
print(cbind(probs,quantiles=qGGcon))
##       probs   quantiles
##  [1,] 0.001    3597.322
##  [2,] 0.005    4552.048
##  [3,] 0.010    5152.217
##  [4,] 0.020    5947.668
##  [5,] 0.050    7512.977
##  [6,] 0.100    9444.735
##  [7,] 0.200   12862.841
##  [8,] 0.300   16486.242
##  [9,] 0.400   20785.789
## [10,] 0.500   26297.458
## [11,] 0.600   33948.842
## [12,] 0.700   45747.276
## [13,] 0.800   67307.652
## [14,] 0.900  124046.755
## [15,] 0.950  221737.401
## [16,] 0.980  467740.997
## [17,] 0.990  816251.786
## [18,] 0.995 1419688.270
## [19,] 0.999 5103330.179

Cociente entre el decil más bajo con el decil más alto

     q10<-qGG(p=c(.10),mu=par[1],sigma=par[2],nu=par[3])
     q90<-qGG(p=c(.90),mu=par[1],sigma=par[2],nu=par[3])
     
     i0_10<-integrate(function(x) x*dGG(x=x, mu=par[1], sigma=par[2],nu=par[3]), 
                    lower= 0, upper=q10)
     i0_90<-integrate(function(x) x*dGG(x=x, mu=par[1], sigma=par[2],nu=par[3]), 
                     lower= 0, upper=q90 )
     num<-ExGGcon-i0_90$v
     numGG<-num     # valor promedio de ingreso para el ultimo decil (90-10)
     denGG<-i0_10$v  # valor promedio de ingreso para el primer decil 
     resGG<-num/denGG   # cociente entre primer decil y ultimo decil 
     mGG<-max(muestra)  #  máximo de la muestra 
     i<-integrate(function(x) x*dGG(x=x, mu=par[1], sigma=par[2],nu=par[3]), 
                  lower= 0, upper=mGG )
     cocGG<-ExGGcon-i$v   ## ingreso promedio mayor al máximo observado en muestra      
     cocienteMax <- cocGG/ExGGcon #Relación ingreso del máximo observado a ingresos totales

Promedios de ingreso para distintos quantiles

   par<-xGGcon
   pes<-c( 0, 0.1,.2,.3,.4,.5,.6,.7,.8,.9,.99,.999999)
   qq<-c()
   ii<-c()
   qq[1]<-0
   for ( i in 2:12)
     qq[i]<-qGG(p=pes[i],mu=par[1],sigma=par[2],nu=par[3])
   
   for ( i in 2:12)
     ii[i]<-integrate(function(x) x*dGG(x=x, mu=par[1],sigma=par[2],nu=par[3]), 
                      lower= qq[i-1], upper=qq[i])$v
   
   
   ii[i+1]<-  ii[i+1]<-integrate(function(x) x*dGG(x=x, mu=par[1],sigma=par[2],nu=par[3]), lower= qq[12], upper=100000000000)$v
   ii<-ii[-1]
   ii[13]<-sum(ii)
   if (ii[13]==ExGGcon) 
     ii[13]<-ii[13] else 
     {ii[12]<-ExGGcon-sum(ii[-c(12,13)])
      ii[13]<-ExGGcon}
   iiGG<-ii
   
   deciles<-c(
     "0---10",
     "10---20",
     "20---30",
     "30---40",
     "40---50",
     "50---60",
     "60---70",
     "70---80",
     "80---90",
     "90---99",
     "99---999999",
     "999999--Inf","suma")
   perGG<-qq[12]  ## percentil .9999999 
   
   promedios<-cbind(deciles,iiGG)
  print(promedios)
##       deciles       iiGG              
##  [1,] "0---10"      "709.316787243632"
##  [2,] "10---20"     "1257.75040111376"
##  [3,] "20---30"     "1779.83314814099"
##  [4,] "30---40"     "2383.20292240231"
##  [5,] "40---50"     "3136.47613665463"
##  [6,] "50---60"     "4142.5868820134" 
##  [7,] "60---70"     "5600.67705316391"
##  [8,] "70---80"     "7989.02273584048"
##  [9,] "80---90"     "12938.4872162444"
## [10,] "90---99"     "30129.0863057142"
## [11,] "99---999999" "21702.8185972577"
## [12,] "999999--Inf" "401.543088942926"
## [13,] "suma"        "92170.8012747323"

Si se desea obtener el gini se puede utilizar la siguiente función. Como entrada necesita un vector de ingresos y de salida se obtiene el gini de esos ingresos

gini2<- function(x, unbiased = TRUE, na.rm = FALSE){
    if (na.rm)
      x <- x[!na.ind]
    n <- length(x)
    mu <- mean(x)
    N <- if (unbiased) n * (n - 1) else n * n
    ox <- x[order(x)] ## quitar la informacion irrelevante DROP 
    ## la funci?n crossprod = productos cruzados 
    dsum <- drop(crossprod(2 * 1:n - n - 1,  ox)) 
    dsum / (mu * N)
  }
 set.seed(10)
   x<-xGGcon
   ma<- rGG(1000, x[1], x[2],x[3])
   giniGG<-gini2(ma)
   print(giniGG)
## [1] 0.73443

Tambien puede usarse la función ya declarada “gini” de la librería laeken

#install.packages("laeken")
library(laeken)
gini(ma)
## Value:
## [1] 73.36956