## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5,
  eval = requireNamespace("did", quietly = TRUE)
)

## ----libs---------------------------------------------------------------------
library(phdid)
library(did)
data(mpdta, package = "did")

## ----first-stage--------------------------------------------------------------
first <- att_gt(
  yname = "lemp", tname = "year", idname = "countyreal",
  gname = "first.treat", control_group = "notyettreated",
  data = mpdta, bstrap = FALSE, cband = FALSE
)

d <- ph_data(first)
d

## ----homogeneity--------------------------------------------------------------
homogeneity_test(d)

## ----l0-----------------------------------------------------------------------
fit_l0 <- l0_ph(d)
fit_l0

## ----l0-plot------------------------------------------------------------------
plot(fit_l0)

## ----l0-agg-------------------------------------------------------------------
rbind(
  flexible = aggregate(flex_twfe(d), "overall"),
  l0 = aggregate(fit_l0, "overall")
)

## ----bayes--------------------------------------------------------------------
fit <- bayes_ph(d, alpha = 1, iters = 6000, burn = 1000, seed = 7,
                progress = FALSE)
fit

## ----bayes-agg----------------------------------------------------------------
aggregate(fit, "overall")

## ----coclust------------------------------------------------------------------
round(coclustering(fit), 2)
plot(fit)

## ----enumerate----------------------------------------------------------------
exact <- enumerate_partitions(d, alpha = 1)
exact

data.frame(
  quantity = c("E[# groups]", "overall effect", "interval lower",
               "interval upper"),
  gibbs = c(fit$m_mean, aggregate(fit, "overall")$estimate,
            aggregate(fit, "overall")$conf.low,
            aggregate(fit, "overall")$conf.high),
  exact = c(exact$m_mean, exact$overall$mean, exact$overall$lower,
            exact$overall$upper)
)

cat("largest co-clustering discrepancy:",
    round(max(abs(fit$coclust - exact$coclust)), 4), "\n")

## ----cov-check----------------------------------------------------------------
cmp <- covariance_check(d, alpha = 1, iters = 1500, burn = 400, seed = 1)
cmp$summary

## ----alpha, fig.height = 7----------------------------------------------------
sens <- alpha_sensitivity(
  d, alpha_grid = c(0.1, 0.5, 1, 5, 14, 100),
  iters = 1500, burn = 400, seed = 3
)
sens
plot_sensitivity(sens)

## ----alpha-cells, fig.height = 5----------------------------------------------
sens_cells <- alpha_sensitivity(
  d, alpha_grid = c(0.1, 0.5, 1, 5, 14, 100),
  type = "cells", iters = 1500, burn = 400, seed = 5
)
plot_sensitivity(sens_cells)

## ----lambda-cells, fig.height = 5---------------------------------------------
lam <- lambda_sensitivity(d, type = "cells", by = "m")
plot_sensitivity(lam)

## ----dynamic------------------------------------------------------------------
aggregate(fit, "dynamic")