QSM_Cleveland_1.0
Archer Thomas
2025-06-28
setwd("C:/Users/arche/Desktop/Summer 25/Site Fund/qsm/qsm_data")current_data <- st_read("cle_metro_shape.shp")## Reading layer `cle_metro_shape' from data source
## `C:\Users\arche\Desktop\Summer 25\Site Fund\qsm\qsm_data\cle_metro_shape.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 628 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 612698 ymin: 146975.8 xmax: 725246.8 ymax: 244076.2
## Projected CRS: NAD83(2011) / Ohio Northnaive_dist_matrix <- st_distance(current_data)
naive_dist_matrix <- units::drop_units(naive_dist_matrix)
naive_dist_matrix <- naive_dist_matrix / 1000
kmh <- 30
naive_dist_matrix <- (naive_dist_matrix / kmh) * 60A Quantitative Spatial Model of the Cleveland-Elyria Metropolitan Region
Introduction
t
Methodology
Data
The ICGities R package requires four spatial datasets covering the entirety of the urban/commuting area in question: the residential population, the workplace population, the area of each spatial unit, and average floorspace prices (per square foot or meter). In general, the spatial unit in question should be as granular as possible while also remaining meaningful in terms of sample size. I therefore chose to employ
Results
current_data[is.na(current_data)] <- 1
current_data[] <- lapply(current_data, function(col) {
if (is.numeric(col)) {
col[col == 0] <- 1
}
return(col)
})
N=as.matrix(nrow(current_data))
L_i=as.matrix(current_data$home_pop)
L_j=as.matrix(current_data$work_pop)
K=as.matrix(current_data$X_area)
t_ij=naive_dist_matrix
Q=as.matrix(current_data$X_idw_mean)inverted_values <- inversionModel(N=N,L_i=L_i,L_j=L_j,Q=Q,K=K,t_ij=t_ij)## Inverting wages...
## Iteration: 10, error: 2.9762945778.
## Iteration: 20, error: 1.418065065.
## Iteration: 30, error: 0.9008252114.
## Iteration: 40, error: 0.8693331511.
## Iteration: 50, error: 0.8373611153.
## Iteration: 60, error: 0.8048169367.
## Iteration: 70, error: 0.7717285332.
## Iteration: 80, error: 0.7381618667.
## Iteration: 90, error: 0.7041680088.
## Iteration: 100, error: 0.6697452643.
## Iteration: 110, error: 0.6348123394.
## Iteration: 120, error: 0.5991923008.
## Iteration: 130, error: 0.5626097764.
## Iteration: 140, error: 0.5247072532.
## Iteration: 150, error: 0.4850903275.
## Iteration: 160, error: 0.4434139818.
## Iteration: 170, error: 0.3995169357.
## Iteration: 180, error: 0.3535914419.
## Iteration: 190, error: 0.3063390324.
## Iteration: 200, error: 0.2590246423.
## Iteration: 210, error: 0.2133431287.
## Iteration: 220, error: 0.1710896294.
## Iteration: 230, error: 0.1337491198.
## Iteration: 240, error: 0.1021881878.
## Iteration: 250, error: 0.0765712606.
## Iteration: 260, error: 0.0564875146.
## Iteration: 270, error: 0.0411802545.
## Iteration: 280, error: 0.0297683389.
## Iteration: 290, error: 0.0269680135.
## Iteration: 300, error: 0.02504048.
## Iteration: 310, error: 0.0232455921.
## Iteration: 320, error: 0.0215737603.
## Iteration: 330, error: 0.0200163951.
## Iteration: 340, error: 0.0185657223.
## Iteration: 350, error: 0.0172146411.
## Iteration: 360, error: 0.0159566175.
## Iteration: 370, error: 0.0147856009.
## Iteration: 380, error: 0.0136959612.
## Iteration: 390, error: 0.0126824394.
## Iteration: 400, error: 0.0117401095.
## Iteration: 410, error: 0.0108643475.
## Iteration: 420, error: 0.0100508077.
## Iteration: 430, error: 0.0092954021.
## Iteration: 440, error: 0.008594284.
## Iteration: 450, error: 0.0079438332.
## Iteration: 460, error: 0.0073406442.
## Iteration: 470, error: 0.006781514.
## Iteration: 480, error: 0.0062634322.
## Iteration: 490, error: 0.0057835715.
## Iteration: 500, error: 0.0053392782.
## Iteration: 510, error: 0.0049280637.
## Iteration: 520, error: 0.0045475963.
## Iteration: 530, error: 0.0041956928.
## Iteration: 540, error: 0.0038703111.
## Iteration: 550, error: 0.0035695427.
## Iteration: 560, error: 0.0032916049.
## Iteration: 570, error: 0.0030348346.
## Iteration: 580, error: 0.0027976809.
## Iteration: 590, error: 0.0025786986.
## Iteration: 600, error: 0.0023765425.
## Iteration: 610, error: 0.0021899605.
## Iteration: 620, error: 0.0020177889.
## Iteration: 630, error: 0.0018589458.
## Iteration: 640, error: 0.0017124269.
## Iteration: 650, error: 0.0015772996.
## Iteration: 660, error: 0.0014526991.
## Iteration: 670, error: 0.0013378233.
## Iteration: 680, error: 0.001231929.
## Iteration: 690, error: 0.0011343276.
## Iteration: 700, error: 0.0010443816.
## Iteration: 710, error: 0.0009615007.
## Iteration: 720, error: 0.0008851391.
## Iteration: 730, error: 0.0008147917.
## Iteration: 740, error: 0.0007499916.
## Iteration: 750, error: 0.0006903071.
## Iteration: 760, error: 0.0006353396.
## Iteration: 770, error: 0.0005847207.
## Iteration: 780, error: 0.0005381101.
## Iteration: 790, error: 0.0004951938.
## Iteration: 800, error: 0.0004556817.
## Iteration: 810, error: 0.0004193065.
## Iteration: 820, error: 0.0003858212.
## Iteration: 830, error: 0.000354998.
## Iteration: 840, error: 0.0003266269.
## Iteration: 850, error: 0.0003005144.
## Iteration: 860, error: 0.0002764816.
## Iteration: 870, error: 0.0002543642.
## Iteration: 880, error: 0.0002340103.
## Iteration: 890, error: 0.0002152801.
## Iteration: 900, error: 0.0001980447.
## Iteration: 910, error: 0.0001821855.
## Iteration: 920, error: 0.0001675931.
## Iteration: 930, error: 0.0001541666.
## Iteration: 940, error: 0.0001418134.
## Iteration: 950, error: 0.000130448.
## Iteration: 960, error: 0.0001199916.
## Iteration: 970, error: 0.0001103718.
## Iteration: 980, error: 0.0001015219.
## Iteration: 990, error: 9.33805e-05.
## Iteration: 1000, error: 8.58909e-05.
## Reached maximum number of iterations (1000). Error=8.58909e-05.a=as.matrix(inverted_values$a)
b=as.matrix(inverted_values$b)
varphi=as.matrix(inverted_values$varphi)
w=as.matrix(inverted_values$w)
u=as.matrix(inverted_values$u)
Q_norm=as.matrix(inverted_values$Q_norm)
ttheta=as.matrix(inverted_values$ttheta)productivity_before <- st_sf(as.data.frame(inverted_values$A), geometry = st_geometry(current_data))
plot(productivity_before)
amenities_before <- st_sf(as.data.frame(inverted_values$B), geometry = st_geometry(current_data))
plot(amenities_before)
amenities_before <- st_sf(as.data.frame(b), geometry = st_geometry(current_data))
plot(amenities_before)
model_closed <- IGCities::solveModel(
N=N,
L_i=L_i,
L_j=L_j,
K=K,
t_ij=t_ij,
a=a,
b=b,
varphi=varphi,
w_eq=w,
u_eq=u,
Q_eq=Q_norm,
ttheta_eq=ttheta
)## Solving model...
## Iteration: 10, error: 2.09073e-05.
## Iteration: 20, error: 7.047e-07.
## Iteration: 30, error: 3.8548e-06.
## Iteration: 40, error: 1.7503e-06.
## Iteration: 50, error: 1.879e-07.
## Iteration: 60, error: 2.105e-07.
## Iteration: 70, error: 1.367e-07.
## Iteration: 80, error: 3.21e-08.
## Iteration: 90, error: 7.4e-09.
## Iteration: 100, error: 8.3e-09.
## Iteration: 110, error: 2.7e-09.
## Iteration: 120, error: 5e-10.
## Iteration: 130, error: 8e-10.
## Iteration: 140, error: 4e-10.
## Iteration: 150, error: 2e-10.
## Converged after 155 iterations. Error=1e-10.#' Array operator to mimic different-dimension-array element-wise operations in
#' MATLAB. It receives as input two arrays of potentially different dimensions,
#' it resizes them to have same dimensions and finally performs the element-wise
#' operation.
#'
#' @param array1 The first array
#' @param array2 The second array
#' @param operation The operation. It can take values: '+', '-', '*',
#' '/' and '^'
#'
#' @return An array with dimensions equal to the "largest" input array. It is
#' the result of applying the operator element-wise to both input arrays.
array_operator = function(array1, array2, operation){
if(is.array(array1) == FALSE){
array1 = array(array1, dim=dim(array2))
} else if(is.array(array2) == FALSE){
array2 = array(array2, dim=dim(array1))
}
dim1 = dim(array1)
dim2 = dim(array2)
if(length(dim1) < length(dim2)){
array1 = kronecker(array1, array(1, dim = array(1, dim=length(dim2))))
dim1 = dim(array1)
} else if(length(dim1) > length(dim2)){
array2 = kronecker(array2, array(1, dim = array(1, dim=length(dim1))))
dim2 = dim(array2)
}
if(min(dim1-dim2) < 0){
array1_reshaped = kronecker(array1, array(1, dim=dim2-dim1+1))
array2_reshaped = array2
} else if(max(dim1-dim2) > 0){
array1_reshaped = array1
array2_reshaped = kronecker(array2, array(1, dim=dim1-dim2+1))
} else{
array1_reshaped = array1
array2_reshaped = array2
}
if(operation == '*'){
array_output = array1_reshaped*array2_reshaped
} else if(operation == '/'){
array_output = array1_reshaped/array2_reshaped
} else if(operation == '+'){
array_output = array1_reshaped+array2_reshaped
} else if(operation == '-'){
array_output = array1_reshaped-array2_reshaped
} else if(operation == '^'){
array_output = array1_reshaped^array2_reshaped
}
return(array_output)
}
#' Computes average income in each location, which is the weighted average of
#' the income of the people living in the location.
#'
#' @param lambda_ij_i NxN matrix - Probability of individuals in each
#' location of working in each location.
#' @param w_tr NxS - Wages in each location in each sector.
av_income_simple = function(lambda_ij_i,
w_tr){
y_bar = (sumDims2(array_operator(lambda_ij_i, w_tr, '*'), 2));
return(list(y_bar=y_bar))
}
#' Function to transform travel times into iceberg commuting costs
#' @param t_ij NxN matrix - Travel time matrix across locations
#' @param epsilon Float - Parameter that transforms travel times to commuting costs
#'
#' @return A NxN matrix of commuting costs
commuting_matrix = function(t_ij,
epsilon){
tau = exp(epsilon*t_ij)
return(list(tau=tau))
}
#' Computes residential and commercial floorspace supply and equilibrium prices.
#'
#' @param Q Nx1 array - Floorspaces prices.
#' @param K Nx1 array - Land supply.
#' @param w NxS - Wages in each location in each sector.
#' @param L_j Nx1 matrix - Number of workers in each location.
#' @param y_bar - Average income in each location.
#' @param L_i Nx1 matrix - Number of residents in each location.
#' @param beta Float - Cobb-Douglas parameter output elasticity wrt labor.
#' @param alpha Float - Utility parameter that determines preferences for
#' consumption.
#' @param mu Float - Floorspace prod function: output elast wrt capita, 1-mu wrt land.
density_development = function(Q,
K,
w,
L_j,
y_bar,
L_i,
beta,
alpha,
mu){
Q_mean = exp(mean(log(Q)));
Q_norm = Q/Q_mean;
FS_f = ((1-beta)/beta)*(array_operator(array_operator(w, L_j, '*'), Q_norm, '/'));
FS_r = (1-alpha)*(array_operator(array_operator(y_bar, L_i, '*'),Q_norm,'/'));
FS = FS_f+FS_r;
varphi = array_operator(FS, K^(1-mu), '/');
return(list(Q_mean = Q_mean, Q_norm = Q_norm, FS_f = FS_f, FS_r = FS_r, FS = FS, varphi=varphi))
}
#' Function to invert model, so amenities, wages, productivities, and development density
#' are chosen to match model to data.
#'
#' @param N Integer - Number of locations.
#' @param L_i Nx1 matrix - Number of residents in each location.
#' @param L_j Nx1 matrix - Number of workers in each location.
#' @param Q Nx1 matrix - Floorspace prices
#' @param K Nx1 matrix - Land area
#' @param t_ij NxN matrix - Travel times across all possible locations.
#' @param alpha Float - Utility parameter that determines preferences for
#' consumption.
#' @param beta Float - Output elasticity wrt labor
#' @param theta Float - Commuting elasticity and migration elasticity.
#' @param delta Float - Decay parameter agglomeration
#' @param rho Float - Decay parameter congestion
#' @param lambda Float - Agglomeration force
#' @param epsilon Float - Parameter that transforms travel times to commuting costs
#' @param mu Float - Floorspace prod function: output elast wrt capital, 1-mu wrt land.
#' @param eta Float - Congestion force
#' @param nu_init Float - Convergence parameter to update wages.
#' Default nu=0.01.
#' @param tol Int - tolerance factor
#' @param maxiter Integer - Maximum number of iterations for convergence.
#' Default maxiter=1000.
#'
#' @return Equilibrium values.
#' @export
#'
#' @examples
#' N=5
#' L_i = c(63, 261, 213, 182, 113)
#' L_j = c(86, 278, 189, 180, 99)
#' Q = c(2123, 1576, 1371, 1931, 1637)
#' K = c(0.44, 1.45, 1.15, 0.87, 0.58)
#' t_ij = rbind(c(0.0, 6.6, 5.5, 5.6, 6.4),
#' c(6.7, 0.0, 3.9, 4.6, 4.4),
#' c(5.5, 3.9, 0.0, 2.8, 3.0),
#' c(5.6, 4.6, 2.8, 0.0, 2.7),
#' c(6.4, 4.4, 3.0, 2.7, 0.0))
#'
#' inversionModel(N=N,
#' L_i=L_i,
#' L_j=L_j,
#' Q=Q,
#' K=K,
#' t_ij=t_ij)
#'
inversionModel = function(N,
L_i,
L_j,
Q,
K,
t_ij,
alpha=0.7,
beta=0.7,
theta=7,
delta=0.3585,
rho=0.9094,
lambda=0.01,
epsilon=0.01,
mu=0.3,
eta=0.1548,
nu_init=0.005,
tol=10^-10,
maxiter=1000){
# Formatting of input data
if(is.data.frame(L_i)){
L_i = array(unlist(L_i), dim(L_i))
} else if(is.null(dim(L_i))){
L_i = array(L_i, dim=c(N,1))
}
if(is.data.frame(L_j)){
L_j = array(unlist(L_j), dim(L_j))
} else if(is.null(dim(L_j))){
L_j = array(L_j, dim=c(N,1))
}
if(is.data.frame(K)){
K = array(unlist(K), dim(K))
} else if(is.null(dim(K))){
K = array(K, dim=c(N,1))
}
if(is.data.frame(Q)){
Q = array(unlist(Q), dim(Q))
} else if(is.null(dim(Q))){
Q = array(Q, dim=c(N,1))
}
t_ij = array(unlist(t_ij), dim(t_ij))
# Normalize L_i to have the same size as L_j
L_i=L_i*sum(L_j)/sum(L_i)
# Initialization
w_init=array(1, dim=c(N,1))
# Transformation of travel times to trade costs
D = commuting_matrix(t_ij=t_ij,
epsilon=epsilon)
tau = D$tau
# Finding the wages that match the data
WI = wages_inversion(N=N,
w_init=w_init,
theta=theta,
tau=tau,
L_i=L_i,
L_j=L_j,
nu_init=nu_init,
tol=tol,
maxiter=maxiter)
# Equilibrium wages
w = WI$w
w_tr = WI$w_tr
W_i = WI$W_i
lambda_ij_i = WI$lambda_ij_i
# Average income
Inc = av_income_simple(lambda_ij_i=lambda_ij_i,
w_tr = w_tr
)
y_bar = Inc$y_bar
#Density of development
DensD = density_development(Q=Q,
K=K,
w=w,
L_j=L_j,
y_bar=y_bar,
L_i=L_i,
beta=beta,
alpha=alpha,
mu=mu
)
Q_mean = DensD$Q_mean
Q_norm = DensD$Q_norm
FS_f = DensD$FS_f
FS_r = DensD$FS_r
FS = DensD$FS
varphi = DensD$varphi
ttheta = FS_f/FS
#Productivities
Prod = productivity(N=N,
Q=Q,
w=w,
L_j=L_j,
K=K,
t_ij = t_ij,
delta=delta,
lambda=lambda,
beta=beta
)
A = Prod$A
a = Prod$a
# Amenities
AM = living_amenities_simple(theta=theta,
N=N,
L_i=L_i,
W_i=W_i,
Q=Q,
K=K,
alpha=alpha,
t_ij=t_ij,
rho=rho,
eta=eta
)
B = AM$B
b = AM$b
# Save and export
Q_alpha = Q_norm^(1-alpha)
u = array_operator(array_operator(W_i,Q_alpha,'/'),B,'*')
U = (sumDims(u,1))^(1/theta)
return(list(A=A, a=a, u=u, B=B, b=b, w=w, varphi=varphi, U=U, Q_norm=Q_norm, ttheta=ttheta))
}
#' Function to estimate amenity parameters of locations where users live.
#'
#' @param theta Float - Parameter that governs the reallocation of workers across
#' locations in the city. This parameter measures how sensible are migration
#' flows within the city to changes in real income.
#' @param N Integer - Number of locations.
#' @param L_i Nx1 matrix - Total residents.
#' @param W_i Nx1 matrix - Market access measure in each location.
#' @param Q Nx1 matrix - Floor space prices.
#' @param K Nx1 matrix - Land area
#' @param alpha Float - Para
#' @param t_ij NxN matrix - Travel times across locations.
#' @param rho Float - decay parameter for amenities.
#' @param eta Float - congestion force
#'
#' @return Matrix with the amenity distribution of living in each location.
living_amenities_simple = function(theta,
N,
L_i,
W_i,
Q,
K,
alpha,
t_ij,
rho,
eta){
Q_mean = exp(mean(log(Q)));
Q_norm = Q/Q_mean;
L_i_mean = exp(mean(log(L_i)));
L_i_norm = L_i/L_i_mean;
W_i_mean = exp(mean(log(W_i)));
W_i_norm = W_i/W_i_mean;
B = array_operator(array_operator(L_i_norm^(1/theta), Q_norm^(1-alpha), '*'), W_i_norm^((-1)), '*');
L_i_dens = (array_operator(L_i, K, '/'));
L_i_dens_per = aperm(array(L_i_dens, dim=c(N,1)), c(2,1));
L_i_dens_rep = kronecker(L_i_dens_per, array(1, dim=c(N, 1)));
Omega = sumDims2(array_operator(exp(-rho*t_ij), L_i_dens_rep, '*'), 2);
b = array_operator(B, Omega^(-eta), "/");
b = b*(L_i>0);
return(list(B = B, b = b))
}
#' Computes productivity levels in each location
#'
#' @param N Float - Number of locations.
#' @param Q Nx1 matrix - Floorspace prices in each location.
#' @param w Nx1 matrix - wages in each location.
#' @param L_j Nx1 matrix - Employment in each location.
#' @param K Nx1 matrix - Land in each location.
#' @param t_ij NxN matrix - Travel times matrix.
#' @param delta Float - decay parameter agglomeration.
#' @param lambda Float - agglomeration force.
#' @param beta Float - Output elasticity wrt labor
productivity = function(N,
Q,
w,
L_j,
K,
t_ij,
delta,
lambda,
beta){
Q_mean = exp(mean(log(Q)));
Q_norm = Q/Q_mean;
beta_tilde = ((1-beta)^(1-beta))*(beta^beta);
A = (1/beta_tilde)*(array_operator(Q_norm^(1-beta), w^beta, '*'));
L_j_dens = (array_operator(L_j, K, '/'));
L_j_dens_per = aperm(array(L_j_dens, dim=c(N,1)), c(2,1));
L_j_dens_rep = kronecker(L_j_dens_per, array(1, dim=c(N, 1)));
Upsilon = sumDims2(array_operator(exp(-delta*t_ij), L_j_dens_rep, '*'), 2);
a = array_operator(A, Upsilon^lambda, "/");
a = a*(L_j>0)
return(list(A = A, a = a))
}
#' Collapse array along one of the dimensions by adding the elements along that
#' dimension.
#'
#' @param array Array to collapse along one dimension.
#' @param dimension Dimension to collapse the array.
#'
#' @return An array that has been collapsed along the given dimension.
sumDims = function(array, dimension){
dim1 = length(dim(array))
keep_dims = setdiff(seq(1, dim1, length.out = dim1), c(dimension))
array_output = apply(array, MARGIN = keep_dims, FUN = 'sum')
array_output = kronecker(array_output, array(1, dim=array(1, dim=dim1)))
dim_reshape = seq(1, dim1, length.out = dim1)
dim_reshape[(dimension+1):length(dim_reshape)] = dim_reshape[dimension:(length(dim_reshape)-1)]
dim_reshape[dimension] = dim1
array_output = aperm(array_output, dim_reshape);
return(array_output)
}
#' Collapse array 2 along one of the dimensions by adding the elements along that
#' dimension.
#'
#' @param array Array to collapse along one dimension.
#' @param dimension Dimension to collapse the array.
#'
#' @return An array that has been collapsed along the given dimension.
sumDims2 = function(array, dimension){
dim1 = length(dim(array))
keep_dims = setdiff(seq(1, dim1, length.out = dim1), c(dimension))
array_output = apply(array, MARGIN = keep_dims, FUN = 'sum')
array_output = kronecker(array_output, array(1, dim=array(1, dim=dim1)))
return(array_output)
}
#' Function to compute equilibrium wages that make the model labor in every
#' location in equal to the observed data. It finds the w's
#' such that equation (3.2) holds.
#'
#' @param N Integer - Number of locations.
#' @param w_init Initial vector of wages.
#' @param theta Float - Commuting elasticity.
#' @param tau NxN matrix - Commuting cost matrix across all locations.
#' @param L_i Nx1 matrix - Number of residents in each location.
#' @param L_j Nx1 matrix - Number of workers in each location.
#' @param nu_init Float - Convergence parameter to update wages.
#' Default nu=0.01.
#' @param tol Float - Maximum tolerable error for estimating total labor.
#' Default tol=10^-10.
#' @param maxiter Integer - Maximum number of iterations for convergence.
#' Default maxiter=10000.
#'
#' @return A list with equilibrium wages and probability of workers in each
#' location working in every other location.
wages_inversion = function(N,
w_init,
theta,
tau,
L_i,
L_j,
nu_init=0.05,
tol=10^-10,
maxiter=10000){
# Settings
outerdiff = Inf
w = w_init
iter = 0
nu = nu_init
cat("Inverting wages...\n")
while(outerdiff>tol & iter<maxiter){
# 1) Labor supply
# Indirect utility
w_tr = aperm(array(w, dim=c(N,1)), c(2,1));
rep_w_tr = kronecker(w_tr^theta, array(1, dim=c(N, 1)));
# Constructing emp` loyment shares
w_tr_tau = array_operator(w_tr^theta, tau^(-theta), '*');
lambda_ij_i = array_operator(w_tr_tau, sumDims2(w_tr_tau,2), '/');
W_i = (sumDims2(w_tr_tau,2))^(1/theta);
# Labor is equal to probabilities * total number of residents * proportion of workers in each sector.
L_ij = array_operator(L_i, lambda_ij_i, '*')
L_j_tr = sumDims2(L_ij, 1)
# L_j_model = aperm(L_j_tr, c(2, 1));
Ratio_supply = array_operator(L_j_tr, w, "/");
w_prime = array_operator(L_j, Ratio_supply, "/");
z_L = array_operator(w, w_prime, '-');
w = array_operator(w*(1-nu), w_prime*nu, '+');
w_mean = exp(mean(log(w)))
w = w/w_mean;
outerdiff = max(abs(z_L))
iter = iter+1;
if(iter %% 10 == 0){
cat(paste0("Iteration: ", iter, ", error: ", round(outerdiff, 10), ".\n"))
}
}
if(outerdiff<=tol){
cat(paste0("Converged after ", iter, " iterations. Error=", round(outerdiff, 10), ".\n"))
} else{
cat(paste0("Reached maximum number of iterations (", iter, "). Error=", round(outerdiff, 10), ".\n"))
}
return(list(w=w, w_tr=w_tr, W_i=W_i, lambda_ij_i=lambda_ij_i))
}solveModel_open = function(N,
L_i,
L_j,
K,
t_ij,
a,
b,
varphi,
w_eq,
u_eq,
Q_eq,
ttheta_eq,
alpha=0.7,
beta=0.7,
theta=7,
mu=0.3,
delta=0.3585,
lambda=0.01,
rho=0.9094,
eta=0.1548,
epsilon=0.01,
zeta=0.95,
tol=10^-10,
maxiter=1000,
U_bar){
# Formatting of input data
if(is.data.frame(L_i)){
L_i = array(unlist(L_i), dim(L_i))
} else if(is.null(dim(L_i))){
L_i = array(L_i, dim=c(N,1))
}
if(is.data.frame(L_j)){
L_j = array(unlist(L_j), dim(L_j))
} else if(is.null(dim(L_j))){
L_j = array(L_j, dim=c(N,1))
}
if(is.data.frame(K)){
K = array(unlist(K), dim(K))
} else if(is.null(dim(K))){
K = array(K, dim=c(N,1))
}
t_ij = array(unlist(t_ij), dim(t_ij))
if(is.null(dim(a))){
a = array(a, dim=c(N,1))
}
if(is.null(dim(b))){
b = array(b, dim=c(N,1))
}
if(is.null(dim(varphi))){
varphi = array(varphi, dim=c(N,1))
}
if(is.null(dim(w_eq))){
w_eq = array(w_eq, dim=c(N,1))
}
if(is.null(dim(u_eq))){
u_eq = array(u_eq, dim=c(N,1))
}
if(is.null(dim(Q_eq))){
Q_eq = array(Q_eq, dim=c(N,1))
}
if(is.null(dim(ttheta_eq))){
ttheta_eq = array(ttheta_eq, dim=c(N,1))
}
# Normalize L_i to have the same size as L_j
L_i=L_i*sum(L_j)/sum(L_i)
D = commuting_matrix(t_ij=t_ij, epsilon = epsilon)
tau = D$tau
L_i = array(unlist(L_i),dim(L_i))
L_j = array(unlist(L_j),dim(L_j))
K = array(unlist(K), dim(K))
# Settings
outerdiff = Inf;
w = w_eq;
u = u_eq;
Q = Q_eq;
ttheta = ttheta_eq
iter = 0;
zeta_init = zeta;
cat("Solving model...\n")
while(outerdiff>tol & iter < maxiter){
# 1) Labor supply equation
w_tr = aperm(array(w, dim=c(N,1)), c(2,1));
rep_w_tr = kronecker(w_tr^theta, array(1, dim=c(N, 1)));
# Constructing employment shares
w_tr_tau = array_operator(w_tr^theta, tau^(-theta), '*');
lambda_ij_i = array_operator(w_tr_tau, sumDims2(w_tr_tau,2), '/');
W_i = (sumDims2(w_tr_tau,2))^(1/theta);
# Labor is equal to probabilities * total number of residents * proportion of workers in each sector.
L_ij = array_operator(L_i, lambda_ij_i, '*')
L_j = sumDims2(L_ij, 1)
L = sum(L_i)
lambda_i = L_i/L
# 2 average income
av_income = av_income_simple(lambda_ij_i=lambda_ij_i,w_tr = w_tr)
ybar = av_income$y_bar
# 3 Total floorspace
FS = array_operator(varphi,K^(1-mu),"*")
# 4 Agglomeration externalities
L_j_dens = (array_operator(L_j, K, '/'));
L_j_dens_per = aperm(array(L_j_dens, dim=c(N,1)), c(2,1));
L_j_dens_rep = kronecker(L_j_dens_per, array(1, dim=c(N, 1)));
Upsilon = sumDims2(array_operator(exp(-delta*t_ij), L_j_dens_rep, '*'), 2);
A = array_operator(a, Upsilon^lambda, '*')
# 5 Amenities
L_i_dens = (array_operator(L_i, K, '/'));
L_i_dens_per = aperm(array(L_i_dens, dim=c(N,1)), c(2,1));
L_i_dens_rep = kronecker(L_i_dens_per, array(1, dim=c(N, 1)));
Omega = sumDims2(array_operator(exp(-rho*t_ij), L_i_dens_rep, '*'), 2);
B = array_operator(b, Omega^(-eta),'*')
# 6 Residents, probabilities, and welfare
u = array_operator(array_operator(W_i, Q^(1-alpha), '/'), B, '*')
U = sum(u^theta)
lambda_i_upd = (u^theta)/U
U = U^(1/theta)
if (U > U_bar) {
L = L+1
} else if (U < U_bar) {
L = L-1
}
# 7 Total output by location
FS_f = array_operator(ttheta,array_operator(varphi, K^(1-mu), '*'), '*')
Y = array_operator(A, array_operator(L_j^beta, FS_f^(1-beta), '*'), '*')
Q_upd1 = (1-beta)*array_operator(Y,FS_f, '/')
w_upd = beta*array_operator(Y, L_j, '/')
# 8 Housing prices
FS_r = array_operator((1-ttheta), array_operator(varphi, K^(1-mu), '*'), '*')
X = array_operator(ybar, L_i, '*')
Q_upd2 = (1-alpha)*array_operator(X, FS_r, '/')
Q_upd = Q_upd1*(a>0) + Q_upd2*(a==0 & b>0)
# 9 Share of commercial floorspace
LP = array_operator(Q_upd1, array_operator(varphi, K^(1-mu), '*'), '*')
ttheta_upd = (1-beta)*array_operator(Y, LP, '/')
ttheta_upd = (b==0)+ttheta_upd*(b>0)
# 10 Calculating the main differences
z_w = array_operator(w, w_upd, '-')
z_L = array_operator(lambda_i, lambda_i_upd, '-')
z_Q = array_operator(Q, Q_upd, '-')
z_theta = array_operator(ttheta, ttheta_upd, '-')
outerdiff = max(c(max(abs(z_w)), max(abs(z_L)), max(abs(z_Q)), max(abs(z_theta))))
iter = iter+1
# 11 New vector of variables
lambda_i = zeta*lambda_i + (1-zeta)*lambda_i_upd
Q = zeta*Q + (1-zeta)*Q_upd
w = zeta*w + (1-zeta)*w_upd
ttheta = zeta*ttheta + (1-zeta)*ttheta_upd
L_i = lambda_i*L
if(iter %% 10 == 0){
cat(paste0("Iteration: ", iter, ", error: ", round(outerdiff, 10), ".\n"))
}
}
if(outerdiff<=tol){
cat(paste0("Converged after ", iter, " iterations. Error=", round(outerdiff, 10), ".\n"))
} else{
cat(paste0("Reached maximum number of iterations (", iter, "). Error=", round(outerdiff, 10), ".\n"))
}
return(list(w=w, W_i=W_i, B=B, A=A, Q=Q, lambda_ij_i=lambda_ij_i, L_i=L_i, L_j=L_j,
ybar=ybar, lambda_i=lambda_i, ttheta=ttheta, u=u, U=U))
}## Solving model...
## Iteration: 10, error: 2.18451e-05.
## Iteration: 20, error: 1.6792e-06.
## Iteration: 30, error: 3.2782e-06.
## Iteration: 40, error: 2.6658e-06.
## Iteration: 50, error: 1.2209e-06.
## Iteration: 60, error: 1.2597e-06.
## Iteration: 70, error: 1.1489e-06.
## Iteration: 80, error: 1.1923e-06.
## Iteration: 90, error: 1.1819e-06.
## Iteration: 100, error: 1.1783e-06.
## Iteration: 110, error: 1.1808e-06.
## Iteration: 120, error: 1.1808e-06.
## Iteration: 130, error: 1.18e-06.
## Iteration: 140, error: 1.1806e-06.
## Iteration: 150, error: 1.1803e-06.
## Iteration: 160, error: 1.1805e-06.
## Iteration: 170, error: 1.1804e-06.
## Iteration: 180, error: 1.1804e-06.
## Iteration: 190, error: 1.1804e-06.
## Iteration: 200, error: 1.1804e-06.
## Iteration: 210, error: 1.1804e-06.
## Iteration: 220, error: 1.1804e-06.
## Iteration: 230, error: 1.1804e-06.
## Iteration: 240, error: 1.1804e-06.
## Iteration: 250, error: 1.1804e-06.
## Iteration: 260, error: 1.1804e-06.
## Iteration: 270, error: 1.1804e-06.
## Iteration: 280, error: 1.1804e-06.
## Iteration: 290, error: 1.1804e-06.
## Iteration: 300, error: 1.1804e-06.
## Iteration: 310, error: 1.1804e-06.
## Iteration: 320, error: 1.1804e-06.
## Iteration: 330, error: 1.1804e-06.
## Iteration: 340, error: 1.1804e-06.
## Iteration: 350, error: 1.1804e-06.
## Iteration: 360, error: 1.1804e-06.
## Iteration: 370, error: 1.1804e-06.
## Iteration: 380, error: 1.1804e-06.
## Iteration: 390, error: 1.1804e-06.
## Iteration: 400, error: 1.1804e-06.
## Iteration: 410, error: 1.1804e-06.
## Iteration: 420, error: 1.1804e-06.
## Iteration: 430, error: 1.1804e-06.
## Iteration: 440, error: 1.1804e-06.
## Iteration: 450, error: 1.1804e-06.
## Iteration: 460, error: 1.1804e-06.
## Iteration: 470, error: 1.1804e-06.
## Iteration: 480, error: 1.1804e-06.
## Iteration: 490, error: 1.1804e-06.
## Iteration: 500, error: 1.1804e-06.
## Iteration: 510, error: 1.1804e-06.
## Iteration: 520, error: 1.1804e-06.
## Iteration: 530, error: 1.1804e-06.
## Iteration: 540, error: 1.1804e-06.
## Iteration: 550, error: 1.1804e-06.
## Iteration: 560, error: 1.1804e-06.
## Iteration: 570, error: 1.1804e-06.
## Iteration: 580, error: 1.1804e-06.
## Iteration: 590, error: 1.1804e-06.
## Iteration: 600, error: 1.1804e-06.
## Iteration: 610, error: 1.1804e-06.
## Iteration: 620, error: 1.1804e-06.
## Iteration: 630, error: 1.1804e-06.
## Iteration: 640, error: 1.1804e-06.
## Iteration: 650, error: 1.1804e-06.
## Iteration: 660, error: 1.1804e-06.
## Iteration: 670, error: 1.1804e-06.
## Iteration: 680, error: 1.1804e-06.
## Iteration: 690, error: 1.1804e-06.
## Iteration: 700, error: 1.1804e-06.
## Iteration: 710, error: 1.1804e-06.
## Iteration: 720, error: 1.1804e-06.
## Iteration: 730, error: 1.1804e-06.
## Iteration: 740, error: 1.1804e-06.
## Iteration: 750, error: 1.1804e-06.
## Iteration: 760, error: 1.1804e-06.
## Iteration: 770, error: 1.1804e-06.
## Iteration: 780, error: 1.1804e-06.
## Iteration: 790, error: 1.1804e-06.
## Iteration: 800, error: 1.1804e-06.
## Iteration: 810, error: 1.1804e-06.
## Iteration: 820, error: 1.1804e-06.
## Iteration: 830, error: 1.1804e-06.
## Iteration: 840, error: 1.1804e-06.
## Iteration: 850, error: 1.1804e-06.
## Iteration: 860, error: 1.1804e-06.
## Iteration: 870, error: 1.1804e-06.
## Iteration: 880, error: 1.1804e-06.
## Iteration: 890, error: 1.1804e-06.
## Iteration: 900, error: 1.1804e-06.
## Iteration: 910, error: 1.1804e-06.
## Iteration: 920, error: 1.1804e-06.
## Iteration: 930, error: 1.1804e-06.
## Iteration: 940, error: 1.1804e-06.
## Iteration: 950, error: 1.1804e-06.
## Iteration: 960, error: 1.1804e-06.
## Iteration: 970, error: 1.1804e-06.
## Iteration: 980, error: 1.1804e-06.
## Iteration: 990, error: 1.1804e-06.
## Iteration: 1000, error: 1.1804e-06.
## Reached maximum number of iterations (1000). Error=1.1804e-06.Including Plots
You can also embed plots, for example:
Note that the echo = FALSE parameter was added to the
code chunk to prevent printing of the R code that generated the
plot.