Edgeworth Box 以前に描いたコードを整理しなおしたと言うか、手計算で微分 して均衡や契約曲線を求めていたのでどちらかと言うと作り直したものですが、applyやwhichやexpressionなどに慣れるのに良さそうだったので。数値解析で均衡点を求めているので、初期財の配分を変えたり、効用関数を他の凸で連続なモノ（e.g. CES型）に変えても、そこの部分だけの変更でプロットできるはずです。
loadLib <- function ( ...){
for ( lname in unlist ( list ( ...))){
if ( ! any ( suppressWarnings ( library ( quietly= TRUE , verbose= FALSE )$ results[, "Package" ] == lname))){
stop ( sprintf ( "Do install.packages( \" %s \" ) before runnning this script." , lname))
}
library ( lname, character.only = TRUE )
}
}
loadLib ( "nleqslv" )
initial_goods <- matrix ( c ( 0.9 , 0.1 , 0.3 , 0.7 ), 2 , 2 ,
dimnames = list ( c ( "person 1" , "person 2" ), c ( "g/s A" , "g/s B" )))
U_pt <- function ( b, gs_a, gs_b){
if ( gs_a <= 0 || gs_b<= 0 ) return ( - Inf )
as.numeric ( b[ 1 ] * log ( gs_a) + b[ 2 ] * log ( gs_b) + b[ 3 ] * log ( gs_a) ^ 2 + b[ 4 ] * log ( gs_b) ^ 2 + b[ 5 ] * log ( gs_a) * log ( gs_b))
}
U_1 <- function ( gs_a, gs_b) U_pt ( c ( 0.4 , 0.6 , 0 , - 0.5 , 0 ), gs_a, gs_b)
U_2 <- function ( gs_a, gs_b) U_pt ( c ( 0.5 , 0.5 , 0 , - 0.4 , 0 ), gs_a, gs_b)
sum_of_goods <- apply ( initial_goods, 2 , sum)
dU <- function ( U, gs_a, gs_b, h= 1e-4 ){
r <- c ( U ( gs_a + h, gs_b) - U ( gs_a - h, gs_b),
U ( gs_a, gs_b + h) - U ( gs_a, gs_b - h)) / ( 2 * h)
names ( r) <- c ( "d(g/s A)" , "d(g/s B)" )
r
}
getPrice <- function ( goods){
dU1 <- dU ( U_1, goods[ "person 1" , "g/s A" ], goods[ "person 1" , "g/s B" ])
dU2 <- dU ( U_2, goods[ "person 2" , "g/s A" ], goods[ "person 2" , "g/s B" ])
r <- c ( dU1[ "d(g/s A)" ] / dU1[ "d(g/s B)" ], dU2[ "d(g/s A)" ] / dU2[ "d(g/s B)" ])
names ( r) <- c ( "person 1" , "person 2" )
r
}
getDistribution <- function ( initial_goods, exchange){
goods <- initial_goods
goods[ "person 1" , ] <- initial_goods[ "person 1" , ] - exchange
goods[ "person 2" , ] <- initial_goods[ "person 2" , ] + exchange
goods
}
getEquilibrium <- function ( initial_goods){
r_nleqslv <- nleqslv ( c ( 0.0 , - 0.0 ), function ( exchange){
goods <- getDistribution ( initial_goods, exchange)
p <- getPrice ( goods)
c ( p[ "person 1" ] - p[ "person 2" ], exchange[ 1 ] * p[ "person 1" ] + exchange[ 2 ])
})
exchange <- r_nleqslv$ x
goods <- getDistribution ( initial_goods, exchange)
price <- getPrice ( goods)
shpe <- substitute ( B0 - a* ( A - A0), list ( B0= goods[ 1 , 2 ], a= price[ 1 ], A0= goods[ 1 , 1 ]))
list ( price= price, exchange= exchange, goods= goods, shp= function ( A){ eval ( shpe) })
}
equilibrium <- getEquilibrium ( initial_goods)
drawIDC <- function ( goods, lwd= c ( 1.5 , 1.5 ), col= c ( "blue" , "red" , "pink" ), names= NULL , IsCore= FALSE , density= 10 ){
UL_1 <- U_1 ( goods[ "person 1" , "g/s A" ], goods[ "person 1" , "g/s B" ])
UL_2 <- U_2 ( goods[ "person 2" , "g/s A" ], goods[ "person 2" , "g/s B" ])
demand_B <- function ( UL, U, A){
B <- numeric ( length ( A))
for ( i in 1 : length ( A)){
r_nlm <- nlm ( function ( b){
if ( b <= 0 ){
return ( 1e+6 )
}
( UL - U ( A[ i], b)) ^ 2
}, 1 )
B[ i] <- with ( r_nlm, ifelse ( code<= 3 , estimate[ 1 ], NA ))
}
B
}
scale <- 1.05
A <- seq ( 1e-6 , sum_of_goods[ "g/s A" ] * scale, length.out= 101 )
IC_1 <- matrix ( c ( A, demand_B ( UL_1, U_1, A)),
length ( A), 2 )
colnames ( IC_1) <- c ( "g/s A" , "g/s B" )
IC_2 <- matrix ( c ( sum_of_goods[ "g/s A" ] - A,
sum_of_goods[ "g/s B" ] - demand_B ( UL_2, U_2, A)),
length ( A), 2 )
colnames ( IC_2) <- c ( "g/s A" , "g/s B" )
lines ( IC_1[, "g/s B" ] ~ IC_1[, "g/s A" ], lwd= lwd[ 1 ], col= col[ 1 ])
lines ( IC_2[, "g/s B" ] ~ IC_2[, "g/s A" ], lwd= lwd[ 2 ], col= col[ 2 ])
if ( ! is.null ( names)){
xadj <- ( par ()$ usr[ 2 ] - par ()$ usr[ 1 ]) / 75
b <- max ( IC_1[ , "g/s B" ][ IC_1[ , "g/s B" ] < sum_of_goods[ "g/s B" ]], na.rm = TRUE )
i <- which ( IC_1[ , "g/s B" ] == b, arr.ind= TRUE )
text ( IC_1[ i, "g/s A" ] + xadj, b, names[ 1 ], col= col[ 1 ], adj= c ( 0 , 1 ))
b <- min ( IC_2[ , "g/s B" ][ IC_2[ , "g/s B" ] > 0 ], na.rm = TRUE )
i <- which ( IC_2[ , "g/s B" ] == b, arr.ind= TRUE )
text ( IC_2[ i, "g/s A" ] - xadj, b, names[ 2 ], col= col[ 2 ], adj= c ( 1 , 0 ))
}
if ( IsCore){
af_1 <- approxfun ( IC_1[ , "g/s A" ], IC_1[ , "g/s B" ])
af_2 <- approxfun ( IC_2[ , "g/s A" ], IC_2[ , "g/s B" ])
A <- seq ( 1e-6 , sum_of_goods[ "g/s A" ] - 1e-6 , length.out= 100 )
A <- A[ af_1 ( A) <= af_2 ( A)]
x <- c ( A, rev ( A))
y <- c ( af_1 ( A), af_2 ( rev ( A)))
polygon ( x, y, density= density, col= col[ 3 ], border= NA )
}
}
par ( bg= "white" , mar= c ( 0 , 0 , 0 , 0 ), cex= 1 )
plot.new ()
margin <- 0.1
xlim <- c ( - margin, ( 1 + margin)) * sum_of_goods[ "g/s A" ]
ylim <- c ( - margin, ( 1 + margin)) * sum_of_goods[ "g/s B" ]
plot.window ( xlim= xlim, ylim= ylim)
col <- c ( "blue" , "red" )
with ( list ( l= 0.125 , lwd= 2 ), {
xlim <- c ( - l/ 2 , ( 1.0 + l/ 2 )) * sum_of_goods[ "g/s A" ]
ylim <- c ( - l/ 2 , ( 1.0 + l/ 2 )) * sum_of_goods[ "g/s B" ]
arrows ( 0.0 , 0.0 , xlim[ 2 ], 0 , col= col[ 1 ], lwd= lwd, length= l)
arrows ( 0.0 , 0.0 , 0 , ylim[ 2 ], col= col[ 1 ], lwd= lwd, length= l)
arrows ( sum_of_goods[ "g/s A" ], sum_of_goods[ "g/s B" ], sum_of_goods[ "g/s A" ], ylim[ 1 ], col= col[ 2 ], lwd= lwd, length= l)
arrows ( sum_of_goods[ "g/s A" ], sum_of_goods[ "g/s B" ], xlim[ 1 ], sum_of_goods[ "g/s B" ], col= col[ 2 ], lwd= lwd, length= l)
text ( 0 , 0 , expression ( O[ 1 ]), col= col[ 1 ], adj= c ( 1 , 1 ))
text ( xlim[ 2 ], 0 , expression ( G[ 1 ] ^ A), col= col[ 1 ], adj= c ( 0 , 1 ))
text ( 0 , ylim[ 2 ], expression ( G[ 1 ] ^ B), col= col[ 1 ], adj= c ( 1 , 0 ))
text ( sum_of_goods[ "g/s A" ], sum_of_goods[ "g/s B" ], expression ( O[ 2 ]), col= col[ 2 ], adj= c ( 0 , 0 ))
text ( xlim[ 1 ], sum_of_goods[ "g/s B" ], expression ( G[ 2 ] ^ A), col= col[ 2 ], adj= c ( 1 , 0 ))
text ( sum_of_goods[ "g/s A" ], ylim[ 1 ], expression ( G[ 2 ] ^ B), col= col[ 2 ], adj= c ( 0 , 1 ))
})
grid <- expand.grid (
"g/s A" = seq ( 1e-2 , 1 - 1e-2 , length.out= 11 ) * sum_of_goods[ "g/s A" ],
"g/s B" = seq ( 1e-2 , 1 - 1e-2 , length.out= 11 ) * sum_of_goods[ "g/s B" ])
r_apply <- apply ( grid, 1 , function ( goodsOfPrsn1){
goods <- matrix ( c ( goodsOfPrsn1, sum_of_goods - goodsOfPrsn1),
2 , 2 , dimnames= dimnames ( initial_goods), byrow= TRUE )
equilibrium <- getEquilibrium ( goods)
with ( equilibrium, c ( goods[ "person 1" , "g/s A" ], goods[ "person 1" , "g/s B" ]))
})
r_apply <- cbind ( r_apply, with ( equilibrium, c ( goods[ "person 1" , "g/s A" ], goods[ "person 1" , "g/s B" ])))
rownames ( r_apply) <- c ( "g/s A" , "g/s B" )
r_apply <- r_apply[, order ( r_apply[ "g/s A" , ])]
lines ( r_apply[ "g/s A" , ], r_apply[ "g/s B" , ], lty= 3 , col= "darkgray" )
drawIDC ( initial_goods, names= c ( expression ( U[ 1 ] ^ I), expression ( U[ 2 ] ^ I)), IsCore= TRUE )
with ( equilibrium, {
drawIDC ( goods, col= col, names= c ( expression ( U[ 1 ] ^ E), expression ( U[ 2 ] ^ E)))
xadj <- ( par ()$ usr[ 2 ] - par ()$ usr[ 1 ]) / 50
pch <- 21
points ( initial_goods[ "person 1" , "g/s A" ], initial_goods[ "person 1" , "g/s B" ], pch= pch, col= "blue" , bg= "blue" )
text ( initial_goods[ "person 1" , "g/s A" ] + xadj, initial_goods[ "person 1" , "g/s B" ], expression ( I), adj= c ( 0 , 0 ))
points ( goods[ "person 1" , "g/s A" ], goods[ "person 1" , "g/s B" ], pch= pch, col= "black" , bg= "black" )
text ( goods[ "person 1" , "g/s A" ] + xadj, goods[ "person 1" , "g/s B" ], expression ( E), adj= c ( 0 , 0 ))
A <- seq ( xlim[ 1 ], xlim[ 2 ], length.out= 2 )
lines ( A, shp ( A))
})
img_path <- "img"
if ( ! dir.exists ( img_path)){
dir.create ( img_path)
}
dev.copy ( png, file= sprintf ( "%s/edgeworth_box.png" , img_path), width= 600 , height= 400 , type= "cairo" ); dev.off ()
