An example of estimations 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

Load the required libraries.

library(diseq)
library(magrittr)

Prepare the data. Normally this step is long and it 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_model_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 <- 2000
tobs <- 5

alpha_d <- -0.1
beta_d0 <- 9.8
beta_d <- c(0.3, -0.2)
eta_d <- c(0.6, -0.1)

alpha_s <- 0.1
beta_s0 <- 5.1
beta_s <- c(0.9)
eta_s <- c(-0.5, 0.2)

gamma <- 1.2
beta_p0 <- 3.1
beta_p <- c(0.8)

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

seed <- 443

stochastic_adjustment_data <- simulate_model_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
)

Initialize the model

The constructor sets the basic parameters for model initialization and constructs a model object. The needed arguments for a construction call are configured as follows:

Set the fields that uniquely identify data rows. For panel data this should be a vector of the entity identifier and the time columns.

key_columns <- c("id", "date")

Set time field. This is needed for the calculation of price’s first differences. The diseq_directional, diseq_deterministic_adjustment, and diseq_stochastic_adjustment models require calculating the first differences.

time_column <- c("date")

Set the quantity variable. The quantity variable is the observable. For the equilibrium models, it is equal with both the demanded and supplied quantities. For the disequilibrium models, the observed quantity represents either a demanded or a supplied quantity. Each disequilibrium model resolves that state of the observation in a different way.

quantity_column <- "Q"

Set the price variable.

price_column <- "P"

Set the remaining demand and supply variables. Simply include the factor variables here as in a usual lm formula. Indicator variables and interaction terms will be created automatically by the constructor. For the diseq_directional model, the price cannot go in both equations. For the rest of the models, the price can go in both equations if treated as exogenous. The diseq_stochastic_adjustment requires also the specification of price dynamics.

demand_specification <- paste0(price_column, " + Xd1 + Xd2 + X1 + X2")
supply_specification <- "Xs1 + X1 + X2"
price_specification <- "Xp1"

Set the verbosity level. This controls the level of messaging. The object displays
- error: always,
- warning: . \(\ge\) 1,
- info: . \(\ge\) 2,
- verbose: . \(\ge\) 3 and
- debug :. \(\ge\) 4.

verbose <- 2

Should the estimation allow for correlated demand and supply shocks?

use_correlated_shocks <- TRUE

Using the above parameterization, construct the model objects. Here, we construct two equilibrium models and and four disequilibrium models, using in all cases the same data simulated by the process based on the stochastic price adjustment model. Of course, this is only done to simplify the exposition of the functionality. The constructors of the models that use price dynamics information in the estimation, i.e. diseq_directional, diseq_deterministic_adjustment, and diseq_stochastic_adjustment, will automatically generate lagged prices and drop one observation per entity.

eq2sls <- new(
  "eq_2sls",
  key_columns,
  quantity_column, price_column,
  demand_specification, paste0(price_column, " + ", supply_specification),
  stochastic_adjustment_data,
  verbose = verbose
)
#> Info: This is 'Equilibrium 2SLS with independent shocks' model
eqfiml <- new(
  "eq_fiml",
  key_columns,
  quantity_column, price_column,
  demand_specification, paste0(price_column, " + ", supply_specification),
  stochastic_adjustment_data,
  use_correlated_shocks = use_correlated_shocks, verbose = verbose
)
#> Info: This is 'Equilibrium FIML with correlated shocks' model
bsmdl <- new(
  "diseq_basic",
  key_columns,
  quantity_column, price_column,
  demand_specification, paste0(price_column, " + ", supply_specification),
  stochastic_adjustment_data,
  use_correlated_shocks = use_correlated_shocks, verbose = verbose
)
#> Info: This is 'Basic with correlated shocks' model
drmdl <- new(
  "diseq_directional",
  key_columns, time_column,
  quantity_column, price_column,
  demand_specification, supply_specification,
  stochastic_adjustment_data,
  use_correlated_shocks = use_correlated_shocks, verbose = verbose
)
#> Info: This is 'Directional with correlated shocks' model
#> Info: Dropping 2000 rows by generating 'LAGGED_P'.
#> Info: Sample separated with 9 rows in excess supply and 7991 in excess demand regime.
damdl <- new(
  "diseq_deterministic_adjustment",
  key_columns, time_column,
  quantity_column, price_column,
  demand_specification, paste0(price_column, " + ", supply_specification),
  stochastic_adjustment_data,
  use_correlated_shocks = use_correlated_shocks, verbose = verbose
)
#> Info: This is 'Deterministic Adjustment with correlated shocks' model
#> Info: Dropping 2000 rows by generating 'LAGGED_P'.
#> Info: Sample separated with 9 rows in excess supply and 7991 in excess demand regime.
samdl <- new(
  "diseq_stochastic_adjustment",
  key_columns, time_column,
  quantity_column, price_column,
  demand_specification, paste0(price_column, " + ", supply_specification),
  price_specification,
  stochastic_adjustment_data,
  use_correlated_shocks = use_correlated_shocks, verbose = verbose
)
#> Info: This is 'Stochastic Adjustment with correlated shocks' model
#> Info: Dropping 2000 rows by generating 'LAGGED_P'.

Estimation

Set the estimation parameters. The only model that is estimated by least squares is eq_2sls. The remaining models are estimated using full information maximum likelihood. First choose an optimization method and the corresponding optimization controls. The available methods are:

Nelder-Mead: Does not require the gradient of the likelihood to be known.
BFGS: Uses the analytically calculated gradients. By default the diseq package uses this method.
L-BFGS-B: Constrained optimization.

optimization_method <- "BFGS"
optimization_controls <- list(REPORT = 10, maxit = 10000, reltol = 1e-6)
use_heteroscedasticity_consistent_errors <- TRUE
cluster_errors_by <- c("id")

Then, estimate the models. See the documentation for more options.

eq2sls <- estimate(eq2sls)
eqfiml_est <- estimate(eqfiml,
  control = optimization_controls, method = optimization_method,
  cluster_errors_by = cluster_errors_by
)
#> Warning: `summarise_each_()` is deprecated as of dplyr 0.7.0.
#> Please use `across()` instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_warnings()` to see where this warning was generated.
#> Warning: `as.tibble()` is deprecated as of tibble 2.0.0.
#> Please use `as_tibble()` instead.
#> The signature and semantics have changed, see `?as_tibble`.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_warnings()` to see where this warning was generated.
#> Warning: `funs()` is deprecated as of dplyr 0.8.0.
#> Please use a list of either functions or lambdas: 
#> 
#>   # Simple named list: 
#>   list(mean = mean, median = median)
#> 
#>   # Auto named with `tibble::lst()`: 
#>   tibble::lst(mean, median)
#> 
#>   # Using lambdas
#>   list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_warnings()` to see where this warning was generated.
bsmdl_est <- estimate(bsmdl,
  control = optimization_controls, method = optimization_method,
  use_heteroscedasticity_consistent_errors = use_heteroscedasticity_consistent_errors
)
drmdl_est <- estimate(drmdl,
  control = optimization_controls, method = optimization_method,
  use_heteroscedasticity_consistent_errors = use_heteroscedasticity_consistent_errors
)
damdl_est <- estimate(damdl,
  control = optimization_controls, method = optimization_method,
  cluster_errors_by = cluster_errors_by
)
samdl_est <- estimate(samdl,
  control = optimization_controls, method = optimization_method,
  cluster_errors_by = cluster_errors_by
)

Post estimation analysis

Marginal effects

Calculate marginal effects on the shortage probabilities. Diseq 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_.

variables <- c(price_column, "Xd1", "Xd2", "X1", "X2", "Xs1")

bsmdl_mme <- sapply(variables,
  function(v) get_mean_marginal_effect(bsmdl, bsmdl_est, v),
  USE.NAMES = FALSE
)
drmdl_mme <- sapply(variables,
  function(v) get_mean_marginal_effect(drmdl, drmdl_est, v),
  USE.NAMES = FALSE
)
damdl_mme <- sapply(variables,
  function(v) get_mean_marginal_effect(damdl, damdl_est, v),
  USE.NAMES = FALSE
)
samdl_mme <- sapply(variables,
  function(v) get_mean_marginal_effect(samdl, samdl_est, v),
  USE.NAMES = FALSE
)
bsmdl_mem <- sapply(variables,
  function(v) get_marginal_effect_at_mean(bsmdl, bsmdl_est, v),
  USE.NAMES = FALSE
)
drmdl_mem <- sapply(variables,
  function(v) get_marginal_effect_at_mean(drmdl, drmdl_est, v),
  USE.NAMES = FALSE
)
damdl_mem <- sapply(variables,
  function(v) get_marginal_effect_at_mean(damdl, damdl_est, v),
  USE.NAMES = FALSE
)
samdl_mem <- sapply(variables,
  function(v) get_marginal_effect_at_mean(samdl, samdl_est, v),
  USE.NAMES = FALSE
)

cbind(
  bsmdl_mme, drmdl_mme, damdl_mme, samdl_mme, bsmdl_mem, drmdl_mem, damdl_mem,
  samdl_mem
)
#>         bsmdl_mme     drmdl_mme     damdl_mme   samdl_mme   bsmdl_mem
#> B_P   -0.03116252 -0.0001311817 -0.0009663191 -0.04037366 -0.04445984
#> D_Xd1  0.06380171  0.0003380489  0.0014376743  0.06073446  0.09102643
#> D_Xd2 -0.03847228  0.0001247948 -0.0009408077 -0.03545453 -0.05488873
#> B_X1   0.22980175  0.0009789019  0.0053218947  0.22223228  0.32786008
#> B_X2  -0.07253721  0.0020569042 -0.0012012494 -0.05206977 -0.10348944
#> S_Xs1 -0.19442845 -0.0006311324 -0.0044621154 -0.19309798 -0.27739271
#>           drmdl_mem     damdl_mem   samdl_mem
#> B_P   -6.251634e-05 -0.0002319415 -0.05602628
#> D_Xd1  1.611016e-04  0.0003450789  0.08428082
#> D_Xd2  5.947257e-05 -0.0002258181 -0.04920002
#> B_X1   4.665084e-04  0.0012773920  0.30839032
#> B_X2   9.802444e-04 -0.0002883308 -0.07225689
#> S_Xs1 -3.007744e-04 -0.0010710228 -0.26796084

Shortages

Copy the disequilibrium model tibble 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 the shocks of demand and supply.
Relative shortages: The point estimates of the shortages are normalized by the estimated supplied quantity.

mdt <- tibble::add_column(
  bsmdl@model_tibble,
  normalized_shortages = c(get_normalized_shortages(bsmdl, bsmdl_est@coef)),
  shortage_probabilities = c(get_shortage_probabilities(bsmdl, bsmdl_est@coef)),
  relative_shortages = c(get_relative_shortages(bsmdl, bsmdl_est@coef))
)

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.

abs_estsep <- c(
  nobs = length(has_shortage(bsmdl, bsmdl_est@coef)),
  nshortages = sum(has_shortage(bsmdl, bsmdl_est@coef)),
  nsurpluses = sum(!has_shortage(bsmdl, bsmdl_est@coef))
)
print(abs_estsep)
#>       nobs nshortages nsurpluses 
#>      10000       6005       3995

rel_estsep <- abs_estsep / abs_estsep["nobs"]
names(rel_estsep) <- c("total", "shortages_share", "surpluses_share")
print(rel_estsep)
#>           total shortages_share surpluses_share 
#>          1.0000          0.6005          0.3995

if (requireNamespace("ggplot2", quietly = TRUE)) {
  ggplot2::ggplot(mdt, ggplot2::aes(normalized_shortages)) +
    ggplot2::geom_density() +
    ggplot2::ggtitle(paste0(
      "Normalized shortages density (",
      get_model_description(bsmdl), ")"
    ))
}

Summaries

All the model estimates support the summary function. The eq_2sls provides also the first stage estimation. The summary output comes from systemfit. The remaining models are estimated using maximum likelihood and the summary functionality is based on bbmle.

summary(eq2sls@first_stage_model)
#> 
#> Call:
#> lm(formula = first_stage_formula, data = object@model_tibble)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -16.6628  -5.8648   0.6677   6.2961  18.0452 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) 20.63325    0.84994  24.276  < 2e-16 ***
#> Xd1          0.35530    0.15154   2.345   0.0191 *  
#> Xd2         -0.06874    0.15005  -0.458   0.6469    
#> X1           0.58717    0.15052   3.901 9.65e-05 ***
#> X2          -0.18630    0.15021  -1.240   0.2149    
#> Xs1         -0.60577    0.15126  -4.005 6.25e-05 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 7.503 on 9994 degrees of freedom
#> Multiple R-squared:  0.003851,   Adjusted R-squared:  0.003353 
#> F-statistic: 7.728 on 5 and 9994 DF,  p-value: 2.898e-07
summary(eq2sls@system_model)
#> 
#> systemfit results 
#> method: 2SLS 
#> 
#>            N    DF     SSR  detRCov   OLS-R2 McElroy-R2
#> system 20000 19989 17839.9 0.000636 0.083066  -0.833223
#> 
#>            N   DF     SSR      MSE    RMSE       R2   Adj R2
#> demand 10000 9994 8916.39 0.892175 0.94455 0.083433 0.082974
#> supply 10000 9995 8923.52 0.892799 0.94488 0.082700 0.082333
#> 
#> The covariance matrix of the residuals
#>          demand   supply
#> demand 0.892175 0.892130
#> supply 0.892130 0.892799
#> 
#> The correlations of the residuals
#>        demand supply
#> demand 1.0000 0.9996
#> supply 0.9996 1.0000
#> 
#> 
#> 2SLS estimates for 'demand' (equation 1)
#> Model Formula: Q ~ P_FITTED + Xd1 + Xd2 + X1 + X2
#> <environment: 0x55cc69aead50>
#> Instruments: ~Xd1 + Xd2 + X1 + X2 + Xs1
#> <environment: 0x55cc69aead50>
#> 
#>               Estimate Std. Error   t value   Pr(>|t|)    
#> (Intercept) 25.1518990  0.6082042  41.35437 < 2.22e-16 ***
#> P_FITTED    -0.8989314  0.0314346 -28.59686 < 2.22e-16 ***
#> Xd1          0.4679311  0.0221114  21.16244 < 2.22e-16 ***
#> Xd2         -0.1449714  0.0190430  -7.61285  2.931e-14 ***
#> X1           0.4682736  0.0264901  17.67734 < 2.22e-16 ***
#> X2          -0.1159127  0.0196993  -5.88409  4.131e-09 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.94455 on 9994 degrees of freedom
#> Number of observations: 10000 Degrees of Freedom: 9994 
#> SSR: 8916.394146 MSE: 0.892175 Root MSE: 0.94455 
#> Multiple R-Squared: 0.083433 Adjusted R-Squared: 0.082974 
#> 
#> 
#> 2SLS estimates for 'supply' (equation 2)
#> Model Formula: Q ~ P_FITTED + Xs1 + X1 + X2
#> <environment: 0x55cc69aead50>
#> Instruments: ~Xd1 + Xd2 + X1 + X2 + Xs1
#> <environment: 0x55cc69aead50>
#> 
#>               Estimate Std. Error  t value   Pr(>|t|)    
#> (Intercept) -2.7883417  1.1282336 -2.47142   0.013474 *  
#> P_FITTED     0.4476886  0.0526790  8.49842 < 2.22e-16 ***
#> Xs1          0.8150032  0.0372171 21.89861 < 2.22e-16 ***
#> X1          -0.3225859  0.0362761 -8.89252 < 2.22e-16 ***
#> X2           0.1349937  0.0213726  6.31620 2.7948e-10 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.94488 on 9995 degrees of freedom
#> Number of observations: 10000 Degrees of Freedom: 9995 
#> SSR: 8923.523074 MSE: 0.892799 Root MSE: 0.94488 
#> Multiple R-Squared: 0.0827 Adjusted R-Squared: 0.082333
bbmle::summary(eqfiml_est)
#> Maximum likelihood estimation
#> 
#> Call:
#> `bbmle::mle2`(list(control = list(REPORT = 10, maxit = 10000, 
#>     reltol = 1e-06), method = "BFGS", start = c(D_P = 0.0231239615913565, 
#> D_CONST = 7.53375095184334, D_Xd1 = 0.140015230368245, D_Xd2 = -0.0737970096786045, 
#> D_X1 = -0.074691417704378, D_X2 = 0.0462502352837326, S_P = 0.0248455188388157, 
#> S_CONST = 6.24286219183355, S_Xs1 = 0.558355600251176, S_X1 = -0.0743221338857006, 
#> S_X2 = 0.055148632519299, D_VARIANCE = 1, S_VARIANCE = 1, RHO = 0
#> ), minuslogl = function (...) 
#> minus_log_likelihood(object, ...), gr = function (...) 
#> gradient(object, ...)))
#> 
#> Coefficients:
#>               Estimate  Std. Error   z value     Pr(z)    
#> D_P         0.31592736  0.00456872   69.1501 < 2.2e-16 ***
#> D_CONST     1.71943548  0.05532016   31.0815 < 2.2e-16 ***
#> D_Xd1       0.04612264  0.00649816    7.0978 1.268e-12 ***
#> D_Xd2      -0.02692379  0.00614201   -4.3835 1.168e-05 ***
#> D_X1       -0.22509044  0.04705933   -4.7831 1.726e-06 ***
#> D_X2        0.12337354  0.04366452    2.8255 0.0047209 ** 
#> S_P         0.46397700  0.00202615  228.9945 < 2.2e-16 ***
#> S_CONST    -0.65225650  0.16817571   -3.8784 0.0001051 ***
#> S_Xs1      -0.22488392  0.01061251  -21.1905 < 2.2e-16 ***
#> S_X1       -0.30169732  0.06588591   -4.5791 4.670e-06 ***
#> S_X2        0.15729886  0.06348024    2.4779 0.0132151 *  
#> D_VARIANCE  5.62349166  0.14559276   38.6248 < 2.2e-16 ***
#> S_VARIANCE 11.71209229  0.05019747  233.3204 < 2.2e-16 ***
#> RHO         0.99196541  0.00062413 1589.3453 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> -2 log L: 95759.78
bbmle::summary(bsmdl_est)
#> Maximum likelihood estimation
#> 
#> Call:
#> `bbmle::mle2`(list(control = list(REPORT = 10, maxit = 10000, 
#>     reltol = 1e-06), method = "BFGS", skip.hessian = FALSE, start = c(D_P = 0.0231239615913565, 
#> D_CONST = 7.53375095184334, D_Xd1 = 0.140015230368245, D_Xd2 = -0.0737970096786045, 
#> D_X1 = -0.074691417704378, D_X2 = 0.0462502352837326, S_P = 0.0248455188388157, 
#> S_CONST = 6.24286219183355, S_Xs1 = 0.558355600251176, S_X1 = -0.0743221338857006, 
#> S_X2 = 0.055148632519299, D_VARIANCE = 1, S_VARIANCE = 1, RHO = 0
#> ), minuslogl = function (...) 
#> minus_log_likelihood(object, ...), gr = function (...) 
#> gradient(object, ...)))
#> 
#> Coefficients:
#>               Estimate  Std. Error   z value     Pr(z)    
#> D_P        -6.7662e-02  5.3178e-05 -1272.351 < 2.2e-16 ***
#> D_CONST     9.0991e+00  2.3479e-04 38753.539 < 2.2e-16 ***
#> D_Xd1       3.0514e-01  4.0574e-04   752.051 < 2.2e-16 ***
#> D_Xd2      -1.8400e-01  7.6105e-04  -241.766 < 2.2e-16 ***
#> D_X1        5.7769e-01  1.3951e-04  4140.954 < 2.2e-16 ***
#> D_X2       -1.4433e-01  1.0246e-03  -140.865 < 2.2e-16 ***
#> S_P         8.1375e-02  2.1225e-03    38.340 < 2.2e-16 ***
#> S_CONST     5.4182e+00  1.6727e-03  3239.296 < 2.2e-16 ***
#> S_Xs1       9.2987e-01  5.2691e-03   176.475 < 2.2e-16 ***
#> S_X1       -5.2135e-01  3.8454e-03  -135.578 < 2.2e-16 ***
#> S_X2        2.0259e-01  5.0656e-03    39.992 < 2.2e-16 ***
#> D_VARIANCE  9.2862e-01  3.6439e-04  2548.418 < 2.2e-16 ***
#> S_VARIANCE  9.9596e-01  4.9565e-05 20094.056 < 2.2e-16 ***
#> RHO         1.7311e-01  6.8595e-05  2523.626 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> -2 log L: 26189.5
bbmle::summary(damdl_est)
#> Maximum likelihood estimation
#> 
#> Call:
#> `bbmle::mle2`(list(control = list(REPORT = 10, maxit = 10000, 
#>     reltol = 1e-06), method = "BFGS", skip.hessian = FALSE, start = c(D_P = -0.00260947112358515, 
#> D_CONST = 7.98130444023755, D_Xd1 = 0.15954441702268, D_Xd2 = -0.0956278531130884, 
#> D_X1 = 0.0454572982786684, D_X2 = 0.0215741337521343, S_P = 0.000425847047005336, 
#> S_CONST = 6.8876653018125, S_Xs1 = 0.467154949725634, S_X1 = 0.0442638269176111, 
#> S_X2 = 0.0274612448510507, P_DIFF = 1, D_VARIANCE = 1, S_VARIANCE = 1, 
#> RHO = 0), minuslogl = function (...) 
#> minus_log_likelihood(object, ...), gr = function (...) 
#> gradient(object, ...)))
#> 
#> Coefficients:
#>              Estimate Std. Error  z value     Pr(z)    
#> D_P        -0.1003664  0.0064966 -15.4492 < 2.2e-16 ***
#> D_CONST    12.3271066  0.2892418  42.6187 < 2.2e-16 ***
#> D_Xd1       0.1511296  0.0237067   6.3750 1.830e-10 ***
#> D_Xd2      -0.0988986  0.0222613  -4.4426 8.887e-06 ***
#> D_X1        0.6028142  0.0464252  12.9846 < 2.2e-16 ***
#> D_X2       -0.0992912  0.0308558  -3.2179  0.001291 ** 
#> S_P         0.0012140  0.0018291   0.6637  0.506887    
#> S_CONST     6.8681300  0.0972483  70.6247 < 2.2e-16 ***
#> S_Xs1       0.4690616  0.0202024  23.2181 < 2.2e-16 ***
#> S_X1        0.0433718  0.0205144   2.1142  0.034497 *  
#> S_X2        0.0269852  0.0199371   1.3535  0.175891    
#> P_DIFF      0.6093889  0.0375820  16.2149 < 2.2e-16 ***
#> D_VARIANCE  1.6865709  0.1125298  14.9878 < 2.2e-16 ***
#> S_VARIANCE  0.7979176  0.0125662  63.4973 < 2.2e-16 ***
#> RHO         0.6725925  0.0247984  27.1224 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> -2 log L: 48414.93
bbmle::summary(samdl_est)
#> Maximum likelihood estimation
#> 
#> Call:
#> `bbmle::mle2`(list(control = list(REPORT = 10, maxit = 10000, 
#>     reltol = 1e-06), method = "BFGS", skip.hessian = FALSE, start = c(D_P = -0.00260947112358515, 
#> D_CONST = 7.98130444023755, D_Xd1 = 0.15954441702268, D_Xd2 = -0.0956278531130884, 
#> D_X1 = 0.0454572982786684, D_X2 = 0.0215741337521343, S_P = 0.000425847047005336, 
#> S_CONST = 6.8876653018125, S_Xs1 = 0.467154949725634, S_X1 = 0.0442638269176111, 
#> S_X2 = 0.0274612448510507, P_DIFF = 1, P_CONST = 8.28009900264685, 
#> P_Xp1 = -0.013207398816623, D_VARIANCE = 1, S_VARIANCE = 1, P_VARIANCE = 1, 
#> RHO_DS = 0, RHO_DP = 0, RHO_SP = 0), minuslogl = function (...) 
#> minus_log_likelihood(object, ...), gr = function (...) 
#> gradient(object, ...)))
#> 
#> Coefficients:
#>              Estimate Std. Error  z value     Pr(z)    
#> D_P        -0.0971316  0.0042175 -23.0303 < 2.2e-16 ***
#> D_CONST     9.7139866  0.1608675  60.3850 < 2.2e-16 ***
#> D_Xd1       0.2950376  0.0264703  11.1460 < 2.2e-16 ***
#> D_Xd2      -0.1722320  0.0265270  -6.4927 8.431e-11 ***
#> D_X1        0.5641836  0.0324945  17.3624 < 2.2e-16 ***
#> D_X2       -0.0869383  0.0266174  -3.2662   0.00109 ** 
#> S_P         0.0989968  0.0052399  18.8930 < 2.2e-16 ***
#> S_CONST     5.1636578  0.1528269  33.7876 < 2.2e-16 ***
#> S_Xs1       0.9380370  0.0347249  27.0133 < 2.2e-16 ***
#> S_X1       -0.5153828  0.0359326 -14.3430 < 2.2e-16 ***
#> S_X2        0.1660077  0.0274347   6.0510 1.439e-09 ***
#> P_DIFF      1.1852181  0.0332512  35.6443 < 2.2e-16 ***
#> P_CONST     3.2090705  0.0998930  32.1251 < 2.2e-16 ***
#> P_Xp1       0.7695691  0.0316378  24.3243 < 2.2e-16 ***
#> D_VARIANCE  0.9769727  0.0360423  27.1063 < 2.2e-16 ***
#> S_VARIANCE  1.0493231  0.0446523  23.4999 < 2.2e-16 ***
#> P_VARIANCE  0.9399843  0.0675561  13.9141 < 2.2e-16 ***
#> RHO_DS      0.0389034  0.0595109   0.6537   0.51329    
#> RHO_DP      0.0029079  0.0463269   0.0628   0.94995    
#> RHO_SP      0.0013893  0.0458459   0.0303   0.97582    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> -2 log L: 46400.96

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

market <- cbind(
  demand = get_demanded_quantities(bsmdl, bsmdl_est@coef)[, 1],
  supply = get_supplied_quantities(bsmdl, bsmdl_est@coef)[, 1]
)
summary(market)
#>      demand           supply      
#>  Min.   : 7.133   Min.   : 6.099  
#>  1st Qu.: 8.625   1st Qu.: 8.078  
#>  Median : 9.050   Median : 8.664  
#>  Mean   : 9.075   Mean   : 8.644  
#>  3rd Qu.: 9.520   3rd Qu.: 9.220  
#>  Max.   :11.087   Max.   :11.221

The package offers also basic aggregation functionality.

aggregates <- c(
  demand = get_aggregate_demand(bsmdl, bsmdl_est@coef),
  supply = get_aggregate_supply(bsmdl, bsmdl_est@coef)
)
aggregates
#>   demand   supply 
#> 90746.69 86439.00