Monocentric city model with two income classes

JOURNAL OF ECONOMIC THEORY 13, 396-413 (1976) Comparative Statics of a Residential Economy with Several Classes JOHN HARTWICK

    library(rootSolve)
## Warning: package 'rootSolve' was built under R version 4.0.3
    a <- 0.35
    t <- 200
    w <- c(5000,12000)
    N <- c(500,500)
    J <- 2
    x <- rep(1,J+1)
    
    u <- rep(151.332,J)


    # Function to make bid-max lot size for agent of type j
    # with income and utility (w_j,u_j). 
    # Bid-max lot size is a function of distance to center d
    # given type specifics (w_j,u_j).

    make_g <- function(u_j,w_j,t,a)
        {
            g <- function(d)
                {
                    nominator <- a*u_j^(1/a)
                    denominator <- (w_j-t*d)^((1-a)/a)*.6*pi*d
                    bid_lot <- nominator/denominator
                    # return should be 1/bid_lot
                    # but denominator can be 0 so do not calculate it
                    # and use it as denominator instead 
                    denominator <- ifelse((w_j-t*d)>0,denominator,0)
                    return(denominator/nominator)
                }
            return(g)
        }


    # Solve the integral of g with respect to r from x_L = x_j to x_U = x_j+1
    # such that the integral equals N_j 

    # Initiate with j=1 so x_l is set to 0 and u_j is some guess
    j <- 1
    x_L <- x[j]
    u_j <- u[j]
        
        # Get exogenous income and population size for relevant segment
        # of population.
        w_j <- w[j]
        N_j <- N[j]

        g <- make_g(u_j,w_j,t,a)
        v <- function(x_U) {return(N_j-integrate(g,x_L,x_U)$value)} 
        x[j+1] <- uniroot(v,c(x_L+0.0001,500))$root

    # Initiate j=2
    j <- 2
        
        # Find u_2 given I_1,I_2 at x and u_1
        u[j] <- u[j-1]*(w[j] - t*x[j])/(w[j-1]-t*x[j])
    
    x_L <- x[j]
    u_j <- u[j]

        # Get exogenous income and population size for relevant segment
        # of population.
        w_j <- w[j]
        N_j <- N[j]

        g <- make_g(u_j,w_j,t,a)
        v <- function(x_U) {return(N_j-integrate(g,x_L,x_U)$value)} 
        x[j+1] <- uniroot(v,c(x_L+0.0001,500))$root

    

    # Calculate price at x[j+1]
    j <- 3  
    p_edge <- ((w[j-1] - t*x[j])/u[j-1])^(1/a)
    x
## [1]  1.000000  8.266838 20.756959
    u
## [1] 151.3320 467.8663
    p_edge
## [1] 3155.438
    p_max <- (w[1]/u[1])^(1/a)+10
    plot(x[1]:(x[J+1]+1),type="n",ylim=c(0,p_max))
    for (j in 1:J)
        {
            x_L <- x[j]
            x_U <- x[j+1]
    
            d <- seq(x_L,x_U,length.out=100)
            u_j <- u[j] 
            p <- ((w[j] - t*d)/u_j)^(1/a)
            points(d,p,type="l")

        }