Stability and bifurcation analysis of $SIQR$ for the COVID-19 epidemic model with time delay

1.   Introduction

Since inception, Data Envelopment Analysis (DEA) developed by Charnes et al. (1978), has been widely used to evaluate the performance of a set of Decision Making Units (DMUs), especially those with multiple inputs and multiple outputs (e.g., Cook et.al, 2009; Cooper et al., 2007; Liang et al., 2008; Thanassoulis et al., 2011 and Wu and Liang, 2010). In the original DEA models, all DMUs are treated as a black box and the internal structure is always ignored; see, e.g., Lewis and Sexton (2004) for more details. However, there are cases where DMUs may consist of two-stage network structures with intermediate measures. In other words, DMUs under evaluation share a common feature found in many two-stage network structures, namely the outputs from the first stage become the inputs to the second stage. This is referred to intermediate measure. For instance, Seiford and Zhu (1999) developed the profitability and marketability of US commercial banks based on a two-stage Network Data Envelopment Analysis (NDEA) model. Kao and Hwang (2008) studied the efficiency of 24 insurance firms based on a two-stage model where operating and insurance expenses are applied to generate premiums in the first stage and the underwriting and investment profits are produced in the second stage by employing the intermediate premiums as a resource. Liang et al. (2008) developed cooperative and non-cooperative models to address the tension between two stages caused by the intermediate measures.

The DEA models are classified into four series: original DEA models; efficiency decomposition models; NDEA models; and game-theoretic models (e.g., Cook et al., 2010). Also, Cook and Zhu (2014) revised various DEA models for measuring efficiency in the two-stage network structures. Li et al. (2018) presented an extended model of Despotis et al.'s model (2016) to generate a Pareto solution and identified the leader stage of a two-stage DEA model. They displayed that the optimal solution for the developed model is also a leader-follower solution and that the global optimal solution can be specified by comparing the efficiency scores difference for the upper and lower bounds of the two stages. In recent years, many researchers also studied two-stage NDEA models with shared resources; see, e.g., Yu and Shi (2014), Moreno et al. (2015), Wu et al. (2016 a, 2016 b), Ang and Chen (2016), Guo et al. (2017), Li et al. (2017 c) and Izadikhah et al. (2018).

In all the above mentioned DEA models, all inputs and outputs parameters are assumed to be exact and the effect of uncertainty is ignored. Research showed that a small perturbation in the problem data can lead to a serious variation in ranking. Robust optimization is a widely used approach to deal with the data uncertainty (e.g., Soyster, 1973; Mulvey et al., 1995; Despotis and Smirlis, 2002). This approach firstly takes the percentage of the perturbation in the data into consideration and then obtains the robust efficiency (e.g., Ben-Tal and Nemirovski, 1999). In the realm of robust DEA, Sadjadi and Omrani (2008) considered uncertainty in output parameters for the performance evaluation of electricity distribution firms. They showed that models presented by Bertsimas and Sim(2003, 2004) and Bertsimas and Thiele (2006) were easier and more practical than the robust DEA model introduced by Ben-Tal and Nemirovski (2000). Salahi et al. (2016) presented the robust counterpart of the CCR model in the envelopment form and showed that it is the same as the optimistic robust counterpart of the multiplier form of the CCR model. Then they computed robust solutions for common set of weights under interval uncertainties using robust efficiency scores of units considering as ideal solutions. In another research, Wu et al. (2017) transformed a robust DEA optimization model into a second-order cone equivalent to immunize against output perturbation in an uncertainty set. As well, in a recent study, Salahi et al. (2018) proposed equivalent formulizations of the robust Russell measure model and its enhanced model for interval and ellipsoidal uncertainties in their best- and worst-cases. The authors indicated that the built formulizations remain convex for both best- and worst-cases under interval uncertainty as well as worst-case with ellipsoidal uncertainty. The effect of uncertainty also is studied in two-stage models. Kao and Liu (2011) developed a two-stage DEA model from deterministic to uncertain situations, where the observations are represented by fuzzy numbers. In their study, the extension principle is utilized to develop a pair of two-level mathematical program to calculate the lower and upper bounds of the α-cut of the fuzzy efficiency. As well, Liu (2014) proposed a fuzzy two-stage NDEA model with assurance region. A new two-stage Stackelberg fuzzy DEA model with utilizing the Monte Carlo simulation is proposed by Tavana and Khalili-Damghani (2014), that discriminately ranked the efficiency scores in each stage for branches of a commercial bank. In another study, Alimohammadi (2016) proposed a robust two-stage DEA model in order to evaluate the efficiency of the electrical networks under uncertainty. Hatefi et al. (2016) applied an integrated forward-reverse logistics network with hybrid facilities under uncertainty in which the impact of random facility disruptions is relieved. Zhou et al. (2017) suggested stochastic NDEA models for two-stage systems under the centralized control organization mechanism and showed that ignoring the internal structure of the whole process is not convenient. Also, Huang et al. (2018) developed a NDEA model to Copula-Based network stochastic frontier analysis on evidence from the U.S. commercial banks in 2009 and evaluated technical efficiencies of the stochastic production and cost frontiers. Their dynamic model recognized multiple banks production processes that was independent of whether, more or less deposits are consistent with higher bank efficiency.

Another approach to deal with uncertainty is the p-robust approach. It was first introduced to deal with uncertainty in facility layout by Kouvelis et al. (1992) and used subsequently in a network design problem; see, e.g., Gutierrez et al. (1996). Next the p-robustness opinion by Mo and Horrison (2005) was first used in a supply-chain network design to show that the relative regret in each scenario must not be more than constant p. After that, Snyder and Daskin (2006) introduced a new approach for optimizing under uncertainty known as a stochastic p-robust optimization approach. In their approach, the objective function minimizes the expected costs, while the p-robustness condition is incorporated into the model as a constraint and the required decisions are partitioned into two stages. It is a combination of the robust optimization and traditional stochastic approaches, each of which has some disadvantages to cope the uncertainty. The stochastic models search to minimize the total expected cost among all scenarios. The optimal solution obtained using it may be very good for some scenarios but very poor for the others. The robust approach usually seeks min-max regret solutions that appear effective no matter which scenario is realized.

In this paper, we utilize stochastic p-robust approach for the Stackelberg (leader-follower) and centralized game-theoretic DEA models of Liang et al. (2008) to achieve robustness against the existing uncertainty. The developed models are based on the discrete robust and stochastic optimization approaches that apply probabilistic scenarios to obtain the effect of imprecise input and output parameters and calculate the efficiency scores for DMUs. It can overcome the difficulties of models under uncertain parameters and give reliable answers. In fact, the aim of this study is to create a robust system for NDEA models through comparison of game-theoretic models in which large reductions in regret are possible with little increases in the expected efficiency of branches. The proposed models are very sensitive to the parameters change. One of the main advantages of these models is that it enables decision makers to make a trade-off between the expected value and the regret value of the efficiency of DMUs under probable scenarios and also, a trade-off between robustness of the solution and the robustness of the model. Also, the uncertainty under different scenarios is formulated with a specified probability for input and output data instead of using point estimates. Like the fuzzy approach, this approach does not require any membership functions and the α-cut variables (i.e., minimizing the computational efforts). Also, it does not use a Monte Carlo simulation procedure to discriminately rank the efficient DUMs. Therefore, we can obtain the scenario ranking based on the p-values on the mind.

The rest of the paper is organized as follows. The next section reviews a short preliminary to the game-theoretic models utilized in two-stage DEA models. Section 3 presents the stochastic p-robust approach for centralized and non-centralized models. In Section 4, we apply the proposed models to a real numerical example to demonstrate its efficiency. Sensitivity analysis is provided in Section 5 and at the end, concluding remarks and some directions for future research are given in the last section.

2.   Preliminaries

2.1.   Structure of two-stage NDEA models

In this section, we briefly review the game-theoretic models used in an original two-stage process. Liang et al. (2006) employed the notions of the Stackelberg game (or leader-follower/the seller and buyer) and the cooperative game to extend models for measuring performance in supply chain settings. In their model, the second stage includes both the outputs from the first stage and some additional inputs. Since some of the inputs to the second stage were not from the first stage, one of the DEA models was nonlinear. Although the model was nonlinear, it could be solved as parametric linear programming problems, and the best solution could be found using a heuristic technique.

Liang et al. (2008) applied the model of Liang et al. (2006) for two-stage processes (which includes only intermediate measures joining the to two stages, as illustrated in Figure 1) and then used the overall efficiency definition of Kao and Hwang (2008) to get linear DEA models. In fact, Liang et al. (2008) presented alternative models to address the conflict between stages caused by the intermediate measures, and at the same time provide efficiency scores for both individual stages and the overall process. Their non-cooperative and cooperative models displayed that the overall efficiency of the two-stage process was the product of efficiencies of the two stages.

In general, the two models derived from game-theoretic axioms are the cooperative (centralized) and the non-cooperative (non-centralized) game models; see, e.g., Cook et al. (2010). The non-centralized model considers the two stages as players in a game and pursues a leader-follower model. This model often mentioned to as a Stackelberg game, requires selecting one of the two stages as the leader and then obtaining multipliers for the inputs and outputs that produce the best feasible score for that stage. The efficiency score for the other stage, namely the follower is then attained by detecting the best potential weights for its inputs and outputs, but with the limitation that the score of the leader is not compromised. The centralized game model gains the best overall efficiency score for the two stages incorporated. Consider an original two-stage process as shown in Figure 1. Suppose we have n DMUs, and each DMU_j (j = 1, …, n) uses m inputs x_ij ^s, (i = 1, …, m) in the first stage to produce D outputs z_dj^s , (d = 1, …, D) under scenario s = {1, …, S}. Then, these D outputs become the inputs to the second stage and will be referred to as intermediate measures. The outputs from the second stage are y_rj ^s, (r = 1, …, t) under scenario s. Therefore, x_ij^s, z_dj^s and y _rj^s show the i ^th input, the d^th intermediate output and the r^th output of DMUj based on the s^th scenario.

For DMUj we denote the efficiency scores of the first stage as e_j^1s and the second as e_j^2s based on the s^th scenario. Based upon the radial constant returns to scale (CRS) DEA model of Charnes et al. (1978), we have the following input-oriented DEA models for each stage:

In models (1), v_i, w_d¹, w_d², and u_r are non-negative decision variables of inputs, intermediate measures and outputs, respectively. Besides, x_ij^s, z_dj^s and y_rj^s show uncertain parameters related to the i^th input, the d^th intermediate measure and the r^th output of DMUj in terms of the s^th scenario, respectively. It is noteworthy that, here, we assume the weights related to the intermediate parameters are the same, i.e., w_d¹ = w_d² regardless of whether they are considered as inputs or outputs; see, e.g., Kao and Hwang (2008) and Liang et al. (2008) for more details on the properties of models. As in Seiford and Zhu (1999), we can use two separate DEA analysis for the two stages, so that if the first stage is efficient and the second stage is not, then the second stage improves its performance by decreasing the inputs z_dj^s, that the reduced z_dj^s may make the first stage inefficient. However, to compute the overall efficiency, Chen et al. (2009) proposed an additive or arithmetic mean methodology for incorporating the two stages, as opposed to the geometric-type model of Kao and Hwang (2008).

Definition 2.1 If e_j^1s and e_j^2s are the efficiency scores of the first and the second stages under the s^th scenario, respectively then the efficiency of the overall two-stage process either is $\frac{1}{2}\left( e_{j}^{1s}+e_{j}^{2s} \right)\text{ or }e_{j}^{1s}\times e_{j}^{2s}$.

Definition 2.2 If e_j^1s≤1 and e_j^2s≤1 are the efficiency measures of the first and the second stages under the s^th scenario, respectively in the input oriented model, then the two-stage process is efficient if and only if e_j^1s = e_j^2s = 1.

It is noted that, if $e_{j}^{s}=\sum\nolimits_{r=1}^{t}{{}}{{u}_{r}}y_{rj}^{s}/\sum\nolimits_{i=1}^{m}{{}}{{v}_{i}}x_{ij}^{S}$ is the two-stage overall efficiency score under the s^th scenario, then at optimality, in the second stage we have $e_{0}^{2s*}=\sum\nolimits_{r=1}^{t}{{}}u_{r}^{*}y_{ro}^{s}/e_{0}^{1s}\ \ \text{ and }\ \ \sum\nolimits_{i=1}^{m}{{}}v_{i}^{*}x_{io}^{s}=1$, under the s^th scenario. On the other hand, $e_{0}^{1s*}=\sum\nolimits_{d=1}^{D}{{}}w_{d}^{1*}z_{do}^{s}$ therefore, $e_{0}^{1s*}\times e_{0}^{2s*}=\sum\nolimits_{r=1}^{t}{{}}u_{r}^{*}y_{ro}^{s}/\sum\nolimits_{i=1}^{m}{{}}v_{i}^{*}x_{io}^{s}$. As a result, $e_{0}^{1s*}\times e_{0}^{2s*}=\sum\nolimits_{r=1}^{t}{{}}u_{r}^{*}y_{ro}^{s}$, and at optimality our proposed models give $e_{j}^{S}=e_{j}^{1s}\times e_{j}^{2s}$.

2.2.   Non-centralized model

As mentioned above, one form of a non-centralized game is the leader-follower model. In the game theory literature, the leader-follower model is also referred to as the Stackelberg model. If we consider the first stage as the leader, then the first stage's performance is more significant, and the efficiency of the second stage (follower) is computed subject to that the leader's efficiency stays constant. By accepting that the first and second stages are corresponding with the leader and the follower, respectively we compute the efficiency for the leader model under the s^th scenario, applying the CCR model (2) as follows:

For computing the second stage's efficiency, namely e₀^2s* the model can be expressed as:

At last, in non-centralized model we state the overall efficiency with different scenarios as following. If e₀^1s* and e₀^2s* are the optimal efficiency scores of models (2) and (3), respectively then the overall efficiency is equal to $e_{0}^{1s*}\times e_{0}^{2s*}=\sum\nolimits_{r=1}^{t}{{}}u_{r}^{*}y_{ro}^{S}$. Accordingly, this illustrates that the overall efficiency in non-centralized model in the s^th scenario is equal to the product of the efficiencies of individual stages.

2.3.   Centralized model

Another model for evaluating the efficiency of the two-stage process is the centralized model that specifies a set of optimal weights on the intermediate factors in order to maximize the overall or agammaegate efficiency score.

Remark1. By assuming that in models (1) w_d¹= w_d², it can be seen that $e_{0}^{1s}\times e_{0}^{2s}=\sum\nolimits_{r=1}^{t}{{}}{{u}_{r}}y_{ro}^{s}/\sum\nolimits_{i=1}^{m}{{}}{{v}_{i}}x_{io}^{s}$. As a result, instead of maximizing the average of e₀^1s and e₀^2s under the s^th scenario, we maximize the product of e₀^1s and e₀^2s of individual stages. Therefore, in the centralized model by letting w_d¹ = w_d², the efficiency scores of both stages are computed simultaneously. As mentioned, we write the centralized model by replacing w_d¹ = w_d² in the models (1), as follows:

It is noted that the e₀^{c-s 1} in model (4) is the optimal efficiency value for the centralized model. Model (4) can be transformed into the linear model (5) as

¹ "c" stands for the centralized model.

Model (5) gives the overall efficiency of the two-stage process. Nonetheless, it is obvious that, if there is uncertainty in data set, models (2), (3) or (5) might be infeasible at optimal solution of nominal problem. Thus, it is essential to choose an alternative model such that a small variation in data set cannot change the rankings. To tackle this case, we apply the stochastic p-robust optimization approach of Snyder and Daskin (2006) that will be illustrated in the next section.

3.   Stochastic p-robust two-stage NDEA models

In this section, first we describe the p-robust concept, and afterwards present the mathematical formulation of two-stage NDEA models under uncertainty based on the stochastic p-robust approach.

3.1.   P-robust concept

Suppose S be a set of scenarios, and let P_s be a deterministic (i.e., single scenario) maximization problem for scenario index s, so that there is a different problem P_s for each scenario s $\epsilon$ S. Let Z_s^* > 0 be the optimal objective value for P_s. Also, let X be the feasible vector in terms of the weights of outputs and inputs, respectively, and Z_s(X) be the objective value of problem P_s under solution X. Then X is called p-robust (p≥0 is a non-negative constant) if for all s $\epsilon$ S the following inequality holds:

The left hand side of formula (6) shows the relative regret under s^th scenario and p≥0 is a parameter that shows the robustness level between different values of each scenario. The relative regret for each scenario is limited by it. We can write inequality (6) as follows:

Finally, in order to control the relative regret related to the scenarios, the p-robust restrictions are added to the model.

Definition 3.1 A DMU_j under evaluation under the different scenarios is stochastic p-robust efficient if its optimal objective function value is one.

3.2.   Stochastic p-robust two-stage NDEA models

As mentioned before, in two-stage NDEA models each DMU is composed of two sub-DMUs in sequences, and the intermediate products of the sub-DMU in the first stage are used as input by the sub-DMU in the second stage. Nonetheless, to cope with uncertainty situations the precise models cannot lead to correct results. In fact, uncertainty can change final results and units rankings. Therefore, the original two-stage NDEA models must be robust upon uncertainty. In this case, we propose two stochastic p-robust two-stage NDEA models to cope with this theme.

3.2.1.   Stochastic-robust non-centralized model

Here we propose stochastic p-robust Stackelberg game models or the leader-follower versions of uncertain DEA models (2) and (3). At first, we formulate stochastic p-robust model for the first stage (i.e., the leader problem) as follows:

The objective function of model (8) computes the expected efficiency value of DMUs according to the data from each scenario in the leader stage. The first constraints represent the p-robust restrictions. Moreover, as retaining the leader's efficiency fixed, the stochastic p-robust model for the second stage (i.e. the follower problem) for all scenarios can be formulated as follows:

The objective function of model (9) maximizes the expected efficiency value of DMUs according to data from each scenario in the follower stage. In the objective function, q_s is the probability that scenario s happens (it is unclear which scenario will happen in the future, in other words, there is no information about the probability of chance of each scenario). In this model, the uncertainty in the parameters is defined by discrete scenarios. The objective function maximizes expected efficiency value of DMUs according to the data from each scenario. The first constraints show the p-robust restrictions and constraints (9.f) measure the efficiency value of the first stage (i.e., leader stage) under scenarios. Moreover, the relation between the overall efficiency and the efficiency score of each stage, i.e., f₀^s*, f₀^1s* and f₀^2s*, respectively are illustrated in the sequel.

Definition 3.2.1 DMU_j under evaluation in the s^th scenario is overall efficient if and only f_j^s = 1, j = 1, …, n.

Definition 3.2.2 The k^th stage sub-DMU_j, under the s^th scenario is efficient if f_j^ks = 1, j = 1, …, n and k = 1, 2.

Proposition 3.2.1 DMU₀ in the s^th scenario is overall efficient if and only if both stages in the s^th scenario is efficient.

Proof. According to Definition 3.2.1, if DMU₀ is overall efficient, then f₀^s = 1. Since f₀^s = f₀^1s ×f₀^2s, and f₀^1s ≤1, f₀^2s ≤1, therefore, the divisional efficiencies f₀^1s and f₀^2s must satisfy f₀^1s = f₀^2s = 1. On the contrary, if both stages are efficient, i.e., f₀^1s = f₀^2s = 1, then since f₀^s = f₀^1s ×f₀^2s, so the efficiency value of f₀^s must be equal to one. Hence, based on to Definition 3.2.1, DMU₀ must be overall efficient.

3.2.2.   Stochastic-robust centralized model

The stochastic p-robust centralized model (5) under uncertainty is as follows:

In model (10), the objective function maximizes expected efficiency value of DMUs according to the data from each scenario. Constraints (10.b) are the p-robust restrictions. This set of restrictions may not allow the scenario efficiency take value more than 100 (1-p%) of the ideal efficiency scores gained by each scenario. The parameter p controls the relative regret between all scenarios. If p = ∞, then the p-robust constraints in models (8), (9) and (10) become redundant. Now we assume x_ko = $\max \left\{ {{x}_{io}}|1 \le i\le m \right\}>0$, and then setting $\left({{w}_{1 }}, \ldots, {{w}_{D}}, {{v}_{1 }}, \ldots, {{v}_{m}} \right) = \left(0, \ldots, 0, \ldots, 1 /{{x}_{ko}}, 0, \ldots \right)$, constraints $\sum\nolimits_{i = 1 }^{m}{{}}{{v}_{i}}x_{io}^{s} = 1 \ \ \text{ and}\ \ \sum\nolimits_{d = 1 }^{D}{{}}\text{ }w_{d}^{1 }z_{dj}^{s}-\sum\nolimits_{i = 1 }^{m}{{}}{{v}_{i}}x_{ij}^{s}\le 0, \text{ imply }\ \ \sum\nolimits_{d = 1 }^{D}{{}}\text{ }w_{d}^{1 }z_{do}^{s}\le 1$. Thus, we get $e_{0 }^{1 s*}\le \frac{1 }{1-p}$. So for very small p's, there may not be p-robust solutions for models (8), (9) and (10), and they may be infeasible. We should note that the p-values can be different for any problem and are determined by the decision-maker. In general, the p-values usually should not be considered smaller than 0.2 and its upper bound is obtained by try and error, and can be increased to one.

Theorem 3.2.1 For a specific DMU₀ under different scenarios, f_o^c-S*≥ f_o^1S*×f_o^2S* where f_o^1S* and f_o^2S* are the optimal values of the non-centralized models (8) and (9), respectively and f_o^c-S* is the optimal value of model (10).

Proof Suppose that the optimal value of model (10) is $f_{o}^{c-S *} = \sum_{s = 1 }^{S} q^{s} \sum_{r = 1 }^{t} u_{r}^{*} y_{r o}^{s}$, and the optimal value of models (8) and (9) are $f_{o}^{1 S *} = \sum_{s = 1 }^{S} q^{s} \sum_{d = 1 }^{D} w_{d}^{1 *} z_{d o}^{s}$ and $f_{o}^{2 S*} = \sum\nolimits_{s = 1 }^{S}{{}}{{q}^{S}}\left(\frac{\sum\nolimits_{r = 1 }^{t}{{}}u_{r}^{*}y_{ro}^{s}}{f_{0 }^{1 s*}} \right)$, respectively. Since $\sum_{i = 1 }^{m} v_{i}^{*} x_{i o}^{s} = 1 \text { and } \sum_{d = 1 }^{D} w_{d}^{1 *} z_{d o}^{s} \leq \sum_{i = 1 }^{m} v_{i}^{*} x_{i o}^{s}$, we further get $\sum_{s = 1 }^{S} q^{s} \sum_{i = 1 }^{m} v_{i}^{*} x_{i o}^{s} = 1 \text { and } \sum_{s = 1 }^{S} q^{s} \sum_{d = 1 }^{D} w_{d}^{1 *} z_{d o}^{s} \leq \sum_{s = 1 }^{S} q^{s} \sum_{i = 1 }^{m} v_{i}^{*} x_{i 0 }^{s}$, respectively. Then

Theorem 3.2.2 The expected efficiency values of the centralized model (5) under different scenarios are greater than those provided by model (10).

Proof. Since the feasible region of model (10) is the subset of the feasible region of model (5) and both models have the same objective function, thus the result follows.

Theorem 3.2.3 The expected regret values of model (5) in different scenarios are greater than those provided by model (10).

Proof. Let the difference between the expected efficiency values of each model with the expected ideal efficiency is

It is clear that the smaller value of α for a model shows that model produces more exact results. According to the Theorem 3.2.2, we have

in which ${\bar{\pi }}$ and ${{\phi }^{*}}$ are the optimal objective values to dual of the expected model (5) and dual of model (10), respectively. From inequality (12), we further can get

The left hand side of (13) is the α-value of expected value model (10) and the right hand side shows the α-value of expected value model (5) and thus the result follows.

4.   An application

The proposed models in the previous section are used to evaluate the technical efficiency of an Iranian commercial bank branches. In the realm of performance measurement, banks' performance analysis has been of primary interest due to its socioeconomic importance. Several researchers studied uncertainty subject in the bank (e.g., Kao and Liu, 2009; Paradi and Zhu, 2013). As well, the study of two-stage NDEA models were relatively active in banking industry in recent years and several important works have been reported; see, e.g., Avkiran (2009), Fukuyama and Weber (2010), and Akther et al. (2013) for more details.

It is notable that, on the one hand, due to the existence of some variables, data obtained are not precise and they are estimated with a specific error level (e.g., Zhou et al., 2017). For instance, the earnings quality as a qualitative variable for a bank system is not exactly available; also, the number of customers is not often reported precisely. On the other hand, for the banks, it may be beneficial to conceal real information and reveal deceptive input and output data. Moreover, real and accurate data about key performance criteria of all banks do not always exist. Therefore, it is important to analyze efficiency of banks under uncertainty. Here, we evaluate the presented models under discrete scenarios, provided by the bank network analyzers (i.e., s₁ = Pessimistic, s₂ = Medium, s₃ = Optimistic). We note that the Snyder and Daskin (2006) assumed that all scenarios have the same occurrence probability (i.e. ${{q}^{s}} = {}^{1 }\!\!\diagup\!\!{}_{|S|}\; ;\ \ \forall s\in S$), which should be q^s = 0.33 for each scenario. But, we considered the average probability of the second and third scenarios, because the data-set in these two scenarios were similar. Therefore, the occurrence probabilities of scenarios are considered as 0.25, 0.5 and 0.25, respectively. Table 1 shows the input, intermediate and output parameters employed in the first and second stages of the proposed models.

As well, the results of the normalized data under three scenarios are provided in Table 2.

First, we obtain the ideal efficiency score of each DMU based on each scenario by applying models (2), (3) and (5). Then we utilize models (8), (9) and (10) for the efficiency analysis of DMUs. The related results are given in Table 3. The columns of Table 3 illustrate the ideal efficiency score according to the values defined under each scenario.

As can be seen, DMUs #3, 5, 8, 9, 11, 13, 14, 16, 17 and 19 gained the efficiency score 1 for all scenarios in the leader division that is 50 % of total branches; in the follower division, DMUs #7 and 20 obtained the efficiency score 1 for all scenarios that is 0.1 % of total branches and finally, in the centralized model none of the DMUs are efficient in all scenarios. Next, we solved models (8), (9) and (10) to get the efficiency scores for different p-values and probabilities for each scenario that are reported in Tables 4 to 6. According to these results, models (8), (9) and (10) give infeasible results for some DMUs when p is small, e.g., p ≤0.58 here, thus we do not report the corresponding results here.

As mentioned before in Section 3.2.2, the proposed models give infeasible results in some scenarios when small values are considered for p. The p-value is different for each problem and is determined by the decision-maker. Here, by try and error, we conclude that the p-value must not be smaller than 0.55. In this study, on the one hand, when p ≤ 0.55 according to the results, our models gives infeasible results for most of the DMUs. As the p-value increases, the efficiency score is improved and the number of infeasible DMUs gradually decline and we observe feasible results. On the other hand, for p ≥ 0.74 the results do not change i.e., the efficiency scores remain fixed. So we do not continue and stop it for the other p-values. Therefore, here, we only consider p ≥ 0.55 and do not show the results of p < 0.55.

For instance, with increasing the p-value from 0.56 to 0.70, the efficiency score of DMU16 changes. This change also can be seen in other DMUs such as 18 and 19. Models (8), (9) and (10) maximize the expected efficiency scores of a DMU (here, a bank branch) based on each scenario while p-robust constraints control the relative difference between its efficiency score generated by the model and ideal efficiency in each scenario.

It should be noted that the overall efficiency for stochastic p-robust Stackelberg game models (8) and (9) can be determined as f₀^s* = f₀^1s* ×f₀^2s* . Table 6 shows the product of the efficiency results of first and second stages (i.e., leader and follower models) for three scenarios.

In order to compare models (8), (9) and (10), we consider p = 0.58 and the related results are reported in Tables 6 and 7. After that, we reported the efficiency scores in Table 8 and also illustrated in Figure 2.

As shown in the last two columns of Table 8, the efficiency scores of the stochastic p-robust centralized model (10) yields better results for the overall efficiency scores than the stochastic p-robust Stackelberg model. So it can be deduced that selection of the centralized/cooperative strategy is preferable over the non-centralized/non-cooperative (or, Stackelberg) strategy.

Also, Figure 2 represents efficiency scores of stochastic p-robust centralized model and Stackelberg model that is model (10) and the product of models (8) and (9) for each DMUs with p = 0.58 in three scenarios.

As can be seen, for all DMUs (branches) the stochastic p-robust centralized efficiency scores are better than the other one. Therefore, the stochastic p-robust centralized model as the preferred model is selected and to get deep insight from stochastic p-robust NDEA models, we compare model (10) with deterministic model. To this end, we propound two scales as $\text{EV} = \sum\nolimits_{s\epsilon S}{{}}{{q}^{s}}{{Z}_{s}}(X)$ and $\operatorname{Reg} = \sum\nolimits_{s\epsilon S}{{}}{{q}^{s}}\left(Z_{s}^{*}-{{Z}_{s}}(X) \right)/Z_{s}^{*}$ that the first scale calculates the expected efficiency score of all DMUs by taking into account the occurrence probabilities of the each abovementioned scenario i.e., 0.25, 0.5 and 0.25, respectively and the second scale computes the expected relative regret for each DMU. Ultimately, the calculated efficiency scores by the robust model (10) are compared with these two scales for three scenarios. The results are reported in Table 9 and illustrated in Figure 2. In this experiment, we let p = 0.58.

As shown in this Table 9, the efficiency scores of the p-robust centralized model (10) for all DMUs are smaller than the efficiency score of the expected centralized model (5). On the other hand, we indicated that the expected relative regret in both models is not equal. Moreover, from Table 9 we find that, DMU5 and DMU15 have better performance compared with other units in both of them, whereas DMU17 has better performance compared with other units in the P-robust centralized model. Further, in expected model (5), some DMUs have earned a similar rank, in other words, there is no discrimination among some DMUs. For example, according to Table 9, the efficiency scores in DMUs #8, 13, 18 and 19 are equal, while this is not the case in model (10). Therefore, the performance of the model (10) is like a ranking model, and this affirms the superiority of it. Although, it seems that the efficiency scores in the expected model (5) are greater than model (10) but it does not provide a more distinguishing efficiency score among some DMUs. As can be seen, in model (10), none of the DMUs have equal efficiency scores.

Figure 2 shows the expected efficiency scores between the two models based on the EV scale in each scenario. We should note that, EV scale is only defined for the centralized model (5), and it is the efficiency scores in the model (10).

As seen from Figure 3. the EV scale values show that the efficiency scores of p-robust centralized model (10) is less than the other one. Also, the Reg scale values of model (10) is less than or equal the other one. As mentioned before, the relative regret amount shows the relative difference between the ideal efficiency obtained from each scenario and efficiency of each model that showed by Theorem 3.2.2 and Theorem 3.2.3.

5.   Sensitivity analysis results

In this section, we perform sensitivity analysis for the expected version of the centralized model (5) and model (10) to different probabilities. We ran these models considering a set of different probability vectors. The related results are summarized in Tables 10 and 11, where the first row shows the value of probabilities as ${{q}^{s}} = \left({{q}^{1 }}, {{q}^{2 }}, {{q}^{3 }} \right)$ and columns 2–9 present the difference between the expected efficiency values of each model with the expected ideal efficiency score, in the other words, i.e. the α values. As a matter of fact, this difference, for both abovementioned models, shows the amount of error with pertinent probabilities in each scenario. It is clear that the smaller the difference value, the better the result. Since the expected efficiency value of that model is closer to the ideal expected efficiency value, thus it produces more precise results.

As it can be seen from Tables 10 and 11, our proposed model (10) generates better results compared to the other one. Comparison of the above-mentioned models with varying probabilities corroborates that the stochastic p-robust NDEA model (i.e., model (10)) gives better results than the expected centralized model (5). This supports the advantage of our stochastic p-robust NDEA model (10) in contrast to other NDEA models.

6.   Concluding remarks

Uncertainty is an inherent part of real performance evaluation problems and protecting versus the worst-case is impossible, intricate, and time-consuming. To tackle this issues, and in order to evaluate the efficiency of the two-stage process, a scenario-based robust optimization model according to the Snyder and Daskin's approach (2006) is designed to develop the two-stage NDEA models dealing with uncertainty. This model incorporates the benefits of original stochastic and robust optimization models, therefore, can dominate the difficulties of models in uncertain parameters and give more acceptable results compared the other existing models. The main contributions of this model are as follows. It is less conservative than the other existence worst-case models and enables decision makers to make a benchmark between the expected efficiency of DMUs in different scenarios, a benchmark between robustness of the solution and the robustness of the model. We presented the stochastic p-robust two-stage NDEA model for both cooperative and non-cooperative forms. Further, in Theorems 3.2.1, 3.2.2 and 3.2.3 in Subsection 3.2 it is proved that the choice of the cooperative model provides better efficiency in comparison with the non-cooperative form.

Moreover, the performance of our proposed cooperative model is like a ranking model, and it provides a more distinguishing efficiency score among DMUs. Also, p-robustness adds a feasibility theme that is not present in most other robustness measures. Sensitivity analysis on different probability vectors can derive a spectrum of solutions that may be helpful for managerial tradeoffs. Finally, the applicability of the proposed models, are discussed for a case study in the banking industry showing their advantages over the other models.

The concept of uncontrollable inputs and undesirable outputs has also widespread applications, thus, including them in the model would be an interesting future research direction. Moreover, the parameters of the proposed model may be changed during the planning. Therefore, we can expand the suggested model into the Malmquist model in dynamic condition. Since the return to scale model is linear, under the uncertainty situations, one can apply the proposed model to analyze the efficiency of DMUs. The combination of stochastic p-robust approach with other NDEA models such as additive and slack based NDEA measures with shared resources could also be considered for future research direction.

Acknowledgments

The authors are grateful to the anonymous referees for their useful comments and suggestions.

Conflict of interest

The authors declare no conflicts of interest in this paper.