msel: Multivariate and multinomial sample selection and endogenous...

In switchSelection: Endogenous Switching and Sample Selection Regression Models

Multivariate and multinomial sample selection and endogenous switching models with multiple outcomes.

Description

This function allows to estimate parameters of the multivariate and multinomial sample selection and endogenous switching models with multiple outcomes. Both maximum-likelihood and two-step estimators are implemented.

Usage

msel(
  formula = NA,
  formula2 = NA,
  formula3 = NA,
  data = NULL,
  groups = NA,
  groups2 = NA,
  groups3 = NA,
  marginal = list(),
  opt_type = "optim",
  opt_args = NA,
  start = NULL,
  estimator = "ml",
  cov_type = "mm",
  degrees = NA,
  degrees3 = NA,
  n_sim = 1000,
  n_cores = 1,
  control = list(),
  regularization = list(),
  type3 = "logit"
)

Arguments

formula	a list which i-th element is an object of class "formula" describing the form of the linear predictor (index) of the i-th ordinal equation. Mean and variance equations should be separated by the '|' symbol.
formula2	a list which i-th element is an object of class "formula" describing the form of the linear predictor (index) of the i-th continuous equation.
formula3	an object of class "formula" describing the form of the linear predictor (index) of the multinomial equation.
data	a data frame containing the variables in the model.
groups	a matrix which (i, j)-th element is the j-th ordinal category (value starting from 0) of the i-th dependent ordinal variable. Each row of this matrix describes observable (in data) combination of categories - values of the ordinal dependent variables i.e., from the left hand side of formula. Special category '-1' means that variable in the j-th column is unobservable when other dependent variables have particular values i.e., given in the same row. See 'Details' for additional information.
groups2	the same as groups argument but for the continuous dependent variables from formula2. See 'Details' for additional information.
groups3	the same as groups argument but for the multinomial dependent variable from formula3. See 'Details' for additional information.
marginal	a list such that marginal[[i]] represents parameters of the marginal distribution of the random error of the i-th ordered equation and names(marginal)[i] is a name of this distribution. Marginal distributions are the same as in pmnorm.
opt_type	a character representing the optimization function to be used. If opt_type = "optim" then optim will be used. If opt_type = "gena" then gena will be applied i.e., a genetic algorithm. If opt_type = "pso" then pso will be used i.e., a particle swarm optimization.
opt_args	a list of input arguments for the optimization function selected via the opt_type argument. See 'Details' for additional information.
start	a numeric vector of initial parameters' values. It will be used as a starting point for the optimization purposes. It is also possible to provide an object of class 'msel' then its 'par' element will be used as a starting point.
estimator	a character determining estimation method. If estimator = "ml" then maximum-likelihood method will be used. If estimator = "2step" then two-step estimation procedure similar to Heckman's method will be applied. See 'Details' for additional information.
cov_type	a character determining the type of the covariance matrix estimate to be returned. First, suppose that estimator = "ml" then the following estimators are available. If cov_type = "hessian" then negative inverse of Hessian matrix will be applied. If cov_type = "gop" then inverse of Jacobian outer products will be used. If cov_type = "sandwich" then sandwich covariance matrix estimator will be applied. if cov_type = "no" then an identity matrix will be used. If cov_type = "mm" (default) then sandwich estimator is used along with the penalized log-likelihood function. Therefore if there is no regularization then "mm" and "sandwhich" estimators are the same. Second, suppose that estimator = "2step" then available estimators of the asymptotic covariance matrix are identity matrix cov_type = "no" and the method of moments cov_type = "mm".
degrees	a vector of non-negative integers such that degrees[i] represents the degree of the polynomial which elements are selectivity terms associated with the i-th ordered equation. See 'Details' for additional information.
degrees3	a vector of non-negative integers such that degrees3[i] represents the degree of the polynomial which elements are selectivity terms associated with the i-th multinomial equation. See 'Details' for additional information.
n_sim	an integer representing the number of GHK draws when there are more than 3 ordered equations. Otherwise alternative (much more efficient) algorithms will be used to calculate multivariate normal probabilities.
n_cores	a positive integer representing the number of CPU cores used for the parallel computing. If possible it is highly recommend to set it equal to the number of available physical cores.
control	a list of control parameters. See 'Details' for additional information.
regularization	a list of control parameters for regularization. Element ridge_ind is a vector of indexes of parameters subject to regularization according to quadratic (ridge) penalty function. These indexes correspond to parameters from par output element. Set show_ind argument of summary.msel to TRUE to see these indexes. Element ridge_scale is a numeric vector of weights of the ridge penalty function. Element ridge_location is a numeric vector of values to be subtracted from the parameters before they pass into the penalty function. Elements lasso_ind, lasso_scale and lasso_location are the same but for the lasso penalty term.
type3	a character determining the type of the multinomial model to be considered. If type3 = "logit" then multinomial logit model will be used. If type3 = "probit" then multinomial probit model will be applied. See 'Details' for additional information.

Details

This function allows to estimate multivariate and multinomial sample selection and endogenous switching models with multiple outcomes. These models are the systems of ordinal, continuous and multinomial equations described by formula, formula2 and formula3 respectively.

Ordinal equations

Argument formula determines the regressors of the multivariate ordered probit model with the heteroscedastic random errors:

z_{ji}^{*} = w_{ji}\gamma_{j} + \sigma_{ji}^{*}u_{ji},

\sigma_{ji}^{*} = \exp(w_{ji}^{*}\gamma_{j}^{*}), \quad u_{i}\sim N\left(\begin{bmatrix}0\\ \vdots\\ 0\end{bmatrix}, \Sigma\right), i.i.d.,

z_{ji}=\begin{cases} 0\text{, if }z_{ji}^{*}\leq c_{j1}\\ 1\text{, if }c_{j1}<z_{ji}^{*}\leq c_{j2}\\ 2\text{, if }c_{j2}<z_{ji}^{*}\leq c_{j3}\\ \vdots\\ m_{J-1}\text{, if }z_{(J-1)i}^{*}>c_{(J-1)m_{j-1}}\\ \end{cases},

z_{i}=(z_{1i},...,z_{Ji})^{T},\quad u_{i} = (u_{1i},u_{2i},...,u_{Ji})^{T},

i\in\{1,2,...,n\},\quad j\in\{1,2,...,J\},

where:

• n - the number of observations. If there are no omitted observations then n equals to nrow(data).

• J - the number of ordinal equations which equals to length(formula).

• z_{ji}^{*} - an unobservable (latent) value of the j-th dependent ordinal variable.

• z_{ji} - an observable (ordinal) value of the j-th dependent ordinal variable.

• m_{j} - the number of categories of z_{ji}\in\{0,1,...,m_{j}\}.

• c_{jk} - the k-th cut (threshold) of the j-th dependent ordinal variable.

• w_{ji} - the regressors of the j-th mean equation which should be described in formula[[j]].

• \gamma_{j} - the regression coefficients of the j-th mean equation.

• w_{ji}\gamma_{j} - the linear predictor (index) of the j-th mean equation.

• u_{i} - a multivariate normal random vector which elements are standard normal random variables.

• \Sigma - a correlation matrix of u_{i} so \Sigma_{t_{1}t_{2}}=\text{corr}\left(u_{it_{1}}, u_{it_{2}}\right).

• \sigma_{ji}^{*} - a heteroscedastic standard deviation.

• \sigma_{ji}^{*}u_{ji} - the heteroscedastic random errors.

• w_{ji}^{*} - the regressors of the j-th variance equation which should be described in formula[[j]] after '|' symbol.

• \gamma_{j}^{*} - the regression coefficients of the j-th variance equation.

• w_{ji}^{*}\gamma_{j}^{*} - the linear predictor (index) of the j-th variance equation.

Constant terms (intercepts) are excluded from the model for identification purposes. If z_{ji} is a binary variable then -c_{j1} may be interpreted as a constant term of the j-th ordinal equation. If all z_{ji} are binary variables then the model becomes a multivariate probit.

By default the joint distribution of u_{i} is multivariate normal. However by using marginal argument it is possible to consider the joint distribution that is determined by the Gaussian copula and possibly non-normal marginal distributions. Specifically, names(marginal)[i] is the name of the marginal distribution of u_{ji} and marginal[[i]] is the number of parameters of this distribution. The marginal distributions are the same as in pmnorm.

Multinomial equation

Argument formula3 determines the regressors of the multinomial equation:

\tilde{z}_{ji}^{*} = \tilde{w}_{i}\tilde{\gamma}_{j} + \tilde{u}_{ji}, \qquad j\in\{0,1,...,\tilde{J} - 1\},

\tilde{z}_{i}=\underset{t\in\{0,...,\tilde{J}-1\}}{\text{argmax}} \text{ }\tilde{z}_{ti}^{*}, \qquad \tilde{u}_{i} = (\tilde{u}_{0i},\tilde{u}_{1i},..., \tilde{u}_{(\tilde{J}-1)i})^{T},

where:

• \tilde{J} - the number of the alternatives i.e., possible values of the dependent variable of the multinomial equation.

• \tilde{z}_{ji}^{*} - an unobservable (latent) value of the j-th alternative. Usually \tilde{z}_{ji}^{*} is interpreted as a utility of the j-th alternative.

• \tilde{z}_{i} - a selected alternative i.e., one which provides the greatest utility \tilde{z}_{ji}^{*}.

• \tilde{w}_{i} - the regressors of the multinomial equation which should be described in formula3 and assumed to be the same for all the alternatives.

• \tilde{\gamma}_{j} - the regression coefficients of the j-th alternative's equation.

• \tilde{w}_{i}\tilde{\gamma}_{j} - the linear predictor (index) of the j-th alternative's equation.

• \tilde{u}_{i} - a vector of random errors.

For the identification purposes it is assumed that the regression coefficients of the last alternative are zero \tilde{\gamma}_{\tilde{J} - 1} = (0,...,0)^{T}.

Joint distribution of \tilde{u}_{i} depends on the value of type3 argument. If type3 = "logit" then multinomial logit model is considered. It assumes that \tilde{u}_{ji} are independent and their marginal distribution is Gumbel (error value distribution). If type3 = "probit" then multinomial probit model is used so it is assumed that joint distribution of \tilde{u}_{i} is multivariate normal with zero mean and the covariance matrix \tilde{\Sigma}. For identification purposes it is also assumed that \tilde{\Sigma}_{11} = 1 so Var(\tilde{u}_{0i}) = 1. In addition \tilde{u}_{(\tilde{J}-1)i}=0 which implies \tilde{\Sigma}_{t(\tilde{J})}=\tilde{\Sigma}_{(\tilde{J})t}=0 for all t\in\{0,...,\tilde{J}-1\}.

Continuous equations

Argument formula3 determines the regressors of the continuous equations:

y_{r^{(v)}i}^{(v)} = x_{i}^{(v)}\beta_{r^{(v)}}^{(v)} + \varepsilon_{r^{(v)}i}^{(v)},

r^{(v)}\in\{0,1,...,R^{(v)} - 1\}, \quad v\in \{1,...,V\},

where:

• V - the number of continuous equations.

• r^{(v)} - the regime of the v-th continuous equation.

• y_{r^{(v)}i}^{(v)} - the r^{(v)}-th potential outcome (in a sense of the Neyman-Rubin framework) of the v-th continuous equation.

• R^{(v)} - the number of the potential outcomes of the v-th continuous equation.

• x_{i}^{(v)} - the regressors of the v-th continuous equation. They are provided via formula2[[v]].

• \beta_{r^{(v)}}^{(v)} - regression coefficients of the v-th continuous equation in the regime r^{(v)}.

• \varepsilon_{r^{(v)}i}^{(v)} - a random error of the v-th continuous equation in the regime r^{(v)}.

Estimation of the model with multivariate ordinal and multinomial equations

If formula2 is not provided then maximum-likelihood estimator is used estimator = "ml" to estimate the parameters of the multivariate ordinal and multinomial equations.

If both formula and formula3 are provided then parameters of the multivariate ordered and multinomial equations are estimated under the assumption that u_{i} is independent of \tilde{u}_{i}. Therefore the estimates are identical to those obtained by separate estimation of these models. We may relax the independence assumption during future updates.

Estimation of the model with continuous equations

If formula and formula3 are not provided then it is assumed that each continuous equation has only one regime so R^{(v)} = 1 for each v\in\{1,...,V\}.

If estimator = "ml" then maximum-likelihood estimator is used under the assumption that the distribution of random errors is multivariate normal. If estimator = "2step" then V-stage least squares estimator is used. In the latter case v-th equation is estimated by the least squares estimator and for every v_{0} such that v_{0} < v if y_{0i}^{(v_{0})} is included in x_{i}^{(v)} then y_{0i}^{(v_{0})} is substituted with its estimate \hat{y}_{0i}^{(v_{0})} obtained on the v_{0}-th step. Therefore if v_{0} < v then if some endogenous variable appears on the left hand side of formula2[[v]] it should not appear on the right hand side of formula2[[v0]].

Estimation of the model with ordinal outcomes and multivariate ordered selection

Suppose that z_{i} represents the ordinal potential outcomes while the observable values are z_{i}^{(o)} = g^{(1)}(z_{i}), where function g^{(1)}(z_{i}) converts unobservable values of z_{i} to -1. Therefore z_{ji}^{(o)} instead of z_{i} appears on the left hand side of the formula[[j]].

For example, consider a binary variable for the employment status (0 - unemployed, 1 - employed) and an ordinal variable z_{2i} (ranging from 0 to 2) for job satisfaction (0 - unsatisfied, 1- satisfied, 2 - highly satisfied). Then z_{2i} is observable only when z_{1i} equals 1 since job satisfaction observable only for working individuals. Consequently z_{2i}^{(o)} should be equal to -1 (minus one) whenever z_{1i} equals to 0:

z_{i}^{(o)}=g^{(1)}(z_{i})= \begin{cases} (1, 2)\text{, if }z_{i} = (1, 2)\\ (1, 1)\text{, if }z_{i} = (1, 1)\\ (1, 0)\text{, if }z_{i} = (1, 0)\\ (0, -1)\text{, otherwise}\\ \end{cases}

Rows of the matrix groups determine all possible values of the function g^{(1)}(z_{i}). In this particular example matrix groups should have the following form:

\text{groups} = \begin{bmatrix} 1 & 2\\ 1 & 1\\ 1 & 0\\ 0 & -1 \end{bmatrix}.

Notice that in this particular example it is necessary to ensure that in data if z_{1i}^{(o)} equals 0 then z_{2i}^{(o)} equals to -1. Also in this case matrix groups will be created automatically so there is no need to provide it manually.

By using groups argument it is straightforward to consider various other models with ordinal outcomes and multivariate ordered selection mechanisms.

If some values of the ordinal variables z_{ji} are missing (by random) i.e., take NA value then the contribution of other ordinal dependent variables (for the i-th observation) still may be included into the likelihood function by manually substituting NA with -1 in the data. However, ensure that this particular (missing) z_{ji} is not a regressor for other dependent variable that may happen in the hierarchical systems.

It is useful to use groups argument to consider causal inference models. For example, suppose that z_{1i} represents a binary treatment variable for the university diploma (0 - no diploma, 1 - has diploma). Also, z_{2i} is a binary potential outcome for the employment status (0 - unemployed, 1 - employed) of the individual if she has university diploma. Finally, z_{3i} is a binary potential outcome for the employment status (0 - unemployed, 1 - employed) if individual does not have university diploma. Since z_{2i} is observable only if z_{1i} = 1 and z_{3i} is observable only when z_{2i} = 0 we get:

z_{i}^{(o)} = \begin{cases} (1, 1, -1)\text{, if }z_{1i} = 1\text{ and }z_{2i}=1\\ (1, 0, -1)\text{, if }z_{1i} = 1\text{ and }z_{2i}=0\\ (0, -1, 1)\text{, if }z_{1i} = 0\text{ and }z_{3i}=1\\ (0, -1, 0)\text{, if }z_{1i} = 0\text{ and }z_{3i}=0 \end{cases}.

Therefore to estimate this model it is necessary to ensure that in data we have z_{2i}^{(o)} = -1 when z_{1i}^{(o)} = 0 and z_{3i}^{(o)} = -1 if z_{1i}^{(o)} = 1. Also the groups argument should include all possible values of z_{i}^{(o)}:

\text{groups} = \begin{bmatrix} 1 & 1 & -1\\ 1 & 0 & -1\\ 0 & -1 & 1\\ 0 & -1 & 0 \end{bmatrix}.

Estimation of the model with continuous outcomes and multivariate ordered selection

To simplify the notations suppose that there is only one continuous equation with multiple regimes:

y_{ri} = x_{i}\beta_{r} + \varepsilon_{ri},

y_{i} = \begin{cases} \text{unobservable, if }g^{(2)}\left(z_{i}^{(o)}\right) = -1\\ y_{0i}\text{, if }g^{(2)}\left(z_{i}^{(o)}\right) = 0\\ y_{1i}\text{, if }g^{(2)}\left(z_{i}^{(o)}\right) = 1\\ \vdots\\ y_{(R^{*}-1)i}\text{, if }g^{(2)}\left(z_{i}^{(o)}\right) = R^{*} - 1\\ \end{cases},

\varepsilon_{i} = (\varepsilon_{0i},\varepsilon_{1i},...,\varepsilon_{(R^{*}-1)i}),\quad r\in\{0,1,...,R^{*} - 1\},

where:

• y_{i} - an observable continuous outcome.

• r - the index of the potential outcome.

• R^{*} - the number of the regimes.

• y_{ri} - a continuous potential outcome i.e., the value of y_{i} in the regime r.

• x_{i} - the vector of regressors provided in formula2.

• \beta_{r} - the regression coefficients in the r-th regime.

• g^{(2)}(z_{i}^{(o)}) - a function determining which potential outcome is observable depending on the observable values of the ordinal variables z_{i}^{(o)}.

• \varepsilon_{i} - a vector of random errors.

The value of groups2[i] argument specifies the value of g^{(2)}(z_{i}^{(o)}) when z_{i}^{(o)} equals to groups[i, ]. The values of y_{i} in data such that g^{(2)}(z_{i}^{(o)})=-1 should be set to NA.

For example, consider a system of equations for wage y_{i}, employment status z_{1i} (0 - unemployed, 1 - employed) and job satisfaction z_{2i} (0 - unsatisfied, 1- satisfied, 2 - highly satisfied). Notice that wage and job satisfaction are observable only for working individuals. Also suppose that wage is unobservable for unsatisfied workers and observable in different regimes for other workers. Namely, for satisfied workers z_{2i} = 1 we observe y_{0i} and for highly satisfied workers z_{2i} = 2 we observe y_{1i}.

To estimate this model it is necessary to manually specify the structure of the equations via groups and groups2 arguments by providing all possible combinations of the ordinal variables and the regimes of the continuous equation:

\text{groups} = \begin{bmatrix} 1 & 2\\ 1 & 1\\ 1 & 0\\ 0 & -1 \end{bmatrix}, \text{groups2} = \begin{bmatrix} 1\\ 0\\ -1\\ -1 \end{bmatrix}.

Notice that groups2[1] = 1 indicates that when groups[1, ] = c(1, 2) i.e. z_{1i} = 1 and z_{2i} = 2 we observe y_{i} in regime 1 corresponding to the wage of highly satisfied workers y_{1i}. Similarly groups2[2] = 0 indicates that when groups[2, ] = c(1, 1) i.e., z_{1i} = 1 and z_{2i} = 1 we observe y_{i} in the regime 0 corresponding to the wage of the satisfied workers y_{0i}. Also, groups3[3] = -1 means that when groups[3, ] = c(1, 0) i.e., z_{1i} = 1 and z_{2i} = 0 we do not observe the wage of the worker y_{i} since he is unsatisfied. Finally, groups3[4] = -1 means that when groups[4, ] = c(0, -1) i.e., z_{1i} = 0 and z_{2i}^{(o)} = -1 we do not observe the wage of the worker y_{i} since he is unemployed.

If the joint distribution of \varepsilon_{i} and u_{i} is multivariate normal then according to Kossova and Potanin (2018):

y_{ri}=x_{i}\beta_{r}+\sum\limits_{j=1}^{J} \Sigma^{(12)}_{j(r+1)}\lambda_{ji}+\varepsilon_{ri}^{*},

where:

\varepsilon_{ri}^{*} = \varepsilon_{ri} - E(\varepsilon_{i}|z_{1i},...z_{Ji})=\varepsilon_{ri} - \sum\limits_{j=1}^{J} \Sigma^{(12)}_{j(r+1)}\lambda_{ji},

\lambda_{ji}= \lambda_{ji}^{(1)} + \lambda_{ji}^{(2)},\qquad \Sigma^{(12)}_{j(r+1)} = cov(u_{ji}, \varepsilon_{r+1}),

\lambda_{ji}^{(1)}= \begin{cases} 0\text{, if }z_{ji} = 0\\ -\frac{\partial \ln P_{i}^{*}}{\partial a_{ji}}\text{, otherwise } \end{cases}, \qquad \lambda_{ji}^{(2)}= \begin{cases} 0\text{, if }z_{ji} = m_{j} - 1\\ -\frac{\partial \ln P_{i}^{*}}{\partial b_{ji}}\text{, otherwise } \end{cases},

P_{i}^{*}(a_{1i},...,a_{Ji};b_{1i},...,b_{Ji})= P(a_{1i}\leq u_{1i}\leq b_{1i},....,a_{Ji}\leq u_{Ji}\leq b_{Ji}),

a_{ji}= \begin{cases} -\infty\text{, if }z_{ji} = 0\\ \frac{c_{jz_{ji}}-w_{ji}\gamma_{j}}{\sigma_{ji}^{*}}\text{, otherwise } \end{cases}, \qquad b_{ji}= \begin{cases} \infty\text{, if }z_{ji} = m_{j} - 1\\ \frac{c_{j(z_{ji}+1)}-w_{ji}\gamma_{j}}{\sigma_{ji}^{*}}\text{, otherwise } \end{cases}.

Notice that the regression equation has selectivity terms \lambda_{ji} which may be correlated with x_{i}. Therefore until random errors u_{i} and \varepsilon_{i} are correlated the least squares estimator of x_{i} on y_{ri} is inconsistent. To get consistent estimates of \beta_{r} it is possible to use maximum-likelihood estimator = "ml" or two-step (method of moments) estimator = "2step" estimator.

If the two-step estimator is used then on the first step \lambda_{ji} are estimated as the functions of the estimates of the multinomial heteroscedastic ordered probit model. On the second step least squares regression of y_{ri} on x_{i} and the first step estimates \hat{\lambda}_{ji} is used to estimate \beta_{r} and \Sigma^{(12)}_{j(r+1)}.

If the joint distribution of random errors \varepsilon_{i}, u_{i} is not multivariate normal then \lambda_{ji} terms may enter continuous equation non-linearly. Following Newey (2009) and E. Kossova, L. Kupriianova, and B. Potanin (2020) it is assumed that:

y_{ri}= x_{i}\beta_{r} + \zeta_{i}\tau_{r} + \varepsilon_{i}^{*},

where \tau_{r} is a n_{\lambda}-dimensional column vector and:

\zeta_{i} = (\zeta_{1}(\lambda_{i}),..., \zeta_{n_{\lambda}}(\lambda_{i})),\qquad \lambda_{i} = (\lambda_{1i},..., \lambda_{Ji}).

Functions \zeta_{t}(\lambda_{i}) are specified manually by the user in the formula2 inside I(). For example, to specify \zeta_{t}(\lambda_{i}) = \lambda_{1i}\times \lambda_{2i} it is sufficient to have a term I(lambda1 * lambda2) in formula2. Notice that to avoid the confusions no variables in data should have the names containing "lambda". Otherwise these variables will be dropped.

It is possible to specify \zeta_{t}(\lambda_{i}) functions as the polynomial without interaction terms by using degrees argument. Specifically, if degrees[j] = t then lambdaj, I(lambdaj ^ 2),...,I(lambdaj ^ t) terms are added to formula2. However, if degrees argument is used then no functions of lambdaj should be provided manually in formula2. Otherwise it will be assumed that degrees is a vector of zeros. Also if estimator = "2step", there is not lambdaj terms in formula2 and degrees is NA then degrees will be converted in a vector of ones.

If there are multiple continuous equations then formula2 should be a list of formulas. Further, if estimator = "2step" then the second step is a V-stage least squares estimator with lambda terms. If they are provided via degrees argument then it should be a matrix which v-th row corresponds to the v-th continuous equation.

For example, consider previous example with additional continuous equation for working hours y_{i}^{(2)} which does not vary with the satisfaction of the workers z_{2i} but observable only for the employed individuals z_{1i} = 1. To estimate the system with this additional continuous equation simply substitute all y_{i}^{(2)} (such that z_{1i}=0) in data with NA and specify:

\text{groups} = \begin{bmatrix} 1 & 2\\ 1 & 1\\ 1 & 0\\ 0 & -1 \end{bmatrix}, \text{groups2} = \begin{bmatrix} 1 & 0\\ 0 & 0\\ -1 & 0\\ -1 & -1 \end{bmatrix}.

Notice that groups2[, 1] describes the regimes of the wage equation y_{i}^{(1)} while groups[, 2] contains the regimes of the hours equation y_{i}^{(2)}. Note that formula of the first equation (wage) should be specified in formula2[[1]] and formula of the second equation should be provided via formula2[[2]] i.e., as the first and the second elements in a formula2 list correspondingly.

If marginal argument is used then aforementioned formula of \lambda_{ji} is slightly modified to address for the non-normal marginal distribution of u_{ji}.

Estimation of the model with continuous outcomes and multinomial selection

The only difference with the previous model is that the observable value of the continuous equation depends on the value of the multinomial equation described in formula3. Therefore g^{(2)}(z_{i}^{(o)}) is substituted with g^{(3)}(\tilde{z}_{i}). Also groups3 argument instead of groups is used. Since there is only a single multinomial equation argument groups3 is a vector. If groups3[k] = q and groups2[k, v] = r then \tilde{z}_{i} = q implies y_{i}^{(v)} = y_{ri}^{(v)}. Remind that a special value r = -1 implies that y_{i}^{(v)} is unobservable.

Only two-step estimator = "2step" estimator of this model is available that is similar to the one described above. The only difference is the formula used to estimate selectivity terms. If type = "logit" then two-step estimator of Dubin-McFadden (1984) is used so selectivity terms are as follows:

\lambda_{ji} = \begin{cases} -\ln P(\tilde{z}_{i} = j)\text{, if }\tilde{z}_{i} = j\\ \frac{P(\tilde{z}_{i} = j)\ln P(\tilde{z}_{i} = j)}{1 - P(\tilde{z}_{i} = j)}\text{, otherwise} \end{cases},

where j\in\{0,...,\tilde{J} - 1\} and:

P(\tilde{z}_{i} = j) = \begin{cases} \frac{e^{(\tilde{w}_{i}\tilde{\gamma}_{j})}}{1 + \sum\limits_{q=0}^{\tilde{J} - 2}e^{(\tilde{w}_{i}\tilde{\gamma}_{q})}} \text{, if }j\ne \tilde{J} - 1\\ \frac{1}{1 + \sum\limits_{q=0}^{\tilde{J} - 2}e^{(\tilde{w}_{i}\tilde{\gamma}_{q})}} \text{, otherwise}\\ \end{cases}.

If type = "probit" then two-step estimator of Kossova and Potanin (2022) is used with the selectivity terms of the form:

\lambda_{ji} = \left(A^{(\tilde{z}_{i})}\lambda_{i}^{*}\right)_{j},

where j\in\{0,...,\tilde{J} - 2\} and:

A^{(j)}_{t_{1}t_{2}} = \begin{cases} 1\text{, if }t_{1} = j + 1\\ -1\text{, if }t_{1} < j + 1\text{ and }t_{1} = t_{2}\\ -1\text{, if }t_{1} > j + 1\text{ and }t_{1} = t_{2} + 1\\ 0\text{, otherwise } \end{cases},

t_{1},t_{2}\in\{1,...,\tilde{J} - 1\},

\lambda_{i}^{*} = \nabla \ln P\left(\tilde{z}_{i}\right) = \nabla \ln F_{\tilde{u}^{(ji)}}\left(\tilde{z}_{1}^{(ji)},..., \tilde{z}_{\tilde{J}-1}^{(ji)}\right)= \left(\lambda_{1i}^{*},...,\lambda_{(\tilde{J}-1)i}^{*} \right)^{T},

\tilde{u}^{(ji)} = \left(\tilde{u}_{0i}-\tilde{u}_{ji},\tilde{u}_{1i}-\tilde{u}_{ji},..., \tilde{u}_{(j-1)i}-\tilde{u}_{ji},\tilde{u}_{(j+1)i}-\tilde{u}_{ji},..., \tilde{u}_{(\tilde{J} - 1)i}-\tilde{u}_{ji}\right),

\tilde{z}^{(ji)} = \tilde{w}_{i}\left( \left(\tilde{\gamma}_{j}-\tilde{\gamma}_{0}\right), \left(\tilde{\gamma}_{j}-\tilde{\gamma}_{1}\right),..., \left(\tilde{\gamma}_{j}-\tilde{\gamma}_{j-1}\right), \left(\tilde{\gamma}_{j}-\tilde{\gamma}_{j+1}\right),..., \left(\tilde{\gamma}_{j}-\tilde{\gamma}_{\tilde{J} - 1}\right)\right).

Probabilities P(z_{i}) are calculated by using a cumulative distribution function of the multivariate normal distribution with zero mean and the covariance matrix of \tilde{u}^{(ji)}.

In formula2 selectivity terms associated with the multinomial equation should be named lambdaj_mn instead of lambdaj. Argument degrees3 is similar to degrees.

Consider a simple example of this model. Suppose that \tilde{z}_{i} is a multinomial variable for the employment status (0 - unemployed, 1 - working in IT sector, 2 - working in other sector). Wage y_{i} is unobservable for unemployed \tilde{z}_{i} = 0, equals to y_{0i} in IT sector \tilde{z}_{i} = 1 and equals to y_{1i} in other sectors \tilde{z}_{i} = 2. To estimate this model set groups3 = (0, 1, 2) and groups2 = (-1, 0, 1). Then substitute all y_{i} such that \tilde{z}_{i} = 0 with NA.

Estimation of the model with continuous outcomes and mixed selection

It is possible to consider the model with continuous outcomes and both multinomial and ordinal selection equations. Remind that it is assumed that random errors of the ordered and multinomial equations are independent. Therefore if formula and formula3 are provided then both lambdaj and lambdaj_mn are included in formula2. Only two-step estimator estimator = "2step" is available for this model.

Missing values

If any of the left hand side variables (regressors) of formula[[j]] is missing then the right hand side variable of formula[[j]] will be set to NA in data. Similar is true for formula2 and formula3.

Additional information

Functions pmnorm and dmnorm are used internally for calculation of multivariate normal probabilities, densities and their derivatives.

Currently control has no input arguments intended for the users. This argument is used for some internal purposes of the package.

Optimization always starts with optim. If opt_type = "gena" or opt_type = "pso" then gena or pso is used to proceed optimization starting from initial point provided by optim. Manual arguments to optim should be provided in a form of a list through opt_args$optim. Similarly opt_args$gena and opt_args$pso provide manual arguments to gena and pso. For example to provide Nelder-Mead optimizer to optim and restrict the number of genetic algorithm iterations to 10 make opt_args = list(optim = list(method = "Nelder-Mead"), gena = list(maxiter = 10)).

If estimator = "2step" then it is possible to precalculate the first step model with msel function and pass it through the formula argument. It allows to experiment with various formula2 and degrees specifications without extra computational burden associated with the first step estimation.

If estimator = "2step" then the method of moments estimator of the asymptotic covariance matrix is used as described in Meijer and Wansbeek (2007).

Value

This function returns an object of class "msel". It is a list which may contain the following elements:

• par - a vector of the estimates of the parameters.

• estimator - the same as the estimator input argument.

• type3 - the same as the type3 input argument.

• formula - the same as formula input argument but all elements are coerced to "formula" type.

• formula2 - the same as formula2 input argument but all elements are coerced to "formula" type.

• formula3 - the same as formula3 input argument but all elements are coerced to "formula" type.

• model1 - an object of class "msel" with the first step estimation results.

• data - the same as data input argument but without missing (by random) values.

• W_mean - a list such that W_mean[[j]] is a matrix of the regressors w_{ji} of the mean equation of z_{ji}^{*}.

• W_var - a list such that W_var[[j]] is a matrix of the regressors w_{ji}^{*} of the variance equation of z_{ji}^{*}.

• X - a list such that X[[v]] is a numeric matrix of regressors x_{i}^{(v)} of the v-th continuous equation y_{i}^{(v)}.

• W_mn - a matrix of the regressors \tilde{w}_{i} of the multinomial equation of \tilde{z}_{ji}^{*}.

• dependent - a numeric matrix which j-th column dependent[, j] is a vector of the ordinal dependent variable z_{ji}^{(o)} values.

• dependent2 - a numeric matrix which v-th column dependent2[, v] is a vector of the continuous dependent variable y_{i}^{(v)} values.

• dependent3 - a numeric vector of values of the multinomial dependent variable \tilde{z}_{j}.

• groups - the same as groups input argument or automatically generated matrix representing the structure of the system of equations. Please, see 'Details' section above for more information.

• groups2 - the same as groups2 input argument or automatically generated matrix representing the structure of the system of equations. Please, see 'Details' section above for more information.

• groups3 - the same as groups3 input argument or automatically generated matrix representing the structure of the system of equations. Please, see 'Details' section above for more information.

• marginal - the same as marginal input argument.

• ind - a list containing some indexes partition of the model (not intended for the users).

• start - the same as the start input argument.

• twostep - a list such that twostep[[v]][[r + 1]] is an object of class "lm" associated with the second step estimates of the v-th equation in the regime r.

• y_pred - a numeric matrix such that y_pred[, v] is a second step prediction of the y_{i}^{(v)}.

• coef - a list which j-th element coef[[j]] is a numeric vector representing \hat{\gamma}_{j}.

• coef_var - a list which j-th element coef_var[[j]] is a numeric vector representing \hat{\gamma}_{j}^{*}.

• cuts - a list which j-th element cuts[[j]] is a numeric vector representing \hat{c}_{j}.

• coef2 - a list of numeric matrices such that coef2[[v]][, r + 1] is a numeric vector representing \hat{\beta}_{r^{(v)}}^{(v)}.

• sigma - a numeric matrix such that sigma[j, t] is a numeric value representing \widehat{Cov}(u_{ji}, u_{ti}).

• var2 - a list of numeric vectors such that var2[[v]][r + 1] represents \widehat{Var}(\varepsilon_{ri}^{(v)}).

• cov2 - a list of numeric matrices which element sigma2[[v]][j, r + 1] represents \widehat{Cov}(u_{ji}, \varepsilon_{ri}^{(v)})).

• sigma2 - a list of numeric matrices representing the estimates of the covariances between random errors of the continuous equations in different regimes \widehat{Cov}(\varepsilon_{r^{(v)}i}^{(v)}, \varepsilon_{r^{(t)}i}^{(t)}).

• marginal_par - a list such that marginal_par[[j]] is a numeric vector of estimates of the parameters of the marginal distribution of u_{ji}.

• coef3 - a numeric matrix such that coef3[j + 1, ] is a numerc vector representing \hat{\tilde{\gamma}}_{j}.

• sigma3 - a numeric matrix such that sigma3[j + 1, t + 1] is a numeric value representing \widehat{Cov}(\tilde{u}_{ji}, \tilde{u}_{ti}).

• lambda - a numeric matrix such that lambda[i, j] corresponds to \hat{\lambda}_{ji} of the ordinal equations.

• lambda_mn - a numeric matrix such that lambda_mn[i, j] corresponds to \hat{\lambda}_{ji} of the multinomial equation.

• H - if estimator = "ml" then H is a Hessian matrix of the log-likelihood function. If estimator = "2step" then H is a numeric matrix of the derivatives of mean sample scores respect to the estimated parameters.

• J - if estimator = "ml" then J is a Jacobian matrix of the log-likelihood function. If estimator = "2step" then J is a numeric matrix such that J[, s] is a vector of sample scores associated with the parameter par[s].

• cov_type - the same as cov_type input argument.

• cov - an estimate of the asymptotic covariance matrix of the parameters' estimator.

• tbl - a list of special tables used to create a summary (not intended for the users).

• se - a numeric vector of standard errors of the estimates.

• p_value - a numeric vectors of the p-values of the tests on the significance of the estimated parameters where the null hypothesis is that corresponding parameter is zero.

• logLik - the value of log-likelihood function at par.

• other - a list of additional variables not intended for the users.

It is highly recommended to get the estimates of the parameters via coef.msel function.

References

W. K. Newey (2009). Two-step series estimation of sample selection models. The Econometrics Journal, vol. 12(1), pages 217-229.

E. Kossova, B. Potanin (2018). Heckman method and switching regression model multivariate generalization. Applied Econometrics, vol. 50, pages 114-143.

E. Kossova, L. Kupriianova, B. Potanin (2020). Parametric and semiparametric multivariate sample selection models estimators' accuracy: Comparative analysis on simulated data. Applied Econometrics, vol. 57, pages 119-139.

E. Kossova, B. Potanin (2022). Estimation of Gaussian multinomial endogenous switching model. Applied Econometrics, vol. 67, pages 121-143.

E. Meijer, T. Wansbeek (2007). The sample selection model from a method of moments perspecrive. Econometric Reviews, vol. 26(1), pages 25-51.

Examples

# ---------------------------------------
# CPS data example
# ---------------------------------------

# Set seed for reproducibility
set.seed(123)

# Upload data
data(cps)

# Prepare the variable for education
cps$educ                    <- NA
cps$educ[cps$basic == 1]    <- 0
cps$educ[cps$bachelor == 1] <- 1
cps$educ[cps$master == 1]   <- 2

# Labor supply (probit) model
f_work <- work ~ age + I(age ^ 2) + bachelor + master + health + 
                 slwage + nchild
model1 <- msel(f_work, data = cps)
summary(model1)
    
# Education choice (ordered probit) model
f_educ <- educ ~ age + I(age ^ 2) + sbachelor + smaster
model2 <- msel(f_educ, data = cps)
summary(model2)   

# Education choice (multinomial logit) model
model3 <- msel(formula3 = f_educ, data = cps, type3 = "logit")
summary(model3)  

# Education choice (multinomial probit) model
model4 <- msel(formula3 = f_educ, data = cps, type3 = "probit")
summary(model4) 

# Labor supply with endogenous ordinal education 
# treatment (recursive or hierarchical ordered probit) model
model5 <- msel(list(f_work, f_educ), data = cps)  
summary(model5)         

# Sample selection (on employment) Heckman's model
f_lwage                  <- lwage ~ age + I(age ^ 2) + 
                                    bachelor + master + health
model6                   <- msel(f_work, f_lwage, data = cps)
summary(model6)

# Ordinal endogenous education treatment with non-random 
# sample selection into employment
model7 <- msel(list(f_work, f_educ), f_lwage, data = cps)
summary(model7)
 
# Ordinal endogenous switching on education model with 
# non-random selection into employment
groups   <- cbind(c(1, 1, 1, 0, 0, 0),
                  c(0, 1, 2, 0, 1, 2))
groups2  <- matrix(c(0, 1, 2, -1, -1, -1), ncol = 1)
f_lwage2 <- lwage ~ age + I(age ^ 2) + health
model8   <- msel(list(f_work, f_educ), f_lwage2,
                      groups = groups, groups2 = groups2,
                      data   = cps)
summary(model8)

# Multinomial endogenous switching on education model with 
# non-random selection into employment
groups  <- matrix(c(1, 1, 1, 0, 0, 0), ncol = 1)
groups2 <- matrix(c(0, 1, 2, -1, -1, -1), ncol = 1)
groups3 <- c(0, 1, 2, 0, 1, 2)
model9  <- msel(f_work, f_lwage2, f_educ,
                groups    = groups,  groups2 = groups2, 
                groups3   = groups3, data    = cps,
                estimator = "2step")
summary(model9)


# ---------------------------------------
# Simulated data example 1
# Ordered probit and other univariate
# ordered choice models
# --------------------------------------

# ---
# Step 1
# Simulation of the data
# ---

# Load required package
library("mnorm")

# Set seed for reproducibility
set.seed(123)

# The number of observations
n <- 1000

# Regressors (covariates)
w1 <- runif(n = n, min = -1, max = 1)
w2 <- runif(n = n, min = -1, max = 1)

# Random errors
u <- rnorm(n = n, mean = 0, sd = 1)

# Coefficients
gamma <- c(-1, 2)

# Linear predictor (index)
li <- gamma[1] * w1 + gamma[2] * w2

# Latent variable
z_star <- li + u

# Cuts
cuts <- c(-1, 0.5, 2)

# Observable ordinal outcome
z                                           <- rep(0, n)
z[(z_star > cuts[1]) & (z_star <= cuts[2])] <- 1
z[(z_star > cuts[2]) & (z_star <= cuts[3])] <- 2
z[z_star > cuts[3]]                         <- 3
table(z)

# Data
data <- data.frame(w1 = w1, w2 = w2, z = z)

# ---
# Step 2
# Estimation of the parameters
# ---

# Estimation
model <- msel(z ~ w1 + w2, data = data)
summary(model)

# Compare the estimates and true values of the parameters
  # regression coefficients
gamma_est <- coef(model, type = "coef", eq = 1)
cbind(true = gamma, estimate = gamma_est)
  # cuts
cuts_est <- coef(model, type = "cuts", eq = 1)
cbind(true = cuts, estimate = cuts_est)   

# ---
# Step 3
# Estimation of the probabilities and marginal effects
# ---

# Predict conditional probability of the dependent variable
# equals to 2 for every observation in a sample.
# P(z = 2 | w)
prob <- predict(model, group = 2, type = "prob")
head(prob)

# Calculate mean marginal effect of w2 on P(z = 1 | w)
mean(predict(model, group = 1, type = "prob", me = "w2"))

# Calculate probabilities and marginal effects
# for manually provided observations.
  # new data
newdata <- data.frame(z = c(1, 1), 
                      w1 = c(0.5, 0.2), 
                      w2 = c(-0.3, 0.8))
  # probability P(z = 2 | w)
predict(model, group = 2, type = "prob", newdata = newdata)
  # linear predictor (index)
predict(model, type = "li", newdata = newdata)  
  # marginal effect of w1 on P(z = 2 | w)
predict(model, group = 2, type = "prob", newdata = newdata, me = "w1")
  # marginal effect of w1 and w2 on P(z = 3 | w)
predict(model, group   = 3,       type = "prob", 
               newdata = newdata, me   = c("w1", "w2"))
  # marginal effect of w2 on the linear predictor (index)
predict(model, group = 2, type = "li", newdata = newdata, me = "w2")
  # discrete marginal effect:
  # P(z = 2 | w1 = 0.5, w2) - P(z = 2 | w1 = 0.2, w2)
predict(model, group = 2,    type = "prob", newdata = newdata, 
               me    = "w2", eps  = c(0.2, 0.5))

# ---
# Step 4
# Ordered logit model
# ---

# Estimate ordered logit model with a unit variance
# that is just a matter of reparametrization i.e.,
# do not affect signs and significance of the coefficients
# and dot not affect at all the marginal effects
logit <- msel(z ~ w1 + w2, data = data, marginal = list("logistic" = 0))
summary(logit)

# Compare ordered probit and ordered logit models
# using Akaike and Bayesian information criteria
  # AIC
c(probit = AIC(model), logit = AIC(logit))
  # BIC
c(probit = BIC(model), logit = BIC(logit))

# Estimate some probabilities and marginal effects
  # probability P(z = 1 | w)
predict(logit, group = 1, type = "prob", newdata = newdata)
  # marginal effect of w2 on P(z = 1 | w)
predict(logit, group = 1, type = "prob", newdata = newdata, me = "w2")

# ---
# Step 5
# Semiparametric ordered choice model with 
# Gallant and Nychka distribution
# ---

# Estimate semiparametric model
pgn <- msel(z ~ w1 + w2, data = data, marginal = list("PGN" = 2))
summary(pgn)

# Estimate some probabilities and marginal effects
  # probability P(z = 3 | w)
predict(pgn, group = 3, type = "prob", newdata = newdata)
  # marginal effect of w2 on P(z = 3 | w)
predict(pgn, group = 3, type = "prob", newdata = newdata, me = "w2")

# Test the normality assumption via the likelihood ratio test
summary(lrtest_msel(model, pgn))

# Test the normality assumption via the Wald test
test_fn <- function(object)
{
  marginal_par <- coef(object, type = "marginal", eq = 1)
  return(marginal_par)
}
test_result <- test_msel(object = pgn, test = "wald", fn = test_fn)
summary(test_result)

 

# ---------------------------------------
# Simulated data example 2
# Heteroscedastic ordered probit model
# ---------------------------------------

# Load required package
library("mnorm")

# ---
# Step 1
# Simulation of the data
# ---

# Set seed for reproducibility
set.seed(123)

# The number of observations
n <- 1000

# Regressors (covariates)
w1 <- runif(n = n, min = -1, max = 1)
w2 <- runif(n = n, min = -1, max = 1)
w3 <- runif(n = n, min = -1, max = 1)

# Random errors
u <- rnorm(n, mean = 0, sd = 1)

# Coefficients of the mean equation
gamma <- c(-1, 2)

# Coefficients of the variance equation
gamma_het <- c(0.5, -1)

# Linear predictor (index) of the mean equation
li <- gamma[1] * w1 + gamma[2] * w2

# Linear predictor (index) of the variance equation
li_het <- gamma_het[1] * w2 + gamma_het[2] * w3

# Heteroscedastic stdandard deviation
# i.e., value of the variance equation
sd_het <- exp(li_het)

# Latent variable
z_star <- li + u * sd_het

# Cuts
cuts <- c(-1, 0.5, 2)

# Observable ordinal outcome
z                                           <- rep(0, n)
z[(z_star > cuts[1]) & (z_star <= cuts[2])] <- 1
z[(z_star > cuts[2]) & (z_star <= cuts[3])] <- 2
z[z_star > cuts[3]]                         <- 3
table(z)

# Data
data <- data.frame(w1 = w1, w2 = w2, w3 = w3, z = z)

# ---
# Step 2
# Estimation of the parameters
# ---

# Estimation
model <- msel(z ~ w1 + w2 | w2 + w3,
                   data = data)
summary(model)

# Compare the estimates and true values of the parameters
  # regression coefficients of the mean equation
gamma_est <- coef(model, type = "coef", eq = 1)
cbind(true = gamma, estimate = gamma_est)
  # regression coefficients of the variance equation
gamma_het_est <- coef(model, type = "coef_var", eq = 1)
cbind(true = gamma_het, estimate = gamma_het_est)  
  # cuts
cuts_est <- coef(model, type = "cuts", eq = 1)
cbind(true = cuts, estimate = cuts_est)   

# Likelihood-ratio test for the homoscedasticity
model0 <- msel(z ~ w1 + w2, data = data)
summary(lrtest_msel(model, model0))

# Wald test for the homoscedasticity
test_fn <- function(object)
{
  val <- coef(object, type = "coef_var", eq = 1)
  return(val)
}
test_result <- test_msel(model, test = "wald", fn = test_fn)
summary(test_result)

# ---
# Step 3
# Estimation of the probabilities and marginal effects
# ---

# Predict probability of the dependent variable
# equals to 2 for every observation in a sample
# P(z = 2 | w)
prob <- predict(model, group = 2, type = "prob")
head(prob)

# Calculate mean marginal effect of w2 on P(z = 1 | w)
mean(predict(model, group = 1, type = "prob", me = "w2"))

# Estimate conditional probabilities, linear predictors (indexes) and 
# heteroscedastic standard deviations for manually 
# provided observations.
  # new data
newdata <- data.frame(z = c(1, 1), 
                      w1 = c(0.5, 0.2), 
                      w2 = c(-0.3, 0.8),
                      w3 = c(0.6, 0.1))
  # probability P(z = 2 | w)
predict(model, group = 2, type = "prob", newdata = newdata)
  # linear predictor (index)
predict(model, type = "li", newdata = newdata)
  # standard deviation
predict(model, type = "sd", newdata = newdata)
  # marginal effect of w3 on P(z = 3 | w)
predict(model, group = 3, type = "prob", newdata = newdata, me = "w3")
  # marginal effect of w2 on the standard error
predict(model, group = 2, type = "sd", newdata = newdata, me = "w2")
  # discrete marginal effect: 
  # P(Z = 2 | w1 = 0.5, w2) - P(Z = 2 | w1 = 0.2, w2)
predict(model, group = 2, type = "prob", newdata = newdata,
        me = "w2", eps = c(0.2, 0.5))



# ---------------------------------------
# Simulated data example 3
# Bivariate ordered probit model with 
# heteroscedastic second equation
# ---------------------------------------

# Load required package
library("mnorm")

# ---
# Step 1
# Simulation of the data
# ---

# Set seed for reproducibility
set.seed(123)

# The number of observations
n <- 1000

# Regressors (covariates)
w1 <- runif(n = n, min = -1, max = 1)
w2 <- runif(n = n, min = -1, max = 1)
w3 <- runif(n = n, min = -1, max = 1)
w4 <- runif(n = n, min = -1, max = 1)

# Covariance matrix of random errors
rho <- 0.5
sigma <- matrix(c(1,   rho,
                  rho, 1), 
                nrow = 2)

# Random errors
u <- mnorm::rmnorm(n = n, mean = c(0, 0), sigma = sigma)

# Coefficients
gamma1 <- c(-1, 2)
gamma2 <- c(1, 1.5)

# Coefficients of the variance equation
gamma2_het <- c(0.5, -1)

# Linear predictors (indexes)
li1 <- gamma1[1] * w1 + gamma1[2] * w2
li2 <- gamma2[1] * w2 + gamma2[2] * w3

# Linear predictor (index) of the variance equation
li2_het <- gamma2_het[1] * w2 + gamma2_het[2] * w4

# Heteroscedastic stdandard deviation
# i.e. value of variance equation
sd2_het <- exp(li2_het)

# Latent variables
z1_star <- li1 + u[, 1]
z2_star <- li2 + u[, 2] * sd2_het

# Cuts
cuts1 <- c(-1, 0.5, 2)
cuts2 <- c(-2, 0)

# Observable ordinal outcome
  # first outcome
z1                                               <- rep(0, n)
z1[(z1_star > cuts1[1]) & (z1_star <= cuts1[2])] <- 1
z1[(z1_star > cuts1[2]) & (z1_star <= cuts1[3])] <- 2
z1[z1_star > cuts1[3]]                           <- 3
  # second outcome
z2                                               <- rep(0, n)
z2[(z2_star > cuts2[1]) & (z2_star <= cuts2[2])] <- 1
z2[z2_star > cuts2[2]]                           <- 2
  # distribution
table(z1, z2)

# Data
data <- data.frame(w1 = w1, w2 = w2, w3 = w3, 
                   w4 = w4, z1 = z1, z2 = z2)

# ---
# Step 2
# Estimation of the parameters
# ---

# Estimation
model <- msel(list(z1 ~ w1 + w2,
                   z2 ~ w2 + w3 | w2 + w4),
              data = data)
summary(model)

# Compare the estimates and true values of parameters
  # regression coefficients of the first equation
gamma1_est <- coef(model, type = "coef", eq = 1)
cbind(true = gamma1, estimate = gamma1_est)
  # regression coefficients of the second equation
gamma2_est <- coef(model, type = "coef", eq = 2)
cbind(true = gamma2, estimate = gamma2_est)
  # cuts of the first equation
cuts1_est <- coef(model, type = "cuts", eq = 1)
cbind(true = cuts1, estimate = cuts1_est)   
  # cuts of the second equation
cuts2_est <- coef(model, type = "cuts", eq = 2)
cbind(true = cuts2, estimate = cuts2_est) 
  # correlation coefficients
rho_est <- coef(model, type = "cov1", eq = c(1, 2))
cbind(true = rho, estimate = rho_est)  
  # regression coefficients of the variance equation
gamma2_het_est <- coef(model, type = "coef_var", eq = 2)
cbind(true = gamma2_het, estimate = gamma2_het_est)

# ---
# Step 3
# Estimation of the probabilities and linear predictors (indexes)
# ---

# Predict probability P(z1 = 2, z2 = 0 | w)
prob <- predict(model, group = c(2, 0), type = "prob")
head(prob)

# Calculate mean marginal effect of w2 on:
  # P(z1 = 1 | w)
mean(predict(model, group = c(1, -1), type = "prob", me = "w2"))
  # P(z1 = 1, z2 = 0 | w)
mean(predict(model, group = c(1, 0), type = "prob", me = "w2"))

# Calculate conditional probabilities and linear predictors (indexes)
# for the manually provided observations.
  # new data
newdata <- data.frame(z1 = c(1, 1), 
                      z2 = c(1, 1),
                      w1 = c(0.5, 0.2), 
                      w2 = c(-0.3, 0.8),
                      w3 = c(0.6, 0.1),
                      w4 = c(0.3, -0.5))
  # probability P(z1 = 2, z2 = 0 | w)
predict(model, group = c(2, 0), type = "prob", newdata = newdata)
  # linear predictor (index)
predict(model, type = "li", newdata = newdata)  
  # marginal probability P(z2 = 1 | w)
predict(model, group = c(-1, 1), type = "prob", newdata = newdata) 
  # marginal probability P(z1 = 3 | w)
predict(model, group = c(3, -1), type = "prob", newdata = newdata)
  # conditional probability P(z1 = 2 | z2 = 0, w)
predict(model, group = c(2, 0), given_ind = 2,
        type = "prob", newdata = newdata) 
  # conditional probability P(z2 = 1 | z1 = 3, w)
predict(model, group = c(3, 1), given_ind = 1,
        type = "prob", newdata = newdata) 
  # marginal effect of w4 on P(Z2 = 2 | w)
predict(model, group = c(-1, 2),
        type = "prob", newdata = newdata, me = "w4")   
  # marginal effect of w4 on P(z1 = 3, Z2 = 2 | w)
predict(model, group = c(3, 2),
        type = "prob", newdata = newdata, me = "w4") 
  # marginal effect of w4 on P(z1 = 3 | z2 = 2, w)
predict(model, group = c(3, 2), given_ind = 2,
        type = "prob", newdata = newdata, me = "w4")         

# ---
# Step 4
# Replication under the non-random sample selection
# ---

# Suppose that z2 is unobservable when z1 = 1 or z1 = 3
z2[(z1 == 1) | (z1 == 3)] <- -1
data$z2 <- z2

# Replicate the estimation procedure
model <- msel(list(z1 ~ w1 + w2,
                   z2 ~ w2 + w3 | w2 + w4),
              cov_type = "gop", data = data)
summary(model)

# Compare estimates and true values of the parameters
  # regression coefficients of the first equation
gamma1_est <- coef(model, type = "coef", eq = 1)
cbind(true = gamma1, estimate = gamma1_est)
  # regression coefficients of the second equation
gamma2_est <- coef(model, type = "coef", eq = 2)  
cbind(true = gamma2, estimate = gamma2_est)
  # cuts of the first equation
cuts1_est <- coef(model, type = "cuts", eq = 1)  
cbind(true = cuts1, estimate = cuts1_est)   
  # cuts of the second equation
cuts2_est <- coef(model, type = "cuts", eq = 2)  
cbind(true = cuts2, estimate = cuts2_est) 
  # correlation coefficient
rho_est <- coef(model, type = "cov1", eq = c(1, 2))
cbind(true = rho, estimate = rho_est)  
  # regression coefficients of the variance equation
gamma2_het_est <- coef(model, type = "coef_var", eq = 2)  
cbind(true = gamma2_het, estimate = gamma2_het_est)

# ---
# Step 5
# Semiparametric model with marginal logistic and PGN distributions
# ---

# Estimate the model
model <- msel(list(z1 ~ w1 + w2,
                   z2 ~ w2 + w3 | w2 + w4),
              data = data, marginal = list(PGN = 3, logistic = NULL))
summary(model)



# ---------------------------------------
# Simulated data example 4
# Heckman model with ordinal 
# selection mechanism
# ---------------------------------------

# Load required package
library("mnorm")

# ---
# Step 1
# Simulation of the data
# ---

# Set seed for reproducibility
set.seed(123)

# The number of observations
n <- 1000

# Regressors (covariates)
w1 <- runif(n = n, min = -1, max = 1)
w2 <- runif(n = n, min = -1, max = 1)
w3 <- runif(n = n, min = -1, max = 1)

# Random errors
rho    <- 0.5
var.y  <- 0.3
sd.y   <- sqrt(var.y)
sigma  <- matrix(c(1,          rho * sd.y,
                   rho * sd.y, var.y),
                 nrow = 2)
errors <- mnorm::rmnorm(n = n, mean = c(0, 0), sigma = sigma)
u      <- errors[, 1]
eps    <- errors[, 2]

# Coefficients
gamma <- c(-1, 2)
beta  <- c(1, -1, 1)

# Linear predictor (index)
li   <- gamma[1] * w1 + gamma[2] * w2
li.y <- beta[1] + beta[2] * w1 + beta[3] * w3

# Latent variable
z_star <- li + u
y_star <- li.y + eps

# Cuts
cuts <- c(-1, 0.5, 2)

# Observable ordered outcome
z                                           <- rep(0, n)
z[(z_star > cuts[1]) & (z_star <= cuts[2])] <- 1
z[(z_star > cuts[2]) & (z_star <= cuts[3])] <- 2
z[z_star > cuts[3]]                         <- 3
table(z)

# Observable continuous outcome such
# that outcome 'y' is observable only 
# when 'z > 1' and unobservable otherwise
# i.e. when 'z <= 1' we code 'y' as 'NA'
y         <- y_star
y[z <= 1] <- NA

# Data
data <- data.frame(w1 = w1, w2 = w2, w3 = w3, 
                   z = z, y = y)

# ---
# Step 2
# Estimation of parameters
# ---

# Estimation
model <- msel(z ~ w1 + w2, y ~ w1 + w3, data = data)
summary(model)

# Compare estimates and true values of the parameters
  # regression coefficients of the ordinal equation
gamma_est <- coef(model, type = "coef", eq = 1)
cbind(true = gamma, estimate = gamma_est)
  # cuts
cuts_est <- coef(model, type = "cuts", eq = 1)
cbind(true = cuts, estimate = cuts_est)   
  # regression coefficients of the continuous equation
beta_est <- coef(model, type = "coef2", eq2 = 1, regime = 0)
cbind(true = beta, estimate = beta_est)  
  # variance
var.y_est <- coef(model, type = "var", eq2 = 1, regime = 0)
cbind(true = var.y, estimate = var.y_est)
  # covariance
cov_est <- coef(model, type = "cov12", eq = 1, eq2 = 1)
cbind(true = rho * sd.y, estimate = cov_est)

# ---
# Step 3
# Estimation of the expectations and marginal effects
# ---

# New data
newdata <- data.frame(z = 1, 
                      y = 1,
                      w1 = 0.1, 
                      w2 = 0.2,
                      w3 = 0.3)

# Predict the unconditional expectation of the continuous outcome
# E(y | w)
predict(model, group = -1, group2 = 0, newdata = newdata)

# Predict the conditional expectations of the continuous outcome
  # E(y | z = 2, w)
predict(model, group = 2, group2 = 0, newdata = newdata)
  # E(y | z = 0, w)
predict(model, group = 0, group2 = 0, newdata = newdata)

# ---
# Step 4
# Classic Heckman's two-step estimation procedure
# ---

# Estimate the model by using the two-step estimator
model_ts <- msel(z ~ w1 + w2, y ~ w1 + w3,
                 data = data, estimator = "2step")
summary(model_ts)

# Check the estimates accuracy
tbl <- cbind(true    = beta, 
             ml      = model$coef2[[1]][1, ],
             twostep = model_ts$coef2[[1]][1, -4])  
print(tbl)

# ---
# Step 5
# Semiparametric estimation procedure
# ---

# Estimate the model using Lee's method 
# assuming logistic distribution of the
# random errors of the selection equation
model_Lee <- msel(z ~ w1 + w2, 
                  y ~ w1 + w3,
                  data = data, estimator = "2step",
                  marginal = list(logistic = NULL))
summary(model_Lee)

# One step estimation is also available as well
# as more complex marginal distributions.
# Consider random errors in selection equation
# following PGN distribution with three parameters.
model_sp <- msel(z ~ w1 + w2, 
                 y ~ w1 + w3,
                 data = data, marginal = list(PGN = 3))
summary(model_sp)

# To additionally relax normality assumption of
# random error of continuous equation it is possible
# to use Newey's two-step procedure.
model_Newey <- msel(z ~ w1 + w2, 
                    y ~ w1 + w3,
                    data      = data,    marginal = list(logistic = 0),
                    estimator = "2step", degrees  = 2)
summary(model_Newey)



# ---------------------------------------
# Simulated data example 5
# Endogenous switching model with 
# heteroscedastic ordered selection 
# mechanism
# ---------------------------------------

# Load required package
library("mnorm")

# ---
# Step 1
# Simulation of the data
# ---

# Set seed for reproducibility
set.seed(123)

# The number of observations
n <- 1000

# Regressors (covariates)
w1 <- runif(n = n, min = -1, max = 1)
w2 <- runif(n = n, min = -1, max = 1)
w3 <- runif(n = n, min = -1, max = 1)

# Random errors
rho_0    <- -0.8
rho_1    <- -0.7
var2_0   <- 0.9
var2_1   <- 1
sd_y_0   <- sqrt(var2_0)
sd_y_1   <- sqrt(var2_1)
cor_y_01 <- 0.7
cov2_01  <- sd_y_0 * sd_y_1 * cor_y_01
cov2_z_0 <- rho_0 * sd_y_0
cov2_z_1 <- rho_1 * sd_y_1
sigma    <- matrix(c(1,        cov2_z_0, cov2_z_1,
                     cov2_z_0, var2_0,   cov2_01,
                     cov2_z_1, cov2_01,  var2_1),
                   nrow = 3)
errors   <- mnorm::rmnorm(n = n, mean = c(0, 0, 0), sigma = sigma)
u        <- errors[, 1]
eps_0    <- errors[, 2]
eps_1    <- errors[, 3]

# Coefficients
gamma     <- c(-1, 2)
gamma_het <- c(0.5, -1)
beta_0    <- c(1, -1, 1)
beta_1    <- c(2, -1.5, 0.5)

# Linear predictor (index) of the ordinal equation
  # mean
li <- gamma[1] * w1 + gamma[2] * w2
  # variance
li_het <- gamma_het[1] * w2 + gamma_het[2] * w3

# Linear predictor (index) of the continuous equation
  # regime 0
li_y_0 <- beta_0[1] + beta_0[2] * w1 + beta_0[3] * w3
  # regime 1
li_y_1 <- beta_1[1] + beta_1[2] * w1 + beta_1[3] * w3

# Latent variable
z_star   <- li + u * exp(li_het)
y_0_star <- li_y_0 + eps_0
y_1_star <- li_y_1 + eps_1

# Cuts
cuts <- c(-1, 0.5, 2)

# Observable ordinal outcome
z                                           <- rep(0, n)
z[(z_star > cuts[1]) & (z_star <= cuts[2])] <- 1
z[(z_star > cuts[2]) & (z_star <= cuts[3])] <- 2
z[z_star > cuts[3]]                         <- 3
table(z)

# Observable continuous outcome such that y' is 
# observable in regime 1 when 'z = 1', 
# observable in regime 0 when 'z <= 1',
# unobservable           when 'z = 0'
y         <- rep(NA, n)
y[z == 0] <- NA
y[z == 1] <- y_0_star[z == 1]
y[z > 1]  <- y_1_star[z > 1]

# Data
data <- data.frame(w1 = w1, w2 = w2, w3 = w3, 
                   z  = z,  y  = y)

# ---
# Step 2
# Estimation of the parameters
# ---

# Assign groups
groups <- matrix(0:3, ncol = 1)
groups2 <- matrix(c(-1, 0, 1, 1), ncol = 1)

# Estimation
model <- msel(list(z ~ w1 + w2 | w2 + w3),
                   list(y ~ w1 + w3),
                   groups = groups, groups2 = groups2,
                   data = data)
summary(model)

# Compare estimates and true values of parameters
  # regression coefficients of the ordinal equation
gamma_est     <- coef(model, type = "coef",     eq = 1)
gamma_het_est <- coef(model, type = "coef_var", eq = 1)
cbind(true = gamma,     estimate = gamma_est)
cbind(true = gamma_het, estimate = gamma_het_est)
  # cuts
cuts_est <- coef(model, type = "cuts", eq = 1)
cbind(true = cuts, estimate = cuts_est)   
  # regression coefficients of the continuous equation
beta_0_test <- coef(model, type = "coef2", eq2 = 1, regime = 0)
beta_1_test <- coef(model, type = "coef2", eq2 = 1, regime = 1)
cbind(true = beta_0, estimate = beta_0_test) 
cbind(true = beta_1, estimate = beta_1_test)
  # variances
var2_0_est <- coef(model, type = "var", eq2 = 1, regime = 0)
var2_1_est <- coef(model, type = "var", eq2 = 1, regime = 1)
cbind(true = c(var2_0, var2_1), estimate = c(var2_0_est, var2_1_est)) 
  # covariances between the random errors
cov2_z_0_est <- coef(model, type = "cov12", eq = 1, eq2 = 1, regime = 0)
cov2_z_1_est <- coef(model, type = "cov12", eq = 1, eq2 = 1, regime = 1)
cbind(true     = c(cov2_z_0,     cov2_z_1), 
      estimate = c(cov2_z_0_est, cov2_z_1_est))

# ---
# Step 3
# Estimation of the expectations and marginal effects
# ---

# New data
newdata <- data.frame(z  = 1,   y  = 1,
                      w1 = 0.1, w2 = 0.2, w3 = 0.3)

# Predict the unconditional expectation of the 
# continuous outcome E(yr | w)
  # regime 0 
predict(model, group = -1, group2 = 0, newdata = newdata)
  # regime 1
predict(model, group = -1, group2 = 1, newdata = newdata)

# Predict the conditional expectations of the continuous outcome
# given condition 'z == 0' for regime 1 i.e., E(y1 | z = 0, w)
predict(model, group = 0, group2 = 1, newdata = newdata)



# ---------------------------------------
# Simulated data example 6
# Endogenous switching model with 
# multivariate heteroscedastic ordered 
# selection mechanism
# ---------------------------------------

# Load required package
library("mnorm")

# ---
# Step 1
# Simulation of the data
# ---

# Set seed for reproducibility
set.seed(123)

# The number of observations
n <- 1000

# Regressors (covariates)
w1 <- runif(n = n, min = -1, max = 1)
w2 <- runif(n = n, min = -1, max = 1)
w3 <- runif(n = n, min = -1, max = 1)
w4 <- runif(n = n, min = -1, max = 1)

# Random errors
rho_z1_z2 <- 0.5
rho_y0_z1 <- 0.6
rho_y0_z2 <- 0.7
rho_y1_z1 <- 0.65
rho_y1_z2 <- 0.75
var20     <- 0.9
var21     <- 1
sd_y0     <- sqrt(var20)
sd_y1     <- sqrt(var21)
cor_y01   <- 0.7
cov201    <- sd_y0 * sd_y1 * cor_y01
cov20_z1  <- rho_y0_z1 * sd_y0
cov21_z1  <- rho_y1_z1 * sd_y1
cov20_z2  <- rho_y0_z2 * sd_y0
cov21_z2  <- rho_y1_z2 * sd_y1
sigma     <- matrix(c(1,         rho_z1_z2,  cov20_z1,  cov21_z1, 
                      rho_z1_z2, 1,          cov20_z2,  cov21_z2,
                      cov20_z1,  cov20_z2,   var20,     cov201,
                      cov21_z1,  cov21_z2,   cov201,    var21),
                    nrow = 4)
errors    <- mnorm::rmnorm(n = n, mean = c(0, 0, 0, 0), sigma = sigma)
u1        <- errors[, 1]
u2        <- errors[, 2]
eps0      <- errors[, 3]
eps1      <- errors[, 4]

# Coefficients
gamma1     <- c(-1, 2)
gamma1_het <- c(0.5, -1)
gamma2     <- c(1, 1)
beta0      <- c(1, -1, 1, -1.2)
beta1      <- c(2, -1.5, 0.5, 1.2)

# Linear index (predictor) of the ordinal equation
  # mean
li1 <- gamma1[1] * w1 + gamma1[2] * w2
li2 <- gamma2[1] * w1 + gamma2[2] * w3
  # variance
li1_het <- gamma1_het[1] * w2 + gamma1_het[2] * w3

# Linear predictor (index) of the continuous equation
  # regime 0
li_y0 <- beta0[1] + beta0[2] * w1 + beta0[3] * w3 + beta0[4] * w4
  # regime 1
li_y1 <- beta1[1] + beta1[2] * w1 + beta1[3] * w3 + beta1[4] * w4

# Latent variables
z1_star <- li1 + u1 * exp(li1_het)
z2_star <- li2 + u2
y0_star <- li_y0 + eps0
y1_star <- li_y1 + eps1

# Cuts
cuts1 <- c(-1, 1)
cuts2 <- c(0)

# Observable ordered outcome
  # first
z1                                               <- rep(0, n)
z1[(z1_star > cuts1[1]) & (z1_star <= cuts1[2])] <- 1
z1[z1_star > cuts1[2]]                           <- 2
  # second
z2                     <- rep(0, n)
z2[z2_star > cuts2[1]] <- 1
table(z1, z2)

# Observable continuous outcome such
# that outcome 'y' is 
# in regime 0  when 'z1 == 1', 
# in regime 1  when 'z1 == 0' or 'z1 == 2',
# unobservable when 'z2 == 0'
y          <- rep(NA, n)
y[z1 == 1] <- y0_star[z1 == 1]
y[z1 != 1] <- y1_star[z1 != 1]
y[z2 == 0] <- NA

# Data
data <- data.frame(w1 = w1, w2 = w2, w3 = w3, w4 = w4,
                   z1 = z1, z2 = z2, y  = y)

# ---
# Step 2
# Estimation of the parameters
# ---

# Assign groups
groups <- matrix(c(0, 0,
                   0, 1,
                   1, 0,
                   1, 1,
                   2, 0,
                   2, 1), 
                 byrow = TRUE, ncol = 2)
groups2 <- matrix(c(-1, 1, -1, 0, -1, 1), ncol = 1)

# Estimation
model <- msel(list(z1 ~ w1 + w2 | w2 + w3,
                   z2 ~ w1 + w3),
              y ~ w1 + w3 + w4,
              groups = groups, groups2 = groups2,
              data   = data)
summary(model)

# Compare estimates and true values of the parameters
  # regression coefficients of the first ordered equation
gamma1_est      <- coef(model, type = "coef",     eq = 1)
gamma1__het_est <- coef(model, type = "coef_var", eq = 1)
cbind(true = gamma1,     estimate = gamma1_est)
cbind(true = gamma1_het, estimate = gamma1__het_est)
  # regression coefficients of the second ordered equation
gamma2_est <- coef(model, type = "coef", eq = 2)
cbind(true = gamma2, estimate = gamma2_est)
  # cuts
cuts1_est <- coef(model, type = "cuts", eq = 1)
cuts2_est <- coef(model, type = "cuts", eq = 2)
cbind(true = cuts1, estimate = cuts1_est)   
cbind(true = cuts2, estimate = cuts2_est)  
  # regression coefficients of the continuous equation
beta0_est <- coef(model, type = "coef2", eq2 = 1, regime = 0)
beta1_est <- coef(model, type = "coef2", eq2 = 1, regime = 1)
cbind(true = beta0, estimate = beta0_est) 
cbind(true = beta1, estimate = beta1_est)
  # variances
var20_est <- coef(model, type = "var", eq2 = 1, regime = 0)
var21_est <- coef(model, type = "var", eq2 = 1, regime = 1)
cbind(true = c(var20, var21), estimate = c(var20_est, var21_est)) 
  # covariances
cov_y0_z1_est <- coef(model, type = "cov12", eq = 1, eq2 = 1, regime = 0)
cov_y0_z2_est <- coef(model, type = "cov12", eq = 2, eq2 = 1, regime = 0)
cov_y1_z1_est <- coef(model, type = "cov12", eq = 1, eq2 = 1, regime = 1)
cov_y1_z2_est <- coef(model, type = "cov12", eq = 2, eq2 = 1, regime = 1)
cbind(true     = c(cov20_z1, cov20_z2), 
      estimate = c(cov_y0_z1_est, cov_y0_z2_est))
cbind(true     = c(cov21_z1, cov21_z2), 
      estimate = c(cov_y1_z1_est,  cov_y1_z2_est))

# ---
# Step 3
# Estimation of the expectations and marginal effects
# ---

# New data
newdata <- data.frame(z1 = 1,   z2 = 1,   y  = 1,
                      w1 = 0.1, w2 = 0.2, w3 = 0.3, w4 = 0.4)

# Predict the unconditional expectation of the continuous outcome
  # regime 0 
predict(model, group = c(-1, -1), group2 = 0, newdata = newdata)
  # regime 1
predict(model, group = c(-1, -1), group2 = 1, newdata = newdata)

# Predict the conditional expectations of the continuous outcome
# E(y1 | z1 = 2, z2 = 1, w)
predict(model, group = c(2, 1), group2 = 1, newdata = newdata)

# Marginal effect of w3 on E(y1 | z1 = 2, z2 = 1, w)
predict(model, group = c(2, 1), group2 = 1, newdata = newdata, me = "w3")

# ---
# Step 4
# Two-step estimation procedure
# ---

# For a comparison reasons let's estimate the model
# via the least squares
model.ls.0 <- lm(y ~ w1 + w3 + w4,
                 data = data[!is.na(data$y) & (data$z1 == 1), ])
model.ls.1 <- lm(y ~ w1 + w3 + w4,
                 data = data[!is.na(data$y) & (data$z1 != 1), ])


# Apply the two-step procedure
model_ts <-  msel(list(z1 ~ w1 + w2 | w2 + w3,
                       z2 ~ w1 + w3),
                  y ~ w1 + w3 + w4,
                  groups    = groups,  groups2 = groups2,
                  estimator = "2step", data    = data)
summary(model_ts)

# Use the two-step procedure with logistic marginal distributions
# that is multivariate generalization of the Lee's method
model_Lee <-  msel(list(z1 ~ w1 + w2 | w2 + w3,
                        z2 ~ w1 + w3),
                   y ~ w1 + w3 + w4,
                   marginal = list(logistic = NULL, logistic = NULL),
                   groups    = groups,  groups2 = groups2,
                   estimator = "2step", data    = data)
                        
# Apply the Newey's method
model_Newey <-  msel(list(z1 ~ w1 + w2 | w2 + w3,
                          z2 ~ w1 + w3),
                     y ~ w1 + w3 + w4,
                     marginal = list(logistic = NULL, logistic = NULL),
                     degrees = c(2, 3), groups = groups, groups2 = groups2,
                     estimator = "2step", data = data)

# Compare accuracy of the methods
  # beta0
tbl0 <- cbind(true    = beta0, 
              ls      = coef(model.ls.0),
              ml      = model$coef2[[1]][1, 1:length(beta0)],
              twostep = model_ts$coef2[[1]][1, 1:length(beta0)],
              Lee     = model_Lee$coef2[[1]][1, 1:length(beta0)],
              Newey   = model_Newey$coef2[[1]][1, 1:length(beta0)])  
print(tbl0)
  # beta1
tbl1 <- cbind(true    = beta1, 
              ls      = coef(model.ls.1),
              ml      = model$coef2[[1]][2, 1:length(beta1)],
              twostep = model_ts$coef2[[1]][2, 1:length(beta1)],
              Lee     = model_Lee$coef2[[1]][2, 1:length(beta1)],
              Newey   = model_Newey$coef2[[1]][2, 1:length(beta1)])  
print(tbl1)



# ---------------------------------------
# Simulated data example 7
# Endogenous multivariate switching model 
# with multivariate heteroscedastic
# ordered selection mechanism
# ---------------------------------------

# Load required package
library("mnorm")

# ---
# Step 1
# Simulation of the data
# ---

# Set seed for reproducibility
set.seed(123)

# The number of observations
n <- 1000

# Regressors (covariates)
w1 <- runif(n = n, min = -1, max = 1)
w2 <- runif(n = n, min = -1, max = 1)
w3 <- runif(n = n, min = -1, max = 1)
w4 <- runif(n = n, min = -1, max = 1)
w5 <- runif(n = n, min = -1, max = 1)

# Random errors
var_y0 <- 0.9
var_y1 <- 1
var_g0 <- 1.1
var_g1 <- 1.2
var_g2 <- 1.3
A      <- rWishart(1, 7, diag(7))[, , 1]
B      <- diag(sqrt(c(1, 1, var_y0, var_y1, 
                      var_g0, var_g1, var_g2)))
sigma  <-  B  %*% cov2cor(A) %*% B
errors <- mnorm::rmnorm(n = n, mean = rep(0, nrow(sigma)), sigma = sigma)
u1     <- errors[, 1]
u2     <- errors[, 2]
eps0_y <- errors[, 3]
eps1_y <- errors[, 4]
eps0_g <- errors[, 5]
eps1_g <- errors[, 6]
eps2_g <- errors[, 7]

# Coefficients
gamma1     <- c(-1, 2)
gamma1_het <- c(0.5, -1)
gamma2     <- c(1, 1)
beta0_y    <- c(1, -1, 1, -1.2)
beta1_y    <- c(2, -1.5, 0.5, 1.2)
beta0_g    <- c(-1, 1, 1, 1)
beta1_g    <- c(1, -1, 1, 1)
beta2_g    <- c(1, 1, -1, 1)

# Linear predictor (index) of the ordinal equation
  # mean
li1 <- gamma1[1] * w1 + gamma1[2] * w2
li2 <- gamma2[1] * w1 + gamma2[2] * w3
  # variance
li1_het <- gamma1_het[1] * w2 + gamma1_het[2] * w3

# Linear predictor (index) of the first continuous equation
  # regime 0
li_y0 <- beta0_y[1] + beta0_y[2] * w1 + beta0_y[3] * w3 + beta0_y[4] * w4
  # regime 1
li_y1 <- beta1_y[1] + beta1_y[2] * w1 + beta1_y[3] * w3 + beta1_y[4] * w4

# Linear predictor (index) of the second continuous equation
  # regime 0
li_g0 <- beta0_g[1] + beta0_g[2] * w2 + beta0_g[3] * w3 + beta0_g[4] * w5
  # regime 1
li_g1 <- beta1_g[1] + beta1_g[2] * w2 + beta1_g[3] * w3 + beta1_g[4] * w5
  # regime 2
li_g2 <- beta2_g[1] + beta2_g[2] * w2 + beta2_g[3] * w3 + beta2_g[4] * w5

# Latent variables
z1_star <- li1 + u1 * exp(li1_het)
z2_star <- li2 + u2
y0_star <- li_y0 + eps0_y
y1_star <- li_y1 + eps1_y
g0_star <- li_g0 + eps0_g
g1_star <- li_g1 + eps1_g
g2_star <- li_g2 + eps2_g

# Cuts
cuts1 <- c(-1, 1)
cuts2 <- c(0)

# Observable ordered outcome
  # first
z1                                               <- rep(0, n)
z1[(z1_star > cuts1[1]) & (z1_star <= cuts1[2])] <- 1
z1[z1_star > cuts1[2]]                           <- 2
  # second
z2                     <- rep(0, n)
z2[z2_star > cuts2[1]] <- 1
table(z1, z2)

# Observable continuous outcome such that outcome 'y' is 
# in regime 0  when 'z1 == 1', 
# in regime 1  when 'z1 == 0' or 'z1 == 2',
# unobservable when 'z2 == 0'
y          <- rep(NA, n)
y[z1 == 1] <- y0_star[z1 == 1]
y[z1 != 1] <- y1_star[z1 != 1]
y[z2 == 0] <- NA

#' # Observable continuous outcome such
# that outcome 'g' is 
# in regime 0 when 'z1 == z2', 
# in regime 1 when 'z1 > z2',
# in regime 2 when 'z1 < z2',
g           <- rep(NA, n)
g[z1 == z2] <- g0_star[z1 == z2]
g[z1 > z2]  <- g1_star[z1 > z2]
g[z1 < z2]  <- g2_star[z1 < z2]

# Data
data <- data.frame(w1 = w1, w2 = w2, w3 = w3, w4 = w4, w5 = w5,
                   z1 = z1, z2 = z2, y  = y,  g  = g)

# ---
# Step 2
# Estimation of the parameters
# ---

# Assign groups
groups <- matrix(c(0, 0,
                   0, 1,
                   1, 0,
                   1, 1,
                   2, 0,
                   2, 1), 
                 byrow = TRUE, ncol = 2)

# Assign groups 2
  # prepare the matrix
groups2 <- matrix(NA, nrow = nrow(groups), ncol = 2)
  # fill the matrix
groups2[groups[, 1] == 1, 1]                        <- 0
groups2[(groups[, 1] == 0) | (groups[, 1] == 2), 1] <- 1
groups2[groups[, 2] == 0, 1]                        <- -1
groups2[groups[, 1] == groups[, 2], 2]              <- 0
groups2[groups[, 1] > groups[, 2], 2]               <- 1
groups2[groups[, 1] < groups[, 2], 2]               <- 2

# The structure of the model
cbind(groups, groups2)

# Estimation
model <- msel(list(z1 ~ w1 + w2 | w2 + w3, z2 ~ w1 + w3),
              list(y ~ w1 + w3 + w4, g ~ w2 + w3 + w5),
              groups = groups, groups2 = groups2, data = data)
summary(model)

# Compare estimates and true values of the parameters
  # regression coefficients of the first ordered equation
gamma1_est     <- coef(model, type = "coef", eq = 1)
gamma1_het_est <- coef(model, type = "coef_var", eq = 1)
cbind(true = gamma1,     estimate = gamma1_est)
cbind(true = gamma1_het, estimate = gamma1_het_est)
  # regression coefficients of the second ordered equation
gamma2_est <- coef(model, type = "coef", eq = 2)
cbind(true = gamma2, estimate = gamma2_est)
  # cuts
cuts1_est <- coef(model, type = "cuts", eq = 1)
cuts2_est <- coef(model, type = "cuts", eq = 2)
cbind(true = cuts1, estimate = cuts1_est)   
cbind(true = cuts2, estimate = cuts2_est)  
  # regression coefficients of the first continuous equation
beta0_y_est <- coef(model, type = "coef2", eq2 = 1, regime = 0)
beta1_y_est <- coef(model, type = "coef2", eq2 = 1, regime = 1)
cbind(true = beta0_y, estimate = beta0_y_est) 
cbind(true = beta1_y, estimate = beta1_y_est)
  # regression coefficients of the second continuous equation
beta0_g_est <- coef(model, type = "coef2", eq2 = 2, regime = 0)
beta1_g_est <- coef(model, type = "coef2", eq2 = 2, regime = 1)
beta2_g_est <- coef(model, type = "coef2", eq2 = 2, regime = 2)
cbind(true = beta0_g, estimate = beta0_g_est) 
cbind(true = beta1_g, estimate = beta1_g_est)
cbind(true = beta2_g, estimate = beta2_g_est)
  # variances of the first continuous equation
var_y0_est <- coef(model, type = "var", eq2 = 1, regime = 0)
var_y1_est <- coef(model, type = "var", eq2 = 1, regime = 1)
cbind(true = c(var_y0, var_y1), estimate = c(var_y0_est, var_y1_est)) 
  # variances of the second continuous equation
var_g0_est <- coef(model, type = "var", eq2 = 2, regime = 0)
var_g1_est <- coef(model, type = "var", eq2 = 2, regime = 1)
var_g2_est <- coef(model, type = "var", eq2 = 2, regime = 2)
cbind(true     = c(var_g0, var_g1, var_g2), 
      estimate = c(var_g0_est, var_g1_est, var_g2_est)) 
  # correlation between the ordinal equations
sigma12_est <- coef(model, type = "cov1", eq = c(1, 2))
cbind(true = c(sigma[1, 2]), estimate = sigma12_est)
  # covariances between the continuous and ordinal equations
cbind(true = sigma[1:2, 3], estimate = model$cov2[[1]][1, ])
cbind(true = sigma[1:2, 4], estimate = model$cov2[[1]][2, ])
cbind(true = sigma[1:2, 5], estimate = model$cov2[[2]][1, ])
cbind(true = sigma[1:2, 6], estimate = model$cov2[[2]][2, ])
cbind(true = sigma[1:2, 7], estimate = model$cov2[[2]][3, ])
  # covariances between the continuous equations
sigma2_est <- coef(model, type = "cov2")[[1]]
cbind(true     = c(sigma[4, 7], sigma[3, 5], sigma[4, 6]), 
      estimate = sigma2_est) 
      
# ---
# Step 3
# Estimation of the expectations and marginal effects
# ---

# New data
newdata <- data.frame(z1 = 1,   z2 = 1,   y  = 1,   g  = 1,
                      w1 = 0.1, w2 = 0.2, w3 = 0.3, w4 = 0.4, w5 = 0.5)

# Predict unconditional expectation of the dependent variable
  # regime 0 for 'y' and regime 1 for 'g' i.e. E(y0 | w), E(g1 | w)
predict(model, group = c(-1, -1), group2 = c(0, 1), newdata = newdata)

# Predict conditional expectations of the dependent variable
# E(y0 | z1 = 2, z2 = 1, w), E(g1 | z1 = 2, z2 = 1, w)
predict(model, group = c(2, 1), group2 = c(0, 1), newdata = newdata)

# Marginal effect of w3 on 
# E(y1 | z1 = 2, z2 = 1, w) and E(g1 | z1 = 2, z2 = 1, w)
predict(model, group   = c(2, 1), group2 = c(0, 1), 
               newdata = newdata, me     = "w3")

# ---
# Step 4
# Two-step estimation procedure
# ---

# Provide manually selectivity terms
model2 <- msel(list(z1 ~ w1 + w2 | w2 + w3, z2 ~ w1 + w3),
               list(y ~ w1 + w3 + w4 + 
                        lambda1 + lambda2 + I(lambda1 * lambda2), 
                    g ~ w2 + w3 + w5 + lambda1 + lambda2),
               groups = groups, groups2   = groups2,
               data   = data,   estimator = "2step")
summary(model2)



# ---------------------------------------
# Simulated data example 8
# Multinomial endogenous switching and 
# selection model (probit)
# ---------------------------------------

# Load required package
library("mnorm")

# ---
# Step 1
# Simulation of the data
# ---

# Set seed for reproducibility
set.seed(123)

# The number of observations
n <- 10000

# Random errors
  # variances and correlations
sd.z2    <- sqrt(0.9)
cor.z    <- 0.3
sd.y0    <- sqrt(2)
cor.z1y0 <- 0.4
cor.z2y0 <- 0.7
sd.y1    <- sqrt(1.8)
cor.z1y1 <- 0.3
cor.z2y1 <- 0.6
cor.y    <- 0.8
  # the covariance matrix
sigma <- matrix(c(
1,                cor.z * sd.z2,            cor.z1y0 * sd.y0,         cor.z1y1 * sd.y1,
cor.z * sd.z2,    sd.z2 ^ 2,                cor.z2y0 * sd.z2 * sd.y0, cor.z2y1 * sd.z2 * sd.y1,
cor.z1y0 * sd.y0, cor.z2y0 * sd.z2 * sd.y0, sd.y0 ^ 2,                cor.y * sd.y0 * sd.y1,
cor.z1y1 * sd.y1, cor.z2y1 * sd.z2 * sd.y1, cor.y * sd.y0 * sd.y1,    sd.y1 ^ 2),
                ncol = 4, byrow = TRUE)
colnames(sigma) <- c("z1", "z2", "y0", "y1")
rownames(sigma) <- colnames(sigma)

# Simulate the random errors
errors <- rmnorm(n, c(0, 0, 0, 0), sigma)
u      <- errors[, 1:2]
eps    <- errors[, 3:4]

# Regressors (covariates)
x1 <- runif(n, -1, 1)
x2 <- runif(n, -1, 1)
x3 <- (x2 + runif(n, -1, 1)) / 2
W  <- cbind(1, x1, x2)
X  <- cbind(1, x1, x3)

# Coefficients
gamma0 <- c(0.1,    1,   1)
gamma1 <- c(0.2, -1.2, 0.8)
beta0  <- c(1,     -3,   4) 
beta1  <- c(1,      4,  -3)

# Linear predictors (indexes)
z1.li <- W %*% gamma0
z2.li <- W %*% gamma1
y0.li <- X %*% beta0
y1.li <- X %*% beta1

# Latent variables
z1.star <- z1.li + u[, 1]
z2.star <- z2.li + u[, 2]
y0.star <- y0.li + eps[, 1]
y1.star <- y1.li + eps[, 2]

# Obvservable variable as a dummy
z1 <- (z1.star > z2.star) & (z1.star > 0)
z2 <- (z2.star > z1.star) & (z2.star > 0)
z3 <- (z1 != 1) & (z2 != 1)

# Observable multinomial variable
z     <- rep(0, n)
z[z1] <- 0
z[z2] <- 1
z[z3] <- 2
table(z)

# Make unobservable values of the continuous outcome
y         <-  rep(NA, n)
y[z == 1] <- y0.star[z == 1]
y[z == 2] <- y1.star[z == 2]

# Data
data <- data.frame(z = z, y = y, x1 = x1, x2 = x2, x3 = x3)
            
# ---
# Step 2
# Estimation of the parameters
# ---

# Define the groups
groups3 <- c(0, 1, 2)
groups2 <- matrix(c(-1, 0, 1), ncol = 1)

# Two-step method
model <- msel(formula3  = z ~ x1 + x2, formula2 = y ~ x1 + x3, 
              groups3   = groups3,     groups2   = groups2,
              data      = data,        estimator = "2step",
              type3     = "probit")
summary(model)

# Compare estimates and true values of parameters
  # regression coefficients of the continuous equation
beta0_est <- coef(model, type = "coef2", eq2 = 1, regime = 0)
beta1_est <- coef(model, type = "coef2", eq2 = 1, regime = 1)
cbind(true = beta0, est = beta0_est[1:length(beta0)])
cbind(true = beta1, est = beta1_est[1:length(beta1)])
  # regression coefficients of the multinomial equations
gamma0_est <- coef(model, type = "coef3", eq3 = 0)
gamma1_est <- coef(model, type = "coef3", eq3 = 1)
cbind(true = gamma0, est = gamma0_est)
cbind(true = gamma1, est = gamma1_est)
  # compare the covariances between
  # z1 and z2
cbind(true = cor.z * sd.z2, 
      est  = coef(model, type = "cov3", eq3 = c(0, 1)))
  # z1 and y0
cbind(true = cor.z1y0 * sd.y0, 
      est  = beta0_est["lambda1_mn"]) 
  # z2 and y0
cbind(true = cor.z2y0 * sd.y0, 
      est  = beta0_est["lambda2_mn"])   
  # z1 and y1
cbind(true = cor.z1y1 * sd.y1, 
      est  = beta1_est["lambda1_mn"]) 
  # z2 and y1
cbind(true = cor.z2y1 * sd.y1, 
      est  = beta1_est["lambda2_mn"]) 
            
# ---
# Step 3
# Predictions and marginal effects
# ---  

# Unconditional expectation E(y1 | w) for every observation in a sample
predict(model, type = "val", group2 = 1, group3 = -1)

# Marginal effect of x1 on conditional expectation E(y0 | z = 1, w) 
# for every observation in a sample
predict(model, type = "val", group2 = 0, group3 = 1, me = "x1")    

# Calculate predictions and marginal effects
# for manually provided observations
# using aforementioned models.
newdata <- data.frame(z  = c(1, 1), 
                      y  = c(1, 1), 
                      x1 = c(0.5, 0.2), 
                      x2 = c(-0.3, 0.8),
                      x3 = c(0.6, -0.7))
                      
# Unconditional expectation E(y0 | w)
predict(model, type = "val", group2 = 0, group3 = -1, newdata = newdata)

# Conditional expectation E(y1 | z=2, w)
predict(model, type = "val", group2 = 1, group3 = 2, newdata = newdata)

# Marginal effect of x2 on E(y0 | z = 1, w)
predict(model, type = "val", group2  = 0, group3 = 1, 
               me   = "x2",  newdata = newdata)
  
# ---
# Step 4
# Multinomial logit selection
# ---                
        
# Two-step method
model2 <- msel(formula3  = z ~ x1 + x2, formula2 = y ~ x1 + x3, 
               groups3   = groups3,     groups2   = groups2,
               data      = data,        estimator = "2step",
               type3     = "logit")
summary(model2)    

# Compare the estimates
beta0_est2 <- coef(model2, type = "coef2", eq2 = 1, regime = 0)[]
beta1_est2 <- coef(model2, type = "coef2", eq2 = 1, regime = 1)

# beta0 coefficients
cbind(true = beta0, probit = beta0_est[1:3], logit = beta0_est2[1:3])
   
# beta1 coefficients
cbind(true = beta1, probit = beta1_est[1:3], logit = beta1_est2[1:3])            

             

# ---------------------------------------
# Simulated data example 9
# Multinomial endogenous switching and 
# selection model (logit)
# ---------------------------------------

# ---
# Step 1
# Simulation of the data
# ---

# Set seed for reproducibility
set.seed(123)

# The number of observations
n <- 10000

# Random errors
u      <- matrix(-log(-log(runif(n * 3))), nrow = n, ncol = 3)
tau0   <- matrix(c(0.6, -0.4, 0.3), ncol = 1)
tau1   <- matrix(c(-0.3, 0.5, 0.2), ncol = 1)
eps0   <- (u - 0.57721656649) %*% tau0 + rnorm(n)
eps1   <- (u - 0.57721656649) %*% tau1 + rnorm(n)

# Regressors (covariates)
x1 <- runif(n = n, min = -1, max = 1)
x2 <- runif(n = n, min = -1, max = 1)
x3 <- runif(n = n, min = -1, max = 1)

# Coefficients
gamma.0 <- c(0.2, -2,  2)
gamma.1 <- c(0.1,  2, -2)
beta.0  <- c(2,    2,  2)
beta.1  <- c(1,   -2,  2)

# Linear predictors (indexes)
z0.li <- gamma.0[1] + gamma.0[2] * x1 + gamma.0[3] * x2
z1.li <- gamma.1[1] + gamma.1[2] * x1 + gamma.1[3] * x2

# Latent variables
z0.star <- z0.li + u[, 1]
z1.star <- z1.li + u[, 2]
z2.star <- u[, 3]
y0.star <- beta.0[1] + beta.0[2] * x1 + beta.0[3] * x3 + eps0
y1.star <- beta.1[1] + beta.1[2] * x1 + beta.1[3] * x3 + eps1

# Observable multinomial variable
z                                            <- rep(2, n)
z[(z0.star > z1.star) & (z0.star > z2.star)] <- 0
z[(z1.star > z0.star) & (z1.star > z2.star)] <- 1
table(z)

# Unobservable values of the continuous outcome
y         <- rep(NA, n)
y[z == 0] <- y0.star[z == 0]
y[z == 1] <- y1.star[z == 1]

# Data
data <- data.frame(x1 = x1, x2 = x2, x3 = x3, z = z, y = y)
            
# ---
# Step 2
# Estimation of the parameters
# ---

# Define the groups
groups3 <- c(0, 1, 2)
groups2 <- c(0, 1, -1)

# Two-step estimator of Dubin-McFadden
model <- msel(formula3  = z ~ x1 + x2, formula2 = y ~ x1 + x3, 
              groups3   = groups3,     groups2   = groups2,
              data      = data,        estimator = "2step",
              type3     = "logit")
summary(model)

# Least squaes estimates (benchmark)
lm0 <- lm(y ~ x1 + x3, data = data[data$z == 0, ])
lm1 <- lm(y ~ x1 + x3, data = data[data$z == 1, ])


# Compare the estimates of beta0
cbind(true = beta.0,
      DMF  = coef(model, type = "coef2", eq2 = 1, regime = 0),
      LS   = coef(lm0))
      
# Compare the estimates of beta1
cbind(true = beta.1,
      DMF  = coef(model, type = "coef2", eq2 = 1, regime = 1),
      LS   = coef(lm1))

             

switchSelection documentation built on Sept. 26, 2024, 5:07 p.m.