# Chance Constrained Radial Models It solves input and output oriented chance constrained radial DEA models under constant (CCR model), variable (BCC model), non-increasing, non-decreasing or generalized returns to scale, based on the model in Cooper et al. (2002). By default, models are solved in a two-stage process (slacks are maximized). We consider $\mathcal{D}=\left\{ \textrm{DMU}_1, \ldots ,\textrm{DMU}_n \right\} $ a set of $n$ DMUs with $m$ stochastic inputs and $s$ stochastic outputs. Matrices $\tilde{X}=(\tilde{x}_{ij})$ and $\tilde{Y}=(\tilde{y}_{rj})$ are the input and output data matrices, respectively, where $\tilde{x}_{ij}$ and $\tilde{y}_{rj}$ represent the $i$-th input and $r$-th output of the $j$-th DMU. Moreover, we denote by $X=(x_{ij})$ and $Y=(y_{rj})$ their expected values. In general, we denote vectors by bold-face letters and they are considered as column vectors unless otherwise stated. The $0$-vector is denoted by $\bm{0}$ and the context determines its dimension. Given $0<\alpha <1$, the first stage program for $\text{DMU}_o$ with constant returns to scale is given by [REMOVE_ME]$$ \min \limits_{\theta, \bm{\lambda}}\quad \theta [REMOVE_ME_2]$$ [REMOVE_ME]$$ \text{s.t.}\quad P\left\{ \left( \theta \tilde{\mathbf{x}}_o-\tilde{X}\bm{\lambda}\right) _i\geq 0\right\} \geq 1-\alpha,\quad i=1,\ldots ,m, [REMOVE_ME_2]$$ [REMOVE_ME]$$ P\left\{ \left( \tilde{Y}\bm{\lambda}-\tilde{\mathbf{y}}_{o}\right) _r\geq 0\right\} \geq 1-\alpha,\quad r=1,\ldots ,s, [REMOVE_ME_2]$$ [REMOVE_ME]$$ \bm{\lambda}\geq \mathbf{0}, [REMOVE_ME_2]$$ where $\bm{\lambda}=(\lambda_1,\ldots,\lambda_n)^\top$, c("$\\tilde{\\mathbf{x}}_o\n$", "$=(\\tilde{x}_{1o},\\ldots,\\tilde{x}_{mo})^\\top$") and c("$\\tilde{\\mathbf{y}}_o=\n$", "$(\\tilde{y}_{1o},\\ldots,\\tilde{y}_{so})^\\top$") are column vectors. Different returns to scale can be easily considered by adding the corresponding constraints: $\mathbf{e}\bm{\lambda}=1$ (VRS), c("$0\\leq \\mathbf{e}\\bm{\\lambda}\n$", "$\\leq 1$") (NIRS), $\mathbf{e}\bm{\lambda}\geq 1$ (NDRS) or c("$L\\leq \\mathbf{e}\n$", "$\\bm{\\lambda}\\leq U$") (GRS), with $0\leq L\leq 1$ and $U\geq 1$, where $\mathbf{e}=(1,\ldots ,1)$ is a row vector. The corresponding second stage program is given by [REMOVE_ME]$$ \max \limits_{\bm{\lambda},\mathbf{s}^-,\mathbf{s}^+}\quad\mathbf{w}^-\mathbf{s}^-+\mathbf{w}^+\mathbf{s}^+ [REMOVE_ME_2]$$ [REMOVE_ME]$$ \text{s.t.}\quad P\left\{ \left( \theta ^*\tilde{\mathbf{x}}_o-\tilde{X}\bm{\lambda}-\mathbf{s}^-\right) _i\geq 0\right\} = 1-\alpha,\quad i=1,\ldots ,m, [REMOVE_ME_2]$$ [REMOVE_ME]$$ P\left\{ \left( \tilde{Y}\bm{\lambda}-\tilde{\mathbf{y}}_{o}-\mathbf{s}^+\right) _r\geq 0\right\} = 1-\alpha,\quad r=1,\ldots ,s, [REMOVE_ME_2]$$ [REMOVE_ME]$$ \bm{\lambda}\geq \mathbf{0},\,\, \mathbf{s}^-\geq \mathbf{0},\,\,\mathbf{s}^+\geq \mathbf{0}, [REMOVE_ME_2]$$ where $\theta ^*$ is the optimal objective function of the first stage program, $\mathbf{s}^-,\mathbf{s}^+$ are column vectors with the slacks, and $\mathbf{w}^-,\mathbf{w}^+$ are positive row vectors with the weights for the slacks. The deterministic equivalents for a multivariate normal distribution of inputs/outputs are given by [REMOVE_ME]$$ \min \limits_{\theta, \bm{\lambda}}\quad \theta [REMOVE_ME_2]$$ [REMOVE_ME]$$ \text{s.t.}\quad X\bm{\lambda}-\Phi ^{-1}(\alpha)\bm{\sigma} ^-(\theta,\bm{\lambda}) \leq \theta \mathbf{x}_o, [REMOVE_ME_2]$$ [REMOVE_ME]$$ Y\bm{\lambda}+\Phi ^{-1}(\alpha)\bm{\sigma} ^+(\bm{\lambda}) \geq \mathbf{y}_o, [REMOVE_ME_2]$$ [REMOVE_ME]$$ \bm{\lambda}\geq \mathbf{0}, [REMOVE_ME_2]$$ and for the second stage, [REMOVE_ME]$$ \max \limits_{\bm{\lambda},\mathbf{s}^-,\mathbf{s}^+}\quad\mathbf{w}^-\mathbf{s}^-+\mathbf{w}^+\mathbf{s}^+ [REMOVE_ME_2]$$ [REMOVE_ME]$$ \text{s.t.}\quad X\bm{\lambda}+\mathbf{s}^--\Phi ^{-1}(\alpha)\bm{\sigma} ^-(\theta ^*,\bm{\lambda}) =\theta ^* \mathbf{x}_o, [REMOVE_ME_2]$$ [REMOVE_ME]$$ Y\bm{\lambda}-\mathbf{s}^++\Phi ^{-1}(\alpha)\bm{\sigma} ^+(\bm{\lambda}) = \mathbf{y}_{o}, [REMOVE_ME_2]$$ [REMOVE_ME]$$ \bm{\lambda}\geq \mathbf{0},\,\, \mathbf{s}^-\geq \mathbf{0},\,\,\mathbf{s}^+\geq \mathbf{0}, [REMOVE_ME_2]$$ where $\Phi $ is the standard normal distribution, and [REMOVE_ME]$$ \displaystyle \left( \sigma ^-_i\left( \theta, \bm{\lambda}\right)\right) ^2 =\sum _{j,q=1}^n\lambda _j\lambda _q\mathrm{Cov}(\tilde{x}_{ij},\tilde{x}_{iq})-2\theta\sum _{j=1}^n\lambda _j\mathrm{Cov}(\tilde{x}_{ij},\tilde{x}_{io}) [REMOVE_ME_2]$$ [REMOVE_ME]$$ +\theta ^2\mathrm{Var}(\tilde{x}_{io}),\quad i=1,\ldots ,m, [REMOVE_ME_2]$$ [REMOVE_ME]$$ \displaystyle \left( \sigma ^+_r\left( \bm{\lambda}\right)\right) ^2 =\sum _{j,q=1}^n\lambda _j\lambda _q\mathrm{Cov}(\tilde{y}_{rj},\tilde{y}_{rq})-2\sum _{j=1}^n\lambda _j\mathrm{Cov}(\tilde{y}_{rj},\tilde{y}_{ro}) [REMOVE_ME_2]$$ [REMOVE_ME]$$ +\mathrm{Var}(\tilde{y}_{ro}),\quad r=1,\ldots ,s. [REMOVE_ME_2]$$ ```r modelstoch_radial(datadea, alpha = 0.05, dmu_eval = NULL, dmu_ref = NULL, orientation = c("io", "oo"), rts = c("crs", "vrs", "nirs", "ndrs", "grs"), L = 1, U = 1, solver = c("alabama", "cccp", "cccp2", "slsqp"), give_X = TRUE, n_attempts_max = 5, maxslack = FALSE, weight_slack_i = 1, weight_slack_o = 1, vtrans_i = NULL, vtrans_o = NULL, compute_target = TRUE, returnqp = FALSE, silent_ud = FALSE, ...) ``` ## Arguments - `datadea`: The data of class `deadata_stoch`, including `n` DMUs, and the expected values of `m` inputs and `s` outputs. - `alpha`: A value for parameter alpha. - `dmu_eval`: A numeric vector containing which DMUs have to be evaluated. If `NULL` (default), all DMUs are considered. - `dmu_ref`: A numeric vector containing which DMUs are the evaluation reference set. If `NULL` (default), all DMUs are considered. - `orientation`: A string, equal to "io" (input oriented) or "oo" (output oriented). - `rts`: A string, determining the type of returns to scale, equal to "crs" (constant), "vrs" (variable), "nirs" (non-increasing), "ndrs" (non-decreasing) or "grs" (generalized). - `L`: Lower bound for the generalized returns to scale (grs). - `U`: Upper bound for the generalized returns to scale (grs). - `solver`: Character string with the name of the solver used by function `solvecop` from package `optiSolve`. - `give_X`: Logical. If it is `TRUE`, it uses an initial vector (given by the evaluated DMU) for the solver, except for "cccp". If it is `FALSE`, the initial vector is given internally by the solver and it is usually randomly generated. - `n_attempts_max`: A value with the maximum number of attempts if the solver does not converge. Each attempt uses a different initial vector. - `maxslack`: Logical. If it is `TRUE`, it computes the max slack solution. - `weight_slack_i`: A value, vector of length `m`, or matrix `m` x `ne` (where `ne` is the length of `dmu_eval`) with the weights of the input slacks for the max slack solution. - `weight_slack_o`: A value, vector of length `s`, or matrix `s` x `ne` (where `ne` is the length of `dmu_eval`) with the weights of the output slacks for the max slack solution. - `vtrans_i`: Numeric vector of translation for undesirable inputs. If `vtrans_i[i]` is `NA`, then it applies the "max + 1" translation to the i-th undesirable input. If `vtrans_i` is a constant, then it applies the same translation to all undesirable inputs. If `vtrans_i` is `NULL`, then it applies the "max + 1" translation to all undesirable inputs. - `vtrans_o`: Numeric vector of translation for undesirable outputs, analogous to `vtrans_i`, but applied to outputs. - `compute_target`: Logical. If it is `TRUE`, it computes targets of the max slack solution. - `returnqp`: Logical. If it is `TRUE`, it returns the quadratic problems (objective function and constraints) of stage 1. - `silent_ud`: Logical, to avoid warnings related with undesirable variables. - `...`: Other parameters, like the initial vector `X`, to be passed to the solver. ## Returns A list with the results for the evaluated DMUs and other parameters for reproducibility. ## Description It solves input and output oriented chance constrained radial DEA models under constant (CCR model), variable (BCC model), non-increasing, non-decreasing or generalized returns to scale, based on the model in Cooper et al. (2002). By default, models are solved in a two-stage process (slacks are maximized). We consider $\mathcal{D}=\left\{ \textrm{DMU}_1, \ldots ,\textrm{DMU}_n \right\} $ a set of $n$ DMUs with $m$ stochastic inputs and $s$ stochastic outputs. Matrices $\tilde{X}=(\tilde{x}_{ij})$ and $\tilde{Y}=(\tilde{y}_{rj})$ are the input and output data matrices, respectively, where $\tilde{x}_{ij}$ and $\tilde{y}_{rj}$ represent the $i$-th input and $r$-th output of the $j$-th DMU. Moreover, we denote by $X=(x_{ij})$ and $Y=(y_{rj})$ their expected values. In general, we denote vectors by bold-face letters and they are considered as column vectors unless otherwise stated. The $0$-vector is denoted by $\bm{0}$ and the context determines its dimension. Given $0<\alpha <1$, the first stage program for $\text{DMU}_o$ with constant returns to scale is given by $$ \min \limits_{\theta, \bm{\lambda}}\quad \theta $$ $$ \text{s.t.}\quad P\left\{ \left( \theta \tilde{\mathbf{x}}_o-\tilde{X}\bm{\lambda}\right) _i\geq 0\right\} \geq 1-\alpha,\quad i=1,\ldots ,m, $$ $$ P\left\{ \left( \tilde{Y}\bm{\lambda}-\tilde{\mathbf{y}}_{o}\right) _r\geq 0\right\} \geq 1-\alpha,\quad r=1,\ldots ,s, $$ $$ \bm{\lambda}\geq \mathbf{0}, $$ where $\bm{\lambda}=(\lambda_1,\ldots,\lambda_n)^\top$, c("$\\tilde{\\mathbf{x}}_o\n$", "$=(\\tilde{x}_{1o},\\ldots,\\tilde{x}_{mo})^\\top$") and c("$\\tilde{\\mathbf{y}}_o=\n$", "$(\\tilde{y}_{1o},\\ldots,\\tilde{y}_{so})^\\top$") are column vectors. Different returns to scale can be easily considered by adding the corresponding constraints: $\mathbf{e}\bm{\lambda}=1$ (VRS), c("$0\\leq \\mathbf{e}\\bm{\\lambda}\n$", "$\\leq 1$") (NIRS), $\mathbf{e}\bm{\lambda}\geq 1$ (NDRS) or c("$L\\leq \\mathbf{e}\n$", "$\\bm{\\lambda}\\leq U$") (GRS), with $0\leq L\leq 1$ and $U\geq 1$, where $\mathbf{e}=(1,\ldots ,1)$ is a row vector. The corresponding second stage program is given by $$ \max \limits_{\bm{\lambda},\mathbf{s}^-,\mathbf{s}^+}\quad\mathbf{w}^-\mathbf{s}^-+\mathbf{w}^+\mathbf{s}^+ $$ $$ \text{s.t.}\quad P\left\{ \left( \theta ^*\tilde{\mathbf{x}}_o-\tilde{X}\bm{\lambda}-\mathbf{s}^-\right) _i\geq 0\right\} = 1-\alpha,\quad i=1,\ldots ,m, $$ $$ P\left\{ \left( \tilde{Y}\bm{\lambda}-\tilde{\mathbf{y}}_{o}-\mathbf{s}^+\right) _r\geq 0\right\} = 1-\alpha,\quad r=1,\ldots ,s, $$ $$ \bm{\lambda}\geq \mathbf{0},\,\, \mathbf{s}^-\geq \mathbf{0},\,\,\mathbf{s}^+\geq \mathbf{0}, $$ where $\theta ^*$ is the optimal objective function of the first stage program, $\mathbf{s}^-,\mathbf{s}^+$ are column vectors with the slacks, and $\mathbf{w}^-,\mathbf{w}^+$ are positive row vectors with the weights for the slacks. The deterministic equivalents for a multivariate normal distribution of inputs/outputs are given by $$ \min \limits_{\theta, \bm{\lambda}}\quad \theta $$ $$ \text{s.t.}\quad X\bm{\lambda}-\Phi ^{-1}(\alpha)\bm{\sigma} ^-(\theta,\bm{\lambda}) \leq \theta \mathbf{x}_o, $$ $$ Y\bm{\lambda}+\Phi ^{-1}(\alpha)\bm{\sigma} ^+(\bm{\lambda}) \geq \mathbf{y}_o, $$ $$ \bm{\lambda}\geq \mathbf{0}, $$ and for the second stage, $$ \max \limits_{\bm{\lambda},\mathbf{s}^-,\mathbf{s}^+}\quad\mathbf{w}^-\mathbf{s}^-+\mathbf{w}^+\mathbf{s}^+ $$ $$ \text{s.t.}\quad X\bm{\lambda}+\mathbf{s}^--\Phi ^{-1}(\alpha)\bm{\sigma} ^-(\theta ^*,\bm{\lambda}) =\theta ^* \mathbf{x}_o, $$ $$ Y\bm{\lambda}-\mathbf{s}^++\Phi ^{-1}(\alpha)\bm{\sigma} ^+(\bm{\lambda}) = \mathbf{y}_{o}, $$ $$ \bm{\lambda}\geq \mathbf{0},\,\, \mathbf{s}^-\geq \mathbf{0},\,\,\mathbf{s}^+\geq \mathbf{0}, $$ where $\Phi $ is the standard normal distribution, and $$ \displaystyle \left( \sigma ^-_i\left( \theta, \bm{\lambda}\right)\right) ^2 =\sum _{j,q=1}^n\lambda _j\lambda _q\mathrm{Cov}(\tilde{x}_{ij},\tilde{x}_{iq})-2\theta\sum _{j=1}^n\lambda _j\mathrm{Cov}(\tilde{x}_{ij},\tilde{x}_{io}) $$ $$ +\theta ^2\mathrm{Var}(\tilde{x}_{io}),\quad i=1,\ldots ,m, $$ $$ \displaystyle \left( \sigma ^+_r\left( \bm{\lambda}\right)\right) ^2 =\sum _{j,q=1}^n\lambda _j\lambda _q\mathrm{Cov}(\tilde{y}_{rj},\tilde{y}_{rq})-2\sum _{j=1}^n\lambda _j\mathrm{Cov}(\tilde{y}_{rj},\tilde{y}_{ro}) $$ $$ +\mathrm{Var}(\tilde{y}_{ro}),\quad r=1,\ldots ,s. $$ ## Examples ```r # Example 1. library(deaR) data("Coll_Blasco_2006") ni <- 2 # number of inputs no <- 2 # number of outputs data_example <- make_deadata(datadea = Coll_Blasco_2006, ni = ni, no = no) nd <- length(data_example$dmunames) # number of DMUs var_input <- matrix(1, nrow = ni, ncol = nd) var_output <- matrix(1, nrow = no, ncol = nd) data_stoch <- make_deadata_stoch(datadea = data_example, var_input = var_input, var_output = var_output) # Evaluate the sixth DMU Collstoch <- modelstoch_radial(data_stoch, dmu_eval = 6) efficiencies(Collstoch) # Example 2. Deterministic data with one stochastic input. library(deaR) dmunames <- c("A", "B", "C") nd <- length(dmunames) # Number of DMUs inputnames <- c("Input_1", "Input_2") ni <- length(inputnames) # Number of Inputs outputnames <- c("Output_1", "Output_2", "Output_3") no <- length(outputnames) # Number of Outputs X <- matrix(c(5, 14, 8, 15, 7, 12), nrow = ni, ncol = nd, dimnames = list(inputnames, dmunames)) Y <- matrix(c(9, 4, 16, 5, 7, 10, 4, 9, 13), nrow = no, ncol = nd, dimnames = list(outputnames, dmunames)) datadea <- make_deadata(inputs = X, outputs = Y) covX <- array(0, dim = c(2, 3, 3)) # The 2nd input is stochastic. # Since the corresponding 3x3 covariances matrix is symmetric, only values # above the diagonal are necessary. covX[2, 1, ] <- c(1.4, 0.9, 0.6) covX[2, 2, 2:3] <- c(1.5, 0.7) covX[2, 3, 3] <- 1.2 # Alternatively (note that values below the diagonal are ignored). covX[2, , ] <- matrix(c(1.4, 0.9, 0.6, 0, 1.5, 0.7, 0, 0, 1.2), nrow = 3, byrow = TRUE) datadea_stoch <- make_deadata_stoch(datadea, cov_input = covX) alpha <- 0.025 res <- modelstoch_radial(datadea_stoch, alpha = alpha, rts = "vrs") efficiencies(res) # Example 3. Replication of results in Land et al. (1993) library(deaR) data("PFT1981") # Selecting DMUs in Program Follow Through (PFT) PFT <- PFT1981[1:49, ] PFT <- make_deadata(PFT, inputs = 2:6, outputs = 7:9) c <- 0.5 var_output <- matrix(c^2, nrow = 3, ncol = 49) PFT_stoch <- make_deadata_stoch(datadea = PFT, var_output = var_output) # Evaluate the second DMU res <- modelstoch_radial(PFT_stoch, dmu_eval = 2) efficiencies(res) ``` ## References Cooper, W.W.; Deng, H.; Huang, Z.; Li, S.X. (2002). “Chance constrained programming approaches to technical efficiencies and inefficiencies in stochastic data envelopment analysis", Journal of the Operational Research Society, 53:12, 1347-1356. tools:::Rd_expr_doi("10.1057/palgrave.jors.2601433") Land, K.C; Lovell, C.A.K.; Thore, S. (1993). "Chance-constrained data envelopment analysis", Managerial and Decision Economics, Vol. 14, No. 6, 541-554. ## Author(s) Vicente Bolós (vicente.bolos@uv.es ). **Department of Business Mathematics** Rafael Benítez (rafael.suarez@uv.es ). **Department of Business Mathematics** Vicente Coll-Serrano (vicente.coll@uv.es ). **Quantitative Methods for Measuring Culture (MC2). Applied Economics.** University of Valencia (Spain)

modelstoch_radial function