library(tidyverse) # for data wrangling
library(texreg) # for charts
library(lmtest) # Likelihood ratio test
library(broom) # model tables
library(gtools) # significance
library(Zelig) # simulation
library(viridis) # colors

# Save the diamonds dataset as an object in my environment
jp <- read_csv("jp_solar.csv") %>%
  select(muni_code, date, year, solar_rate, disaster, pop) %>%
  arrange(muni_code, date)

# View first 3 rows of dataset
jp %>% head(3)

# Let's count up in how many rows the outcome is zero, vs. not zero
jp %>%
  group_by(values = if_else(solar_rate == 0, "0 Units", "1 Unit or More")) %>%
  summarize(count = n()) %>%
  ungroup() %>%
  mutate(percent = count / sum(count) * 100,
         percent = paste(round(percent,1), "%", sep = "")) %>%
  # Visualize!
  ggplot(mapping = aes(x = values, y = count, 
                       fill = values, 
                       label = percent)) +
  geom_col(color = "black") +
  geom_text(nudge_y = -200, size = 7, color = "white") +
  scale_fill_viridis(option = "mako", discrete = TRUE, begin = 0.8, end = 0.2) +
  guides(fill = "none") +
  theme_classic(base_size = 14) +
  labs(x = "Municipality-Month Observations by Level of Rooftop Solar Adoption",
       y = "Observations (n = 6,321)")

# Let's a vector of my unique observations, sorted from smallest to largest.
myoutcome <- jp$solar_rate %>% unique() %>% sort()
# Let's view the first 3
myoutcome %>% head(3)

# Get your constant
myconstant <- myoutcome[2] / 2
# View your constnat
myconstant

# Get a dataset comparing our different transformation options
compare <- bind_rows(
  # Let's first get a version of the dataset as it is originally
  jp %>%
    mutate(type = "Original\nNo Adjustment",
           finding = "Excludes 38.4%",
           position = 0),
  # Next, let's get a version where we impute zeros with our realistic constnat
  jp %>%
    mutate(solar_rate = if_else(
      condition = solar_rate == 0, 
      true = myconstant, false = solar_rate),
      type = "Realistic\nAdjustment",
      finding = "Shows frequent\nlow values\npretty well",
      position = 0),
  # Then, let's get a version where we impute zeros with our unrealistic constant
  jp %>%
    mutate(solar_rate = if_else(
      condition = solar_rate == 0, 
      true = 0.000000000001, false = solar_rate),
           type = "Unrealistic\nAdjustment",
           finding = "Creates\nvery strong\nartificial outliers",
           position = -10)
)

# Let's map their distributions
compare %>%
  ggplot(mapping = aes(x = log(solar_rate))) +
  geom_histogram(mapping = aes(y = ..density..), bins = 30, color = "white") +
  geom_density(fill = "steelblue", alpha = 0.5) +
  geom_text(mapping = aes(x = position, y = 0.50, 
                          label = finding), size = 5) +
  facet_wrap(~type, scales = "free_x") +
  theme_classic(base_size = 14) +
  theme(panel.border = element_rect(fill = NA, color = "black")) +
  labs(x = "Logged Solar Rate", y = "(%) Frequency (n = 6321)")

# Let's a vector of my unique observations, sorted from smallest to largest.
myoutcome <- jp$solar_rate %>% unique() %>% sort()

# Get your constant
myconstant <- myoutcome[2] / 2

# Update our data.frame with new value
jp2 <- jp %>%
    mutate(solar_log = if_else(
      condition = solar_rate == 0, 
      true = myconstant,
      false = solar_rate) %>%
        log())

jp3 <- jp2 %>%
  # Must order by city and time
  arrange(muni_code, date) %>%
  # Must group by municipality - IMPORTANT
  group_by(muni_code) %>%
  # Then lag the solar rate, using n to specify how many steps.
  mutate(
    # First let's lag by 1 row,
    solar_log_lag1 = lag(solar_log, n = 1),
    # Then let's lag by 2 rows
    solar_log_lag2 = lag(solar_log, n = 2),
    # Then let's lag by 3 rows
    solar_log_lag3 = lag(solar_log, n = 3)) %>%
  ungroup() %>%
  # Turn date into a factor for easy modeling
  mutate(date = factor(date))

jp3 %>%
  select(muni_code, date, solar_log, contains("lag") ) %>%
  head(5)

m1 <- jp3 %>%
  lm(formula = solar_log ~ disaster + date)

m2 <- jp3 %>%
  lm(formula = solar_log ~ disaster + date +
       solar_log_lag1)

m3 <- jp3 %>%
  lm(formula = solar_log ~ disaster + date +
       solar_log_lag1 + solar_log_lag2)

m4 <- jp3 %>%
  lm(formula = solar_log ~ disaster + date +
       solar_log_lag1 + solar_log_lag2 + solar_log_lag3)

library(texreg)
htmlreg(
  list(m1,m2,m3,m4), 
  # We'll omit our many date-fixed effects for readability
  omit.coef = "date",
  caption = "*** p < 0.001; ** p < 0.01; * p < 0.05.<br>With monthly fixed effects, omitted from table.")

z1 <- jp3 %>%
  drop_na() %>%
  lm(formula = solar_log ~ disaster + date)

z2 <- jp3 %>%
  drop_na() %>%
  lm(formula = solar_log ~ disaster + date +
       solar_log_lag1)

z3 <- jp3 %>%
  drop_na() %>%
  lm(formula = solar_log ~ disaster + date +
       solar_log_lag1 + solar_log_lag2)

z4 <- jp3 %>%
  drop_na() %>%
  lm(formula = solar_log ~ disaster + date +
       solar_log_lag1 + solar_log_lag2 + solar_log_lag3)

htmlreg(
  list(z1,z2,z3,z4), 
  # We'll omit our many date-fixed effects for readability
  omit.coef = "date",
  caption = "*** p < 0.001; ** p < 0.01; * p < 0.05.<br>With monthly fixed effects, omitted from table.")

library(lmtest)
library(broom)
library(gtools)

result <- lrtest(z1,z2,z3,z4) %>%
  tidy() %>%
  mutate(sig = stars.pval(p.value))

result

library(Zelig)

# Let's get the interquartile range for logged solar adoption
myrange <- jp3$solar_log %>% 
  quantile(probs = c(0.25, 0.75), na.rm = TRUE) 

# Simulate effects for high solar for 3-months prior
simlag3 <- z4 %>%
  # Convert to Zelig object
  to_zelig() %>%
  # Set range of predictor values
  setx(solar_log_lag3 = seq(from = myrange[1], to = myrange[2], length.out = 10),
       solar_log_lag2 = seq(from = myrange[1], to = myrange[2], length.out = 10),
       solar_log_lag1 = seq(from = myrange[1], to = myrange[2], length.out = 10)) %>%
  # Simulate expected value
  sim() %>%
  # Extract expected value
  zelig_qi_to_df() %>%
  select(x = solar_log_lag3, ev = expected_value)

# Simulate effects for high solar for 2-months prior
simlag2 <- z4 %>%
  # Convert to Zelig object
  to_zelig() %>%
  # Set range of predictor values
  setx(solar_log_lag2 = seq(from = myrange[1], to = myrange[2], length.out = 10),
       solar_log_lag1 = seq(from = myrange[1], to = myrange[2], length.out = 10)) %>%
  # Simulate expected value
  sim() %>%
  # Extract expected value
  zelig_qi_to_df() %>%
  select(x = solar_log_lag2, ev = expected_value)

# Simulate effects for high solar for 1 months prior
simlag1 <- z4 %>%
  # Convert to Zelig object
  to_zelig() %>%
  # Set range of predictor values
  setx(solar_log_lag1 = seq(from = myrange[1], to = myrange[2], length.out = 10)) %>%
  # Simulate expected value
  sim() %>%
  # Extract expected value
  zelig_qi_to_df() %>%
  select(x = solar_log_lag1, ev = expected_value)

mysim <- bind_rows(
  simlag1 %>% 
    mutate(type = "1 month"),
  simlag2 %>% 
    mutate(type = "2 months"),
  simlag3 %>% 
    mutate(type = "3 months")) %>%
  # exponentiate the predictor and expected value
  # To bring them back into their original range
  mutate(ev = exp(ev),
         x = exp(x))

mystats <- mysim %>%
  # Get a 95% confidence interval for each type and predictor
  group_by(type, x) %>%
  summarize(lower = quantile(ev, probs = 0.025),
            median = quantile(ev, probs = 0.50),
            upper = quantile(ev, probs = 0.975))

mystats %>% 
  ggplot(mapping = aes(x = x, y = median, ymin = lower, ymax = upper, fill = type)) +
  geom_ribbon(alpha = 0.5, color = "white") +
  scale_fill_viridis(option = "plasma", begin = 0.2, end = 0.8, discrete = TRUE) +
  labs(fill = "Increased\nSolar Adoption\nfor X months\nprior",
       x = "Rate of Solar Adoption in Prior Month(s)",
       y = "Expected Rate of Solar Adoption\n(with 95% simulated confidence interval)",
       subtitle = "Effect of Past Solar Adoption\non Future Solar Adoption",
       caption = "Based on 1000 simulations in the Zelig package in R,\nfor an otherwise average Japanese city.") +
  theme_bw()