#------------------------------------------------------------------------------
#            1. Merge to-gpinter-ize do-file
#------------------------------------------------------------------------------
  # -------- 0. Packages
library(readxl)
library(haven)
library(dplyr)
library(tidyr)
library(writexl)
library(stringr)

project_dir <- "/Users/manuelestebanarias/Desktop/2025_replication_extension_nonmain_mer"

#--------- 1.0 Import current macrodata 

aggregates_hist <- download_wid (
  indicators = c("ahweal","anninc"),
  areas      = c("FR", "ES", "IT", "SE", "OC-PPP", "OK-PPP", "BR", "OE-PPP", "ZA", "OJ-PPP", "CN", "JP", "OB-PPP"),
  year       = 1820:1910,
  perc       = "p0p100",
  ages       = 992,
  pop        = "i")

aggregates_hist <- aggregates_hist %>%
  pivot_wider(
    names_from  = variable,
    values_from = value
  ) %>%
  select(-percentile) %>%
  rename(iso=country) 




#--------- 1.1 Import current data 
current_wid <- download_wid (
  indicators = c("sptinc","aptinc","tptinc"),
  areas      = c("FR", "ES", "IT", "SE", "OC-PPP", "OK-PPP", "BR", "OE-PPP", "ZA", "OJ-PPP", "CN", "JP", "OB-PPP"),
  year       = 1820:1910,
  ages       = 992, 
  pop = "j")

current_wid$variable <- substr(current_wid$variable, 1, 1)


current_wid <-current_wid %>%
  rename(
    p       = percentile,
    iso     = country
  )

current_wid_top10 <- current_wid %>%
  filter(p=="p90p100" & variable=="s")

current_wid <- current_wid %>%
  # Parsing the percentiles
  mutate(
    p_min = NA_real_,
    p_max = NA_real_
  ) %>%
  mutate(
    p_min = if_else(
      str_detect(p, "^p([0-9\\.]+)p([0-9\\.]+)$"),
      round(1000 * as.numeric(str_match(p, "^p([0-9\\.]+)p([0-9\\.]+)$")[,2])),
      p_min
    ),
    p_max = if_else(
      str_detect(p, "^p([0-9\\.]+)p([0-9\\.]+)$"),
      round(1000 * as.numeric(str_match(p, "^p([0-9\\.]+)p([0-9\\.]+)$")[,3])),
      p_max
    )
  ) %>%
  mutate(
    p_min = if_else(
      str_detect(p, "^p([0-9\\.]+)$"),
      round(1000 * as.numeric(str_match(p, "^p([0-9\\.]+)$")[,2])),
      p_min
    ),
    p_max = if_else(is.na(p_max), 1000*100, p_max)
  ) %>%
  mutate(
    p_max = case_when(
      is.na(p_max) & between(p_min, 0, 98000)   ~ p_min + 1000,
      is.na(p_max) & between(p_min, 99000, 99800) ~ p_min + 100,
      is.na(p_max) & between(p_min, 99900, 99980) ~ p_min + 10,
      is.na(p_max) & between(p_min, 99990, 99999) ~ p_min + 1,
      TRUE ~ p_max
    )
  ) %>%
  # Reconstruct p 
  mutate(
    p = if_else(
      !is.na(p_max),
      paste0("p", round(p_min / 1e3, 3), "p", round(p_max / 1e3, 3)),
      p
    )
  ) %>%
  # Keeping the g-percentiles
  mutate(
    n = round(p_max - p_min, 1)
  )%>%
  filter(n %in% c(1, 10, 100, 1000)) %>%
  filter(!(n == 1000 & p_min >= 99000)) %>%
  filter(!(n == 100  & p_min >= 99900)) %>%
  filter(!(n == 10   & p_min >= 99990)) %>%
  # Clean columns
  select(-p_max, -p) %>%
  rename(p = p_min)

current_wid <- current_wid %>%
  pivot_wider(
    names_from  = variable,
    values_from = value
  )%>%
  mutate(name = "historical_sptinc992j")


current_wid <- current_wid %>%
  left_join(aggregates_hist, by = c("iso", "year")) 

current_wid <- current_wid %>%
  mutate(
    a = a/anninc992i,
    t = t/anninc992i
  ) %>%
  select(- ahweal992i,- anninc992i) 

#--------- 2. Calculate coefficient of change
# Extract  currenct data of the WID
input_data <- current_wid_top10%>%
    filter(p == 90000)%>%
    rename(current = value)

index <- read_excel(file.path(project_dir,'Corrections_round1.xlsx'),
                      sheet = 'Corrections', range = 'A4:AA9')

# Rename columns: prefix with 'k'
colnames(index)[1] <- 'year'  

# Reshape from wide to long format
index <- index %>%
  select( year, contains("corrected"))%>%
  pivot_longer(
    cols = - year,
    names_to = "iso",
    values_to = "pc2021"
  )%>%
  mutate(iso = str_replace(iso, " corrected", "")) %>%
  filter(iso != "year") %>%  
  mutate(
    iso = ifelse(iso == "QM", "OK", iso)
  )%>%
  arrange(iso, year)

# Bring the current data
index <- index %>%
  left_join(input_data, by = c("iso", "year"))%>%
  mutate(
    k = pc2021/current
  )%>%
  select(iso, year, pc2021, k )

rm(input_data)

# -------- 3.  Prepare data for gpinter
# Create fractile groups
# Fractile breaks: from 0 to 75000 by 25000, then 90000, 99000, 100000
breaks <- c(seq(0, 75000, 25000), 90000, 99000, 100000)

data <- current_wid %>%
  mutate(
    ftile_index = cut(p, breaks = breaks, include.lowest = TRUE, right = FALSE, labels = FALSE),
    ftile = breaks[ftile_index]  # this gives you the left edge of the interval
  )

# Calculate population
data <- data %>%
  group_by(iso, year) %>%
  arrange(p) %>%
  mutate(
    lead_p = lead(p),
    pop = (lead_p - p) / 1e5,
    pop = ifelse(is.na(pop), 1 / 1e5, pop)
  ) %>%
  ungroup() %>%
  select(-lead_p)

# Collapse data by iso-year-ftile
data <- data %>%
  group_by(iso, year, ftile) %>%
  summarise(
    t = min(t, na.rm = TRUE),
    s = sum(s, na.rm = TRUE),       
    a = sum(a*pop, na.rm = TRUE),
    pop = sum(pop, na.rm = TRUE),
    .groups = "drop"
  )

data <- data %>%
  mutate(
    a = a/pop
  ) 
  
# Verify consistency
data_org <- data %>%
  group_by(iso, year) %>%
  mutate(
    shr = sum(s, na.rm = TRUE),
    apop = a * pop,
    avg = sum(apop, na.rm = TRUE)
  ) %>%
  ungroup()

# Check: average and share should be close to 1
stopifnot(all(abs(data_org$shr - 1) < 1e-2))
stopifnot(all(abs(data_org$avg - 1) < 1e-2))

#  Cleanup
data_org <- data_org %>%
  select(-shr, -avg, -apop) %>%
  rename(p = ftile) %>%
  mutate(
    p =p,
    name = "historical_sptinc992j"
  )

#---------- 5. Finding and applying a new coeff 
# Note: The k corrections dis-adjust the series (the K would keep the series 
# balance only if the share of the top 10 was 50 but given that not, its application
# will make the sum(bracket_avg*pop)!=1). As so correction to this index needs to 
# be applied as well

# Merge correction index 
data <- data_org %>%
  left_join(index, by = c("iso", "year"))

data <- data %>%
  mutate(alpha=NA_real_)

countries <- unique(data$iso)
years     <- unique(data$year)
for ( x in countries) {
  for ( y in years) {
    image <- data%>%
      filter(iso  == x)%>%
      filter(year == y)
    Sopt<- image$pc2021[1]
    S0  <- image$s[1]
    S25 <- image$s[2]
    S50 <- image$s[3]
    S75 <- image$s[4]
    S90 <- image$s[5]
    S99 <- image$s[6]
    k   <- image$k[1]
    
    A <- S0 + S25 + S50 + S75
    B <- S90 + S99
    
    f <- function(lph) {
      (Sopt * ((A / (k + lph)) + B * (k + lph)) / B) - k - lph
    }
    
    lph<-uniroot(f, lower = -k, upper = 2)$root
    
    data <- data %>%
      mutate(alpha = if_else(iso == x & year ==y, lph, alpha))
  }
}
rm(image, A , B, countries, years, k, lph, S0, S25, S50, S75,S90,S99, Sopt, x, y)

data <- data%>%
  mutate(
    k_adj = k+alpha
  )%>%
  select(iso, name, year,p, pop, t, s, a, k_adj)
  

# 2. Multiply variables by k if p > 90000 or >= 90000
data <- data %>%
  mutate(
    t = if_else(p > 90000, t * k_adj, t),
    s = if_else(p >= 90000, s * k_adj, s),
    a = if_else(p >= 90000, a * k_adj, a)
  )

# 3. Divide variables by k if p < 90000
data <- data %>%
  mutate(
    t = if_else(p < 90000, t / k_adj, t),
    s = if_else(p < 90000, s / k_adj, s),
    a = if_else(p < 90000, a / k_adj, a)
  )

# Verify consistency
data <- data %>%
  group_by(iso, year) %>%
  mutate(
    shr = sum(s, na.rm = TRUE),
    apop = a * pop,
    avg = sum(apop, na.rm = TRUE)
  ) %>%
  ungroup()


# adjust
data <- data %>%
  mutate(
    a= a/avg,
    t= t/avg,
    s=s/shr
  )
#  Cleanup
data <- data %>%
  select(-shr, -avg, -apop, -k_adj)


# Verify 2nd consistency
data <- data %>%
  group_by(iso, year) %>%
  mutate(
    shr = sum(s, na.rm = TRUE),
    apop = a * pop,
    avg = sum(apop, na.rm = TRUE)
  ) %>%
  ungroup()

# Check: average and share should be close to 1
stopifnot(all(abs(data$shr - 1) < 1e-2))
stopifnot(all(abs(data$avg - 1) < 1e-2))

#  Cleanup
data <- data %>%
  select(-shr, -avg, -apop)


#------------------------------------------------------------------------------
#            2. Run the Gpinter
#------------------------------------------------------------------------------
library(tidyverse)
library(haven)
library(gpinter)
library(here)
library(progressr)

handlers(handler_progress(
  format = ":spin :current/:total [:bar] :percent in :elapsedfull ETA: :eta"
))
options(dplyr.summarise.inform = FALSE)

gperc <- c(
  seq(0, 99000, 1000), seq(99100, 99900, 100),
  seq(99910, 99990, 10), seq(99991, 99999, 1)
)


groups <- data %>%
  select(iso, year, pop ,p, name, a, t) %>%
  drop_na() %>%
  group_by(iso, year, name) %>%
  mutate(average = sum(a * pop, na.rm = TRUE)) %>%
  arrange(iso, year, name, p) %>%
  group_by(iso, year, name) %>%
  filter(n() > 2)%>%
  group_split()


with_progress({
  p <- progressor(along = groups)
  gpinterized <- groups %>% map_dfr(~ {
    dist <- tabulation_fit(p = .x$p/1e5, threshold = .x$t,  bracketavg = .x$a, average = .x$average[1], fast = TRUE)
    tab <- generate_tabulation(dist, gperc/1e5) %>% as_tibble() %>% transmute(
      iso = first(.x$iso),
      year = first(.x$year),
      name = first(.x$name),
      p = round(fractile*1e5),
      t = threshold,
      ts = top_share,
      bs = bottom_share,
      s = bracket_share,
      a = bracket_average
    )
    p()
    return(tab)
  })
})


write_dta( gpinterized, here(file.path(project_dir,"/output1_round1-gpinterized.dta")))

