A Backward Simulation Method for Stochastic Optimal Control Problems
Zhiyi Shen, Chengguo Weng
TL;DR
This paper introduces a generalized backward simulation algorithm for stochastic optimal control problems, improving computational efficiency and accuracy over traditional methods, demonstrated through an insurance product pricing application.
Contribution
The paper extends the Least-Squares Monte Carlo algorithm to a broader class of stochastic control models with novel simulation and estimation techniques.
Findings
The algorithm bypasses forward simulation and control randomization.
It avoids extrapolating the value function.
It reduces computational burden in parameter tuning.
Abstract
A number of optimal decision problems with uncertainty can be formulated into a stochastic optimal control framework. The Least-Squares Monte Carlo (LSMC) algorithm is a popular numerical method to approach solutions of such stochastic control problems as analytical solutions are not tractable in general. This paper generalizes the LSMC algorithm proposed in Shen and Weng (2017) to solve a wide class of stochastic optimal control models. Our algorithm has three pillars: a construction of auxiliary stochastic control model, an artificial simulation of the post-action value of state process, and a shape-preserving sieve estimation method which equip the algorithm with a number of merits including bypassing forward simulation and control randomization, evading extrapolating the value function, and alleviating computational burden of the tuning parameter selection. The efficacy of the…
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6| Parameter | Value |
|---|---|
| Volatility rate | |
| Risk-free rate | |
| Insurance fee rate | |
| Number of time steps | 12 |
| Length of time interval | 1/12 |
| Discouting factor | 0.9975 |
| Initial purchase payment | 1 |
| Withdrawal penalty | 0.8 |
| Guaranteed payment percentage | |
| Setting | SPSE | RSE | |||
|---|---|---|---|---|---|
| Mean | S.d. | Mean | S.d. | ||
| 0 | |||||
| 1 | |||||
| 2 | |||||
| 3 | |||||
| 4 | |||||
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsInsurance, Mortality, Demography, Risk Management · Stochastic processes and financial applications · Financial Risk and Volatility Modeling
A Backward Simulation Method for
Stochastic Optimal Control Problems 111The authors are thankful to Alexander Schied, Pengfei Li, Degui Li, and Zhaoxing Gao for helpful discussions.
Zhiyi Shen 222Department of Statistics and Actuarial Science, University of Waterloo. Email Address: [email protected], and Chengguo Weng 333Department of Statistics and Actuarial Science, University of Waterloo. Email Address: [email protected]
Abstract
A number of optimal decision problems with uncertainty can be formulated into a stochastic optimal control framework. The Least-Squares Monte Carlo (LSMC) algorithm is a popular numerical method to approach solutions of such stochastic control problems as analytical solutions are not tractable in general. This paper generalizes the LSMC algorithm proposed in Shen and Weng (2017) to solve a wide class of stochastic optimal control models. Our algorithm has three pillars: a construction of auxiliary stochastic control model, an artificial simulation of post-action value of state process, and a shape-preserving sieve estimation method which equip the algorithm with a number of merits including bypassing forward simulation and control randomization, evading extrapolating the value function, and alleviating computational burden of the tuning parameter selection. The efficacy of the algorithm is corroborated by an application to pricing equity-linked insurance products.
1 Introduction
The stochastic optimal control model is a prevalent paradigm for solving optimal decision problems with uncertainty in a variety of fields, particularly, financial engineering. In solving discrete-time stochastic optimal control problems, the Dynamic Programming Principle (DPP) is a prevailing tool which characterizes the optimal value function as the solution to a backward recursive equation system, often known as the Bellman equation. This reduces the stochastic optimization problem into two separate problems:
- solving a sequence of deterministic optimization problems and 2) evaluating the conditional expectation terms in the Bellman equation. In spite of the theoretical appealingness of the DPP, there generally does not exist closed-form solution of the Bellman equation, which impedes the application of stochastic optimal control models to complicated real-world problems. Recently, a number of numerical methods have been proposed in the literature to approach the optimal or suboptimal solutions to various stochastic optimal control problems by combining Monte Carlo Simulation with nonparametric regression methods.
In a statistical setting, the typical goal of nonparametric regression methods is to estimate the functional form of the expectation of a response variable conditioning on a covariate variable. This naturally motivates one to use certain nonparametric regression methods to evaluate the conditional expectation (also known as the continuation value in the context of pricing Bermudan option) involved in the Bellman equation where the value function at the next time step and the state variable at the current time step are taken as the response and covariate variables, respectively. Such a ground-breaking idea was incubated in a series of papers including Carriere (1996), Longstaff and Schwartz (2001), and Tsitsiklis and Van Roy (2001), and the corresponding numerical algorithms are often referred to as the Least-Squares Monte Carlo (LSMC) algorithms. Since then, the LSMC algorithm has witnessed remarkable popularity in solving optimal stopping problems, a special class of stochastic control problems; see, e.g., Clément et al. (2002), Stentoft (2004), Egloff (2005), Glasserman and Yu (2004), Glasserman et al. (2004), Egloff et al. (2007), Zanger (2009), Zanger (2013), Belomestny et al. (2009), Belomestny (2011), and the references therein.
The problem of solving general stochastic control problems by resorting to the LSMC algorithm is considerably more involved. To understand the crux, let us note that the LSMC method has two building blocks: 1) a forward simulation of the state process and 2) a backward updating procedure which employs the nonparametric regression to estimate the continuation value. In an optimal stopping problem, the evolution of the state process is independent of decision maker’s (DM’s) action and therefore, the forward simulation of the sample paths of the state process is relatively straightforward. In stark contrast to this, in a general stochastic optimal control setting, the state process is influenced by the DM’s action and accordingly, its simulation is unattainable without specifying the DM’s action. Ideally, one may expect to simulate the state process driven by the optimal action of the DM. However, the optimal action should be determined by solving the Bellman equation in a backward recursion manner, which is incongruous with the need of forward simulation in an LSMC algorithm. To circumvent this, Kharroubi et al. (2014) proposes to first draw the DM’s action from a random distribution and then simulate the sample paths of the state process based on the initialized action. This method has been referred to as the control randomization method in the literature and applied in the LSMC algorithm to solve many specific stochastic control problems; see, e.g., Cong and Oosterlee (2016), Zhang et al. (2018), and Huang and Kwok (2016), among others. Despite the wide usage of the control randomization method, the accuracy of the numerical estimate is impaired over the region with sparse sample points (Zhang et al., 2018, Section 3.3) and the spread of the sample paths is sensitive to the specific way of initializing the action. In an extreme case, the LSMC algorithm might even miss the optimal solution under a dismal choice of random distribution from which the PH’s action are drawn; see Shen and Weng (2017), for instance. It is also notable that most literature bind together the control randomization and the forward simulation of the state process. However, this paper will show that the forward simulation is not imperative in an LSMC algorithm and the merits of abjuring the forward simulation are extant in several aspects. The limitations of the control randomization and the consequential forward simulation will be elaborated in Section 2.2 of this paper.
Besides the simulation of the state process, the approximation of the conditional expectation term in a Bellman equation is also taxing for several reasons. Firstly, the prevalent regression methods only warrant the accuracy of the regression estimate over a compact support, see, e.g., Newey (1997), Stentoft (2004), and Zanger (2013), whereas the state variable generally takes value in an unbounded set. Some literature compromise to first truncate the domain of the continuation function and then use extrapolation techniques when the knowledge of the function outside the truncated region is required. It is worth noting that this problem is not acute in the context of optimal stopping problem but is severe in a general stochastic control setting. This is because, in the latter case, one has to traverse all admissible actions, which calls for the values of the continuation function over a domain that is wider than spreading range of sample paths. Secondly, in order to avoid overfitting or underfitting, most nonparametric regression methods thirst for an appropriate choice of the tuning parameter, e.g., the number of basis functions in a linear sieve estimation method (see the sequel Section 3.3). This is often resolved by computationally expensive cross-validation methods, see, e.g., Li (1987). However, in view of the extraordinarily large number of simulated paths, such a tuning parameter selection procedure is intolerable in implementing the LSMC algorithm. The aforementioned challenges will be investigated in details in the sequel section.
The contribution of this paper is summarized as follows. Firstly, we restrain the value set of the state process into a compact set, which evades the undesirable extrapolating value function estimate during the backward recursion of the LSMC algorithm. The value function accompanying the truncated state process is shown to be a legitimate approximation for the primal value function under a suitable choice of the truncation parameter. Secondly, we generalize the idea of Shen and Weng (2017) to simulate the post-action value of the state process from an artificial probability distribution. This eliminates the need for the forward simulation and is consistent with the backward induction nature of the Bellman equation. The memory as well as time costs of the artificial simulation method are considerably less than those of the control-randomization-based forward simulation method. Thirdly, we introduce a shape-preserving sieve estimation method to approximate the conditional expectation term involved in the Bellman equation. By exploiting certain shape information of the continuation function, the sieve estimate is insensitive to the tuning parameter and accordingly reduces the computational cost of the tuning parameter selection. We refer to the proposed LSMC algorithm as the Backward Simulation and Backward Updating (BSBU) algorithm. Finally, we establish the convergence result of BSBU algorithm which sheds light on how the numerical error propagates over the backward recursion procedure.
This paper is organized as follows. Section 2 gives a tour through the LSMC algorithm and shows its challenges in solving general stochastic optimal control problems. Section 3 gives the main results of the paper: a construction of auxiliary stochastic optimal control model, the BSBU algorithm, and the associated convergence analysis. Section 4 applies the BSBU algorithm to the pricing problem of an equity-linked insurance product and Section 5 conducts the corresponding numerical experiments. Finally, Section 6 concludes the paper.
2 Basic Framework and Motivations
2.1 Stochastic Optimal Control Model
We restrict our attention to a collection of consecutive time points labeled by on which a decision maker (DM) may take action. The uncertainty faced by the DM is formulated by a probability space equipped with a filtration \mathbb{F}=\big{\{}\mathcal{F}_{t}\big{\}}_{t\in\mathcal{T}}. The DM’s action is described by a discrete-time stochastic process with . Let be a certain state process valued in with . Starting from an initial state , it evolves recursively according to the following transition equation:
[TABLE]
where is a sequence of independent random variables valued in with . reflects the uncertainty faced by the DM at time step and is referred to as random innovation in what follows. For brevity of notation, in what follows, we compress the dependency of the state process on the action and the readers should always bear in mind that implicitly depends on the DM’s action up to time . We give the formal definition of the action as follows.
Definition 1** (Admissible Action).**
We call a discrete-time process an admissible action if it satisfies:
(i)
* is -measurable for ;*
(ii)
* for , where is some function valued as a subset of with .*
The function in the preceding definition corresponds to a certain state constraint on the action taken at time . Denote by the set of the DM’s admissible actions. Consider a discrete-time stochastic optimal control problem in the following form:
[TABLE]
where is a certain discounting factor, and are the intermediate and terminal reward functions, respectively. In order to ensure the well-posedness of the stochastic control problem (2), we impose the following assumption which is conventional in literature, see Rogers (2007) and Belomestny et al. (2010) for instance.
Assumption 1**.**
[TABLE]
The Dynamic Programming Principle states that the value function can be solved recursively:
[TABLE]
where
[TABLE]
We proceed by rewriting the transition equation (1) into the following form:
[TABLE]
where and are some measurable functions with . It is worth stressing that any transition function can be rewritten into the above form since one may choose as identity function (i.e., ) and the above equation holds trivially. Nevertheless, it is instructive to introduce the function as it brings the benefit of dimension reduction. We will explain this more in the sequel. Combing Eqs. (4) and (5), we get
[TABLE]
Hereafter, we call the post-action value of the state process at time . It constitutes an essential component in the LSMC algorithm proposed in Section 3.2. Define function
[TABLE]
We observe the following relationship between and :
[TABLE]
The crucial implication of the above relation is that it suffices to recover the functional form of in order to evaluate since is known at the first hand. The motivation of rewriting the transition equation into Eq. (5) is now clear: maps a -dimensional vector into a -dimensional vector, which compresses the dimension if and it is more efficient to recover the function than due to such a dimension reduction. It is also worth noting that is solely determined by the probability distribution of according to Eq. (6). This implies that it is not necessary to know the exact distribution of in the evaluation of the function . In view of the relation (7), the Bellman equation (3) can be equivalently written as
[TABLE]
The above equation system states that, given the value function at time step , one may first evaluate continuation function according to Eqs. (6) and (7) and then obtain the value function at time step via solving an optimization problem in the second line of Eq. (8). The information propagation behind the above recursive procedure is illustrated in Figure 1.
2.2 A Tour Through LSMC Algorithm
We proceed by briefly reviewing the Least-squares Monte Carlo (LSMC) algorithm. We will show its limitations in several aspects which motivate the algorithm we will propose in the subsequent sections.
2.2.1 “Forward simulation and backward updating” (FSBU) algorithm
There has been voluminous literature on the LSMC for optimal stopping problem, while the literature on the LSMC for general stochastic optimal control problem is thin. Most literature addresses the LSMC for the stochastic control problems arising in some specific applications, see, e.g., Carmona and Ludkovski (2010), Barrera-Esteve et al. (2006), Huang and Kwok (2016), Shen and Weng (2017), Cong and Oosterlee (2016), and Zhang et al. (2018), among others. An LSMC algorithm for a class of stochastic control problem is developed in Belomestny et al. (2010).
For most variants of the LSMC algorithm, they can be decomposed into two pillars: (i) a forward simulation of the state process and (ii) a backward updating of control policies. We review these algorithms in a unified paradigm as follows.
Initiation: Set . For , do the two steps below. 2. 2.
Forward Simulation:
- 2.1
Control randomization Generate a random sample of the DM’s action up to time step :
[TABLE]
with each generated by a certain heuristic rule. 2. 2.2
Simulation of state process Simulate a random sample of the random innovations:
[TABLE]
The sample of the state process up to time step is given by
[TABLE]
where for . 3. 3.
Backward Updating
- 3.1
Regression Given a numerical estimate of value function at time step , denoted by , construct the random sample
[TABLE]
Further construct a random sample of post-action value of the state process as follows:
[TABLE]
Take and as the samples of response variable and regressor, respectively, and employ a certain non-parametric regression to obtain a regression estimate for . 2. 3.2
Optimization An estimate for the value function at time step is given by
[TABLE]
We henceforth call the above algorithm as the Forward Simulation and Backward Updating (FSBU) algorithm.
Remark 1** (Randomness of and ).**
The superscript E in and stresses that they are numerical estimates of the true value function and continuation function, respectively. Since a certain regression technique is employed to get such numerical estimates, they essentially depend on the random samples and and hence on all previously generated random samples going from step down to step , i.e., and for . Such dependency is suppressed in notation for brevity, but the readers should keep in mind that both and are random functions.
2.2.2 Challenges
There are several challenges in implementing the previously introduced FSBU algorithm to solve a stochastic control problem. Below, we make some comments on the challenges from three aspects.
(I)
Limitation of control randomization As the DM’s optimal action is not tractable at the very first but should be solved in the backward updating stage of the algorithm, Step 1.1 randomly generates a feasible action, which is referred to as control randomization method; see Kharroubi et al. (2014). For some selected action , the accuracy of the regression estimate can be warranted only over the support of the resulting sampling points , say , which might be smaller than those for other actions. However, in order to solve the optimization problem in Step 2.2 (see Eq. (10)), one requires the knowledge of over the range of post-action value for all feasible actions because all possible values of the action should be traversed and taken as the input of the function in evaluating ; see Figure 2 for a graphical illustration. As a compromise, one may use certain extrapolation methods to infer the value of outside the region , which incurs extra error and is hard to justify its legitimacy.
(II)
Cost of forward simulation It is notable that, at time step of the above FSBU algorithm, a new random sample of the state process that is independent of the sample at time step is simulated; see Figure 3 for a graphical illustration. This is required in order to apply the nonparametric regression theory to establish the convergence result, see Section 2.3 of (Zanger, 2013, p. 511) for instance. On the contrary, using a single sample causes in-sample bias because the numerical estimate of value function obtained at time step is correlated with ; see, e.g., Section 3.1 of Choi et al. (2018) and the earlier Remark 1. The total time cost in a forward simulation procedure of the above LSMC algorithm is of 444Suppose the time cost of simulating a path over each time interval is . Then the time cost of simulating a whole path up to time step is about , and the forward simulation in the whole LSMC algorithm has approximate time cost of , accordingly.. Simulating the whole path of state process can be time-consuming especially when one uses some approximation schemes to simulate general stochastic differential equations555It is the authors’ experience that a single simulation of paths of the Heston model over a -year period takes seconds by using the R package “msde” on a MacBook Pro (2.8 GHz Intel Core i7).. Besides the issue of time cost, the memory cost in a single simulation is of with and being the number of time steps and dimensionality of the state process, respectively, which is sizable for a large .
(III)
Choice of regression technique Despite the voluminous literature on nonparametric regression, the choice of the nonparametric regression method in Step 2.1 should be meticulous. In the above FSBU algorithm, the sample size in the regression problem corresponds to the number of simulated paths and is generally recommended in the literature be chosen larger than one hundred thousand, which makes most regression methods computationally prohibitive. To be specific, the local methods such as local-polynomial regression are clearly not wise choices as they require running a regression at each sample point. It is worthy to point out that even computing a single point in the sample is fairly time-consuming as it involves a local optimization problem (see Eq. (10)). Furthermore, the nuisance of high memory cost also burdens most nonparametric regression methods. For example, the kernel regression and Isotonic regression methods require storing all sample points in order to recover the functional form of the regression function over some support, and the memory cost is extraordinarily large， accordingly. The above two thorny issues are escalated by the fact that almost all nonparametric regression techniques involve a computationally-intensive cross-validation procedure to determine the tuning parameter (e.g., the bandwidth in local regression methods and the number of basis functions in global regression methods) in order to avoid overfitting or underfitting.
2.2.3 Motivations
In view of the previous items (I)–(III), the thrust behind this paper is to explore possible answers to the following questions:
(Q1)
How to avoid theoretically shaky extrapolation?
(Q2)
Is it possible to bypass the forward simulation in an LSMC algorithm?
(Q3)
Is there a regression method that is insensitive to tuning parameter?
In terms of (Q1), in the sequel section, we will construct an auxiliary stochastic control problem where the accompanying state process only takes values in a bounded set. This construction sidesteps extrapolating the regression function outside the region where the sample distributes. In response to (Q2), we will propose to directly simulate the post-action value of state process. For (Q3), we will introduce a shape-preserving sieve estimation method to infer the continuation function. The resulting sieve estimate, on the one hand, is insensitive to the tuning parameter, and on the other hand, preserves certain shape properties of the continuation function.
3 Main Results
3.1 Localization and Error Bound
As commented in the item “Limitation of control randomization” in the last section, it is necessary to know the value of the continuation function over the whole range of post-action value of the state process which is wider than the set where the regression sample suffuse. It is notable that the range of post-action value is unbounded if the state process takes value in an unbounded set, which is particularly the case in many finance applications. Therefore, it is generally inevitable to infer the continuation function outside the support of the sample and the error incurred by extrapolating the regression estimate is hard to quantify.
The aim of this subsection is to find a certain way to circumvent the unsound extrapolation in the implementation of an LSMC algorithm. The pivotal idea is to first construct an auxiliary stochastic optimal control problem where the accompanying state process only takes values in a bounded set and then show the discrepancy between the auxiliary problem and the primal one is marginal in a certain sense. To formalize the idea, we let be a bounded subset of the set where the subscript denotes a certain truncation parameter. Further denote (resp. ) as the interior (resp. boundary) of . Given the initial state , define the following stopping time:
[TABLE]
with the convention: if for all . Let be the closure of the set and assume it to be strictly convex.
We recursively define an auxiliary state process as follows:
[TABLE]
where with denoting the Euclidean -norm. Since is a compact and strictly convex set, is unique and lies on the boundary set for .
Below we give some interpretations regarding the auxiliary state process defined in the above Eq. (12). The primal state process coincides with the auxiliary state process until the stopping time . Once the primal state process passes through the interior of the truncated domain, the auxiliary state process freezes at a certain point in the boundary set thereafter. The evolution mechanisms of the primal and auxiliary state processes are illustrated in Figure 4. The following proposition gives the transition equation of .
Proposition 1**.**
The auxiliary state process defined by Eq. (12) admits the following transition equation across each time point: and
[TABLE]
where
[TABLE]
and is the transition equation relating the pre-action and post-action values of the primal state process defined in Eq. (5).
The proof of the above proposition is relegated to Appendix B.1 for clarity of presentation. The preceding Eq. (13) essentially states that is a Markov chain by itself, and accordingly, it is the sole state process of the auxiliary stochastic control model defined in the sequel.
Let be the set of all admissible actions for the auxiliary state process which is defined as:
[TABLE]
Relative to the primal stochastic optimal control problem (2), we consider the following auxiliary problem:
[TABLE]
where is defined recursively by Eq. (13) for any given action . Since the state process freezes once it reaches the boundary set , the value function in Eq. (15) is given by
[TABLE]
with . Over the interior of the truncated domain, the above value function can be solved in a similar backward recursion way as does, that is,
[TABLE]
where is defined in line with Eq. (6) with replaced by . It is worth noting that, in evaluating \tilde{C}_{t}\big{(}K(x,a)\big{)}, the knowledge of over might be in need, and in such a situation, Eq. (16) is invoked.
We make some comparisons between Eq. (17) and the Bellman equation (8) associated with the primal stochastic control model. Firstly, in both equations, the state constraint , transition equation between pre-action and post-action values , and reward functions are exactly the same. Secondly, the value function is solely defined on a bounded set , whilst is defined on the set which might be unbounded in many financial applications as the primal state process may correspond to a certain risky asset valued on the whole positive real line.
We will characterize the discrepancy between the value functions and in a certain sense. To this end, it is necessary to impose some assumptions on the state process and the reward functions.
Assumption 2**.**
Let . There exists a measurable function satisfying
[TABLE]
in Eq. (18) gives an upper bound for the probability that the auxiliary state process disagrees with the primal at some time before maturity regardless of the DM’s action. Since the primary difference between the auxiliary and primal value functions stems from the disparity between the associated state processes, it is not surprising that the above inequality (18) plays an important role in characterizing the approximation error of as we will see later in the proof of Theorem 1. The expression of should be specified for each specific application.
Assumption 3**.**
(i)
There exists a measurable function and a generic constant independent of and such that \big{|}f_{T}(x)\big{|}^{2}\leq B(x), \sup_{a\in A_{t}(x)}\big{|}f_{t}(x,a)\big{|}^{2}\leq B(x),
[TABLE]
(ii)
There exists a measurable function satisfying
[TABLE]
In many applications has a polynomial form and in such a situation, the above assumption states that the reward functions are bounded by a certain polynomial from the above uniformly in . In the context of pricing financial products, this assumption says that the policy payoffs have a polynomial growth rate.
The following theorem quantifies the error stemming from using the auxiliary problem (15) as a proxy for the primal stochastic control model (2).
Theorem 1** (Truncation Error Estimate).**
Suppose Assumptions 1, 2, and 3 hold. Then
[TABLE]
The proof of the above theorem is relegated to Appendix B.2. The inequality (20) can be understood as follows. The term \big{(}\xi(R)+\zeta\big{)}\mathcal{E}(X_{0},R) corresponds to an upper bound for the discrepancy between the reward functions of the two stochastic control models (2) and (15) at each time step. Since such a difference primarily stems from replacing the primal state process by , it is not surprising that the term appears in the error estimate. Furthermore, the two terms and correspond to certain upper bounds of the magnitudes of the reward terms and , respectively, and therefore a square root arises in the inequality (20). Finally, the discrepancy between the two value functions is amplified as the time horizon is prolonged, which is reflected by the existence of factor in the above error estimate.
3.2 A Backward Simulation and Backward Updating Algorithm
3.2.1 Simulation of post-action value
In this subsection we propose an LSMC algorithm which simulates the state process without referring to the optimal action. Recall from Step 2.1 of the FSBU algorithm in Section 2.2 that the ultimate goal of simulating the state process is generating a random sample of the post-action value of the state process which acts as a crucial input for the regression step. This naturally inspires us to directly simulate the post-action value from an artificial probability distribution. The term “artificial” stresses the fact that such a distribution might not coincide with the distribution of under the optimal action process.
Since the value function is explicitly given by Eq. (16) over , the primary goal of our proposed LSMC algorithm is to get a numerical estimate for the value function over the open set . In view of this, we may circumscribe the support of the artificial probability distribution that the post-action values are simulated from. First note that the range of post-action value of the auxiliary state process denoted by is given by
[TABLE]
Consider the following subset:
[TABLE]
for , where is the set of all values the random innovation might take and is the transition equation relating the post-action value at time step to the state variable at time step which is given in Eq. (14).
The preceding equation states that will stop at a certain point in the boundary set if lies in the set ; see Figure 5 for a graphical illustration. To make the matter more concrete, let us consider the example of pricing variable annuities (see, e.g. Huang and Kwok (2016) and Shen and Weng (2017)) where corresponds the post-withdrawal value of the investment account. If the investment account is depleted after the policyholder’s withdrawal (i.e., ), it remains exhausted forever (i.e., for ). In such an example, is a singleton . In view of the above discussion and Eq. (16), for any , we observe
[TABLE]
which has a value independent of and is given by Eq. (16). Therefore, at time step , it suffices to get a regression estimate for the continuation function on the set .
3.2.2 The algorithm
Now we present the Backward Simulation and Backward Updating (BSBU) algorithm as follows.
Initiation: Set for . For , do the two steps below. 2. 2.
Backward Simulation:
- 2.1
Simulation of post-action value Generate a sample of the post-action values denoted by
[TABLE]
from a probability distribution with support . 2. 2.2
Simulation of the state process Construct the sample of the state process at time step according to
[TABLE]
with being a sample of the random innovations. 3. 3.
Backward Updating:
- 3.1
Data preparation Given a numerical estimate of value function at time step , denoted by , construct the sample
[TABLE] 2. 3.2
Regression Take and as the samples of response variable and regresssor, respectively, and employ a certain non-parametric regression to obtain a regression estimate over the set . For , we set with given by Eq. (22). 3. 3.3
Optimization An estimate for the value function at time step is given by:
[TABLE]
For , we set with given by Eq. (16).
In Step 3.2, we prescribe for the value of when because might fall in the set . Similarly, in Step 3.3, Eq. (16) is invoked to evaluate for as generated by Eq. (23) may lie on , the boundary set of the truncated domain. The backward information propagation in the above BSBU algorithm is illustrated in Figure 6.
3.2.3 Discussions
Comparing the above BSBU algorithm and the FSBU counterpart in Section 2.2, we have the following observations.
Firstly, the primary difference of the two algorithms lies in how to generate the post-action values of state process, i.e., . The control randomization method is a forward simulation scheme while the BSBU algorithm directly generates post-action value from a certain prior distribution. Indeed, the FSBU algorithm can be viewed as a special BSBU algorithm if is chosen as the probability distribution of the post-action value from a control randomization procedure. In general, both methods do not yield the distribution of driven by the optimal action, and thus, there is no loss to directly generate from a prior distribution . 2. 2.
Secondly, the BSBU method has the advantage of reducing memory and time costs. On the one hand, one does not need to store the sample of whole trajectories at each time step in the BSBU algorithm. On the other hand, the total time cost of simulating the state process is of in the BSBU algorithm, while it is of in the FSBU counterpart; see the item “Cost of forward simulation” in Section 2.2. 3. 3.
Thirdly, the BSBU algorithm circumvents extrapolating the numerical estimates of the continuation function and the value function. It is notable that and are obtained over the sets and , respectively, at time step ; see Steps 2.2-2.3 of the above BSBU algorithm. At the time step , the BSBU algorithm does not require the knowledge of the value function (resp. continuation function) outside (resp. ) in computing the regression data ; see Figure 7 for a graphical illustration. This nice property inherits from the construction of the auxiliary state process whose values are confined to a bounded set. In the FSBU algorithm, however, the state process is not restrained and might lie outside the regression domain for . In such a situation, extrapolating the numerical solution causes extra error which is hard to control.
3.3 Sieve Estimation Method
In this subsection, we discuss the details of the regression method used to estimate the continuation function in our BSBU algorithm.
3.3.1 Selection criteria for regression method
In Section 2.2 we have discussed the potential issue which may be associated with a regression method when used for the estimation of the continuation function of a stochastic control problem; see the item “Choice of regression method”. Based on the discussion, we propose the following criteria for the choice of regression method in estimating continuation function.
(C1)
Small memory cost The regression problem embedded in an LSMC algorithm usually exhibits extraordinarily large sample size. Thus, an appropriate regression method should have small memory requirement. This criterion excludes the kernel method (Nadaraya (1964) and Watson (1964)), local-polynomial regression method (Fan and Gijbels (1996)), and Isotonic regression method (Robertson et al. (1988)) which require storing all sample points in the memory in order to compute the regression function at any point in the domain.
(C2)
Computationally cheap In almost all nonparametric regression methods, a certain parameter (referred to as tuning parameter in statistical literature) is used to avoid undesirable overfitting or underfitting of the regression model. Determining the optimal value of such a tuning parameter is usually computationally intensive. Therefore, an ideal regression method should be insensitive to the tuning parameter.
In view of the above two criteria, we have a limited number of suitable choices despite the voluminous nonparametric regression methods in the literature. In the following, we discuss a class of regression methods referred to as the sieve estimation method which include the least-squares method of Longstaff and Schwartz (2001) as a special case.
3.3.2 Shape-preserving sieve estimation
We give a brief introduction to the sieve estimation method; refer to Chen (2007) for a comprehensive review. Suppose we have a sample of independent and identically distributed (i.i.d.) random pairs where is a -valued random vector with compact support and is a univariate random variable. Define the function as
[TABLE]
which is independent of . In the context of our BSBU algorithm, and correspond to and , respectively, and the parallel function is the continuation function .
The sieve estimation method strives to estimate the functional form of by solving the following optimization problem:
[TABLE]
where is a finite-dimensional functional space depending on a certain parameter and is called as sieve space. Intuitively, the ampler the sieve space is, the smaller the “gap” between the and the function would be. The price to pay is that larger estimation error is incurred for a richer sieve space due to limited sample size . Therefore, one has to balance such a trade-off by controlling the complexity of the sieve space and this is achieved by tuning the parameter . To make the matter more concrete, we consider two examples of the sieve space in the sequel.
Example 1** (Linear Sieve Space).**
Let be a sequence of basis functions indexed by . Consider the sieve space defined by
[TABLE]
The above set is essentially a linear span of finitely many basis functions and is referred to as linear sieve space in the statistical literature.
In the present context of stochastic control, the regression function corresponds to the continuation function and it exhibits some shape properties such as monotonicity in many applications; see Del Moral et al. (2012) for pricing American option and Huang and Kwok (2016) for valuing equity-linked insurance product, among others. In view of this, it is natural to expect the element in the sieve space satisfies such shape constraints, which in turn preserves the financial interpretations of the numerical result. This can be achieved by considering a special linear sieve space in the following example.
Example 2** (Shape-Preserving Sieve Space).**
Let be a sequence of basis functions indexed by . Denote with . Consider the sieve space defined by
[TABLE]
where is some integer-valued function, is a -by- matrix, and is a -by-1 null vector.
Wang and Ghosh (2012a, b) show that each element in the sieve space in Eq. (29) is a convex, concave, or monotone function (with respect to each coordinate) with a special the choice of matrix given that are Bernstein polynomials. Some forms of are relegated to Appendix A.1.
For a linear sieve space defined either in Eq. (28) or Eq. (29), the solution of the preceding optimization problem (27) is given by the following form:
[TABLE]
where and is the optimizer of the following optimization problem:
[TABLE]
The dependency of and on is suppressed for brevity. In general, one has to solve a constrained quadratic programming problem to obtain .
3.3.3 Discussions
One clear merit of the above linear sieve estimation is that one only needs to store the vector for future evaluation of the regression function at any point in the domain because basis functions are explicitly known at the first hand. This makes the linear sieve estimation method tailored to our present problem in terms of the criterion (C1).
For the criterion (C2), it is well documented in statistical literature that when the true regression function satisfies certain shape constraints, the shape-preserving estimate obtained by (31) with given by Eq. (29) is insensitive to the tuning parameter ; see, e.g., Meyer (2008) and Wang and Ghosh (2012a, b). However, this is legitimate only when the true conditional mean function exhibits such convexity, concavity or monotonicity property. For the general case when there is no prior shape information of , one has to use the sieve space (28) and the regression estimate is sensitive to the choice of . Under such a situation, should be determined in a data-driven manner. In Appendix A.2, we present some common methods of selecting discussed in the statistical literature.
Finally, the convergence of the sieve estimate to the conditional mean function is ensured under some technical conditions. These conditions are summarized in Assumption 4 which is relegated to Appendix A.3 for the clarity of presentation.
3.4 Convergence Analysis of BSBU Algorithm
Now we are ready to conduct convergence analysis of the BSBU algorithm proposed in Section 3.2. For the regression method employed in the algorithm, we restrict our attention to the linear sieve estimator given by Eqs. (30) and (31) in the last subsection.
A complete convergence analysis of the BSBU algorithm should take account of three types of errors:
(E1)
Truncation Error The truncation error is caused by taking as a proxy for .
(E2)
Sieve Estimation Error At each step of the BSBU algorithm, the sieve estimation method is employed to get an estimate for the continuation function. The associated sieve estimation error stems from two resources: (a) the bias caused by using a finite-dimensional sieve space to approximate continuation function; and (b) the statistical error in estimating coefficients of basis functions under a limited sample size of .
(E3)
Accumulation Error The primal goal of the regression step in the BSBU algorithm is to estimate the continuation function of the auxiliary stochastic control problem, i.e., \tilde{C}_{t}(\cdot)={\mathbb{E}}\left[\tilde{V}_{t+1}(X_{t+1})\big{|}X_{t^{+}}=\cdot\ \right]. Thus, in principle, one should generate a random sample
[TABLE]
based on which the sieve estimation method can be employed to get a regression estimate. However, is not known exact at each time step of the algorithm and is replaced by its numerical estimate ; see Step 2.2 of the BSBU algorithm in Section 3.2. Therefore, the the algorithm error accumulated from time step down to triggers a new type of error in addition to (E1) and (E2).
(E1) has been investigated in Theorem 1. The discrepancy between and is contributed by (E2) and (E3). Distinguishing these two types of error plays a crucial role in our convergence analysis and this is inspired by Belomestny et al. (2010). Our main convergence result is summarized in the following theorem.
Theorem 2** (BSBU Algorithm Error).**
Suppose that
(i)
Assumptions 1–3 and Assumption 5 in Appendix B.3 hold;
(ii)
Assumption 4 in Appendix A holds for and uniformly in , where and are given in Steps 2.1 and 2.2 of the BSBU algorithm.
Then, there exists a constant such that
[TABLE]
with “Big O p” notation defined in Definition 2 of Appendix B.3.
The above theorem basically states that the numerical solution converges to in probability as both the number of basis functions and number of simulated paths approach infinity at the rate specified by Condition (v) in Assumption 4. Since Theorem 1 shows that the discrepancy between and shrinks as increases, the numerical estimate is a legitimate approximation for when , , and are considerable. The R.H.S. of Eq. (32) reveals that the overall BSBU algorithm error arises from the two resources discussed in the previous item (E2), which are indicated by the terms and , respectively. Furthermore, Eq. (32) also shows that such a regression error is escalated by a factor at each time step, which reflects the error accumulation from time step down to time step [math] and is in line with the earlier discussion in the item (E3).
4 Application: Pricing Equity-linked Insurance Products
In this section, we apply the BSBU algorithm to the pricing of equity-linked insurance products. This pricing problem is an appropriate example to show the limitations of the FSBU algorithm commented in Section 2.2. For the convenience of illustration, the contract we study here is a simplified version of variable annuities (VAs); for the discussions on more generic policies, we refer to Azimzadeh and Forsyth (2015), Huang and Kwok (2016), Huang et al. (2017), and Shen and Weng (2017), among others.
4.1 Contract Description
We give a brief introduction to the VA. VAs are equity-linked insurance products issued by insurance companies. At the inception of the contract, the policyholder (PH) pays a lump sum to the insurer which is invested into a certain risky asset. The PH is entitled to withdraw any portion of the investment before the maturity. She also enjoys certain guaranteed payments regardless of the performance of the investment account. Therefore, the insurer provides downside protection for a potential market decline. As a compensation, the insurer deducts insurance fees from the investment account and trade available securities to hedge his risk exposure. Thus, no-arbitrage pricing has been the dominating paradigm for pricing VAs in the literature. The primary challenge of this pricing problem stems from the uncertainty of the PH’s withdrawal behavior. This is conventionally resolved by studying the optimal withdrawal strategy of the PH, which naturally leads to a stochastic control problem; see Dai et al. (2008), Chen et al. (2008), Huang and Kwok (2016), and many others.
4.2 Model Setup
In the following, we exemplify the model setup of Section 2 in the present pricing problem. The lattice corresponds to the collection of all available withdrawal dates. The first decision variable represents the PH’s decision to initialize the withdrawal or not by taking values and [math], respectively. As we will see later, the payoff functions depend on the timing of the first withdrawal of the PH. Therefore, a state variable is introduced to record the first-withdrawal-time, and its evolution mechanism is prescribed as follows: , and
[TABLE]
for . The feasible set of is a singleton if the withdrawal has been initialized, i.e., ; otherwise, it is .
Denote and . The second state variable corresponds to the investment account and it evolves according to
[TABLE]
where is the withdrawal amount of the PH at time , is the absolute return of the underlying asset over , and is a certain percentage depending on the first-withdrawal-time . The above equation implies that the PH can withdraw up to the amount of even if the investment account is depleted, i.e., . The jump mechanism of the investment account across each withdrawal date is illustrated in Figure 8.
Now, the state process and the DM’s action are and , respectively, with the superscript “” denoting vector transpose. In accordance with Eqs. (33) and (34), the accompanying transition equation is X_{t+1}=H\big{(}K(X_{t},a_{t}),\varepsilon_{t+1}\big{)}, where
[TABLE]
with . The dependency of on is suppressed for notational brevity.
Next, we discuss the feasible set of the PH’s action. In principle, the withdrawal amount takes values in a continum \big{[}0,W_{t}\vee G(I_{t})P_{0}\big{]} (see Eq. (34)). However, it can be shown that the optimal withdrawal amount is limited to three choices:
- ,
- , and 3) under certain contract specifications; see Azimzadeh and Forsyth (2015), Huang and Kwok (2016), Huang et al. (2017), and Shen and Weng (2017). Via a similar argument adopted by the above references, one may show that this conclusion still holds for the contract considered here. Therefore, we restrict the feasible set of action into the following discrete set:
[TABLE]
for . As a convention, the PH is not allowed to withdraw at inception, and thus .
We proceed by specifying the reward functions which corresponds to the policy payoffs in the present context. Before maturity, the cash inflow of the PH is her withdrawal amount subject to some penalty:
[TABLE]
with being the penalty rate. In other words, the withdrawal amount in excess of the guaranteed amount is subject to a proportional penalty. Conventionally, . At maturity, the policy payoff is the remaining value of the investment account, i.e., .
Finally, we give the interpretation of value function in the present context. with (resp., ) corresponds to the no-arbitrage price of the contract at withdrawal date given that the investment account has a value of and the first withdrawal is triggered at withdrawal date (resp., no withdrawal has been taken).
4.3 A BSBU Algorithm for the Pricing Problem
The state process generally takes value in the unbounded set . We consider a truncated domain: with . Consequently, we may define the auxiliary state process as in Eq. (13). The range of the post-action value is given by , where and , respectively. This is in line with Eq. (21).
Now we are almost ready to employ the BSBU algorithm developed in Section 3.2 to solve the present pricing problem. It is worth noting that a discrete state variable appears in the present context and the continuation function, in general, is not continuous with respect to the post-action value accompanying this state variable, i.e., ; see Eq. (35). Consequently, Condition (iii) of Assumption 4 might not hold here; see Appendix A.3. However, for each given value of , the continuation function is still continuous with respect to , the post-action value associated with the investment account value. And therefore, one may repeat Step 3.2 of the BSBU algorithm for every distinct value of . It is easy to see the convergence of the resulting BSBU algorithm is not influenced by this modification.
Finally, it remains to specify how to simulate the post-action value of the state process in order to pave the way to implementing the BSBU algorithm. In the sequel section, we will address this issue in details and, in particular, we will compare the control randomization method with our artificial simulation method.
5 Numerical Experiments
This section devotes to conducting numerical experiments to show the merits of the BSBU algorithm in the context of pricing the variable annuity product addressed in the last section.
5.1 Parameter Setting
We first present the parameter setting for our numerical experiments. We consider time steps and the time interval between two consecutive withdrawal date is assumed to be . This corresponds to a contract with one-year maturity and monthly withdrawal frequency. The PH’s initial investment is assumed to be one unit, i.e., . The guaranteed payment percentage is prescribed as follows:
[TABLE]
In other words, the PH enjoys a larger amount of guaranteed payment if she postpones the initiation of the withdrawal. As a result, the value function/continuation is not continuous with respect to the state variable .
Let and be the annualized risk-free rate and insurance fee rate, respectively. We assume the absolutely return of the underlying fund follows a log-normal distribution with {\mathbb{E}}[\log\varepsilon_{t+1}]=\big{(}r-q-\sigma^{2}/2\big{)}\delta and under a risk-neutral pricing measure. This implicitly assumes the underlying fund evolves according to a Geometric Brownian Motion with annualized volatility rate . Finally, the discounting rate is given by . All of the above market and contract parameters are summarized in Table 1.
Finally, we discuss the choice of truncation parameter . Under the present context, it is easy to see that the function in Assumption 2 is bounded from above by the tail probability of the continuous running maximum of a geometric Brownian Motion. To be specific, we have
[TABLE]
with , where is a standard Brownian Motion, is the cumulative distribution function of a standard normal distribution, and the last equality follows by the Reflection Principle (see, e.g., Shreve (2004, pp. 297)). Let . Then the R.H.S. of the above inequality approximately equals to under the parameter setting in Table 1. It is also easy to see that is quadratic in in the present example, and therefore, the truncation error is marginal according to the error bound in Eq. (20). In view of this, we fix in all subsequent numerical experiments.
5.2 Forward Simulation v.s. Artificial Simulation
Next, we would like to show the limitations of the forward simulation based on control randomization in generating random samples of the state process. Below, we present some prevalent control randomization methods.
(CR0)
Given a simulated , the PH’s action is simulated from a degenerated distribution with one single point mass at .
(CR1)
Given , the DM’s action is simulated from a discrete uniform distribution with support set \big{\{}(0,0)^{{\intercal}},\left(G(I_{t})P_{0},1\right)^{{\intercal}},\left(W_{t},1\right)^{{\intercal}}\big{\}} if ; and \big{\{}(0,0)^{{\intercal}},\left(G(I_{t})P_{0},1\right)^{{\intercal}},\left(W_{t},1\right)^{{\intercal}}\big{\}}, otherwise.
(CR2)
Given , the DM’s action is simulated from a discrete uniform distribution with support set \big{\{}(0,0)^{{\intercal}},\left(G(I_{t})P_{0},1\right)^{{\intercal}}\big{\}} if ; and \big{\{}(0,0)^{{\intercal}},\left(G(I_{t})P_{0},1\right)^{{\intercal}}\big{\}}, otherwise.
Given the above rules of generating the PH’s action, one may simulate the state process in a forward manner in accordance with Steps 2.1 and 2.2 of the FSBU algorithm; see Section 2.2.
(CR0) is first proposed by Huang and Kwok (2016) in the context of pricing Guaranteed Lifelong Withdrawal Benefit, a particular type of variable annuity policy. It initializes the withdrawal at and the resulting simulated the state variable (resp., its accompanying post-action value ) equals a fixed value for all although (resp., ), in principle, can take any value in (resp., ). A consequential annoying issue is that the obtained estimate for the value function/continuation function is invariant to the first-withdrawal-time , which is not sensible since the later the PH initializes the withdrawal the larger guaranteed amount she could enjoy in remaining contract life.
(CR1) uniformly simulates the PH’s action from its feasible set ; see Eq. (36). By virtue of this, there always exists some paths with which corresponds to the scenario that the withdrawal has not been initialized. This in turn guarantees that, in principle, (resp., ) can take any value in (resp., ). However, this strategy is also not satisfactory: an overwhelming portion of paths are absorbed by the state , i.e., the depletion of investment account, and very sparse sample points of the investment account are positive. This is graphically illustrated in the top panel of Figure 9 where sample paths are plotted for the clarity of presentation. So it is not hard to expect that the accuracy of the regression estimate is severely impaired over .
To alleviate the serious problem mentioned above, (CR2) discards the strategy of depleting the investment account, i.e., , in simulating the PH’s action. And therefore, the simulated investment account value can spread over a wider range than that accompanying (CR1); see the bottom panel of Figure 9. This phenomenon is more palpable from the histograms of which are collected by Figure 10. Nevertheless, (CR2)’s performance in simulating the is undesirable: Figure 11 shows that a substantial portion of sample points of the first-withdrawal-time are concentrated in first few values that can take. To understand the crux, we note that at the first possible withdrawal date, one-half of sample paths exhibit the initiation of the withdrawal; among the remaining paths, one-half of them witness the withdrawal in the consecutive withdrawal date. Therefore, the portion of positive declines at an exponential rate as increases which is in line with Figure 11. In view of this, it can be expected that the consequential numerical estimate for the value function sustains significant error at state with a large .
Overall, none of the above rules (CR0)–(CR2) gives agreeable performance. It is hard to figure out an ideal way to randomize the PH’s action which can sidestep the thorny issues mentioned above. This shows one drawback of binding together control randomization and forward simulation in addition to the issue of computational cost; see also the item “Limitation of control randomization” of Section 2.2.
To circumvent the annoying problems mentioned above, in the sequel numerical experiments, we simulate the post-action value of the state process at each time step as follows: where and are simulated from two independent uniform distributions with support sets and , respectively. This ensures the post-action value evenly distributed over .
5.3 Raw Sieve Estimation v.s. Shape-Preserving Sieve Estimation
In the sequel, we conduct several numerical experiments to compare the regression estimates for the continuation function produced by two regression methods: the raw sieve estimation (RSE) method and the shape-preserving sieve estimation (SPSE) method. The RSE and SPSE are essentially the linear sieve estimation method discussed in Section 3.3 with sieve spaces (28) and (29), respectively. It is easy to show that the continuation function is monotone and therefore we incorporate this shape constraint in the SPSE method. The expression of the accompanying matrix is given in Appendix A.1. The RSE method is equivalent to the least-squares method commonly adopted in the literature; see, for instance, Longstaff and Schwartz (2001). The aim of this subsection is to show the advantage of incorporating shape constraints in the regression step of an LSMC algorithm.
In the first numerical experiment, we compare the SPSE and RSE for the continuation function at each time step. For the fairness of the comparison, for both methods, we take \bm{\phi}(\cdot)=\big{(}\phi_{0}(\cdot),\dots,\phi_{J}(\cdot)\big{)}^{{\intercal}} as a vector of univariate Bernstein polynomials up to order in both sieve estimation methods. Figure 12 collects the plots of regression estimates as a function of at odd time steps with . To better show the subtle difference between the estimates produced by SPSE and RSE, the plots are restricted on the interval . From Figure 12, we can see that the discrepancy between the regression estimates accompanying SPSE and RSE is not conspicuous at large time step but becomes more significant as the time step goes down. Despite the continuation function, in principle, is a monotone function with respect to , the dotted lines in Figure 12 show that its regression estimate produced by the RSE does not inherit this monotonicity and loses certain economic interpretations, accordingly. This issue is more serious at smaller time steps as shown by the bottom panel of Figure 12. This is not surprising because once the monotonicity is lost at a certain time step, the regression estimate obtained in the consecutive time step will be influenced, which in turn exaggerates the problem. In contrast, as depicted by the solid lines in Figure 12, the SPSE method always preserves the monotonicity of the continuation function and therefore the corresponding regression estimates are more economically sensible. This shows the first advantage of the SPSE method in terms of preserving certain shape properties of the continuation function.
Next, we implement the BSBU algorithm to compute with where the SPSE and RSE are employed, respectively. These estimates approximate the no-arbitrage price of the VA policy at the inception and therefore are of most interest in the present context; see the last paragraph of Section 4.2. It is worth noting that is random due to the randomness of the simulated sample; see also Remark 1. And therefore, we repeat the BSBU algorithm 40 times in order to study the stability of under a finite sample size. Table 2 summarizes the mean and standard deviation of under different pairs of and . The “S.d.” column of the table discloses that the standard deviation accompanying the SPSE is nearly one half of that associated with the RSE under all numerical settings. For the numerical settings 0-2 of Table 2, Figure 13 delineates the corresponding density plots of the 40 estimates. By comparing the left and right panels of Figure 13, we can easily perceive that the numerical estimates accompanying SPSE are less volatile as reflected by the more spiked shape of the corresponding density plots. This observation is consistent with Table 2. To sum up, the SPSE surpasses the RSE in terms of smaller standard deviation of the resulting numerical estimate for the optimal value function at the initial state.
From the settings 0-2 of Table 2, we also observe that for both methods, the change of the mean of the numerical estimate is not substantial as the number of basis functions hikes from to . In the settings 1, 3 and 4 of Table 2, we fix and increase the number of simulated paths from to . We witness that standard deviation decreases as the number of simulated paths climbs. This descending trend is also confirmed by the box plots depicted in Figure 14: the height of the box shrinks as the number of simulated paths hikes. All of these show the convergence of the BSBU algorithm which is in line with the convergence result established in Section 3.4.
Overall, the advantages of the SPSE over the RSE are extant at least in two-fold. Firstly, the SPSE produces economically sensible regression estimates by inheriting certain shape properties of the true continuation function. Secondly, the consequential estimate for the optimal value function accompanying the SPSE method is less volatile than that produced by the RSE method under a finite number of simulated sample paths.
6 Conclusion
This paper develops a novel LSMC algorithm, referred to as Backward Simulation and Backward Updating (BSBU) algorithm, to solve discrete-time stochastic optimal control problems. We first introduce an auxiliary stochastic control problem where the state process only takes value in a compact set. This enables the BSBU algorithm to successfully sidestep extrapolating value function estimate. We further show the optimal value function of the auxiliary problem is a legitimate approximation for that of the original problem with an appropriate choice of the truncation parameter. To circumvent the drawbacks of forward simulation and control randomization, we propose to directly simulate the post-action value of the state process from an artificial probability distribution. The pivotal idea behind this artificial simulation method is that the continuation function is solely determined by the distribution of random innovation term. Moreover, motivated by the shape information of the continuation function, we introduce a shape-preserving sieve estimation technique to alleviate the computational burden of tuning parameter selection involved in the regression step of an LSMC algorithm. Furthermore, convergence result of the BSBU algorithm is established by resorting to the theory of nonparametric sieve estimation. Finally, we confirm the merits of the BSBU algorithm through an application to pricing equity-linked insurance products and the corresponding numerical experiments.
Appendix A Supplements for Sieve Estimation Method
A.1 Forms of Matrix
To make the paper self-contained, we collect several forms of the constraint matrix in (29) which ensures monotonicity, convexity, or concavity of the sieve estimate (30); see Wang and Ghosh (2012a, b) for a justification. For the simplicity of notation, we only address the case .
Monotonicity Suppose the conditional mean defined in Eq. (26) is monotone. Then the corresponding monotonicity-preserved sieve estimate is obtained from (30) with given by Eq. (29) and
[TABLE]
Convexity/Concavity If is convex, we choose the matrix as
[TABLE]
Moreover, the matrix accompanying a concave is obtained by taking negative of the above matrix.
Convexity and Monotonicity If is convex and monotone, the corresponding is given by
[TABLE]
A.2 A Data-driven Choice of
Below we present some common methods of choosing the number of basis functions in a sieve estimation method; see, e.g., Li (1987).
Mallows’s For a discrete set , is determined by solving the following minimization problem:
[TABLE]
where is given in (30) and which is an estimate for the variance of residual term.
Generalized cross-validation is determined by
[TABLE]
with given in (30).
Leave-one-out cross-validation Select to minimize
[TABLE]
where is similarly obtained by Eq. (30) with the sample point removed.
It is worth stressing that, among the above three selection methods, the leave-one-out cross-validation method is most computationally expensive as one has to compute the regression estimate times in a single evaluation of . This is clearly computationally prohibitive when is considerable, which is particularly the case in the present context of the LSMC algorithm. The Mallows’s criterion and generalized cross-validation method are relatively less cumbersome but still undesirable under a sizable . These show the importance of avoiding such a tuning parameter selection procedure which is one important thrust behind proposing shape-preserved sieve estimation method.
A.3 Technical Assumption of Sieve Estimation Method
We impose the following assumption accompanying the sieve estimation method discussed in Section 3.3 which follows from Newey (1997).
Assumption 4**.**
(i)
* are i.i.d. and has compact support . Furthermore, is bounded over .*
(ii)
There exists a sequence such that with denoting the supremum norm of a continuous function over .
(iii)
For the sieve space defined either in Eq. (28) or Eq. (29), there exists a -by-1 vector and a sequence such that as , and
[TABLE]
where we remind that .
(iv)
Let . There exists a positive constant independent of such that with and denoting the smallest and largest eigenvalues of , respectively.
(v)
As , , and .
We give some comments on the above technical conditions.
The i.i.d. condition in Part (i) of the above assumption discloses the necessity of generating an independent sample at each time step in an LSMC algorithm; see also the discussion in the earlier item “Cost of forward simulation” of Section 2.2. Part (i) further requires has a compact support, which is conventional in literature see, e.g., Newey (1997) and Chen (2007). In the context of BSBU algorithm, this shows that restraining the state process into a bounded domain is not only beneficial in eliminating undesirable extrapolation but also indispensable in guaranteeing the convergence of the regression estimate to the continuation function. This has also been pointed out in the literature, see, e.g., Stentoft (2004) and Zanger (2013). 2. 2.
Part (ii) specifies how the magnitude of is amplified as the number of basis function functions grows up. In particular, Newey (1997) shows that for B-splines and for power series; for the cases of other types of basis functions, we refer to Chen (2007). 3. 3.
Part (iii) states that there exists a function in the sieve space that “best” approximates the conditional mean function under the supremum norm; see Figure 15 for a graphical illustration. The existence of vector (referred to as oracle) is guaranteed by the convexity of sieve space . For the sieve space (29), the existence of relies on the convexity, concavity or monotonicity of the function which follows by the Property 3.2 of Wang and Ghosh (2012b). Figure 15 depicts the relationship between and : their discrepancy vanishes as increases and, for a fixed , the sieve estimate converges to as the sample size approaches infinity. Therefore, one may view the sieve estimation as a two-stage approximation for the conditional mean function . 4. 4.
The condition in Part (iv) ensures the design matrix of the regression problem is nonsingular with a high probability and does not blow up as approaches infinity. Finally, Part (v) prescribes the growth rates of and in order to avoid overfitting or underfitting.
Appendix B Proofs of Statements
B.1 Proof of Proposition 1
B.1.1 Preliminary
Lemma 1**.**
For any -adapted process , the following statements hold:
(i)
* for ;*
(ii)
* for ,*
where and are defined in Eqs. (11) and (12), respectively.
Proof of Lemma 1.
(i)
According to Eq. (12), we observe
[TABLE]
where the second identity is by the definition of the stopping time and the fact that . To show the statement in Part (i) of Lemma 1, it suffices to prove . Indeed, implies , and thus .
(ii)
In view of Part (i) and Eq. (12), we obtain
[TABLE]
Therefore, we obtain
[TABLE]
This proves Part (ii) of Lemma 1.
∎
B.1.2 Proof of the Main Result
Proof of Proposition 1.
By exploiting Lemma 1 and Eq. (12), we get
[TABLE]
where the fourth equality follows from Eq. (12) and the last equality follows from Lemma 1.
The above equation together with Eqs. (5) and (14) yields Eq. (13). This completes the proof. ∎
B.2 Proof of Theorem 1
B.2.1 Preliminary
Recall that and implicitly depend on certain actions ; see Eqs. (1) and (13), respectively. In the sequel, we sometimes stress such dependency by writing (resp. ) and (resp. ).
Lemma 2**.**
For the state process defined through Eq. (13), the following statements hold.
(i)
For any , there exists such that for all almost surely.
(ii)
For any , there exists such that for all almost surely.
Proof of Lemma 2.
(i)
Given and , we construct as follows: , and
[TABLE]
where for and .
It is easy to see from the above construction that is -adapted. It remains to show that
[TABLE]
Firstly, we observe and . As induction hypothesis, we assume the statement (B.2) holds for time step . For time step , we split the discussions into two cases.
- 1.
If , then
[TABLE]
where the first and third equalities follow by Eq. (13) and the second equality is due to the induction hypothesis. 2. 2.
In the second case that , we apply Eq. (13) to get
[TABLE]
where the second equality follows by Eq. (B.1) and the induction hypothesis (B.2).
In either of the above cases, we have . This combined with Eq. (B.1) implies
[TABLE]
Otherwise, . This proves the statement (B.2) for time step . The proof of Part (i) is complete.
(ii)
Given and , we construct as follows: , and
[TABLE]
where is any measurable function satisfying for and .
It is easy to see that is -adapted. Next, we use a forward induction argument to show that
[TABLE]
The above statement holds trivially for . As induction hypothesis, we assume it holds for time step . For time step , we consider two separate cases.
- 1.
If , Eq. (13) in combined with (B.4) implies
[TABLE] 2. 2.
In the second case that , applying Eq. (13) gives
[TABLE]
where the second equality follows by Eq. (B.3) and the induction hypothesis (B.4).
Overall, we always observe . To prove the statement (B.4) holds for time step , it remains to show We split the discussion into two separate cases.
- 1.
Firstly, suppose Eq. (B.3) implies
[TABLE] 2. 2.
Secondly, suppose By Part (i) of Lemma 1, and thus, it follows from Eq. (12) that . Consequently, we apply Eq. (B.3) to get
[TABLE]
The proof of Part (ii) is complete.
∎
A direct consequence of the preceding lemma is the following corollary.
Corollary 1**.**
The value function defined in Eq. (15) exhibits:
[TABLE]
It is worth noting that the optimization problems in Eq. (B.5) and Eq. (15) are taken over the set and , respectively. The above corollary states that the optimal values of these two optimization problems are exactly the same as given by .
B.2.2 Proof of the Main Result
Proof of Theorem 1.
In view of Eqs. (2) and (B.5), we obtain
[TABLE]
Below we establish upper bounds for the and defined in the above display, respectively. Let . Note that
[TABLE]
Accordingly, we get
[TABLE]
with , where the first inequality is by triangular inequality and Assumption 3 and the third inequality is due to Part (ii) of Assumption 3. Applying Cauchy–Schwarz inequality twice gives
[TABLE]
where the third inequality follows because for two real numbers and . In view of Assumption 3, we get
[TABLE]
Combing the last two displays with Assumption 2 implies
[TABLE]
A similar argument gives
[TABLE]
Combining (B.6), (B.8), and (B.9) together implies
[TABLE]
The proof is complete. ∎
B.3 Proof of Theorem 2
B.3.1 Preliminary lemmas
We first give the definitions of “Big O p” and “Small O p” notations which are commonplaces in statistical literature.
Definition 2**.**
(i)
For two sequences of random variables and indexed by , we say if
(ii)
Moreover, we say if for all .
Some Matrices Let \mathbf{h}_{t}(x)=\left(\sup\limits_{a\in A_{t}(x)}\phi_{1}\big{(}K(x,a)\big{)},\dots,\sup\limits_{a\in A_{t}(x)}\phi_{J}\big{(}K(x,a)\big{)}\right)^{{\intercal}}, for , and we suppress its dependency on . Define matrices
[TABLE]
for with the superscript denoting vector transpose. It is palpable that is a finite-sample estimate for . In the sequel, we denote (resp. ) as the largest (resp. smallest) eigenvalue of a square matrix . We impose the following Assumption on the eigenvalues of .
Assumption 5**.**
(i)
For any fixed and , is a compact set. Moreover, and are continuous functions for .
(ii)
There exists a positive constant independent of and such that .
Part (i) of the preceding assumption guarantees that the function is well-defined for . The continuity requirement of can be removed if is a lattice (discrete set), which is particularly the case when the stochastic optimal control problem exhibits the Bang-bang solution, see, e.g., Azimzadeh and Forsyth (2015) and Huang and Kwok (2016). Part (ii) requires the largest eigenvalue of the matrix does not blow up as and approach infinity. This condition ensures the sample eigenvalue converges to the non-sample counterpart as approaches infinity as shown in the sequel Lemma 3.
Moreover, we define matrices
[TABLE]
The following lemma relates the eigenvalues of and to those of and .
Lemma 3**.**
(i)
Suppose Condition (ii) of Theorem 2 is satisfied. Then,
[TABLE]
and
[TABLE]
for .
(ii)
Suppose Assumption 5 holds. In addition, Condition (v) of Assumption 4 is satisfied. Then, for .
Proof of Lemma 3.
Lemma 3 can be proved by a similar argument as that used in the proof of Eq. (A.1) in Newey (1997). ∎
The above lemma shows the sample eigenvalues converge to the non-sample counterparts as approaches infinity. In view of Condition (iv) of Assumption 4, Lemma 3 also implies the largest (resp., smallest) eigenvalue of is bounded from above (resp., below) with probability approaching 1 as . This fact is exploited in the proofs of sequel Lemmas 4 and 5.
Pseudo Estimate, Oracle, and True Estimate Next, we introduce the concept of pseudo estimate. Let (resp. ) be the solution to the optimization problem in Eq. (31) with (resp. ) and . Given and , denote the associated regression estimates by and , respectively. is essentially the sieve estimate for the continuation function when the true value function is employed in the regression. We further define function for by substituting in Eq. (25) with . For , we set with given by Eq. (16).
Admittedly, in the implementation of the BSBU algorithm, is not tractable because the true value function is unknown and should be replaced by the numerical estimate obtained inductively. For this reason, following Belomestny et al. (2010), we call the pseudo estimate. Despite this, the pseudo estimate plays an indispensable role in establishing the convergence result of Theorem 2. In addition to the two estimates and defined in the above, we further define the oracle as the solution to the optimization problem (A.1) with replaced by .
The following lemma discloses that the gap between pseudo estimate and the oracle vanishes when both and increase at a certain rate.
Lemma 4**.**
Suppose the conditions of Theorem 2 are satisfied. Then,
[TABLE]
Proof of Lemma 4.
Recall that solves the optimization problem:
[TABLE]
On the other hand, is a suboptimal solution to the above optimization problem. Therefore, we get
[TABLE]
where is a -by- matrix with -th row being and is a -by-1 vector with -th element given by .
By adding and subtracting the term in the L.H.S. of the above inequality, we get
[TABLE]
where we use the shorthand notations and . Expanding both sides of the above inequality gives
[TABLE]
where the second inequality is by Hölder’s inequality. For the L.H.S. of the above inequality, it follows from the definition of the smallest eigenvalue that
[TABLE]
Combing the last two inequalities together implies
[TABLE]
It follows from Lemma 3 that the event holds with probability approaching as . And therefore,
[TABLE]
holds with probability approaching as .
It follows as in Eq. (A.2) of Newey (1997, pp. 163) that This in conjunction with the last display proves the desired result. The proof is complete. ∎
The next lemma relates the discrepancy between the pseudo estimate and the true estimate to the estimation error of the value function at the previous time step.
Lemma 5**.**
Suppose the conditions of Theorem 2 are satisfied. Then, for , there exists a constant independent of , and such that
[TABLE]
holds with probability approaching 1 as , where and are two -by-1 vectors with -th element given by and , respectively.
Proof of Lemma 5.
Using the argument as in the proof of inequality (B.10), we obtain
[TABLE]
holds with probability approaching 1 as , where we adopt shorthand notations and On the other hand, it follows from Lemma 3 that
[TABLE]
holds with probability approaching 1 as .
Combing the above two inequalities implies
[TABLE]
By adding and subtracting the term in the R.H.S. of the above inequality, we get
[TABLE]
where the last equality is guaranteed by Part (ii) of Assumption 4. Combing the last two inequalities implies
[TABLE]
holds with probability approaching 1 as . This proves Lemma 5. ∎
The statement of the above lemma is not hard to expect because the primary difference between the pseudo estimate and the true estimate stems from the the estimation error of value function.
The final lemma quantifies the discrepancy between the value function and its numerical estimate under the empirical norm.
Lemma 6**.**
Let be the probability distribution function of for . Suppose the assumptions of Theorem 2 hold. Then
[TABLE]
where and are two -by- vectors with -th element being and , respectively.
Proof of Lemma 6.
We use a backward induction procedure to prove the statement of Lemma 6. For , we note that is in agreement with because for . We get for , accordingly. Furthermore, we observe that
[TABLE]
where the third inequality is by the definition of function and the last equality is guaranteed by Assumption 4.
Consequently, we obtain
[TABLE]
where the second inequality follows from the definition of the largest eigenvalue of a matrix and the last equality is guaranteed by Lemma 4 and Lemma 3. In view of the above display, Eq. (B.11) holds for .
As induction hypothesis, we assume (B.11) holds for . Note that, for ,
[TABLE]
An argument similar to the one used in establishing (B.13) shows that
[TABLE]
Next, we investigate the term . Observe that
[TABLE]
We adopt the same argument as in the proof of (B.13) to get
[TABLE]
Applying Lemma 5 yields
[TABLE]
where the last equality is due to induction hypothesis (B.11) and (see Lemma 3). The above display in conjunction with (B.14) and (B.15) implies
[TABLE]
This completes the proof. ∎
B.3.2 Proof of the Main Result
Proof of Theorem 2.
Following the arguments used to prove (B.12), we get
[TABLE]
where the last equality is by Lemma 4 and Part (ii) of Assumption 5.
On the other hand, an argument similar to the one used in deriving (B.16) shows
[TABLE]
The above two displays in conjunction with (B.11) implies
[TABLE]
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Azimzadeh and Forsyth (2015) Parsiad Azimzadeh and Peter A Forsyth. The existence of optimal bang-bang controls for gmxb contracts. SIAM Journal on Financial Mathematics , 6(1):117–139, 2015.
- 2Barrera-Esteve et al. (2006) Christophe Barrera-Esteve, Florent Bergeret, Charles Dossal, Emmanuel Gobet, Asma Meziou, Rémi Munos, and Damien Reboul-Salze. Numerical methods for the pricing of swing options: a stochastic control approach. Methodology and computing in applied probability , 8(4):517–540, 2006.
- 3Belomestny (2011) Denis Belomestny. Pricing bermudan options by nonparametric regression: optimal rates of convergence for lower estimates. Finance and Stochastics , 15(4):655–683, 2011.
- 4Belomestny et al. (2009) Denis Belomestny, Grigori Milstein, and Vladimir Spokoiny. Regression methods in pricing american and bermudan options using consumption processes. Quantitative Finance , 9(3):315–327, 2009.
- 5Belomestny et al. (2010) Denis Belomestny, Anastasia Kolodko, and John Schoenmakers. Regression methods for stochastic control problems and their convergence analysis. SIAM Journal on Control and Optimization , 48(5):3562–3588, 2010.
- 6Carmona and Ludkovski (2010) René Carmona and Michael Ludkovski. Valuation of energy storage: An optimal switching approach. Quantitative finance , 10(4):359–374, 2010.
- 7Carriere (1996) Jacques F Carriere. Valuation of the early-exercise price for options using simulations and nonparametric regression. Insurance: mathematics and Economics , 19(1):19–30, 1996.
- 8Chen (2007) Xiaohong Chen. Large sample sieve estimation of semi-nonparametric models. Handbook of econometrics , 6:5549–5632, 2007.