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.
Load the required libraries.
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
)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 + X2Set the estimation parameters.
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.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: 0x557cd48e7d68>
#> Instruments: ~Xd1 + Xd2 + X1 + X2 + Xs1
#> <environment: 0x557cd48e7d68>
#>
#> 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: 0x557cd48e7d68>
#> Instruments: ~Xd1 + Xd2 + X1 + X2 + Xs1
#> <environment: 0x557cd48e7d68>
#>
#> 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.67The 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.0192Since 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.03684207Moreover, 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