R code for http://grantmcdermott.github.io/2010/10/25/why-we-need-maths-in-economics-2b_25/

maths-economics-2b-schelling

rm(list=ls())

library(dplyr)
library(ggplot2)
library(cowplot)
library(extrafont)

greens <- 
  data_frame(pop = seq(0, 50, length = 51)) %>%
  mutate(tol = 2 - pop/(50/2)) %>%
  mutate(race = "Greens")
reds <- 
  data_frame(pop = seq(0, 100, length = 101)) %>%
  mutate(tol = 2 - pop/(100/2)) %>%
  mutate(race = "Reds")

schelling <- bind_rows(greens, reds)

schelling %>%
  ggplot(aes(x = pop, y = tol, col = race)) +
  geom_line() +
  scale_colour_manual(values = c("Greens" = "Green", "Reds" = "Red"),
                      labels = expression(y[G]==2-x/25, 
                                          y[R]==2-x/50),
                      name = ""
                      ) +
  labs(x = "Population member (x)", y = "Ratio (y:1)",
       title = "Tolerance schedules") +
  scale_x_continuous(breaks = seq(0, 100, 10)) +
  scale_y_continuous(breaks = seq(0, 2, .2)) +
  theme(text = element_text(family = "Lato"),
        panel.grid.major = element_line(color = "gray80", size = 0.2)
        ) +
  ggsave("schelling-schedule.png", width = 8, height = 6)

# library(ggvis)
# schelling %>% 
#   mutate(cols = ifelse(race == "Greens", "lawngreen", "red")) %>%
#   group_by(race) %>%
#   ggvis(~pop, ~tol, stroke:= ~cols, strokeWidth:= 2) %>% 
#   layer_lines()


schelling <- 
  schelling %>%
  mutate(lvl = pop * tol)

motion <-
  bind_rows(
    data_frame(Reds = filter(schelling, race == "Reds")$pop, 
               Greens = filter(schelling, race == "Reds")$lvl) %>%
      mutate(race = "Reds"),
    data_frame(Reds = filter(schelling, race == "Greens")$lvl,
               Greens = filter(schelling, race == "Greens")$pop) %>%
      mutate(race = "Greens")
  )


y_r <- function(x) x*(2 - x/50) ## Reds
x_g <- function(y) y*(2 - y/25) ## Greens
y_g <- function(x) {
  ## To map Greens on to the Red axis
  ## Solve for the quadratic equation in x_g above to get roots:
  ## x = y*(2 - y/25) ==> y^2 - 50y + 25x = 0 y = [50 +/- sqrt((-50)^2 - 100*x)]/2
  x_1 = (50 + sqrt((-50)^2 - 100*x))/2
  x_2 = (50 - sqrt((-50)^2 - 100*x))/2
  return(c(x_1,x_2))
  } ## Greens

## Equilibria (intersection points)
intrsct <- 
  data_frame(
    Reds = c(0, 100, optimize(function(x) abs(y_r(x) - max(y_g(x))), c(0, 25))$minimum)
    ) %>%
  mutate(Greens = ifelse(Reds == 0, y_g(Reds), y_r(Reds))) %>%
  mutate(eqms = ifelse(Reds %in% c(0, 100), "Stable", "Unstable"))

## Arrows
phase <- data_frame(Reds = c(12, 10, 30, 40), 
                    Greens = c(15, 35, 50, 25), 
                    vReds = c(5, -10, -5, 10), 
                    vGreens = c(6.25, 5, -6.25, -5)) 

motion %>%
  ggplot(aes(x = Reds, y = Greens)) +
  geom_path(aes(col = race)) +
  geom_point(data = intrsct, aes(shape = eqms), size = 4, stroke = 1) +
  geom_segment(data = phase, aes(xend = Reds+vReds, yend = Greens+vGreens),
               arrow = arrow()) +
  labs(title = "Dynamics of motion") +
  scale_colour_manual(values = c("Greens" = "Green", "Reds" = "Red"),
                      guide = F) +
  scale_shape_manual(name = "Equilbrium",
                     values = c("Stable" = 21, "Unstable" = 24)) +
  scale_x_continuous(breaks = seq(0, 100, 10)) +
  scale_y_continuous(breaks = seq(0, 50, 5)) +
  theme(text = element_text(family = "Lato"),
        panel.grid.major = element_line(color = "gray80", size = 0.2)
        ) +
  ggsave("schelling-dynamics.png", width = 8, height = 6)

# motion %>% 
#   mutate(cols = ifelse(race == "Greens", "lawngreen", "red")) %>%
#   group_by(race) %>%
#   ggvis(~Reds, ~Greens, stroke:= ~cols, strokeWidth:= 2) %>% 
#   layer_paths()

##################
## Second example
# 100 members of a population with linear, downward-sloping preferences and a 
# median tolerance ratio of 2.5 implies a curve of y = 5 - 50/20

greens2 <- 
  data_frame(pop = seq(0, 100, length = 101)) %>%
  mutate(tol = 5 - pop/20) %>%
  mutate(race = "Greens")
reds2 <- 
  data_frame(pop = seq(0, 100, length = 101)) %>%
  mutate(tol = 5 - pop/20) %>%
  mutate(race = "Reds")

schelling2 <- bind_rows(greens2, reds2)

schelling2 %>%
  ggplot(aes(x = pop, y = tol, col = race)) +
  geom_line() +
  scale_colour_manual(values = c("Greens" = "Green", "Reds" = "Red"),
                      labels = expression(y[G]==5-x/20, 
                                          y[R]==5-x/20),
                      name = ""
                      ) +
  labs(x = "Population member (x)", y = "Ratio (y:1)",
       title = "Tolerance schedules (Example 2") +
  scale_x_continuous(breaks = seq(0, 100, 10)) +
  scale_y_continuous(breaks = seq(0, 5, .5)) +
  theme(text = element_text(family = "Lato"),
        panel.grid.major = element_line(color = "gray80", size = 0.2)
        ) +
  ggsave("schelling-schedule-2.png", width = 8, height = 6)

schelling2 <- 
  schelling2 %>%
  mutate(lvl = pop * tol)

motion2 <-
  bind_rows(
    data_frame(Reds = filter(schelling2, race == "Reds")$pop, 
               Greens = filter(schelling2, race == "Reds")$lvl) %>%
      mutate(race = "Reds"),
    data_frame(Reds = filter(schelling2, race == "Greens")$lvl,
               Greens = filter(schelling2, race == "Greens")$pop) %>%
      mutate(race = "Greens")
    )


y_r2 <- function(x) x*(5 - x/20) ## Reds
x_g2 <- function(y) y*(5 - y/20) ## Greens
y_g2 <- function(x) {
  ## To map Greens onto the Red axis
  ## Solve for the quadratic equation in x_g above to get roots:
  ## x = y*(5 - y/20) ==> y^2 - 100y + 20x = 0 ==> y = [100 +/- sqrt((-100)^2.5 - 80*x)]/2
  x_1 = (100 + sqrt((-100)^2 - 80*x))/2
  x_2 = (100 - sqrt((-100)^2 - 80*x))/2
  return(c(x_1,x_2))
} ## Greens

## Equilibria (intersection points)
intrsct2 <- 
  data_frame(
    Reds = c(0, 100, 
             optimize(function(x) abs(y_r2(x) - max(y_g2(x))), c(10, 30))$minimum,
             optimize(function(x) abs(y_r2(x) - max(y_g2(x))), c(30, 90))$minimum,
             optimize(function(x) abs(y_r2(x) - min(y_g2(x))), c(90, 120))$minimum
             )
    ) %>%
  mutate(Greens = ifelse(Reds == 0, y_g2(Reds), y_r2(Reds))) %>%
  mutate(eqms = ifelse(Reds %in% c(0, 100, median(Reds)), "Stable", "Unstable"))

## Arrows
phase2 <- data_frame(Reds = c(50, 110, 50, 110, 10, 20, 80, 120), 
                    Greens = c(50, 110, 110, 50, 70, 120, 15, 15), 
                    vReds = c(10, -10, 10, -10, -10, -10, 10, -10), 
                    vGreens = c(10, -10, -10, 10, 10, -10, -10, -10)) 

motion2 %>%
  ggplot(aes(x = Reds, y = Greens)) +
  geom_path(aes(col = race)) +
  geom_point(data = intrsct2, aes(shape = eqms), size = 4, stroke = 1) +
  geom_segment(data = phase2, aes(xend = Reds+vReds, yend = Greens+vGreens),
               arrow = arrow()) +
  labs(title = "Dynamics of motion (Example 2)") +
  scale_colour_manual(values = c("Greens" = "Green", "Reds" = "Red"),
                      guide = F) +
  scale_shape_manual(name = "Equilbrium",
                     values = c("Stable" = 21, "Unstable" = 24)) +
  scale_x_continuous(breaks = seq(0, 130, 10)) +
  scale_y_continuous(breaks = seq(0, 130, 10)) +
  theme(text = element_text(family = "Lato"),
        panel.grid.major = element_line(color = "gray80", size = 0.2)
        ) +
  ggsave("schelling-dynamics-2.png", width = 8, height = 6)