###############################################
# Día 3: estimación de medidas de desigualad
###############################################

# Práctica 1:

ingreso.a <- c(13, 25, 15, 7, 12, 38, 42, 53, 7)
ingreso.b <- c(4, 23, 36, 18, 39, 20, 9, 45, 12)

# Curva de lorenz de la localidad A

sort(ingreso.a)
length(ingreso.a)
sum(ingreso.a)

prop.renta.a <- sort(ingreso.a)/sum(ingreso.a)
prop.renta.a
prop.pop.a <- rep(1/length(ingreso.a), length(ingreso.a))
prop.pop.a

cum.renta.a <- c(0, cumsum(prop.renta.a))
cum.pop.a <- c(0, cumsum(prop.pop.a))
cbind(cum.renta.a, cum.pop.a)


prop.renta.b <- sort(ingreso.b)/sum(ingreso.b)
prop.pop.b <- rep(1/length(ingreso.b), length(ingreso.b))
cum.renta.b <- c(0, cumsum(prop.renta.b))
cum.pop.b <- c(0, cumsum(prop.pop.b))


# Representación gráfica de la curva de Lorenz
plot(cum.pop.a, cum.renta.a, type = "o",
     xlab = "Porcentaje de población", 
     ylab = "Proporción de renta", 
     main = "Curva de Lorenz de las localidades A y B", 
     pch = 20, col= "Red")
abline(0, 1, lwd = 2, lty = 1)
points(cum.pop.b, cum.renta.b, type = "o", col= "blue", pch=20)


# Práctica 2:
library(ineq)

# Curva de Lorenz con ineq
lc.a <- Lc(ingreso.a)
lc.a
# Representación gráfica de la curva de Lorenz
plot(lc.a$p, lc.a$L, type = "o",
     xlab = "Porcentaje de población", 
     ylab = "Proporción de renta", 
     main = "Curva de Lorenz de las localidades A y B", 
     pch = 20, col= "Red")
abline(0, 1, lwd = 2, lty = 1)
lc.b <- Lc(ingreso.b)
points(lc.b$p, lc.b$L, type = "o", col= "blue", pch=20)

#Localidad A
# Índice de Gini
Gini(ingreso.a)
# Índice de Atkinson
Atkinson(ingreso.a, parameter = 1.25)
#Índice de Theil
entropy(ingreso.a, parameter = 1)
# MLD
entropy(ingreso.a, parameter = 0)

#Localidad B
# Índice de Gini
Gini(ingreso.b)
# Índice de Atkinson
Atkinson(ingreso.b, parameter = 1.25)
#Índice de Theil
entropy(ingreso.b, parameter = 1)
# MLD
entropy(ingreso.b, parameter = 0)

# La práctica no lo pide
Atkinson(ingreso.a, parameter = 2)
Atkinson(ingreso.b, parameter = 2)

# Efecto de impuesto fijo frente a proporcional
# Ingreso neto con impuesto fijo de 1000$
ingreso.a1 <- ingreso.a - 1

# Ingreso neto con impuesto proporcional del 10%
ingreso.a2 <- ingreso.a * 0.9

Gini(ingreso.a)
Gini(ingreso.a1)
Gini(ingreso.a2)


####################################################
# Estimación de medidas de desigualdad con 
# información limitada
#####################################################

shares <- c(3.2, 8.3, 14.4, 23, 51)
shares.c <- cumsum(shares) / sum(shares)
shares.c <- c(0, shares.c)
shares.c
pop <- seq(0, 1, 0.2)

# Representación gráfica de la curva de Lorenz
plot(pop, shares.c, type = "o",
     xlab = "Porcentaje de población", 
     ylab = "Proporción de renta", 
     main = "Curva de Lorenz de USA (2012)", 
     pch = 20, col= "Red")
abline(0, 1, lwd = 2, lty = 1)

# Cota baja del índice de Gini

cbind(shares.c, pop)

gini.emp <- 1 - sum((shares.c[-1] + shares.c[-length(shares.c)]) * 
                      (pop[-1] - pop[-length(pop)])) 
gini.emp

##################################################
# Práctica 5
##################################################
library(GB2group)
library(ineq)
library(GB2)

fitgroup.gb2(shares, gini.e = 0.477, pc.inc = 50626, rescale = 1)

# Estimación de medidas de desigualdad por simulación paramétrica (Monte Carlo)
ajuste.gb2 <- fitgroup.gb2(shares, gini.e = 0.477, pc.inc = 50605, rescale = 1)
param.gb2 <- ajuste.gb2$nls.estimation[1, ]

# Paso 1 - simulo una muestra de mi modelo paramétrico
# Tamaño de muestra grande
N <- 10^6
set.seed(7)
sim.sam <- rgb2(N, param.gb2[1], param.gb2[2], param.gb2[3], param.gb2[4])

# Paso 2: Calcular las medidas de desigualdad deseadas
Gini(sim.sam)
Atkinson(sim.sam, parameter = 1.25)
entropy(sim.sam, parameter = 1) # Índice de Theil
entropy(sim.sam, parameter = 0) # MLD
var(sim.sam) * (N - 1) / N # Varianza
Gini(sim.sam) * mean(sim.sam)

# Distribución Beta 2
# Estimación de medidas de desigualdad por simulación paramétrica (Monte Carlo)
ajuste.b2 <- fitgroup.b2(shares, gini.e = 0.477, pc.inc = 50605, rescale = 1)
param.b2 <- ajuste.b2$nls.estimation[1, ]

# Paso 1 - simulo una muestra de mi modelo paramétrico (BETA 2)
# Tamaño de muestra grande
N <- 10^6
set.seed(7)
sim.sam <- rgb2(N, 1, param.b2[1], param.b2[2], param.b2[3])

# Paso 2: Calcular las medidas de desigualdad deseadas
Gini(sim.sam)
Atkinson(sim.sam, parameter = 1.25)
entropy(sim.sam, parameter = 1) # Índice de Theil
entropy(sim.sam, parameter = 0) # MLD
var(sim.sam) * (N - 1) / N # Varianza
Gini(sim.sam) * mean(sim.sam)


# Distribución Singh- Maddala
# Estimación de medidas de desigualdad por simulación paramétrica (Monte Carlo)
ajuste.sm <- fitgroup.sm(shares, gini.e = 0.477, pc.inc = 50626, rescale = 1)
param.sm <- ajuste.sm$nls.estimation[1, ]

# Paso 1 - simulo una muestra de mi modelo paramétrico (BETA 2)
# Tamaño de muestra grande
N <- 10^6
set.seed(7)
sim.sam <- rgb2(N, param.sm[1], param.sm[2], 1, param.sm[3])

# Paso 2: Calcular las medidas de desigualdad deseadas
Gini(sim.sam)
Atkinson(sim.sam, parameter = 1.25)
entropy(sim.sam, parameter = 1) # Índice de Theil
entropy(sim.sam, parameter = 0) # MLD
var(sim.sam) * (N - 1) / N # Varianza
Gini(sim.sam) * mean(sim.sam)

# Distribución Dagum
# Estimación de medidas de desigualdad por simulación paramétrica (Monte Carlo)
ajuste.da <- fitgroup.da(shares, gini.e = 0.477, pc.inc = 50626, rescale = 1)
param.da <- ajuste.da$nls.estimation[1, ]

# Paso 1 - simulo una muestra de mi modelo paramétrico (dagum)
# Tamaño de muestra grande
N <- 10^6
set.seed(7)
sim.sam <- rgb2(N, param.da[1], param.da[2], param.da[3], 1)

# Paso 2: Calcular las medidas de desigualdad deseadas
Gini(sim.sam)
Atkinson(sim.sam, parameter = 1.25)
entropy(sim.sam, parameter = 1) # Índice de Theil
entropy(sim.sam, parameter = 0) # MLD
var(sim.sam) * (N - 1) / N # Varianza
Gini(sim.sam) * mean(sim.sam)


# Representación gráfica de los ajustes
fit.plot(c("ajuste.gb2", "ajuste.b2"), fit.type = 1, fit.legend = TRUE,
         l.size = 0.7)
abline(0, 1)

# Calculamos medidas de bondad del ajuste: 
# Vamos a calcular la desviación cuadrática media

sse.gb2 <- ajuste.gb2$nls.rss/(length(shares) - 1)
sse.gb2

sse.b2 <- ajuste.b2$nls.rss/(length(shares) - 1)
sse.b2

sse.sm <- ajuste.sm$nls.rss/(length(shares) - 1)
sse.sm

sse.da <- ajuste.da$nls.rss/(length(shares) - 1)
sse.da

# Vamos a calcular AIC

aic.gb2 <- sse.gb2 * exp(2 * (length(param.gb2) - 1) / (length(shares) - 1))
aic.gb2

aic.b2 <- sse.b2 * exp(2 * (length(param.b2) - 1) / (length(shares) - 1))
aic.b2

aic.sm <- sse.sm * exp(2 * (length(param.sm) - 1) / (length(shares) - 1))
aic.sm

aic.da <- sse.da * exp(2 * (length(param.da) - 1) / (length(shares) - 1))
aic.da

#######################################################
# Estimación de medidas de desigualdad a nivel regional
#########################################################

shares <- c(3.2, 8.3, 14.4, 23, 51)
#Ajusto GB2 en USA
ajuste.gb2 <- fitgroup.gb2(shares, gini.e = 0.477, pc.inc = 50605, rescale = 1)
param.gb2 <- ajuste.gb2$nls.estimation[1, ]

renta.c <- c(2.2, 4, 5.2, 6.5, 7.9, 9.2, 10.7, 12.6 ,15.4, 26.1)
#Ajusto GB2 en Canadá
ajuste.gb2.c <- fitgroup.gb2(renta.c, gini.e = 0.354, pc.inc = 41622, rescale = 1)
param.gb2.c <- ajuste.gb2.c$nls.estimation[1, ]


#Estimamos las medidas de entropía para los distintos valores del parámetro de sensibilidad
# En USA
N <- 10^6
set.seed(7)
sim.sam <- rgb2(N, param.gb2[1], param.gb2[2], param.gb2[3], param.gb2[4])
ge1.usa <- entropy(sim.sam, parameter = 1) # Índice de Theil
ge0.usa <- entropy(sim.sam, parameter = 0) # MLD
ge15.usa <- entropy(sim.sam, parameter = 1.5) # Índice de Theil

# En Canada
N <- 10^6
set.seed(7)
sim.sam <- rgb2(N, param.gb2.c[1], param.gb2.c[2], param.gb2.c[3], param.gb2.c[4])
ge1.can <- entropy(sim.sam, parameter = 1) # Índice de Theil
ge0.can <- entropy(sim.sam, parameter = 0) # MLD
ge15.can <- entropy(sim.sam, parameter = 1.5) # Índice de Theil


# Creamos las proporciones de población
poblacion <- c(314043872, 34900705)
li <- poblacion/sum(poblacion)
li

# Creamos las proporciones de renta
medias <- c(50605, 41622)
media.r <- sum(li * medias)
media.r
si <- li * medias / media.r
si

# Una vez que tenemos todos los componentes, calculamos las medidas regionales
# Entropía generalizada (theta = 1.5)
ge15 <- c(ge15.usa, ge15.can)
theta <- 1.5
GEw <- sum(li^(1 - theta) * si^theta * ge15)
GEw

GEb <- 1/(theta * (theta  - 1)) * (sum(li * (medias / media.r)^theta) - 1)
GEb

# medida de entropía generalizada en Norte América
GEw + GEb

# Índice de Theil
ge1 <- c(ge1.usa, ge1.can)
Tw <- sum(si * ge1)
Tb <- sum(si * log(medias / media.r))
Tw; Tb

# Índice de Theil regional
Tw + Tb

# MLD
ge0 <- c(ge0.usa, ge0.can)
Lw <- sum(li * ge0)
Lb <- sum(li * log(media.r / medias))
# MLD regional
Lw +Lb