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 markets. More usage details can be found in the documentation of the package.

Setup the environment

Load the required 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.

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.

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(
  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.

summary(eq_reg)
#> Equilibrium Model for Markets in Equilibrium:
#>   Demand RHS        :   D_P + D_Xd1 + D_Xd2 + D_X1 + D_X2
#>   Supply RHS        :   S_P + S_Xs1 + S_X1 + S_X2
#>   Market Clearing   : Q = D_Q = S_Q
#>   Shocks            : Correlated
#>   Nobs              : 5000
#>   Sample Separation : Not Separated
#>   Quantity Var      : Q
#>   Price Var         : P
#>   Key Var(s)        : id, date
#>   Time Var          : date
#> 
#> Least squares estimation:
#>   Method              : 2SLS
#> 
#> Shocks:
#>   D_VARIANCE          : 5.66427
#>   S_VARIANCE          : 5.70121
#>   RHO                 : -0.396832
#> 
#> First Stage:
#> 
#> Call:
#> lm(formula = first_stage_formula, data = object@data)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -3.9385 -0.9419 -0.0076  0.8862  4.5895 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)  4.49078    0.01802  249.14   <2e-16 ***
#> Xd1          0.72818    0.01799   40.48   <2e-16 ***
#> Xd2         -0.33706    0.01797  -18.75   <2e-16 ***
#> X1           1.31269    0.01791   73.29   <2e-16 ***
#> X2          -1.36787    0.01829  -74.78   <2e-16 ***
#> Xs1         -0.93267    0.01792  -52.05   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 1.274 on 4994 degrees of freedom
#> Multiple R-squared:  0.7616, Adjusted R-squared:  0.7614 
#> F-statistic:  3191 on 5 and 4994 DF,  p-value: < 2.2e-16
#> 
#> 
#> Demand Equation:
#> 
#> Call:
#> lm(formula = demand_formula, data = object@data)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -7.3458 -0.8461  0.1607  0.9935  3.9047 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) 26.05365    0.09618  270.89   <2e-16 ***
#> P_FITTED    -1.92560    0.02095  -91.90   <2e-16 ***
#> Xd1          2.29771    0.02498   92.00   <2e-16 ***
#> Xd2         -1.01303    0.02077  -48.78   <2e-16 ***
#> X1           2.59564    0.03410   76.11   <2e-16 ***
#> X2          -1.07589    0.03494  -30.79   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 1.39 on 4994 degrees of freedom
#> Multiple R-squared:  0.7722, Adjusted R-squared:  0.772 
#> F-statistic:  3386 on 5 and 4994 DF,  p-value: < 2.2e-16
#> 
#> 
#> Supply Equation:
#> 
#> Call:
#> lm(formula = supply_formula, data = object@data)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -7.3439 -0.8408  0.1576  0.9897  3.8530 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) 12.00283    0.11136  107.78   <2e-16 ***
#> P_FITTED     1.20302    0.02439   49.33   <2e-16 ***
#> Xs1          2.91737    0.03006   97.04   <2e-16 ***
#> X1          -1.51101    0.03777  -40.00   <2e-16 ***
#> X2           3.20439    0.03895   82.27   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 1.39 on 4995 degrees of freedom
#> Multiple R-squared:  0.772,  Adjusted R-squared:  0.7718 
#> F-statistic:  4227 on 4 and 4995 DF,  p-value: < 2.2e-16
summary(eq_fit)
#> Equilibrium Model for Markets in Equilibrium:
#>   Demand RHS        :   D_P + D_Xd1 + D_Xd2 + D_X1 + D_X2
#>   Supply RHS        :   S_P + S_Xs1 + S_X1 + S_X2
#>   Market Clearing   : Q = D_Q = S_Q
#>   Shocks            : Correlated
#>   Nobs              : 5000
#>   Sample Separation : Not Separated
#>   Quantity Var      : Q
#>   Price Var         : P
#>   Key Var(s)        : id, date
#>   Time Var          : date
#> 
#> Maximum likelihood estimation:
#>   Method              : BFGS
#>   Convergence Status  : success
#>   Starting Values     :
#>        D_P    D_CONST      D_Xd1      D_Xd2       D_X1       D_X2        S_P 
#>    -1.9256    26.0536     2.2977    -1.0130     2.5956    -1.0759     1.2030 
#>    S_CONST      S_Xs1       S_X1       S_X2 D_VARIANCE S_VARIANCE        RHO 
#>    12.0028     2.9174    -1.5110     3.2044     5.6643     5.7012    -0.3968 
#> 
#> Coefficients:
#>              Estimate Std. Error   z value      Pr(>|z|) 
#>  D_P        -1.9256975 0.03725261 -51.69296  0.000000e+00 ***
#>  D_CONST    26.0542576 0.16556392 157.36676  0.000000e+00 ***
#>  D_Xd1       2.3055160 0.04361180  52.86451  0.000000e+00 ***
#>  D_Xd2      -0.9963758 0.03348349 -29.75723 1.398342e-194 ***
#>  D_X1        2.5956115 0.05978398  43.41651  0.000000e+00 ***
#>  D_X2       -1.0762772 0.06019306 -17.88042  1.675703e-71 ***
#>  S_P         1.2032908 0.04260645  28.24199 1.785188e-175 ***
#>  S_CONST    12.0015431 0.19885727  60.35255  0.000000e+00 ***
#>  S_Xs1       2.9177581 0.05258619  55.48525  0.000000e+00 ***
#>  S_X1       -1.5113726 0.06625113 -22.81278 3.423737e-115 ***
#>  S_X2        3.2047666 0.06723523  47.66499  0.000000e+00 ***
#>  D_VARIANCE  5.6639426 0.22577277  25.08692 6.909728e-139 ***
#>  S_VARIANCE  5.7012019 0.25913008  22.00131 2.797632e-107 ***
#>  RHO        -0.3970413 0.01730149 -22.94838 1.529232e-116 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> -2 log L: 33488
summary(bs_fit)
#> Basic Model for Markets in Disequilibrium:
#>   Demand RHS        :   D_P + D_Xd1 + D_Xd2 + D_X1 + D_X2
#>   Supply RHS        :   S_P + S_Xs1 + S_X1 + S_X2
#>   Short Side Rule   : Q = min(D_Q, S_Q)
#>   Shocks            : Correlated
#>   Nobs              : 5000
#>   Sample Separation : Not Separated
#>   Quantity Var      : Q
#>   Price Var         : P
#>   Key Var(s)        : id, date
#>   Time Var          : date
#> 
#> Maximum likelihood estimation:
#>   Method              : Nelder-Mead
#>   Max Iterations      : 1e+05
#>   Convergence Status  : success
#>   Starting Values     :
#>        D_P    D_CONST      D_Xd1      D_Xd2       D_X1       D_X2        S_P 
#>   -0.91370   21.50697    1.55130   -0.68159    1.24540    0.30929    0.08391 
#>    S_CONST      S_Xs1       S_X1       S_X2 D_VARIANCE S_VARIANCE        RHO 
#>   17.03305    1.86928   -0.02750    1.66942    3.10289    2.85599    0.00000 
#> 
#> Coefficients:
#>              Estimate Std. Error    z value     Pr(>|z|) 
#>  D_P        -0.7791371 0.07931915  -9.822812 8.980303e-23 ***
#>  D_CONST    22.6344592 0.29751814  76.077578 0.000000e+00 ***
#>  D_Xd1       2.0117793 0.03523165  57.101478 0.000000e+00 ***
#>  D_Xd2      -0.6992671 0.04193071 -16.676727 1.935354e-62 ***
#>  D_X1        1.7158904 0.14576979  11.771234 5.491502e-32 ***
#>  D_X2       -0.5892051 0.07792261  -7.561413 3.987138e-14 ***
#>  S_P         0.3290079 0.04549028   7.232488 4.742254e-13 ***
#>  S_CONST    17.5777710 0.23544440  74.657843 0.000000e+00 ***
#>  S_Xs1       2.7156153 0.03555906  76.369148 0.000000e+00 ***
#>  S_X1       -0.9950805 0.09697323 -10.261394 1.051773e-24 ***
#>  S_X2        2.8737782 0.06496191  44.237895 0.000000e+00 ***
#>  D_VARIANCE  2.5069085 0.57009934   4.397319 1.095963e-05 ***
#>  S_VARIANCE  1.7364128 0.22626449   7.674261 1.663745e-14 ***
#>  RHO        -0.6685024 0.06631271 -10.081061 6.699888e-24 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> -2 log L: 14692.09
summary(da_fit)
#> Deterministic Adjustment Model for Markets in Disequilibrium:
#>   Demand RHS        :   D_P + D_Xd1 + D_Xd2 + D_X1 + D_X2
#>   Supply RHS        :   S_P + S_Xs1 + S_X1 + S_X2
#>   Short Side Rule   : Q = min(D_Q, S_Q)
#>   Separation Rule   : P_DIFF analogous to (D_Q - S_Q)
#>   Shocks            : Correlated
#>   Nobs              : 4000
#>   Sample Separation : Demand Obs = 1758, Supply Obs = 2242
#>   Quantity Var      : Q
#>   Price Var         : P
#>   Key Var(s)        : id, date
#>   Time Var          : date
#> 
#> Maximum likelihood estimation:
#>   Method              : BFGS
#>   Convergence Status  : success
#>   Starting Values     :
#>        D_P    D_CONST      D_Xd1      D_Xd2       D_X1       D_X2        S_P 
#>   -1.00158   22.19624    1.67280   -0.73584    1.48376    0.04390    0.10972 
#>    S_CONST      S_Xs1       S_X1       S_X2     P_DIFF D_VARIANCE S_VARIANCE 
#>   16.83321    1.80322    0.05562    1.57367    0.92085    2.58729    2.98552 
#>        RHO 
#>    0.00000 
#> 
#> Coefficients:
#>              Estimate Std. Error    z value      Pr(>|z|) 
#>  D_P        -1.3431249 0.01644569 -81.670328  0.000000e+00 ***
#>  D_CONST    25.3875755 0.08826183 287.639344  0.000000e+00 ***
#>  D_Xd1       2.3709604 0.02657810  89.207292  0.000000e+00 ***
#>  D_Xd2      -1.0175771 0.02196629 -46.324492  0.000000e+00 ***
#>  D_X1        2.7721752 0.03298164  84.052075  0.000000e+00 ***
#>  D_X2       -1.3064350 0.03369226 -38.775526  0.000000e+00 ***
#>  S_P         0.5271217 0.01507390  34.969158 6.624101e-268 ***
#>  S_CONST    15.8320910 0.07487340 211.451463  0.000000e+00 ***
#>  S_Xs1       2.6868156 0.02499791 107.481594  0.000000e+00 ***
#>  S_X1       -1.1388234 0.02977150 -38.252129  0.000000e+00 ***
#>  S_X2        2.8187413 0.02949716  95.559769  0.000000e+00 ***
#>  P_DIFF      1.1021707 0.01308521  84.230261  0.000000e+00 ***
#>  D_VARIANCE  1.7590874 0.04799147  36.654166 3.927494e-294 ***
#>  S_VARIANCE  1.5970307 0.04168354  38.313223  0.000000e+00 ***
#>  RHO         0.0600957 0.01783407   3.369713  7.524658e-04 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> -2 log L: 18105.45
summary(sa_fit)
#> Stochastic Adjustment Model for Markets in Disequilibrium:
#>   Demand RHS        :   D_P + D_Xd1 + D_Xd2 + D_X1 + D_X2
#>   Supply RHS        :   S_P + S_Xs1 + S_X1 + S_X2
#>   Price Dynamics RHS:   I(D_Q - S_Q) + Xp1
#>   Short Side Rule   : Q = min(D_Q, S_Q)
#>   Shocks            : Correlated
#>   Nobs              : 4000
#>   Sample Separation : Not Separated
#>   Quantity Var      : Q
#>   Price Var         : P
#>   Key Var(s)        : id, date
#>   Time Var          : date
#> 
#> Maximum likelihood estimation:
#>   Method              : BFGS
#>   Max Iterations      : 1e+05
#>   Convergence Status  : success
#>   Starting Values     :
#>        D_P    D_CONST      D_Xd1      D_Xd2       D_X1       D_X2        S_P 
#>   -1.00158   22.19624    1.67280   -0.73584    1.48376    0.04390    0.10972 
#>    S_CONST      S_Xs1       S_X1       S_X2     P_DIFF    P_CONST      P_Xp1 
#>   16.83321    1.80322    0.05562    1.57367    0.92252    0.42843   -0.07090 
#> D_VARIANCE S_VARIANCE P_VARIANCE     RHO_DS     RHO_DP     RHO_SP 
#>    2.58729    2.98552    5.60450    0.00000    0.00000    0.00000 
#> 
#> Coefficients:
#>               Estimate Std. Error     z value      Pr(>|z|) 
#>  D_P        -1.31678069 0.01374289 -95.8154072  0.000000e+00 ***
#>  D_CONST    24.88425787 0.07604569 327.2277235  0.000000e+00 ***
#>  D_Xd1       2.29596931 0.02228707 103.0179994  0.000000e+00 ***
#>  D_Xd2      -0.98502218 0.01939144 -50.7967401  0.000000e+00 ***
#>  D_X1        2.62337313 0.02912280  90.0796982  0.000000e+00 ***
#>  D_X2       -1.14372398 0.03011240 -37.9818330  0.000000e+00 ***
#>  S_P         0.63216311 0.01611099  39.2380068  0.000000e+00 ***
#>  S_CONST    15.98170066 0.07430136 215.0929829  0.000000e+00 ***
#>  S_Xs1       2.97894891 0.02786543 106.9048104  0.000000e+00 ***
#>  S_X1       -1.55155068 0.03410057 -45.4992563  0.000000e+00 ***
#>  S_X2        3.24779056 0.03521129  92.2371880  0.000000e+00 ***
#>  P_DIFF      1.21068401 0.01331229  90.9448118  0.000000e+00 ***
#>  P_CONST     0.89786543 0.03672925  24.4455156 5.616513e-132 ***
#>  P_Xp1      -0.10822251 0.02145413  -5.0443670  4.550253e-07 ***
#>  D_VARIANCE  1.04693421 0.03435125  30.4773222 5.206602e-204 ***
#>  S_VARIANCE  0.98237039 0.03766173  26.0840478 5.531219e-150 ***
#>  P_VARIANCE  0.99622390 0.09802695  10.1627555  2.907411e-24 ***
#>  RHO_DS     -0.02631395 0.06405949  -0.4107736  6.812385e-01  
#>  RHO_DP     -0.04298693 0.05583361  -0.7699113  4.413525e-01  
#>  RHO_SP     -0.01373388 0.06060979  -0.2265951  8.207386e-01  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> -2 log L: 16932.03

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_.

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
#>                bs         dr         da          sa
#> B_P   -0.10160382 -0.1366526 -0.1919618 -0.19683597
#> D_Xd1  0.18445642  0.2518488  0.2433550  0.23188424
#> D_Xd2 -0.06411454 -0.1076561 -0.1044440 -0.09948352
#> B_X1   0.24856404  0.3359714  0.4014243  0.42165155
#> B_X2  -0.31751469 -0.3534875 -0.4234075 -0.44352639
#> S_Xs1 -0.24898988 -0.2633837 -0.2757743 -0.30086260

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
#>               bs         dr         da         sa
#> B_P   -0.1666732 -0.3001412 -0.4054235 -0.4984543
#> D_Xd1  0.3025865  0.5531559  0.5139659  0.5872082
#> D_Xd2 -0.1051750 -0.2364539 -0.2205857 -0.2519255
#> B_X1   0.4077501  0.7379211  0.8478083  1.0677624
#> B_X2  -0.5208583 -0.7763931 -0.8942368 -1.1231568
#> S_Xs1 -0.4084487 -0.5784910 -0.5824356 -0.7618845

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.

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.

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))])
}
#> Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
#> 3.5.0.
#> ℹ Please use the `legend.position.inside` argument of `theme()` instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.

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.   : 8.45   Min.   : 4.754  
#>  1st Qu.:16.69   1st Qu.:16.746  
#>  Median :18.57   Median :19.093  
#>  Mean   :18.55   Mean   :19.115  
#>  3rd Qu.:20.36   3rd Qu.:21.454  
#>  Max.   :29.74   Max.   :30.763

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