High-Order Portfolio Optimization Problem with Background Risk

Xiaolu Zhou

Jinan University, Guangzhou, China.

DOI: 10.4236/ojbm.2021.93052                

Abstract

After Markowitz proposed the mean-variance model, the research on portfolio problems has been a hot topic for many investors. The research on portfolio optimization is becoming more and more perfect. The investment theory changes from second-order moment to high-order moment, and from single-stage to multi-stage. More and more factors affecting portfolio optimization are taken into consideration. In this paper, a high-order portfolio optimization problem considering background risks is studied. Firstly, an optimization model of high-order moments including background risks is established, and the genetic algorithm is used to solve the model. Finally, the effects of background risks and high-order moments on the portfolio optimization model are analyzed empirically.

Keywords

Background Risk, Higher Moment, Genetic Algorithm

1. Introduction

After Markowitz proposed the mean-variance model ( Markowitz, 1952), many researchers and investors found some deficiencies in mean-variance model. Some limitations of the mean-variance model are as follows: in order to make the mean-variance model consistent with the expected utility principle, it is necessary and sufficient for the investor’s utility function to be quadratic or to ensure that the return rate of risk assets obeys the normal distribution ( Liu & Liping, 2004; Maringer & Parpas, 2009), but this is not feasible in reality. A large number of studies have shown ( Zhang, He, Zhang, & Pandey, 2017; Zhang, 2017) that the return rate on risk assets is characterized by obvious asymmetry and non-normal distribution with sharp peaks and thick tails. This makes skewness and kurtosis a factor for investors to consider.

From the research of investment theories in recent decades, kurtosis and skewness have become an important factor to consider in portfolio selection ( Naqvi, Mirza, Naqvi, & Rizvi, 2017), and the inclusion of kurtosis is emphasized in portfolio selection and optimization. Including kurtosis and skewness in portfolio selection makes it non-convex ( Kellner, Lienland, & Utz, 2019). Since the inclusion of higher-order moments is non-convexity, various objectives, such as revenue maximization and positive skewness, can be obtained, while the variance and skewness are minimized. FAZ, Gökgöz and Atmaca (2017) think that the addition of skewness not only improves the efficiency of the mean variance portfolio, but also has an impact on the optimization of the portfolio and its selection. Ioannis Oikonomou’s study further confirms the influence of skewness and kurtosis on the portfolio problem ( Oikonomou, Platanakis, & Sutcliffe, 2018). The research results of Khan, Naqvi and Ghafoor (2020) show that the optimal combination including skewness and kurtosis is sustainable and significantly different from the mean-variance optimal combination with asymmetry and fat-tail risk. Therefore, it is not enough to consider only the second-order moment of the mean variance. The consideration of higher-order moments has become an inevitable trend in the study of portfolio problems.

On the other hand, scholars have done a lot of research on portfolio selection under background risk. Background risk usually refers to the risk that cannot be traded in the financial market, such as real estate investment, health condition, etc., and beyond the investor’s control in the short term. Therefore this risk can be considered exogenous. In the general portfolio optimization model, investors usually only consider the endogenous risks such as stock price risk, and less consider the background risks. But there is the risk that the stock price generates, although it is the primary risk, it is not the total risk. Many empirical results show that background risks affect investment decisions. Atatella et al. (2012) show that health risk as a background risk can influence portfolio selection, but this effect can be moderated by a protective and universal national health service. There is also a lot of theoretical research on portfolio selection under background risk. Heaton and Lucas (2000) proposed a theoretical model of portfolio selection decision with background risk, and reached a conclusion consistent with their previous empirical findings, that is, background risk plays an important role in portfolio decision making. Jiang et al. (2010) analyzed the hedging behavior of investors in the face of background risk. Guo et al. (2018) studied the behavior of banking enterprises in the mean-variance model and the risk-taking behavior under the background risk. Then Guo, Chan, Wong, and Zhu (2019) discussed the influence of TIGATING background on the portfolio selection of venture investors, Mean-VaR Mean-CVaR and Mean-Variance framework-work, and analyzed the existence of background risks in the mean-VaR boundary and Mean-VaR effective boundary of characteristics. Xiaoxia Huang and Tingting Yang Donlinks (2020) studied the influence of background risk on individual investment decisions under the framework of uncertainty theory. An uncertain mean-variance model with background risk is proposed, and the optimal solution is given when the return on stock and the return on background assets obey the normal uncertain distribution. All the above studies show that background risk greatly affects portfolio selection. Decisions that do not take into account background risk will lead to investment failure.

2. The Establishment of Model

In order to compare high-order portfolio optimization models with background risks, a general high-order portfolio optimization model with high-order moments is presented.

This model is constructed by considering kurtosis, skewness, variance and expectation. Suppose that in a frictionless market where short selling is not allowed. There are n risky assets with return vector $R=\left(R_1,R_2,\cdots ,R_n\right)$, where the vector $\mu ={\left({\mu }_1,{\mu }_2,\cdots ,{\mu }_n\right)}^{\text{T}}$ represents the expected return vector of risky assets, $x={\left(x_1,x_2,\cdots ,x_n\right)}^{\text{T}}$ is the investment weight vector of risk assets in the portfolio, ${\displaystyle \underset{i=1}{\overset{n+1}{\sum }}{x}_i}=x\left[1\right]=1$. Where is the unit column vector of order n. The return rate of investment portfolio is $R_P=x^{\text{T}}R={\displaystyle \underset{i=1}{\overset{n}{\sum }}{x}_i}{r}_i$, TB asically, In order to solve the aspect problem, the method proposed by Athayde and Flores (2003) is used to represent the model with the tensor in physics. we compute the first four moments of R by.

$\text{Mean}={\mu }_P=E\left(R\right)=x^{\text{T}}\mu ={\displaystyle \underset{i=1}{\overset{n}{\sum }}{x}_{i}E\left(r_i\right)}$.

$\text{Variance}=V_P={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{x}_i{x}_j}}{\sigma }_{ij}=x^{\text{T}}\left(\begin{array}{ccc}{\sigma }_{11}& \cdots & {\sigma }_{1n}\\ ⋮& \ddots & ⋮\\ {\sigma }_{n1}& \cdots & {\sigma }_{nn}\end{array}\right)x=x^{\text{T}}{M}_{2}x$.

$M_2=E\left[\left(R-\mu \right){\left(R-\mu \right)}^{\text{T}}\right]={\left\{{\sigma }_{ij}\right\}}_{n\times n}$ is the $n\times n$ covariance matrix

$\text{Skewness}=S_P={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\displaystyle \underset{m=1}{\overset{n}{\sum }}{\omega }_i{\omega }_j}}}{\omega }_m{s}_{ijm}\left(i\ne j\right)=x^{\text{T}}{M}_3\left(x\otimes x\right)$.

$M_3=E\left[\left(R-\mu \right){\left(R-\mu \right)}^{\text{T}}\otimes {\left(R-\mu \right)}^{\text{T}}\right]={\left\{s_{ijm}\right\}}_{n\times n^2}$ is the $n\times n^2$ Coskewness matrix, $\otimes$ is the Kronecker product.

$\text{Kurtosis}=K_P={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\displaystyle \underset{m=1}{\overset{n}{\sum }}{\displaystyle \underset{l}{\overset{n}{\sum }}{\omega }_i{\omega }_j{\omega }_m{\omega }_l{k}_{ijml}\left(i\ne j\right)}}}}=x^{\text{T}}{M}_4\left(x\otimes x\otimes x\right)$.

$M_4=E\left[\left(R-\mu \right){\left(R-\mu \right)}^{\text{T}}\otimes {\left(R-\mu \right)}^{\text{T}}\otimes {\left(R-\mu \right)}^{\text{T}}\right]={\left\{k_{ijml}\right\}}_{n\times n^3}$ is the $n\times n^3$ Cokurtosis matrix.

The introduction of tensor can make the solution of the problem more convenient. After each order moment of the portfolio is expressed by tensor, the model with higher order moment can be expressed as the following (P1):

$\min K_P=x^{\text{T}}{M}_4\left(x\otimes x\otimes x\right)$

$\max S_P=x^{\text{T}}{M}_3\left(x\otimes x\right)$

$\min V_p=x^{\text{T}}{M}_{2}x$

$\max {\mu }_p=x^{\text{T}}\mu$

$\text{s}\text{.t}\text{.}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{x}_i}=1$

$x_i\ge 0$

Due to the existence of background risks, the background risks should be taken into account in the model construction to achieve a better purpose of risk dispersion. Therefore, the next step is to build a model considering the background risk on the basis of the high-order moment model.

So let’s assume that the total wealth of the investors is zero W₀, meanwhile $\omega \left(0\le \omega \le 1\right)$ Allocating wealth $\omega \left(0\le \omega \le 1\right)$ to risky financial assets. Thus, It means that the wealth available $1-\omega$ is allocated to the background risky assets, $\omega$ is background risk preference ratio. $r_b$ is represent the return on background risk. so, the total return received by the investor at the end of the investment is as follows:

$R_p=W_0\left(\omega X^{\text{T}}r+\left(1-\omega \right)r_b\right)$

And then the expected rate of return at this point is:

${\mu }_p=E\left(R\right)=W_0\left(\omega X^{\text{T}}E\left(r\right)+\left(1-\omega \right)E\left(r_b\right)\right)$

$E\left(r_b\right)$ is represents the expected rate of return of background risk. Minimization kurtosis is still used as the objective function, then the high-order portfolio optimization model with background risks based on this model can be expressed as the following model (P2):

$\min K_P=x^{\text{T}}{M}_4\left(x\otimes x\otimes x\right)$

$\max S_P=x^{\text{T}}{M}_3\left(x\otimes x\right)$

$\min V_p=x^{\text{T}}{M}_{2}x$

$\max {\mu }_p=W_0\left(\omega X^{\text{T}}E\left(r\right)+1-\omega E\left(r_b\right)\right)$

$\text{s}\text{.t}\text{.}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{x}_i}=1$

$x_i\ge 0$

3. Solving Model

Due to the complexity of solving the problem of the established model, this paper adopts genetic algorithm to solve the model ( Yang & Lin, 2014). Genetic algorithm is widely used in intelligent optimization algorithm. Many researchers use intelligent optimization algorithm to solve problems related to investment portfolio. Apply genetic algorithms. The solving steps of model are as follows:

Step 1: We have N chromosomes as the initial population, generates Nthnx1’s 0 - 1 distributed random numbers, So if I take the first random number vector, NX1, and I transform it, it’s going to be a portfolio. That gives you the initial population, then we have N chromosomes.

Step 2: Determination of fitness function, in the model (P3), the multi-objective model was transformed into the single-objective optimization model (P4) ( Yang & Lin, 2014) as the fitness function. as follows:

$\min Z={\left|\frac{d_1}{{\mu }_P}\right|}^{{\lambda }_1}+{\left|\frac{d_2}{V_p}\right|}^{{\lambda }_2}+{\left|\frac{d_3}{S_P}\right|}^{{\lambda }_3}+{\left|\frac{d_4}{K_p}\right|}^{{\lambda }_4}$

$\text{s}\text{.t}\text{.}{\mu }_p=W_0\left(\omega X^{\text{T}}E\left(r\right)+\left(1-\omega \right)E\left(r_b\right)\right)+d_2$

$V_p=x^{\text{T}}{M}_{2}x-d_2$, $S_P=x^{\text{T}}{M}_3\left(x\otimes x\right)+d_3$, $K_P=x^{\text{T}}{M}_4\left(x\otimes x\otimes x\right)-d_4$

${\displaystyle \underset{i=1}{\overset{n}{\sum }}{x}_i}=1$, $x_i\ge 0$, $d_i\ge 0,i=1,2,3,4$

where, is ${\lambda }_i$ the preference coefficient of each moment, $d_i$ is the target variable, Used to evaluate the deviation between the actual moment and the expected value.

Step 3: Selection, crossover, and mutation are performed on chromosomes

Using roulette to select operations: The fitness function max Z obtained according to Step 2 Fitness values for each chromosome $eval\left(f_i\right)=Z\left(F_i\right)$ in count N paternal chromosomes $F=\left(F_1,F_2,\cdots ,F_n\right)$, So for the chromosome $F_i$, the probability of selection is $p_i=F_i/{\displaystyle \underset{i=1}{\overset{N}{\sum }}{F}_i}$, the roulette wheel method was used to select the new population according to the different fitness values of each individual $F^{\left(1\right)}=\left(F_1^{\left(1\right)},F_2^{\left(1\right)},\cdots ,F_N^{\left(1\right)}\right)$.

A crossover operation is performed for crossover probability $P_c$ using an arithmetic crossover algorithm: For the new population that we got using the roulette wheel algorithm $F^{\left(1\right)}=\left(F_1^{\left(1\right)},F_2^{\left(1\right)},\cdots ,F_N^{\left(1\right)}\right)$, $P_c$ as the probability to select chromosome $F_i^{\left(1\right)}$, $F_j^{\left(1\right)}$, A linear combination of them as follow:

$\{\begin{array}{l}{F}_i^{\left(2\right)}=\xi F_i^{\left(1\right)}+\left(1-\xi \right)F_j^{\left(1\right)}\\ F_j^{\left(2\right)}=\xi F_i^{\left(1\right)}+\left(1-\xi \right)F_i^{(\; 1\; )}\end{array}$

So there are produce two new individuals $F_i^{\left(2\right)}$, $F_j^{\left(2\right)}$, apparently vector $F_i^{\left(2\right)}$ add $F_j^{\left(2\right)}$ is one, It’s also an investment portfolio, $\xi$ is arithmetic crossover operator, The new population is obtained by crossover operation $F^{\left(2\right)}=\left(F_1^{\left(2\right)},F_2^{\left(2\right)},\cdots ,F_N^{\left(2\right)}\right)$.

By exchanging the basic bits, the mutation operation is carried out for the mutation probability $P_M$ : If you choose the chromosome $F_i^{\left(2\right)}$ and ${\beta }_i^{\left(2\right)}$ ’s gene of population $F^{\left(2\right)}=\left(F_1^{\left(2\right)},F_2^{\left(2\right)},\cdots ,F_N^{\left(2\right)}\right)$ to process mutate operation, only need to swap ${\beta }_i^{\left(2\right)}$ gene and ${\beta }_{i-1}^{\left(2\right)}$ gene’g in $F_i^{\left(2\right)}$. If you choose the ${\beta }_1^{\left(2\right)}$ gene in $F_i^{\left(2\right)}$, need to swap the ${\beta }_1^{\left(2\right)}$ and ${\beta }_2^{\left(2\right)}$ in $F_i^{\left(2\right)}$. then $F_i^{\left(2\right)}=\left({\beta }_1^{\left(2\right)},{\beta }_2^{\left(2\right)},\cdots ,{\beta }_{i-1}^{\left(2\right)},{\beta }_i^{\left(2\right)},\cdots ,{\beta }_n^{\left(2\right)}\right)$ Mutate thus to $F_i^{\left(3\right)}=\left({\beta }_1^{\left(3\right)},{\beta }_2^{\left(3\right)},\cdots ,{\beta }_{i-1}^{\left(3\right)},{\beta }_i^{\left(3\right)},\cdots ,{\beta }_n^{\left(3\right)}\right)$, by mutating, the new population becomes $F^{\left(3\right)}=\left(F_1^{\left(3\right)},F_2^{\left(3\right)},\cdots ,F_n^{(\; 3\; )}\right)$

Step 4: Repeat Step 3 for M times, and select the chromosome with the highest fitness in each generation, that is, the optimal investment portfolio.

4. Empirical Analysis

This paper randomly selects stocks from 19 different industries in the market by using get_data_yahoo function in Python, and calculates the average annual return of 19 stocks from January 1, 2018 to December 30, 2019. Then, five stocks with higher returns and less correlation coefficients were selected for analysis. Table 1 shows the ticker symbols and expected annualized returns for 19 stocks.

According to the foregoingselects five stocks in China’s stock market for research, which are: Gree (000651), Hengrui Pharmaceutical (600276), Wuliangye (000858), Ping An (601318), Ping An Bank (000001).

Sample data analysis

As can be seen from Table 2, all stock returns are positive, while all skewness is less than 0, indicating that the probability density function of the return rate of portfolio risky assets has the feature of left bias. The kurtosis is less than 3, indicating that the data distribution presents a dumpy feature. The results of the Jarque-Bera statistics show that the data are far from obeying the assumption of a normal distribution. In conclusion, the influence of skewness and kurtosis

Table 1. Code and annualized returns of 19 stocks.

Table 2. Expectation, variance, skewness, kurtosis and Jarque-Bera statistics of stocks.

should be considered when studying portfolio optimization.

Matlab was used to write genetic algorithm, and each parameter in the model was set as: initial wealth $W_0=1$, Population size $N=200$, Crossover probability $P_c=0.5$, The arithmetic crossover operator $\xi =0.7$ ; Mutation probability $P_M=0.05$ Generarions $M=1000$. Three groups of different parameter preferences were selected for the study: 1) ${\lambda }_1=1$, ${\lambda }_2=1$, ${\lambda }_3=0$, ${\lambda }_4=0$ ; 2) ${\lambda }_1=1$, ${\lambda }_2=1$, ${\lambda }_3=3$, ${\lambda }_4=1$ ; 3) ${\lambda }_1=1$, ${\lambda }_2=1$, ${\lambda }_3=1$, ${\lambda }_4=3$.

1) When background risk $\omega =1$, That is, when there is no background risk, The resulting portfolio situation is as follows (Table 3).

The first case (1, 1, 0, 0), At this time, Keede investors will consider investing 24.06% of the amount in Gree, 29.14% in Hengrui Pharmaceutical, 19.92% in Wuliang Liquor, 1.08% in Ping An, and 25.79% in Ping An Bank.

The second case (1, 1, 3, 1). It shows that when higher order moments are considered, that the investors are willing to pursue more excess returns regardless of risk level, At this time, it can be calculated that investors consider investing 26.64% of their funds in Gree, 19.36% in Hengrui Pharmaceutical, 21.61% in Wuliangye, 21.20% in Ping An of China, and 11.19% in Ping An Bank.

The third case (1, 1, 1, 3) Indicates that investors attach more emphasis on risk control, investors consider to invest 35.60% of the funds in Gree, 9.30% in Hengrui Pharmaceutical, 38.34% in Wuliang Liquor, 1.73% in China Ping’an, and 15.04% in Ping An Bank.

The above results show that different preference portfolios have a great impact on investors’ choice of different assets.

2) When there is a background risk, hypothesis $\omega =1/2$, The resulting portfolio situation is as follows (Table 4):

Contrast that with the situation without background risk

Investment in Gree 36.97%, Hengrui Pharmaceutical 1.34%, Wuliangye 17.07%, Ping An 14.84% and Ping An Bank 29.78% are required while preference parameters is (1, 1, 0, 0).

Table 3. Portfolio of preference moments without background risk.

Table 4. Portfolios with different preference degrees when the background risk preference $\omega =1/2$.

We need to invest 31.02% in Gree, 7.29% in Hengrui Pharmaceutical, 19.10% in Wuliangye, 17.38% in Ping An of China and 25.21% in Ping An Bank while preference parameters are (1, 1, 3, 1).

We need to invest 23.07% in Gree, 24.54% in Hengrui Pharmaceutical, 8.98% in Wuliangye, 22.11% in Ping An of China and 21.30% in Ping An Bank while preference parameters are (1, 1, 1, 3).

By comparison, it can be found that when the background risk is included, the investment proportion of each asset has changed greatly. Therefore, considering the background has a great impact on the portfolio optimization and should not be ignored.

5. Conclusion

In this paper, background risk is added into the high-order moment model considering portfolio optimization to form a new model containing background risk. Due to the difficulty of analyzing and solving the model, this paper adopts the genetic algorithm in the intelligent optimization algorithm to solve the problem, and transforms the multi-objective optimization problem into the single-objective optimization problem.

In the end, this paper selects five stocks in Chinese stock market: Gree (000651), Hengrui Pharmaceutical (600276), Wuliangye (000858), Ping An (601318) and Ping An Bank (000001) for research. The comparative analysis shows that the investment proportion of the five stocks changed differently after the inclusion of background risk compared with that without background risk. Therefore, the consideration of background risk should be considered in the portfolio optimization problem.

Because this paper establishes a static portfolio model, and in the real economic market, portfolio returns are time-varying, and different time points will present different distributions. Therefore, in the next step, the static model can be transformed into a dynamic model for research, and the portfolio can be studied from a dynamic perspective to make it more suitable for the market.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Atatella, V., Brunetti, M., & Maestas, N. (2012). Household Portfolio Choices, Health Status and Health Care Systems: A Cross-Country Analysis Based on Share. Journal of Banking and Finance, 36, 1320-1335.
[2]	Athayde, G. M., & Flores Jr., R. G. (2003). Incorporating Skewness and Kurtosis in Portfolio Optimization: A Multidimensional Efficient Set (pp. 243-257). Oxford: Butterworth-Heinemann.
[3]	Gokgoz, F., & Atmaca, M. E. (2017). Portfolio Optimization under Lower Partial Moments in Emerging Electricity Markets: Evidence from Turkey. Renewable and Sustainable Energy Reviews, 67, 437-449.
[4]	Guo, X., Chan, R. H., Wong, W. K., & Zhu, L. (2019). Mean-Variance, Mean-VaR, Mean-CVaR Models for Portfolio Selection with Background Risk. Risk Management, 21, 73-98.
[5]	Guo, X., Wagener, A., Wong, W. K., & Zhu, L. (2018). The Two-Moment Decision Model with Additive Risks. Risk Management, 20, 77-94.
[6]	Heaton, J., & Lucas, D. (2000). Portfolio Choice in the Presence of Background Risk. The Economic Journal, 110, 1-26.
[7]	Huang, X. X., & Yang, T. T. (2020). How Does Background Risk Affect Portfolio Choice: An Analysis Based on Uncertain Mean-Variance Model with Background Risk. Journal of Banking and Finance, 111, Article ID: 105726. https://doi.org/10.1016/j.jbankfin.2019.105726
[8]	Jiang, C., Ma, Y., & An, Y. (2010). An Analysis of Portfolio Selection with Background Risk. Journal of Computational and Applied Mathematics, 330, 59-69.
[9]	Kellner, F., Lienland, B., & Utz, S. (2019). An a Posteriori Decision Support Methodology for Solving the Multi-Criteria Supplier Selection Problem. European Journal of Operational Research, 272, 505-522.
[10]	Khan, K. I., Naqvi, S. M. W. A., Ghafoor, M. M. et al. (2020). Sustainable Portfolio Optimization with Higher-Order Moments of Risk. Sustainability, 12, 2006.
[11]	Liu, L. P. (2004). A New Foundation for the Mean-Variance Analysis. European Journal of Operational Research, 158, 229-242.
[12]	Maringer, D., & Parpas, P. (2009). Global Optimization of Higher Order Moments in Portfolio Selection. Journal of Global Optimization, 43, 219-230. https://doi.org/10.1007/s10898-007-9224-3
[13]	Markowitz, H. M. (1952). Portfolio Selection. Journal of Finance, 7, 77-91. https://doi.org/10.1111/j.1540-6261.1952.tb01525.x
[14]	Naqvi, B., Mirza, N., Naqvi, W. A., & Rizvi, S. (2017). Portfolio Optimisation with Higher Moments of Risk at the Pakistan Stock Exchange. Economic Research-Ekonomska Istrazivanja, 30, 1594-1610. https://doi.org/10.1080/1331677X.2017.1340182
[15]	Oikonomou, I., Platanakis, E., & Sutcliffe, C. (2018). Socially Responsible Investment Portfolios: Does the Optimization Process Matter? The British Accounting Review, 50, 379-401. https://doi.org/10.1016/j.bar.2017.10.003
[16]	Yang, J. H., & Lin, R. J. (2014). The Higher Moment Portfolio Model Based on Genetic Algorithm. Henan Science, 32, 697-702.
[17]	Zhang, T. X. (2017). An Improved High-Moment Method for Reliability Analysis. Structural and Multidisciplinary Optimization, 56, 1225-1232. https://doi.org/10.1007/s00158-017-1715-3
[18]	Zhang, X. F., He, W., Zhang, Y. M., & Pandey, M. D. (2017). An Effective Approach for Probabilistic Lifetime Modelling Based on the Principle of Maximum Entropy with Fractional Moments. Applied Mathematical Modelling, 51, 626-642. https://doi.org/10.1016/j.apm.2017.07.036