Basic usage

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:

key_columns <- c("id", "date")
time_column <- c("date")
quantity_column <- "Q"
price_column <- "P"
demand_specification <- paste0(price_column, " + Xd1 + Xd2 + X1 + X2")
supply_specification <- "Xs1 + X1 + X2"
price_specification <- "Xp1"
verbose <- 2
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:

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

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:

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 = [email protected]_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: 0x55ee9875c8e0>
#> Instruments: ~Xd1 + Xd2 + X1 + X2 + Xs1
#> <environment: 0x55ee9875c8e0>
#> 
#>               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: 0x55ee9875c8e0>
#> Instruments: ~Xd1 + Xd2 + X1 + X2 + Xs1
#> <environment: 0x55ee9875c8e0>
#> 
#>               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

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