A partir de la tabla de indicadores económicos para 2017:
Ajuste de la distribución a partir del método de máxima verosimilitud, reportar la estimación puntual de los parámetros y la distribución seleccionada.
Mediante el método Bootstrap, (a) calcule intervalo de confianza, (b) error estándar y (c) sesgo. Reporte en una tabla las cantidades para cada indicador.
A continuación se leen los datos, cada indicador se normaliza dividiendo en 100.
competitividad = read_excel("Competitividad.xlsx", sheet = "2017")
colnames(competitividad) <- c("departamento", "global", "fei", "il", "bsch", "cti", "igp", "region")
competitividad$global = competitividad$global/100
competitividad$fei = competitividad$fei /100
competitividad$il = competitividad$il /100
competitividad$bsch = competitividad$bsch /100
competitividad$cti = competitividad$cti /100
competitividad$igp = competitividad$igp /100
attach(competitividad)
summary(competitividad)
## departamento global fei il
## Length:32 Min. :0.1530 Min. :0.1750 Min. :0.0300
## Class :character 1st Qu.:0.3653 1st Qu.:0.3382 1st Qu.:0.4323
## Mode :character Median :0.4865 Median :0.4925 Median :0.5495
## Mean :0.4799 Mean :0.4863 Mean :0.5129
## 3rd Qu.:0.6038 3rd Qu.:0.6025 3rd Qu.:0.6610
## Max. :0.8750 Max. :0.9190 Max. :0.7860
## bsch cti igp region
## Min. :0.1510 Min. :0.0490 Min. :0.3300 Length:32
## 1st Qu.:0.4465 1st Qu.:0.1050 1st Qu.:0.5335 Class :character
## Median :0.6245 Median :0.1770 Median :0.6175 Mode :character
## Mean :0.5670 Mean :0.2553 Mean :0.6086
## 3rd Qu.:0.7163 3rd Qu.:0.3500 3rd Qu.:0.7010
## Max. :0.8730 Max. :0.9540 Max. :0.8360
En seguida se presentan las distribuciones ajustadas para cada indicador con sus correspondientes parámetros:
x = seq(0, 1, 0.01)
fitdist2 = function(v, d, m='mle', s=NULL, f=NULL) fitdist(v, d, m, s, f)$estimate
par(mfrow = c(3,2))
p = fitdist2(global, 'norm')
titulo = paste('dist = normal, media =',round(p[1],2), 'std =', round(p[2],2))
hist(global, freq = F, main=titulo); lines(x, dnorm(x, p[1], p[2]))
p = fitdist2(fei, 'gamma')
titulo = paste('dist = gamma, media =',round(p[1]/p[2],2), 'gamma =', round(p[2],2))
hist(fei, freq = F, main=titulo); lines(x, dgamma(x, p[1], p[2]))
p = fitdist2(il, 'beta')
titulo = paste('dist = beta, media =',round(p[1]/(p[1]+p[2]),2), 'beta =', round(p[2],2))
hist(il, freq = F, main=titulo);lines(x, dbeta(x,p[1], p[2]))
p = fitdist2(bsch, 'beta')
titulo = paste('dist = beta, media =',round(p[1]/(p[1]+p[2]),2), 'beta =', round(p[2],2))
hist(bsch, freq = F, main=titulo);lines(x, dbeta(x,p[1], p[2]))
p = fitdist2(cti, 'gamma')
titulo = paste('dist = gamma, media =',round(p[1]/p[2],2), 'gamma =', round(p[2],2))
hist(cti, freq = F, main=titulo);lines(x, dgamma(x,p[1], p[2]))
p = fitdist2(igp, 'beta')
titulo = paste('dist = beta, media =',round(p[1]/(p[1]+p[2]),2), 'beta =', round(p[2],2))
hist(igp, freq = F, main=titulo);lines(x, dbeta(x,p[1], p[2]))
Ahora se realizan las estimaciones Bootstrap para cada media.
est = function(est, variable) {
b.est = boot(variable, est, 1000, sim = 'ordinary')
t0 = b.est$t0
error = sd(b.est$t)
sesgo = b.est$t0 - mean(b.est$t)
inter = c(boot.ci(b.est, type = 'perc')[4]$percent[4],
boot.ci(b.est, type = 'perc')[4]$percent[5])
return(list('t0'=t0, 'error'=error, 'sesgo'=sesgo, 'ci'=inter))
}
est_norm_mean = function(x, i) {fitdist(x[i], 'norm' )$estimate[1]}
est_beta_mean = function(x, i) {p = fitdist(x[i], 'beta' )$estimate; p[1]/(p[1]+p[2])}
est_gamma_mean = function(x, i) {p = fitdist(x[i], 'gamma')$estimate; p[1]/p[2]}
est_global = est(est_norm_mean, global)
est_fei = est(est_gamma_mean, fei)
est_il = est(est_beta_mean, il)
est_bsch = est(est_beta_mean, bsch)
est_cti = est(est_gamma_mean, cti)
est_igp = est(est_beta_mean, igp)
media = c(round(est_global$t0, 3),
round(est_fei$t0, 3),
round(est_il$t0, 3),
round(est_bsch$t0, 3),
round(est_cti$t0, 3),
round(est_igp$t0, 3))
errores = c(round(est_global$error, 3),
round(est_fei$error, 3),
round(est_il$error, 3),
round(est_bsch$error, 3),
round(est_cti$error, 3),
round(est_igp$error, 3))
sesgos = c(round(est_global$sesgo, 5),
round(est_fei$sesgo, 5),
round(est_il$sesgo, 5),
round(est_bsch$sesgo, 5),
round(est_cti$sesgo, 5),
round(est_igp$sesgo, 5))
interv = c(paste(round(est_global$ci[1],3),'-',round(est_global$ci[2],3)),
paste(round(est_fei$ci[1], 3),'-',round(est_fei$ci[2], 3)),
paste(round(est_il$ci[1], 3),'-',round(est_il$ci[2], 3)),
paste(round(est_bsch$ci[1], 3),'-',round(est_bsch$ci[2], 3)),
paste(round(est_cti$ci[1], 3),'-',round(est_cti$ci[2], 3)),
paste(round(est_igp$ci[1], 3),'-',round(est_igp$ci[2], 3)))
data = data.frame('media' = media,
'error_std' = errores,
'sesgo' = sesgos,
'interv_95' = interv)
rownames(data) = c('global', 'fei', 'il', 'bsch', 'cti', 'igp')
kable(data)
| media | error_std | sesgo | interv_95 | |
|---|---|---|---|---|
| global | 0.480 | 0.032 | 0.00030 | 0.42 - 0.545 |
| fei | 0.486 | 0.032 | -0.00076 | 0.425 - 0.55 |
| il | 0.493 | 0.037 | -0.00056 | 0.418 - 0.564 |
| bsch | 0.558 | 0.036 | 0.00066 | 0.486 - 0.625 |
| cti | 0.255 | 0.034 | 0.00042 | 0.192 - 0.329 |
| igp | 0.608 | 0.021 | 0.00036 | 0.568 - 0.649 |
Y en seguida para los parámetros de escala:
| escala | error_std | sesgo | interv_95 | |
|---|---|---|---|---|
| global_std | 0.174 | 0.018 | 0.00426 | 0.135 - 0.205 |
| fei_gamma | 14.041 | 3.705 | -1.03045 | 9.743 - 25.282 |
| il_beta | 2.362 | 1.067 | -0.37288 | 1.619 - 5.694 |
| bsch_beta | 2.140 | 0.484 | -0.17276 | 1.705 - 3.355 |
| cti_gamma | 7.471 | 2.397 | -0.79710 | 4.996 - 14.715 |
| igp_beta | 6.297 | 1.707 | -0.52347 | 4.499 - 11.334 |