# Some examples of use-cases using the market models of the package

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

## Setup the environment

library(diseq)
library(magrittr)
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 diseq 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.

nobs <- 1000
tobs <- 10

alpha_d <- -0.3
beta_d0 <- 6.8
beta_d <- c(0.3, -0.02)
eta_d <- c(0.6, -0.1)

alpha_s <- 0.6
beta_s0 <- 4.1
beta_s <- c(0.9)
eta_s <- c(-0.5, 0.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

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.

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.

eq_reg <- equilibrium_model(
estimation_options = list(method = "2SLS")
)
eq_fit <- equilibrium_model(
estimation_options = list(standard_errors = c("id"))
)
bs_fit <- diseq_basic(
estimation_options = list(
method = "Nelder-Mead", control = list(maxit = 1e+5),
standard_errors = "heteroscedastic"
)
)
dr_fit <- diseq_directional(
formula(update(Formula(market_spec), . ~ . | . - P)),
estimation_options = list(standard_errors = "heteroscedastic")
)
estimation_options = list(standard_errors = c("id"))
)
formula(update(Formula(market_spec), . ~ . | . | Xp1)),
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.

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 (", model_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. market <- cbind( demand = demanded_quantities(fit)[, 1], supply = supplied_quantities(fit)[, 1] ) summary(market) #> demand supply #> Min. :3.175 Min. : 2.529 #> 1st Qu.:5.241 1st Qu.: 5.670 #> Median :5.684 Median : 6.366 #> Mean :5.683 Mean : 6.362 #> 3rd Qu.:6.129 3rd Qu.: 7.054 #> Max. :8.148 Max. :10.747 The package also offers basic aggregation functionality. 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 <- tibble::tibble( 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  (", model_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
}