### Inequality Analysis of House Sizes ###

# Load required libraries
library(ggplot2)
library(dplyr)
library(ineq)

# ----------------------------------------
Define the dataset (house sizes)
# ----------------------------------------
house_sizes <- c(
  15.57, 13.87, 18.81, 23.76, 19.25, 21.97, 18.88, 25.24, 21.41, 21.61,
  21.81, 25.7, 25.26, 27.41, 24.52, 36.71, 27.71, 27.32, 32.44, 28.49,
  31.04, 29.95, 21.36, 33.76, 34.9, 22.22, 22.81, 38.95, 25.4, 33.95,
  34.64, 27.69, 19.72, 33.58, 29.91, 36.77, 38.6, 30.24, 44.01, 26.06,
  46.23, 35.82, 43.93, 32.29, 34.56, 31.08, 36.81, 24.61, 35.99, 22.96,
  28.4, 38.36, 33.47, 42.97, 32.76, 36.00, 35.6, 36.56, 35.23, 45.89,
  40.05, 23.67, 45.92, 42.34, 31.32, 33.03, 40.68, 44.87, 47.04, 30.22,
  58.83, 46.29, 42.27, 56.97, 63.6, 36.77, 37.65, 33.98, 42.58, 72.13,
  34.98, 24.84, 42.51, 49.59, 36.96, 45.98, 37.63, 39.6, 42.98, 39.43,
  54.03, 43.16, 37.76, 37.17, 36.58, 36.36, 37.94, 35.19, 31.78, 37.94,
  52.59, 41.02, 53.34, 42.75, 40.16, 35.96, 49.75, 64.08, 50.15, 55.51,
  54.45, 57.44, 65.91, 40.94, 46.18, 39.47, 53.34, 52.94, 44.62, 45.23,
  40.69, 43.96, 41.18, 34.19, 44.19, 53.09, 54.34, 46.84, 37.74, 45.05,
  52.66, 71.09, 69.06, 43.01, 49.21, 56.29, 54.52, 50.13, 45.22, 62.08,
  55.68, 36.26, 61.39, 80.43, 55.25, 38.59, 62.75, 40.55, 43.44, 51.87,
  52.94, 53.32, 50.1, 44.17, 48.61, 46.69, 55.02, 32.33, 73.87, 46.12,
  47.73, 40.92, 45.23, 39.58, 48.52, 65.51, 62.55, 70.54, 56.84, 40.06,
  53.83, 88.41, 60.27, 53.48, 42.41, 65.2, 68.79, 49.18, 57.25, 65.05,
  70.26, 66.28, 59.42, 46.65, 42.09, 68.32, 88.5, 48.9, 53.64, 51.89,
  82.48, 56.91, 57.33, 66.9, 83.91, 92.93, 100.79, 70.26, 49.22, 97.27,
  54.00, 111.41, 104.71, 84.37, 61.72, 69.59, 46.04, 56.55, 72.55, 68.14,
  92.68, 66.42, 123.53, 104.79, 110.97, 118.99, 90.97
)

# ----------------------------------------
Step 1: Estimate required number of runs for stable Gini estimate
# ----------------------------------------

set.seed(123)
min_size <- 10
max_size <- 100
max_runs <- 10000
tolerance <- 0.01
n_repeats <- 1000

runs_needed <- numeric(n_repeats)
montecarlo_data <- list()

for (r in 1:n_repeats) {
  metric_values <- numeric(max_runs)
  running_mean <- numeric(max_runs)
  running_se <- numeric(max_runs)
  
  for (i in 1:max_runs) {
    sample_size <- sample(min_size:max_size, 1)
    sample_i <- sample(house_sizes, sample_size, replace = TRUE)
    metric_values[i] <- Gini(sample_i)
    
    running_mean[i] <- mean(metric_values[1:i])
    running_se[i] <- sd(metric_values[1:i]) / sqrt(i)
    
    if (i > 30) {
      rel_error <- running_se[i] / running_mean[i]
      if (rel_error < tolerance) {
        runs_needed[r] <- i
        break
      }
    }
    if (i == max_runs) {
      runs_needed[r] <- max_runs
    }
  }
  
  # Store data for first 20 repetitions for plotting
  if (r <= 20) {
    df <- data.frame(run = 1:runs_needed[r],
                     mean_gini = running_mean[1:runs_needed[r]],
                     repetition = as.factor(r))
    montecarlo_data[[r]] <- df
  }
}

# Combine data into a single data frame
montecarlo_df <- bind_rows(montecarlo_data)

# ----------------------------------------
Plot Monte Carlo analysis
# ----------------------------------------

# Plot 1: Histogram of how many runs were needed for convergence
mean_runs <- mean(runs_needed)

ggplot(data.frame(runs_needed), aes(x = runs_needed)) +
  geom_histogram(binwidth = 10, fill = "grey", color = "black") +
  geom_vline(xintercept = mean_runs, color = "red", linetype = "dashed", size = 1) +
  labs(title = "",
       x = "Number of Runs", y = "Frequency") +
  theme_bw()

# Plot 2: Running mean Gini index for 20 example repetitions
ggplot(montecarlo_df, aes(x = run, y = mean_gini, group = repetition)) +
  geom_line(color = "grey50", size = 0.7) +
  labs(title = "",
       x = "Number of Runs", y = "Running Mean Gini Index") +
  theme_bw()

# ----------------------------------------
Step 2: Monte Carlo Simulation of 500 Random Gini Samples
# ----------------------------------------

n_samples <- 500
gini_values <- numeric(n_samples)

for (i in 1:n_samples) {
  sample_size <- sample(10:100, 1)
  sample_i <- sample(house_sizes, sample_size, replace = FALSE)
  gini_values[i] <- Gini(sample_i)
}

# Plot 3: Histogram of Gini values
ggplot(data.frame(Gini = gini_values), aes(x = Gini)) +
  geom_histogram(binwidth = 0.005, fill = "grey", color = "black") +
  labs(title = "Distribution of Gini Indices (Random Samples)",
       x = "Gini Index", y = "Frequency") +
  theme_bw()

# Output 5-number summary of Gini results
summary(gini_values)