# An example of market-clearing assessment

This short tutorial gives an example of statistically assessing whether a market is in an equilibrium state. The tutorial assumes some familiarity with the concepts and the functionality of the package. The basic_usage vignette can be helpful in acquiring this familiarity.

## Setup the environment

Load the required libraries.

    library(diseq)
    library(magrittr)

Prepare the data. Here, we simply simulate data using a data generating process for a market in equilibrium.

    nobs <- 1000
    tobs <- 5

    alpha_d <- -3.9
    beta_d0 <- 28.9
    beta_d <- c(2.1, -0.7)
    eta_d <- c(3.5, 6.25)

    alpha_s <- 2.8
    beta_s0 <- 26.2
    beta_s <- c(2.65)
    eta_s <- c(1.15, 4.2)

    sigma_d <- 0.8
    sigma_s <- 1.1
    rho_ds <- 0.0

    seed <- 42

    eq_data <- simulate_data(
      "equilibrium_model", nobs, tobs,
      alpha_d, beta_d0, beta_d, eta_d,
      alpha_s, beta_s0, beta_s, eta_s,
      NA, NA, c(NA),
      sigma_d = sigma_d, sigma_s = sigma_s, rho_ds = rho_ds,
      seed = seed
    )

## Estimate the models

Prepare the basic parameters for model initialization.

    verbose <- 2
    correlated_shocks <- TRUE
    formula <-   Q | P | id | date ~ P + Xd1 + Xd2 + X1 + X2 | P + Xs1 + X1 + X2

Set the estimation parameters.

    optimization_method <- "BFGS"
    optimization_options <- list(maxit = 10000, reltol = 1e-8)

Using the above parameterization, construct and estimate the model objects. Here we estimate two equilibrium models and four disequilibrium models. All the models are constructed using the simulated data from a model of market in equilibrium.

    eq_reg <- equilibrium_model(
      formula, eq_data[eq_data$date != 1, ],
      correlated_shocks = correlated_shocks, verbose = verbose,
      estimation_options = list(method = "2SLS")
    )
    #> Info: This is 'Equilibrium with correlated shocks' model
    #> Warning: Removing unobserved '1' level(s).
    eq_fit <- equilibrium_model(
      formula, eq_data[eq_data$date != 1, ],
      correlated_shocks = correlated_shocks, verbose = verbose,
      estimation_options = list(
        control = optimization_options, method = optimization_method
      )
    )
    #> Info: This is 'Equilibrium with correlated shocks' model
    #> Warning: Removing unobserved '1' level(s).
    bs_fit <- diseq_basic(
      formula, eq_data[eq_data$date != 1, ],
      correlated_shocks = correlated_shocks, verbose = verbose,
      estimation_options = list(
        control = optimization_options, method = optimization_method
      )
    )
    #> Info: This is 'Basic with correlated shocks' model
    #> Warning: Removing unobserved '1' level(s).
    da_fit <- diseq_deterministic_adjustment(
      formula, eq_data,
      correlated_shocks = correlated_shocks, verbose = verbose,
      estimation_options = list(
        control = optimization_options, method = optimization_method
      )
    )
    #> Info: This is 'Deterministic Adjustment with correlated shocks' model
    #> Info: Dropping 1000 rows to generate 'LAGGED_P'.
    #> Info: Sample separated with 2000 rows in excess supply and 2000 in excess demand state.

## Post estimation analysis

### Summaries

All the models provide estimates for the simulated data. Even with simulated data, it is difficult to assess which model performs better by examining only the summaries in separation or collectively.

    summary(eq_reg@fit[[1]]$first_stage_model)
    #>
    #> Call:
    #> lm(formula = first_stage_formula, data = object@model_tibble)
    #>
    #> Residuals:
    #>      Min       1Q   Median       3Q      Max
    #> -0.88388 -0.14215  0.00305  0.14088  0.81967
    #>
    #> Coefficients:
    #>              Estimate Std. Error t value Pr(>|t|)
    #> (Intercept)  0.402566   0.003302  121.91   <2e-16 ***
    #> Xd1          0.313982   0.003288   95.50   <2e-16 ***
    #> Xd2         -0.103262   0.003274  -31.54   <2e-16 ***
    #> X1           0.348932   0.003262  106.98   <2e-16 ***
    #> X2           0.306357   0.003222   95.07   <2e-16 ***
    #> Xs1         -0.398531   0.003284 -121.36   <2e-16 ***
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #>
    #> Residual standard error: 0.2087 on 3994 degrees of freedom
    #> Multiple R-squared:  0.9194, Adjusted R-squared:  0.9193
    #> F-statistic:  9110 on 5 and 3994 DF,  p-value: < 2.2e-16
    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              : 4000
    #>   Sample Separation : Not Separated
    #>   Quantity Var      : Q
    #>   Price Var         : P
    #>   Key Var(s)        : id, date
    #>   Time Var          : date
    #>
    #> systemfit results
    #> method: 2SLS
    #>
    #>           N   DF     SSR  detRCov   OLS-R2 McElroy-R2
    #> system 8000 7989 7718.08 0.833283 0.972727   0.976083
    #>
    #>           N   DF     SSR      MSE     RMSE       R2   Adj R2
    #> demand 4000 3994 2597.52 0.650356 0.806447 0.981643 0.981620
    #> supply 4000 3995 5120.56 1.281742 1.132141 0.963812 0.963776
    #>
    #> The covariance matrix of the residuals
    #>            demand     supply
    #> demand  0.6503560 -0.0174889
    #> supply -0.0174889  1.2817422
    #>
    #> The correlations of the residuals
    #>            demand     supply
    #> demand  1.0000000 -0.0191552
    #> supply -0.0191552  1.0000000
    #>
    #>
    #> 2SLS estimates for 'demand' (equation 1)
    #> Model Formula: Q ~ P + Xd1 + Xd2 + X1 + X2
    #> <environment: 0x55f54e57b118>
    #> Instruments: ~Xd1 + Xd2 + X1 + X2 + Xs1
    #> <environment: 0x55f54e57b118>
    #>
    #>               Estimate Std. Error  t value   Pr(>|t|)
    #> (Intercept) 28.8908936  0.0182726 1581.108 < 2.22e-16 ***
    #> P           -3.8765576  0.0318337 -121.775 < 2.22e-16 ***
    #> Xd1          2.1050973  0.0162669  129.410 < 2.22e-16 ***
    #> Xd2         -0.7077763  0.0131230  -53.934 < 2.22e-16 ***
    #> X1           3.5074743  0.0167043  209.975 < 2.22e-16 ***
    #> X2           6.2465504  0.0158419  394.306 < 2.22e-16 ***
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #>
    #> Residual standard error: 0.806447 on 3994 degrees of freedom
    #> Number of observations: 4000 Degrees of Freedom: 3994
    #> SSR: 2597.521998 MSE: 0.650356 Root MSE: 0.806447
    #> Multiple R-Squared: 0.981643 Adjusted R-Squared: 0.98162
    #>
    #>
    #> 2SLS estimates for 'supply' (equation 2)
    #> Model Formula: Q ~ P + Xs1 + X1 + X2
    #> <environment: 0x55f54e57b118>
    #> Instruments: ~Xd1 + Xd2 + X1 + X2 + Xs1
    #> <environment: 0x55f54e57b118>
    #>
    #>               Estimate Std. Error  t value   Pr(>|t|)
    #> (Intercept) 26.1862350  0.0279220 937.8368 < 2.22e-16 ***
    #> P            2.8426868  0.0539250  52.7156 < 2.22e-16 ***
    #> Xs1          2.6776568  0.0281531  95.1107 < 2.22e-16 ***
    #> X1           1.1627870  0.0256741  45.2902 < 2.22e-16 ***
    #> X2           4.1882859  0.0242913 172.4191 < 2.22e-16 ***
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #>
    #> Residual standard error: 1.132141 on 3995 degrees of freedom
    #> Number of observations: 4000 Degrees of Freedom: 3995
    #> SSR: 5120.560239 MSE: 1.281742 Root MSE: 1.132141
    #> Multiple R-Squared: 0.963812 Adjusted R-Squared: 0.963776
    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              : 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      : 10000
    #>   Relative Tolerance  : 1e-08
    #>   Convergence Status  : success
    #>   Starting Values     :
    #>        D_P    D_CONST      D_Xd1      D_Xd2       D_X1       D_X2        S_P
    #>   -3.87656   28.89089    2.10510   -0.70778    3.50747    6.24655    2.84269
    #>    S_CONST      S_Xs1       S_X1       S_X2 D_VARIANCE S_VARIANCE        RHO
    #>   26.18624    2.67766    1.16279    4.18829    0.65036    1.28174   -0.01916
    #>
    #> Coefficients
    #>            Estimate Std. Error  z value      Pr(z)
    #> D_P         -3.8766    0.03181 -121.857  0.000e+00
    #> D_CONST     28.8909    0.01826 1582.297  0.000e+00
    #> D_Xd1        2.1050    0.01624  129.610  0.000e+00
    #> D_Xd2       -0.7079    0.01313  -53.918  0.000e+00
    #> D_X1         3.5075    0.01669  210.118  0.000e+00
    #> D_X2         6.2465    0.01583  394.602  0.000e+00
    #> S_P          2.8430    0.05390   52.750  0.000e+00
    #> S_CONST     26.1861    0.02791  938.343  0.000e+00
    #> S_Xs1        2.6778    0.02814   95.168  0.000e+00
    #> S_X1         1.1627    0.02566   45.310  0.000e+00
    #> S_X2         4.1882    0.02428  172.510  0.000e+00
    #> D_VARIANCE   0.6494    0.01584   41.010  0.000e+00
    #> S_VARIANCE   1.2803    0.03540   36.169 1.881e-286
    #> RHO         -0.0192    0.01807   -1.062  2.881e-01
    #>
    #> -2 log L: 6722.67
    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              : 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      : 10000
    #>   Relative Tolerance  : 1e-08
    #>   Convergence Status  : success
    #>   Starting Values     :
    #>        D_P    D_CONST      D_Xd1      D_Xd2       D_X1       D_X2        S_P
    #>    -3.3899    28.6909     1.9497    -0.6543     3.3398     6.0968     1.5864
    #>    S_CONST      S_Xs1       S_X1       S_X2 D_VARIANCE S_VARIANCE        RHO
    #>    26.6854     2.1697     1.5963     4.5814     0.6012     1.0378     0.0000
    #>
    #> Coefficients
    #>            Estimate Std. Error z value      Pr(z)
    #> D_P         -3.5088    0.04544  -77.22  0.000e+00
    #> D_CONST     29.1229    0.03412  853.59  0.000e+00
    #> D_Xd1        2.0719    0.02544   81.44  0.000e+00
    #> D_Xd2       -0.7088    0.01662  -42.65  0.000e+00
    #> D_X1         3.3653    0.02546  132.20  0.000e+00
    #> D_X2         6.1435    0.02472  248.49  0.000e+00
    #> S_P          2.0706    0.10797   19.18  5.585e-82
    #> S_CONST     27.8542    0.05857  475.60  0.000e+00
    #> S_Xs1        2.5828    0.06431   40.16  0.000e+00
    #> S_X1         1.4476    0.05402   26.80 3.340e-158
    #> S_X2         4.3986    0.05373   81.86  0.000e+00
    #> D_VARIANCE   0.6386    0.02336   27.34 1.359e-164
    #> S_VARIANCE   1.1249    0.07243   15.53  2.131e-54
    #> RHO         -0.4709    0.04036  -11.67  1.872e-31
    #>
    #> -2 log L: 8338.91
    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 = 2000, Supply Obs = 2000
    #>   Quantity Var      : Q
    #>   Price Var         : P
    #>   Key Var(s)        : id, date
    #>   Time Var          : date
    #>
    #> Maximum likelihood estimation
    #>   Method              : BFGS
    #>   Max Iterations      : 10000
    #>   Relative Tolerance  : 1e-08
    #>   Convergence Status  : success
    #>   Starting Values     :
    #>        D_P    D_CONST      D_Xd1      D_Xd2       D_X1       D_X2        S_P
    #> -3.3899420 28.6908902  1.9497497 -0.6543147  3.3398395  6.0967924  1.5864480
    #>    S_CONST      S_Xs1       S_X1       S_X2     P_DIFF D_VARIANCE S_VARIANCE
    #> 26.6853581  2.1696502  1.5962938  4.5814415 -0.0007971  0.6012371  1.0378385
    #>        RHO
    #>  0.0000000
    #>
    #> Coefficients
    #>             Estimate Std. Error    z value      Pr(z)
    #> D_P        -3.875716    0.03396 -114.13950  0.000e+00
    #> D_CONST    28.891196    0.01875 1540.93773  0.000e+00
    #> D_Xd1       2.105016    0.01624  129.59540  0.000e+00
    #> D_Xd2      -0.707967    0.01313  -53.92060  0.000e+00
    #> D_X1        3.507454    0.01670  210.07124  0.000e+00
    #> D_X2        6.246540    0.01583  394.60023  0.000e+00
    #> S_P         2.842225    0.05490   51.76647  0.000e+00
    #> S_CONST    26.187074    0.03095  846.01070  0.000e+00
    #> S_Xs1       2.677787    0.02814   95.16504  0.000e+00
    #> S_X1        1.162670    0.02566   45.30858  0.000e+00
    #> S_X2        4.188175    0.02428  172.49527  0.000e+00
    #> P_DIFF      0.001553    0.02176    0.07137  9.431e-01
    #> D_VARIANCE  0.649373    0.01584   41.00771  0.000e+00
    #> S_VARIANCE  1.280272    0.03540   36.16826 1.922e-286
    #> RHO        -0.019191    0.01807   -1.06181  2.883e-01
    #>
    #> -2 log L: 6722.67

### Model selection

The deterministic adjustment model has price dynamics that are analogous to excess demand and estimates one extra parameter. The directional model estimates one parameter less as the model does not have enough equations to identify prices in both demand and supply equations. The estimated parameters are summarized as follows.

    da_coef <- coef(da_fit)
    coef_names <- names(da_coef)

    sim_coef <- c(
      alpha_d, beta_d0, beta_d, eta_d,
      alpha_s, beta_s0, beta_s, eta_s,
      NA,
      sigma_d, sigma_s,
      rho_ds
    )

    coef_tbl <- function(fit) {
        tibble::tibble(coef = names(coef(fit)),!!substitute(fit) := coef(fit))
    }

    comp <- coef_tbl(da_fit) %>%
        dplyr::left_join(coef_tbl(bs_fit), by = "coef") %>%
        dplyr::left_join(coef_tbl(eq_reg), by = "coef") %>%
        dplyr::left_join(coef_tbl(eq_fit), by = "coef") %>%
        dplyr::mutate(sim = sim_coef) %>%
        dplyr::mutate(sim = sim_coef,
                      da_fit_err = abs(da_fit - sim),
                      bs_fit_err = abs(bs_fit - sim),
                      eq_fit_err = abs(eq_fit - sim))

    comp
    #> # A tibble: 15 × 9
    #>    coef     da_fit bs_fit  eq_reg  eq_fit   sim da_fit_err bs_fit_err eq_fit_err
    #>    <chr>     <dbl>  <dbl>   <dbl>   <dbl> <dbl>      <dbl>      <dbl>      <dbl>
    #>  1 D_P    -3.88    -3.51  -3.88   -3.88   -3.9     0.0243     0.391      0.0234
    #>  2 D_CON… 28.9     29.1   28.9    28.9    28.9     0.00880    0.223      0.00910
    #>  3 D_Xd1   2.11     2.07   2.11    2.11    2.1     0.00502    0.0281     0.00505
    #>  4 D_Xd2  -0.708   -0.709 -0.708  -0.708  -0.7     0.00797    0.00883    0.00795
    #>  5 D_X1    3.51     3.37   3.51    3.51    3.5     0.00745    0.135      0.00749
    #>  6 D_X2    6.25     6.14   6.25    6.25    6.25    0.00346    0.106      0.00345
    #>  7 S_P     2.84     2.07   2.84    2.84    2.8     0.0422     0.729      0.0430
    #>  8 S_CON… 26.2     27.9   26.2    26.2    26.2     0.0129     1.65       0.0139
    #>  9 S_Xs1   2.68     2.58   2.68    2.68    2.65    0.0278     0.0672     0.0278
    #> 10 S_X1    1.16     1.45   1.16    1.16    1.15    0.0127     0.298      0.0127
    #> 11 S_X2    4.19     4.40   4.19    4.19    4.2     0.0118     0.199      0.0118
    #> 12 P_DIFF  0.00155 NA     NA      NA      NA      NA         NA         NA
    #> 13 D_VAR…  0.649    0.639  0.650   0.649   0.8     0.151      0.161      0.151
    #> 14 S_VAR…  1.28     1.12   1.28    1.28    1.1     0.180      0.0249     0.180
    #> 15 RHO    -0.0192  -0.471 -0.0192 -0.0192  0       0.0192     0.471      0.0192

Since we have used simulated data, we can calculate the average absolute error of the parameter estimation for each of the models. The population values are unknown in practice, and this calculation is impossible.

    model_errors <- colMeans(comp[, grep("err", colnames(comp))], na.rm = TRUE)
    model_errors
    #> da_fit_err bs_fit_err eq_fit_err
    #> 0.03675054 0.32116445 0.03684207

Moreover, the average absolute error cannot provide an overall estimation assessment as the market models have different parameter spaces. To assess the overall model performance, one can instead use an information criterion.

    fits <- c(da_fit, bs_fit, eq_fit)
    model_names <- sapply(fits, function(m) m@model_type_string)
    model_obs <- sapply(fits, nobs)
    aic <- sapply(fits, AIC)
    df <- sapply(fits, function(m) attr(logLik(m), "df"))
    seltbl <- tibble::tibble(Model = model_names, AIC = aic,
                             D.F = df, Obs. = model_obs,
                             Abs.Error = model_errors) %>%
      dplyr::arrange(AIC)
    seltbl
    #> # A tibble: 3 × 5
    #>   Model                      AIC   D.F  Obs. Abs.Error
    #>   <chr>                    <dbl> <int> <int>     <dbl>
    #> 1 Equilibrium              6751.    14  4000    0.0368
    #> 2 Deterministic Adjustment 6753.    15  4000    0.0368
    #> 3 Basic                    8367.    14  4000    0.321