# Load required libraries
library(fixest)     # For fixed effects regression
library(haven)      # For loading Stata .dta files
library(dplyr)      # For data manipulation

# Load data
url <- "https://github.com/scunning1975/mixtape/raw/master/Giffen.dta"
data <- read_dta(url)

# Filter data to Hunan province and the minimum person_id (mimicking Stata's "if province=="Hunan" & person_id==min_p")
data_hunan <- data %>%
  filter(province == "Hunan", person_id == min_p)

# Column 1: Full sample
model_full <- feols(pct_ch_hh_rice ~ pct_ch_sub_rice_arc + pct_ch_hh_pay + pct_ch_hh_nonwage + pct_ch_hh_people | county_time, 
                    data = data_hunan, cluster = ~hhid)

# Column 2: Calorie share <= 0.8
model_cal_share_leq_08 <- feols(pct_ch_hh_rice ~ pct_ch_sub_rice_arc + pct_ch_hh_pay + pct_ch_hh_nonwage + pct_ch_hh_people | county_time, 
                                data = filter(data_hunan, hh_staple_calorie_share_1 <= 0.8), cluster = ~hhid)

# Column 3: Calorie share > 0.8
model_cal_share_gt_08 <- feols(pct_ch_hh_rice ~ pct_ch_sub_rice_arc + pct_ch_hh_pay + pct_ch_hh_nonwage + pct_ch_hh_people | county_time, 
                               data = filter(data_hunan, hh_staple_calorie_share_1 > 0.8), cluster = ~hhid)

# Column 4: Calorie share between 0.6 and 0.8
model_cal_share_06_08 <- feols(pct_ch_hh_rice ~ pct_ch_sub_rice_arc + pct_ch_hh_pay + pct_ch_hh_nonwage + pct_ch_hh_people | county_time, 
                               data = filter(data_hunan, hh_staple_calorie_share_1 >= 0.6 & hh_staple_calorie_share_1 <= 0.8), 
                               cluster = ~hhid)

# Displaying results
etable(model_full, model_cal_share_leq_08, model_cal_share_gt_08, model_cal_share_06_08,
       dict = c(pct_ch_hh_rice = "Percent Change in Rice Consumption",
                pct_ch_sub_rice_arc = "Percent Change in Subsidy on Rice",
                pct_ch_hh_pay = "Percent Change in Household Pay",
                pct_ch_hh_nonwage = "Percent Change in Non-Wage Income",
                pct_ch_hh_people = "Percent Change in Household Size"),
       headers = c("Full Sample", "Calorie Share <= 0.8", "Calorie Share > 0.8", "Calorie Share 0.6-0.8"))