---
title: "Exercise template for 'Earnings management'"
author: Your name
format: html
bibliography: book.bib
---

```{r}
#| include: false
#| message: false
library(tidyverse)
library(DBI)
library(dbplyr)
library(farr)
library(broom)
library(furrr)
library(knitr)
```

## Measuring earnings management

### Discussion questions

1. @Jones:1991vx focuses on a small number of firms.
Why does @Jones:1991vx have such a small sample?
What are the disadvantages of a small sample?
Are there advantages of a smaller sample or narrower focus?

2. What are the primary conclusions of @Jones:1991vx?
Which table presents the main results of @Jones:1991vx?
Describe the empirical test used in that table.
Can you suggest an alternative approach?
What do you see as the primary challenges to the conclusions of @Jones:1991vx?

3. Can you think of refinements to the broad research question?
What tests might you use to examine these?

4. @McNichols:1988vq state at the outset that their paper "examines whether managers manipulate earnings."
Is this a good statement of the main research question of @McNichols:1988vq?
If not, suggest an alternative summary of the research questions of @McNichols:1988vq.

5. What do @McNichols:1988vq mean by "nondiscretionary accruals"?
How "operationalizable" is this concept?^[See <https://en.wiktionary.org/wiki/operationalizable#English>.]

6. @McNichols:1988vq say "if $\mathit{DA}$ were observable, accrual-based tests of earnings management would be expressed in terms of the following regression:

$$ \mathit{DA} = \alpha + \beta \textit{PART} + \epsilon $$
where $\textit{PART}$ is a dummy variable that partitions the data into two groups for which earnings management predictions are specified".
@Healy:1985jg points out that bonus plans can give managers incentives to increase earnings or decrease earnings depending on the situation.
How is this problematic for the formulation of @McNichols:1988vq above?
How might a researcher address this?

7. What are the benefits and costs of focusing on a single item (bad debt expense) in a study of earnings management?

8. The main results of @McNichols:1988vq are in Tables 6 and 7.
How persuasive do you find the evidence of earnings management found in the "residual provision" columns of those tables?

9. How well does the $\mathit{PART}$ framework apply to @Jones:1991vx?
Does the framework require modification for this paper?
In which periods would  $\mathit{PART}$ be set to one in @Jones:1991vx?

## Evaluating measures of earnings management

```{r}
#| include: false
#| label: model-functions
calc_accruals <- function(df) {
  df |>
    group_by(gvkey) |>
    arrange(datadate) |>
    mutate(lag_at = lag(at),
           d_ca = act - lag(act),
           d_cash = che - lag(che),
           d_cl = lct - lag(lct),
           d_std = dlc - lag(dlc),
           d_rev = sale - lag(sale),
           d_rec = rect - lag(rect)) |>
    ungroup() |>
    mutate(acc_raw =  (d_ca - d_cash - d_cl + d_std) - dp)
}

fit_jones <- function(df) {
  fm <- lm(acc_at ~ one_at + d_rev_at + ppe_at - 1,
           data = df, model = FALSE, subset = !part)

  df |>
    mutate(nda_jones = predict(fm, newdata = df),
           da_jones = acc_at - nda_jones) |>
    select(fyear, nda_jones, da_jones)
}

fit_mod_jones <- function(df) {
  fm <- lm(acc_at ~ one_at + d_rev_alt_at + ppe_at - 1,
           data = df, model = FALSE, subset = !part)
  df |>
    mutate(nda_mod_jones = predict(fm, newdata = df),
           da_mod_jones = acc_at - nda_mod_jones) |>
    select(fyear, nda_mod_jones, da_mod_jones)
}

get_nda <- function(df) {

  df_mod <-
    df |>
    calc_accruals() |>
    mutate(sic2 = str_pad(as.character(sic %/% 100), width = 2, pad = "0"),
           acc_at = acc_raw / lag_at,
           one_at = 1 / lag_at,
           d_rev_at = d_rev / lag_at,
           d_rev_alt_at = (d_rev - d_rec) / lag_at,
           ppe_at = ppegt / lag_at) |>
    group_by(sic2, fyear) |>
    mutate(acc_ind = median(if_else(part, NA, acc_at), na.rm = TRUE)) |>
    ungroup()

  da_healy <-
    df_mod |>
    group_by(gvkey) |>
    arrange(fyear) |>
    mutate(nda_healy = mean(if_else(part, NA, acc_at), na.rm = TRUE),
           da_healy = acc_at - nda_healy,
           nda_deangelo = lag(acc_at),
           da_deangelo = acc_at - nda_deangelo) |>
    ungroup() |>
    select(gvkey, fyear, part, nda_healy, da_healy, nda_deangelo,
           da_deangelo)

  df_jones <-
    df_mod |>
    nest_by(gvkey) |>
    reframe(fit_jones(data))

  df_mod_jones <-
    df_mod |>
    nest_by(gvkey) |>
    reframe(fit_mod_jones(data))

  fit_industry <- function(df) {
    fm <- lm(acc_at ~ acc_ind, data = df, model = FALSE, subset = !part)

    df |>
      mutate(nda_industry = suppressWarnings(predict(fm, newdata = df)),
             da_industry = acc_at - nda_industry) |>
      select(fyear, nda_industry, da_industry)
  }

  df_industry <-
    df_mod |>
    nest_by(gvkey) |>
    reframe(fit_industry(data))

  da_healy |>
    left_join(df_jones, by = c("gvkey", "fyear")) |>
    left_join(df_mod_jones, by = c("gvkey", "fyear")) |>
    left_join(df_industry, by = c("gvkey", "fyear"))
}

fit_model <- function(df, measure = "healy") {
  df |>
    nest_by(gvkey) |>
    summarize(model = list(lm(as.formula(str_c("da_", measure, " ~ part")),
                              model = FALSE, data = data)),
              .groups = "drop") |>
    mutate(measure = !!measure)
}

multi_fit <- function(df) {
  models <- c("healy", "deangelo", "jones", "mod_jones", "industry")
  models |>
    map(\(x) fit_model(df, x)) |>
    list_rbind()
}

get_stats <- function(fm) {
  fm |>
    tidy() |>
    filter(term == "partTRUE") |>
    select(-term)
}

table_1_stats <- function(x) {
  tibble(mean = mean(x, na.rm = TRUE),
         sd = sd(x, na.rm = TRUE),
         q1 = quantile(x, p = 0.25, na.rm = TRUE),
         median = median(x, na.rm = TRUE),
         q3 = quantile(x, p = 0.75, na.rm = TRUE))
}

h_test <- function(fm) {
  coefs <- coef(summary(fm))

  if ("partTRUE" %in% rownames(coefs)) {
    t_stat <- coefs["partTRUE", "t value"]
    df <- fm$df.residual

    return(tibble(neg_p01 = pt(t_stat, df, lower = TRUE) < 0.01,
                  neg_p05 = pt(t_stat, df, lower = TRUE) < 0.05,
                  pos_p01 = pt(t_stat, df, lower = FALSE) < 0.01,
                  pos_p05 = pt(t_stat, df, lower = FALSE) < 0.05))
  } else {
    return(tibble(neg_p01 = NA, neg_p05 = NA, pos_p01 = NA, pos_p05 = NA))
  }
}

binom_test <- function(x, p) {
  x <- x[!is.na(x)]
  binom.test(sum(x), length(x), p = p)$p.value
}
```

```{r}
#| label: database-connection
#| include: false
#| cache: false
db <- dbConnect(RPostgres::Postgres(),
                bigint = "integer", check_interrupts = TRUE)
funda <- tbl(db, Id(schema = "comp", table = "funda"))
company <- tbl(db, Id(schema = "comp", table = "company"))
```

```{r}
#| include: false
#| label: shared-data
#| cache: true
#| cache.vars: [acc_data_raw, test_sample, merged_sample_1, earn_deciles, merged_sample_2]
#| cache.extra: !expr digest::digest(list(formals(calc_accruals), body(calc_accruals)))
sics <-
  company |>
  select(gvkey, sic) |>
  mutate(sic = as.integer(sic))

funda_mod <-
  funda |>
  filter(indfmt == "INDL", datafmt == "STD",
         consol == "C", popsrc == "D") |>
  left_join(sics, by = "gvkey") |>
  mutate(sic = coalesce(sich, sic))

acc_data_raw <-
  funda_mod |>
  filter(!is.na(at),
         pddur == 12,
         !between(sic, 6000, 6999)) |>
  mutate(across(c(che, dlc, sale, rect), \(x) coalesce(x, 0))) |>
  select(gvkey, datadate, fyear, at, ib, dp, rect, ppegt, ni, sale,
         act, che, lct, dlc, sic) |>
  filter(between(fyear, 1950, 1991)) |>
  arrange(gvkey, fyear) |>
  collect()

test_sample <-
  acc_data_raw |>
  calc_accruals() |>
  filter(lag_at > 0, sale > 0, ppegt > 0, !is.na(acc_raw),
         !is.na(d_rev), !is.na(d_rec), !is.na(ppegt)) |>
  group_by(gvkey) |>
  filter(n() >= 13) |>
  ungroup() |>
  arrange(gvkey, fyear) |>
  select(gvkey, fyear)

set.seed(2022)
sample_1_firm_years <-
  test_sample |>
  group_by(gvkey) |>
  arrange(fyear, .by_group = TRUE) |>
  filter(row_number() > 2) |>
  ungroup() |>
  mutate(rand = rnorm(n())) |>
  arrange(rand) |>
  distinct(gvkey, .keep_all = TRUE) |>
  slice_head(n = 1000) |>
  select(gvkey, fyear) |>
  mutate(part = TRUE)

sample_1 <-
  test_sample |>
  semi_join(sample_1_firm_years, by = "gvkey") |>
  left_join(sample_1_firm_years, by = c("gvkey", "fyear")) |>
  mutate(part = coalesce(part, FALSE))

merged_sample_1 <-
  sample_1 |>
  inner_join(acc_data_raw, by = c("gvkey", "fyear"))

earn_deciles <-
  acc_data_raw |>
  calc_accruals() |>
  mutate(earn = ib / lag_at) |>
  semi_join(test_sample, by = c("gvkey", "fyear")) |>
  arrange(gvkey, fyear) |>
  mutate(earn_dec = ntile(earn, 10)) |>
  select(gvkey, fyear, earn_dec)

set.seed(2023)
sample_2_firm_years <-
  test_sample |>
  group_by(gvkey) |>
  arrange(fyear, .by_group = TRUE) |>
  mutate(event_ok = row_number() > 2) |>
  ungroup() |>
  inner_join(earn_deciles, by = c("gvkey", "fyear")) |>
  filter(earn_dec == 10, event_ok) |>
  mutate(rand = rnorm(n())) |>
  arrange(rand) |>
  distinct(gvkey, .keep_all = TRUE) |>
  slice_head(n = 1000) |>
  select(gvkey, fyear) |>
  mutate(part = TRUE)

sample_2 <-
  test_sample |>
  semi_join(sample_2_firm_years, by = "gvkey") |>
  left_join(sample_2_firm_years, by = c("gvkey", "fyear")) |>
  mutate(part = coalesce(part, FALSE))

merged_sample_2 <-
  sample_2 |>
  inner_join(acc_data_raw, by = c("gvkey", "fyear"))
```

```{r}
#| include: false
#| cache: false
dbDisconnect(db)
```

```{r}
#| include: false
#| label: baseline-models
#| cache: true
#| cache.vars: [results]
#| cache.extra: !expr digest::digest(lapply(list(calc_accruals, fit_jones, fit_mod_jones, get_nda, fit_model, multi_fit), function(f) list(formals(f), body(f))))
reg_data <- get_nda(merged_sample_1)
results <- multi_fit(reg_data)
```

```{r}
#| echo: false
head(results)
```

```{r}
#| label: tbl-em-1
#| tbl-cap: "Results of tests of earnings management: Sample 1"
#| echo: false
results |>
  mutate(stats = map(model, get_stats)) |>
  unnest_wider(stats) |>
  pivot_longer(estimate:statistic, names_to = "stat") |>
  group_by(measure, stat) |>
  summarize(table_1_stats(value), .groups = "drop") |>
  kable(digits = 4)
```

```{r}
#| include: false
test_results <- results |> mutate(map_dfr(model, h_test))
```

```{r}
#| label: tbl-em-2
#| tbl-cap: "Type I error rates"
#| echo: false
test_results |>
  group_by(measure) |>
  summarize(tests = sum(!is.na(neg_p05)),
            across(matches("p0"), \(x) mean(x, na.rm = TRUE))) |>
  kable(digits = 4)
```

```{r}
#| label: tbl-em-2-sig
#| tbl-cap: "p-values for null that Type I error rates equal size of tests"
#| echo: false
test_results |>
  group_by(measure) |>
  summarize(neg_p01 = binom_test(neg_p01, p = 0.01),
            neg_p05 = binom_test(neg_p05, p = 0.05),
            pos_p01 = binom_test(pos_p01, p = 0.01),
            pos_p05 = binom_test(pos_p05, p = 0.05)) |>
  kable(digits = 4)
```

The tables count available tests; sparse industry-year benchmarks can leave some Industry Model event predictions unavailable.

```{r}
#| include: false
one_gvkey <- merged_sample_1 |>
  count(gvkey) |>
  filter(n == 20) |>
  arrange(gvkey) |>
  slice_head(n = 1) |>
  pull(gvkey)
```

```{r}
df_test <- merged_sample_1 |> filter(gvkey == one_gvkey)
```

```{r}
df_mod <-
  df_test |>
  calc_accruals() |>
  mutate(acc_at = acc_raw / lag_at,
         one_at = 1 / lag_at,
         d_rev_at = d_rev / lag_at,
         d_rev_alt_at = (d_rev - d_rec) / lag_at,
         ppe_at = ppegt / lag_at) |>
  ungroup()
```

```{r}
fm1a <- lm(acc_at ~ one_at + d_rev_at + ppe_at,
          data = df_mod, subset = !part)
```

```{r}
res1 <-
  df_mod |>
  mutate(nda_jones = predict(fm1a, newdata = pick(everything()))) |>
  select(fyear, part, acc_at, nda_jones) |>
  mutate(da_jones = acc_at - nda_jones)
```

```{r}
fm2a <- lm(da_jones ~ part, data = res1)
```

```{r}
fm2 <- lm(acc_at ~ one_at + d_rev_at + ppe_at + part,
          data = df_mod)
```

### Discussion questions and exercises

1. What interpretation do @Dechow:1995wr provide for their Table 1 results?

2. Compare the results in @tbl-em-1 with those in Table 1 of @Dechow:1995wr.
What differences appear to be significant?

3. Compare the values in the standard deviation column of Table 1 of @Dechow:1995wr with other statistics.
Do these differences make sense?
Or do they suggest anomalies in the underlying data?

4. Compare the values in the standard deviation column of the "earnings management" rows of Table 1 of @Dechow:1995wr with the values in the mean column of the standard error rows.
What is the relationship between these values?
What would you expect the relationship between these values to be?
Do you observe similar relations in @tbl-em-1?

5. Focusing on the Healy Model, DeAngelo Model, and the Industry Model, compare the rejection rates in @tbl-em-2 with those presented in Table 2 of @Dechow:1995wr.
What might explain any differences?
Could these be attributed to differences between our results in @tbl-em-1 and those reported in Table 1 of @Dechow:1995wr?
Or do you expect that these differences have another cause?

6. How do you interpret the results from `binom_test()` reported in @tbl-em-2-sig?
Does it make sense to interpret each of the columns independent of the others?

7. Confirm that the coefficient on $\textit{PART}$ from the regression in `fm2a` can be recovered from the regression in `fm2`.
  Are the standard errors the same?

8. Modify the code above to check that the same holds for the Modified Jones Model.

9. We described the Jones Model above as "a (differently) modified Jones Model".
In what way is the model different from the Jones Model estimated in `fit_jones()` above?
Does the @Salkever:1976ue equivalence hold if we use the Jones Model from `fit_jones()`?
If so, why?
If not, how might this affect how you would use the Jones Model and the @Salkever:1976ue approach?
(For example, do we expect the "(differently) modified Jones Model" to produce materially different results from the Jones Model?)

10. Do the issues related to a first and second stage apply to either the Healy Model or the DeAngelo Model or both?
If so, could we apply the @Salkever:1976ue approach to address these issues?
If not, are there "one-stage" equivalents to the Healy Model and DeAngelo Model approaches as implemented above?

11. Produce an equivalent of Table 3 from @Dechow:1995wr by adapting the code used above to create `merged_sample_2` and @tbl-em-2.
(*Challenge version*: Implement the approach of @Salkever:1976ue in doing so.)

12. Produce an equivalent of Table 4 from @Dechow:1995wr by adapting the code used above to create `merged_sample_2` and @tbl-em-2.

## Power of tests of earnings management

```{r}
#| include: false
#| label: power-functions
manipulate <- function(df, level = 0, type) {
  df <-
    df |>
    group_by(gvkey) |>
    arrange(datadate) |>
    mutate(ni_ratio = median(if_else(part, NA, ni / sale), na.rm = TRUE),
           lag_at = lag(at),
           manip_amt = lag_at * level,
           margin_ok = is.finite(ni_ratio) & ni_ratio > 0,
           manip_amt_gross = if_else(margin_ok, manip_amt / ni_ratio, 0))

  if (type == "expense") {
    df |>
      mutate(lct = if_else(part, lct - manip_amt, lct)) |>
      ungroup()
  } else if (type == "revenue") {
    df |>
      mutate(sale = case_when(part ~ sale + manip_amt,
                              lag(part) ~ sale - lag(manip_amt),
                              .default = sale),
             rect = if_else(part, rect + manip_amt, rect),
             act = if_else(part, act + manip_amt, act)) |>
      ungroup()
  } else if (type == "margin") {
    df |>
      mutate(sale = case_when(part & margin_ok ~
                                sale + manip_amt_gross,
                              lag(part) & margin_ok ~
                                sale - lag(manip_amt_gross),
                              .default = sale),
             rect = if_else(part & margin_ok,
                            rect + manip_amt_gross, rect),
             act = if_else(part & margin_ok,
                           act + manip_amt_gross, act),
             lct = if_else(part & margin_ok,
                           lct + manip_amt_gross - manip_amt, lct)) |>
      ungroup()
  } else {
    df |>
      ungroup()
  }
}

h_test_5 <- function(fm) {
  coefs <- coef(summary(fm))

  if ("partTRUE" %in% rownames(coefs)) {
    t_stat <- coefs["partTRUE", "t value"]
    df <- fm$df.residual
    return(pt(t_stat, df, lower = FALSE) < 0.05)
  } else {
    return(NA)
  }
}

calc_power <- function(df, level, type) {
  df |>
    manipulate(level = level, type = type) |>
    get_nda() |>
    multi_fit() |>
    mutate(rej_null = map_lgl(model, h_test_5)) |>
    group_by(measure) |>
    summarize(tests = sum(!is.na(rej_null)),
              prop_reject = mean(rej_null, na.rm = TRUE),
              .groups = "drop")
}
```

```{r}
#| include: false
#| label: power-results
#| cache: true
#| cache.vars: [power_plot_data]
#| cache.extra: !expr digest::digest(lapply(list(manipulate, get_nda, fit_jones, fit_mod_jones, multi_fit, h_test_5), function(f) list(formals(f), body(f))))
plan(multisession, workers = 2)
manip_df <-
  expand_grid(level = seq(from = 0, to = 1, by = 0.1),
              manip_type = c("expense", "revenue", "margin")) |>
  mutate(results = future_map2(
    level, manip_type,
    \(x, y) calc_power(merged_sample_1, x, y),
    .options = furrr_options(packages = c("dplyr", "tidyr", "purrr",
                                          "broom", "stringr"))
  ))
plan(sequential)

power_plot_data <- manip_df |> unnest(results)
```

```{r}
#| label: fig-power
#| fig-cap: Power functions for tests of earnings management
#| fig-height: 5
#| fig-pos: H
#| echo: false
power_plot_data |>
  ggplot(aes(x = level, y = prop_reject)) +
  geom_line() +
  facet_grid(measure ~ manip_type)
```

### Discussion questions

1. How do the results in @fig-power compare with those in Figure 4 of @Dechow:1995wr?

2. According to the [SEC's filing](https://www.sec.gov/litigation/admin/3439329.txt) referenced above related to B&L, "B&L recognized, in contravention of GAAP and the Company's own revenue recognition policies, \$42.1 million of revenue, resulting in at least a \$17.6 million, or 11%, overstatement of the net income originally reported for its 1993 fiscal year."
According to [a subsequent SEC filing,](https://www.sec.gov/Archives/edgar/data/10427/0000010427-95-000006.txt) B&L's total assets for 1994 were \$2,457,731,000 (it seems reasonable to assume that the 1993 value was not radically different from this).
Based on this information (plus any information in the SEC's filing), which of @Dechow:1995wr's three categories did B&L's earnings management fall into?
What is the approximate magnitude relative to the $x$-axes of the plots in Figure 4 of @Dechow:1995wr (or the equivalent above)?
Based on these data points, what is the approximate estimated probability of the various models detecting earnings management of this magnitude?

3. What do you view as the implications of the power analysis conducted above for research on earnings management?
Are these implications consistent with the extensive literature on earnings management subsequent to @Dechow:1995wr?
If so, explain why.
If not, how would you reconcile the inconsistencies?

4. Does each of the three forms of earnings management implemented in `manipulate()` above agree precisely with the corresponding description in @Dechow:1995wr [pp. 201--202]?
If not, does one approach seem more correct than the other?
(Note that *one* issue arises with negative or zero net income ratio.
How are such cases handled by @Dechow:1995wr and by `manipulate()`?)