#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
#            1. Correct data 
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
  
# Objective: Find the motive for doing this code at the bottom*

# -------- 0. Packages ---------------------------------------------------------
library(readxl)
library(haven)
library(dplyr)
library(tidyr)
library(writexl)
library(stringr)
library(purrr)

sum_or_na <- function(x) {
  if (all(is.na(x))) NA_real_ else sum(x, na.rm = TRUE)
}


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("OJ-PPP", "OE-PPP", "US","SE", "GB", "JP", "FR", "JP","CA", "CN"),
  year       = 1820:1970,
  perc       = "p0p100",
  ages       = 992,
  pop        = "i")

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

#--------- 1. Import current data (already gpintereised) -----------------------
# Keep iso-year-widcode to be changed
current_wid <- download_wid (
  indicators = c("sptinc","aptinc","tptinc"),
  areas      = c("OJ-PPP", "OE-PPP", "US","SE", "GB", "JP", "FR", "JP","CA", "CN"),
  year       = 1820:1970,
  ages       = 992, 
  pop = "j") %>%
  rename(iso=country) %>%
  filter(!(iso == "US" & year >= 1920))

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

# Extract  currenct data of the WID for the section 2
input_data <- current_wid %>%
  filter(percentile %in% c(
    "p0p10", "p0p25", "p0p50", "p0p75", "p0p90", "p0p99",
    "p10p100", "p25p100", "p50p100", "p75p100", "p90p100", "p99p100"
  )) %>%
  mutate(
   variable = case_when(
      str_detect(percentile, "p100") & variable!="t" ~ paste0("t", variable),
      str_detect(percentile, "^p0")  & variable!="t" ~ paste0("b", variable),
      TRUE                           ~ variable
    )
  ) %>%
  pivot_wider(
    names_from = variable,
    values_from = value
  ) %>%
  rename(p=percentile) %>%
  select(iso, year, p, t, ts, bs) %>%
  rename_with(~ paste0("curr_", .x), c(t, ts, bs)) %>%
  mutate(
    p = case_when(
      p == "p0p10"   ~ 10000,
      p == "p0p25"   ~ 15000,
      p == "p0p50"   ~ 50000,
      p == "p0p75"   ~ 75000,
      p == "p0p90"   ~ 90000,
      p == "p0p99"   ~ 99000,
      p == "p10p100" ~ 10000,
      p == "p25p100" ~ 25000,
      p == "p50p100" ~ 50000,
      p == "p75p100" ~ 75000,
      p == "p90p100" ~ 90000,
      p == "p99p100" ~ 99000,
      TRUE ~ NA_real_
    )
  )



current_wid <-current_wid %>%
  rename(
    p       = percentile
  )

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  correction indexes
index <- read_excel(file.path(project_dir,'Corrections_round2.xlsx'),
                      sheet = 'CorrectionsSeptember2025fix', range = 'A1:Bk13')

# 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 = "corrected"
  )%>%
  mutate(iso = str_replace(iso, " corrected", "")) %>%
  filter(iso != "year") %>%  
  mutate(
    iso = ifelse(iso == "QM", "OK", iso)
  )%>%
  mutate(iso = str_replace_all(iso, " ", ""))%>%
  mutate(
    p   = substr(iso, 3, nchar(iso)),
    iso = substr(iso, 1, 2)
  )%>%
  mutate(
    onelet = if_else(str_detect(p, "P"), "corr_t",if_else(str_detect(p, "B"), "corr_bs","corr_ts"))
  )%>%
  arrange(iso, year, p, onelet)

# Set up the percentiles
index <- index %>%
  pivot_wider(
    names_from  = onelet,
    values_from = corrected
    ) %>%
  mutate(
    kind = if_else(str_detect(p, "B"), "B", if_else(str_detect(p, "T"), "T",""))
  )%>%
  mutate(
    p = as.numeric(str_extract(p, "\\d+")) 
  )%>%
  mutate(
    p = if_else(kind == "T", (100 - p)*1000, p*1000)
  )%>%
  select(-kind)

# Bring the current data
index <- index %>%
  left_join(input_data, by = c("iso", "year","p"))%>%
  mutate(
    k_bs = corr_bs / curr_bs,
    k_ts = corr_ts / curr_ts,
    k_t = corr_t / curr_t
  )%>%
  select(iso, year, p, corr_bs,corr_ts, corr_t , k_bs, k_ts, k_t )%>%  # k_t
  filter(!if_all(c(corr_bs, corr_ts, corr_t), is.na))%>% #k_t
  mutate(
    kind = if_else(!is.na(corr_ts), "ts", if_else(!is.na(corr_bs), "bs", "t"))
  )%>%
  arrange(iso, year, p)

# generate information for the cases
index <- index %>%
  filter(!if_all(c(corr_bs, corr_ts, corr_t), is.na)) %>%
  add_count(iso, year, name = "count") 

index_aux <- index %>%
  group_by(iso, year) %>%
  summarise(
    kinds_available = paste(sort(unique(kind)), collapse = ", "),
    .groups = "drop"
  )
index <- index %>%
    inner_join( index_aux, by = c("iso", "year"))

p_t_available <- unique(index$p[index$kinds_available == "t"])


rm(index_aux)

# -------- 3.  Prepare data for gpinter ----------------------------------------
# Create fractile groups
breaks <- c(0,10000,25000,50000,75000, 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 %>%
  rename(p = ftile) %>%
  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) %>%
  mutate(
    name = "historical_sptinc992j"
  )


#---------- 4. Method 1: Single corrections of inner shares -------------------------------------------------------------------------------------
# objetive :Finding and applying a new coeff  when only one correction is needed (indexone)
# 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

# Generate index for files with Only 1 correction
  indexone <- index %>%
  filter(count == 1 & p!=99000 & p!=0) %>% # We cannot apply the methods to the border
  filter(kinds_available!="t") %>%
  select(iso, year, p, corr_ts, k_ts, corr_bs, k_bs) %>%
  mutate(
    corr_ts = if_else(is.na(corr_ts), corr_bs, corr_ts),
    kind    = if_else(is.na(k_ts),  "bs",  "ts"),
    k_ts    = if_else(is.na(k_ts),  k_bs,  k_ts), ### in the bottom correction we want to adjust the other side of the benchmark p
    
  ) %>%
  select(iso, year, p, k_ts, corr_ts, kind) %>%
  rename(
    k         = k_ts,
    corrected = corr_ts
  ) %>%
  rename(p_k=p)

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

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

countries <- unique(data$iso)

for ( x in countries) {
  image0 <- data%>%
    filter(iso  == x ) 
  
  years     <- unique(image0$year)
  
  for ( y in years) {
    image <- data%>%
      filter( iso  == x & year == y)

    Sopt <- image$corrected[1]
    k    <- image$k[1]
    p_k  <- image$p_k[1]
    kind <- image$kind[1]
    
    A <- sum(image$s[ image$p <  p_k ])   # all s with p <  p_k
    B <- sum(image$s[ image$p >= p_k ])   # all s with p >= p_k
    
    if (kind == "ts") {
      f <- function(lph) {
        (Sopt * ((A / (k + lph)) + B * (k + lph)) / B) - k - lph
      }
    } else {
      f <- function(lph) {
        (Sopt * ((A * (k + lph)) + B / (k + lph)) / A) - k - lph
      }
    }
    lph<-uniroot(f, lower = -k, upper = 2)$root
    
    data <- data %>%
      mutate(alpha = if_else(iso == x & year ==y, lph, alpha))
  }
}
rm(image0, image, A , B, countries, years, k, lph, Sopt, x, y)

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


# 2. Multiply variables by k if p > 90000 or >= 90000 and Divide variables by k if p < 90000
data <- data %>%
  mutate(
    t = case_when(
      kind == "ts" & p >  p_k ~ t * k_adj,
      kind == "ts" & p <  p_k ~ t / k_adj,
      kind == "bs" & p >  p_k ~ t / k_adj,
      kind == "bs" & p <  p_k ~ t * k_adj,
      TRUE ~ t
    ),
    s = case_when(
      kind == "ts" & p >= p_k ~ s * k_adj,
      kind == "ts" & p <  p_k ~ s / k_adj,
      kind == "bs" & p >= p_k ~ s / k_adj,
      kind == "bs" & p <  p_k ~ s * k_adj,
      TRUE ~ s
    ),
    a = case_when(
      kind == "ts" & p >= p_k ~ a * k_adj,
      kind == "ts" & p <  p_k ~ a / k_adj,
      kind == "bs" & p >= p_k ~ a / k_adj,
      kind == "bs" & p <  p_k ~ a * k_adj,
      TRUE ~ 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-6))
stopifnot(all(abs(data$avg - 1) < 1e-6))

#  Cleanup
dataone <- data %>%
  select(iso, name, year, p, pop, t, s, a)
 # select(-shr, -avg, -apop)

#---------- 5. Method2: single and compiled border changes ----------------------------------------------------------------------------

# Generate index kinds
indextwo <- index %>%
  filter(kinds_available!="t") %>%
  filter(count!=1 | p==99000 | p==10000) %>%
  mutate(
    p=if_else(p==10000,     0,  p), # bottom shares belong to the break below
    p=if_else(p==25000, 10000,  p),
    p=if_else(p==50000, 25000,  p)
  )

#---------- 51. Method2: single border changes 
# Single border shares
indextwo1 <- indextwo %>%
  filter(count==1 &( p==99000 | p==10000)) %>%
  rename( p_k=p) %>%
  mutate(
    kind = if_else(!is.na(corr_ts), "ts", "bs")
  )

# Merge correction index 
data <- data_org %>%
  inner_join(indextwo1, by = c("iso", "year")) %>%
  mutate(
    s = if_else(!is.na(corr_bs) & p==p_k, corr_bs, s),
    s = if_else(!is.na(corr_ts) & p==p_k, corr_ts, s),
    t = if_else(!is.na(corr_bs) & p==p_k, t*k_bs, if_else(!is.na(corr_ts) & p==p_k, t*k_ts, t))
  ) 

# re-weight the rest of the shares and correct under the same value the  trehsholds
data <- data %>%
  group_by(iso, year) %>%      # <-- adjust grouping as needed
  mutate(
    rest_s = if_else(
      kind == "ts",
      sum(s[p != 99000], na.rm = TRUE),   # when kind == "ts" (so the correction at top1)
      sum(s[p != 0],    na.rm = TRUE)     # when kind != "ts" (so the correction at bot10)
    )
  ) %>%
  ungroup() %>%
  mutate(
    adj_ratio=  if_else(!is.na(corr_bs), rest_s/(1-corr_bs), rest_s/(1-corr_ts)),
    s = if_else(
      kind == "ts"  & p != 99000,
      s / adj_ratio,
      if_else(kind != "ts" & p != 0,
              s / adj_ratio,
              s)    # otherwise leave s unchanged
    ),
    t = if_else(
      kind == "ts"  & p != 99000,
      t / adj_ratio,
      if_else(kind != "ts" & p != 0,
              t / adj_ratio,
              t)    # otherwise leave s unchanged
    )
  )
  
# Correct the rest of the distribution
data <- data %>%
  mutate(
    a= (s*1)/pop
  ) %>%
  select(-p_k, -corr_bs, -corr_ts, -corr_t, -kind, -rest_s, -adj_ratio, -k_bs, -k_ts, -count)
  
# 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()

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

#  Cleanup
datatwo1 <- data %>%
  select(iso, name, year, p, pop, t, s, a)
  #select(-shr, -avg, -apop)
  
  
#---------- 52. Method2: Compiled border changes 
# compiled border shares
indextwo2 <- indextwo %>%
  filter(count!=1) %>%
  select(-kind) %>%
  group_by(iso, year, p) %>%
  summarise(
    corr_ts = sum_or_na(corr_ts),   
    corr_bs  = sum_or_na(corr_bs),
    corr_t   = sum_or_na(corr_t),
    k_ts   = sum_or_na(k_ts),
    k_bs   = sum_or_na(k_bs),
    .groups = "drop"
  )  %>%
  add_count(iso, year, name = "count")

concernedtwo<- indextwo2 %>%
  select(iso,year) %>%
  distinct(iso, year, .keep_all = TRUE)
  
indextwo2t <- indextwo2 %>%
  select(iso, year, p, corr_t)%>%
  mutate(
    p = case_when(
      p == 25000 ~ 50000,
      p == 10000 ~ 25000,
      p == 0     ~ 10000,
      TRUE       ~ p               # keep original if no match
    )
  ) %>%
  filter(!is.na(corr_t))


indextwo2s <- indextwo2 %>%
  select(-corr_t) %>%
  filter(!if_all(c(corr_bs, corr_ts), is.na)) %>%
  add_count(iso, year, name = "count") %>%
  arrange(iso, year, p) %>%
  group_by(iso, year) %>%
  mutate(
    corr_bs2 = corr_bs, 
    corr_ts2 = corr_ts,
    corr_bs = corr_bs - if_else(!is.na( lag(corr_bs2)), lag(corr_bs2),  0),
    corr_ts = corr_ts - if_else(!is.na(lead(corr_ts2)), lead(corr_ts2), 0)
  ) %>%
  ungroup() %>%
  select(-corr_bs2, -corr_ts2)

# Keep only concerned iso year
data <- data_org %>%
 inner_join(concernedtwo, by = c("iso", "year"))
  

# Correct the thresholds
data <- data %>%
  left_join(indextwo2t, by = c("iso", "year","p"))%>%
  mutate(
    t=if_else(!is.na(corr_t),corr_t, t),
    changed_t=if_else(!is.na(corr_t),1, 0)
  )%>%
  select(-corr_t)

# Correct the target shares
data <- data %>%
  left_join(indextwo2s, by = c("iso", "year","p"))%>%
  mutate(
    s = if_else(!is.na(corr_bs), corr_bs, s),
    s = if_else(!is.na(corr_ts), corr_ts, s),
    changed_s= if_else(!is.na(corr_bs) | !is.na(corr_ts), 1, 0) ) %>%
  select(-corr_ts, -corr_bs, -k_ts, -k_bs, -count)

# re-balance the rest of the series
data <- data %>%
  group_by(iso, year) %>%
  mutate(
    s_rest = sum(s[changed_s == 0], na.rm = TRUE),
    s_modified = sum(s[changed_s == 1], na.rm = TRUE)
  ) %>%
  ungroup() %>%
  mutate(
    adj_ratio=  s_rest/(1-s_modified),
    s = if_else(changed_s == 0, s / adj_ratio, s),
    t = if_else(changed_s == 0 & changed_t == 0, t / adj_ratio, t),
    )

# Correct the rest of the distribution
data <- data %>%
  mutate(
    a_adj= a/((s*1)/pop),
  ) %>%
  mutate(
    a= (s*1)/pop,
    t=if_else(changed_s == 1 & changed_t==0, t/a_adj, t)) %>%
  select(-s_modified, -changed_t, -changed_s, -s_rest, -adj_ratio, -a_adj)

# 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()

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

#  Cleanup
datatwo2 <- data %>%
  select(iso, name, year, p, pop, t, s, a)
  #select(-shr, -avg, -apop)

#---------- 5. Method3: changes  exclusively to the thresholds ----------------------------------------------------------------------------
## Some corrections just imply thresholds targets. while the correction can be done directly, as in the precedent threshold correction, 
## this can result in a dis-balanced distribution. For the requested correction, the only-thresholds corrections happen only in the middle (50)
## and in one of the middle border ( 90 & 25). The poposed method is: 1 to correct the threhsolds directly. 2 to correct the border share by 
## the shame ratio that corrected the "middle border". 3. to estimated the diasjustment of the whol series and correct the other shares by this ratio.

# Generate index for files with threshold correction
indexthree <- index %>%
  filter(kinds_available=="t") %>% # Keep corrections that are just thresholds
  select (iso, year , p , k_t) %>%
  rename(p_t=p) %>%
  pivot_wider(
    names_from = p_t,       # the column whose unique values become column names
    values_from = k_t,      # the column that fills the cells
    names_prefix = "k_t_"   # optional: prefix for new columns
  )

## Correct thresholds
data <- data_org %>%
  inner_join(indexthree, by = c("iso", "year"), relationship = "many-to-many") %>%
  mutate(
    modif_t =0,
    modif_s =0
  )

#for (p_val in p_t_available) {
#  colname <- paste0("k_t_", p_val)
#  data <- data %>%
#    mutate(
#      t = if_else(p == p_val & !is.na(.data[[colname]]),
#                      .data[[colname]]*t, t),
#      modif_t = if_else(p == p_val & !is.na(.data[[colname]]),
#                  1, modif_t )
#      )
#}

# if assertion fails, the next code block may need to be  changed
expected <- c(10000, 50000, 90000)

if(length(p_t_available) == 3 && all(sort(p_t_available) == sort(expected))) {
  print("Check passed: vector p_t_available contains exactly the expected elements")
} else {
  stop("Check failed: Thesholds to be corrected in p_t_available are not only 50 and 90 or10, so the code bellow may need to be changed")
}

## correct the data
data <- data %>%
  mutate(
    t = case_when(
      #correction of border sides
      p <= 10000 & !is.na(k_t_10000) ~ t * k_t_10000,
      p >= 90000 & !is.na(k_t_90000) ~ t * k_t_90000,
      
      # correction of the mid sides
      (p == 25000 | p==50000) & !is.na(k_t_10000) ~ t * k_t_50000,
      (p == 75000 | p==50000)  & !is.na(k_t_90000) ~ t * k_t_50000,
      TRUE ~ t
    ),
    s = case_when(
      # correction of the border sided
      p <= 10000 & !is.na(k_t_10000) ~ s * k_t_10000,
      p >= 90000 & !is.na(k_t_90000) ~ s * k_t_90000,
      # correction of the mid sides
      (p == 25000 ) & !is.na(k_t_10000) ~ s * k_t_50000,
      (p == 75000 | p==50000)  & !is.na(k_t_90000) ~ s * k_t_50000,
      TRUE ~ s
    ),
    modif_t = case_when(
      # correction of the mid sides
      (p<=50000) & !is.na(k_t_10000) ~ 1,
      (p>=50000)  & !is.na(k_t_90000) ~ 1,
      TRUE ~ modif_s
    ),
    modif_s = case_when(
      # correction of the mid sides
      (p<50000) & !is.na(k_t_10000) ~ 1,
      (p>=50000)  & !is.na(k_t_90000) ~ 1,
      TRUE ~ modif_s
    )
  )


# Correct the rest of the series
data<- data %>%
  group_by(iso, year) %>%
  mutate(
    observed = sum(ifelse(modif_s == 0, s, 0), na.rm = TRUE),
    expected = sum(ifelse(modif_s == 1, s, 0), na.rm = TRUE),
  ) %>%
  ungroup() %>%
  mutate (
    expected=1-expected,
    adj=expected/observed,
    s=if_else(modif_s==0, s*adj,s),
    t=if_else(modif_t==0, t*adj,t)
  ) %>%
  select(iso, name, year, p, pop, t, s, a)

# readjust the average in fucniton of the new shares
data <- data %>%
  mutate(
    a= (s*1)/pop
  )

# 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()

# 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
datathree <- data %>%
  select(-shr, -avg, -apop)


#------------ . Append calcualted datasets
data_modif <- bind_rows(dataone, datatwo1, datatwo2, datathree)

rm (index, indextwo, indextwo2, input_data, data, concernedtwo, current_wid)

#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
#            2. Call the rest of the historical Series.
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------

# -------- 1.  Call current data non-modified ----------------------------------
aggregates_hist <- download_wid (
  indicators = c("ahweal","anninc"),
  year       = 1820:1970,
  perc       = "p0p100",
  ages       = 992,
  pop        = "i")

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

#--------- 1. Import current data (already gpintereised) -----------------------
# Keep iso-year-widcode to be changed
current_wid <- download_wid (
  indicators = c("sptinc","aptinc","tptinc"),
  year       = 1820:1970,
  ages       = 992, 
  pop = "j") %>%
  rename(iso=country) %>%
  filter(!(iso == "US" & year >= 1920))

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

# Extract  currenct data of the WID for the section 2
input_data <- current_wid %>%
  filter(percentile %in% c(
    "p0p10", "p0p25", "p0p50", "p0p75", "p0p90", "p0p99",
    "p10p100", "p25p100", "p50p100", "p75p100", "p90p100", "p99p100"
  )) %>%
  mutate(
    variable = case_when(
      str_detect(percentile, "p100") & variable!="t" ~ paste0("t", variable),
      str_detect(percentile, "^p0")  & variable!="t" ~ paste0("b", variable),
      TRUE                           ~ variable
    )
  ) %>%
  pivot_wider(
    names_from = variable,
    values_from = value
  ) %>%
  rename(p=percentile) %>%
  select(iso, year, p, t, ts, bs) %>%
  rename_with(~ paste0("curr_", .x), c(t, ts, bs)) %>%
  mutate(
    p = case_when(
      p == "p0p10"   ~ 10000,
      p == "p0p25"   ~ 15000,
      p == "p0p50"   ~ 50000,
      p == "p0p75"   ~ 75000,
      p == "p0p90"   ~ 90000,
      p == "p0p99"   ~ 99000,
      p == "p10p100" ~ 10000,
      p == "p25p100" ~ 25000,
      p == "p50p100" ~ 50000,
      p == "p75p100" ~ 75000,
      p == "p90p100" ~ 90000,
      p == "p99p100" ~ 99000,
      TRUE ~ NA_real_
    )
  )



current_wid <-current_wid %>%
  rename(
    p       = percentile
  )

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) 


  
# exlcude modifed iso-years
iso_year_to_exclude <- data_modif %>% 
  distinct(iso, year)

current_wid <- current_wid %>%
  anti_join(iso_year_to_exclude, by = c("iso", "year"))


# -------- 2.  Prepare data for gpinter ----------------------------------------
# Create fractile groups
breaks <- c(0,10000,25000,50000,75000, 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 %>%
  rename(p = ftile) %>%
  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_nomodif <- data_org %>%
  select(-shr, -avg, -apop) %>%
  mutate(
    name = "historical_sptinc992j"
  )

rm(iso_year_to_exclude, data)

#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
#            3. 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)
)

data <- bind_rows(data_modif, data_nomodif) %>%
  arrange(iso, year, p)


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,"/output2_round2-gpinterized.dta")))