A Monte Carlo-Based Framework for Two-Stage Stochastic Programming: Application to Bond Portfolio Optimization

TL;DR

This paper introduces a Monte Carlo method to optimize bond portfolios under uncertain market conditions, ensuring cost efficiency and liability fulfillment.

Contribution

A novel Monte Carlo-based framework for two-stage stochastic programming applied to bond portfolio optimization.

Findings

The approach successfully minimizes bond portfolio costs under random market conditions.

The method meets liabilities effectively, providing robust portfolio solutions.

The algorithm determines the required number of scenarios to convert stochastic problems into deterministic ones.

Abstract

This paper presents a Monte Carlo simulation-based approach for solving stochastic two-stage bond portfolio optimization problems. The main objective is to optimize the cost of the bond portfolio while making decisions on bond purchases, holdings, and sales under random market conditions such as interest rate fluctuations and liabilities. The proposed algorithm identifies the number of randomly generated scenarios required to convert the stochastic problem into a deterministic one, subsequently solving it as a Mixed-Integer Linear Program. The practical relevance of this research is shown through an application of the proposed method to a real-world bond market. The results indicate that the proposed approach successfully minimizes costs and meets liabilities, providing a robust solution for bond portfolio optimization.

Figures

Full text

1. Introduction

This paper focuses on the use of two-stage stochastic programming in bond portfolio optimization. In the realm of financial decision-making, uncertainty plays a critical role, particularly in bond portfolio optimization where market parameters such as interest rates and liabilities are inherently volatile. Two-stage stochastic programming provides a powerful modeling framework for bond portfolio optimization problems by enabling early-stage decisions followed by adaptive responses to realized outcomes. Bond portfolios, even those constructed conservatively, are subject to substantial risks due to fluctuating interest rates, bond prices, cash requirements, liabilities, and potential defaults. Effective management of these portfolios requires robust optimization techniques capable of handling such uncertainties.

The portfolio selection problem is one of the most prominent research topics in modern finance (Ardakani 2024 [1]; Yan et al., 2023 [2]; Novais et al., 2022 [3]). Finding a combination of assets that satisfies an investor’s needs is the main goal of portfolio selection problem theory. The theory of portfolio optimization was developed by Markowitz (1952) [4]. Bond portfolios are subject to major uncertainties and risks even when they are conservatively assembled. Some of the uncertainties associated with the bond portfolio model can be seen as uncertainties in the interest rates and their effects on bond prices, including the cash requirements or liabilities of the portfolio and the possibility of bond re-calls or defaults due to bankruptcy, to name just a few (Shapiro 1988 [5]). The dynamic model of the bond portfolio problem proposed by Bradley and Crane (1972) [6] can be viewed as a multi-stage decision problem in which actions such as “buying, holding, or selling” bonds are taken at successive time points.

In this paper, we propose a Monte Carlo simulation-based methodology to solve a two-stage stochastic bond portfolio optimization problem, incorporating concepts from uncertainty quantification to guide scenario generation and improve decision quality. Specifically, we use the sample average approximation (SAA) approach to reduce the uncertainty space by generating a representative set of market scenarios (Hu et al. 2020 [7]; Jiang and Li 2021 [8]). SAA leverages iterative random sampling from input distributions to estimate objective function values, thereby directing the optimization process (Bertsimas et al., 2018 [9]). Through the systematic generation of random scenarios, the methodology transforms the stochastic framework into a deterministic problem structure, enabling resolution via Mixed-Integer Linear Programming (MILP) techniques implemented through commercial optimization software. An entropy-based approach can also be used for portfolio optimization (Novais et al., 2022 [3]; Mercurio et al., 2020 [10]).

Let x and y designate the first- and second-stage decision vectors, respectively, and let $[eqn]$ be the random vector to be observed. Then, the decision-observation scheme (Prékopa 1995 [11]) is

$$[eqn]$$

The second-stage problem is formulated under the assumption that the values of x and $[eqn]$ are known:

$$[eqn]$$

Assume that the first-stage decision vector x satisfies the deterministic constraints: $[eqn]$ , and let F be the set of all those x vectors for which problem (1) has a feasible solution for every possible value of the random vector $[eqn]$ . Matrices $[eqn]$ associated with vector x and $[eqn]$ associated with vector y are called the technology matrix and the recourse matrix, respectively.

Let $[eqn]$ designate the optimum value of problem (1), and  $[eqn]$ be the recourse function. Then, the first-stage problem can be formulated as:

$$[eqn]$$

equivalently,

$$[eqn]$$

 where

$[eqn]$ : first-stage decision vector, $[eqn]$ : second-stage decision vector, $[eqn]$ : a random vector to be observed,F: set of all $[eqn]$ for which problem (1) has feasible solutions for every possible value of the random vector $[eqn]$ , $[eqn]$ : cost vector corresponding to the first-stage decision vector x, $[eqn]$ : cost vector corresponding to the second-stage decision vector y, $[eqn]$ : coefficient matrix, and  $[eqn]$ : right-hand side vector of the equality constraints.

Let $[eqn]$ and $[eqn]$ denote the support of the random vector $[eqn]$ . Then, problem (3) has a finite optimum if, and only if, there exists a vector $[eqn]$ such that $[eqn]$ . This assertion follows from the duality theory in linear programming (Prékopa 1995 [11]). If such vector z exists (z is dual feasible) and the expected value of the random variable $[eqn]$ exists, then for every $[eqn]$ , the expectation $[eqn]$ exists and is a convex function on the convex polyhedron F (Prékopa 1995 [11]).

One approach to solve problem (3) is to approximate $[eqn]$ using the Monte Carlo simulation-based method (Shapiro 2001 [12]). Monte Carlo sampling constitutes a widely implemented methodology for addressing stochastic optimization challenges through iterative random sampling procedures. The approach estimates objective function values via repeated sampling from specified input distributions, utilizing these estimates to direct the optimization trajectory. This computational framework operates by generating substantial quantities of random samples from appropriate probability distributions—whether uniform, normal, exponential, or others—contingent upon problem characteristics. Once the samples have been generated, an appropriate optimization algorithm can be used to estimate the objective function and search for an optimal solution (Brandimarte 2014 [13]).

The motivation of this research stems from challenges in bond portfolio optimization under uncertainty, particularly in the Saudi Sukuk (bond) market. To address these challenges, we adopt and extend the two-stage stochastic programming framework by introducing a novel adaptive sampling rule, thereby combining methodological innovation with practical relevance. While two-stage stochastic programming combined with Monte Carlo simulation and sample average approximation (SAA) is well established in the literature, the originality of our study lies in extending these methods both methodologically and in application. Methodologically, we develop a novel adaptive sampling procedure that introduces a dual-condition stopping criterion based simultaneously on hypothesis tests of the sample mean and the sample variance across successive iterations. If both tests confirm equality, the algorithm terminates since enlarging the sample would not yield further improvement; otherwise, the sample size is expanded and the procedure is repeated. This ensures that the sample size is increased only when statistically significant improvements are observed, thereby enhancing computational efficiency without compromising solution quality. To the best of our knowledge, such a dual statistical stopping mechanism has not been applied before in stochastic programming. From an application perspective, although the SAA approach is common in stochastic programming in general, its direct use in bond portfolio optimization remains virtually unexplored. Existing works on bond portfolios under uncertainty (Alreshidi et al., 2020 [14]; Alkailany et al., 2022 [15]) typically rely on fixed scenario trees or deterministic recourse models. By adopting the SAA in this context, and further innovating with the proposed dual statistical stopping criterion, our paper extends a general methodology into a new and underdeveloped application area, thereby contributing both methodological novelty and practical relevance.

The remainder of this paper is organized as follows: Section 2 reviews the relevant literature on bond portfolio optimization and stochastic programming. Section 3 presents the basic mathematical models for the studied problem. Section 4 details the methodology of the Monte Carlo simulation approach. Section 5 presents a case study applying the proposed method to the Saudi Sukuk market and discusses its results and implications. Finally, Section 6 concludes the paper with a summary of findings and potential directions for future research.

2. Literature Review

The portfolio selection problem has been a crucial research topic in modern finance, with the primary goal of identifying asset combinations that fulfill investor needs. The foundational theory of portfolio optimization was introduced by Markowitz (1952) [4], establishing a basis for understanding and managing the uncertainties and risks inherent in bond portfolios. These risks include fluctuating interest rates, bond prices, cash requirements, liabilities, and the potential for bond recalls or defaults (Shapiro 1988 [5]).

Bradley and Crane (1972) [6] pioneered a dynamic bond portfolio model, conceptualizing it as a multi-stage decision process involving sequential actions such as buying, holding, or selling bonds. This model was enhanced by Hodges and Schaefer (1977) [16], demonstrating that a deterministic optimization approach could reduce bond portfolio cost and initial investments while maintaining constraints and reinvesting remaining cash. Ronn (1987) [17] further developed this model by integrating bid and ask prices to reflect the investor’s yield curve.

In addressing the complexities of cash flows, interest rates, and obligations, Bradley and Crane (1972) [6] introduced a dynamic multi-period bond portfolio model utilizing a decomposition technique. Sequential decision theory and two-stage programming under uncertainty were employed by Shapiro (1988) [5]. Shapiro’s model incorporated multiple scenarios to address uncertainties in fixed-income portfolio selection.

Korn and Koziol (2006) [18] expanded on the bond portfolio optimization problem, further exploring stochastic programming approaches. Stoyan and Kwon (2011) [19] applied a stochastic-goal mixed-integer programming method to integrate stock and bond portfolios. He and Qu (2014) [20] proposed a two-stage mixed-integer programming model that addressed trading constraints and market uncertainty.

Recent contributions to the field include Maggioni et al. (2020) [21], that provided three mathematical formulations for multi-horizon stochastic programs, and Alreshidi et al. (2020) [14], that introduced new mixed-integer stochastic programming models for bond portfolio optimization in the Saudi bond market. Alkailany et al. (2022) [15] applied the two-stage stochastic programming model developed by Alreshidid et al. (2020) [14] to the U.S. bond market with discrete random variables. A recent comprehensive survey by Salo et al. (2024) [22] highlights both the enduring relevance and the contemporary advances of portfolio optimization, providing context for our contribution to the bond portfolio setting.

Stochastic programming, also known as optimization under uncertainty, is a crucial branch of optimization that addresses problems involving uncertain parameters. In such scenarios, some of the input parameters in the objective function or constraints are random variables, which significantly complicate the decision-making process (see, e.g., Prékopa, 1995 [11]). To manage this uncertainty, two primary approaches are typically employed: analyzing the statistical characteristics of the random parameters or reformulating the problem to account for the joint probability distribution of these parameters. The resulting numerical problem is often complex and requires sophisticated mathematical programming techniques to solve (Prékopa 1995 [11]).

Stochastic programming represents a prominent approach among uncertainty modeling techniques in optimization, offering broad adoptibility and implementation potential across scientific and engineering domains. The integration of probability theory, optimization, statistics, and functional analysis forms the mathematical foundation of stochastic programming methodology. Contemporary computational capabilities and algorithmic developments have facilitated the field’s evolution, enabling increasingly sophisticated modeling frameworks and solution methodologies, as demonstrated through extensive research (Zhang et al., 2023 [23]; Herion and Römisch 2022 [24]; Chan et al. 2025 [25]).

The evolution of stochastic programming in finance traces back to Markowitz’s (1952) [4] seminal work on portfolio selection, which laid the groundwork for optimization under uncertainty. Dantzig and Infanger (1993) [26] extended these concepts to multi-stage stochastic programs for portfolio management. A stochastic network programming approach for financial planning problems was presented by Mulvey and Vladimirou (1992) [27]. Some other applications include fixed-income securities management and asset-liability modeling by Zenios (1995) [28], management of real-world investment portfolios by Carino et al. (1998) [29], optimal portfolio selection by Dentcheva and Ruszczynski (2006) [30], portfolio management by Gülpınar et al. (2016) [31], and methods for incorporating risk-averse approach by Shapiro and Ugurlu (2016) [32]. Following the 2008 financial crisis, Jobst et al. (2009) [33] developed integrated liquidity and market risk models using multi-stage stochastic programming, and Fábán et al. (2011) [34] introduced risk-averse portfolio optimization approaches. Recently, Bertsimas and Kallus (2020) [35] incorporated machine learning techniques within stochastic programming frameworks to enhance predictive accuracy in portfolio optimization.

Hodges and Schaefer (1977) [16] presented a deterministic optimization model to improve the bond portfolio selection performance efficiently. The model minimizes the cost and initial investment of a bond portfolio while constraints are satisfied at each period with the same cash flow pattern, assuming any remaining cash is used for reinvestment in the current period. Later, Ronn (1987) [17] developed a model of bond portfolio based on Hodges and Schaefer’s (1977) [16] findings, introducing the bid and ask price of bond into the model. The proposed model may display an investors yield curve. Bradley et al. (1972) [6] introduced a dynamic model for multi-period bond portfolio optimization and proposed a decomposition technique to address the model’s computational complexity. They explored two major approaches, sequential decision theory and two-stage stochastic programming, to solve multi-asset, two-stage portfolio problems, explicitly incorporating uncertainty in cash flows, interest rates, and financial obligations. As they noted, normative models for such decision problems tend to become significantly large, particularly when the model’s dynamic structure and the stochastic nature of future interest rates and cash flows are fully incorporated. The second approach allows bond portfolio managers to minimize the cost of a bond portfolio taking into account various constraints on the bonds using the historical data (Shapiro 1988 [5]).

A stochastic programming model for fixed income portfolio optimization was developed by Shapiro (1988) [5] to treat model uncertainties. The idea is to cover most of the future expected cases by considering a different number of scenarios. The optimal solution strategy for selling, buying, and/or holding for the bond portfolio is calculated by the proposed model, where cash flows, interest rates, and/or liabilities may have uncertainties (Shapiro 1988 [5]). For more information on the bond portfolio optimization problem, the reader is referred to Korn and Koziol (2006) [18], Alreshidi et al. (2020) [14] and the references therein.

A Stochastic-Goal Mixed-Integer Programming approach for an integrated stock and bond portfolio problem was used by Stoyan and Kwon (2011) [19]. He and Qu (2014) [20] proposed a two-stage mixed-integer programming model for a portfolio selection issue that takes into account a wide range of actual trading constraints in addition to market random uncertainty over asset prices. Complementary to these approaches, Consigli et al. (2025) [36] developed asset-liability management models with sequential stochastic dominance, illustrating the growing focus on risk-aware extensions of stochastic programming. The reader is referred to a recent article by Maggioni et al. (2019) [21] in which three general mathematical formulations of a multi-horizon stochastic program are presented, and definitions of the classical Expected Value problem and Wait-and-See problem are expanded in a multi-horizon framework. Alreshidi et al. (2020) [14] proposed three new mixed-integer stochastic programming models for two-stage bond portfolio optimization and its application on the Saudi Sukuk (bond) market, where an investor may optimize the cost of a bond portfolio under different scenarios. Alkailany et al. (2022) [15] presented a two-stage bond portfolio optimization model and its application in the U.S. market, assuming the random input parameters are independent, identically distributed discrete random variables. Related research has also applied multi-stage stochastic programming to corporate bond issuance, underscoring the growing importance of stochastic methods in fixed-income markets (Yang et al., 2022 [37]). Recent studies emphasize the role of adaptive two-stage stochastic models in large-scale planning contexts, such as capacity expansion (Basciftci et al., 2024 [38]).

In this research, we present a stochastic two-stage bond portfolio model as a modification of a model introduced by Alreshidi et al. (2020) [14]. Unlike the previous study, where the underlying distribution is discrete, we assume that the random variables associated with the interest rate and second-stage obligations follow a continuous distribution. As a result, the problem involves an infinite number of scenarios, making it more complex and realistic in modeling financial uncertainties.

While the model proposed by Alreshidi et al. (2020) [14] serves as a foundational reference, we propose a Monte Carlo Simulation framework, which can effectively handle the computational complexity and is particularly well suited for problems characterized by continuous uncertainty.

The proposed algorithmic framework offers several advantages over prior studies. It enhances the accuracy of scenario generation through continuous probability distributions, leading to more reliable decision-making outcomes. Additionally, the use of Monte Carlo simulation improves computational efficiency in handling high-dimensional uncertainty, making the model both theoretically significant and practically applicable in real-world financial decision-making, particularly in bond portfolio optimization.

The exact origin of using Monte Carlo simulation methods to solve stochastic programming problems is difficult to pinpoint. Rubinstein and Shapiro (1993) [39] covered the stochastic counterparts algorithm, laying foundational concepts for optimization using Monte Carlo methods. Similar ideas were employed by Geyer and Thompson (1992) [40] to compute maximum likelihood estimators through Gibbs sampling-based Monte Carlo approaches. Infanger (1992) [41] developed algorithms combining classical decomposition and Monte Carlo sampling to effectively solve large-scale stochastic linear programs.

It is worth noting that the term “sample average approximation (SAA)” is not uniformly used in the literature. Plambeck et al. (1996) [42] referred to it as “sample-path optimization.” This approach, also known as sample average approximation, is based on Monte Carlo concepts (Gurkan et al., 1994 [43]). Shapiro and Wardi (1996) [44] provided a framework ensuring the convergence of a Monte Carlo-based approximation method with a probability one. Later on, Shapiro and Homem-de-Mello (1998) [45] illustrated the methodology of using a Monte Carlo simulation-based approach on a two-stage stochastic program with recourse, discussing algorithmic details, fundamental principles, and providing several numerical examples. Moreover, the error estimation, stopping rules, and validation analysis were developed by incorporating statistical inference into numerical algorithms (Shapiro 1988 [5]).

Sampling-based approaches have been effectively applied in various stochastic optimization contexts. Notable applications include vehicle routing (Verweij et al., 2023 [46]), engineering design (Royset and Polak 2004 [47]), supply chain network design (Santoso et al., 2005 [48]), generation and transmission (Jirutitijaroen and Singh 2008 [49]), scheduling (Keller and Bayraksan 2009 [50]), and asset management (Hilli et al., 2007 [51]). The most recent survey on Monte Carlo sampling-based approaches for stochastic optimization is discussed by Homem-de-Mello and Bayraksan (2014) [52]. Their work provides a comprehensive overview, offering practical guidance for tackling stochastic optimization problems using sampling methods. Recent advances such as adaptive sequential SAA methods (Pasupathy and Song 2021 [53]) demonstrate how iterative control of sample sizes can improve efficiency. Extensions of the SAA approach have increasingly focused on integrating auxiliary information, such as covariates, to improve sample efficiency (Kannan et al. 2025 [54]), highlighting the ongoing evolution of adaptive sampling methods. Further, a recent work by Lew et al. (2022) [55] examined limitations of the SAA under equality constraints and non-convex problems, reinforcing the need for adaptive strategies.

In the following section, we provide the two mathematical models that will serve as basis for our Monte Carlo simulation approach.

3. Basic Mathematical Models for Bond Portfolio Optimization

3.1. Bond Portfolio Optimization Problem for Stochastic Programming with Recourse

Our modeling approach builds on the classical two-stage recourse stochastic programming framework (Shapiro 1988 [5]). In this framework, a set of first-stage (here-and-now) decisions is taken before the uncertainty is revealed. After the random outcomes are realized, second-stage (recourse) decisions are made to adjust the initial plan at an additional cost. Formally, the generic two-stage stochastic programming model can be written as

$$[eqn]$$

where x represents the first-stage decision vector, c is the associated cost, $[eqn]$ is a random vector of uncertain parameters, and  $[eqn]$ denotes the optimal value of the recourse problem given x and realization $[eqn]$ . This structure provides a flexible foundation for modeling portfolio optimization problems under uncertainty.

In the context of bond portfolio optimization, Shapiro’s recourse stochastic programming model provides a foundation for addressing uncertainties. The model employs the following notations (Shapiro 1988 [5]):

• B: number of bonds
• T: number of time periods
• $[eqn]$ : amount of bond j purchased for all $[eqn]$
• $[eqn]$ : cost of bond j for all $[eqn]$
• $[eqn]$ : last time period when all parameters are known with certainty
• $[eqn]$ : cash surplus accumulated at the end of the period t for all $[eqn]$
• $[eqn]$ : cash flows obtained from bond j at time t for $[eqn]$ and $[eqn]$
• K: number of scenarios of the uncertain future
• $[eqn]$ : re-investment rate under scenario k for period t, where $[eqn]$ is the interest rate for all $[eqn]$ and $[eqn]$
• $[eqn]$ : cash liability to be met in period t under scenario k for all $[eqn]$ and $[eqn]$
• $[eqn]$ : probability that scenario k will occur for all $[eqn]$
• $[eqn]$ : additional cash required in period t under scenario k for all $[eqn]$ and $[eqn]$
• $[eqn]$ : discount factor for cash flows in period t under scenario k for all $[eqn]$ and $[eqn]$

Then, problem (5) minimizes the cost of bonds purchased in the first stage and the expected cost of bonds purchased in the second stage under various scenarios, while maximizing the expected gain from selling bonds that have not yet matured. The constraints of problem (5) ensure that obligations are met by cash flows from bonds and reinvested surplus (Shapiro 1988 [5]).

$$[eqn]$$

3.2. Two-Stage Mixed-Integer Stochastic Programming Models for Bond Portfolio Optimization

In this section, we consider an investor who aims to secure the required liabilities over a fixed time period through investments in bonds. The investor seeks to develop an optimal investment plan that ensures these liabilities are met through the returns generated from the bond portfolio. To achieve this, the investor needs to make strategic decisions regarding the selection and allocation of bonds, considering both immediate investment opportunities (first-stage) and future adjustments based on market uncertainties (second-stage). This problem naturally leads to the formulation of a two-stage bond portfolio optimization model, where the initial investment decisions are made before the uncertainty is realized, and recourse actions are taken in the second stage to adjust the portfolio in response to observed financial conditions.

Alreshidi et al. (2020) [14] formulated a new bond portfolio optimization model as a two-stage stochastic programming problem. The model’s objective is to maximize the expected gain from selling bonds under various scenarios, defined by interest rates in each period, while minimizing the cost of purchasing bonds in both stages. Alreshidi et al. (2020) [14] formulated the following two-stage bond portfolio optimization model:

$$[eqn]$$

where the input parameters and decision variables are as defined in Table 1 and Table 2, respectively.

The objective function of problem (6) aims to minimize the cost of bonds purchased in the first-stage, $[eqn]$ , where $[eqn]$ is the buying price of the bond j in Stage 1.

Similarly, it seeks to reduce the expected cost of bonds purchased in the second-stage under scenario k, $[eqn]$ , where $[eqn]$ is the buying price of bond j purchased under scenario k as defined before, and  $[eqn]$ is the probability that scenario k will occur. Additionally, it aims to maximize the expected gain from selling bonds not yet matured in Stage 2, $[eqn]$ , where $[eqn]$ is the selling price of bond j under scenario k at time t, and  $[eqn]$ is the probability that scenario k will occur.

Overall, the objective is to minimize the initial investment, $[eqn]$ , by balancing the costs of bond purchases in both stages and maximizing the expected returns from bond sales under various scenarios.

Note that the first set of equality constraints in problem (6) ensures that for each period $[eqn]$ , the obligation requirements are met by the cash flow from bonds j (which have not yet matured) purchased in Stage 1, combined with the surplus from the previous period reinvested in the current period. These constraints are referred to as the liability constraints for Stage 1.

The second set of equality constraints represents the liability constraint at the transition stage $[eqn]$ , moving from Stage 1 to Stage 2, where $[eqn]$ is the last time period in which all parameters are known with certainty under scenario $[eqn]$ . For this period, the liability requirement will be met by holding or selling the bonds owned or buying new bonds only in the first period of the second stage. Note that the first term is the total cash flow obtained from bonds purchased in Stage 1 that have not yet reached maturity. The second term represents the cash flow for bonds sold in the first period in Stage 2. The third and fourth terms are the interest and cash from the sales of the bonds in Stage 2 and cash flow from the bonds purchased in Stage 2. Note also that the cash flow for the transition period uses the cash surplus $[eqn]$ , which is accumulated at the end of the final period of Stage 1 and may be used for the purchase of new bonds in the first period of Stage 2 under any scenario k.

The third set of equality constraints is similar to the transition stage constraints. It represents the remaining periods, $[eqn]$ . The main difference is the calculation of cash surplus that will be used to purchase new bonds in the current period of Stage 2, in this case, will be the cash accumulated at the end of the previous period of Stage 2.

The fourth set of constraints ensures that the amount of bond j purchased in Stage 1 is greater than the total amount of bond j for any $[eqn]$ (that has not yet reached the maturity) to be sold in Stage 2 under any scenario k at the beginning of the period $[eqn]$ . The integrality constraints in problem (6) have been added for the amount of bonds that may be purchased or sold within the planning horizon in addition to the liability constraints. We set an upper bound for the amount of bond j that can be purchased in Stage 1 or Stage 2.

A thorough examination of problem (6) proposed by Alreshidi et al. (2020) [14] leads to the observation that the amount of funds resulting from selling bonds in Stage 2 is subtracted from the objective function (cost of bond portfolio); on the other hand, those funds are used to meet liabilities constraints. Therefore, to respect accounting laws, we propose to modify problem (6) by replacing the second and third set of liability constraints with the following constraints, where reinvestment rates, $[eqn]$ , are replaced with interest rates, $[eqn]$ , under scenario k and for periods $[eqn]$ in Stage 2:

$$[eqn]$$ $$[eqn]$$

This results in the following two-stage bond portfolio optimization model:

$$[eqn]$$

where input parameters and decision variables are as defined in Table 1 and Table 2, respectively.

4. A Monte Carlo Simulation Approach to Solve Problem (9)

The Monte Carlo simulation method is a powerful tool for solving stochastic programming problems. This method involves using random sampling to approximate mathematical functions and simulate the behavior of complex systems. Its flexibility and robustness make it particularly well suited for optimization problems that involve uncertain or random elements. The basic idea of Monte Carlo simulation in the context of stochastic optimization is to generate a large number of random samples from the input distribution, evaluate the objective function for each sample, and use these estimates to guide the search for the optimal solution. This approach is especially useful when the expectation function in a stochastic programming model cannot be calculated accurately or efficiently.

We propose a Monte Carlo simulation-based approach to solve problem (9). Our method is motivated by Algorithm 1 developed by Homem-de-Mello and Bayraksan (2014) [52] that employs sampling methods in stochastic optimization as outlined below: Algorithm 1 Sample Average Approximation Framework (Homem-de-Mello and Bayraksan 2014 [52])

1. Initialization: Choose an initial guess $[eqn]$ , and set iteration counter $[eqn]$ .
2. Sampling: Obtain a realization $[eqn]$ of $[eqn]$ .
3. Optimization: Perform optimization steps on the function

$$[eqn]$$

  possibly using information from previous iterations to obtain $[eqn]$ . 4. Convergence: Check the stopping criteria; if not satisfied, then set $[eqn]$ and return to Step 2.

Algorithm 1 can be viewed as a general framework that allows for numerous variations of the sample average approximation (SAA) approach (Homem-de-Mello and Bayraksan 2014 [52]). To develop a variety of SAA algorithms, several key questions must be addressed, including determining the sample size for each iteration, generation of realizations from previous iterations, selection of a method for the solution of the underlying optimization problem, and determining the stopping criteria. To approximate the expectation in the stochastic objective function, we apply the SAA technique. From an information-theoretic perspective, each scenario sampled represents a discrete state of nature that contributes to the overall uncertainty in the model. By generating a sufficiently large and diverse set of scenarios, we aim to approximate the true probability distribution of market behaviors.

To solve problem (9), we propose to modify Algorithm 1 as described in Algorithm 2: Algorithm 2 Modified Sample Average Approximation Framework to Solve Two Stage Problems

1. Initialization: Set $[eqn]$ , and  $[eqn]$ , where $[eqn]$ and $[eqn]$ are initial values of sample size and number of realizations, respectively.
2. Sampling: Generate R different realizations, each with k different samples

$$[eqn]$$

3. Optimization: Compute $[eqn]$ , …, $[eqn]$ and solve the underlying integer linear program using deterministic optimization techniques.
4. Estimation: Compute the average optimal cost and its variance:

$$[eqn]$$ $$[eqn]$$

5. First Convergence Check: If $[eqn]$ is accepted, then go to Step 6. Otherwise, go to Step 7.
6. Second Convergence Check: If $[eqn]$ is accepted, then go to Step 8. Otherwise, go to Step 7.
7. Looping: Set $[eqn]$ , where c is the number of increments for each iteration, and return to Step 2.
8. Upper and Lower Bounds: Compute an upper bound (UB) and a lower bound (LB) for each candidate solution, and select the best solution based on the smallest optimality gap between UB and LB.

The main idea behind determining the sample size is to choose the minimum number of samples that provide a highly accurate approximation. Although larger sample sizes yield better approximations, they also increase computational costs. Our proposed methodology computes the sample size $[eqn]$ at each iteration as $[eqn]$ for each $[eqn]$ , where $[eqn]$ and c are constants (input parameters). More details about the sample size and an investigation of the rate of convergence can be found in Homem-de-Mello and Bayraksan (2014) [52].

For Algorithm 2, we propose to generate a new set of realizations $[eqn]$ at each iteration instead of extending the realizations from the previous iteration.

As for the stopping criteria, we propose two methods, one related to the average cost and the other to the standard deviation of the cost.

For the first criterion, Algorithm 2 stops based on the result of the following statistical test:

$$[eqn]$$

For the second criterion, related to the standard deviation of the cost, Algorithm 2 uses the following statistical test:

$$[eqn]$$

Algorithm 2 continues iterating until both null hypotheses are accepted, indicating that further iterations are unlikely to result in significant improvements.

To select the best solution among R candidate solutions delivered in the last iteration of the proposed algorithm, we assess an upper bound and a lower bound for each candidate solution and consider the solution with the smallest optimality gap. For the selection of optimal solutions, we adopt a framework similar to the one used by Alnaggar et al. (2020) [56] for solving a two-stage transportation problem.

Let ( $[eqn]$ ) denote the vector solution of the first-stage problem for the $[eqn]$ candidate solution and $[eqn]$ denote its objective function value. To determine a lower bound of the objective function value, we start by computing the average $[eqn]$ and the variance $[eqn]$ of the objective function values $[eqn]$ .

$$[eqn]$$ $$[eqn]$$

Then the lower bound can be obtained as:

$$[eqn]$$

where $[eqn]$ is the $[eqn]$ critical value of the standard normal distribution, given that $[eqn]$ .

The upper bound of the true objective function for each potential solution is determined by examining the solution using a large scenario tree of size $[eqn]$ , assumed to accurately represent the actual cost distribution. Since each scenario $[eqn]$ is an independent, identically distributed (i.i.d.) random sample, the task of assessing a candidate solution breaks down into $[eqn]$ sub-problems. The size of the scenario tree $[eqn]$ significantly exceeds the size of the scenario tree handled in the last iteration of Algorithm 2, $[eqn]$ , where K is the index of the last iteration of the algorithm.

We represent the objective function value, $[eqn]$ , of a specific sub-problem i as:

$$[eqn]$$

It is important to highlight that, since each subproblem is addressed independently, a large value for $[eqn]$ does not impose a substantial computational load.

The estimation of the objective function value of the second-stage problem, denoted as $[eqn]$ , is calculated by:

$$[eqn]$$

The value of the objective function $[eqn]$ and its variance $[eqn]$ for candidate solution ( $[eqn]$ ) are computed as follows:

$$[eqn]$$ $$[eqn]$$

Finally, the upper bound, $[eqn]$ , of a candidate solution can be computed as:

$$[eqn]$$

where $[eqn]$ is the $[eqn]$ critical value of the standard normal distribution. The upper bound, $[eqn]$ , delivered by Algorithm 2 is the smallest value, $[eqn]$ for all $[eqn]$ , as defined in Equation (18). The solution that results in the smallest optimality gap, $[eqn]$ for all candidate solutions, $[eqn]$ , is the final solution, ( $[eqn]$ ), of the SAA, which corresponds to the solution with the smallest upper bound $[eqn]$ , is as shown in Equation (19).

$$[eqn]$$ $$[eqn]$$

5. Case Study: Monte Carlo Simulation for Two-Stage Bond Portfolio Optimization in the Saudi Sukuk Market

In this section, we apply Algorithm 2 to solve problem (9) using data from the Saudi Sukuk (bond) market. This case study demonstrates the practical implementation and effectiveness of the proposed methodology in a real-world financial context.

The data for this case study are sourced from the Saudi Stock Exchange (Tadawul) Annual Report of 2011, which provides detailed information on various tradable Sukuk and bonds (Afshar 2013 [57]; Alswaidan 2013 [58]; Alreshidi et al., 2020 [14]).

We select nine bonds for analysis, including SABIC 1, SABIC 2, SABIC 3, SIPCHEM, Saudi Electricity 1, Saudi Electricity 2, Saudi Electricity 3, SATORP, and Saudi Hollandi Bank 2. Key input parameters are summarized below (see Alreshidi et al., 2020 [14] for more details):

• −Issue size ranges from Saudi Riyal (SAR) 725 million for Saudi Hollandi Bank 2 to SAR 8000 million for SABIC 2, with a total of SAR 41,274 million.
• −Par value varies between SAR 10 thousand for SABIC 2, SABIC 3, SE 3 and SAR 500,000 for Saudi Electricity 1.
• −Maturity date varies, earliest with 6 July 2016 for SIPCHEM and latest 10 May 2030 for Saudi Electricity 3.
• −Coupon rates spread between SIBOR for Saudi Hollandi Bank 2 and SIBOR + 1.75% for SIPCHEM.
• −There is a total of 49 transactions, ranging from 2 for SABIC 1 and SABIC 3 to 14 for Saudi Electricity 2, while SATORP reported no transactions.
• −Value traded ranges from SAR 1097.4 thousand for SABIC 3 to about SAR 1.5 billion for Saudi Electricity 3, with a total of over SAR 1.8 billion.
• −Nominal value traded ranges from SAR 1.1 million for SABIC 3 to SAR 1.45 milion for Saudi Electricity 3, with a total of over SAR 1.8 billion.
• −Stage 1 interest rates at time periods t are 0.0148, 0.0317, 0.0112, 0.0112 for $[eqn]$ , respectively. The following assumptions are made for the bond portfolio optimization model and calculations of the input parameters are as presented by Alreshidi et al. (2020) [14]:

The total planning horizon is $[eqn]$ years, where the first stage covers $[eqn]$ years, while the second stage covers the remaining 6 years. We start with $[eqn]$ scenarios at each realization and increment by 30 at each iteration until the stopping criteria are met.

The interest rate, $[eqn]$ at time t, in both stages is calculated as $[eqn]$ , where SIBOR, the Saudi Interbank Offered Rate, is the benchmark interest rate used by banks in Saudi Arabia for lending to each other and is published daily by the Saudi Arabian Monetary Authority (SAMA), and A is a random variable generated uniformly in the range $[eqn]$ .

SABIC 1, SABIC 2, SABIC 3, SIPCHEM are considered to be Stage 1 bonds, and Saudi Electricity 1, Saudi Electricity 2, Saudi Electricity 3, SATORP, and Saudi Hollandi Bank 2 are considered to be Stage 2 bonds.

The cost, $[eqn]$ , of a bond j in Stage 1 and Stage 2 is considered the same as the par value in the Saudi data. Costs for SABIC 1, SABIC 2, SABIC 3, and SIPCHEM bonds are SAR 50,000, SAR 10,000, SAR 10,000, and SAR 100,000, respectively. Costs for Saudi Electricity 1, Saudi Electricity 2, Saudi Electricity 3, SATORP, and Saudi Hollandi Bank 2 bonds are SAR 500,000, SAR 100,000, SAR 10,000, SAR 100,000, and SAR 100,000, respectively.

The buying price, $[eqn]$ , of a bond $[eqn]$ under scenario $[eqn]$ in the second-stage is calculated as $[eqn]$ , where $[eqn]$ represents the purchase price of bond $[eqn]$ , which corresponds to its nominal buying price. Additionally, $[eqn]$ serves as an impact factor for each scenario k, where $[eqn]$ , in Stage 2.

The selling price, $[eqn]$ , of a bond $[eqn]$ in Stage 2 under scenario $[eqn]$ for period $[eqn]$ is obtained as follows:

$$[eqn]$$

where all parameters are as defined in Table 1.

The cash flow, $[eqn]$ , for bond $[eqn]$ in period $[eqn]$ in Stage 1 is calculated as below:

$$[eqn]$$

where $[eqn]$ denotes the purchase price of bond $[eqn]$ in Stage 1, $[eqn]$ represents the coupon rate, and $[eqn]$ signify the maturity date of any bond $[eqn]$ . The cash flows for Stage 1 bonds, calculated using Equation (21), indicate that SABIC 1, SABIC 2, SABIC 3, and SIPCHEM generate annual amounts of SAR 1575, SAR 313, SAR 323, and SAR 4500, respectively, which remain constant over the four-year horizon.

The cash flows, $[eqn]$ , for each bond $[eqn]$ in Stage 2 under scenario $[eqn]$ are calculated using the equation below:

$$[eqn]$$

where $[eqn]$ represent the purchase price of bond j acquired in Stage 2 under scenario k and $[eqn]$ and $[eqn]$ represent the coupon rate and the maturity date of bond j in Stage 2, respectively.

Stage 1 liabilities, $[eqn]$ , $[eqn]$ , are a percentage (randomly chosen in the interval [2%, 5%]) of the total cash flows from available bonds at the beginning of Stage 1, where $[eqn]$ = 1,009,895, $[eqn]$ = 2,275,167, $[eqn]$ = 3,064,511, and $[eqn]$ = 4,086,015.

Stage 2 liabilities, $[eqn]$ , are randomly chosen in the interval [0.9B, 1.1B] for the second stage, where B is a percentage (randomly chosen in the interval [0.2%, 2%] of the total cash flows from available bonds at the beginning of the second stage.

After solving the two-stage bond-portfolio optimization problem (9) using the Monte Carlo simulation-based approach described in Algortihm 2, we obtain the optimal costs for the Saudi Sukuk (bond) data. This is achieved by implementing the algorithm in C++ with CPLEX. The optimal costs of the two-stage bond-portfolio optimization problems for all realizations and different iterations of the algorithm are presented in Table 3, Table 4 and Table 5.

Table 6 presents the average cost and variance across different iterations. For each iteration, $[eqn]$ realizations were considered, and the average cost and variance were computed based on these 30 realizations.

Figure 1 illustrates the average cost for different sample sizes (number of scenarios) $[eqn]$ to $[eqn]$ . We observe that, as the sample size (number of scenarios) increases, the average costs stabilize, with the values converging closer to each other. Additionally, Figure 2 shows the variance of the average cost, revealing that variance decreases as the sample size (number of scenarios) increases.

At each iteration, the statistical tests for both hypotheses were conducted based on the computed sample averages and variance estimates obtained through Algorithm 2. Table 6 presents the average cost (in Saudi Riyal SAR) and variance for different sample sizes (number of scenarios) ranging from $[eqn]$ to $[eqn]$ scenarios for problem (9).

Based on the results obtained for problem (9) shown in Table 7, we observe that as the number of scenarios, K, increases, the average cost initially fluctuates, with noticeable differences in earlier iterations. This variation suggests that the optimization process is still adjusting to increase sample sizes. However, as the number of scenarios, K, continues to grow, the differences in mean and variance become statistically insignificant since $[eqn]$ , especially after $[eqn]$ scenarios. The F-statistic confirms variance stabilization, indicating that cost fluctuations across realizations have reached a consistent state. At $[eqn]$ the stopping criterion is met, as both the T-statistic and F-statistic show that further increases in the number of scenarios, K, will not yield statistically significant improvements in cost reduction or stability. Hence, additional iterations would not be computationally efficient, and therefore, $[eqn]$ , is the optimal stopping sample size.

To investigate the optimal cost, we focus on the optimal sample size, $[eqn]$ , representing the optimal number of scenarios. To select the best solution among the 30 realizations, we calculate the lower bound, $[eqn]$ , in Equation (12) and the upper bound, $[eqn]$ , in Equation (18) to identify the smallest optimality gap, as described in Section 4. Table 8 indicates that the smallest gap among the 30 realizations is 1.1438, with the optimality range between [SAR 22,652,406 and SAR 22,914,500], where SAR is Saudi Riyal.

Based on the chosen solution, we determine the bonds to be purchased in Stage 1. Specifically, 3979 units of SIPCHEM and 478 units of SABIC 3 are purchased in period 1 of Stage 1, as shown in Table 9.

Table 10 presents the optimal investment plan for Stage 1, based on the modified Model. The total optimal cost of the investment plan, according to the modified Model 1 is SAR 22,637,904.

Table 11 illustrates the computational efficiency of Algorithm 2 based on varying sample sizes (number of scenarios). All computational experiments were carried out on a desktop with Intel(R) Core (TM) i74930K CPU 3.4 GHz processor with 16 GB of memory under a Windows environment. As expected, the computation time increases as the number of scenarios increases. We observe that the computational cost has increased by 50 times when the sample size increases from 1000 scenarios to 2000 scenarios, whereas it significantly decreases for 5000 scenarios. We remark that optimization solvers dynamically adjust their internal methods based on the input data. The decrease in computational cost for the problem with 2000 scenarios to the problem with 5000 scenarios may be a result of more efficient numerical methods or multithreading optimization that are triggered at $[eqn]$ scenarios but were not activated at $[eqn]$ scenarios.

In summary, the results of applying Algorithm 2 along with the proposed stopping criteria to the Saudi Sukuk (bond) data demonstrate the proposed method’s effectiveness and computational efficiency in solving large-scale stochastic optimization problems. The proposed algorithm, Algorithm 2, showcases remarkable computational efficiency.

6. Conclusions

This paper presents a Monte Carlo simulation-based approach for solving two-stage stochastic programming problems, specifically applied to bond portfolio optimization. By leveraging random sampling techniques, the proposed method effectively transforms the stochastic problem into a deterministic one, enabling the use of standard optimization techniques. The application of this approach to the Saudi Sukuk (bond) market provides valuable insights and demonstrates the practicality and robustness of the methodology.

The key contributions of this research are many-fold: (1) The integration of Monte Carlo simulation with two-stage stochastic programming offers a novel and efficient way to address uncertainty in bond portfolio optimization. This approach allows for accurate approximation of expectation functions and effective scenario management. (2) The robustness of the approach in dealing with market uncertainties ensures that optimized portfolios remain resilient across different scenarios. (3) The proposed algorithm showcases remarkable computational efficiency and manages to handle large scenario trees and deliver solutions in a reasonable time frame. This efficiency is vital for practical applications, where timely decision-making is crucial. (4) The optimization results indicate substantial cost reductions and effective liability management. Algorithm 2 consistently identifies optimal bond portfolios across various realizations and iterations. By employing the Monte Carlo simulation approach, Algorithm 2 converts the stochastic problem into a deterministic one, making it feasible to solve with standard optimization techniques. (5) One of the significant strengths of the Monte Carlo simulation approach is its ability to handle market uncertainties. The scenario analysis reveals that Algorithm 2 can effectively optimize bond portfolios under varying interest rates and liability conditions. The robustness of the proposed method ensures that the optimized portfolios remain resilient against different market scenarios, thereby minimizing risks and maximizing returns. (6) To ensure the reliability of the results, statistical tests are employed to validate the convergence of the algorithm. The use of sample average approximation (SAA) allows for accurate approximation of the expectation function. The convergence criteria, based on statistical significance tests for average costs and variances, confirm the stability and precision of the solutions. The upper and lower bounds computed for each candidate solution further are validated by the optimality of the results. (7) The practical application of this methodology to the Saudi Sukuk (bond) market demonstrates its real-world relevance. The case study highlights the adaptability of the Monte Carlo simulation approach to different financial instruments and market conditions. The results provide valuable insights for investors and portfolio managers, enabling them to make informed decisions that optimize bond portfolios while effectively managing risks. (8) Incorporating uncertainties such as fluctuating interest rates and liabilities in bond portfolio optimization, the proposed approach allows for better risk management and more resilient investment strategies. The methodology can be applied to other financial markets and instruments, making it a versatile tool for global investment strategies.