# Multivariate and multinomial sample selection and endogenous switching models with multiple outcomes. 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. ```r 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. ## Returns 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 c("$\\widehat{Cov}(\\varepsilon_{r^{(v)}i}^{(v)}, \n$", "$\\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. ## 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}^{*}), \quadu_{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}c_{(J-1)m_{j-1}}\\\end{cases}, $$ $$ z_{i}=(z_{1i},...,z_{Ji})^{T},\quadu_{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 multinomialequations 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 orderedselection 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 multivariateordered 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}),\quadr\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},\qquadb_{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 multinomialselection 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). ## Examples ```r # --------------------------------------- # 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)) ``` ## 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.

msel function