Stochastic control problems in finance often involve complex controls at discrete times. As a result, numerically solving such problems using, for example, methods based on partial differential or integrodifferential equations inevitably gives rise to low-order (usually at most second-order) accuracy. In many cases, Fourier methods can be used to efficiently advance solutions between control monitoring dates, and numerical optimization methods can then be applied across decision times. However, Fourier methods are not monotone, and as a result they give rise to possible violations of arbitrage inequalities. This is problematic in the context of control problems, where the control is determined by comparing value functions. In this paper, we give a preprocessing step for Fourier methods that involves projecting the Green’s function onto the set of linear basis functions. The resulting algorithm is guaranteed to be monotone (to within a tolerance), ℓ ∞-stable and satisfies an ε-discrete comparison principle. In addition, the algorithm has the same complexity per step as a standard Fourier method and second-order accuracy for smooth problems.
Optimal stochastic control problems in finance often involve monitoring or making decisions at discrete points in time. These monitoring times typically cause difficulties when solving optimal stochastic control problems numerically, both for efficiency and correctness. This is because, for efficiency, numerical methods are typically applied from one monitoring time to the next; correctness arises as an issue when the decision is determined by comparing value functions, which is somewhat problematic when discrete approximations are not monotone. These optimal stochastic problems arise in many important financial applications and include problems such as asset allocation (Li and Ng 2000; Huang 2010; Forsyth and Vetzal 2017; Cong and Oosterlee 2016), pricing of variable annuities (Bauer et al 2008; Dai et al 2008; Chen et al 2008; Ignatieva et al 2018; Alonso-Garcia et al 2018; Huang et al 2017) and hedging in discrete time (Remillard and Rubenthaler 2013; Angelini and Herzel 2014).
These optimal control problems are typically modeled as the solutions of partial integrodifferential equations (PIDEs), which can be solved via numerical, finite-difference (Chen et al 2008) or Monte Carlo (Cong and Oosterlee 2016) methods. When cast into dynamic programming form, the optimal control problem reduces to solving a PIDE backward in time between each decision point, and then determining the optimal control at each such point. In many cases, including those mentioned above, the models are based on fairly simple stochastic processes, with the main interest being the behavior of the optimal controls. These parsimonious stochastic models can be justified if we are looking at long-term problems, eg, variable annuities or saving for retirement, where the time scales are of the order of ten to thirty years. In these situations, it is reasonable to use a parsimonious stochastic process model.
In these (and many other) situations, the characteristic function of the associated stochastic process is known in closed form. For the PIDE types that appear in financial problems, knowing the characteristic function implies that the Fourier transform of the solution is also known in closed form. By discretizing these Fourier transforms, we obtain an approximation to the solution, which can be used for effective numerical computation. A natural approach in this case is to use a Fourier scheme to advance the solution in a single time step between decision times, and then to apply a numerical optimization approach to advance the solution across the decision time. This technique is repeated until the current time is reached (Ruijter et al 2013; Lippa 2013). These methods are based on Fourier space time-stepping (FST) (Jackson et al 2008), the convolution (CONV) technique (Lord et al 2008) or the Fourier cosine (COS) algorithm (Fang and Oosterlee 2008). Fourier methods have been applied to, among other things, the pricing of exotic variance products and volatility derivatives (Zheng and Kwok 2014), guaranteed minimum withdrawal benefits (Ignatieva et al 2018; Alonso-Garcia et al 2018; Huang et al 2017) and equity-indexed annuities (Deng et al 2017).
Fourier methods have a number of advantages over finite-difference and other methods; primarily, there are no time-stepping errors between decision dates. These methods also provide easy handling for stochastic processes involving jump diffusion (Lippa 2013) and regime switching (Jackson et al 2008). Although Fourier methods typically need a large number of discretization points, the algorithms reduce to using finite fast Fourier transforms (FFTs) that operate efficiently on most platforms (including graphics processing units). The algorithms are also quite easy to implement. For example, using Fourier methods for the pricing of variable annuities reduces to the use of discrete FFTs and local optimization. A detailed knowledge of partial differential equation algorithms is not actually required in this case. Fourier methods also easily extend to multifactor stochastic processes, where finite-difference methods have difficulties because of cross-derivative terms. Of course, Fourier methods suffer from the curse of dimensionality, and hence are restricted, except in special cases, to problems of dimension 3 or less. Finally, Fourier methods have good convergence properties for problems with noncomplex controls. For example, for European option pricing, in cases where the characteristic function of the underlying stochastic process is known, the COS method achieves exponential convergence (with regard to the number of terms in the Fourier series) (Fang and Oosterlee 2008).
A major drawback with current Fourier methods is that they are not monotone. In the contingent claims context, monotone methods preserve arbitrage inequalities, or discrete comparison properties, independent of any discretization errors. As a concrete example, consider the case of a variable annuity contract, with ratchet features and withdrawal controls at each decision date. Suppose contract has a larger payoff at the terminal time than contract . Then a monotone scheme generates a value for contract that is always larger than the value of contract , at all points in time and space, regardless of the accuracy of the numerical scheme. In a sense, the arbitrage inequality (discrete comparison) condition is the financial equivalence of conservation of mass in engineering computations. Use of nonmonotone methods is especially problematic in the context of control problems, where the control is determined by comparing value functions.
Monotonicity is also relevant for the convergence of numerical schemes. In general, optimal control problems posed as PIDEs are nonlinear and do not have unique solutions. The financially relevant solution is the viscosity solution of the PIDEs, and it is well known that a discretization of a PIDE converges to the viscosity solution if it is monotone, consistent and stable (Barles and Souganidis 1991). There are examples where nonmonotone discretizations fail to converge (Obermann 2006) and also examples where there is convergence (Pooley et al 2003) but not to the financially correct viscosity solution. In addition, in cases where the Green’s function has a thin peak, existing nonmonotonic Fourier methods require a very small space step that often results in numerical issues. Finally, monotone schemes are more reliable for the numerical computation of Greeks (ie, derivatives of the solution), which is often important for financial instruments.
The starting point for this paper is the assumption that we have a closed-form representation of the Fourier transform of the Green’s function of the stochastic process PIDE. From a practical point of view, we also assume that a spatial shift property holds. The latter assumption can be removed but at the cost of increasing the computational complexity of our method. We will discuss these assumptions further below.
In this paper, we present a new Fourier algorithm in which monotonicity can be guaranteed to within a user-specified numerical tolerance. The algorithm is for use with general optimal control problems in finance. In these general control problems, the objective function may be complex and nonsmooth; hence the optimal control at each step must be determined by a numerical optimization procedure. Indeed, in many cases, this is done by discretizing the control and using exhaustive search. Reconstructing the Fourier coefficients is typically done by assuming the control is constant over discretized intervals of the physical space, by numerically determining the control at the midpoint of these intervals and, finally, by reconstructing the Fourier coefficients by quadrature. This is equivalent to using a type of trapezoidal rule to reconstruct the Fourier coefficients, and hence this can be at most second-order accurate (in terms of the physical domain mesh size).
In fact, we show how the FST or CONV schemes can be modified to get new schemes in which monotonicity can be guaranteed to within a user-specified numerical tolerance. Our approach is similar to that used in these schemes, which first approximate the solution of a linear PIDE by a Green’s function convolution, then discretize the convolution and finally carry out the dense matrix–vector multiplication efficiently using an FFT. In our case, we discretize the value function and generate a continuous approximation of the function by assuming linear basis (or alternatively piecewise constant) functions. Given this approximation, we carry out an exact integration of the convolution integral and then truncate the series approximation of this integral so that monotonicity holds to within a certain tolerance. Consequently, we prove that our algorithm has an -discrete comparison property, that is, given a tolerance , a discrete comparison (also called arbitrage inequality) holds to , independent of the number of discretization nodes and time steps. This is similar in spirit to the -monotone schemes discussed in, for example, Bokanowski et al (2018). Typically, the convergence to the integral is exponential in the series truncation parameter, so it is inexpensive to specify a small tolerance. The key idea here is that the number of terms required to accurately determine the projection of the Green’s function onto linear basis functions can be larger than the number of basis functions. After an initial setup cost, the complexity per step is the same as for the standard FST or CONV methods. This requires only a small change to existing codes in order to guarantee monotonicity. The desirable property of our method is that monotonicity can be guaranteed (to within a certain tolerance) independent of the number of FST (CONV) grid nodes or the time-step size.
While Fourier methods have good convergence properties for vanilla contracts or problems where controls are smooth, it is a different story for general optimal control problems. For example, if the COS method is applied to optimal control problems, then it is challenging to maintain exponential convergence, as the optimal control must be determined in the physical space. Hence, a highly accurate recursive expression for the Fourier coefficients must be found after application of the optimal controls in order to maintain exponential convergence. In the case of bang-bang controls, it is often possible to separate the physical domain into regions where the control is constant. If these regions are determined to high accuracy, then an accurate algorithm for recursive generation of the Fourier coefficients can be developed (Ruijter et al 2013). However, even for the case of an American option, this requires careful analysis and implementation (Fang and Oosterlee 2008). Our interest is in general problems, where the control may not be of the bang-bang type, and we expect that such good convergence properties will not hold. In addition, in the path-dependent case, the problem is usually converted to Markovian form through additional state variables. The dynamics of these state variables are typically represented by a deterministic equation (between monitoring dates). At monitoring dates, the state variable may have nonsmooth jumps (eg, cashflows), and hence the standard approach would be to discretize this state variable and then to interpolate the value function across the monitoring dates. If linear interpolation is used, this also implies that the solution is at most second-order accurate at a monitoring date.
While monotone schemes have good numerical properties, they appear to be inherently low-order methods. However, it would seem that in the most general case it is difficult to develop high-order schemes for control problems. For example, in the COS method, this difficulty can be traced to that of reconstructing the Fourier coefficients after numerically determining the optimal control at discrete points in the physical space. Consequently, in this paper we focus on FST or CONV techniques, which use straightforward procedures to move between Fourier space and the physical space (and vice versa).
We illustrate the behavior of our algorithm by comparing various implementations of FST/CONV on some model option-pricing examples, particularly European and Bermudan options. In addition, we demonstrate the use of the monotone scheme methods on a realistic asset allocation problem. Our main conclusion is that for problems with complex controls, where we can expect fairly low-order convergence to the solution, a small change to standard FST or CONV methods can be made that guarantees monotonicity, at least to within a user-specified tolerance. This does not alter the order of convergence in this case; hence, we can ensure a monotone scheme with only a slightly increased setup cost. After the initialization, the complexity per step of the monotone method is the same as that of the standard FST/CONV algorithm.
The remainder of this paper is as follows. In the next section, we describe our optimal control problem in a general setting. Section 3 describes existing Fourier methods and contrasts them with our new monotone Fourier method presented in Section 4. The monotone algorithm for solving optimal control problems is then given in Section 5, with properties of the algorithm and proofs given in Section 6. Wraparound is an important issue for Fourier methods, particularly in the case of our control problems. Our method of minimizing such an error is described in Section 7. Section 8 presents two numerical examples used to stress test the monotone algorithm. This is followed by an application of our algorithm to the multiperiod mean–variance optimal asset allocation problem, a general optimal control problem well suited to our monotone methods. The paper ends with our conclusions and topics for future research.
2 General control formulation
In this section, we describe our optimal control problem in a general setting. Consider a set of intervention (or monitoring) times
with the inception time of the investment and the terminal time. For simplicity, we specify the set of intervention times to be equidistant, that is, for each .
Let and (with ) denote the instants before and after the th monitoring time . We define a value function with domain (we restrict our attention to one-dimensional problems for ease of exposition) that satisfies
with a partial integrodifferential operator. At we find an optimal control via
where is an intervention operator and is the set of admissible controls.
It is more natural to rewrite these equations going backward in time, , that is, in terms of time to completion. In this case, the value function is and satisfies
Here the control , and now refers to the set of backward intervention times
A typical intervention operator has the form
As an example, in the context of portfolio allocation, we can interpret as a rebalancing rule. In general, there can also be cashflows associated with the decision process, as in the case of variable annuities. However, for simplicity, we will ignore such a generalization in this paper, and we will instead assume that the intervention operator has the form (2.6). In our asset allocation example (described later), the cashflows are modeled by updating a path-dependent variable.
3 Convolution and Fourier space time-stepping methods
In this section, we derive the FST and the closely related CONV technique in an intuitive fashion. This will allow us to contrast these methods with the monotone technique developed in the next section. For ease of exposition, we will continue to restrict our attention to one-dimensional problems. However, there is no difficulty generalizing this approach to the multidimensional case. In a financial context, it is often the case that the variable , where is an asset price.
3.1 Green’s functions
A solution of the PIDE (2.4),
can be represented in terms of the Green’s function of the PIDE, a function typically of the form . However, in many cases this will have the form ; we will assume this to hold in our problems. More formally, we make the following assumptions, which we assume hold in the rest of this work.
Assumption 3.1 (The form of Green’s function).
The Green’s function can be written as
where is known in closed form, and is independent of .
If we view the Green’s function in Assumption 3.1 as a scaled conditional probability density , then our assumption is that only depends on and via their difference . This assumption holds for Lévy processes (independent and stationary increments), but it does not hold, for example, for a Heston stochastic volatility model or for mean-reverting Ornstein–Uhlenbeck processes (but see Surkov (2010), Zhang et al (2012) and Shan (2014) for possible workarounds). The -monotonicity modifications described in this paper also hold when we do not have , but at the price of reduced efficiency. This is discussed later, in Section 4.2. The second assumption, that we know the Fourier transform of our Green’s function in closed form, holds, for example, in situations where the characteristic function of the underlying stochastic process is known. In the case of Lévy processes, the Lévy–Khintchine formula provides such an explicit representation of the characteristic function.
From Assumption 3.1, the exact solution of our PIDE is then
The Green’s function has a number of important properties (Garroni and Menaldi 1992). For this work, the two properties
are particularly important.11 1 For the examples in this paper the constant is explicitly given (in each example) in Section 3.1.1. is a constant independent of . These properties are formally proven in Garroni and Menaldi (1992), but they can also be deduced from the interpretation of the Green’s function as a scaled probability density.
We define the Fourier transform pair for the Green’s function as
with a closed-form expression for being available.
As is typically the case, we assume that the Green’s function decays to zero as , ie, is negligible outside a region . Choosing and , we localize the computational domain for the integral in (3.2) so that . We can therefore replace the Fourier transform pair (3.4) by its Fourier series equivalent:
Note that the Fourier series (3.5) implies a periodic extension of , that is, . The localization assumption also then implies that is periodically extended.
Let and choose points , by
Then, the integral in (3.7) can be approximated by a quadrature rule with weights , giving
is the discrete Fourier transform (DFT) of . In the following, we will consider two cases for the weights : the trapezoidal rule and Simpson’s quadrature. Substituting (3.8) and (3.9) into (3.7) and truncating the infinite sum to then gives
Thus, is the inverse DFT of the product .
In summary, we can obtain a discrete set of values for the solution by first going to the Fourier domain by constructing its Fourier transform using a set of quadrature weights and then returning to the physical domain by convolution of with the Fourier transform of the Green’s function. The cost is then that of doing a single FFT and inverse FFT (iFFT).
There are four significant approximations in these steps: localization of the computational domain; representation of the Green’s function by a truncated Fourier series; a periodic extension of the solution; and approximation of the integral in (3.7) by a quadrature rule. For details of the effect of the errors from these approximations we refer the reader to the discussion in Lord et al (2008).
3.1.1 Examples of Green’s functions
We consider the generic PIDE (which would be typical of a problem where the underlying asset follows a jump diffusion):
where , , , , are constants, and is the jump-size density. If, for example, , where is the risk-free rate, then this is the option-pricing equation, while if , then this PIDE arises in asset allocation.
with being the complex conjugate of . Integrating (3.13) gives
from which we can deduce that the Fourier transform of the Green’s function is
Two common choices for jump-size density are the double exponential (with , , constants)
and the lognormal (with , constants)
In the case of a double exponential jump distribution, we have
while in the case of a lognormal jump-size distribution, we have
which means that in these cases becomes
3.2 The FST/CONV algorithms
The FST and CONV algorithms are described using the previous approximations. Let be the vector of solution values just after , and let be the vector of quadrature weights:
Furthermore, let , , be a linear interpolation operator:
The full FST/CONV algorithm applied to a control problem is illustrated in Algorithm 3.3. We refer the reader to Lippa (2013), Ignatieva et al (2018) and Huang et al (2017) for applications in finance.
Algorithm 3.3 (FST/CONV Fourier method).
is the Hadamard product of vectors .
Input: number of time steps and initial solution
for do time-step loop
, optimal control
In Jackson et al (2008), the authors describe their FST method in slightly different terms. There, they use a continuous Fourier transform to convert the PIDE into Fourier space. The PIDE in physical space then reduces to a linear first-order differential equation in Fourier space that can be solved in closed form (as in Section 3.1.1). In this way, the method is able to produce exact pricing results between monitoring dates (if any) of an option, using a continuous domain. In practice, using a discrete computational domain leads to approximations, as a discrete Fourier transform is used to approximate the continuous Fourier transform.
4 An -monotone Fourier method
Our monotone Fourier method proceeds in a similar fashion to the previous section, but it is based on a slightly different philosophy. We begin by discretizing the value function, and then generate a continuous approximation of the value function by assuming linear basis functions. Given this approximation, we carry out an exact integration of the convolution integral. We can then truncate the series approximation of this integral so that monotonicity holds to within a certain tolerance (using a truncation parameter to keep track of the number of terms). Typically, the convergence to the integral is exponential in the series truncation parameter, so it is inexpensive to specify a small tolerance. The key idea here is that the number of terms required to accurately determine the projection of the Green’s function onto a given set of linear basis functions can be larger than the number of basis functions.
An additional important point is that, after the initial setup cost, the complexity per step is the same as for the standard FST and CONV methods. This requires only a small change to existing FST or CONV codes in order to guarantee monotonicity. The desirable property of this method is that monotonicity can be guaranteed (to within a small tolerance) independently of the number of FST (CONV) grid nodes or time-step size.
4.1 A monotone scheme
We proceed as follows. As before, we assume a localized computational domain:
and discretize this problem on the grid :
where and with and . Setting , we can now represent the solution as a linear combination:
where the are piecewise linear basis functions, ie,
Here we have used , a property that follows from the properties of linear basis functions. Setting , for gives
as the averaged projection of the Green’s function onto the basis functions . Note that, for this projection, , since the exact Green’s function has for all , and of course . Therefore, the scheme (4.4) is monotone for any .
Remark 4.1 (Green’s function availability in closed form).
If the Green’s function is available in closed form, rather than just its Fourier transform, then (4.6) can be used to compute the terms directly, as, for example, in Tanskanen and Lukkarinen (2003). However, in general this will require a numerical integration. If the Fourier transform of the Green’s function is known, we will derive a technique to efficiently compute to an arbitrary level of accuracy.
4.2 Approximating the monotone scheme
The scheme (4.4) is monotone, since the weights given in (4.6) are nonnegative. However it is only possible for us to approximate these weights, and this prevents us from guaranteeing monotonicity. In this subsection, we show how we overcome this issue.
Recall that our starting point is that , the Fourier series of the Green’s function, is known in closed form. We then replaced our Green’s function by its localized, periodic approximation,
This is then approximated by truncating the series.
A key point is that the truncation of the projection of the Green’s function does not have to use the same number of terms as the number of basis functions. That is, set , with defined in (4.2) and for . Suppose we now truncate the Fourier series for the projected linear basis form for to terms. Let denote the use of a truncated Fourier series with truncation parameter for a fixed value of so that the Fourier series (4.8) truncates to
Using the notation , we have
so our sequence is periodic.
Remark 4.2 (Efficient computation of the projections).
For , we define as the DFT of the :
that is, the basis function is integrated over a grid of size , and so is larger than the grid spacing on the grid. As , there is no error in evaluating these integrals (projections) for a fixed value of . For any finite , there is an error due to the use of a truncated Fourier series.
Again, we emphasize that the truncation for the Fourier series representation of the projection of the Green’s function in (4.9) does not have to use the same number of terms () as the discrete convolution (). Instead, we can take a very accurate expansion of the Green’s function projection and then translate this back to the coarse grid using (4.12). There is no further loss of information in this last step. As remarked above, we only use the Fourier representation of to carry out the discrete convolution, ie, a dense matrix–vector multiplication, efficiently. The discrete convolution in Fourier space is exactly equivalent to the discrete convolution in physical space, assuming periodic extensions.
Remark 4.3 (Assumption 3.1 revisited).
The assumption permits fast computation of a dense matrix–vector multiplication using an FFT. As mentioned earlier, this assumption holds for Lévy processes, but it does not hold, for example, for a Heston stochastic volatility model. However, the basic idea of projection of the Green’s function onto linear basis functions can be used even if Assumption 3.1 does not hold. The cost, in this case, is a loss of computational efficiency. As an example, in the case of the Heston stochastic volatility model, we have a closed form for the characteristic function, but here the Green’s function has the form , where is the variance and , where is the asset price. In this case, we can use an FFT effectively in the -direction, but not in the -direction.
Remark 4.4 (Relation to the COS method).
In the COS method, the solution is also expanded in a Fourier series. This gives an exponential convergence of the entire algorithm for smooth , which in turn requires that we have a highly accurate Fourier representation of . However, suppose is obtained by applying an impulse control using a numerical optimization method at discrete points on a previous step, using linear interpolation (the only interpolation method that is monotone in general). In that case, we will not have an accurate representation of the Fourier series of . In addition, it does not seem possible to ensure monotonicity for the COS method. So far, we have only assumed that the can be expanded in terms of piecewise linear basis functions. This property can be used to guarantee monotonicity. However, convergence will be slower than the COS method if the solution is smooth.
Remark 4.5 (Piecewise constant basis functions).
The equations and the previous discussion in this section also hold if our basis functions are piecewise constant functions, ie, basis functions that are nonzero over . In this case, computing the integral in (4.7) gives
with the subsequent equations also requiring slight modifications.
4.3 Computing the monotone scheme
In order to ensure our monotone approach is effective, it remains to compute the discrete convolution (4.4) efficiently. For the DFT pair for and , we recall that and so
Suppose we write as a DFT:
where the last equation follows from the classical orthogonality properties of th roots of unity.
From (4.14), we have
which we recognize as the inverse DFT of .
Remark 4.6 (Monotonicity).
Equations (4.18) and (4.4) are algebraic identities (assuming periodic extensions). Hence, if we use (4.18) to advance the solution, then this is algebraically identical to using (4.4) to advance the solution. Thus, we can analyze the properties of (4.18) by analyzing (4.4). In particular, if , then the scheme is monotone.
5 A monotone algorithm for solution of the control problem
In this section, we describe our monotone algorithm for the control problem (2.4), (2.5). Let be the vector of values of our solution just after , as defined in (3.20), and let be the linear interpolation operator defined as in (3.21). Let
Let us assume that our Green’s function is not an explicit function of , and that we instead have and the time steps are all constant, ie,
In this case, we can compute only once. If these two assumptions do not hold, then would have to be recomputed frequently, and hence our algorithm for ensuring monotonicity becomes more costly.
Algorithm 5.1 describes the computation of . Here, we test for monotonicity (up to a small tolerance) by minimizing the effect of any negative weights, by requiring that
The test for accuracy of the projection occurs by the comparison
Both monotonicity and convergence tests are scaled by so that these quantities are bounded as for all (the Green’s function becomes unbounded as , but the integral of the Green’s function is bounded by unity). In addition, the monotonicity test scales by in order to eliminate the number of time steps from our monotonicity bounds. This is discussed further in Section 6.
Algorithm 5.1 (Initialization of the monotone Fourier method).
Require: closed-form expression for , the Fourier transform of the Green’s function
Input: , ,
Let and compute .
for , , until convergence do construct accurate
if () and () then
break from for loop convergence test
end for end accurate loop
Output: weights , in Fourier domain.
In Algorithm 5.1, the test on line 5 will ensure that monotonicity holds to a user-specified tolerance and the test on line 6 ensures accuracy of the projections. The complete monotone algorithm for the control problem is given in Algorithm 5.4.
Remark 5.2 (Convergence of Algorithm 5.1).
Remark 5.3 (Complexity).
The complexity of using (4.18) to advance the time (excluding the cost of determining an optimal control) is operations, roughly the same as the usual FST/CONV methods.
Algorithm 5.4 (Monotone Fourier method).
Require: Weights , for in Fourier domain (from Algorithm 5.1)
Input: number of time steps , initial solution
for do time-step loop
end for end time-step loop
5.1 Convergence of truncated Fourier series for the projected Green’s functions
Since the Green’s function for (3.11) is a smooth function for any finite , we can expect uniform convergence of the Fourier series to the exact Green’s function, assuming that . This can also be seen from the exponential decay of the Fourier coefficients, which we demonstrate in this section. Since the exact Green’s function is nonnegative, the projected Green’s function (4.8) converges to a nonnegative value at every point . Consider the case of the truncated projection on linear basis functions
and . The error in the truncated series is then
Bounding the sum gives
Noting that , and from (5.4) and , we have
so that this test is usually satisfied to within round-off for .
6 Properties of the monotone Fourier method
In this section, we prove a number of properties satisfied by our -monotone Fourier algorithm. The main properties include stability and a type of -discrete comparison principle.
Let be a constant such that the exact Green’s function satisfies . Then, for all ,
For , , we have
Theorem 6.2 ( stability).
Assume that is computed using Algorithm 5.1, that is computed from
Then, for every , we have
From (6.1), we obtain
which then implies
From Lemma 6.1, we get
Remark 6.3 (Jump condition).
Lemma 6.5 (Minimum value of solution).
and, since this is valid for any , using (6.6) we obtain
Theorem 6.7 (-discrete comparison principle).
Let , . Then,
we then have
Hence, using the definition of the intervention operator (2.6), we obtain
where . Since and , the result follows. ∎
Remark 6.9 (Continuously observed impulse control problems).
By determining the optimal control at each time step, we can apply our monotone Fourier method to the continuously observed impulse control problem:
This is effectively a method whereby the optimal control is applied explicitly, as in Chen and Forsyth (2008). Using the methods developed in this paper combined with those from Chen and Forsyth (2008), it is straightforward to show that the -monotone Fourier technique is -stable and consistent in the viscosity sense as . The -monotone Fourier method is also monotone to , where is the discretization parameter. Thus, it is possible to show convergence to the viscosity solution using the results in Barles and Souganidis (1991) extended as in Azimzadeh et al (2018), using the -monotonicity property as in Bokanowski et al (2018).
7 Minimization of wraparound error
The use of the convolution form for our solution (4.18) is rigorously correct for a periodic extension of the solution and the Green’s function. In normal option-pricing applications, the wraparound error due to periodic extension causes little error. However, in control applications, the values used in the optimization step (2.5) may be near the ends of the grid, and thus large errors may result (Lippa 2013; Ruijter et al 2013; Ignatieva et al 2018). Hence, we need to consider methods to reduce errors associated with wraparound.
In order to minimize the effect of wraparound, we proceed in the following manner. Given the localized problem on with nodes, we construct an auxiliary grid with nodes, on the domain , where
with . We construct and store the DFT of the projection of the Green’s function , , on this auxiliary grid. We then replace line 4 in Algorithm 5.4 by applying the DFT to the solution on the auxiliary grid
where is an asymptotic form of the solution, which we assume to be available by financial reasoning. On the auxiliary grid near , we simply extend the solution by the constant value at , which is expected to generate a small error, since the grid spacing (in terms of ) is very small. We then carry out lines 4 and 5 of Algorithm 5.4 on the auxiliary grid and generate by discarding all the values on the auxiliary grid that are not on the original grid (as these are contaminated by wraparound errors). The errors incurred by using extensions (7.2a) and (7.2b) can be made small by choosing and sufficiently large.
(Use of asymptotic form to reduce wraparound error) Use of the above technique necessitates some changes to the proof of Theorem 6.7. However, the main result is the same, with adjustments to some of the constants in the bounds. This is a tedious algebraic exercise, which we omit.
(Additional complexity to reduce wraparound) For a one-dimensional problem, the complexity for one time step is
where is the number of nodes in the original grid. In the case of the path-dependent problem in Section 9, if there are nodes in the -direction and nodes in the bond direction, then the complexity for one time step is .
8 Numerical examples
8.1 European option
|Expiry time||0.25 years|
|Initial asset price||100|
|(a) Monotone methods|
|Piecewise linear||Piecewise constant|
Consider a European option written on an underlying stock whose price follows a jump–diffusion process. Denote by the random number representing the jump multiplier, so that when a jump occurs we have . The risk-neutral process followed by is
where denotes the expectation operator. Here, is the increment of a Wiener process, is the risk-free rate, is the volatility, is a Poisson process with positive intensity parameter , and are independent and identically distributed positive random variables. The density function , , is assumed to be double exponential (Kou and Wang 2004):
with the expectation
Given that a jump occurs, is the probability of an upward jump and is the probability of a downward jump.
The price of a European call option with is then given as the solution to
The Green’s function for this problem is given in Section 3.1.1.
|(a) Monotone methods|
|Piecewise linear||Piecewise constant|
The particular parameters for this test are given in Table 1, with the results appearing in Table 2. All methods obtain smooth second-order convergence, with the exception of the FST/CONV Simpson rule method, which gives fourth-order convergence, due to the higher-order quadrature method. This is to be expected in this case, since there is a node at the strike. Increasing and alters the last two digits of the results in the table. This is due to the effects of both localizing the problem to and FFT wraparound.
In order to stress these Fourier methods, we repeat this example using an expiry time of . Since the Green’s function in the physical space converges to a delta function as , we can expect this to be challenging for Fourier methods, as a large number of terms will be required in the Fourier series in order to get an accurate representation of the Green’s function in the physical space. The results for this test are shown in Table 3. The monotone method with piecewise linear basis functions gives reasonable results for all grid sizes. The standard FST/CONV methods are quite poor, except for very large numbers of nodes. Indeed, using Simpson’s rule on coarse grids even results in values larger than at , which violates the provable bound for a call option.
This phenomenon can be explained by examining Figure 1, which shows the projection of the Green’s functions for the monotone method (piecewise linear basis function) and the truncated Green’s function for the FST/CONV method. The projection of the Green’s function for the monotone method in Figure 1(b) clearly has the expected properties: very peaked near and nonnegative for all . In contrast, the FST/CONV numerical Green’s function is oscillatory and negative for some values of . Figure 2 compares the FST/CONV (trapezoidal) solution with the monotone (piecewise linear) solution, on a coarse grid with 512 nodes. The monotone solution can never produce a value less than zero (to within the tolerance). Note that monotonicity is clearly violated for the FST/CONV solution, with negative values for a call option. The oscillations are even more pronounced if Simpson’s quadrature is used for the FST/CONV method.
An estimate of the error in computing the projected Green’s function is given in (5.6). We can see that a very small time step affects the exponent in (5.6). For the extreme case of , (the problem in Table 1), we observe that, for , in Algorithm 5.1 is approximately , indicating that a very high accuracy projection can be achieved under extreme situations. For the same problem (2056 nodes) with , we find that in Algorithm 5.1 is approximately for .
From these tests, we can conclude both that the monotone method is robust for all time-step sizes and that for smooth problems and large time steps the monotone method exhibits a slower rate of convergence than high-order techniques, as expected.
8.2 Bermudan option with nonproportional discrete dividends
|Expiry time||10 years|
|Initial asset price||100|
|Monitoring frequency||1.0 years|
|(a) Monotone methods|
|Piecewise linear||Piecewise constant|
Let us now assume that we have the same underlying process (8.1) as in the previous subsection, except that the density function for is assumed to be normal:
with expectation . Rather than a European option, we will now consider a Bermudan put option that can be early exercised at fixed monitoring times . In addition, the underlying asset pays a fixed dividend amount at , ie, immediately after the early exercise opportunity in forward time. Between monitoring dates, the option price is given by (8.4). At monitoring dates, we have the condition
The expression in (8.6) ensures that the no-arbitrage condition holds, ie, the dividend cannot be larger than the stock price, taking into account the localized grid. Linear interpolation is used to evaluate the option value in (8.6). The parameters for this problem are listed in Table 4, with numerical results given in Table 5. All methods perform similarly, with second-order convergence. We can see that once we use a linear interpolation to impose the control there is no benefit, in terms of convergence order, to using a high-order method.
9 Multiperiod mean–variance optimal asset allocation problem
In this section, we give an example of a realistic problem with complex controls: the multiperiod mean–variance optimal asset allocation problem. Here we consider the case of an investor with a portfolio consisting of a bond index and a stock index. The amount invested in the stock index follows the process under the objective measure
with the double exponential jump-size distribution (8.2), while the amount in the bond index follows
The investor injects cash at time , with total wealth at time being . Let be the total wealth before cash injection. It turns out that, in the multiperiod mean–variance case, in some circumstances, it is optimal to withdraw cash from the portfolio (Cui et al 2014; Dang and Forsyth 2016). Denote this optimal cash withdrawal as . The total wealth after cash injection and withdrawal is then
We then select an amount to invest in the bond, so that
Since only cash withdrawals are allowed, we have . The control at rebalancing time consists of the pair . That is, after withdrawing from the portfolio we rebalance to a portfolio with in stock and in bonds. A no-leverage and no-shorting constraint is enforced by
In order to determine the mean–variance optimal solution to this asset allocation problem, we make use of the embedding result (Li and Ng 2000; Zhou and Li 2000). The mean–variance optimal strategy can be posed as
where can be viewed as a parameter that traces out the efficient frontier.
be the discounted future contributions to the portfolio at time . If
then the optimal strategy is to withdraw cash from the portfolio, and invest the remainder in the risk-free asset. This is optimal in this case since (Cui et al 2012; Dang and Forsyth 2016), which is the minimum of problem (9.6).
In the following, we will refer to any cash withdrawn from the portfolio as a surplus or free cashflow (Bauerle and Grether 2015). For the sake of argument, we will assume that the surplus cash is invested in a risk-free asset but does not contribute to the computation of the terminal mean and variance. Other possibilities are discussed in Dang and Forsyth (2016).
The solution of (9.6) is the so-called pre-commitment solution. We can interpret the pre-commitment solution in the following way. At , we decide which Pareto point is desirable (ie, a point on the efficient frontier). This fixes the value of . At any time , we can regard the optimal policy as the time-consistent solution to the problem of minimizing the expected quadratic loss with respect to the fixed target wealth , which can be viewed as a useful practical objective function (Vigna 2014; Menoncin and Vigna 2017).
9.1 Optimal control problem
A brief overview of the PIDE for the solution of the mean–variance optimal control problem is given below (we refer the reader to Dang and Forsyth (2014) for additional details).
Let the value function with be defined as
Let the set of observation times backward in time be . For , satisfies
on the localized domain .
9.2 Computational details
We solve problem (9.9) combined with the optimal control (9.12) on the localized domain . We discretize in the -direction using an equally spaced grid with nodes, and in the -direction using an unequally spaced grid with nodes. Set and denote the discrete solution at by
Let be a two-dimensional linear interpolation operator acting on the discrete solution values . Given the solution at , we use Algorithm 5.4 to advance the solution to . For the mean–variance problem, we extend this algorithm to approximate (9.11), as described in Algorithm 9.1.
Algorithm 9.1 (Advance time ).
Require: , , (from Algorithm 5.1)
for do advance time loop
end for end advance time loop
In order to advance the solution from to , we approximate the solution to the optimal control problem (9.12). The optimal control is approximated by discretizing the candidate control using the discretized grid and exhaustive search:
This algorithm converges to the solution of the original control problem as . This can be proved using similar steps to the finite-difference case (Dang and Forsyth 2014). For brevity, we omit the proof.
(Practical implementation enhancements) As noted by several authors, since the Green’s function and the solution are real, the Fourier coefficients satisfy symmetry relations. Hence, and need to be computed and stored only for . It is also possible to arrange the step in line 2 of Algorithm 9.1 and the optimal control step of (9.13) so that only a single interpolation error is introduced at each node. Note that the Fourier series representation of the Green’s function is only used to compute the projection of the Green’s function onto linear basis functions. After this initial step, we use FFTs only to efficiently carry out a dense matrix–vector multiplication (the convolution) at each step. Use of the FFT here is algebraically identical to carrying out the convolution in the physical space. The only approximation being used in this step is the periodic extension of the solution.
9.3 Numerical example
The data for this problem is given in Table 6. It was determined by fitting to the monthly returns from the Center for Research in Security Prices (CRSP) through Wharton Research Data Services, for the period January 1926 to December 2015. We use the monthly CRSP value-weighted (capitalization weighted) total return index (“vwretd”), which includes all distributions for all domestic stocks trading on major US exchanges, and the monthly ninety-day Treasury bill return index from CRSP. Both this index and the equity index are in nominal terms, so we adjust them for inflation by using the US Consumer Price Index (also supplied by CRSP). We use real indexes, since investors saving for retirement are focused on real (not nominal) wealth goals.
|Expiry time||30 years|
|Real interest rate||0.00827|
|Asymptotic form ,|
As a first test, we fix , and then increase the number of nodes in the -direction () and in the -direction (). We use the monotone scheme, with linear basis functions. In Table 7, we show the value function and the mean and standard deviation of the final wealth, which are of practical importance. The value function shows smooth second-order convergence, which is to be expected. Even though the optimal control is correct only to order (since we optimize by discretizing the controls and using exhaustive search), the value function is correct to (since it is an extreme point).
We expect that the derived quantities , , which are based on the controls computed as a by-product of computing the value function, should show a lower-order convergence. Recall that these quantities are evaluated by storing the controls and then solving a linear PIDE. In fact, we do see somewhat erratic convergence for these quantities. As an independent check, we used the stored controls from solving for the value function (on the finest grid), and then carried out Monte Carlo simulations to directly compute the mean and standard deviation of the final wealth. The results are shown in Table 8.
Of more practical interest is the following computation. In Table 9, we show the results obtained by rebalancing to a constant weight in equities at each monitoring date. We specify that the portfolio is rebalanced to 0.60 in stocks and 0.40 in bonds (a common default recommendation). We then solve for the value function using the monotone Fourier method, allowing to vary, but fixing the expected value so that is the same as for the 60:40 constant proportion strategy. This is done by using a Newton iteration, where each evaluation of the residual function requires a solution for the value function and the expected value equation. The results of this test are shown in Table 10. In this case, fixing the mean and allowing to vary results in smooth convergence of the standard deviation. From a practical point of view, we can see that the optimal strategy has the same expected value as the constant proportion strategy, but the standard deviation is reduced from 512 to 241, and the median of the optimal strategy is 936 compared with a median of for the constant proportion strategy. A heat map of the optimal strategy is shown in Figure 3.
Many problems in finance give rise to discretely monitored complex control problems. In many cases, the optimal controls are not of a simple bang-bang type. A numerical procedure must then be used to determine the optimal control at discrete points in the physical domain. In these situations, there is little hope of obtaining a high-order accurate solution after the control is applied. If we desire a monotone scheme, which increases robustness and reliability for our computations, then we are limited to the use of linear interpolation; hence, we can obtain at most second-order accuracy.
Traditional FST/CONV methods assume knowledge of the Fourier transform of the Green’s function, but then approximate this function by a truncated Fourier series. As a result, these methods are not monotone. Instead, when the Fourier transform of the Green’s function is known, we carry out a preprocessing step by projecting the Green’s function (in the physical space) onto a set of linear basis functions. These integrals can then be computed to within a specified tolerance, allowing us to guarantee a monotone scheme to within the tolerance. This monotone scheme is robust to small time steps, which is observably not the case for the standard FST/CONV methods, and indeed this lack of robustness is a major pitfall of the latter.
When the Green’s function depends on time only through the time-step size and the monitoring dates for the control are equally spaced (which is typically the case), the final monotone algorithm has the same complexity per step as the original FST/CONV algorithms, and the same order of convergence for smooth control problems. It is a simple process to add this preprocessing step to existing FST/CONV software. This results in more robust and more reliable algorithms for optimal stochastic control problems.
Declaration of interest
The authors report no conflicts of interest. The authors alone are responsible for the content and writing of the paper.
This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC).
The results presented in Section 9 were calculated based on data from Historical Indexes, © 2015 Center for Research in Security Prices (CRSP), The University of Chicago Booth School of Business. Wharton Research Data Services (WRDS) was used in preparing this article. This service and the data available thereon constitute valuable intellectual property and trade secrets of WRDS and/or its third-party suppliers.
- Alonso-Garcia, J., Wood, O., and Ziveyi, J. (2018). Pricing and hedging guaranteed minimum withdrawal benefits under a general Lévy framework using the COS method. Quantitative Finance 18(6), 1049–1075 (https://doi.org/10.1080/14697688.2017.1357832).
- Angelini, F., and Herzel, S. (2014). Delta hedging in discrete time under stochastic interest rate. Journal of Computational and Applied Mathematics 259, 385–393 (https://doi.org/10.1016/j.cam.2013.06.022).
- Azimzadeh, P., Bayraktar, E., and Labahn, G. (2018). Convergence of implicit schemes for Hamilton–Jacobi–Bellman quasi-variational inequalities. SIAM Journal on Control and Optimization 56(6), 3994–4016.
- Barles, G., and Souganidis, P. (1991). Convergence of approximation schemes for fully nonlinear equations. Asymptotic Analysis 4, 271–283.
- Bauer, D., Kling, A., and Russ, J. (2008). A universal pricing framework for guaranteed minimum benefits in variable annuities. ASTIN Bulletin 38, 621–651 (https://doi.org/10.1017/S0515036100015312).
- Bauerle, N., and Grether, S. (2015). Complete markets do not allow free cash flow streams. Mathematical Methods of Operations Research 81, 137–146 (https://doi.org/10.1007/s00186-014-0489-2).
- Bokanowski, O., Picarelli, A., and Reisinger, C. (2018). High-order filtered schemes for time dependent second order HJB equations. ESAIM: Mathematical Modelling and Numerical Analysis 52(1), 69–97 (https://doi.org/10.1051/m2an/2017039).
- Chen, Z., and Forsyth, P. A. (2008). A numerical scheme for the impulse control formulation for pricing variable annuities with a guaranteed minimum withdrawal benefit (GMWB). Numerische Mathematik 109, 535–569 (https://doi.org/10.1007/s00211-008-0152-z).
- Chen, Z., Vetzal, K., and Forsyth, P. A. (2008). The effect of modelling parameters on the value of GMWB guarantees. Insurance: Mathematics and Economics 43(1), 165–173 (https://doi.org/10.1016/j.insmatheco.2008.04.003).
- Cong, F., and Oosterlee, C. W. (2016). Multi-period mean–variance portfolio optimization based on Monte-Carlo simulation. Journal of Economic Dynamics and Control 64, 23–38 (https://doi.org/10.1016/j.jedc.2016.01.001).
- Cui, X., Li, D., Wang, S., and Zhu, S. (2012). Better than dynamic mean–variance: time inconsistency and free cash flow stream. Mathematical Finance 22(2), 346–378 (https://doi.org/10.1111/j.1467-9965.2010.00461.x).
- Cui, X., Gao, J., Li, X., and Li, D. (2014). Optimal multi-period mean–variance policy under no-shorting constraint. European Journal of Operational Research 234, 459–468 (https://doi.org/10.1016/j.ejor.2013.02.040).
- Dai, M., Kwok, Y., and Zong, J. (2008). Guaranteed minimum withdrawal benefit in variable annuities. Mathematical Finance 184, 595–611 (https://doi.org/10.1111/j.1467-9965.2008.00349.x).
- Dang, D.-M., and Forsyth, P. A. (2014). Continuous time mean–variance optimal portfolio allocation under jump diffusion: a numerical impulse control approach. Numerical Methods for Partial Differential Equations 30, 664–698 (https://doi.org/10.1002/num.21836).
- Dang, D.-M., and Forsyth, P. A. (2016). Better than pre-commitment mean–variance portfolio allocation strategies: a semi-self-financing Hamilton–Jacobi–Bellman equation approach. European Journal of Operational Research 250, 827–841 (https://doi.org/10.1016/j.ejor.2015.10.015).
- Deng, G., Dulaney, T., McCann, C., and Yan, M. (2017). Efficient valuation of equity-indexed annuities under Lévy processes using Fourier cosine series. The Journal of Computational Finance 21(2), 1–27 (https://doi.org/10.21314/JCF.2017.331).
- Fang, F., and Oosterlee, C. (2008). A novel pricing method for European options based on Fourier cosine series expansions. SIAM Journal on Scientific Computing 31, 826–848 (https://doi.org/10.1137/080718061).
- Forsyth, P., and Vetzal, K. (2017). Robust asset allocation for long-term target-based investing. International Journal of Theoretical and Applied Finance 20(3), 1750017 (https://doi.org/10.1142/S0219024917500170).
- Garroni, M. G., and Menaldi, J. L. (1992). Green Functions for Second Order Parabolic Integro-Differential Problems. Longman Scientific, New York.
- Huang, H.-C. (2010). Optimal multiperiod asset allocation: matching assets to liabilities in a discrete model. Journal of Risk and Insurance 77, 451–472 (https://doi.org/10.1111/j.1539-6975.2009.01350.x).
- Huang, Y., Zeng, P., and Kwok, Y. (2017). Optimal initiation of a guaranteed lifelong withdrawal with dynamic controls. SIAM Journal on Financial Mathematics 8, 804–840 (https://doi.org/10.1137/16M1089575).
- Ignatieva, K., Song, A., and Ziveyi, J. (2018). Fourier space time-stepping algorithm for valuing guaranteed minimum withdrawal benefits in variable annuities under regime-switching and stochastic mortality. ASTIN Bulletin 48, 139–169 (https://doi.org/10.1017/asb.2017.23).
- Jackson, K., Jaimungal, S., and Surkov, V. (2008). Fourier space time-stepping for option pricing with Lévy models. The Journal of Computational Finance 12(2), 1–29 (https://doi.org/10.21314/JCF.2008.178).
- Kou, S., and Wang, H. (2004). Option pricing under a double exponential jump diffusion model. Management Science 50(9), 1178–1192 (https://doi.org/10.1287/mnsc.1030.0163).
- Li, D., and Ng, W.-L. (2000). Optimal dynamic portfolio selection: multiperiod mean–variance formulation. Mathematical Finance 10, 387–406 (https://doi.org/10.1111/1467-9965.00100).
- Lippa, J. (2013). A Fourier space time-stepping approach applied to problems in finance. MMath Thesis, University of Waterloo.
- Lord, R., Fang, F., Bervoets, R., and Oosterlee, C. (2008). A fast and accurate FFT based method for pricing early-exercise options under Lévy processes. SIAM Journal on Scientific Computing 30, 1678–1705 (https://doi.org/10.1137/070683878).
- Menoncin, F., and Vigna, E. (2017). Mean–variance target based optimisation for defined contribution pension schemes in a stochastic framework. Insurance: Mathematics and Economics 76, 172–184 (https://doi.org/10.1016/j.insmatheco.2017.08.002).
- Obermann, A. (2006). Convergent difference schemes for degenerate elliptic and parabolic equations: Hamilton–Jacobi equations and free boundary problems. SIAM Journal of Numerical Analysis 44(2), 879–895 (https://doi.org/10.1137/S0036142903435235).
- Pooley, D., Forsyth, P., and Vetzal, K. (2003). Numerical convergence properties of option pricing PDEs with uncertain volatility. IMA Journal of Numerical Analysis 23, 241–267 (https://doi.org/10.1093/imanum/23.2.241).
- Remillard, B., and Rubenthaler, S. (2013). Optimal hedging in discrete time. Quantitative Finance 13, 819–825 (https://doi.org/10.1080/14697688.2012.745012).
- Ruijter, M., Oosterlee, C., and Albers, R. (2013). On the Fourier cosine series expansion method for stochastic control problems. Numerical Linear Algebra with Applications 20, 598–625 (https://doi.org/10.1002/nla.1866).
- Shan, C. (2014). Commodity options pricing by the Fourier space time-stepping method. MMath Essay, University of Waterloo.
- Surkov, A. (2010). Option pricing using Fourier space time-stepping framework. PhD Thesis, University of Toronto.
- Tanskanen, A., and Lukkarinen, J. (2003). Fair valuation of path-dependent participating life insurance contracts. Insurance: Mathematics and Economics 33, 595–609 (https://doi.org/10.1016/j.insmatheco.2003.08.003).
- Vigna, E. (2014). On efficiency of mean–variance based portfolio selection in defined contribution pension schemes. Quantitative Finance 14, 237–258 (https://doi.org/10.1080/14697688.2012.708778).
- Zhang, B., Grezelak, L., and Oosterlee, C. (2012). Efficient pricing of commodity options with early-exercise under the Ornstein–Uhlenbeck process. Applied Numerical Mathematics 62, 91–111 (https://doi.org/10.1016/j.apnum.2011.10.005).
- Zheng, W., and Kwok, Y. (2014). Fourier transform algorithms for pricing and hedging discretely sampled exotic variance products and volatility derivatives under additive processes. The Journal of Computational Finance 18(2), 3–30 (https://doi.org/10.21314/JCF.2014.278).
- Zhou, X. Y., and Li, D. (2000). Continuous-time mean–variance portfolio selection: a stochastic LQ framework. Applied Mathematics and Optimization 42, 19–33 (https://doi.org/10.1007/s002450010003).
Only users who have a paid subscription or are part of a corporate subscription are able to print or copy content.
You are currently unable to print this content. Please contact [email protected] to find out more.
You are currently unable to copy this content. Please contact [email protected] to find out more.
Copyright Infopro Digital Limited. All rights reserved.
You may share this content using our article tools. Printing this content is for the sole use of the Authorised User (named subscriber), as outlined in our terms and conditions - https://www.infopro-insight.com/terms-conditions/insight-subscriptions/
If you would like to purchase additional rights please email [email protected]
Copyright Infopro Digital Limited. All rights reserved.
You may share this content using our article tools. Copying this content is for the sole use of the Authorised User (named subscriber), as outlined in our terms and conditions - https://www.infopro-insight.com/terms-conditions/insight-subscriptions/
If you would like to purchase additional rights please email [email protected]