---
title: "Some examples of use-cases using the market models of the package"
output:
  rmarkdown::html_vignette:
    toc: false
vignette: >
  %\VignetteIndexEntry{Some examples of use-cases using the market models of the package}
  %\VignetteEngine{knitr::rmarkdown}
---

This short tutorial covers the very basic use cases to get you started with *markets*. More usage details can be found in the documentation of the package.

```{r, include = FALSE}
if (requireNamespace("knitr", quietly = TRUE)) {
  knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
}
```

## Setup the environment

Load the required libraries.

```{r setup.libraries}
library(markets)
library(Formula)
```

Prepare the data. Normally this step is long and depends on the nature of the data and the considered market. For this example, we will use simulated data. Although we could simulate data independently from the package, we will use the top-level simulation functionality of *markets* to simplify the process. See the documentation of `simulate_data` for more information on the simulation functionality. Here, we simulate data using a data generating process for a market in disequilibrium with stochastic price dynamics.

```{r setup.data}
nobs <- 1000
tobs <- 5

alpha_d <- -1.3
beta_d0 <- 24.8
beta_d <- c(2.3, -1.02)
eta_d <- c(2.6, -1.1)

alpha_s <- 0.6
beta_s0 <- 16.1
beta_s <- c(2.9)
eta_s <- c(-1.5, 3.2)

gamma <- 1.2
beta_p0 <- 0.9
beta_p <- c(-0.1)

sigma_d <- 1
sigma_s <- 1
sigma_p <- 1
rho_ds <- 0.0
rho_dp <- 0.0
rho_sp <- 0.0

seed <- 4430

stochastic_adjustment_data <- simulate_data(
  "diseq_stochastic_adjustment", nobs, tobs,
  alpha_d, beta_d0, beta_d, eta_d,
  alpha_s, beta_s0, beta_s, eta_s,
  gamma, beta_p0, beta_p,
  sigma_d = sigma_d, sigma_s = sigma_s, sigma_p = sigma_p,
  rho_ds = rho_ds, rho_dp = rho_dp, rho_sp = rho_sp,
  seed = seed
)
```

## Estimate the models


Prepare the basic parameters for model initialization. The `simulate_data` call uses `Q` for the simulated traded quantity, `P` for the simulated prices, `id` for subject identification, and `date` for time identification. It automatically creates the demand-specific variables `Xd1` and `Xd2`, the supply-specific variable `Xs1`, the common (i.e., both demand and supply) variables `X1` and `X2`, and the price dynamics' variable `Xp1`. 

```{r model.parameters}
market_spec <- Q | P | id | date ~ P + Xd1 + Xd2 + X1 + X2 | P + Xs1 + X1 + X2
```

The market specification has to be modified in two cases. For the `diseq_directional`, the price variable is removed from the supply equation because the separation rule of the model can only be used for markets with exclusively either inelastic demand or supply. For the `diseq_stochastic_adjustment`, the right-hand side of the price dynamics equation is appended in the market specification.

By default, the models are estimated by allowing the demand, supply, and price equations to have correlated error shocks. The default verbosity behavior is to display errors and warnings that might occur when estimating the models. 

By default, all models are estimated using full information maximum likelihood based on the `"BFGS"` optimization algorithm. The first `equilibrium_model` call modifies the estimation behavior and estimates the model using two stage least squares. The `diseq_basic` call modifies the default optimization behavior and estimates the model using the `"Nelder-Mead"` optimization methods. 

Standard errors are by default assumed to be homoscedastic. The second `equilibrium_model` and `diseq_deterministic_adjustment` calls modify this behavior by calculating clustered standard errors based on the subject identifier, while the `diseq_basic` and `diseq_directional` calls modify it by calculating heteroscedastic standard errors via the sandwich estimator.


```{r model.constructor}
eq_reg <- equilibrium_model(
  market_spec, stochastic_adjustment_data,
  estimation_options = list(method = "2SLS")
)
eq_fit <- equilibrium_model(
  market_spec, stochastic_adjustment_data,
  estimation_options = list(standard_errors = c("id"))
)
bs_fit <- diseq_basic(
  market_spec, stochastic_adjustment_data,
  estimation_options = list(
    method = "Nelder-Mead", control = list(maxit = 1e+5),
    standard_errors = "heteroscedastic"
  )
)
dr_fit <- diseq_directional(
  formula(update(Formula(market_spec), . ~ . | . - P)),
  stochastic_adjustment_data,
  estimation_options = list(standard_errors = "heteroscedastic")
)
da_fit <- diseq_deterministic_adjustment(
  market_spec, stochastic_adjustment_data,
  estimation_options = list(standard_errors = c("id"))
)
sa_fit <- diseq_stochastic_adjustment(
  formula(update(Formula(market_spec), . ~ . | . | Xp1)),
  stochastic_adjustment_data,
  estimation_options = list(control = list(maxit = 1e+5))
)
```


## Post estimation analysis 

### Summaries

All the model estimates support the `summary` function. The `eq_2sls` also provides the first-stage estimation, but it is not included in the summary and has to be explicitly asked. 

```{r analysis.summaries}
summary(eq_reg)
summary(eq_fit)
summary(bs_fit)
summary(da_fit)
summary(sa_fit)
```

### Marginal effects

Calculate marginal effects on the shortage probabilities. *Markets* offers two marginal effect calls out of the box. The mean marginal effects and the marginal effects ate the mean. Marginal effects on the shortage probabilities are state-dependent. If the variable is only in the demand equation, the output name of the marginal effect is the variable name prefixed by `D_`. If the variable is only in the supply equation, the name of the marginal effect is the variable name prefixed by `S_`. If the variable is in both equations, then it is prefixed by `B_`.

```{r analysis.effects}
diseq_abbrs <- c("bs", "dr", "da", "sa")
diseq_fits <- c(bs_fit, dr_fit, da_fit, sa_fit)
variables <- c("P", "Xd1", "Xd2", "X1", "X2", "Xs1")

apply_marginal <- function(fnc, ...) {
  function(fit) {
    sapply(variables, function(v) fnc(fit, v, ...), USE.NAMES = FALSE)
  }
}

mspm <- sapply(diseq_fits, apply_marginal(shortage_probability_marginal))
colnames(mspm) <- diseq_abbrs
# Mean Shortage Probabilities' Marginal Effects
mspm

spmm <- sapply(
  diseq_fits,
  apply_marginal(shortage_probability_marginal, aggregate = "at_the_mean")
)
colnames(spmm) <- diseq_abbrs
# Shortage Probabilities' Marginal Effects at the Mean
spmm
```

### Shortages

Copy the disequilibrium model data frame and augment it with post-estimation data. The disequilibrium models can be used to estimate:

 * Shortage probabilities. These are the probabilities that the disequilibrium models assign to observing a particular extent of excess demand.
 
 * Normalized shortages. The point estimates of the shortages are normalized by the variance of the difference of demand and supply shocks.
 
 * Relative shortages: The point estimates of the shortages are normalized by the estimated supplied quantity.

```{r analysis.estimates}
fit <- sa_fit
mdt <- cbind(
  fit@model@data,
  shortage_indicators = c(shortage_indicators(fit)),
  normalized_shortages = c(normalized_shortages(fit)),
  shortage_probabilities = c(shortage_probabilities(fit)),
  relative_shortages = c(relative_shortages(fit))
)
```

How is the sample separated post-estimation? The indices of the observations for which the estimated demand is greater than the estimated supply are easily obtained.

```{r analysis.shortages}
if (requireNamespace("ggplot2", quietly = TRUE)) {
  pdt <- data.frame(
    Shortage = c(mdt$normalized_shortages, mdt$relative_shortages),
    Type = c(rep("Normalized", nrow(mdt)), rep("Relative", nrow(mdt))),
    xpos = c(rep(-1.0, nrow(mdt)), rep(1.0, nrow(mdt))),
    ypos = c(
      rep(0.8 * max(mdt$normalized_shortages), nrow(mdt)),
      rep(0.8 * max(mdt$relative_shortages), nrow(mdt))
    )
  )
  ggplot2::ggplot(pdt) +
    ggplot2::stat_density(ggplot2::aes(Shortage,
      linetype = Type,
      color = Type
    ), geom = "line") +
    ggplot2::ggtitle(paste0("Estimated shortages densities (", name(fit), ")")) +
    ggplot2::theme(
      panel.background = ggplot2::element_rect(fill = "transparent"),
      plot.background = ggplot2::element_rect(
        fill = "transparent",
        color = NA
      ),
      legend.background = ggplot2::element_rect(fill = "transparent"),
      legend.box.background = ggplot2::element_rect(
        fill = "transparent",
        color = NA
      ),
      legend.position = c(0.8, 0.8)
    )
} else {
    summary(mdt[, grep("shortage", colnames(mdt))])
}
```

### Fitted values and aggregation

The estimated demanded and supplied quantities can be calculated per observation. 

```{r analysis.market_forces}
market <- cbind(
  demand = demanded_quantities(fit)[, 1],
  supply = supplied_quantities(fit)[, 1]
)
summary(market)
```

The package also offers basic aggregation functionality.

```{r analysis.aggregation}
aggregates <- aggregate_demand(fit) |>
  dplyr::left_join(aggregate_supply(fit), by = "date") |>
  dplyr::mutate(date = as.numeric(date)) |>
  dplyr::rename(demand = D_Q, supply = S_Q)
if (requireNamespace("ggplot2", quietly = TRUE)) {
  pdt <- data.frame(
    Date = c(aggregates$date, aggregates$date),
    Quantity = c(aggregates$demand, aggregates$supply),
    Side = c(rep("Demand", nrow(aggregates)), rep("Supply", nrow(aggregates)))
  )
  ggplot2::ggplot(pdt, ggplot2::aes(x = Date)) +
    ggplot2::geom_line(ggplot2::aes(y = Quantity, linetype = Side, color = Side)) +
    ggplot2::ggtitle(paste0(
      "Aggregate estimated demand and supply  (", name(fit), ")"
    )) +
    ggplot2::theme(
      panel.background = ggplot2::element_rect(fill = "transparent"),
      plot.background = ggplot2::element_rect(
        fill = "transparent", color = NA
      ),
      legend.background = ggplot2::element_rect(fill = "transparent"),
      legend.box.background = ggplot2::element_rect(
        fill = "transparent", color = NA
      ),
      legend.position = c(0.8, 0.5)
    )
} else {
    aggregates
}
```