# Load necessary packages
library(sf); library(ggplot2); library(dplyr)
# Load in the Netherlands Roman Empire data
netherlands_re <- st_read("Data/netherlands_roman_updated.geojson")## Reading layer `netherlands_roman_updated' from data source
## `C:\Users\Milan\OneDrive - Universiteit Utrecht\FINAL thesis code\Data\netherlands_roman_updated.geojson'
## using driver `GeoJSON'
## Simple feature collection with 355 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 3.360782 ymin: 50.72349 xmax: 7.227095 ymax: 53.55458
## Geodetic CRS: WGS 84# Convert the 'is_in_re' column to a factor
netherlands_re$is_in_re <- as.factor(netherlands_re$is_in_re)
# Filter out water observations from municipalities
excluded_areas <- c("IJsselmeer", "Zeeuwse meren", "Zuid Hollandse Meren")
netherlands_re <- netherlands_re %>%
filter(!(NAME_1 %in% excluded_areas))
# Create the map
titlepage_plot <- ggplot(netherlands_re) +
geom_sf(aes(fill = is_in_re), color = "black", size = 0.1) +
scale_fill_manual(
values = c("1" = "#800000", "0" = "white"),
labels = c("0" = "No", "1" = "Yes"),
name = "Roman Empire"
) +
labs(title = "Dutch Municipalities & the Roman Empire") +
theme_minimal() +
theme(legend.position = "bottom",
panel.background = element_rect(fill = "lightblue"),
plot.background = element_rect(fill = "lightblue"))
# Save the plot
ggsave("Figures/titlepage.png", plot = titlepage_plot, width = 6, height = 7, dpi = 300)# Load additional packages
library(tidyverse); library(geodata)
# Roman Empire
re <- cawd::awmc.roman.empire.117.sp |>
st_as_sf()
# Funda house prices with coordinates
funda_data <- read.csv("Data/houses_funda.csv", sep = ";")
funda_data$postcode_4 <- as.character(funda_data$postcode_4)
geopostcode <- st_read("Data/georef-netherlands-postcode-pc4.geojson")## Reading layer `georef-netherlands-postcode-pc4' from data source
## `C:\Users\Milan\OneDrive - Universiteit Utrecht\FINAL thesis code\Data\georef-netherlands-postcode-pc4.geojson'
## using driver `GeoJSON'
## Simple feature collection with 4068 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 3.358524 ymin: 50.75119 xmax: 7.228207 ymax: 53.5186
## Geodetic CRS: WGS 84## Merge the datasets based on postal codes
funda_geo <- funda_data |>
left_join(geopostcode, by = c("postcode_4" = "pc4_code"))
# Convert funda_geo to an sf object
funda_geo <- funda_geo |>
st_as_sf()
# Is the house inside the roman empire?
is_in_re <- st_intersects(funda_geo, re)
is_in_re <- is_in_re |> as.numeric()
# Compute whether house is in empire
funda_geo <- funda_geo |>
mutate(is_in_re = if_else(is.na(is_in_re), 0, is_in_re))
# Compute border
boundary <- re |> st_boundary() |> st_intersection(funda_geo)
# Define the coordinates of the bounding box
coords <- matrix(c(4.6, 52.5,
4.6, 52.6,
6.4, 52.3,
6.4, 51.7,
4.6, 52.5), ncol = 2, byrow = TRUE)
# Convert LINESTRING to POLYGON and create an sf object
bounding_box_polygon <- st_polygon(list(coords))
bounding_box_sf <- st_sf(geometry = st_sfc(bounding_box_polygon))
# Set the Coordinate Reference System (CRS)
st_crs(bounding_box_sf) <- 4326
# Create final boundary
final_boundary <- st_intersection(boundary |> st_union(), bounding_box_sf)
# Distance to border:
distances <- funda_geo |>
st_centroid() |>
st_distance(final_boundary)
# Add to df
funda_geo <- funda_geo |>
mutate(distance_to_border = as.numeric(if_else(is_in_re == 1, distances, - distances)))
# Add geographical control variables
waterbodiesnl <- st_read("Data/KRW_waterbodies.geojson")## Reading layer `KRW_waterbodies' from data source
## `C:\Users\Milan\OneDrive - Universiteit Utrecht\FINAL thesis code\Data\KRW_waterbodies.geojson'
## using driver `GeoJSON'
## Simple feature collection with 54 features and 40 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -4962.57 ymin: 307330.2 xmax: 276801.2 ymax: 642207
## Projected CRS: Amersfoort / RD New# Ensure both datasets are in the same CRS
waterbodiesnl <- st_transform(waterbodiesnl, st_crs(funda_geo))
# Calculate the distance from each house to the nearest water body
dist_to_water <- st_distance(st_centroid(funda_geo), waterbodiesnl)
# Extract the minimum distance for each house
min_dist_to_water <- apply(dist_to_water, 1, min)
# Add the distance to the funda_geo dataset as a new column
funda_geo <- funda_geo |>
mutate(distance_to_water = as.numeric(min_dist_to_water))
# Load ancient city data
ancientcities <- read.csv("Data/Hanson2016_Cities_OxREP.csv")
ancientcities_sf <- st_as_sf(ancientcities, coords = c("Longitude..X.", "Latitude..Y."), crs = 4326)
# Ensure both datasets are in the same CRS
ancientcities_sf <- st_transform(ancientcities_sf, st_crs(funda_geo))
# Calculate the distance from each house to the nearest ancient city
dist_to_ancient_city <- st_distance(st_centroid(funda_geo), ancientcities_sf)
# Extract the minimum distance for each house
min_dist_to_ancient_city <- apply(dist_to_ancient_city, 1, min)
# Add the distance to the funda_geo dataset as a new column
funda_geo <- funda_geo |>
mutate(distance_to_ancient_city = as.numeric(min_dist_to_ancient_city))
# Clean column names
names(funda_geo) <- make.names(names(funda_geo))
# Export to geojson
funda_geo |> write_sf("main_data.geojson")## Reading layer `main_data' from data source
## `C:\Users\Milan\OneDrive - Universiteit Utrecht\FINAL thesis code\main_data.geojson'
## using driver `GeoJSON'
## Simple feature collection with 40020 features and 21 fields (with 4 geometries empty)
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 3.358524 ymin: 50.75559 xmax: 7.228207 ymax: 53.47275
## Geodetic CRS: WGS 84# Remove observations with NA values
main_data <- na.omit(main_data)
# Remove unnecessary variables
main_data <- main_data %>% select(-prices_per_sq_m, -geo_point_2d, -year, -gem_code, -prov_code)
# Remove the municipality of Amsterdam and province of Flevoland observations
main_data <- main_data %>% filter(gem_name != "Amsterdam", prov_name != "Flevoland")# Convert energy labels to numeric for the analysis. Making the highest energy label
# also the highest number
energy_levels <- c("G", "F", "E", "D", "C", "B", "A", "A+", "A++", "A+++", "A++++")
main_data$energy_label <- factor(main_data$energy_label, levels = energy_levels, ordered = TRUE)
main_data$energy_label <- as.numeric(main_data$energy_label)# Convert prices to numeric and correcting values of 1 million + that are incorrectly
# extracted from funda as 1,000 instead of 1,000,000
main_data$prices <- as.numeric(gsub(",", "", main_data$prices))
main_data$prices <- sapply(main_data$prices, function(x) {
if (x < 20000) {
return(x * 1000)
} else {
return(x)
}
})# Converting the footage_interior and footage_exterior variables to numeric and square meters
# Also correcting values of 1000 m^2 + that are incorrectly extracted from funda
main_data$footage_interior <- as.numeric(gsub(",", "", main_data$footage_interior))
main_data$footage_exterior <- as.numeric(gsub(",", "", main_data$footage_exterior))
# Removing 1 observation that would incorrectly be corrected along with
# observations with 1000 m^2+ values
min_footage_interior <- min(main_data$footage_interior, na.rm = TRUE)
main_data <- main_data %>% filter(footage_interior != min_footage_interior)
# Correcting values
main_data$footage_exterior <- sapply(main_data$footage_exterior, function(x) {
if (x < 25000) {
return(x * 1000)
} else {
return(x)
}
})
# Convert to square meters and rename variables
main_data <- main_data %>%
mutate(
footage_interior = footage_interior / 1000,
footage_exterior = footage_exterior / 1000
) %>%
rename(
sqm_interior = footage_interior,
sqm_exterior = footage_exterior
)# Load additional packages
library(stargazer);library(knitr);library(kableExtra);library(summarytools)
library(data.table); library(webshot); library(scales); library(gt)# Convert data frame to data.table
main_data_dt <- as.data.table(main_data)
# Compute descriptive statistics
descriptive_stats <- main_data_dt[, .(
N = .N,
Mean = mean(prices, na.rm = TRUE),
SD = sd(prices, na.rm = TRUE),
Min = min(prices, na.rm = TRUE),
Max = max(prices, na.rm = TRUE)
), by = is_in_re]
# Format the numeric values with commas for thousands
descriptive_stats[, Mean := comma(Mean)]
descriptive_stats[, SD := comma(SD)]
descriptive_stats[, Min := comma(Min)]
descriptive_stats[, Max := comma(Max)]
descriptive_stats[, N := comma(N)]
# Create the table
descriptive_table_gt <- gt(descriptive_stats) %>%
tab_header(
title = "Table 1 - Descriptive Statistics"
) %>%
cols_label(
is_in_re = "Is in Roman Empire",
N = "N",
Mean = "Mean",
SD = "SD",
Min = "Min",
Max = "Max"
) %>%
fmt_number(
columns = vars(N, Mean, SD, Min, Max),
decimals = 2,
use_seps = TRUE
) %>%
tab_source_note(
source_note = "Prices are in euros"
) %>%
cols_align(
align = "center",
columns = everything()
) %>%
tab_options(
table.font.size = "small",
table.font.names = "Arial"
)
# Save the table as an HTML file
gtsave(descriptive_table_gt, "Tables/descstats/desc_stats_table_gt.html")# Linear regression model 1 with only house characteristics as control variables
lrm1 <- lm(prices ~ is_in_re + sqm_interior + sqm_exterior + bedrooms + energy_label, data = main_data)
# Adding geographical control variables to linear regression model 1
lrm2 <- lm(prices ~ is_in_re + sqm_interior + sqm_exterior + bedrooms + energy_label
+ distance_to_border + distance_to_water + distance_to_ancient_city, data = main_data)
# Combine both models into one table and save as HTML
stargazer(lrm1, lrm2, type = "html", title = "Table 2 - Linear Regression",
column.labels = c("LR Model 1", "LR Model 2"),
covariate.labels = c("In Roman Empire", "Interior Size (sqm)", "Exterior Size (sqm)",
"Number of Bedrooms", "Energy Label", "Distance to Border",
"Distance to Water", "Distance to Ancient City"),
dep.var.caption = "Dependent Variable: Housing Prices",
out = "Tables/LRM/lrm1.2_table.html")# Checking for multicolinearity because of negative "bedrooms" coefficient
library(car)
vif2_values <- vif(lrm2)
print(vif2_values)## is_in_re sqm_interior sqm_exterior
## 3.073957 1.881289 1.207393
## bedrooms energy_label distance_to_border
## 1.636848 1.025462 3.167364
## distance_to_water distance_to_ancient_city
## 1.609852 1.958498# Adding interaction terms to lrm 2
lrm3 <- lm(prices ~ is_in_re * sqm_interior + sqm_exterior + bedrooms + energy_label
+ distance_to_border + distance_to_water + distance_to_ancient_city, data = main_data)
# Adding a log transformation for the dependent variable "prices"
lrm4 <- lm(log(prices) ~ is_in_re * sqm_interior + sqm_exterior + bedrooms + energy_label
+ distance_to_border + distance_to_water + distance_to_ancient_city, data = main_data)
# Correct labels based on model summaries
lrmrbc_covariate_labels <- c("In Roman Empire", "Interior Size (sqm)",
"Exterior Size (sqm)", "Number of Bedrooms", "Energy Label",
"Distance to Border", "Distance to Water", "Distance to Ancient City",
"In Roman Empire * Interior Size")
# Create a table for both models and save as HTML
stargazer(lrm3, lrm4, type = "html", title = "Table 3 - Linear Regression Robustness Checks",
column.labels = c("LR Model 3", "LR Model 4"),
covariate.labels = lrmrbc_covariate_labels,
dep.var.labels.include = FALSE,
model.names = FALSE,
header = FALSE,
add.lines = list(c("Dependent Variable", "Housing Prices", "log(Housing Prices)")),
out = "Tables/LRM/lrm3_4_table.html")# Run the balance checks
# Size Interior
rdd_balance1 <- rdrobust(y = main_data$sqm_interior, x = main_data$distance_to_border, c = 0)
# Size Exterior
rdd_balance2 <- rdrobust(y = main_data$sqm_exterior, x = main_data$distance_to_border, c = 0)
# Bedrooms
rdd_balance3 <- rdrobust(y = main_data$bedrooms, x = main_data$distance_to_border, c = 0)
# Energy Label
rdd_balance4 <- rdrobust(y = main_data$energy_label, x = main_data$distance_to_border, c = 0)
# Function to extract p-value and coefficients
balance_coefs <- function(rdd_result) {
coef <- rdd_result$coef[1, 1]
se <- rdd_result$se[1, 1]
p_value <- rdd_result$pv[1, 1]
stars <- if (p_value < 0.01) {
"***"
} else if (p_value < 0.05) {
"**"
} else if (p_value < 0.1) {
"*"
} else {
""
}
return(list(coef = paste0(round(coef, 3), stars), se = round(se, 3)))
}
# Extracting the necessary statistics
balance_results <- data.frame(
Variable = c("Size Interior", "Size Exterior", "Bedrooms", "Energy Label"),
Coefficient = c(balance_coefs(rdd_balance1)$coef,
balance_coefs(rdd_balance2)$coef,
balance_coefs(rdd_balance3)$coef,
balance_coefs(rdd_balance4)$coef),
SE = c(balance_coefs(rdd_balance1)$se,
balance_coefs(rdd_balance2)$se,
balance_coefs(rdd_balance3)$se,
balance_coefs(rdd_balance4)$se),
Control = c(rdd_balance1$N_h[1], rdd_balance2$N_h[1], rdd_balance3$N_h[1], rdd_balance4$N_h[1]),
Treated = c(rdd_balance1$N_h[2], rdd_balance2$N_h[2], rdd_balance3$N_h[2], rdd_balance4$N_h[2]),
Bandwidth = c(rdd_balance1$bws[1, 1], rdd_balance2$bws[1, 1], rdd_balance3$bws[1, 1], rdd_balance4$bws[1, 1])
)
# Creating the table
balance_table <- gt(balance_results) %>%
tab_header(title = "Table 4 - Balance Check Results") %>%
tab_source_note(
source_note = "Statistical Significance: *** p<0.01, ** p<0.05, * p<0.1"
)
# Saving the table as an HTML file
gtsave(balance_table, "Tables/RDD/balance_check_results.html")# RDD Model 1
rdd1 <- rdrobust(y = main_data$prices, x = main_data$distance_to_border, c = 0)
# RDD Model 2
rdd2 <- rdrobust(y = main_data$prices, x = main_data$distance_to_border, c = 0,
covs=cbind(main_data$sqm_interior, main_data$sqm_exterior,
main_data$bedrooms, main_data$energy_label))
# RDD Model 3
rdd3 <- rdrobust(y = main_data$prices, x = main_data$distance_to_border, c = 0,
covs=cbind(main_data$sqm_interior, main_data$sqm_exterior,
main_data$bedrooms, main_data$energy_label,
main_data$distance_to_water, main_data$distance_to_ancient_city))
# Function to extract p-value and coefficients
rdd_coefs <- function(rdd_result) {
coef <- rdd_result$coef[1, 1]
se <- rdd_result$se[1, 1]
p_value <- rdd_result$pv[1, 1]
stars <- if (p_value < 0.01) {
"***"
} else if (p_value < 0.05) {
"**"
} else if (p_value < 0.1) {
"*"
} else {
""
}
return(list(coef = paste0(round(coef, 3), stars), se = round(se, 3)))
}
# Extracting the necessary statistics
rdd_results <- data.frame(
Model = c("RDD Model 1", "RDD Model 2", "RDD Model 3"),
Coefficient = c(rdd_coefs(rdd1)$coef,
rdd_coefs(rdd2)$coef,
rdd_coefs(rdd3)$coef),
SE = c(rdd_coefs(rdd1)$se,
rdd_coefs(rdd2)$se,
rdd_coefs(rdd3)$se),
Control = c(rdd1$N_h[1], rdd2$N_h[1], rdd3$N_h[1]),
Treated = c(rdd1$N_h[2], rdd2$N_h[2], rdd3$N_h[2]),
Bandwidth = c(rdd1$bws[1, 1], rdd2$bws[1, 1], rdd3$bws[1, 1])
)
# Creating the table
rdd_table <- gt(rdd_results) %>%
tab_header(title = "Table 5 - RDD Results") %>%
tab_source_note(
source_note = "Statistical Significance: *** p<0.01, ** p<0.05, * p<0.1"
)
# Saving the table as an HTML file
gtsave(rdd_table, "Tables/RDD/rdd_123.html")# RDD Model 4
rdd4 <- rdrobust(y = log(main_data$prices), x = main_data$distance_to_border, c = 0)
# RDD Model 5
rdd5 <- rdrobust(y = log(main_data$prices), x = main_data$distance_to_border, c = 0,
covs=cbind(main_data$sqm_interior, main_data$sqm_exterior,
main_data$bedrooms, main_data$energy_label))
# RDD Model 6
rdd6 <- rdrobust(y = log(main_data$prices), x = main_data$distance_to_border, c = 0,
covs=cbind(main_data$sqm_interior, main_data$sqm_exterior,
main_data$bedrooms, main_data$energy_label,
main_data$distance_to_water, main_data$distance_to_ancient_city))
# Extracting the necessary statistics for RDD results
rdd_log <- data.frame(
Model = c("RDD Model 4", "RDD Model 5", "RDD Model 6"),
Coefficient = c(rdd_coefs(rdd4)$coef,
rdd_coefs(rdd5)$coef,
rdd_coefs(rdd6)$coef),
SE = c(rdd_coefs(rdd4)$se,
rdd_coefs(rdd5)$se,
rdd_coefs(rdd6)$se),
Control = c(rdd4$N_h[1], rdd5$N_h[1], rdd6$N_h[1]),
Treated = c(rdd4$N_h[2], rdd5$N_h[2], rdd6$N_h[2]),
Bandwidth = c(rdd4$bws[1, 1], rdd5$bws[1, 1], rdd6$bws[1, 1])
)
# Creating the table for RDD results
rdd_table2 <- gt(rdd_log) %>%
tab_header(title = "Table 6 - RDD Results (Log-Transformed Prices)") %>%
tab_source_note(
source_note = "Statistical Significance: *** p<0.01, ** p<0.05, * p<0.1"
)
# Saving the table as an HTML file
gtsave(rdd_table2, "Tables/RDD/rdd_456.html")# RDD models 4, 5 & 6 with cerrd bandwidth selection
rdd4_cerrd <- rdrobust(y = log(main_data$prices), x = main_data$distance_to_border, c = 0, bwselect = "cerrd")
rdd5_cerrd <- rdrobust(y = log(main_data$prices), x = main_data$distance_to_border, c = 0,
covs = cbind(main_data$sqm_interior, main_data$sqm_exterior,
main_data$bedrooms, main_data$energy_label), bwselect = "cerrd")
rdd6_cerrd <- rdrobust(y = log(main_data$prices), x = main_data$distance_to_border, c = 0,
covs = cbind(main_data$sqm_interior, main_data$sqm_exterior,
main_data$bedrooms, main_data$energy_label,
main_data$distance_to_water, main_data$distance_to_ancient_city), bwselect = "cerrd")
# Extracting the necessary statistics for RDD results with cerrd bandwidth selection
rdd_results_cerrd <- data.frame(
Model = c("RDD Model 4 (cerrd)", "RDD Model 5 (cerrd)", "RDD Model 6 (cerrd)"),
Coefficient = c(rdd_coefs(rdd4_cerrd)$coef,
rdd_coefs(rdd5_cerrd)$coef,
rdd_coefs(rdd6_cerrd)$coef),
SE = c(rdd_coefs(rdd4_cerrd)$se,
rdd_coefs(rdd5_cerrd)$se,
rdd_coefs(rdd6_cerrd)$se),
Control = c(rdd4_cerrd$N_h[1], rdd5_cerrd$N_h[1], rdd6_cerrd$N_h[1]),
Treated = c(rdd4_cerrd$N_h[2], rdd5_cerrd$N_h[2], rdd6_cerrd$N_h[2]),
Bandwidth = c(rdd4_cerrd$bws[1, 1], rdd5_cerrd$bws[1, 1], rdd6_cerrd$bws[1, 1])
)
# Creating the table for RDD results with cerrd bandwidth selection
rdd_table_cerrd <- gt(rdd_results_cerrd) %>%
tab_header(title = "Table 7 - RDD Results (With cerrd Bandwidth Selection)") %>%
tab_source_note(
source_note = "Statistical Significance: *** p<0.01, ** p<0.05, * p<0.1"
)
# Saving the table as an HTML file
gtsave(rdd_table_cerrd, "Tables/RDD/rdd_cerrd_456.html")