Distribución por Ingreso y Edad de la población.

Se obtuvieron datos del INEGI (distribución de edades de la población y deciles de ingreso familiar), así como del US Census (distribución de edades e ingreso). Con esto se estima una distribución de edades e ingreso para México.

La idea de hacer esto es poder simular de forma relativamente realista una distribución demográfica de México en términos de ingreso y edad, con el fin de hacer un muestreo con características similares que nos sirva como base de nuestras estrategias, en particular en el sector de seguros.

Los datos se obtuvieron de:

• https://www.census.gov/hhes/www/cpstables/032013/perinc/pinc01_000.htm
• http://www3.inegi.org.mx/sistemas/descarga/

Se hizo una primera limpieza de los datos y se colocó en el archivo: “CombinaIncomeEdad.csv”

Notas:

• Sólo se incluye la población con ingresos menores a US$100,000 anuales y edades entre 25 y 80 años. En particular cortar el nivel de ingreso implica cortar la distribución de cola ancha, pero el modelo es independiente de esto al ser no paramétrico y se puede aplicar igual.
• MÁS INFORMACIÓN Y CLARIFICACIÓN CON ENRIQUE

Preprocesamiento

Primero haremos el preprocesamiento de los datos y un muestreo adecuado. Tres variables de la simulación son:

• N es el tamaño de la muestra
• smt es un factor de smoothing para la simulacion de la poblacion
• FactorMex es el factor estimado que lleva el ingreso anual en dolares en Estados unidos al ingreso mensual en pesos en México
• Recordar que solo se consideran edades entre 25 y 80 años e ingresos mensuales menores a (100K * FactorMex)

N <- 1e+05
smt <- 2
FactorMex <- 0.8

data <- read.csv("CombinaIncomeEdad.csv")

rango <- function(x, y, z) {
    if (z == 0) 
        return(min(which(x <= cumsum(y))))
    if (z > 0) 
        return(min(which(x <= cumsum(unlist(y[z, ])))))
    return(0)
}

aleapob <- runif(N, 1, sum(data$Total.INEGI))
pob <- data.frame(aleapob)
pob$sedad <- sapply(pob$aleapob, FUN = rango, y = data$Total.INEGI, z = 0)

pob$toting <- data$Total.US[pob$sedad]

pob$aleaing <- ceiling(runif(N, 1, pob$toting))
pob$sing <- mapply(FUN = rango, x = pob$aleaing, z = pob$sedad, MoreArgs = list(y = data[1:10, 
    7:47]))

pob$edad <- ceiling(sapply(22.5 + 5 * pob$sedad, FUN = rnorm, n = 1, sd = 5/smt))
pob$ingdol <- ceiling(sapply(-1250 + 2500 * (pob$sing - 1), FUN = rnorm, n = 1, 
    sd = 2500/smt))
pob$ingdol <- ceiling(pob$ingdol * FactorMex)
maxcero <- function(x) {
    return(max(x, 0))
}
pob$edad <- sapply(pob$edad, maxcero)
pob$ingdol <- sapply(pob$ingdol, maxcero)

pob$fedad <- data[pob$sedad, 1]
pob$fedadn <- 22.5 + 5 * (pob$sedad)
pob$fing <- sapply(-1250 + 2500 * (pob$sing - 1), maxcero)

Hasta este momento tenemos simplemente una tabla pob representando la población de clientes cuyas columnas importantes son:

• edad –la edad simulada
• ingdol – el ingreso mensual simulado
• fedad – el grupo de edad como factor
• fedadn – la edad media del grupo como factor
• fing – el ingreso medio del grupo como factor

Visualización

Esta es la distribucion por edad:

pob$fedadn <- 20 + 5 * (pob$sedad)
par(mfrow = c(1, 2))
hist(unlist(pob$fedadn), main = "Frecuencia x edad", xlab = "Edad", col = rgb(0.5, 
    0.5, 0), breaks = as.double(levels(factor(pob$fedadn))))
hist(pob$fing, main = "Frecuencia x ingreso", xlab = "Ingreso mensual $", col = rgb(0, 
    0.5, 0.5))

Estas son scatter plots de la poblacion en el espacio de Edad vs Ingreso:

par(mfrow = c(1, 2))
smoothScatter(pob$edad, pob$ingdol, ylab = "Ingreso mensual $", xlab = "Edad")

## KernSmooth 2.23 loaded
## Copyright M. P. Wand 1997-2009

plot(pob$edad, pob$ingdol, col = rgb(0.5, 0, 0.5, 1000/N), pch = 19, ylab = "Ingreso mensual $", 
    xlab = "Edad")

Ingreso por deciles:

deciles <- quantile(pob$ingdol, seq(0, 1, 0.1))
deciles

##    0%   10%   20%   30%   40%   50%   60%   70%   80%   90%  100% 
##     0   332  6336 11116 16296 21406 27016 33562 41690 54692 81984

Simular asignación de productos

Coche A , Coche B , Vida A , Vida B , Vida C

Ver archivo “Notas Simulacion de Productos.jpg” con diagramas y notas (o fórmulas subsiguientes) para significado

Asignación de seguro de coche

M <- FactorMex * 1e+05
# parametros seguro coche
Pen1 <- 0.9  ##penetracion coche
b1 <- 20000
h1 <- Pen1/(M - b1/2) * M
Pen2 <- 0.9  ##penetracion seguro de coche
b2 <- 20000
h2 <- Pen2/(M - b2/2) * M
ba <- 25000
bb <- 45000

ing <- pob$ingdol
tienecoche <- logical(length(ing))
tieneseguro <- logical(length(ing))
segurocoche <- character(length(ing))

for (i in 1:length(ing)) {
    tienecoche[i] <- runif(1) < (if (ing[i] < b1) 
        ing[i] * h1/b1 else h1)
    tieneseguro[i] <- tienecoche[i] & runif(1) < if (ing[i] < b2) 
        ing[i] * h2/b2 else h2
    if (tieneseguro[i]) {
        umb <- if (ing[i] < ba) {
            1
        } else if (ing[i] < bb) {
            1 - (ing[i] - ba)/(bb - ba)
        } else {
            0
        }
        if (runif(1) <= umb) 
            segurocoche[i] <- "Auto A" else segurocoche[i] <- "Auto B"
    } else segurocoche[i] <- NA
}

Asignación de seguro de vida

# parametros seguro vida
Pen3 <- 0.7  ##penetracion seguro vida
b3 <- 20000
b4 <- 30000
b5 <- 50
h3 <- Pen3/(M - b3/2) * M
h41 <- 1
h42 <- 0.3
h43 <- 0.9
h44 <- 0.3

ing <- pob$ingdol
tienesegvida <- logical(length(ing))
segurovida <- character(length(ing))

for (i in 1:length(ing)) {
    ingmod <- if (pob$edad[i] < b5) {
        ing[i]
    } else {
        ing[i] * (pob$edad[i] - b5)/(80 - b5) + ing[i]
    }
    tienesegvida[i] <- runif(1) < (if (ing[i] < b3) 
        ing[i] * h3/b3 else h3)
    if (tienesegvida[i]) {
        umb1 <- if (ingmod < b4) 
            h41 + ingmod * (h42 - h41)/b4 else h42 - (ingmod - b4) * h42/(M - b4)
        umb2 <- if (ingmod < b4) 
            1 + ingmod * (h43 - 1)/b4 else h43 - (ingmod - b4) * (h43 - h44)/(M - b4)
        random <- runif(1)
        if (random < umb1) 
            segurovida[i] <- "Vida A" else if (random < umb2) 
            segurovida[i] <- "Vida B" else segurovida[i] <- "Vida C"
    } else segurovida[i] <- NA
}

pob$tcoche <- tienecoche
pob$tsegcoche <- tieneseguro
pob$segcoche <- segurocoche
pob$tsegvida <- tienesegvida
pob$segvida <- segurovida

clientes <- data.frame(id = 1:length(pob$edad), edad = pob$edad, ing = pob$ingdol, 
    tcoche = pob$tcoche, tsegcoche = pob$tsegcoche, segcoche = factor(pob$segcoche), 
    tsegvida = pob$tsegvida, segvida = factor(pob$segvida))

Visualización de distribución de productos

ylim = N/20

summary(clientes)

##        id              edad         ing          tcoche       
##  Min.   :     1   Min.   :18   Min.   :    0   Mode :logical  
##  1st Qu.: 25001   1st Qu.:33   1st Qu.: 8698   FALSE:27308    
##  Median : 50000   Median :42   Median :21406   TRUE :72692    
##  Mean   : 50000   Mean   :44   Mean   :24932   NA's :0        
##  3rd Qu.: 75000   3rd Qu.:53   3rd Qu.:37732                  
##  Max.   :100000   Max.   :82   Max.   :81984                  
##  tsegcoche         segcoche      tsegvida         segvida     
##  Mode :logical   Auto A:37348   Mode :logical   Vida A:20107  
##  FALSE:33585     Auto B:29067   FALSE:42336     Vida B:25200  
##  TRUE :66415     NA's  :33585   TRUE :57664     Vida C:12357  
##  NA's :0                        NA's :0         NA's  :42336  
##                                                               
## 

par(mfrow = c(1, 2))

hist(clientes$ing[clientes$segcoche == "Auto A"], n = 10, xlim = c(1, 80000), 
    xlab = "Ingreso", ylim = c(0, ylim * 1.5), ylab = "Frec", main = "Auto A", 
    col = "blue")
hist(clientes$ing[clientes$segcoche == "Auto B"], n = 10, xlim = c(1, 80000), 
    xlab = "Ingreso", ylim = c(0, ylim * 1.5), ylab = "Frec", main = "Auto B", 
    col = "blue")

par(mfrow = c(2, 3))

hist(clientes$ing[clientes$segvida == "Vida A"], n = 10, xlim = c(1, 80000), 
    xlab = "Ingreso", ylim = c(0, ylim), ylab = "Frec", main = "Vida A", col = "orange")
hist(clientes$ing[clientes$segvida == "Vida B"], n = 10, xlim = c(1, 80000), 
    xlab = "Ingreso", ylim = c(0, ylim), ylab = "Frec", main = "Vida B", col = "orange")
hist(clientes$ing[clientes$segvida == "Vida C"], n = 10, xlim = c(1, 80000), 
    xlab = "Ingreso", ylim = c(0, ylim), ylab = "Frec", main = "Vida C", col = "orange")

hist(clientes$edad[clientes$segvida == "Vida A"], n = 10, xlim = c(1, 80), xlab = "Edad", 
    ylim = c(0, ylim), ylab = "Frec", main = "Vida A", col = "orange")
hist(clientes$edad[clientes$segvida == "Vida B"], n = 10, xlim = c(1, 80), xlab = "Edad", 
    ylim = c(0, ylim), ylab = "Frec", main = "Vida B", col = "orange")
hist(clientes$edad[clientes$segvida == "Vida C"], n = 10, xlim = c(1, 80), xlab = "Edad", 
    ylim = c(0, ylim), ylab = "Frec", main = "Vida C", col = "orange")

Exportar tabla

# write.csv(clientes,'clientes.csv')

Tablas segmentadas

Hacer segmentos por edad y cajitas

segedad <- numeric(length(clientes$edad))
seging <- numeric(length(clientes$ing))

pasoedad <- 5
pasoing <- 2500
for (i in 1:length(clientes$edad)) {
    edad <- clientes$edad[i]
    ing <- clientes$ing[i]
    segedad[i] <- edad - if ((edad%%pasoedad) < pasoedad/2) 
        edad%%pasoedad else edad%%pasoedad - pasoedad
    if (ing) 
        seging[i] <- ing - if ((ing%%pasoing) < pasoing/2) 
            ing%%pasoing else ing%%pasoing - pasoing else seging[i] <- 0

}

clientes$sedad <- segedad
clientes$sing <- seging

source("myImagePlot.R")

tca <- xtabs(segcoche == "Auto A" ~ sing + sedad, data = clientes)
myImagePlot(tca[dim(tca)[1]:2, ], title = "Auto A", ylab = "Ingreso $", xlab = "Edad")

tcb <- xtabs(segcoche == "Auto B" ~ sing + sedad, data = clientes)
myImagePlot(tcb[dim(tcb)[1]:2, ], title = "Auto B", ylab = "Ingreso $", xlab = "Edad")

tva <- xtabs(segvida == "Vida A" ~ sing + sedad, data = clientes)
myImagePlot(tva[dim(tva)[1]:2, ], title = "Vida A", ylab = "Ingreso $", xlab = "Edad")

tvb <- xtabs(segvida == "Vida B" ~ sing + sedad, data = clientes)
myImagePlot(tvb[dim(tvb)[1]:2, ], title = "Vida B", ylab = "Ingreso $", xlab = "Edad")

tvc <- xtabs(segvida == "Vida C" ~ sing + sedad, data = clientes)
myImagePlot(tvc[dim(tvc)[1]:2, ], title = "Vida C", ylab = "Ingreso $", xlab = "Edad")

tt <- table(clientes$sing, clientes$sedad)

Probabilidades por segmento

pca <- tca/(tt + 1)
myImagePlot(pca[dim(tca)[1]:2, ], title = "Prob Auto A", ylab = "Ingreso $", 
    xlab = "Edad")

pcb <- tcb/(tt + 1)
myImagePlot(pcb[dim(tcb)[1]:2, ], title = "Prob Auto B", ylab = "Ingreso $", 
    xlab = "Edad")

pva <- tva/(tt + 1)
myImagePlot(pva[dim(tva)[1]:2, ], title = "Prob Vida A", ylab = "Ingreso $", 
    xlab = "Edad")

pvb <- tvb/(tt + 1)
myImagePlot(pvb[dim(tvb)[1]:2, ], title = "Prob Vida B", ylab = "Ingreso $", 
    xlab = "Edad")

pvc <- tvc/(tt + 1)
myImagePlot(pvc[dim(tvc)[1]:2, ], title = "Prob Vida C", ylab = "Ingreso $", 
    xlab = "Edad")