23  Beyond OLS

This part of the book explores some additional topics that are too advanced to be covered in Part I and that do not quite fit in Part III. These include generalized linear models (this chapter), extreme values and sensitivity analysis (Chapter 24), matching (Chapter 25), and a chapter on prediction that also introduces machine learning approaches (Chapter 26). Apart from Chapter 24 building on this chapter, the material in each of these chapters is largely independent of that in the others in this section. Material in this chapter and in Chapters 24 and 25 does build on ideas in Part III to some extent, but likely could be covered before Part III. The material in Chapter 26 is fairly discrete and could be tackled by students with the level of understanding provided by the first few chapters of Part I.

In Chapter 3, we studied the fundamentals of regression analysis with a focus on OLS regression. While we have used OLS as the foundation for almost all regression analyses so far, OLS will at times appear not to be appropriate. For example, if the dependent variable is binary, some aspects of OLS may be difficult to interpret and OLS may be econometrically inefficient. Examples of binary dependent variables include an indicator variable for accounting fraud (see Chapter 26) or for using a Big Four auditor (see Chapter 25).

Similar issues apply to inherently non-negative dependent variables, such as total assets or count variables, such as the number of management forecasts, the number of 8-K filings, or the number of press releases, as seen in Guay et al. (2016), the focus of this chapter.

In this chapter, we cover generalized linear models (GLMs), a class of models that includes probit and logit models—which assume a binary dependent variable—and Poisson regression—which is most naturally motivated by count variables.

A second theme of this chapter is the introduction of additional data skills relevant for applied researchers. In previous chapters we have eschewed larger data sets beyond those available through WRDS. However, most researchers need to know how to manage larger data sets, so we introduce some skills for doing so here.

In particular, we collect metadata for every filing made on the SEC EDGAR system since 1993, producing a data frame that occupies 3 GB if loaded into R. While 3 GB certainly is not “big data”, it is large enough that more efficient ways of working have a significant impact on the facility of manipulating and analysing the data. In the analysis below, we process the data into parquet files that we can access using the DuckDB database package.

Tip

The code in this chapter uses the packages listed below. For instructions on how to set up your computer to use the code found in this book, see Section 1.2. Quarto templates for the exercises below are available at https://github.com/iangow/far_templates/blob/main/README.md.

library(tidyverse)
library(DBI)
library(farr)
library(modelsummary)
library(arrow)         # write_parquet(), read_parquet()
library(googledrive)   # drive_download(), drive_deauth()
library(fixest)        # feols()
library(lmtest)        # coeftest()
library(sandwich)      # vcovHC()
library(httr2)         # request(), req_*(), resp_body_html()
library(rvest)         # html_elements(), html_table()

23.1 Complexity and voluntary disclosure

Like Li et al. (2018) (covered in Chapter 21), Guay et al. (2016) examine the phenomenon of voluntary disclosure. A core idea of Guay et al. (2016) is that users of a more complex firm’s disclosures require more information to understand the business sufficiently well to invest in it. Guay et al. (2016, p. 235) posit that “economic theory suggests that managers will use other disclosure channels to improve the information environment” when their reporting is more complex. A number of “other disclosure channels” are explored in Guay et al. (2016), including the issuance of management forecasts (the primary measure), filing of 8-Ks, and firm-initiated press releases.

Guay et al. (2016) find evidence of greater voluntary disclosure by firms with more complex financial statements. For example, Table 3 of Guay et al. (2016) presents results of OLS regressions with two different measures of financial statement complexity, and with and without firm fixed effects.

23.1.1 Discussion questions

1. Guay et al. (2016, p. 252) argue that the “collective results from our two quasi-natural experiments … validate that our text-based measures of financial statement complexity reflect, at least in part, the complexity of the underlying accounting rules.” Which specific table or figure of Guay et al. (2016) provides the most direct evidence in support of this claim? Can you suggest any alternative tests or presentations to support this claim?

2. What is the causal diagram implied by the regression specification in Table 3 of Guay et al. (2016)? Provide arguments for and against the inclusion of ROA, SpecialItems and Loss as controls.

3. Assuming that the causal effect of FS_complexity implied by the results in Table 3 exists, provide an explanation for the changes in the coefficients when firm fixed effects are added to the model. When might you expect these changes to have the opposite sign?

4. Guay et al. (2016, p. 234) conclude that “collectively, these findings suggest managers use voluntary disclosure to mitigate the negative effects of complex financial statements on the information environment.” Suppose you wanted to design a field experiment to test the propositions embedded in this sentence and were given the support of a regulator to do so (e.g., you can randomize firms to different regimes, as was done with Reg SHO discussed in Section 19.2). How would you implement this experiment? What empirical tests would you use? Are some elements of the hypothesis more difficult than others to test?

5. Clearly Guay et al. (2016) did not have the luxury of running a field experiment. Do you see any differences between the analyses you propose for your (hypothetical) field experiment and those in Guay et al. (2016)? What do you think best explains any differences?

23.2 Generalized linear models

In this section, we introduce the idea of generalized linear models (GLMs), which provides a powerful framework for analysing a wide class of settings beyond OLS. We begin by discussing maximum likelihood as an organizing principle for estimation before applying that to understand regression analysis with binary dependent variable.

After discussing the concept of marginal effects, we finish up this section with a discussion of GLM regressions with count variables.

23.2.1 Maximum likelihood

The concept of maximum likelihood estimation provides a useful framework for thinking about regression, including GLMs. As the name “ordinary least squares” suggests, OLS regression is often motivated as the linear estimator that minimizes the squared residuals in a setting in which \(y = X \beta\).

However, we can take another perspective. Suppose that \(y_i = X_i \beta + \epsilon_i\) where \(\epsilon_i \sim N(0, \sigma^2)\) and \(\epsilon_i\) is independent across individuals (\(i\)). As such, the density for \(y_i\) is given by

\[ f(y_i | X_i \beta) = \frac{1}{\sigma \sqrt{2 \pi}} \exp \left(-\frac{(y_i - X_i \beta)^2}{2\sigma^2} \right) \] But mathematically, the \(y_i\) and \(X_i\) can be seen as given and this function can be viewed instead as a function of \(\beta\) and relabelled as a likelihood function.

\[ \mathcal{L}(\beta; X_i, y_i) = \frac{1}{\sigma \sqrt{2 \pi}} \exp \left(-\frac{(y_i - X_i \beta)^2}{2\sigma^2} \right) \] Given independence, the density for the full sample \(y = (y_1, \dots, y_N)\) is simply the product of the densities of each observation and we can consider the likelihood function for the full sample in the same way.

\[ \mathcal{L}(\beta; X, y) = \prod_{i=1}^N \frac{1}{\sigma \sqrt{2 \pi}} \exp \left(-\frac{(y_i - X_i \beta)^2}{2\sigma^2} \right) \]

We can also write the log-likelihood as follows

\[ \begin{aligned} \log \mathcal{L}(\beta; X, y) &= \sum_{i=1}^N \left( \log \left(\frac{1}{\sigma \sqrt{2 \pi}} \right) + \left(-\frac{(y_i - X_i \beta)^2}{2\sigma^2} \right) \right) \\ &= \sum_{i=1}^N \log \left(\frac{1}{\sigma \sqrt{2 \pi}} \right) + \sum_{i=1}^N \left(-\frac{(y_i - X_i \beta)^2}{2\sigma^2} \right) \\ \end{aligned} \]

The idea of maximum likelihood estimation (MLE) is to choose the values \(\hat{\beta}\) and \(\hat{\sigma}\) that maximize the likelihood function. But note that maximizing the log-likelihood will yield the same solution as maximizing the likelihood (because \(\log(x) > \log(y) \Leftrightarrow x > y\)) and that only the second term in the sum above is affected by the value of \(\beta\).

This means that maximizing the log-likelihood with respect to \(\beta\) is achieved by maximizing the second term. But the second term is nothing but the negative of the sum of square residuals. So the least-squares estimator is the maximum likelihood estimator in this setting.

So why introduce MLE in the first place? First, this equivalence will not hold in general. Merely assuming heteroskedasticity (i.e., that \(\sigma_i^2 \neq \sigma_j^2\) for \(i \neq j\)) is enough to create a wedge between least-squares and MLE.

Second, MLE can be shown to have desirable econometric properties in a broad class of settings. According to Wooldridge (2010, p. 385), “MLE has some desirable efficiency properties: it is generally the most efficient estimation procedure in the class of estimators that use information on the distribution of the endogenous variables given the exogenous variables.” Additionally, MLE has “a central role” (Wooldridge, 2010, p. 386) in estimation of models with limited dependent variables, which we cover in the next section.

23.2.2 Binary dependent variables

If we have a setting in which \(y_i \in \{0, 1\}\) (i.e., the dependent variable is binary), then a generic framework has the conditional density of \(y_i\) as

\[ g(y_i | X_i; \beta ) = \begin{cases} h(X_i \cdot \beta), & y_i = 1 \\ 1 - h(X_i \cdot \beta), &y_i = 0 \\ 0, &\text{otherwise} \\ \end{cases} \]

where \(X = (1, x_{i1}, x_{i2}, \dots, x_{ik})\), \(\beta = (\beta_0, \beta_1, \dots, \beta_k)\), and \(h(\cdot)\) is a known function \(\mathbb{R} \to [0, 1]\).

With this setup, we define the logit model in terms of \(h(\cdot)\):

\[ h(X\beta) = \frac{e^{X\beta}}{1 + e^{X\beta}} \] The probit model is defined using the standard normal distribution function \(\Phi(\cdot)\):

\[ h(X\beta) = \Phi(X\beta) \] If we assume that one of the above functions describes the underlying data-generating process, then we can construct the corresponding maximum-likelihood estimator of \(\beta\). For example, the probit model implies that

\[ f_{Y|X} = g(Y, X; \beta) = \Phi(X \beta)^{Y} (1 - \Phi(X \beta))^{1- Y} \] \[ \mathcal{L}(\beta | Y, X) = \prod_{i = 1}^n \Phi(X_i \beta)^{Y_i} (1 - \Phi(X_i \beta))^{1- Y_i} \] \[ \begin{aligned} \log \mathcal{L}(\beta | Y, X) &= \log \left( \prod_{i = 1}^n \Phi(X_i \beta)^{Y_i} (1 - \Phi(X_i \beta))^{1- Y_i} \right) \\ &= \sum_{i = 1}^n \log \left( \Phi(X_i \beta)^{Y_i} (1 - \Phi(X_i \beta))^{1- Y_i} \right) \\ &= \sum_{i = 1}^n {Y_i} \log \Phi(X_i \beta) + (1 - Y_i) \log (1 - \Phi(X_i \beta)) \\ \end{aligned} \] Thus, given \(Y\) and \(X\), we can obtain the ML estimate of \(\beta\) using a computer to find

\[ \hat{\beta}_{ML} = \underset{\beta \in \mathbb{R^{k+1}}}{\arg \max} \sum_{i = 1}^n {Y_i} \log \Phi(X_i \beta) + (1 - Y_i) \log (1 - \Phi(X_i \beta)) \]

A similar analysis can be conducted for logit regression.

But rarely would a researcher have good reasons to believe strongly that the correct model is the probit model, rather than the logit model, or vice versa. This raises questions about the consequences of estimation using MLE, but with the wrong model. In other words, what happens if we use probit regression when the correct model is logit? Or logit regression when the correct model is probit? And what happens if we simply use OLS to fit a model? (The term linear probability model is used to describe the fitting of OLS with a binary dependent variable, as the fitted values can be interpreted as probabilities that \(y=1\), even though these fitted values can fall outside the \([0, 1]\) interval.)

A natural way to examine these questions is with simulated data. We make a small function (get_probit_data()) that generates data using the assumptions underlying the probit model. We generate random \(x\) and then calculate \(h(X\beta) = \Phi(X\beta)\), which gives us a vector h representing the probability that \(y_i = 1\) for each \(i\). We then compare \(h(X_i \beta)\) with \(u_i\), a draw from the uniform distribution over \([0, 1]\).¹ If \(h(X_i \beta) > u_i\), then \(y = 1\); otherwise \(y = 0\). If \(h(X_i \beta) = 0.9\) then we have a 90% chance that \(h(X_i \beta) > u_i\). Finally, the get_probit_data() function returns the generated data.

get_probit_data <- function(beta_0 = -2, beta_1, n = 1000) {
  x <- rnorm(n)
  h <- pnorm(beta_0 + x * beta_1)
  y <- runif(n = n) < h
  tibble(x, y)
}

23.2.3 Marginal effects

An important concept with non-linear models such as those examined in this chapter is the marginal effect. Given a causal model \(\mathbb{E}[y|X] = f(X)\), the partial derivative of \(x_j\) at \(X = x\)—that is,

\[\left. \frac{\partial \mathbb{E}[y|X]}{\partial x_i} \right|_{X = x}\]

for \(j \in (1, \dots, k)\)—can be viewed as the marginal effect of a change in \(x_j\) on the value of \(\mathbb{E}[y|X]\) evaluated at \(X = x\).

If \(y\) is a binary variable, then \(\mathbb{E}[y_i|X_i] = \mathbb{P}[y_i = 1|X_i]\). If we have a probit model, then \(\mathbb{P}[y_i = 1|X_i] = \Phi(X_i \beta)\). Assuming that \(X_i\) is a continuous variable, we then have the following:²

\[ \begin{aligned} \frac{\partial \mathbb{E}[y_i|X_i]}{\partial x_{ij}} &= \frac{\partial \mathbb{P}[y_i = 1|X_i]}{\partial x_{ij}} \\ &= \frac{\partial \Phi(X_i \beta)}{\partial x_{ij}} \\ &= \frac{\partial \Phi(X_i \beta)}{\partial (X_i \beta)} \frac{\partial (X_i \beta)}{\partial x_{ij}} \\ &= \phi(X_i \beta) \beta_j \\ \end{aligned} \] That is, the marginal effect is the product of the probability density function and the coefficient \(\beta_i\). Note that we used the facts that the first derivative of the cumulative normal distribution \(\Phi(\cdot)\) is the normal density function \(\phi(\cdot)\) and that the partial derivative of \(X_i \beta = \beta_0 + x_{i1} \beta_1 + \dots + x_{ij} \beta_j + \dots + x_{ik} \beta_k\) with respect to \(x_{ij}\) is \(\beta_j\).

With linear models such as those fitted using OLS where \(x_j\) only enters the model in a linear fashion, the marginal effect of \(x_j\) is simply \(\beta_j\) and does not depend on the value \(x_j\) at which it is evaluated. But in models in which \(x_j\) interacts with other variables or enters the model in a non-linear fashion, the marginal effect of \(x_j\) can vary. For example, in the following model, the marginal effect of \(x_1\) will depend on the values of \(x_1\) and \(x_2\) at which it is evaluated.³

\[ y_i = \beta_0 + x_{i1} \beta_1 + x_{i2} \beta_2 + x_{i1} x_{i2} \beta_3 + x^2_{i2} \beta_4 + \epsilon_i \]

With non-linear models the marginal effect of \(x_j\) will generally depend on the value \(x\) at which it is evaluated. One common approach evaluates the marginal effect at the sample mean of \(\overline{X} = (\overline{x}_1, \dots, \overline{x}_k)\). Another approach would be to report the mean of the marginal effects evaluated for each \(X\) in the sample, but this approach may be difficult to interpret if the signs of the marginal effects vary across the sample distribution.

\[ \begin{aligned} \mathbb{E} \left[\frac{\partial \mathbb{E}[y_i|X_i]}{\partial x_{ij}} \right] &= \mathbb{E} \left[\phi(X_i \beta) \right] \beta_j \\ \end{aligned} \]

The following function implements the latter computation for probit, logit, and OLS models. While a more sophisticated version of this function would infer the model’s type from its attributes, this function is adequate for our needs here. The mfx package offers estimates of marginal effects for a variety of models.

get_mfx <- function(fm, type) {
  pdf <- 
    case_when(type == "probit" ~ mean(dnorm(predict(fm, type = "link"))),
              type == "logit" ~ mean(dlogis(predict(fm, type = "link"))),
              type == "Poisson" ~ mean(predict(fm, type = "response")),
              type == "OLS" ~ 1)
  mfx <- pdf * coef(fm)
  mfx
}

The following fit_model() convenience function fits a model, extracts statistics from that model, then returns a data frame containing those statistics.⁴

fit_model <- function(data = df, type = "OLS") {
  family <- ifelse(type %in% c("probit", "logit"), binomial,
                    ifelse(type == "Poisson", poisson, gaussian))
  
  link <- case_when(type %in% c("probit", "logit") ~ type, 
                    type == "Poisson" ~ "log", 
                    .default = "identity")
  
  fm <- glm(y ~ x, family(link = link), data = data)
  coefs <- summary(fm)$coefficients
  
  tibble(type, 
         coef = coefs[2, 1],
         se = coefs[2, 2],
         t_stat = coefs[2, 3],
         p = coefs[2, 4],
         mfx = (get_mfx(fm, type))[["x"]])
}

We then create a function to run an iteration of the simulation: it gets the data using get_probit_data(), extracts statistics using fit_model() for three different model types, then returns a data frame with the results.

run_probit_sim <- function(i, beta_1) {
  df <- get_probit_data(beta_1 = beta_1)
  
  fms <- list(fit_model(df, type = "probit"),
              fit_model(df, type = "logit"),
              fit_model(df, type = "OLS"))
  
  bind_cols(tibble(i), 
            beta_1_true = beta_1, 
            list_rbind(fms))
}

Finally, we run the simulation n_iters times for each of two values of \(\beta_1\). We consider \(\beta_1 = 0\), so that the null hypothesis of no relation between \(x_1\) and \(y\) is true. We also consider \(\beta_1 = 0.3\), so that the null hypothesis of no relation between \(x_1\) and \(y\) is false. The former allows us to check that the rejection rate matches the size of the test (\(\alpha = 0.05\) or \(\alpha = 0.01\)). The latter provides insights into the power of the tests.

set.seed(2024)
n_iters <- 1000L

results <- 
  list(
    map(1:n_iters, \(x) run_probit_sim(x, beta_1 = 0.0)),
    map(1:n_iters, \(x) run_probit_sim(x, beta_1 = 0.3))) |>
  bind_rows() |>
  system_time()

   user  system elapsed 
 15.644   0.752  16.558 

To understand the relationships between the various approaches, we focus on the case where \(\beta_1 \neq 0\) (so there is an underlying relationship between \(x\) and \(y\) to estimate).

rel_df <- 
  results |> 
  filter(beta_1_true != 0) |> 
  pivot_wider(names_from = type, values_from = coef:mfx)

We examine the relation between the estimated coefficients for the various models using the probit model as the benchmark (after all, we know that the underlying data-generating process is a probit model).

fms <- list(
  "Logit" = lm(coef_logit ~ coef_probit, data = rel_df),
  "OLS" = lm(coef_OLS ~ coef_probit, data = rel_df),
  "mfx_logit" = lm(mfx_logit ~ mfx_probit, data = rel_df),
  "mfx_OLS" = lm(coef_OLS ~ mfx_probit, data = rel_df))

The first thing to note from our simulation results is that there are strong relationships between the coefficients and marginal effects across the various models. In Table 23.1, we can see a high correlation (measured using the \(R^2\)) between probit, logit, and OLS coefficients and that this correlation is even higher when we move to estimates of the marginal effects. This strong relationship is also evident in Figure 23.1.

modelsummary(fms,
             estimate = "{estimate}{stars}",
             statistic = "statistic",
             gof_map = c("nobs", "r.squared"),
             stars = c('*' = .1, '**' = 0.05, '***' = .01))

Table 23.1: Relationships between model coefficients

	Logit	OLS	mfx_logit	mfx_OLS
(Intercept)	0.034***	0.001***	0.000	0.000
	(9.937)	(2.868)	(-0.191)	(-0.653)
coef_probit	2.172***	0.057***
	(202.626)	(54.862)
mfx_probit			0.997***	0.998***
			(309.056)	(235.526)
Num.Obs.	1000	1000	1000	1000
R2	0.976	0.751	0.990	0.982

rel_df |> 
  ggplot(aes(x = coef_probit, y = coef_logit)) +
  geom_point() + 
  geom_smooth(formula = "y ~ x", method = "lm", se = FALSE)

Figure 23.1: Relationship between probit and logit model coefficients

Of course, very little attention is paid to the magnitude of the coefficients from estimated probit or logit models in accounting research. While this partly reflects the focus on the statistical significance and sign of the coefficients in the NHST paradigm, it also reflects the difficulty of attaching economic meaning to the values of the coefficients. In this regard, more meaning can be placed on the marginal effect estimates.

Table 23.2 reveals that in our simulation analysis there is little difference in the estimated marginal effect sizes from the three approaches, even though only one (probit) is the correct model. The strong relationship between probit and logit fixed effects is also clear in Figure 23.2.

rel_df |> 
  ggplot(aes(x = mfx_probit, y = mfx_logit)) +
  geom_point() + 
  geom_smooth(formula = "y ~ x", method = "lm", se = FALSE)

Figure 23.2: Relationship between probit and logit marginal effect estimates

But even marginal effect sizes receive little attention in most accounting research papers, as the focus is generally on the coefficients’ signs and statistical significance. Table 23.2 provides insights on these issues too. When \(\beta_1 = 0\) (i.e., the null hypothesis is true), we see that the rejection rates are quite close to the size of the tests (either 5% or 1%) for OLS, logit, and probit and very close to each other. When \(\beta_1 = 0.3\) (i.e., the null hypothesis is false), we see that the rejection rates for OLS, logit, and probit are very similar. We have power over 90% at the 5% level and around 80% at the 1% level. As such, there is arguably little value in this setting in using probit or logit rather than OLS regression, which yields coefficients that are much easier to interpret.

results |> 
  group_by(type, beta_1_true) |>
  summarize(mean_coef = mean(coef),
            mean_mfx = mean(mfx),
            rej_rate_5 = mean(p < 0.05), 
            rej_rate_1 = mean(p < 0.01),
            .groups = "drop") |>
  arrange(beta_1_true, type)

Table 23.2: Probit simulation results

type	beta_1_true	mean_coef	mean_mfx	rej_rate_5	rej_rate_1
OLS	0.0	0.000	0.000	0.035	0.008
logit	0.0	0.006	0.000	0.033	0.008
probit	0.0	0.002	0.000	0.035	0.008
OLS	0.3	0.018	0.018	0.938	0.805
logit	0.3	0.696	0.018	0.938	0.800
probit	0.3	0.305	0.018	0.934	0.801

Of course, the analysis above should not be assumed to apply in all settings. But it seems reasonable for deviations from straightforward OLS analysis in specific settings to be justified. Fortunately, simulation analysis like that above is relatively straightforward to conduct.

23.2.4 Count variables

Another case in which OLS may be regarded as “inadequate” is when the dependent variable is a count variable. One member of the class of GLMs that models such a dependent variable is Poisson regression. In this model, we have \(y \sim \mathrm{Pois}(\mu)\) where \(\mu = \exp(\beta_0 + \beta_1 x)\), where \(\mathrm{Pois}(\mu)\) is the Poisson distribution with parameters \(\mu\).

As above, we start by making a small function (get_pois_data()) that generates data using the assumptions underlying the Poisson model.

get_pois_data <- function(beta_0 = -2, beta_1, n = 1000) {
  x <- rnorm(n)
  h <- exp(beta_0 + x * beta_1)
  y <- rpois(n = n, lambda = h)
  tibble(x, y)
}

We then create a function to run an iteration of the simulation: get the data using get_pois_data(), extract statistics using fit_model() for both Poisson and OLS, then return a data frame with the results.

run_pois_sim <- function(i, beta_1) {
  df <- get_pois_data(beta_1 = beta_1)
  
  fm_pois <- fit_model(df, type = "Poisson")
  fm_ols  <- fit_model(df, type = "OLS")
  
  bind_cols(tibble(i), 
            beta_1_true = beta_1, 
            bind_rows(fm_pois, fm_ols))
}

Finally, we run the simulation n_iters times for each of two values of \(\beta_1\). As above, we consider \(\beta_1 = 0\) and \(\beta_1 = 0.3\).

set.seed(2024)
n_iters <- 1000L

results_pois <- 
  list(
    map(1:n_iters, \(x) run_pois_sim(x, beta_1 = 0.0)),
    map(1:n_iters, \(x) run_pois_sim(x, beta_1 = 0.3))) |>
  bind_rows() |>
  system_time()

   user  system elapsed 
 10.256   0.242  10.624 

Table 23.3 reveals that in our simulation analysis there is little difference in the estimated marginal effect sizes between Poisson and OLS regression. When \(\beta_1 = 0\) (i.e., the null hypothesis is true), we see that the rejection rates are quite close to the size of the tests (either 5% or 1%) and to each other. When \(\beta_1 = 0.3\) (i.e., the null hypothesis is false), we see that the rejection rates for Poisson and OLS are very similar. Again, we have power over 90% at the 5% level and over 80% at the 1% level. As such, similar to what we saw in the analysis of probit analysis above, there is arguably little value in this setting in using Poisson regression rather than OLS regression. This seems especially true given that OLS yields coefficients that are much easier to interpret.

results_pois |> 
  group_by(type, beta_1_true) |>
  summarize(mean_coef = mean(coef),
            mean_mfx = mean(mfx),
            rej_rate_5 = mean(p < 0.05), 
            rej_rate_1 = mean(p < 0.01),
            .groups = "drop") |>
  arrange(beta_1_true, type)

Table 23.3: Poisson simulation results

type	beta_1_true	mean_coef	mean_mfx	rej_rate_5	rej_rate_1
OLS	0.0	0.000	0.000	0.059	0.011
Poisson	0.0	0.001	0.000	0.063	0.010
OLS	0.3	0.043	0.043	0.942	0.849
Poisson	0.3	0.302	0.043	0.940	0.850

23.3 Application: Complexity and voluntary disclosure

To enrich our understanding of GLMs and to examine some additional data skills, we return to a setting similar to that considered in Guay et al. (2016). Like Guay et al. (2016), we will examine the association between complexity of financial reporting and subsequent voluntary disclosures. Our purpose here is not to conduct a close replication of Guay et al. (2016). In place of the measures of financial reporting complexity and of voluntary disclosure used in Guay et al. (2016), we use measures that are readily available. Additionally, we do not restrict ourselves to the sample period used in Guay et al. (2016).

This is the first chapter in which we need to collect a substantial amount of data from sources other than WRDS. To facilitate this, we will perform a number of data steps as quite discrete steps. This allows the re-running of code for one data set without requiring us to rerun code for other data sets.

To make this approach easier, we create save_parquet(), which takes a (possibly remote) data frame and uses the write_parquet() function from the arrow library to save it as a parquet file. The parquet format is described in R for Data Science (Wickham et al., 2023, p. 393) as “an open standards-based format widely used by big data systems.” Parquet files provide a format optimized for data analysis, with a rich type system. More details on the parquet format can be found in Chapter 22 of R for Data Science. We use write_parquet() from the arrow package to store the data in sec_index_2022q2 in a file named sec_index_2022q2.parquet.

save_parquet <- function(df, name, schema = "", path = data_dir) {
  file_path <- file.path(path, schema, str_c(name, ".parquet"))
  write_parquet(collect(df), sink = file_path)
}

23.3.1 Measures of reporting complexity

Guay et al. (2016) consider a number of measures of reporting complexity, including measures based on the “readability”, measured using the Fog index (Guay et al., 2016, p. 239), and length of a firm’s 10-K, measured as the natural logarithm of the number of words. Loughran and McDonald (2014) find “that 10-K document file size provides a simple readability proxy that outperforms the Fog Index, does not require document parsing, facilitates replication, and is correlated with alternative readability constructs.” Thus, we use 10-K document file size here.

Data on 10-K document file size can be collected from SEC EDGAR, but gathering these directly from SEC EDGAR would involve costly web scraping. Thus, we simply get these data from the website provided by Loughran and McDonald (2014).

In this book so far, we have avoided substantial data collection efforts of the kind that are often central in conducting research. From a pedagogical perspective, such data collection efforts typically involve messy details and complications that could distract more than they enlighten. However, we feel that it is useful to readers to provide some examples of more substantial data-collection tasks.

We first specify a directory in which we will store the data we will collect and create this directory if it doesn’t already exist. While the use of the DATA_DIR environment variable does not add anything in the analysis here, it may have a role if we are using local file storage across projects (see Appendix E for more on this point).⁵

Tip

The following line of code below is only indicative. If you have already set DATA_DIR (e.g., if you have set up a data repository as discussed in Appendix E), then you many not need (or want) to run this line. If you are within an RStudio project when running the following line, then it should set DATA_DIR to a data folder inside that project, which is a reasonable choice for current purposes.

Sys.setenv(DATA_DIR = "data")

data_dir <- Sys.getenv("DATA_DIR")
if (!dir.exists(data_dir)) dir.create(data_dir)

Using an approach discussed in Appendix E, we will organize the data in this directory into subdirectories, which we can think of as schemas, much as data on WRDS are organized into schemas.

Specifically, we imagine that this code is being used for a hypothetical project named “glms” and we put all data specific to this project into a directory with that name. We store the full path to this directory, a subdirectory of data_dir, in project_dir.

project_dir <- file.path(data_dir, "glms")
if (!dir.exists(project_dir)) dir.create(project_dir)

It turns out that the Loughran and McDonald (2014) website provides a link to a CSV file stored in Google Drive and we can use the googledrive package to download this. This requires supplying the Google Drive ID (stored in id below) and a local file name to use (stored in lm_data below). The code below only downloads the file if it is not already present in the location stored in lm_data.⁶

drive_deauth()
lm_data <- file.path(project_dir, "lm_10x_summary.csv")
id <- "1puReWu4AMuV0jfWTrrf8IbzNNEU6kfpo"
if (!file.exists(lm_data)) drive_download(as_id(id), lm_data)

Having downloaded this CSV file, we can read it into R. While one way to read this file would use read_csv() from the readr package, we continue the approach of using database-backed dplyr used elsewhere in the book. However, rather than using a PostgreSQL database, we will use a local DuckDB database. DuckDB works much like PostgreSQL, but does not require a server to be set up and comes with excellent facilities for reading common file formats, such as CSV and parquet files (discussed below).

In fact, we can set up a temporary database with a single line of code:

db <- dbConnect(duckdb::duckdb())

Having created this database connection, we can use DuckDB’s native ability to read CSV files provided by its read_csv() function to create a database table as follows. We convert all variable names to lower case, convert cik to integer type and filing_date and cpr to date types, and handle the -99 value used to represent missing values in cpr. (According to the source website, cpr refers to the “conformed period of report” of the filing, which is the last date of the reporting period.)

Note that we force the table to be materialized as a remote temporary table by appending the command compute(). This takes less than a second, but means that the conversion from CSV and processing of dates and missing values and the like happens only once rather than each time we call on lm_10x_summary in later analyses. Judicious use of compute() can also enhance query performance even when a table is not used multiple times by making it easier for the query optimizer in our database to understand the query. Unfortunately, WRDS does not grant its users sufficient privileges to create temporary tables and thus we have not been able to use compute() in previous chapters.⁷

lm_10x_summary <- 
  tbl(db, str_c("read_csv('", lm_data, "')")) |> 
  rename_with(tolower) |>
  mutate(cik = as.integer(cik),
         filing_date = as.character(filing_date),
         cpr = if_else(cpr == -99, NA, as.character(cpr))) |>
  mutate(filing_date = as.Date(strptime(filing_date, '%Y%m%d')),
         cpr = as.Date(strptime(cpr, '%Y%m%d'))) |>
  compute()

This remote table lm_10x_summary is just like the remote tables we first saw in Section 6.1.

lm_10x_summary

# Source:   table<dbplyr_if0zhyOp7u> [?? x 25]
# Database: DuckDB v1.1.3 [iangow@Darwin 24.2.0:R 4.4.2/:memory:]
    cik filing_date acc_num cpr        form_type coname   sic ffind n_words
  <int> <date>      <chr>   <date>     <chr>     <chr>  <dbl> <dbl>   <dbl>
1 60512 1993-08-13  000006… 1993-06-30 10-Q      LOUIS…  1311    30    4068
2 66740 1993-08-13  000006… 1993-06-30 10-Q      MINNE…  2670    38    4389
3 60512 1993-10-07  000006… 1992-12-31 10-K-A    LOUIS…  1311    30    8719
4 60512 1993-11-10  000006… 1993-09-30 10-Q      LOUIS…  1311    30    4938
5 11860 1993-11-12  000001… 1993-09-30 10-Q      BETHL…  3312    19    3823
6 20762 1993-11-12  000095… 1993-09-30 10-Q      CLARK…  2911    30    4136
# ℹ more rows
# ℹ 16 more variables: n_unique_words <dbl>, n_negative <dbl>,
#   n_positive <dbl>, n_uncertainty <dbl>, n_litigious <dbl>,
#   n_strongmodal <dbl>, n_weakmodal <dbl>, n_constraining <dbl>,
#   n_negation <dbl>, grossfilesize <dbl>, netfilesize <dbl>,
#   nontextdoctypechars <dbl>, htmlchars <dbl>, xbrlchars <dbl>,
#   xmlchars <dbl>, n_exhibits <dbl>

From a quick inspection, we can see that the data begin in 1993, but only appear to be comprehensive from some time in 1995 or 1996.

lm_10x_summary |> 
  mutate(year = year(filing_date)) |>
  count(year) |>
  arrange(year) |>
  collect()

# A tibble: 29 × 2
   year     n
  <dbl> <dbl>
1  1993    13
2  1994  9752
3  1995 21880
4  1996 44531
5  1997 55484
6  1998 55735
# ℹ 23 more rows

Finally, we save the data to a parquet file in our glms project directory and disconnect from the in-memory DuckDB database.

lm_10x_summary |>
  save_parquet(name = "lm_10x_summary", schema = "glms")

dbDisconnect(db)

23.3.2 Data on EDGAR filings

As discussed above, Guay et al. (2016) consider a number of measures of voluntary disclosure. While measures based on management forecasts require data from IBES and measures of firm-initiated press releases require data from RavenPack, the measure based on the number of 8-Ks filed in a period can be collected from SEC EDGAR fairly easily and we focus on that measure for this reason.

Here we collect filings data from SEC EDGAR. Specifically, we collect index files for each quarter. These index files include details on each filing in the applicable quarter.

For example, if you visit https://www.sec.gov/Archives/edgar/full-index/2022/QTR2/, you can see a number of files available to download. As these files are mostly different ways of presenting the same information in a variety of formats, we will focus on the company.gz files.

The first practical challenge we face is that if we try to use download.file() to get company.gz, we will get a “403 Forbidden” error message because the SEC does not allow anonymous downloads. To avoid this, we set the HTTPUserAgent option using our email address (obviously, use your actual email address here).

options(HTTPUserAgent = "your_name@some_email.com")

The file company.gz is a compressed fixed-width format text file of the kind examined in Chapter 9.

url <- str_c("https://www.sec.gov/Archives/edgar/full-index/",
             "2022/QTR2/company.gz")
t <- tempfile(fileext = ".gz")
download.file(url, t)

sec_index_2022q2 <-
  read_fwf(t, 
           fwf_cols(company_name = c(1, 62),
                    form_type = c(63, 77),
                    cik = c(78, 89),
                    date_filed = c(90, 101),
                    file_name = c(102, NA)),
           col_types = "cciDc",
           skip = 10,
           locale = locale(encoding = "macintosh")) 

sec_index_2022q2 |> select(company_name, cik)

# A tibble: 291,592 × 2
  company_name                                           cik
  <chr>                                                <int>
1 006 - Series of IPOSharks Venture Master Fund, LLC 1894348
2 008 - Series of IPOSharks Venture Master Fund, LLC 1907507
3 01 Advisors 03, L.P.                               1921173
4 01Fintech LP                                       1928541
5 07 Carbon Capture, LP                              1921663
6 09 Carbon Capture, LP                              1936048
# ℹ 291,586 more rows

sec_index_2022q2 |> select(form_type, date_filed, file_name)

# A tibble: 291,592 × 3
  form_type date_filed file_name                                  
  <chr>     <date>     <chr>                                      
1 D         2022-04-08 edgar/data/1894348/0001892986-22-000003.txt
2 D         2022-06-07 edgar/data/1907507/0001897878-22-000003.txt
3 D         2022-05-20 edgar/data/1921173/0001921173-22-000001.txt
4 D         2022-05-13 edgar/data/1928541/0001928541-22-000001.txt
5 D         2022-04-05 edgar/data/1921663/0001921663-22-000001.txt
6 D         2022-06-30 edgar/data/1936048/0001936048-22-000001.txt
# ℹ 291,586 more rows

As the file we downloaded above is over 4 MB in size and we will need to download dozens of such files, we likely want to store these data in some permanent fashion. While one approach would be to use a PostgreSQL database to which we have write access, we do not assume here that you have such a database.⁸

Another approach would be to simply download all the company.gz files and process them each time we use them. The problem with this approach is that this processing is time-consuming and thus best done once.

Here we use the parquet file format to store the data. As it is easy to imagine using these on projects other than our hypothetical "glms" project, we will store the SEC EDGAR filings data discussed in edgar, a subdirectory of DATA_DIR defined above named.

edgar_dir <- file.path(data_dir, "edgar")
if (!dir.exists(edgar_dir)) dir.create(edgar_dir)
sec_index_2022q2 |> 
  save_parquet("sec_index_2022q2", schema = "edgar")

This parquet file can be read into R in many ways. For example, we could use the open_dataset() function from the arrow library. However, we will use a DuckDB database connection db—created using a single line of code below—and the DuckDB-specific load_parquet() function provided by the farr package to create a database table.⁹

db <- dbConnect(duckdb::duckdb())

sec_index <- load_parquet(db, "sec_index_2022q2", schema = "edgar")

If we were using just a single index file, the benefits of using DuckDB rather than the open_dataset() function might be less clear. However, we will soon see how easy it is to work with parquet files using DuckDB and also look to combine these index files with the data in the lm_10x_summary data frame we created above and data from Compustat and CRSP.

Our next step is to embed the code we used above to download and process company.gz for a single quarter in a function that can download an SEC index file and save it as a parquet file. Armed with this function, we can then process multiple index files in one step. Naturally, the function we create below (get_sec_index()) addresses possibilities ignored above, such as download failures.

get_sec_index <- function(year, quarter, overwrite = FALSE) {

  pq_path <- str_c(edgar_dir, "/sec_index_",
                    year, "q", quarter, ".parquet")
  if (file.exists(pq_path) & !overwrite) return(TRUE)
  
  # Download the zipped index file from the SEC website
  url <- str_c("https://www.sec.gov/Archives/edgar/full-index/",
               year,"/QTR", quarter, "/company.gz")
  
  t <- tempfile(fileext = ".gz")
  result <- try(download.file(url, t))

  # If we didn't encounter an error downloading the file, parse it
  # and save as a parquet file
  if (!inherits(result, "try-error")) {
    temp <-
      read_fwf(t, fwf_cols(company_name = c(1, 62),
                           form_type = c(63, 77),
                           cik = c(78, 89),
                           date_filed = c(90, 101),
                           file_name = c(102, NA)),
               col_types = "ccicc", skip = 10,
               locale = locale(encoding = "macintosh")) |>
      mutate(date_filed = as.Date(date_filed))

    write_parquet(temp, sink = pq_path)
    return(TRUE)
  } else {
    return(FALSE)
  }
}

We next create a data frame containing the years and quarters for which we want to collect SEC index files. SEC EDGAR extends back to 1993 and index files for the current quarter will appear after the first day of the quarter. We can use the crossing() function from the tidyr package to create the Cartesian product of the quarters and years.

now <- now(tz = 'America/New_York') - days(1)
current_year <- as.integer(year(now))
current_qtr <- quarter(now)
year <- 1993L:current_year
quarter <- 1:4L

index_files_to_get <-
  crossing(year, quarter) |>
  filter(year < current_year |
           (year == current_year & quarter <= current_qtr)) 

Finally, we can download the SEC index files using the following few lines of code.

Tip

The following code takes several minutes to run. This code took us around 4 minutes to run, but the time required will depend on things such as the speed of your internet connection.

index_files_downloaded <-
  index_files_to_get |>
  mutate(available = map2(year, quarter, get_sec_index))

After running the code above, we have 125 index files saved as parquet files in edgar_dir occupying 468 MB of drive space. This represents much more than 468 MB of data because parquet files are highly compressed.

Reading these data into our database is simple and quick. We can use a wildcard to combine several parquet files ("sec_index*.parquet" here) into a single table.

sec_index <- load_parquet(db, "sec_index*", "edgar")

We can check that we have data for 125 quarters. As we saw above with the data on 10-Ks, there are few observations in 1993 and a significant increase from 1994 to 1995.

sec_index |>
  mutate(year = year(date_filed),
         quarter = quarter(date_filed)) |>
  count(year, quarter) |>
  arrange(year, quarter) |>
  collect()

# A tibble: 129 × 3
   year quarter     n
  <dbl>   <dbl> <dbl>
1  1993       1     4
2  1993       2     4
3  1993       3     7
4  1993       4    20
5  1994       1 20879
6  1994       2 16500
# ℹ 123 more rows

Each filing on SEC EDGAR has a form_type. As can be seen in Table 23.4, the most common form type is Form 4, which is filed whenever there is a material change in the holdings of company insiders. The second most common form type is Form 8-K, which the SEC describes as “the ‘current report’ companies must file with the SEC to announce major events that shareholders should know about.”

sec_index |> 
  count(form_type, sort = TRUE) |>
  collect(n = 10)

Table 23.4: Most common form types on SEC EDGAR

form_type	n
4	9,052,950
8-K	1,974,529
SC 13G/A	867,478
3	814,247
424B2	797,075
10-Q	689,052
6-K	531,027
497	515,442
SC 13G	456,345
D	420,508

As mentioned above, we will construct a measure of disclosure based on the number of 8-K filings. While 8-K filings represent the second-most common form of filing, they are still a small share of all filings. Figure 23.3 shows some clear trends in the number of filings over time. Form 8-K filings appear to have peaked in number in 2005 and to have been in slow decline in number since then.

sec_index |>
  mutate(year = year(date_filed),
         form_type = if_else(form_type %in% c("4", "8-K"), 
                             form_type, "Other")) |>
  filter(year < max(year, na.rm = TRUE)) |>
  count(year, form_type) |>
  mutate(last_year =  year == max(year, na.rm = TRUE),
         label = if_else(last_year, form_type, NA)) |>
  ggplot(aes(x = year, y = n / 1000, 
             colour = form_type, linetype = form_type)) +
  geom_line() +
  geom_label(aes(label = label), na.rm = TRUE, 
             position = position_nudge(y = -0.4)) + 
  scale_y_continuous(name = "Number of filings (thousands)") +
  scale_x_continuous(breaks = seq(1993, current_year, 2)) + 
  theme(axis.text.x = element_text(angle = 90),
        legend.position = "none")