In this document, I replicate Tables 1, 2, and 3 in Coate & Ravallion (1993). There is one very minor typo in the paper: the equation (5) should be \(\theta_{ij}^*(t + 1) = \min \{ \widehat{\theta}_{ij}, y_i - [y_i^{1 - \rho} - (v(\theta^*(t)) - \bar{v}) (1 - \rho) / r ]^{1 / (1 - \rho)} \}\) instead of \(\theta_{ij}^*(t + 1) = \min \{ \widehat{\theta}_{ij}, y_i - [y_i^{1 - \rho} - (v(\theta^*(t)) - \bar{v}) / r ]^{1 / (1 - \rho)} \}\). This is fixed in the code below.
rm(list = ls())
library(knitr)
# Common settings ------------------
y_vals <- c(1, 2, 3)
theta_hat_vals <- c((y_vals[2] - y_vals[1]) / 2,
(y_vals[3] - y_vals[1]) / 2,
(y_vals[3] - y_vals[2]) / 2)
# v(theta)
v_theta <- function(theta, y_vals, pi, rho){
output <- 0
for (i in 1:3){
if (i > 1){
for (j in 1:(i-1)){
output <- output + pi[i,j] * ((y_vals[i] - theta[i+j-2]) ^ (1 - rho) / (1 - rho) +
(y_vals[j] + theta[i+j-2]) ^ (1 - rho) / (1 - rho))
}
} else {
output <- output
}
output <- output + pi[i, i] * (y_vals[i]) ^ (1 - rho) / (1 - rho)
}
return(output)
}
# Iteration for theta_star
iteration <- function(theta_star, theta_hat, y_val, y_vals, pi, rho, r, v_bar){
new_theta_star_ij <- min(theta_hat,
y_val - (y_val ^ (1 - rho) -
(v_theta(theta_star, y_vals, pi, rho) - v_bar)
* (1 - rho) / r) ^ (1 / (1 - rho)))
return(new_theta_star_ij)
}
# Iteration loop for theta_star
theta_star_loop <- function(init_theta_star, theta_hat_vals, y_vals, pi, rho, r, v_bar){
iter <- 0
tol <- 1e-16
error <- tol + 1
theta_star <- init_theta_star
new_theta_star <- rep(0, 3)
while (iter < 500 & error > tol){
new_theta_star[1] <- iteration(theta_star, theta_hat_vals[1],
y_vals[2], y_vals, pi, rho, r, v_bar)
new_theta_star[2] <- iteration(theta_star, theta_hat_vals[2],
y_vals[3], y_vals, pi, rho, r, v_bar)
new_theta_star[3] <- iteration(theta_star, theta_hat_vals[3],
y_vals[3], y_vals, pi, rho, r, v_bar)
error = max(abs(new_theta_star - theta_star))
theta_star <- new_theta_star
iter <- iter + 1
}
return(theta_star)
}
Replication of Table 1: The “highly covariate” income stream
table_1 <- matrix(rep(0, (5 + 11 + 14) * 6), ncol = 6)
pi <- matrix(c(0.2, 0.05, 0.05, 0.05, 0.3, 0.05, 0.05, 0.05, 0.2), nrow = 3)
r_vals <- c(0.05, 0.15, 0.25)
rho_vals <- c(0.1, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5, 1.7, 1.9, 2.1, 2.3, 2.5, 2.7)
# r = 0.05
r <- r_vals[1]
for (i in 1:5){
# r
table_1[i, 1] <- r
# rho
table_1[i, 2] <- rho_vals[i]
# v_bar
v_bar <- 0
for (j in 1:3){
for (k in 1:3){
v_bar <- v_bar + pi[j,k] * (y_vals[j] ^ (1 - rho_vals[i]) / (1 - rho_vals[i]))
}
}
# theta_star iteration
init_theta_star <- theta_hat_vals
theta_star <- theta_star_loop(init_theta_star, theta_hat_vals,
y_vals, pi, rho_vals[i], r, v_bar)
# theta_star_21
table_1[i, 3] <- theta_star[1]
# theta_star_31
table_1[i, 4] <- theta_star[2]
# theta_star_32
table_1[i, 5] <- theta_star[3]
# gamma
table_1[i, 6] <- (v_theta(theta_star, y_vals, pi, rho_vals[i]) - v_bar) /
(v_theta(theta_hat_vals, y_vals, pi, rho_vals[i]) - v_bar)
}
# r = 0.15
r <- r_vals[2]
for (i in 1:11){
# r
table_1[i+5, 1] <- r
# rho
table_1[i+5, 2] <- rho_vals[i]
# v_bar
v_bar <- 0
for (j in 1:3){
for (k in 1:3){
v_bar <- v_bar + pi[j,k] * (y_vals[j] ^ (1 - rho_vals[i]) / (1 - rho_vals[i]))
}
}
# theta_star iteration
init_theta_star <- theta_hat_vals
theta_star <- theta_star_loop(init_theta_star, theta_hat_vals,
y_vals, pi, rho_vals[i], r, v_bar)
# theta_star_21
table_1[i+5, 3] <- theta_star[1]
# theta_star_31
table_1[i+5, 4] <- theta_star[2]
# theta_star_32
table_1[i+5, 5] <- theta_star[3]
# gamma
table_1[i+5, 6] <- (v_theta(theta_star, y_vals, pi, rho_vals[i]) - v_bar) /
(v_theta(theta_hat_vals, y_vals, pi, rho_vals[i]) - v_bar)
}
# r = 0.25
r <- r_vals[3]
for (i in 1:14){
# r
table_1[i+16, 1] <- r
# rho
table_1[i+16, 2] <- rho_vals[i]
# v_bar
v_bar <- 0
for (j in 1:3){
for (k in 1:3){
v_bar <- v_bar + pi[j,k] * (y_vals[j] ^ (1 - rho_vals[i]) / (1 - rho_vals[i]))
}
}
# theta_star iteration
init_theta_star <- theta_hat_vals
theta_star <- theta_star_loop(init_theta_star, theta_hat_vals,
y_vals, pi, rho_vals[i], r, v_bar)
# theta_star_21
table_1[i+16, 3] <- theta_star[1]
# theta_star_31
table_1[i+16, 4] <- theta_star[2]
# theta_star_32
table_1[i+16, 5] <- theta_star[3]
# gamma
table_1[i+16, 6] <- (v_theta(theta_star, y_vals, pi, rho_vals[i]) - v_bar) /
(v_theta(theta_hat_vals, y_vals, pi, rho_vals[i]) - v_bar)
}
colnames(table_1) <- c("$r$", "$\\rho$", "$\\theta_{21}^*$",
"$\\theta_{31}^*$", "$\\theta_{32}^*$", "$\\gamma$")
kable(table_1, digits = 3, caption = "Highly covariate income streams")
Highly covariate income streams
| 0.05 |
0.1 |
0.000 |
0.000 |
0.000 |
0.000 |
| 0.05 |
0.3 |
0.000 |
0.000 |
0.000 |
0.000 |
| 0.05 |
0.5 |
0.262 |
0.323 |
0.323 |
0.649 |
| 0.05 |
0.7 |
0.500 |
0.686 |
0.500 |
0.942 |
| 0.05 |
0.9 |
0.500 |
0.963 |
0.500 |
0.999 |
| 0.15 |
0.1 |
0.000 |
0.000 |
0.000 |
0.000 |
| 0.15 |
0.3 |
0.000 |
0.000 |
0.000 |
0.000 |
| 0.15 |
0.5 |
0.000 |
0.000 |
0.000 |
0.000 |
| 0.15 |
0.7 |
0.000 |
0.000 |
0.000 |
0.000 |
| 0.15 |
0.9 |
0.001 |
0.001 |
0.001 |
0.003 |
| 0.15 |
1.1 |
0.179 |
0.279 |
0.279 |
0.568 |
| 0.15 |
1.3 |
0.299 |
0.500 |
0.500 |
0.819 |
| 0.15 |
1.5 |
0.386 |
0.685 |
0.500 |
0.934 |
| 0.15 |
1.7 |
0.452 |
0.846 |
0.500 |
0.986 |
| 0.15 |
1.9 |
0.500 |
0.982 |
0.500 |
1.000 |
| 0.15 |
2.1 |
0.500 |
1.000 |
0.500 |
1.000 |
| 0.25 |
0.1 |
0.000 |
0.000 |
0.000 |
0.000 |
| 0.25 |
0.3 |
0.000 |
0.000 |
0.000 |
0.000 |
| 0.25 |
0.5 |
0.000 |
0.000 |
0.000 |
0.000 |
| 0.25 |
0.7 |
0.000 |
0.000 |
0.000 |
0.000 |
| 0.25 |
0.9 |
0.000 |
0.000 |
0.000 |
0.000 |
| 0.25 |
1.1 |
0.000 |
0.000 |
0.000 |
0.000 |
| 0.25 |
1.3 |
0.034 |
0.057 |
0.057 |
0.140 |
| 0.25 |
1.5 |
0.149 |
0.270 |
0.270 |
0.543 |
| 0.25 |
1.7 |
0.234 |
0.451 |
0.451 |
0.761 |
| 0.25 |
1.9 |
0.299 |
0.608 |
0.500 |
0.879 |
| 0.25 |
2.1 |
0.353 |
0.750 |
0.500 |
0.946 |
| 0.25 |
2.3 |
0.398 |
0.875 |
0.500 |
0.981 |
| 0.25 |
2.5 |
0.435 |
0.985 |
0.500 |
0.995 |
| 0.25 |
2.7 |
0.468 |
1.000 |
0.500 |
0.999 |
Replication of Table 2: The “moderately covariate” income stream
table_2 <- matrix(rep(0, (5 + 7 + 9) * 6), ncol = 6)
pi <- matrix(c(0.1, 0.1, 0.1, 0.1, 0.2, 0.1, 0.1, 0.1, 0.1), nrow = 3)
r_vals <- c(0.05, 0.15, 0.25)
rho_vals <- c(0.1, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5, 1.7, 1.9, 2.1, 2.3, 2.5, 2.7)
# r = 0.05
r <- r_vals[1]
for (i in 1:5){
# r
table_2[i, 1] <- r
# rho
table_2[i, 2] <- rho_vals[i]
# v_bar
v_bar <- 0
for (j in 1:3){
for (k in 1:3){
v_bar <- v_bar + pi[j,k] * (y_vals[j] ^ (1 - rho_vals[i]) / (1 - rho_vals[i]))
}
}
# theta_star iteration
init_theta_star <- theta_hat_vals
theta_star <- theta_star_loop(init_theta_star, theta_hat_vals,
y_vals, pi, rho_vals[i], r, v_bar)
# theta_star_21
table_2[i, 3] <- theta_star[1]
# theta_star_31
table_2[i, 4] <- theta_star[2]
# theta_star_32
table_2[i, 5] <- theta_star[3]
# gamma
table_2[i, 6] <- (v_theta(theta_star, y_vals, pi, rho_vals[i]) - v_bar) /
(v_theta(theta_hat_vals, y_vals, pi, rho_vals[i]) - v_bar)
}
# r = 0.15
r <- r_vals[2]
for (i in 1:7){
# r
table_2[i+5, 1] <- r
# rho
table_2[i+5, 2] <- rho_vals[i]
# v_bar
v_bar <- 0
for (j in 1:3){
for (k in 1:3){
v_bar <- v_bar + pi[j,k] * (y_vals[j] ^ (1 - rho_vals[i]) / (1 - rho_vals[i]))
}
}
# theta_star iteration
init_theta_star <- theta_hat_vals
theta_star <- theta_star_loop(init_theta_star, theta_hat_vals,
y_vals, pi, rho_vals[i], r, v_bar)
# theta_star_21
table_2[i+5, 3] <- theta_star[1]
# theta_star_31
table_2[i+5, 4] <- theta_star[2]
# theta_star_32
table_2[i+5, 5] <- theta_star[3]
# gamma
table_2[i+5, 6] <- (v_theta(theta_star, y_vals, pi, rho_vals[i]) - v_bar) /
(v_theta(theta_hat_vals, y_vals, pi, rho_vals[i]) - v_bar)
}
# r = 0.25
r <- r_vals[3]
for (i in 1:9){
# r
table_2[i+12, 1] <- r
# rho
table_2[i+12, 2] <- rho_vals[i]
# v_bar
v_bar <- 0
for (j in 1:3){
for (k in 1:3){
v_bar <- v_bar + pi[j,k] * (y_vals[j] ^ (1 - rho_vals[i]) / (1 - rho_vals[i]))
}
}
# theta_star iteration
init_theta_star <- theta_hat_vals
theta_star <- theta_star_loop(init_theta_star, theta_hat_vals,
y_vals, pi, rho_vals[i], r, v_bar)
# theta_star_21
table_2[i+12, 3] <- theta_star[1]
# theta_star_31
table_2[i+12, 4] <- theta_star[2]
# theta_star_32
table_2[i+12, 5] <- theta_star[3]
# gamma
table_2[i+12, 6] <- (v_theta(theta_star, y_vals, pi, rho_vals[i]) - v_bar) /
(v_theta(theta_hat_vals, y_vals, pi, rho_vals[i]) - v_bar)
}
colnames(table_2) <- c("$r$", "$\\rho$", "$\\theta_{21}^*$",
"$\\theta_{31}^*$", "$\\theta_{32}^*$", "$\\gamma$")
kable(table_2, digits = 3, caption = "Moderately covariate income streams")
Moderately covariate income streams
| 0.05 |
0.1 |
0.000 |
0.000 |
0.000 |
0.000 |
| 0.05 |
0.3 |
0.365 |
0.415 |
0.415 |
0.766 |
| 0.05 |
0.5 |
0.500 |
0.934 |
0.500 |
0.997 |
| 0.05 |
0.7 |
0.500 |
1.000 |
0.500 |
1.000 |
| 0.05 |
0.9 |
0.500 |
1.000 |
0.500 |
1.000 |
| 0.15 |
0.1 |
0.000 |
0.000 |
0.000 |
0.000 |
| 0.15 |
0.3 |
0.000 |
0.000 |
0.000 |
0.000 |
| 0.15 |
0.5 |
0.000 |
0.000 |
0.000 |
0.000 |
| 0.15 |
0.7 |
0.253 |
0.338 |
0.338 |
0.663 |
| 0.15 |
0.9 |
0.434 |
0.628 |
0.500 |
0.916 |
| 0.15 |
1.1 |
0.500 |
0.857 |
0.500 |
0.989 |
| 0.15 |
1.3 |
0.500 |
1.000 |
0.500 |
1.000 |
| 0.25 |
0.1 |
0.000 |
0.000 |
0.000 |
0.000 |
| 0.25 |
0.3 |
0.000 |
0.000 |
0.000 |
0.000 |
| 0.25 |
0.5 |
0.000 |
0.000 |
0.000 |
0.000 |
| 0.25 |
0.7 |
0.000 |
0.000 |
0.000 |
0.000 |
| 0.25 |
0.9 |
0.122 |
0.176 |
0.176 |
0.397 |
| 0.25 |
1.1 |
0.282 |
0.440 |
0.440 |
0.769 |
| 0.25 |
1.3 |
0.390 |
0.649 |
0.500 |
0.920 |
| 0.25 |
1.5 |
0.469 |
0.827 |
0.500 |
0.983 |
| 0.25 |
1.7 |
0.500 |
0.977 |
0.500 |
1.000 |
Replication of Table 3: The “highly non-covariate” income stream
table_3 <- matrix(rep(0, (5 + 6 + 8) * 6), ncol = 6)
pi <- matrix(c(0.1, 0.1, 0.15, 0.1, 0.1, 0.1, 0.15, 0.1, 0.1), nrow = 3)
r_vals <- c(0.05, 0.15, 0.25)
rho_vals <- c(0.1, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5, 1.7, 1.9, 2.1, 2.3, 2.5, 2.7)
# r = 0.05
r <- r_vals[1]
for (i in 1:5){
# r
table_3[i, 1] <- r
# rho
table_3[i, 2] <- rho_vals[i]
# v_bar
v_bar <- 0
for (j in 1:3){
for (k in 1:3){
v_bar <- v_bar + pi[j,k] * (y_vals[j] ^ (1 - rho_vals[i]) / (1 - rho_vals[i]))
}
}
# theta_star iteration
init_theta_star <- theta_hat_vals
theta_star <- theta_star_loop(init_theta_star, theta_hat_vals,
y_vals, pi, rho_vals[i], r, v_bar)
# theta_star_21
table_3[i, 3] <- theta_star[1]
# theta_star_31
table_3[i, 4] <- theta_star[2]
# theta_star_32
table_3[i, 5] <- theta_star[3]
# gamma
table_3[i, 6] <- (v_theta(theta_star, y_vals, pi, rho_vals[i]) - v_bar) /
(v_theta(theta_hat_vals, y_vals, pi, rho_vals[i]) - v_bar)
}
# r = 0.15
r <- r_vals[2]
for (i in 1:6){
# r
table_3[i+5, 1] <- r
# rho
table_3[i+5, 2] <- rho_vals[i]
# v_bar
v_bar <- 0
for (j in 1:3){
for (k in 1:3){
v_bar <- v_bar + pi[j,k] * (y_vals[j] ^ (1 - rho_vals[i]) / (1 - rho_vals[i]))
}
}
# theta_star iteration
init_theta_star <- theta_hat_vals
theta_star <- theta_star_loop(init_theta_star, theta_hat_vals,
y_vals, pi, rho_vals[i], r, v_bar)
# theta_star_21
table_3[i+5, 3] <- theta_star[1]
# theta_star_31
table_3[i+5, 4] <- theta_star[2]
# theta_star_32
table_3[i+5, 5] <- theta_star[3]
# gamma
table_3[i+5, 6] <- (v_theta(theta_star, y_vals, pi, rho_vals[i]) - v_bar) /
(v_theta(theta_hat_vals, y_vals, pi, rho_vals[i]) - v_bar)
}
# r = 0.25
r <- r_vals[3]
for (i in 1:8){
# r
table_3[i+11, 1] <- r
# rho
table_3[i+11, 2] <- rho_vals[i]
# v_bar
v_bar <- 0
for (j in 1:3){
for (k in 1:3){
v_bar <- v_bar + pi[j,k] * (y_vals[j] ^ (1 - rho_vals[i]) / (1 - rho_vals[i]))
}
}
# theta_star iteration
init_theta_star <- theta_hat_vals
theta_star <- theta_star_loop(init_theta_star, theta_hat_vals,
y_vals, pi, rho_vals[i], r, v_bar)
# theta_star_21
table_3[i+11, 3] <- theta_star[1]
# theta_star_31
table_3[i+11, 4] <- theta_star[2]
# theta_star_32
table_3[i+11, 5] <- theta_star[3]
# gamma
table_3[i+11, 6] <- (v_theta(theta_star, y_vals, pi, rho_vals[i]) - v_bar) /
(v_theta(theta_hat_vals, y_vals, pi, rho_vals[i]) - v_bar)
}
colnames(table_3) <- c("$r$", "$\\rho$", "$\\theta_{21}^*$",
"$\\theta_{31}^*$", "$\\theta_{32}^*$", "$\\gamma$")
kable(table_3, digits = 3, caption = "Highly non-covariate income streams")
Highly non-covariate income streams
| 0.05 |
0.1 |
0.000 |
0.000 |
0.000 |
0.000 |
| 0.05 |
0.3 |
0.500 |
0.648 |
0.500 |
0.914 |
| 0.05 |
0.5 |
0.500 |
1.000 |
0.500 |
1.000 |
| 0.05 |
0.7 |
0.500 |
1.000 |
0.500 |
1.000 |
| 0.05 |
0.9 |
0.500 |
1.000 |
0.500 |
1.000 |
| 0.15 |
0.1 |
0.000 |
0.000 |
0.000 |
0.000 |
| 0.15 |
0.3 |
0.000 |
0.000 |
0.000 |
0.000 |
| 0.15 |
0.5 |
0.131 |
0.161 |
0.161 |
0.360 |
| 0.15 |
0.7 |
0.425 |
0.570 |
0.500 |
0.872 |
| 0.15 |
0.9 |
0.500 |
0.857 |
0.500 |
0.987 |
| 0.15 |
1.1 |
0.500 |
1.000 |
0.500 |
1.000 |
| 0.25 |
0.1 |
0.000 |
0.000 |
0.000 |
0.000 |
| 0.25 |
0.3 |
0.000 |
0.000 |
0.000 |
0.000 |
| 0.25 |
0.5 |
0.000 |
0.000 |
0.000 |
0.000 |
| 0.25 |
0.7 |
0.058 |
0.077 |
0.077 |
0.184 |
| 0.25 |
0.9 |
0.287 |
0.415 |
0.415 |
0.734 |
| 0.25 |
1.1 |
0.428 |
0.665 |
0.500 |
0.924 |
| 0.25 |
1.3 |
0.500 |
0.866 |
0.500 |
0.989 |
| 0.25 |
1.5 |
0.500 |
1.000 |
0.500 |
1.000 |