This tutorial examines model initialization and estimation in some detail. Models can be initialized and estimated with a single function call (see basic_usage), which is the recommended approach for most usage cases. However, it is convenient to separate model estimation and initialization on some occasions. This is particularly relevant when estimating the same model using different methods and/or options without re-initializing.
The operations of this vignette cover the many but not all use initialization cases. More usage details can be found in the documentation of the package.
Load the required libraries.
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 <- 2000
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 <- 443
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
)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:
The fields that uniquely identify simulated data records are
id (for subjects) and date (for time). These
variable names are automatically set for the data that
simulate_data generates.
The quantity variable is automatically named Q by
the simulate_data function. The quantity variable is
observable. For the equilibrium models, it is equal to both the demanded
and supplied quantities. The observed quantity represents either a
demanded or a supplied quantity for the disequilibrium models. Each
disequilibrium model resolves that state of the observation in a
different way.
The price variable is set as P by the
simulate_data call.
The right-hand sides of the demand and supply equations. 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 also requires the specification
of price dynamics. The simulate_data call generates 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.
Set the verbosity level. This controls the level of messaging. The verbosity level here is set so that the constructed objects display basic information in addition to errors and warnings.
Using the above parameterization, construct the model objects. Here,
we construct an equilibrium model 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 subject.
eq <- new(
"equilibrium_model",
quantity = Q, price = P,
demand = P + Xd1 + Xd2 + X1 + X2,
supply = P + Xs1 + X1 + X2,
subject = id, time = date,
data = stochastic_adjustment_data,
correlated_shocks = correlated_shocks, verbose = verbose
)
#> Info: This is Equilibrium model.
bs <- new(
"diseq_basic",
quantity = Q, price = P,
demand = P + Xd1 + Xd2 + X1 + X2,
supply = P + Xs1 + X1 + X2,
subject = id, time = date,
data = stochastic_adjustment_data,
correlated_shocks = correlated_shocks, verbose = verbose
)
#> Info: This is Basic model.
dr <- new(
"diseq_directional",
quantity = Q, price = P,
demand = P + Xd1 + Xd2 + X1 + X2,
supply = Xs1 + X1 + X2,
subject = id, time = date,
data = stochastic_adjustment_data,
correlated_shocks = correlated_shocks, verbose = verbose
)
#> Info: This is Directional model.
#> Info: Dropping 2000 rows to generate 'LAGGED_P'.
#> Info: Sample separated with 3473 rows in excess supply and 4527 rows in
#> excess demand states.
da <- new(
"diseq_deterministic_adjustment",
quantity = Q, price = P,
demand = P + Xd1 + Xd2 + X1 + X2,
supply = P + Xs1 + X1 + X2,
subject = id, time = date,
data = stochastic_adjustment_data,
correlated_shocks = correlated_shocks, verbose = verbose
)
#> Info: This is Deterministic Adjustment model.
#> Info: Dropping 2000 rows to generate 'LAGGED_P'.
#> Info: Sample separated with 3473 rows in excess supply and 4527 rows in
#> excess demand states.
sa <- new(
"diseq_stochastic_adjustment",
quantity = Q, price = P,
demand = P + Xd1 + Xd2 + X1 + X2,
supply = P + Xs1 + X1 + X2,
price_dynamics = Xp1,
subject = id, time = date,
data = stochastic_adjustment_data,
correlated_shocks = correlated_shocks, verbose = verbose
)
#> Info: This is Stochastic Adjustment model.
#> Info: Dropping 2000 rows to generate 'LAGGED_P'.First, we need to set the estimation parameters and choose an
estimation method. The only model that can be estimated by least squares
is the equilibrium_model. To estimate the model with this
methodology call markets::estimate with
method = 2SLS set. The equilibrium_model can
also be estimated using full information maximum likelihood, as it is
the case for all the disequilibrium models. One may 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 markets package uses this method.
"L-BFGS-B": Constrained optimization.
optimization_method <- "BFGS"
optimization_options <- list(REPORT = 10, maxit = 10000, reltol = 1e-6)Then, estimate the models. The eq model is estimated
with two different methods, namely two stage least squares and full
information maximum likelihood. Moreover, the bs is
estimated using two different optimization options; these are the
gradient-based "BFGS" method and the simplex-based
"Nelder-Mead" methods. Lastly, the models estimated with
maximal likelihood use different estimation options regarding the
calculation of standard errors. See the documentation for more
options.
estimate(eq, method = "2SLS")
#> 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
#>
#> Least squares estimation:
#> Method : 2SLS
estimate(eq,
control = optimization_options, method = optimization_method,
standard_errors = c("id")
)
#> 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
#>
#> Maximum likelihood estimation:
#> Method : BFGS
#> Convergence Status : success
estimate(bs,
control = optimization_options, method = optimization_method,
standard_errors = "heteroscedastic"
)
#> 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
#>
#> Maximum likelihood estimation:
#> Method : BFGS
#> Convergence Status : success
estimate(bs,
control = optimization_options, method = "Nelder-Mead",
standard_errors = "heteroscedastic"
)
#> 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
#>
#> Maximum likelihood estimation:
#> Method : Nelder-Mead
#> Convergence Status : success
estimate(dr,
control = optimization_options, method = optimization_method,
standard_errors = "heteroscedastic"
)
#> Directional Model for Markets in Disequilibrium:
#> Demand RHS : D_P + D_Xd1 + D_Xd2 + D_X1 + D_X2
#> Supply RHS : S_Xs1 + S_X1 + S_X2
#> Short Side Rule : Q = min(D_Q, S_Q)
#> Separation Rule : P_DIFF >= 0 then D_Q >= S_Q
#> Shocks : Correlated
#>
#> Maximum likelihood estimation:
#> Method : BFGS
#> Convergence Status : success
estimate(da,
control = optimization_options, method = optimization_method,
standard_errors = c("id")
)
#> 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
#>
#> Maximum likelihood estimation:
#> Method : BFGS
#> Convergence Status : success
estimate(sa,
control = optimization_options, method = optimization_method
)
#> 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
#>
#> Maximum likelihood estimation:
#> Method : BFGS
#> Convergence Status : success