New Economic Geography model with R

New economic geography model with R

We first need to read in some packages

################################################################
# Read in libraries
################################################################

library(nleqslv)  # for solving system of nonlinear equations
library(ggplot2)  # for structurally making plots
library(ggthemes) # for using economist theme
library(dplyr)    # for data wrangling
library(cowplot)  # for combining plots

Then we define constants

################################################################
# Define parameters
################################################################

L       <- 2.0  # Total labor force
phi1    <- 0.48 # fraction of food works living in region 1
gam     <- 0.3  # fraction that works in manufacturing
eps     <- 5.0  # elasticity of demand
rho     <- 0.8  # substitution parameter of variety
bet     <- 0.8  # variable costs
alp     <- 0.08 # fixed costs
delta   <- 0.4  # budget share manufacturing

################################################################
# Define iterations and stepsize for transporation and lambda
################################################################

iter_l <- 999
step_l <- 0.001
start_l <- 0.001

iter_t <- 51
start_t <- 1.5
step_t <- 0.01

Our key optimalisation procedure looks as follows:

################################################################
# Definite optimal function
################################################################

equilibrium <- function(x){
  
    Y1 <- x[1]
    Y2 <- x[2]
    W1 <- x[3]
    W2 <- x[4]
    I1 <- x[5]
    I2 <- x[6]
    
    y <- rep(NA, length(x))
      
    y[1] <- Y1-phi1*(1-gam)*L-lam*gam*L*W1
    y[2] <- Y2-(1-phi1)*(1-gam)*L-(1-lam)*gam*L*W2
    y[3] <- W1-rho*bet^(-rho)*(delta/(alp*(eps-1)))^(1/eps)*(Y1*I1^(eps-1)+T^(1-eps)*Y2*I2^(eps-1))^(1/eps)
    y[4] <- W2-rho*bet^(-rho)*(delta/(alp*(eps-1)))^(1/eps)*(T^(1-eps)*Y1*I1^(eps-1)+Y2*I2^(eps-1))^(1/eps)
    y[5] <- I1-(gam*L/(alp*eps))^(1/(1-eps))*(bet/rho)*(lam*W1^(1-eps)+(1-lam)*T^(1-eps)*W2^(1-eps))^(1/(1-eps))
    y[6] <- I2-(gam*L/(alp*eps))^(1/(1-eps))*(bet/rho)*(lam*T^(1-eps)*W1^(1-eps)+(1-lam)*W2^(1-eps))^(1/(1-eps))
    
    return(y)
}

And finally we need the loop below to create the figures

################################################################
# Create the vector where the output is stored
# This is faster than using append
# we will only append the equilibrium dataframe to find the
# stable and unstable equiliria (do that in the slower (outer) 
# loop)
################################################################

rel       <- vector(length = iter_l*iter_t)
lambda    <- vector(length = iter_l*iter_t)
transport <- vector(length = iter_l*iter_t)
welfare   <- vector(length = iter_l*iter_t)
w_man_h   <- vector(length = iter_l*iter_t)
w_man_f   <- vector(length = iter_l*iter_t)
w_farm_h  <- vector(length = iter_l*iter_t)
w_farm_f  <- vector(length = iter_l*iter_t)

################################################################
# Set the double loop for the  optimal solution using the  
# package nleqslv. 
# The fast (inner) loop is over gamma, The slow (outer) loop is 
# over the transportation costs
################################################################

# Completely parameterized
loop_transport <- seq( start_t, start_t + iter_t * step_t - step_t, by = step_t)
loop_gamma <- seq( start_l, start_l + iter_l * step_l - step_l, by = step_l )
equilibria <- data.frame(T = numeric(0), gamma = numeric(0), stable = numeric(0))

# Create intial starting values
start <- c(1,1,1,1,1,1)

iteration <- 0 # General counter
for (T in loop_transport){
  iter_eq <- 0 # Counter to find the equilibria for lambda
  lam_vec <- vector(length = iter_l) # initialize lambda vector
  t_vec   <- vector(length = iter_l) # initialize transport vector
  rel_vec <- vector(length = iter_l) # initialize relative real wage diff vector
  for (lam in loop_gamma){
    iteration <- iteration + 1
    iter_eq   <-  iter_eq + 1
    opt <- nleqslv(start, equilibrium)
    Y1 <- opt$x[1]
    Y2 <- opt$x[2]
    W1 <- opt$x[3]
    W2 <- opt$x[4]
    I1 <- opt$x[5]
    I2 <- opt$x[6]
    
    # Fill the various vectors
    rel[iteration]       <- (W1/I1^delta)/(W2/I2^delta)
    welfare[iteration]   <- Y1/(I1^delta)+Y2/(I2^delta)
    w_man_h[iteration]   <- W1/I1^delta 
    w_man_f[iteration]   <- W2/I2^delta
    w_farm_h[iteration]  <- 1/I1^delta 
    w_farm_f[iteration]  <- 1/I2^delta 
    lambda[iteration]    <- lam
    transport[iteration] <- T
    
    # Needed to find the equilibria (a bit redundant but more readible so)
    lam_vec[iter_eq]     <- lam
    t_vec[iter_eq]       <- T
    rel_vec[iter_eq]     <- (W1/I1^delta)/(W2/I2^delta)
  }
  eq <- data.frame(t_vec, lam_vec, rel_vec)
  eq <- eq %>% 
    mutate(
           dpos = ifelse( ( (rel_vec - 1) >= 0 & ( lag(rel_vec) - 1)  < 0), 1, 0 ),
           dneg = ifelse( ( (rel_vec - 1) <= 0 & ( lag(rel_vec) - 1)  > 0), 1, 0 )
           )
  stable <- eq %>% 
    filter(dneg == 1) %>%
    mutate(stable =1) %>%
    select(-dpos)
  unstable <- eq %>% 
    filter(dpos == 1) %>%
    mutate(stable =0) %>%
    select(-dneg)
  if (nrow(stable) > 0 ) {
      equilibria <- rbind(equilibria, data.frame(stable[1], stable[2], stable[5]))
      } 
  if (nrow(unstable) > 0) {
      equilibria <- rbind(equilibria, data.frame(unstable[1], unstable[2], unstable[5]))
      }
  if (nrow(unstable) == 1){
    equilibria <- rbind(equilibria, c(unstable[1,1], 0, 1))
    equilibria <- rbind(equilibria, c(unstable[1,1], 1, 1))
  }
  if ( (nrow(unstable) == 1) & (nrow(stable) == 1) ) {
    if (stable$lam_vec[1] > unstable$lam_vec[1] ){
      equilibria <- rbind(equilibria, c(unstable[1,1], 0, 1))
    }  
    if (stable$lam_vec[1] < unstable$lam_vec[1] ){
      equilibria <- rbind(equilibria, c(unstable[1,1], 1, 1))
    }
  }
  if ((nrow(unstable) + nrow(stable)) == 3) {
      equilibria <- rbind(equilibria, c(unstable[1,1], 1, 1))
      equilibria <- rbind(equilibria, c(unstable[1,1], 0, 1))
  }
}

################################################################
# Create the dataframe called neg_data
################################################################

neg_data <- data.frame(transport, lambda, rel, welfare,
                       w_man_h, w_man_f, w_farm_h, w_farm_f)

################################################################
# Create the plots
################################################################

#Indicate which lines should be highlighted
top_line <- neg_data[neg_data$transport == "1.5", ]
bottom_line <- neg_data[neg_data$transport == "2", ]
mid_line <- neg_data[neg_data$transport == "1.75", ]

p_wiggle <- ggplot(neg_data) + aes(lambda, rel, group = transport) + geom_line(size = 0.5, colour="grey", alpha = 0.5) +
  geom_line(data = top_line, aes(x = lambda, y = rel, group = transport, colour = "steelblue"), size = 1) +
  geom_line(data = bottom_line, aes(x = lambda, y = rel, group = transport, colour = "black"), size = 1) +
  geom_line(data = mid_line, aes(x = lambda, y = rel, group = transport, colour = "red"), size = 1) +
  scale_colour_discrete(name = "Transportation costs", labels = c("High", "Medium", "Low")) +
  geom_hline(yintercept = 1, size = 1, colour = "red", linetype = 4) +
  theme_economist() +  
  labs(title ="Wiggle diagram", y = "Relative real wage",
       subtitle = "Changes in relative real wage with varying lambda and transportation costs")

p_w_man_h <- ggplot(neg_data) + aes(lambda, w_man_h, group = transport) + 
  geom_line(size = 0.5, colour="grey", alpha = 0.5) +
  geom_line(data = top_line, aes(x = lambda, y = w_man_h, group = transport, color = "steelblue"), size = 1) +
  geom_line(data = bottom_line, aes(x = lambda, y = w_man_h, group = transport, color = "black"), size = 1) +
  geom_line(data = mid_line, aes(x = lambda, y = w_man_h, group = transport, color = "red"), size = 1) +
  scale_colour_discrete(name = "Transportation costs", labels = c("High", "Medium", "Low")) +
  theme_economist() + 
  labs(title ="Wage manufacturing workers home", y = "Wage manufacturers home",
       subtitle = "Changes in wages of workers in manufacturing at home \n with varying lambda and transportation costs")

p_w_man_f <- ggplot(neg_data) + aes(lambda, w_man_f, group = transport) + 
  geom_line(size = 0.5, colour="grey", alpha = 0.5) +
  geom_line(data = top_line, aes(x = lambda, y = w_man_f, group = transport, color = "steelblue"), size = 1) +
  geom_line(data = bottom_line, aes(x = lambda, y = w_man_f, group = transport, color = "black"), size = 1) +
  geom_line(data = mid_line, aes(x = lambda, y = w_man_f, group = transport, color = "red"), size = 1) +
  scale_colour_discrete(name = "Transportation costs", labels = c("High", "Medium", "Low")) +
  theme_economist() + 
  labs(title ="Wage manufacturing workers foreign", y = "Wage manufacturers foreign",
       subtitle = "Changes in wages of workers in manufacturing foreign \n with varying lambda and transportation costs")

p_welfare <- ggplot(neg_data) + aes(lambda, welfare, group = transport) + 
  geom_line(size = 0.5, colour="grey", alpha = 0.5) +
  geom_line(data = top_line, aes(x = lambda, y = welfare, group = transport, color = "steelblue"), size = 1) +
  geom_line(data = bottom_line, aes(x = lambda, y = welfare, group = transport, color = "black"), size = 1) +
  geom_line(data = mid_line, aes(x = lambda, y = welfare, group = transport, color = "red"), size = 1) +
  scale_colour_discrete(name = "Transportation costs", labels = c("High", "Medium", "Low")) +
  theme_economist() +
  labs(title ="Welfare", y = "Total welfare",
       subtitle = "Changes in total welfare with varying lambda and transportation costs")

p_combined <- plot_grid(p_wiggle + theme(legend.position="bottom"), 
                        p_w_man_h + theme(legend.position="bottom"), 
                        p_welfare + theme(legend.position="bottom"), 
                        p_w_man_f + theme(legend.position="bottom"), 
                        labels= c("A", "B", "C", "D"), 
                        ncol=2
                        ) 

p_tomahawk <- ggplot(equilibria) + aes(t_vec, lam_vec) + 
  geom_point(aes(colour = factor(stable))) + 
  theme_economist() +
  theme(legend.title=element_blank()) +
  scale_colour_discrete(breaks = c("0", "1"), labels=c("Unstable equilibrium", "Stable equilibrium")) +
  labs(title ="Tomahawk", y = "lambda", x = "transportation costs")

So we can get our tomahawk picture we want:

ggplot(equilibria) + aes(t_vec, lam_vec) + 
  geom_point(aes(colour = factor(stable))) + 
  theme_economist() +
  theme(legend.title=element_blank()) +
  scale_colour_discrete(breaks = c("0", "1"), labels=c("Unstable equilibrium", "Stable equilibrium")) +
  labs(title ="Tomahawk", y = "lambda", x = "transportation costs")

comments powered by Disqus