drawEdgeworthBox <- function(init_A=0.9, init_B=0.3, CC=TRUE, SHP=TRUE, I=TRUE, E=TRUE, CORE=TRUE){
    # 経済の財は合計1に固定し、かつそれぞれの配分の最小値を0.01とする
    modifyParam <- function(g){
        g <- max(0.01, g)
        g <- min(0.99, g)
        g
    }
    init_A <- modifyParam(init_A)
    init_B <- modifyParam(init_B)

    # 財の配分を表す行列
    initial_goods <- matrix(c(init_A, 1 - init_A, init_B, 1 - init_B), 2, 2,
        dimnames = list(c("person 1", "person 2"), c("g/s A", "g/s B")))

    # トランスログ型効用関数の雛形を定義
    U_pt <- function(b, gs_a, gs_b){
        # b[3]～b[5]≠0だと局所的に凸なだけなので、安易に財の総量を増やしたり、下手にパラメーターをいじると、計算がデタラメに。
        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))
    }

    # Person 1と2の効用関数を定義
    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){
        # 財配分を更新
        # 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)
            # Person 1とPerson 2の限界代替率が一致し、Person 1の支出と収入が一致する点が均衡
            # p*(x0 - x) + (y0 - y) = 0
            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
        }

        # 財Aの量が横軸になる，0は効用関数の定義域から外れるので調整している
        # 無差別曲線を2つ、片方ひっくり返して重ねている図のため、線が枠をはみ出る方がよいので、上限はあえてずらしてる
        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")

        # Person 2の無差別曲線は180度回転してプロットするため、経済全体の財の量からPerson 2の財の量を引いて、対応するPerson 1の財の量を求めてプロットしている
        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財の量を求める
            b <- max(IC_1[ , "g/s B"][IC_1[ , "g/s B"] < sum_of_goods["g/s B"]], na.rm = TRUE)
            # 最大のB財に対応する行を求める
            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財の量を求める
            b <- min(IC_2[ , "g/s B"][IC_2[ , "g/s B"] > 0], na.rm = TRUE)
            # 最小のB財に対応する行を求める
            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))
        }

        # 純粋交換経済のコア
        # Person 2の無差別曲線（の180度回転）のB財の量よりも、Person 1のB財の量が同じか少ない領域を求める
        if(IsCore){
            # 横軸のIC_1の財Aの量とIC_2の財Aの量の対応がずれているので、線形近似関数をつくる
            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))) # 行きはPerson 1の無差別曲線に沿って、帰りはPerson 2の〃
            # 面を書く関数で色を塗る
            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))
    })

    # 契約曲線の計算と描画
    if(CC){
        # 初期時点の組み合わせを格子状に大量に作って計算する
        grid <- expand.grid(
            # 0は効用関数の定義域から外れるので調整している
            "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"])

        # apply関数で行列の列ごとに操作する
        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")

        # A財の大小で並び替えておく
        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")
    }

    # 初期点を通る無差別曲線を描画
    if(I) drawIDC(initial_goods, names=c(expression(U[1]^I), expression(U[2]^I)), IsCore=CORE)

    with(equilibrium, {
        # 均衡を通る無差別曲線を描画
        if(E) drawIDC(goods, col=col, names=c(expression(U[1]^E), expression(U[2]^E)))

        xadj <- (par()$usr[2] - par()$usr[1])/50 # 文字表示の位置調整
        pch <- 21 # 点の形状（21:丸, 22:四角, 23:菱形, 24:三角，25:逆三角, etc...）

        # 初期配分を描画
        if(I){
            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))
        }

        # 競争均衡点を描画
        if(E){
            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))
        }

        # 分離超平面を描画
        if(SHP){
            A <- seq(xlim[1], xlim[2], length.out=2)
            lines(A, shp(A))
        }
    })
}