Basic usage

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)

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

eq2sls <- estimate(eq2sls)
eqfiml_est <- estimate(eqfiml,
  control = optimization_controls,
  method = optimization_method
)
bsmdl_est <- estimate(bsmdl,
  control = optimization_controls,
  method = optimization_method
)
drmdl_est <- estimate(drmdl,
  control = optimization_controls,
  method = optimization_method
)
damdl_est <- estimate(damdl,
  control = optimization_controls,
  method = optimization_method
)
samdl_est <- estimate(samdl,
  control = optimization_controls,
  method = optimization_method
)

Post estimation analysis

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

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

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), ")"
    ))
}

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: 0x55fc92ecf938>
#> Instruments: ~Xd1 + Xd2 + X1 + X2 + Xs1
#> <environment: 0x55fc92ecf938>
#> 
#>               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: 0x55fc92ecf938>
#> Instruments: ~Xd1 + Xd2 + X1 + X2 + Xs1
#> <environment: 0x55fc92ecf938>
#> 
#>               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.3159274  0.0091012  34.7128 < 2.2e-16 ***
#> D_CONST     1.7194355  0.2471955   6.9558 3.506e-12 ***
#> D_Xd1       0.0461226  0.0066524   6.9332 4.114e-12 ***
#> D_Xd2      -0.0269238  0.0062343  -4.3187 1.570e-05 ***
#> D_X1       -0.2250904  0.0478912  -4.7000 2.601e-06 ***
#> D_X2        0.1233735  0.0474833   2.5983  0.009370 ** 
#> S_P         0.4639770  0.0033303 139.3208 < 2.2e-16 ***
#> S_CONST    -0.6522565  0.2319501  -2.8121  0.004923 ** 
#> S_Xs1      -0.2248839  0.0201892 -11.1388 < 2.2e-16 ***
#> S_X1       -0.3016973  0.0691032  -4.3659 1.266e-05 ***
#> S_X2        0.1572989  0.0682277   2.3055  0.021139 *  
#> D_VARIANCE  5.6234917  0.2832408  19.8541 < 2.2e-16 ***
#> S_VARIANCE 11.7120923  0.0845867 138.4626 < 2.2e-16 ***
#> RHO         0.9919654  0.0013962 710.4539 < 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        -0.0676615  0.0058233 -11.6191 < 2.2e-16 ***
#> D_CONST     9.0991117  0.2237501  40.6664 < 2.2e-16 ***
#> D_Xd1       0.3051360  0.0358039   8.5224 < 2.2e-16 ***
#> D_Xd2      -0.1839963  0.0350303  -5.2525 1.501e-07 ***
#> D_X1        0.5776875  0.0517207  11.1694 < 2.2e-16 ***
#> D_X2       -0.1443290  0.0395242  -3.6517 0.0002606 ***
#> S_P         0.0813754  0.0039040  20.8443 < 2.2e-16 ***
#> S_CONST     5.4182457  0.1509925  35.8842 < 2.2e-16 ***
#> S_Xs1       0.9298673  0.0389533  23.8713 < 2.2e-16 ***
#> S_X1       -0.5213550  0.0385037 -13.5404 < 2.2e-16 ***
#> S_X2        0.2025852  0.0327272   6.1901 6.012e-10 ***
#> D_VARIANCE  0.9286177  0.0432209  21.4854 < 2.2e-16 ***
#> S_VARIANCE  0.9959581  0.0321953  30.9348 < 2.2e-16 ***
#> RHO         0.1731088  0.0779616   2.2204 0.0263891 *  
#> ---
#> 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.0066225 -15.1553 < 2.2e-16 ***
#> D_CONST    12.3271066  0.2957220  41.6848 < 2.2e-16 ***
#> D_Xd1       0.1511296  0.0237260   6.3698 1.893e-10 ***
#> D_Xd2      -0.0988986  0.0224765  -4.4001 1.082e-05 ***
#> D_X1        0.6028142  0.0452146  13.3323 < 2.2e-16 ***
#> D_X2       -0.0992912  0.0300319  -3.3062 0.0009457 ***
#> S_P         0.0012140  0.0018120   0.6700 0.5028852    
#> S_CONST     6.8681300  0.0979114  70.1464 < 2.2e-16 ***
#> S_Xs1       0.4690616  0.0200247  23.4242 < 2.2e-16 ***
#> S_X1        0.0433718  0.0200539   2.1628 0.0305596 *  
#> S_X2        0.0269852  0.0200151   1.3482 0.1775806    
#> P_DIFF      0.6093889  0.0376155  16.2005 < 2.2e-16 ***
#> D_VARIANCE  1.6865709  0.1144416  14.7374 < 2.2e-16 ***
#> S_VARIANCE  0.7979176  0.0126176  63.2385 < 2.2e-16 ***
#> RHO         0.6725925  0.0242021  27.7906 < 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.0042877 -22.6535 < 2.2e-16 ***
#> D_CONST     9.7139866  0.1620684  59.9376 < 2.2e-16 ***
#> D_Xd1       0.2950376  0.0265935  11.0943 < 2.2e-16 ***
#> D_Xd2      -0.1722320  0.0262833  -6.5529 5.643e-11 ***
#> D_X1        0.5641836  0.0321914  17.5259 < 2.2e-16 ***
#> D_X2       -0.0869383  0.0262036  -3.3178 0.0009073 ***
#> S_P         0.0989968  0.0052260  18.9433 < 2.2e-16 ***
#> S_CONST     5.1636578  0.1499040  34.4464 < 2.2e-16 ***
#> S_Xs1       0.9380370  0.0346449  27.0758 < 2.2e-16 ***
#> S_X1       -0.5153828  0.0363245 -14.1883 < 2.2e-16 ***
#> S_X2        0.1660077  0.0277562   5.9809 2.219e-09 ***
#> P_DIFF      1.1852181  0.0347836  34.0741 < 2.2e-16 ***
#> P_CONST     3.2090705  0.0968038  33.1502 < 2.2e-16 ***
#> P_Xp1       0.7695691  0.0310035  24.8220 < 2.2e-16 ***
#> D_VARIANCE  0.9769727  0.0373513  26.1563 < 2.2e-16 ***
#> S_VARIANCE  1.0493231  0.0457338  22.9442 < 2.2e-16 ***
#> P_VARIANCE  0.9399843  0.0731163  12.8560 < 2.2e-16 ***
#> RHO_DS      0.0389034  0.0630691   0.6168 0.5373422    
#> RHO_DP      0.0029079  0.0492433   0.0591 0.9529109    
#> RHO_SP      0.0013893  0.0476685   0.0291 0.9767481    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> -2 log L: 46400.96

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

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