Journal of Computational Finance
ISSN:
14601559 (print)
17552850 (online)
Editorinchief: Christoph Reisinger
εmonotone Fourier methods for optimal stochastic control in finance
Peter A. Forsyth and George Labahn
Need to know
 Current Fourier methods (FST/CONV/COS) are not necessarily monotone
 We devise a preprocessing step for FST/CONV methods which are monotone.to a user specified tolerance
 The resulting methods can be used safely for optimal control problems in finance.
Abstract
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 loworder (usually at most secondorder) 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 secondorder accuracy for smooth problems.
Introduction
1 Introduction
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; AlonsoGarcia 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, finitedifference (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 longterm 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 timestepping (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; AlonsoGarcia et al 2018; Huang et al 2017) and equityindexed annuities (Deng et al 2017).
Fourier methods have a number of advantages over finitedifference and other methods; primarily, there are no timestepping 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 finitedifference methods have difficulties because of crossderivative 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 $A$ has a larger payoff at the terminal time than contract $B$. Then a monotone scheme generates a value for contract $A$ that is always larger than the value of contract $B$, 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 closedform 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 userspecified 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 secondorder 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 userspecified 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 $\epsilon $discrete comparison property, that is, given a tolerance $\epsilon $, a discrete comparison (also called arbitrage inequality) holds to $O(\epsilon )$, independent of the number of discretization nodes and time steps. This is similar in spirit to the $\epsilon $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 timestep 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 bangbang 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 bangbang type, and we expect that such good convergence properties will not hold. In addition, in the pathdependent 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 secondorder accurate at a monitoring date.
While monotone schemes have good numerical properties, they appear to be inherently loworder methods. However, it would seem that in the most general case it is difficult to develop highorder 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 optionpricing 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 loworder convergence to the solution, a small change to standard FST or CONV methods can be made that guarantees monotonicity, at least to within a userspecified 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 ${t}_{n}$
$$\widehat{?}\equiv \{{t}_{0}\le \mathrm{\cdots}\le {t}_{M}\},$$  (2.1) 
with ${t}_{0}=0$ the inception time of the investment and ${t}_{M}=T$ the terminal time. For simplicity, we specify the set of intervention times to be equidistant, that is, ${t}_{n}{t}_{n1}=\mathrm{\Delta}t=T/M$ for each $n$.
Let ${t}_{n}^{}={t}_{n}\epsilon $ and ${t}_{n}^{+}={t}_{n}+\epsilon $ (with $\epsilon \to {0}^{+}$) denote the instants before and after the $n$th monitoring time ${t}_{n}$. We define a value function $\widehat{v}(x,t)$ with domain $x\in \mathbb{R}$ (we restrict our attention to onedimensional problems for ease of exposition) that satisfies
${\widehat{v}}_{t}+\mathcal{L}\widehat{v}=0,t\in ({t}_{n}^{+},{t}_{n+1}^{}),$  (2.2) 
with $\mathcal{L}$ a partial integrodifferential operator. At ${t}_{n}\in \widehat{?}$ we find an optimal control $\widehat{c}(x,{t}_{n})$ via
$\widehat{v}(x,{t}_{n}^{})=\underset{\widehat{c}\in \mathbb{Z}}{inf}\mathcal{M}(\widehat{c})\widehat{v}(x,{t}_{n}^{+}),$  (2.3) 
where $\mathcal{M}(\widehat{c})$ is an intervention operator and $\mathbb{Z}$ is the set of admissible controls.
It is more natural to rewrite these equations going backward in time, $\tau =Tt$, that is, in terms of time to completion. In this case, the value function is $v(x,\tau )=\widehat{v}(x,Tt)$ and satisfies
$${v}_{\tau}\mathcal{L}v=0,\tau \in ({\tau}_{n}^{+},{\tau}_{n+1}^{}),$$  (2.4)  
$$v(x,{\tau}_{n}^{+})=\underset{c}{inf}\mathcal{M}(c)v(x,{\tau}_{n}^{}),{\tau}_{n}\in ?.$$  (2.5) 
Here the control $c(x,\tau )=\widehat{c}(x,T\tau )$, and $?$ now refers to the set of backward intervention times
$?\equiv \{{\tau}_{0}\le \mathrm{\cdots}\le {\tau}_{M}\}\mathit{\hspace{1em}}\text{with}{\tau}_{0}=0,{\tau}_{M}=T\mathit{\hspace{1em}}\text{and}\mathit{\hspace{1em}}{\tau}_{n}=T{t}_{Mn}.$ 
A typical intervention operator has the form
$\mathcal{M}(c)v(x,{\tau}_{n}^{})=v(x+\mathrm{\Gamma}(x,{\tau}_{n}^{},c),{\tau}_{n}^{}).$  (2.6) 
As an example, in the context of portfolio allocation, we can interpret $\mathrm{\Gamma}(x,{\tau}_{i}^{},c)$ 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 pathdependent variable.
3 Convolution and Fourier space timestepping 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 onedimensional 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 $x=\mathrm{log}(S)\in (\mathrm{\infty},\mathrm{\infty})$, where $S$ is an asset price.
3.1 Green’s functions
A solution of the PIDE (2.4),
${v}_{\tau}\mathcal{L}v=0,\tau \in ({\tau}_{n},{\tau}_{n+1}],$ 
can be represented in terms of the Green’s function of the PIDE, a function typically of the form $g=g(x,{x}^{\prime},\mathrm{\Delta}\tau )$. However, in many cases this will have the form $g=g(x{x}^{\prime},\mathrm{\Delta}\tau )$; 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
$g(x,{x}^{\prime},\mathrm{\Delta}\tau )$  $=g(x{x}^{\prime},\mathrm{\Delta}\tau )$  
$={\displaystyle {\int}_{\mathrm{\infty}}^{\mathrm{\infty}}}G(\omega ,\mathrm{\Delta}\tau ){\mathrm{e}}^{2\pi \mathrm{i}\omega (x{x}^{\prime})}d\omega ,$  (3.1) 
where $G\mathit{}\mathrm{(}\omega \mathrm{,}\mathrm{\Delta}\mathit{}\tau \mathrm{)}$ is known in closed form, and $G\mathit{}\mathrm{(}\omega \mathrm{,}\mathrm{\Delta}\mathit{}\tau \mathrm{)}$ is independent of $\mathrm{(}x\mathrm{,}{x}^{\mathrm{\prime}}\mathrm{)}$.
Remark 3.2.
If we view the Green’s function in Assumption 3.1 as a scaled conditional probability density $f$, then our assumption is that $f(y\mid x)$ only depends on $x$ and $y$ via their difference $f(y\mid x)=f(yx)$. 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 meanreverting Ornstein–Uhlenbeck processes (but see Surkov (2010), Zhang et al (2012) and Shan (2014) for possible workarounds). The $\epsilon $monotonicity modifications described in this paper also hold when we do not have $g(x,{x}^{\prime},\mathrm{\Delta}\tau )=g(x{x}^{\prime},\mathrm{\Delta}\tau )$, 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
$v(x,\tau +\mathrm{\Delta}\tau )={\displaystyle {\int}_{\mathbb{R}}}g(x{x}^{\prime},\mathrm{\Delta}\tau )v({x}^{\prime},\tau )d{x}^{\prime}.$  (3.2) 
The Green’s function has a number of important properties (Garroni and Menaldi 1992). For this work, the two properties
${\int}_{\mathbb{R}}}g(x,\mathrm{\Delta}\tau )dx={C}_{1}\le 1\mathit{\hspace{1em}}\text{and}\mathit{\hspace{1em}}g(x,\mathrm{\Delta}\tau )\ge 0$  (3.3) 
are particularly important.^{1}^{1} 1 For the examples in this paper the constant ${C}_{1}$ is explicitly given (in each example) in Section 3.1.1. ${C}_{1}$ is a constant independent of $x$. 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
$$\begin{array}{cc}\hfill G(\omega ,\mathrm{\Delta}\tau )& ={\int}_{\mathrm{\infty}}^{\mathrm{\infty}}g(x,\mathrm{\Delta}\tau ){\mathrm{e}}^{2\pi \mathrm{i}\omega x}dx,\hfill \\ \hfill g(x,\mathrm{\Delta}\tau )& ={\int}_{\mathrm{\infty}}^{\mathrm{\infty}}G(\omega ,\mathrm{\Delta}\tau ){\mathrm{e}}^{2\pi \mathrm{i}\omega x}d\omega ,\hfill \end{array}\}$$  (3.4) 
with a closedform expression for $G(\omega ,\mathrm{\Delta}\tau )$ being available.
As is typically the case, we assume that the Green’s function $g(x,\mathrm{\Delta}\tau )$ decays to zero as $x\to \mathrm{\infty}$, ie, $g(x,\mathrm{\Delta}\tau )$ is negligible outside a region $x\in [A,A]$. Choosing $$ and ${x}_{\mathrm{max}}>A$, we localize the computational domain for the integral in (3.2) so that $x\in [{x}_{\mathrm{min}},{x}_{\mathrm{max}}]$. We can therefore replace the Fourier transform pair (3.4) by its Fourier series equivalent:
$$\begin{array}{cc}\hfill G({\omega}_{k},\mathrm{\Delta}\tau )& \simeq {\int}_{{x}_{\mathrm{min}}}^{{x}_{\mathrm{max}}}g(x,\mathrm{\Delta}\tau ){\mathrm{e}}^{2\pi \mathrm{i}{\omega}_{k}x}dx,\hfill \\ \hfill \widehat{g}(x,\mathrm{\Delta}\tau )& =\frac{1}{P}\sum _{k=\mathrm{\infty}}^{\mathrm{\infty}}G({\omega}_{k},\mathrm{\Delta}\tau ){\mathrm{e}}^{2\pi \mathrm{i}{\omega}_{k}x},\hfill \end{array}\}$$  (3.5) 
with $P={x}_{\mathrm{max}}{x}_{\mathrm{min}}$ and ${\omega}_{k}=k/P$. Here the scaling factors in (3.5) are selected to be consistent with the scaling in (3.4). The solution of the PIDE (3.2) is then approximated as
$v(x,\tau +\mathrm{\Delta}\tau )\simeq {\displaystyle {\int}_{{x}_{\mathrm{min}}}^{{x}_{\mathrm{max}}}}\widehat{g}(x{x}^{\prime},\mathrm{\Delta}\tau )v({x}^{\prime},\tau )d{x}^{\prime}.$  (3.6) 
Note that the Fourier series (3.5) implies a periodic extension of $\widehat{g}$, that is, $\widehat{g}(x+P,\tau )=\widehat{g}(x,\tau )$. The localization assumption also then implies that $v(x,\tau )$ is periodically extended.
Substituting the Fourier series (3.5) into (3.6) gives our approximate solution as
$v(x,\tau +\mathrm{\Delta}\tau )\simeq {\displaystyle \frac{1}{P}}{\displaystyle \sum _{k=\mathrm{\infty}}^{\mathrm{\infty}}}G({\omega}_{k},\mathrm{\Delta}\tau ){\mathrm{e}}^{2\pi \mathrm{i}{\omega}_{k}x}{\displaystyle {\int}_{{x}_{\mathrm{min}}}^{{x}_{\mathrm{max}}}}v({x}^{\prime},\tau ){\mathrm{e}}^{2\pi \mathrm{i}{\omega}_{k}{x}^{\prime}}d{x}^{\prime}.$  (3.7) 
Let $\mathrm{\Delta}x=P/N$ and choose points $\{{x}_{j}\}$, $\{{x}_{j}^{\prime}\}$ by
${x}_{j}={\widehat{x}}_{0}+j\mathrm{\Delta}x,{x}_{j}^{\prime}={\widehat{x}}_{0}+j\mathrm{\Delta}x\mathit{\hspace{1em}}\text{for}j=\frac{1}{2}N,\mathrm{\dots},\frac{1}{2}N1.$ 
Then, the integral in (3.7) can be approximated by a quadrature rule with weights ${w}_{\mathrm{\ell}}$, giving
${\int}_{{x}_{\mathrm{min}}}^{{x}_{\mathrm{max}}}}v({x}^{\prime},\tau ){\mathrm{e}}^{2\pi \mathrm{i}{\omega}_{k}{x}^{\prime}}d{x}^{\prime$  $\simeq {\displaystyle \sum _{\mathrm{\ell}=N/2}^{N/21}}{w}_{\mathrm{\ell}}v({x}_{\mathrm{\ell}}^{\prime},\tau )\mathrm{exp}\left(2\pi \mathrm{i}{\displaystyle \frac{k}{P}}{x}_{\mathrm{\ell}}^{\prime}\right)\mathrm{\Delta}x$  
$=P\mathrm{exp}\left(2\pi \mathrm{i}{\displaystyle \frac{k}{P}}{\widehat{x}}_{0}\right)V({\omega}_{k},\tau ),$  (3.8) 
where
$V({\omega}_{k},\tau )={\displaystyle \frac{1}{N}}{\displaystyle \sum _{\mathrm{\ell}=N/2}^{N/21}}{w}_{\mathrm{\ell}}v({x}_{\mathrm{\ell}}^{\prime},\tau ){\mathrm{e}}^{2\pi \mathrm{i}k\mathrm{\ell}/N}$  (3.9) 
is the discrete Fourier transform (DFT) of $\{{w}_{j}v({x}_{j}^{\prime},\tau )\}$. In the following, we will consider two cases for the weights ${w}_{\mathrm{\ell}}$: the trapezoidal rule and Simpson’s quadrature. Substituting (3.8) and (3.9) into (3.7) and truncating the infinite sum to $k\in [\frac{1}{2}N,\frac{1}{2}N1]$ then gives
$v({x}_{j},\tau +\mathrm{\Delta}\tau )$  $\simeq {\displaystyle \frac{1}{P}}{\displaystyle \sum _{k=N/2}^{N/21}}\mathrm{exp}\left(2\pi \mathrm{i}{\displaystyle \frac{k}{P}}{\widehat{x}}_{0}\right)G({\omega}_{k},\mathrm{\Delta}\tau )\mathrm{exp}\left({\displaystyle \frac{2\pi \mathrm{i}kj}{N}}\right)$  
$\mathrm{\hspace{1em}\hspace{1em}}\times P\mathrm{exp}\left(2\pi \mathrm{i}{\displaystyle \frac{k}{P}}{\widehat{x}}_{0}\right)V({\omega}_{k},\tau )$  
$={\displaystyle \sum _{k=N/2}^{N/21}}G({\omega}_{k},\mathrm{\Delta}\tau )V({\omega}_{k},\tau )\mathrm{exp}\left({\displaystyle \frac{2\pi \mathrm{i}kj}{N}}\right).$  (3.10) 
Thus, $\{v({x}_{j},\tau +\mathrm{\Delta}\tau )\}$ is the inverse DFT of the product $\{G({\omega}_{k},\mathrm{\Delta}\tau )V({\omega}_{k},\tau )\}$.
In summary, we can obtain a discrete set of values for the solution $v$ by first going to the Fourier domain by constructing its Fourier transform $V$ using a set of quadrature weights and then returning to the physical domain by convolution of $V$ 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):
${v}_{\tau}=\frac{1}{2}{\sigma}^{2}{v}_{xx}+(\mu \frac{1}{2}{\sigma}^{2}\lambda \kappa ){v}_{x}(\rho +\lambda )v+\lambda {\displaystyle {\int}_{\mathrm{\infty}}^{+\mathrm{\infty}}}v(x+y)f(y)dy,$  (3.11) 
where $\sigma $, $\mu $, $\lambda $, $\kappa $, $\rho $ are constants, and $f(y)$ is the jumpsize density. If, for example, $\rho =\mu =r$, where $r$ is the riskfree rate, then this is the optionpricing equation, while if $\rho =0$, then this PIDE arises in asset allocation.
Let
$$\begin{array}{cc}\hfill v(x,\tau )& ={\int}_{\mathrm{\infty}}^{\mathrm{\infty}}V(\omega ,\tau ){\mathrm{e}}^{2\pi \mathrm{i}\omega x}d\omega ,\hfill \\ \hfill f(y,\tau )& ={\int}_{\mathrm{\infty}}^{\mathrm{\infty}}F(\omega ,\tau ){\mathrm{e}}^{2\pi \mathrm{i}\omega y}d\omega .\hfill \end{array}\}$$  (3.12) 
Substituting (3.12) into (3.11) gives
$$V{(\omega ,\tau )}_{\tau}=\mathrm{\Psi}(\omega )V(\omega ,\tau ),$$  (3.13) 
where
$$\mathrm{\Psi}(\omega )=(\frac{1}{2}{\sigma}^{2}{(2\pi \omega )}^{2}+(\mu \lambda \kappa \frac{1}{2}{\sigma}^{2})(2\pi \mathrm{i}\omega )(\rho +\lambda )+\lambda \overline{F}(\omega )),$$  (3.14) 
with $\overline{F}(\omega )$ being the complex conjugate of $F(\omega )$. Integrating (3.13) gives
$V(\omega ,\tau +\mathrm{\Delta}\tau )={\mathrm{e}}^{\mathrm{\Psi}(\omega )\mathrm{\Delta}\tau}V(\omega ,\tau ),$ 
from which we can deduce that the Fourier transform of the Green’s function $G(\omega ,\mathrm{\Delta}\tau )$ is
$G(\omega ,\mathrm{\Delta}\tau )={\mathrm{e}}^{\mathrm{\Psi}(\omega )\mathrm{\Delta}\tau}.$  (3.15) 
Two common choices for jumpsize density are the double exponential (with ${p}_{u}$, ${\eta}_{1}$, ${\eta}_{2}$ constants)
$$  (3.16) 
and the lognormal (with $\gamma $, $\nu $ constants)
$f(y)={\displaystyle \frac{1}{\sqrt{2\pi}\gamma}}\mathrm{exp}\left({\displaystyle \frac{{(y\nu )}^{2}}{2{\gamma}^{2}}}\right).$  (3.17) 
In the case of a double exponential jump distribution, we have
$\overline{F}(\omega )={\displaystyle \frac{{p}_{u}}{12\pi \mathrm{i}\omega /{\eta}_{1}}}+{\displaystyle \frac{1{p}_{u}}{1+2\pi \mathrm{i}\omega /{\eta}_{2}}},$  (3.18) 
while in the case of a lognormal jumpsize distribution, we have
$\overline{F}(\omega )=\mathrm{exp}(2(\pi \mathrm{i}\omega \nu {(\pi \omega \gamma )}^{2})).$  (3.19) 
From (3.13) and (3.15), we obtain
$G(0,\mathrm{\Delta}\tau )={\mathrm{e}}^{\rho \mathrm{\Delta}\tau},$ 
which means that in these cases ${C}_{1}={\int}_{\mathbb{R}}g(x,\mathrm{\Delta}\tau )dx$ becomes
${C}_{1}=\{\begin{array}{cc}{\mathrm{e}}^{r\mathrm{\Delta}\tau},\hfill & \text{option pricing},\hfill \\ 1,\hfill & \text{mean\u2013variance asset allocation}.\hfill \end{array}$ 
3.2 The FST/CONV algorithms
The FST and CONV algorithms are described using the previous approximations. Let ${({v}^{n})}^{+}$ be the vector of solution values just after ${\tau}_{n}$, and let ${w}_{\mathrm{quad}}$ be the vector of quadrature weights:
$$\begin{array}{cc}\hfill {({v}^{n})}^{+}& =[v({x}_{N/2},{\tau}_{n}^{+}),\mathrm{\dots},v({x}_{N/21},{\tau}_{n}^{+})],\hfill \\ \hfill {w}_{\mathrm{quad}}& =[w({x}_{N/2}),\mathrm{\dots},w({x}_{N/21})].\hfill \end{array}\}$$  (3.20) 
Furthermore, let ${?}_{\mathrm{\Delta}x}(x)$, ${x}_{k}\le x\le {x}_{k+1}$, be a linear interpolation operator:
${?}_{\mathrm{\Delta}x}(x){({v}^{n})}^{+}=\theta v({x}_{k},{\tau}_{n}^{+})+(1\theta )v({x}_{k+1},{\tau}_{n}^{+}),\theta ={\displaystyle \frac{({x}_{k+1}x)}{\mathrm{\Delta}x}}.$  (3.21) 
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).
$x\circ y$ is the Hadamard product of vectors $x,y$.
Require $G=\{G({\omega}_{j},\mathrm{\Delta}\tau )\}$, $j=\frac{1}{2}N,\mathrm{\dots},\frac{1}{2}N1$
 (1)
Input: number of time steps $M$ and initial solution ${({v}^{0})}^{}$
 (2)
${({v}^{0})}^{+}={inf}_{c}\mathcal{M}(c)({?}_{\mathrm{\Delta}x}(x){({v}^{0})}^{})$
 (3)
for $m=1,\mathrm{\dots},M$ do $\mathrm{\{}$timestep loop$\}$
 (4)
${V}^{m1}=\mathrm{FFT}[{w}_{\mathrm{quad}}\circ {({v}^{m1})}^{+}]$ $\{$frequency domain$\}$
 (5)
${({v}^{m})}^{}=\mathrm{iFFT}[{V}^{m1}\circ G]$ $\{$physical domain$\}$
 (6)
$v({x}_{j},{\tau}_{m}^{+})={inf}_{c}\mathcal{M}(c)({?}_{\mathrm{\Delta}x}({x}_{j}){({v}^{m})}^{})\mathit{\hspace{1em}}j=\frac{1}{2}N,\mathrm{\dots},\frac{1}{2}N1$, $\{$optimal control$\}$
 (7)
end for
Remark 3.4.
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 firstorder 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 $\epsilon $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 timestep size.
4.1 A monotone scheme
We proceed as follows. As before, we assume a localized computational domain:
$v(x,\tau +\mathrm{\Delta}\tau )={\displaystyle {\int}_{{x}_{\mathrm{min}}}^{{x}_{\mathrm{max}}}}g(x{x}^{\prime},\mathrm{\Delta}\tau )v({x}^{\prime},\tau )d{x}^{\prime}$  (4.1) 
and discretize this problem on the grid $\{{x}_{j}\},\{{x}_{j}^{\prime}\}$:
${x}_{j}={\widehat{x}}_{0}+j\mathrm{\Delta}x,{x}_{j}^{\prime}={\widehat{x}}_{0}+j\mathrm{\Delta}x,j=\frac{1}{2}N,\mathrm{\dots},\frac{1}{2}N1,$ 
where $P={x}_{\mathrm{max}}{x}_{\mathrm{min}}$ and $\mathrm{\Delta}x=P/N$ with ${x}_{\mathrm{min}}={\widehat{x}}_{0}\frac{1}{2}N\mathrm{\Delta}x$ and ${x}_{\mathrm{max}}={\widehat{x}}_{0}+\frac{1}{2}N\mathrm{\Delta}x$. Setting ${v}_{j}(\tau )=v({x}_{j},\tau )$, we can now represent the solution as a linear combination:
$v(x,\tau )\simeq {\displaystyle \sum _{j=N/2}^{N/21}}{\varphi}_{j}(x)v({x}_{j},\tau )={\displaystyle \sum _{j=N/2}^{N/21}}{\varphi}_{j}(x){v}_{j}(\tau ),$  (4.2) 
where the ${\varphi}_{j}$ are piecewise linear basis functions, ie,
${\varphi}_{j}(x)=\{\begin{array}{cc}{\displaystyle \frac{({x}_{j+1}x)}{\mathrm{\Delta}x}},\hfill & {x}_{j}\le x\le {x}_{j+1},\hfill \\ {\displaystyle \frac{(x{x}_{j1})}{\mathrm{\Delta}x}},\hfill & {x}_{j1}\le x\le {x}_{j},\hfill \\ 0,\hfill & \text{otherwise}.\hfill \end{array}$  (4.3) 
Substituting representation (4.2) into (4.1) gives
${v}_{k}(\tau +\mathrm{\Delta}\tau )$  $={\displaystyle {\int}_{{x}_{\mathrm{min}}}^{{x}_{\mathrm{max}}}}g({x}_{k}x,\mathrm{\Delta}\tau )v(x,\tau )dx$  
$={\displaystyle \sum _{j=N/2}^{N/21}}{v}_{j}(\tau ){\displaystyle {\int}_{{x}_{\mathrm{min}}}^{{x}_{\mathrm{max}}}}{\varphi}_{j}(x)g({x}_{k}x,\mathrm{\Delta}\tau )dx$  
$={\displaystyle \sum _{j=N/2}^{N/21}}{v}_{j}(\tau )\stackrel{~}{g}({x}_{k}{x}_{j},\mathrm{\Delta}\tau )\mathrm{\Delta}x,$  (4.4) 
where
$\stackrel{~}{g}({x}_{k}{x}_{j},\mathrm{\Delta}\tau )={\displaystyle \frac{1}{\mathrm{\Delta}x}}{\displaystyle {\int}_{{x}_{k}{x}_{j}\mathrm{\Delta}x}^{{x}_{k}{x}_{j}+\mathrm{\Delta}x}}{\varphi}_{kj}(x)g(x,\mathrm{\Delta}\tau )dx.$  (4.5) 
Here we have used ${\varphi}_{j}({x}_{k}x)={\varphi}_{kj}(x)$, a property that follows from the properties of linear basis functions. Setting $\mathrm{\ell}=kj$, ${y}_{\mathrm{\ell}}={x}_{k}{x}_{j}=\mathrm{\ell}\mathrm{\Delta}x$ for $\mathrm{\ell}=\frac{1}{2}N,\mathrm{\dots},\frac{1}{2}N1$ gives
$\stackrel{~}{g}({y}_{\mathrm{\ell}},\mathrm{\Delta}\tau )={\displaystyle \frac{1}{\mathrm{\Delta}x}}{\displaystyle {\int}_{{y}_{\mathrm{\ell}}\mathrm{\Delta}x}^{{y}_{\mathrm{\ell}}+\mathrm{\Delta}x}}{\varphi}_{\mathrm{\ell}}(x)g(x,\mathrm{\Delta}\tau )dx$  (4.6) 
as the averaged projection of the Green’s function onto the basis functions ${\varphi}_{\mathrm{\ell}}$. Note that, for this projection, $\stackrel{~}{g}({y}_{\mathrm{\ell}},\mathrm{\Delta}\tau )\ge 0$, since the exact Green’s function has $g(x)\ge 0$ for all $x$, and of course ${\varphi}_{\mathrm{\ell}}(x)\ge 0$. Therefore, the scheme (4.4) is monotone for any $N$.
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 $\stackrel{~}{g}({y}_{\mathrm{\ell}},\mathrm{\Delta}\tau )$ 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 $\stackrel{~}{g}({y}_{\mathrm{\ell}},\mathrm{\Delta}\tau )$ to an arbitrary level of accuracy.
4.2 Approximating the monotone scheme
The scheme (4.4) is monotone, since the weights $\stackrel{~}{g}({y}_{\mathrm{\ell}},\mathrm{\Delta}\tau )$ 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 $G$, the Fourier series of the Green’s function, is known in closed form. We then replaced our Green’s function $g(x,\mathrm{\Delta}\tau )$ by its localized, periodic approximation,
$$\widehat{g}(x,\mathrm{\Delta}\tau )=\frac{1}{P}\sum _{k=\mathrm{\infty}}^{\mathrm{\infty}}{\mathrm{e}}^{2\pi \mathrm{i}{\omega}_{k}x}G({\omega}_{k},\mathrm{\Delta}\tau ),$$  
$${\omega}_{k}=\frac{k}{P},P={x}_{\mathrm{max}}{x}_{\mathrm{min}},$$ 
and then projected the Green’s function onto the linear basis functions. Replacing $g(x,\mathrm{\Delta}\tau )$ by $\widehat{g}(x,\mathrm{\Delta}\tau )$ in (4.6), and assuming uniform convergence of the Fourier series, we integrate (4.6) term by term, resulting in
$\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau )={\displaystyle \frac{1}{P}}{\displaystyle \sum _{k=\mathrm{\infty}}^{\mathrm{\infty}}}\left({\displaystyle \frac{1}{\mathrm{\Delta}x}}{\displaystyle {\int}_{{y}_{j}\mathrm{\Delta}x}^{{y}_{j}+\mathrm{\Delta}x}}{\mathrm{e}}^{2\pi \mathrm{i}{\omega}_{k}x}{\varphi}_{j}(x)dx\right)G({\omega}_{k},\mathrm{\Delta}\tau ).$  (4.7) 
In the case of linear basis functions (4.3), we convert the complex exponential in (4.7) into trigonometric functions, with the resulting integration giving^{2}^{2} 2 For ${\omega}_{k}=0$, we take the limit ${\omega}_{k}\to 0$.
$\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau )={\displaystyle \frac{1}{P}}{\displaystyle \sum _{k=\mathrm{\infty}}^{\mathrm{\infty}}}{\mathrm{e}}^{2\pi \mathrm{i}{\omega}_{k}{y}_{j}}\left({\displaystyle \frac{{\mathrm{sin}}^{2}\pi {\omega}_{k}\mathrm{\Delta}x}{{(\pi {\omega}_{k}\mathrm{\Delta}x)}^{2}}}\right)G({\omega}_{k},\mathrm{\Delta}\tau ).$  (4.8) 
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 ${N}^{\prime}=\alpha N$, with $N$ defined in (4.2) and $\alpha ={2}^{k}$ for $k=1,2,\mathrm{\dots}$. Suppose we now truncate the Fourier series for the projected linear basis form for $\stackrel{~}{g}$ to ${N}^{\prime}$ terms. Let $\stackrel{~}{g}({y}_{k},\mathrm{\Delta}\tau ,\alpha )$ denote the use of a truncated Fourier series with truncation parameter $\alpha $ for a fixed value of $N$ so that the Fourier series (4.8) truncates to
$\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau ,\alpha )={\displaystyle \frac{1}{P}}{\displaystyle \sum _{k=\alpha N/2}^{\alpha N/21}}{\mathrm{e}}^{2\pi \mathrm{i}{\omega}_{k}j\mathrm{\Delta}x}\left({\displaystyle \frac{{\mathrm{sin}}^{2}\pi {\omega}_{k}\mathrm{\Delta}x}{{(\pi {\omega}_{k}\mathrm{\Delta}x)}^{2}}}\right)G({\omega}_{k},\mathrm{\Delta}\tau ).$  (4.9) 
Using the notation ${\stackrel{~}{g}}_{j}(\mathrm{\Delta}\tau ,\alpha )=\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau ,\alpha )$, we have
${\stackrel{~}{g}}_{j+N}(\mathrm{\Delta}\tau ,\alpha )={\stackrel{~}{g}}_{j}(\mathrm{\Delta}\tau ,\alpha ),$ 
so our sequence $\{{\stackrel{~}{g}}_{N/2}(\mathrm{\Delta}\tau ,\alpha ),\mathrm{\dots},{\stackrel{~}{g}}_{N/21}(\mathrm{\Delta}\tau ,\alpha )\}$ is periodic.
Remark 4.2 (Efficient computation of the projections).
It remains to compute the projections. For this, we need to determine the discrete convolution (4.9). Let
${Y}_{k}=\left({\displaystyle \frac{{\mathrm{sin}}^{2}\pi {\omega}_{k}\mathrm{\Delta}x}{{(\pi {\omega}_{k}\mathrm{\Delta}x)}^{2}}}\right)G({\omega}_{k},\mathrm{\Delta}\tau ),k={\displaystyle \frac{\alpha N}{2}},\mathrm{\dots},{\displaystyle \frac{\alpha N}{2}}1.$ 
Then rewriting ${\mathrm{e}}^{2\pi \mathrm{i}{\omega}_{k}j\mathrm{\Delta}x}={\mathrm{e}}^{2\pi \mathrm{i}k\mathrm{\ell}/(\alpha N)}$ with $\mathrm{\ell}=j\alpha $ and defining
${?}_{\mathrm{\ell}}={\displaystyle \frac{1}{P}}{\displaystyle \sum _{k=\alpha N/2}^{\alpha N/21}}\mathrm{exp}\left({\displaystyle \frac{2\pi \mathrm{i}k\mathrm{\ell}}{\alpha N}}\right){Y}_{k},\mathrm{\ell}={\displaystyle \frac{\alpha N}{2}},\mathrm{\dots},{\displaystyle \frac{\alpha N}{2}}1,$  (4.10) 
gives $\{{?}_{\mathrm{\ell}}\}$ as the DFT of the $\{{Y}_{k}\}$ (of size ${N}^{\prime}=\alpha N$). Consequently, using (4.9) and (4.10) yields
$\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau ,\alpha )={?}_{\mathrm{\ell}},\mathrm{\ell}=j\alpha ,j=\frac{1}{2}N,\mathrm{\dots},\frac{1}{2}N1.$  (4.11) 
Thus, the projections $\{\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau ,\alpha )\}$ are computed via a single FFT of size ${N}^{\prime}$.
For $k=\frac{1}{2}N,\mathrm{\dots},\frac{1}{2}N1$, we define $\stackrel{~}{G}({\omega}_{k},\mathrm{\Delta}\tau ,\alpha )$ as the DFT of the $\{\stackrel{~}{g}({y}_{m},\mathrm{\Delta}\tau ,\alpha )\}$:
$\stackrel{~}{G}({\omega}_{k},\mathrm{\Delta}\tau ,\alpha )={\displaystyle \frac{P}{N}}{\displaystyle \sum _{m=N/2}^{N/21}}\mathrm{exp}\left({\displaystyle \frac{2\pi \mathrm{i}mk}{N}}\right)\stackrel{~}{g}({y}_{m},\mathrm{\Delta}\tau ,\alpha ).$  (4.12) 
Note that
$\mathrm{\Delta}x={\displaystyle \frac{P}{N}}>{\displaystyle \frac{P}{{N}^{\prime}}}={\displaystyle \frac{P}{\alpha N}},$ 
that is, the basis function is integrated over a grid of size $\mathrm{\Delta}x>P/{N}^{\prime}$, and so is larger than the grid spacing on the ${N}^{\prime}$ grid. As $\alpha \to \mathrm{\infty}$, there is no error in evaluating these integrals (projections) for a fixed value of $N$. For any finite $\alpha $, 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 ($\alpha N$) as the discrete convolution ($N$). 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 $\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau ,\alpha )$ 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 $g(x,{x}^{\prime},\mathrm{\Delta}\tau )=g(x{x}^{\prime},\mathrm{\Delta}\tau )$ 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 $g=g(\nu ,{\nu}^{\prime},x{x}^{\prime},\mathrm{\Delta}\tau )$, where $\nu $ is the variance and $x=\mathrm{log}S$, where $S$ is the asset price. In this case, we can use an FFT effectively in the $x$direction, but not in the $\nu $direction.
Remark 4.4 (Relation to the COS method).
In the COS method, the solution $v(x,\tau )$ is also expanded in a Fourier series. This gives an exponential convergence of the entire algorithm for smooth $v(x,\tau )$, which in turn requires that we have a highly accurate Fourier representation of $v(x,\tau )$. However, suppose $v(x,\tau )$ 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 $v(x,\tau )$. In addition, it does not seem possible to ensure monotonicity for the COS method. So far, we have only assumed that the $v(x,\tau )$ 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 ${\varphi}_{j}$ that are nonzero over $[{x}_{j}\frac{1}{2}\mathrm{\Delta}x,{x}_{j}+\frac{1}{2}\mathrm{\Delta}x]$. In this case, computing the integral in (4.7) gives
$\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau )={\displaystyle \frac{1}{P}}{\displaystyle \sum _{k=\mathrm{\infty}}^{\mathrm{\infty}}}{\mathrm{e}}^{2\pi \mathrm{i}{\omega}_{k}{y}_{j}}\left({\displaystyle \frac{\mathrm{sin}\pi {\omega}_{k}\mathrm{\Delta}x}{\pi {\omega}_{k}\mathrm{\Delta}x}}\right)G({\omega}_{k},\mathrm{\Delta}\tau )$  (4.13) 
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 ${v}_{j}(\tau )$ and $V({\omega}_{p},\tau )$, we recall that ${x}_{j}={\widehat{x}}_{0}+j\mathrm{\Delta}x$ and so
$$\begin{array}{cc}\hfill {v}_{j}(\tau )& =\sum _{\mathrm{\ell}=N/2}^{N/21}V({\omega}_{\mathrm{\ell}},\tau ){\mathrm{e}}^{2\pi \mathrm{i}{\omega}_{l}{x}_{j}}\hfill \\ & =\sum _{\mathrm{\ell}=N/2}^{N/21}({\mathrm{e}}^{2\pi \mathrm{i}{\omega}_{\mathrm{\ell}}{\widehat{x}}_{0}})V({\omega}_{\mathrm{\ell}},\tau ){\mathrm{e}}^{2\pi \mathrm{i}j\mathrm{\ell}/N},\hfill \\ \hfill V({\omega}_{p},\tau )& =\frac{1}{N}\sum _{\mathrm{\ell}=N/2}^{N/21}{\mathrm{e}}^{2\pi \mathrm{i}{\omega}_{p}{x}_{\mathrm{\ell}}}{v}_{\mathrm{\ell}}(\tau )\hfill \\ & =\frac{1}{N}({\mathrm{e}}^{2\pi \mathrm{i}{\omega}_{p}{\widehat{x}}_{0}})\sum _{\mathrm{\ell}=N/2}^{N/21}\mathrm{exp}\left(\frac{2\pi \mathrm{i}p\mathrm{\ell}}{N}\right){v}_{\mathrm{\ell}}(\tau ).\hfill \end{array}\}$$  (4.14) 
Suppose we write $\stackrel{~}{g}({x}_{k}{x}_{j},\mathrm{\Delta}\tau )$ as a DFT:
${\stackrel{~}{g}}_{kj}(\mathrm{\Delta}\tau ,\alpha )$  $={\displaystyle \frac{1}{P}}{\displaystyle \sum _{p=N/2}^{N/21}}\stackrel{~}{G}({\omega}_{p},\mathrm{\Delta}\tau ,\alpha ){\mathrm{e}}^{2\pi \mathrm{i}(kj)p/N},$  (4.15) 
where we use (4.12) to determine $\stackrel{~}{G}({\omega}_{p},\mathrm{\Delta}\tau ,\alpha )$. Substituting (4.15) and (4.14) into (4.4), we then get
$v({x}_{k},\tau +\mathrm{\Delta}\tau )$  $=\mathrm{\Delta}x{\displaystyle \sum _{j=N/2}^{N/21}}{v}_{j}^{n}\stackrel{~}{g}({x}_{k}{x}_{j},\mathrm{\Delta}\tau ,\alpha )$  
$={\displaystyle \frac{1}{N}}{\displaystyle \sum _{p=N/2}^{N/21}}{\displaystyle \sum _{\mathrm{\ell}=N/2}^{N/21}}({\mathrm{e}}^{2\pi \mathrm{i}{\omega}_{\mathrm{\ell}}{x}_{0}})\stackrel{~}{G}({\omega}_{p},\mathrm{\Delta}\tau ,\alpha )$  
$\mathrm{\hspace{1em}\hspace{1em}}\times V({\omega}_{\mathrm{\ell}},\tau ){\mathrm{e}}^{2\pi \mathrm{i}kp/N}{\displaystyle \sum _{j=N/2}^{N/21}}\mathrm{exp}\left({\displaystyle \frac{2\pi \mathrm{i}j(\mathrm{\ell}p)}{N}}\right)$  
$={\displaystyle \sum _{p=N/2}^{N/21}}({\mathrm{e}}^{2\pi \mathrm{i}{\omega}_{p}{x}_{0}})V({\omega}_{p},\tau )\stackrel{~}{G}({\omega}_{p},\mathrm{\Delta}\tau ,\alpha )\mathrm{exp}\left({\displaystyle \frac{2\pi \mathrm{i}kp}{N}}\right),$  (4.16) 
where the last equation follows from the classical orthogonality properties of $N$th roots of unity.
From (4.14), we have
$V({\omega}_{p},\tau )={\displaystyle \frac{1}{N}}({\mathrm{e}}^{2\pi \mathrm{i}{\omega}_{p}{\widehat{x}}_{0}}){\displaystyle \sum _{\mathrm{\ell}=N/2}^{N/21}}\mathrm{exp}\left({\displaystyle \frac{2\pi \mathrm{i}p\mathrm{\ell}}{N}}\right){v}_{\mathrm{\ell}}(\tau )=({\mathrm{e}}^{2\pi \mathrm{i}{\omega}_{p}{\widehat{x}}_{0}})\stackrel{~}{V}({\omega}_{p},\tau )$  (4.17) 
with
$\stackrel{~}{V}({\omega}_{p},\tau )={\displaystyle \frac{1}{N}}{\displaystyle \sum _{\mathrm{\ell}=N/2}^{N/21}}\mathrm{exp}\left({\displaystyle \frac{2\pi \mathrm{i}p\mathrm{\ell}}{N}}\right){v}_{\mathrm{\ell}}(\tau ),$ 
the DFT of $\{{v}_{\mathrm{\ell}}(\tau )\}$. Finally, substituting (4.17) into (4.16) gives
$v({x}_{k},\tau +\mathrm{\Delta}\tau )={\displaystyle \sum _{p=N/2}^{N/21}}\stackrel{~}{V}({\omega}_{p},\tau )\stackrel{~}{G}({\omega}_{p},\mathrm{\Delta}\tau ,\alpha )\mathrm{exp}\left({\displaystyle \frac{2\pi \mathrm{i}p\mathrm{\ell}}{N}}\right),$  (4.18) 
which we recognize as the inverse DFT of $\{\stackrel{~}{V}({\omega}_{p},\tau )\stackrel{~}{G}({\omega}_{p},\mathrm{\Delta}\tau ,\alpha )\}$.
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 $\stackrel{~}{g}({x}_{k},\mathrm{\Delta}\tau ,\alpha )\ge 0$, then the scheme is monotone.
Remark 4.7 (Converting FST or CONV to monotone form).
Equation (4.18) is formally identical to (3.10). This has the practical result that any FST or CONV software can be converted to a monotone form by a preprocessing step that computes $\stackrel{~}{G}({\omega}_{p},\mathrm{\Delta}\tau ,\alpha )$ and choosing a trapezoidal rule for the integral in (3.8).
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 ${({v}^{n})}^{+}$ be the vector of values of our solution just after ${\tau}_{n}$, as defined in (3.20), and let ${?}_{\mathrm{\Delta}x}(x)$ be the linear interpolation operator defined as in (3.21). Let
$${\stackrel{~}{V}}^{n}=[\stackrel{~}{V}({\omega}_{N/2},{\tau}_{n}),\mathrm{\dots},\stackrel{~}{V}({\omega}_{N/21},{\tau}_{n})]=\mathrm{DFT}[{({v}^{n})}^{}]$$  
and  
$$\stackrel{~}{G}=[\stackrel{~}{G}({\omega}_{N/2},\mathrm{\Delta}\tau ,\alpha ),\mathrm{\dots},\stackrel{~}{G}({\omega}_{N/21},\mathrm{\Delta}\tau ,\alpha )].$$ 
Let us assume that our Green’s function is not an explicit function of $\tau $, and that we instead have $g=g(x{x}^{\prime},\mathrm{\Delta}\tau )$ and the time steps are all constant, ie, ${\tau}_{n+1}{\tau}_{n}=\mathrm{\Delta}\tau =\mathrm{const}.$
In this case, we can compute $\stackrel{~}{G}({\omega}_{k},\mathrm{\Delta}\tau ,\alpha )$ only once. If these two assumptions do not hold, then $\stackrel{~}{G}(\cdot )$ would have to be recomputed frequently, and hence our algorithm for ensuring monotonicity becomes more costly.
Algorithm 5.1 describes the computation of $\stackrel{~}{G}(\cdot )$. 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 $\mathrm{\Delta}x$ so that these quantities are bounded as $\mathrm{\Delta}\tau \to 0$ for all $\mathrm{\Delta}x$ (the Green’s function becomes unbounded as $\mathrm{\Delta}\tau \to 0$, but the integral of the Green’s function is bounded by unity). In addition, the monotonicity test scales ${\epsilon}_{1}$ by $\mathrm{\Delta}\tau /T$ 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: closedform expression for $G(\omega ,\mathrm{\Delta}\tau )$, the Fourier transform of the Green’s function
 (1)
Input: $N$, $\mathrm{\Delta}x$, $\mathrm{\Delta}\tau $
 (2)
Let $\alpha =1$ and compute $\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau ,1)$.
 (3)
for $\alpha ={2}^{k}$, $k=1,2,\mathrm{\dots}$, until convergence do $\{$construct accurate $\stackrel{~}{g}$$\}$
 (4)
 (5)
${\text{test}}_{1}={\sum}_{j}\mathrm{\Delta}x\mathrm{min}(\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau ,\alpha ),0)$ $\{$monotonicity test$\}$
 (6)
${\text{test}}_{2}={\mathrm{max}}_{j}\mathrm{\Delta}x\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau ,\alpha )\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau ,\frac{1}{2}\alpha )$ $\{$accuracy test$\}$
 (7)
if ($$) and ($$) then
 (8)
break from for loop $\{$convergence test$\}$
 (9)
end if
 (10)
end for $\{$end accurate $\stackrel{~}{g}$ loop$\}$
 (11)
Output: weights $\stackrel{~}{G}({\omega}_{j},\mathrm{\Delta}\tau ,\alpha )$, $j=\frac{1}{2}N,\mathrm{\dots},\frac{1}{2}N1$ in Fourier domain.
In Algorithm 5.1, the test on line 5 will ensure that monotonicity holds to a userspecified 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 $O(N\mathrm{log}N)$ operations, roughly the same as the usual FST/CONV methods.
Algorithm 5.4 (Monotone Fourier method).
Require: Weights $\stackrel{~}{G}=\{\stackrel{~}{G}({\omega}_{j},\mathrm{\Delta}\tau ,\alpha )\}$, for $j=\frac{1}{2}N,\mathrm{\dots},\frac{1}{2}N1$ in Fourier domain (from Algorithm 5.1)
 (1)
Input: number of time steps $M$, initial solution ${({v}^{0})}^{}$
 (2)
${({v}^{0})}^{+}={inf}_{c}\mathcal{M}(c)({?}_{\mathrm{\Delta}x}(x){({v}^{0})}^{})$
 (3)
for $m=1,\mathrm{\dots},M$ do $\{$timestep loop$\}$
 (4)
${\stackrel{~}{V}}^{m1}=\text{FFT}[{({v}^{m1})}^{+}]$ $\{$frequency domain$\}$
 (5)
${({v}^{m})}^{}=\text{iFFT}[{\stackrel{~}{V}}^{m1}\circ \stackrel{~}{G}]$ $\{$physical domain$\}$
 (6)
$v({x}_{j},{\tau}_{m}^{+})={inf}_{c}\mathcal{M}(c)({?}_{\mathrm{\Delta}x}({x}_{j}){({v}^{m})}^{})$, $j=\frac{1}{2}N,\mathrm{\dots},\frac{1}{2}N1$

$\{$optimal control$\}$
 (7)
end for $\{$end timestep 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 $\mathrm{\Delta}\tau $, we can expect uniform convergence of the Fourier series to the exact Green’s function, assuming that $\sigma >0$. 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 ${y}_{j}$. Consider the case of the truncated projection on linear basis functions
$\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau ,\alpha )={\displaystyle \frac{1}{P}}{\displaystyle \sum _{k=\alpha N/2}^{\alpha N/21}}{\mathrm{e}}^{2\pi \mathrm{i}{\omega}_{k}{y}_{j}}\left({\displaystyle \frac{{\mathrm{sin}}^{2}\pi {\omega}_{k}\mathrm{\Delta}x}{{(\pi {\omega}_{k}\mathrm{\Delta}x)}^{2}}}\right)G({\omega}_{k},\mathrm{\Delta}\tau )\mathit{\hspace{1em}}\text{with}{\omega}_{k}={\displaystyle \frac{k}{P}}$ 
and $\mathrm{\Delta}x=P/N$. The error in the truncated series is then
$\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau ,\alpha )\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau ,\mathrm{\infty})$  $={\displaystyle \frac{1}{P}}{\displaystyle \sum _{k=\alpha N/2}^{\mathrm{\infty}}}{\mathrm{e}}^{2\pi \mathrm{i}{\omega}_{k}{y}_{j}}\left({\displaystyle \frac{{\mathrm{sin}}^{2}\pi {\omega}_{k}\mathrm{\Delta}x}{{(\pi {\omega}_{k}\mathrm{\Delta}x)}^{2}}}\right)G({\omega}_{k},\mathrm{\Delta}\tau )$  
$\mathrm{\hspace{1em}\hspace{1em}}+{\displaystyle \frac{1}{P}}{\displaystyle \sum _{k=\mathrm{\infty}}^{\alpha N/21}}{\mathrm{e}}^{2\pi \mathrm{i}{\omega}_{k}{y}_{j}}\left({\displaystyle \frac{{\mathrm{sin}}^{2}\pi {\omega}_{k}\mathrm{\Delta}x}{{(\pi {\omega}_{k}\mathrm{\Delta}x)}^{2}}}\right)G({\omega}_{k},\mathrm{\Delta}\tau )$  
$\le {\displaystyle \frac{2}{P}}{\displaystyle \sum _{k=\alpha N/2}^{\mathrm{\infty}}}{\displaystyle \frac{1}{{(\pi {\omega}_{k}\mathrm{\Delta}x)}^{2}}}G({\omega}_{k},\mathrm{\Delta}\tau )$  
$\le {\displaystyle \frac{2}{P}}{\displaystyle \frac{4}{{\pi}^{2}{\alpha}^{2}}}{\displaystyle \sum _{k=\alpha N/2}^{\mathrm{\infty}}}G({\omega}_{k},\mathrm{\Delta}\tau ).$  (5.1) 
From (3.13), and noting that $\mathrm{Re}(\overline{F}(\omega ))\le 1$ (this holds since $\overline{F}(\omega )$ is the Fourier transform of a density function, and can also be seen directly from (3.18), (3.19)), we then have
$\mathrm{Re}(\mathrm{\Psi}(\omega ))$  $=\frac{1}{2}{\sigma}^{2}{(2\pi \omega )}^{2}(\rho +\lambda )+\lambda \mathrm{Re}(\overline{F}(\omega ))$  
$\le \frac{1}{2}{\sigma}^{2}{(2\pi \omega )}^{2}(\rho +\lambda )+\lambda $  
$\le \frac{1}{2}{\sigma}^{2}{(2\pi \omega )}^{2},$  (5.2) 
since $\rho \ge 0$. From (3.15) and (5.2), we have
$G(\omega ,\mathrm{\Delta}\tau )={\mathrm{e}}^{\mathrm{\Psi}(\omega )\mathrm{\Delta}\tau}\le \mathrm{exp}(\frac{1}{2}{\sigma}^{2}{(2\pi \omega )}^{2}\mathrm{\Delta}\tau ).$  (5.3) 
If we let ${C}_{4}=2{\sigma}^{2}{\pi}^{2}\mathrm{\Delta}\tau /{P}^{2}$, then (5.3) and (5.1) imply
$\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau ,\alpha )\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau ,\mathrm{\infty})\le {\displaystyle \frac{8}{P{\pi}^{2}{\alpha}^{2}}}{\displaystyle \sum _{k=\alpha N/2}^{\mathrm{\infty}}}{\mathrm{e}}^{{C}_{4}{k}^{2}}.$ 
Bounding the sum gives
$\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau ,\alpha )\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau ,\mathrm{\infty})$  $\le {\displaystyle \frac{8}{P{\pi}^{2}{\alpha}^{2}}}{\displaystyle \frac{\mathrm{exp}({C}_{4}{N}^{2}{\alpha}^{2}/4)}{1{\mathrm{e}}^{{C}_{4}N\alpha}}}.$  (5.4) 
Consider the monotonicity test in Algorithm 5.1, line 5, given by
$${\text{test}}_{1}=\sum _{j}\mathrm{\Delta}x\mathrm{min}(\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau ,\alpha ),0).$$ 
Noting that $\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau ,\mathrm{\infty})\ge 0$, and from (5.4) and ${\sum}_{j}\mathrm{\Delta}x=P$, we have
${\text{test}}_{1}\le {\displaystyle \frac{8}{{\pi}^{2}{\alpha}^{2}}}{\displaystyle \frac{\mathrm{exp}({C}_{4}{N}^{2}{\alpha}^{2}/4)}{1{\mathrm{e}}^{{C}_{4}N\alpha}}},$  (5.5) 
so that this test is usually satisfied to within roundoff for $\alpha =2,4$.
Consider now the accuracy test on line 6 of Algorithm 5.1, given by
${\text{test}}_{2}=\underset{j}{\mathrm{max}}\mathrm{\Delta}x\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau ,\alpha )\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau ,\frac{1}{2}\alpha ),$ 
which we see is bounded by
${\text{test}}_{2}$  $\le \mathrm{\Delta}x\underset{j}{\mathrm{max}}(\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau ,\alpha )\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau ,\mathrm{\infty})$  
$\mathrm{\hspace{1em}\hspace{1em}}\mathit{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}}+\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau ,\frac{1}{2}\alpha )\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau ,\mathrm{\infty}))$  
$\le \mathrm{\Delta}x\underset{j}{\mathrm{max}}(2\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau ,\frac{1}{2}\alpha )\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau ,\mathrm{\infty}))$  
$\le {\displaystyle \frac{64}{{\pi}^{2}{\alpha}^{2}}}{\displaystyle \frac{\mathrm{\Delta}x}{P}}{\displaystyle \frac{\mathrm{exp}({C}_{4}{N}^{2}{\alpha}^{2}/16)}{1{\mathrm{e}}^{{C}_{4}N\alpha /2}}}.$  (5.6) 
This test will also be satisfied for small values of $\alpha $, although these $\alpha $ should be larger than for the monotonicity test (5.5).
6 Properties of the monotone Fourier method
In this section, we prove a number of properties satisfied by our $\epsilon $monotone Fourier algorithm. The main properties include ${\mathrm{\ell}}_{\mathrm{\infty}}$ stability and a type of $\epsilon $discrete comparison principle.
Lemma 6.1.
Let ${C}_{\mathrm{1}}$ be a constant such that the exact Green’s function satisfies ${C}_{\mathrm{1}}\mathrm{=}{\mathrm{\int}}_{\mathrm{R}}g\mathit{}\mathrm{(}x\mathrm{,}\mathrm{\Delta}\mathit{}\tau \mathrm{)}\mathit{}\mathrm{d}x$. Then, for all $k$,
$\mathrm{\Delta}x{\displaystyle \sum _{j=N/2}^{N/21}}\stackrel{~}{g}({x}_{k}{x}_{j},\mathrm{\Delta}x,\alpha )={C}_{1}\mathit{\hspace{1em}}\mathit{\text{with}}\mathrm{\Delta}x={\displaystyle \frac{P}{N}}.$ 
Proof.
For ${y}_{\mathrm{\ell}}={x}_{k}{x}_{j}$, $\mathrm{\ell}=kj$, we have
$\mathrm{\Delta}x{\displaystyle \sum _{j=N/2}^{N/21}}\stackrel{~}{g}({x}_{k}{x}_{j},\mathrm{\Delta}x,\alpha )$  $={\displaystyle \frac{P}{N}}{\displaystyle \sum _{\mathrm{\ell}=N/2}^{N/21}}\stackrel{~}{g}({y}_{\mathrm{\ell}},\mathrm{\Delta}\tau ,\alpha )$  
$={\displaystyle \frac{P}{N}}{\displaystyle \sum _{\mathrm{\ell}=N/2}^{N/21}}{\displaystyle \frac{1}{P}}{\displaystyle \sum _{k=\alpha N/2}^{\alpha N/21}}{\mathrm{e}}^{2\pi \mathrm{i}{\omega}_{k}\mathrm{\ell}\mathrm{\Delta}x}\left({\displaystyle \frac{{\mathrm{sin}}^{2}\pi {\omega}_{k}\mathrm{\Delta}x}{{(\pi {\omega}_{k}\mathrm{\Delta}x)}^{2}}}\right)G({\omega}_{k},\mathrm{\Delta}\tau )$  
$={\displaystyle \frac{1}{N}}{\displaystyle \sum _{k=\alpha N/2}^{\alpha N/21}}\left({\displaystyle \frac{{\mathrm{sin}}^{2}\pi {\omega}_{k}\mathrm{\Delta}x}{{(\pi {\omega}_{k}\mathrm{\Delta}x)}^{2}}}\right)G({\omega}_{k},\mathrm{\Delta}\tau ){\displaystyle \sum _{\mathrm{\ell}=N/2}^{N/21}}\mathrm{exp}\left({\displaystyle \frac{2\pi \mathrm{i}\mathrm{\ell}k}{N}}\right)$  
$=G(0,\mathrm{\Delta}\tau )={\displaystyle {\int}_{\mathrm{\infty}}^{\mathrm{\infty}}}g(x,\mathrm{\Delta}\tau )dx$  
$={C}_{1}.$ 
∎
Theorem 6.2 (${\mathrm{\ell}}_{\mathrm{\infty}}$ stability).
Assume that $\stackrel{\mathrm{~}}{G}$ is computed using Algorithm 5.1, that ${\mathrm{(}{v}^{n}\mathrm{)}}^{\mathrm{}}$ is computed from
${({v}_{k}^{n})}^{}={\displaystyle \sum _{j=N/2}^{N/21}}\mathrm{\Delta}x{\stackrel{~}{g}}_{kj}{({v}_{j}^{n1})}^{+}$  (6.1) 
and that
${\parallel {({v}^{n})}^{+}\parallel}_{\mathrm{\infty}}\le {\parallel {({v}^{n})}^{}\parallel}_{\mathrm{\infty}}.$  (6.2) 
Then, for every $\mathrm{0}\mathrm{\le}n\mathrm{\le}M$, we have
${\parallel {({v}^{n})}^{+}\parallel}_{\mathrm{\infty}}\le {C}_{2}={\mathrm{e}}^{2{\epsilon}_{1}}{\parallel {({v}^{0})}^{}\parallel}_{\mathrm{\infty}}.$ 
Proof.
From (6.1), we obtain
${({v}_{k}^{n})}^{}$  $={\displaystyle \sum _{j=N/2}^{N/21}}\mathrm{\Delta}x{\stackrel{~}{g}}_{kj}{({v}_{j}^{n1})}^{+}$  
$={\displaystyle \sum _{j=N/2}^{N/21}}\mathrm{\Delta}x\mathrm{max}({\stackrel{~}{g}}_{kj},0){({v}_{j}^{n1})}^{+}+{\displaystyle \sum _{j=N/2}^{N/21}}\mathrm{\Delta}x\mathrm{min}({\stackrel{~}{g}}_{kj},0){({v}_{j}^{n1})}^{+},$  (6.3) 
which then implies
${({v}_{k}^{n})}^{}$  $\le {\parallel {({v}^{n1})}^{+}\parallel}_{\mathrm{\infty}}{\displaystyle \sum _{j=N/2}^{N/21}}\mathrm{\Delta}x\mathrm{max}({\stackrel{~}{g}}_{kj},0)$  
$\mathrm{\hspace{1em}\hspace{1em}}+{\parallel {({v}^{n1})}^{+}\parallel}_{\mathrm{\infty}}{\displaystyle \sum _{j=N/2}^{N/21}}\mathrm{\Delta}x\mathrm{min}({\stackrel{~}{g}}_{kj},0).$ 
From Lemma 6.1, we get
$\sum _{j=N/2}^{N/21}}\mathrm{\Delta}x\mathrm{max}({\stackrel{~}{g}}_{kj},0)={C}_{1}+{\displaystyle \sum _{j=N/2}^{N/21}}\mathrm{\Delta}x\mathrm{min}({\stackrel{~}{g}}_{kj},0),$  (6.4) 
and so, noting lines 5 and 7 in Algorithm 5.1, we have
${({v}_{k}^{n})}^{}$  $\le {\parallel {({v}^{n1})}^{+}\parallel}_{\mathrm{\infty}}\left({C}_{1}+2{\displaystyle \sum _{j=N/2}^{N/21}}\mathrm{\Delta}x\mathrm{min}({\stackrel{~}{g}}_{kj},0)\right)$  
$\le {\parallel {({v}^{n1})}^{+}\parallel}_{\mathrm{\infty}}\left({C}_{1}+2{\epsilon}_{1}{\displaystyle \frac{\mathrm{\Delta}\tau}{T}}\right).$  (6.5) 
Since (6.5) is true for any $k$, we have that
${\parallel {({v}^{n})}^{}\parallel}_{\mathrm{\infty}}\le {\parallel {({v}^{n1})}^{+}\parallel}_{\mathrm{\infty}}\left({C}_{1}+2{\epsilon}_{1}{\displaystyle \frac{\mathrm{\Delta}\tau}{T}}\right),$ 
which, combined with (6.2) and using ${C}_{1}\le 1$, gives
${\parallel {({v}^{n})}^{+}\parallel}_{\mathrm{\infty}}\le {\parallel {({v}^{n1})}^{+}\parallel}_{\mathrm{\infty}}\left(1+2{\epsilon}_{1}{\displaystyle \frac{\mathrm{\Delta}\tau}{T}}\right).$ 
Iterating the above bound and using (6.2) at $n=0$ gives
${\parallel {({v}^{n})}^{+}\parallel}_{\mathrm{\infty}}$  $\le {\parallel {({v}^{0})}^{}\parallel}_{\mathrm{\infty}}{\left(1+2{\epsilon}_{1}{\displaystyle \frac{\mathrm{\Delta}\tau}{T}}\right)}^{n}$  
$\le {\parallel {({v}^{0})}^{}\parallel}_{\mathrm{\infty}}\mathrm{exp}\left(2{\epsilon}_{1}n{\displaystyle \frac{\mathrm{\Delta}\tau}{T}}\right)$  
$\le {\parallel {({v}^{0})}^{}\parallel}_{\mathrm{\infty}}{\mathrm{e}}^{2{\epsilon}_{1}}$  
$={C}_{2}.$ 
∎
Remark 6.3 (Jump condition).
Corollary 6.4 (Stability of Algorithm 5.4).
Algorithm 5.4 is ${\mathrm{\ell}}_{\mathrm{\infty}}$stable.
Lemma 6.5 (Minimum value of solution).
Let ${\mathrm{(}{v}^{n}\mathrm{)}}^{\mathrm{+}}$ be generated using (6.1), and set
${({v}^{n})}_{\mathrm{min}}^{+}=\underset{k}{\mathrm{min}}{({v}_{k}^{n})}^{+}.$ 
If the conditions for Lemma 6.1 are satisfied and
${({v}^{n})}_{\mathrm{min}}^{+}\ge {({v}^{n})}_{\mathrm{min}}^{},$  (6.6) 
then
${({v}^{n})}_{\mathrm{min}}^{+}\ge {({v}^{0})}_{\mathrm{min}}^{}{({C}_{3})}^{n}{C}_{2}({\mathrm{e}}^{{\epsilon}_{1}}1),$ 
where ${C}_{\mathrm{2}}\mathrm{=}{\mathrm{\parallel}{\mathrm{(}{v}^{\mathrm{0}}\mathrm{)}}^{\mathrm{}}\mathrm{\parallel}}_{\mathrm{\infty}}\mathit{}{\mathrm{e}}^{\mathrm{2}\mathit{}{\epsilon}_{\mathrm{1}}}$ is given in Lemma 6.1 and
$${C}_{3}=\sum _{j=N/2}^{N/21}\mathrm{\Delta}x\mathrm{max}({\stackrel{~}{g}}_{kj},0).$$ 
Proof.
From (6.1), and using (6.3) along with the definition of ${C}_{3}$, we obtain
${({v}_{k}^{n})}^{}$  $\ge {({v}^{n1})}_{\mathrm{min}}^{+}{\displaystyle \sum _{j=N/2}^{N/21}}\mathrm{\Delta}x\mathrm{max}({\stackrel{~}{g}}_{kj},0)+{\displaystyle \sum _{j=N/2}^{N/21}}\mathrm{\Delta}x\mathrm{min}({\stackrel{~}{g}}_{kj},0){({v}_{j}^{n})}^{+}$  
$\ge {({v}^{n1})}_{\mathrm{min}}^{+}{\displaystyle \sum _{j=N/2}^{N/21}}\mathrm{\Delta}x\mathrm{max}({\stackrel{~}{g}}_{kj},0)$  
$\mathrm{\hspace{1em}\hspace{1em}}{\parallel {({v}^{n1})}^{+}\parallel}_{\mathrm{\infty}}{\displaystyle \sum _{j=N/2}^{N/21}}\mathrm{\Delta}x\mathrm{min}({\stackrel{~}{g}}_{kj},0)$  
$={({v}^{n1})}_{\mathrm{min}}^{+}{C}_{3}{\parallel {({v}^{n1})}^{+}\parallel}_{\mathrm{\infty}}\left({\displaystyle \sum _{j=N/2}^{N/21}}\mathrm{\Delta}x\mathrm{min}({\stackrel{~}{g}}_{kj},0)\right).$ 
Using Lemma 6.1 and lines 5 and 7 in Algorithm 5.1 then gives
${({v}_{k}^{n})}^{}\ge {({v}^{n1})}_{\mathrm{min}}^{+}{C}_{3}{C}_{2}{\epsilon}_{1}{\displaystyle \frac{\mathrm{\Delta}\tau}{T}}$ 
and, since this is valid for any $k$, using (6.6) we obtain
${({v}^{n})}_{\mathrm{min}}^{+}\ge {({v}^{n1})}_{\mathrm{min}}^{+}{C}_{3}{C}_{2}{\epsilon}_{1}{\displaystyle \frac{\mathrm{\Delta}\tau}{T}}.$ 
Iterating implies
${({v}^{n})}_{\mathrm{min}}^{+}$  $\ge {({v}^{0})}_{\mathrm{min}}^{+}{C}_{3}^{n}{C}_{2}{\epsilon}_{1}{\displaystyle \frac{\mathrm{\Delta}\tau}{T}}\left({\displaystyle \frac{1{C}_{3}^{n}}{1{C}_{3}}}\right)$  
$\ge {({v}^{0})}_{\mathrm{min}}^{}{C}_{3}^{n}{C}_{2}{\epsilon}_{1}{\displaystyle \frac{\mathrm{\Delta}\tau}{T}}\left({\displaystyle \frac{1{C}_{3}^{n}}{1{C}_{3}}}\right),$  (6.7) 
where we again use (6.6) in the last line. From (6.4) and the definition of ${C}_{3}$, we have
${C}_{3}={C}_{1}+{\displaystyle \sum _{j=N/2}^{N/21}}\mathrm{\Delta}x\mathrm{min}({\stackrel{~}{g}}_{kj},0)\le 1+{\epsilon}_{1}{\displaystyle \frac{\mathrm{\Delta}\tau}{T}},$  (6.8) 
where the last inequality follows from lines 5 and 7 in Algorithm 5.1 (and recalling that ${C}_{1}\le 1$). Combining (6.6)–(6.8), and noting that $n\mathrm{\Delta}\tau \le T$, gives
${({v}^{n})}_{\mathrm{min}}^{+}\ge {({v}^{0})}_{\mathrm{min}}^{+}{C}_{3}^{n}{C}_{2}({\mathrm{e}}^{{\epsilon}_{1}}1)\ge {({v}^{0})}_{\mathrm{min}}^{}{C}_{3}^{n}{C}_{2}({\mathrm{e}}^{{\epsilon}_{1}}1).$ 
∎
Remark 6.6.
Theorem 6.7 ($\epsilon $discrete comparison principle).
Suppose we have two independent discrete solutions:
$$\begin{array}{cc}\hfill {({u}^{n})}^{+}& =[u({x}_{N/2},{\tau}_{n}^{+}),\mathrm{\dots},u({x}_{+N/21},{\tau}_{n}^{+})],\hfill \\ \hfill {({w}^{n})}^{+}& =[w({x}_{N/2},{\tau}_{n}^{+}),\mathrm{\dots},w({x}_{+N/21},{\tau}_{n}^{+})],\hfill \end{array}\}$$  (6.9) 
with
${({u}^{0})}^{}\ge {({w}^{0})}^{},$ 
where the inequality is understood in the componentwise sense, and ${\mathrm{(}{u}^{n}\mathrm{)}}^{\mathrm{+}}\mathrm{,}{\mathrm{(}{w}^{n}\mathrm{)}}^{\mathrm{+}}$ are computed using Algorithm 5.4. If $\stackrel{\mathrm{~}}{G}$ is computed using Algorithm 5.1 and ${\mathrm{I}}_{\mathrm{\Delta}\mathit{}x}\mathit{}\mathrm{(}x\mathrm{)}$ is a linear interpolant, then
${({u}^{n})}^{+}{({w}^{n})}^{+}\ge {\epsilon}_{1}{\parallel {({u}^{0}{w}^{0})}^{}\parallel}_{\mathrm{\infty}}+O({\epsilon}_{1}^{2}),{\epsilon}_{1}\to 0.$  (6.10) 
Proof.
Let ${({z}^{n})}^{+}={({u}^{n})}^{+}{({w}^{n})}^{+}$, ${({z}^{n})}^{}={({u}^{n})}^{}{({w}^{n})}^{}$. Then,
${({z}_{k}^{n})}^{}={\displaystyle \sum _{j=N/2}^{N/21}}\mathrm{\Delta}x{\stackrel{~}{g}}_{kj}{({z}_{j}^{n1})}^{+}.$ 
Noting that
${z}_{j}({\tau}_{n}^{+})=\underset{c}{inf}\mathcal{M}(c)({?}_{\mathrm{\Delta}x}({x}_{j}){({u}^{m})}^{})\underset{c}{inf}\mathcal{M}(c)({?}_{\mathrm{\Delta}x}({x}_{j}){({w}^{m})}^{}),$  (6.11) 
we then have
${z}_{j}({\tau}_{n}^{+})\le \underset{c}{sup}\mathcal{M}(c){?}_{\mathrm{\Delta}x}({x}_{j})({({u}^{m})}^{}{({w}^{m})}^{}).$  (6.12) 
Hence, using the definition of the intervention operator (2.6), we obtain
${\parallel {({z}^{n})}^{+}\parallel}_{\mathrm{\infty}}\le {\parallel {({z}^{n})}^{}\parallel}_{\mathrm{\infty}}.$  (6.13) 
Similarly,
${({z}^{n})}_{\mathrm{min}}^{+}$  $=\underset{j}{\mathrm{min}}{z}_{j}({\tau}_{n}^{+})$  
$\ge \underset{j}{\mathrm{min}}\underset{c}{inf}\mathcal{M}(c){?}_{\mathrm{\Delta}x}({x}_{j})({({u}^{m})}^{}{({w}^{m})}^{})$  
$\ge {({z}^{n})}_{\mathrm{min}}^{}.$  (6.14) 
Hence, condition (6.2) of Lemma 6.1 and condition (6.6) of Lemma 6.5 are satisfied. Applying Lemma 6.5 to ${({z}^{n})}^{+}$, ${({z}^{n})}^{}$, we get
${({z}^{n})}_{\mathrm{min}}^{+}\ge {({z}^{0})}_{\mathrm{min}}^{}{({C}_{3})}^{n}{\mathrm{e}}^{2{\epsilon}_{1}}{\parallel {({u}^{0}{w}^{0})}^{}\parallel}_{\mathrm{\infty}}({\mathrm{e}}^{{\epsilon}_{1}}1),$  (6.15) 
where ${C}_{3}={\sum}_{j=N/2}^{N/21}\mathrm{\Delta}x\mathrm{max}({\stackrel{~}{g}}_{kj},0)$. Since ${({z}^{0})}_{\mathrm{min}}^{}\ge 0$ and $0\le {C}_{3}^{n}\le {\mathrm{e}}^{{\epsilon}_{1}}$, the result follows. ∎
Remark 6.8.
If Algorithm 5.1 is used to construct $\stackrel{~}{G}$ for use in Algorithm 5.4, then the $\epsilon $discrete comparison property is satisfied for any $N$, $\mathrm{\Delta}\tau $, $M$ up to order ${\epsilon}_{1}$. Since typically $\stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau ,\alpha )\to \stackrel{~}{g}({y}_{j},\mathrm{\Delta}\tau ,\mathrm{\infty})\ge 0$ exponentially in $\alpha $, in practice it is very inexpensive to make ${\epsilon}_{1}$ as small as desired.
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:
$\mathrm{max}[{v}_{\tau}\mathcal{L}v,v\underset{c}{inf}\mathcal{M}(c)v]=0.$  (6.16) 
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 $\epsilon $monotone Fourier technique is ${\mathrm{\ell}}_{\mathrm{\infty}}$stable and consistent in the viscosity sense as $\mathrm{\Delta}\tau ,\mathrm{\Delta}x\to 0$. The $\epsilon $monotone Fourier method is also monotone to $O(h)$, where $h=O(\mathrm{\Delta}x)=O(\mathrm{\Delta}\tau )$ 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 $\epsilon $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 optionpricing 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 $[{x}_{\mathrm{min}},{x}_{\mathrm{max}}]$ with $N$ nodes, we construct an auxiliary grid with ${N}^{a}=2N$ nodes, on the domain $[{x}_{\mathrm{min}}^{a},{x}_{\mathrm{max}}^{a}]$, where
${x}_{\mathrm{min}}^{a}={x}_{\mathrm{min}}\frac{1}{2}({x}_{\mathrm{max}}{x}_{\mathrm{min}})\mathit{\hspace{1em}}\text{and}\mathit{\hspace{1em}}{x}_{\mathrm{max}}^{a}={x}_{\mathrm{max}}+\frac{1}{2}({x}_{\mathrm{max}}{x}_{\mathrm{min}})$  (7.1) 
with $({x}_{\mathrm{max}}^{a}{x}_{\mathrm{min}}^{a})=2({x}_{\mathrm{max}}{x}_{\mathrm{min}})$. We construct and store the DFT of the projection of the Green’s function $\stackrel{~}{G}({\omega}_{p},\mathrm{\Delta}\tau ,\alpha )$, $p=\frac{1}{2}{N}^{a},\mathrm{\dots},\frac{1}{2}{N}^{a}1$, on this auxiliary grid. We then replace line 4 in Algorithm 5.4 by applying the DFT to the solution $v$ on the auxiliary grid
$v{({x}_{k},{\tau}_{n}^{+})}^{a}$  $=v({x}_{k},{\tau}_{n}^{+})$  $(k$  $=\frac{1}{2}N,\mathrm{\dots},\frac{1}{2}N1)$  
$=v({x}_{N/2},{\tau}_{n}^{+})$  $(k$  $=\frac{1}{2}{N}^{a},\mathrm{\dots},\frac{1}{2}N1)$  (7.2a)  
$=A({x}_{k},{\tau}_{n}^{+})$  $(k$  $=\frac{1}{2}N,\mathrm{\dots},\frac{1}{2}{N}^{a}1),$  (7.2b) 
where $A(x,\tau )$ is an asymptotic form of the solution, which we assume to be available by financial reasoning. On the auxiliary grid near $x\to \mathrm{\infty}$, we simply extend the solution by the constant value at $x={x}_{\mathrm{min}}$, which is expected to generate a small error, since the grid spacing (in terms of $S={\mathrm{e}}^{x}$) is very small. We then carry out lines 4 and 5 of Algorithm 5.4 on the auxiliary grid and generate ${({v}^{n})}^{}$ 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 ${x}_{\mathrm{min}}$ and ${x}_{\mathrm{max}}$ sufficiently large.
Remark 7.1.
(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.
Remark 7.2.
(Additional complexity to reduce wraparound) For a onedimensional problem, the complexity for one time step is
$$O({N}^{a}\mathrm{log}{N}^{a})=O(2N\mathrm{log}(2N)),$$ 
where $N$ is the number of nodes in the original grid. In the case of the pathdependent problem in Section 9, if there are ${N}_{x}$ nodes in the $\mathrm{log}S$direction and ${N}_{b}$ nodes in the bond direction, then the complexity for one time step is $O(2{N}_{b}{N}_{x}\mathrm{log}(2{N}_{x}))$.
8 Numerical examples
8.1 European option
Expiry time  0.25 years 

Strike $K$  100 
Payoff  Call 
Initial asset price ${S}_{\text{0}}$  100 
Riskfree rate $r$  0.05 
Volatility $\sigma $  0.15 
$\lambda $  0.1 
${\eta}_{\text{1}}$  3.0465 
${\eta}_{\text{2}}$  3.0775 
${p}_{u}$  0.3445 
${x}_{\mathrm{max}}$  $\mathrm{log}({S}_{\text{0}})+\text{10}$ 
${x}_{\mathrm{min}}$  $\mathrm{log}({S}_{\text{0}})\text{10}$ 
${\epsilon}_{\text{1}},{\epsilon}_{\text{2}}$  ${\text{10}}^{\text{6}}$ 
Asymptotic form $x\to \mathrm{\infty}$  $A(x)={\mathrm{e}}^{x}$ 
(a) Monotone methods  

Piecewise linear  Piecewise constant  
$?$  Value  Ratio  Value  Ratio 
${\text{2}}^{\text{9}}$  3.9808516210  3.9443958729  
${\text{2}}^{\text{10}}$  3.9753205007  3.9662547470  
${\text{2}}^{\text{11}}$  3.9739391670  4.0  3.9716756819  4.0 
${\text{2}}^{\text{12}}$  3.9735939225  4.0  3.9730282349  4.0 
${\text{2}}^{\text{13}}$  3.9735076171  4.0  3.9733662066  4.0 
${\text{2}}^{\text{14}}$  3.9734860412  4.0  3.9734506895  4.0 
(b) FST/CONV  
Trapezoidal  Simpson  
$?$  Value  Ratio  Value  Ratio 
${\text{2}}^{\text{9}}$  3.9075619850  3.9784907318  
${\text{2}}^{\text{10}}$  3.9571661688  3.9737010716  
${\text{2}}^{\text{11}}$  3.9694107823  4.1  3.9734923202  23.0 
${\text{2}}^{\text{12}}$  3.9724624589  4.0  3.9734796846  17.0 
${\text{2}}^{\text{13}}$  3.9732247908  4.0  3.9734789013  16.0 
${\text{2}}^{\text{14}}$  3.9734153372  4.0  3.9734788524  16.0 
Consider a European option written on an underlying stock whose price $S$ follows a jump–diffusion process. Denote by $\xi $ the random number representing the jump multiplier, so that when a jump occurs we have ${S}_{t}=\xi {S}_{{t}^{}}$. The riskneutral process followed by ${S}_{t}$ is
$\frac{\mathrm{d}{S}_{t}}{{S}_{{t}^{}}}}=(r\lambda \kappa )\mathrm{d}t+\sigma \mathrm{d}Z+\mathrm{d}\left({\displaystyle \sum _{i=1}^{{\pi}_{t}}}({\xi}_{i}1)\right)\mathit{\hspace{1em}}\text{with}\kappa =E[\xi ]1,$  (8.1) 
where $E[\cdot ]$ denotes the expectation operator. Here, $\mathrm{d}Z$ is the increment of a Wiener process, $r$ is the riskfree rate, $\sigma $ is the volatility, ${\pi}_{t}$ is a Poisson process with positive intensity parameter $\lambda $, and ${\xi}_{i}$ are independent and identically distributed positive random variables. The density function $f(y)$, $y=\mathrm{log}(\xi )$, is assumed to be double exponential (Kou and Wang 2004):
$$  (8.2) 
with the expectation
$E[\xi ]={\displaystyle \frac{{p}_{u}{\eta}_{1}}{{\eta}_{1}1}}+{\displaystyle \frac{(1{p}_{u}){\eta}_{2}}{{\eta}_{2}+1}}.$  (8.3) 
Given that a jump occurs, ${p}_{u}$ is the probability of an upward jump and $(1{p}_{u})$ is the probability of a downward jump.
The price of a European call option $v(x,\tau )$ with $x=\mathrm{log}S$ is then given as the solution to
$${v}_{\tau}=\frac{1}{2}{\sigma}^{2}{v}_{xx}+(r\frac{1}{2}{\sigma}^{2}\lambda \kappa ){v}_{x}(r+\lambda )v+\lambda {\int}_{\mathrm{\infty}}^{+\mathrm{\infty}}v(x+y)f(y)dy\mathit{\hspace{1em}}\text{with}v(x,0)=\mathrm{max}({\mathrm{e}}^{x}K,0).$$  (8.4) 
The Green’s function for this problem is given in Section 3.1.1.
(a) Monotone methods  

Piecewise linear  Piecewise constant  
$?$  Value  Ratio  Value  Ratio 
${\text{2}}^{\text{9}}$  0.19662316859  0.94284763015  
${\text{2}}^{\text{10}}$  0.19467436458  0.041410269769  
${\text{2}}^{\text{11}}$  0.19376651687  2.1  0.15335986938  $$8.0 
${\text{2}}^{\text{12}}$  0.19346709107  3.0  0.18477993505  3.6 
${\text{2}}^{\text{13}}$  0.19339179620  4.0  0.19127438852  4.8 
${\text{2}}^{\text{14}}$  0.19337297842  4.0  0.19284673379  4.1 
(b) FST/CONV  
Trapezoidal  Simpson  
$?$  Value  Ratio  Value  Ratio 
${\text{2}}^{\text{9}}$  0.24774086499  319.45747026  
${\text{2}}^{\text{10}}$  0.21909081933  521.62802838  
${\text{2}}^{\text{11}}$  0.18611676723  0.87  439.13444172  $$2.5 
${\text{2}}^{\text{12}}$  0.17728640855  3.7  27.002978049  0.2 
${\text{2}}^{\text{13}}$  0.18913280108  $$0.75  0.19367805822  15.0 
${\text{2}}^{\text{14}}$  0.19231903134  3.7  0.19338110881  $\text{9.0}\times {\text{10}}^{\text{4}}$ 
The particular parameters for this test are given in Table 1, with the results appearing in Table 2. All methods obtain smooth secondorder convergence, with the exception of the FST/CONV Simpson rule method, which gives fourthorder convergence, due to the higherorder quadrature method. This is to be expected in this case, since there is a node at the strike. Increasing ${x}_{\mathrm{max}}$ and ${x}_{\mathrm{min}}$ alters the last two digits of the results in the table. This is due to the effects of both localizing the problem to $[{x}_{\mathrm{min}},{x}_{\mathrm{max}}]$ and FFT wraparound.
In order to stress these Fourier methods, we repeat this example using an expiry time of $T=0.001$. Since the Green’s function in the physical space converges to a delta function as $T\to 0$, 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 ${S}_{0}$ at $S={S}_{0}=100$, 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 $x=0$ and nonnegative for all $x$. In contrast, the FST/CONV numerical Green’s function is oscillatory and negative for some values of $x$. 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.
Remark 8.1 (Error in approximating (4.6) using (4.11)).
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 $T=0.001$, $N=2056$ (the problem in Table 1), we observe that, for $\alpha =8$, ${\text{test}}_{2}$ in Algorithm 5.1 is approximately ${10}^{12}$, indicating that a very high accuracy projection can be achieved under extreme situations. For the same problem (2056 nodes) with $T=0.25$, we find that ${\text{test}}_{2}$ in Algorithm 5.1 is approximately ${10}^{16}$ for $\alpha =2$.
From these tests, we can conclude both that the monotone method is robust for all timestep sizes and that for smooth problems and large time steps the monotone method exhibits a slower rate of convergence than highorder techniques, as expected.
8.2 Bermudan option with nonproportional discrete dividends
Expiry time  10 years 

Strike $K$  100 
Payoff  Put 
Initial asset price ${S}_{\text{0}}$  100 
Riskfree rate $r$  0.05 
Volatility $\sigma $  0.15 
Dividend $D$  1.00 
Monitoring frequency $\mathrm{\Delta}\tau $  1.0 years 
$\lambda $  0.1 
$\nu $  $$1.08 
$\gamma $  0.4 
${x}_{\mathrm{max}}$  $\mathrm{log}({S}_{\text{0}})+\text{10}$ 
${x}_{\mathrm{min}}$  $\mathrm{log}({S}_{\text{0}})\text{10}$ 
${\epsilon}_{\text{1}},{\epsilon}_{\text{2}}$  ${\text{10}}^{\text{6}}$ 
Asymptotic form $x\to \mathrm{\infty}$  $A(x)=\text{0}$ 
(a) Monotone methods  

Piecewise linear  Piecewise constant  
$?$  Value  Ratio  Value  Ratio 
${\text{2}}^{\text{9}}$  24.811127744  24.806532754  
${\text{2}}^{\text{10}}$  24.789931363  24.788800257  
${\text{2}}^{\text{11}}$  24.782264461  2.8  2.4.781982815  2.6 
${\text{2}}^{\text{12}}$  24.781134292  6.8  24.781063962  7.4 
${\text{2}}^{\text{13}}$  24.780822977  3.6  24.780805394  3.6 
${\text{2}}^{\text{14}}$  24.780744620  4.0  24.780740225  4.0 
(b) FST/CONV  
Trapezoidal  Simpson  
$?$  Value  Ratio  Value  Ratio 
${\text{2}}^{\text{9}}$  24.801967268  24.802639420  
${\text{2}}^{\text{10}}$  24.787670043  24.787731820  
${\text{2}}^{\text{11}}$  24.781701225  2.4  24.781787212  2.5 
${\text{2}}^{\text{12}}$  24.780993635  8.4  24.781007785  7.6 
${\text{2}}^{\text{13}}$  24.780787811  3.4  24.780788678  3.6 
${\text{2}}^{\text{14}}$  24.780735831  4.0  24.780737159  4.3 
Let us now assume that we have the same underlying process (8.1) as in the previous subsection, except that the density function for $y=\mathrm{log}(\xi )$ is assumed to be normal:
$f(y)={\displaystyle \frac{1}{\sqrt{2\pi}\gamma}}\mathrm{exp}\left({\displaystyle \frac{{(y\nu )}^{2}}{2{\gamma}^{2}}}\right),$  (8.5) 
with expectation $E[\xi ]={\mathrm{e}}^{\nu +{\gamma}^{2}/2}$. Rather than a European option, we will now consider a Bermudan put option that can be early exercised at fixed monitoring times ${\tau}_{n}$. In addition, the underlying asset pays a fixed dividend amount $D$ at ${\tau}_{n}^{}$, 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
$$v(x,{\tau}_{n}^{+})=\mathrm{max}(v(\mathrm{log}(\mathrm{max}({\mathrm{e}}^{x}D,{\mathrm{e}}^{{x}_{\mathrm{min}}})),{\tau}_{n}^{}),P(x)),\text{with}P(x)=\mathrm{payoff}=\mathrm{max}(K{\mathrm{e}}^{x},0).$$  (8.6) 
The expression $\mathrm{max}({\mathrm{e}}^{x}D,{\mathrm{e}}^{{x}_{\mathrm{min}}})$ in (8.6) ensures that the noarbitrage 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 secondorder 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 highorder 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
$\frac{\mathrm{d}{S}_{t}}{{S}_{{t}^{}}}}=(\mu \lambda \kappa )\mathrm{d}t+\sigma \mathrm{d}Z+\mathrm{d}\left({\displaystyle \sum _{i=1}^{{\pi}_{t}}}({\xi}_{i}1)\right)$  (9.1) 
with the double exponential jumpsize distribution (8.2), while the amount in the bond index follows
$\mathrm{d}{B}_{t}=r{B}_{t}\mathrm{d}t.$  (9.2) 
The investor injects cash ${q}_{n}$ at time ${t}_{n}\in \widehat{?}$, with total wealth at time $t$ being ${W}_{t}={S}_{t}+{B}_{t}$. Let ${W}_{n}^{}={S}_{n}^{}+{B}_{n}^{}$ 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 ${c}_{n}^{*}$. The total wealth after cash injection and withdrawal is then
${W}_{n}^{+}={W}_{n}^{}+{q}_{n}{c}_{n}^{*}.$  (9.3) 
We then select an amount ${b}_{n}^{*}$ to invest in the bond, so that
${B}_{n}^{+}$  $={b}_{n}^{*}\mathit{\hspace{1em}}\text{and}\mathit{\hspace{1em}}{S}_{n}^{+}={W}_{n}^{+}{b}_{n}^{*}.$  (9.4) 
Since only cash withdrawals are allowed, we have ${c}_{n}^{*}\ge 0$. The control at rebalancing time ${t}_{n}$ consists of the pair $({b}_{n}^{*},{c}_{n}^{*})$. That is, after withdrawing ${c}_{n}^{*}$ from the portfolio we rebalance to a portfolio with ${S}_{n}^{+}$ in stock and ${B}_{n}^{+}$ in bonds. A noleverage and noshorting constraint is enforced by
$0\le {b}_{n}^{*}\le {W}_{n}^{+}.$  (9.5) 
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
$$\underset{\{({b}_{0}^{*},{c}_{0}^{*}),\mathrm{\dots},({b}_{M1}^{*},{c}_{M1}^{*})\}}{\mathrm{min}}\mathit{\hspace{1em}}E[{({W}^{*}{W}_{T})}^{2}]$$  
$$\text{subject to}\{\begin{array}{cc}({S}_{t},{B}_{t})\text{follow processes(9.1),(9.2)};t\notin \widehat{?},\hfill & \\ {W}_{n}^{+}={S}_{n}^{}+{B}_{n}^{}+{q}_{n}{c}_{n}^{*},\hfill & \\ {S}_{n}^{+}={W}_{n}^{+}{b}_{n}^{*},{B}_{n}^{+}={b}_{n}^{*};t\in \widehat{?},\hfill & \\ 0\le {b}_{n}^{*}\le {W}_{n}^{+},\hfill & \\ {c}_{n}^{*}\ge 0,\hfill & \end{array}$$  (9.6) 
where ${W}^{*}$ can be viewed as a parameter that traces out the efficient frontier.
Let
${Q}_{\mathrm{\ell}}={\displaystyle \sum _{j=\mathrm{\ell}+1}^{M1}}{\mathrm{e}}^{r({t}_{j}{t}_{\mathrm{\ell}})}{q}_{j}$  (9.7) 
be the discounted future contributions to the portfolio at time ${t}_{\mathrm{\ell}}$. If
$({W}_{n}^{}+{q}_{n})>{W}^{*}{\mathrm{e}}^{r(T{t}_{n})}{Q}_{n},$  (9.8) 
then the optimal strategy is to withdraw cash ${c}_{n}^{*}={W}_{n}^{}+{q}_{n}({W}^{*}{\mathrm{e}}^{r(T{t}_{n})}{Q}_{n})$ from the portfolio, and invest the remainder $({W}^{*}{\mathrm{e}}^{r(T{t}_{n})}{Q}_{n})$ in the riskfree asset. This is optimal in this case since $E[{({W}^{*}{W}_{T})}^{2}]=0$ (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 riskfree 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 socalled precommitment solution. We can interpret the precommitment solution in the following way. At $t=0$, we decide which Pareto point is desirable (ie, a point on the efficient frontier). This fixes the value of ${W}^{*}$. At any time $t>0$, we can regard the optimal policy as the timeconsistent solution to the problem of minimizing the expected quadratic loss with respect to the fixed target wealth ${W}^{*}$, 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 $v(x,b,\tau )$ with $\tau =Tt$ be defined as
$$v(x,b,\tau )=\underset{\{({b}_{0}^{*},{c}_{0}^{*}),\mathrm{\dots},({b}_{M}^{*},{c}_{M}^{*})\}}{inf}\{E[{(\mathrm{min}({W}_{T}{W}^{*},0))}^{2}\mid \mathrm{log}S(t)=x,B(t)=b]\}.$$  (9.9) 
Let the set of observation times backward in time be $?=\{{\tau}_{0},{\tau}_{1},\mathrm{\dots},{\tau}_{M}\}$. For $\tau \notin ?$, $v$ satisfies
$$\begin{array}{cc}\hfill {v}_{\tau}& =\mathcal{L}v+rb{v}_{b},\hfill \\ \hfill \mathcal{L}v& \equiv \frac{1}{2}{\sigma}^{2}{v}_{xx}+(\mu \frac{1}{2}{\sigma}^{2}\lambda \kappa ){v}_{x}(\mu +\lambda )v\hfill \\ & +\lambda {\int}_{\mathrm{\infty}}^{\mathrm{\infty}}v(x+y)f(y)dy,\hfill \\ \hfill v(x,b,0)& ={(\mathrm{min}({\mathrm{e}}^{x}+b{W}^{*},0))}^{2},\hfill \end{array}\}$$  (9.10) 
on the localized domain $(x,b)\in [{x}_{\mathrm{min}},{x}_{\mathrm{max}}]\times [0,{b}_{\mathrm{max}}]$.
If $g(x,\tau )$ is the Green’s function of ${v}_{\tau}=\mathcal{L}v$, then the solution of (9.10) at ${\tau}_{n+1}^{}$, given the solution at ${\tau}_{n}^{+}$, ${\tau}_{n}\in ?$, is
$v(x,b,{\tau}_{n+1}^{})={\displaystyle {\int}_{{x}_{\mathrm{min}}}^{{x}_{\mathrm{max}}}}g(x{x}^{\prime},\mathrm{\Delta}\tau )v({x}^{\prime},b{\mathrm{e}}^{rb\mathrm{\Delta}\tau},{\tau}_{n}^{+})\mathit{\hspace{1em}}\text{with}\mathrm{\Delta}\tau ={\tau}_{n+1}{\tau}_{n}.$  (9.11) 
Equation (9.11) can be regarded as a combination of a Green’s function step for the PIDE ${v}_{\tau}=\mathcal{L}v$ and a characteristic technique to handle the $rb{v}_{b}$ term. At rebalancing times ${\tau}_{n}\in ?$,
$$v(x,b,{\tau}_{n}^{+})=\underset{({b}^{*},{c}^{*})}{\mathrm{min}}v({x}^{\prime},{b}^{*},{\tau}_{n}^{})$$  
$$\text{subject to}\{\begin{array}{cc}{c}^{*}=\mathrm{max}({\mathrm{e}}^{x}+b+{q}_{Mn}{Q}_{Mn},0),\hfill & \\ {W}^{\prime}={\mathrm{e}}^{x}+b+{q}_{Mn}{c}^{*},\hfill & \\ 0\le {b}^{*}\le {W}^{\prime},\hfill & \\ {x}^{\prime}=\mathrm{log}(\mathrm{max}({W}^{\prime}{b}^{*},{\mathrm{e}}^{{x}_{\mathrm{min}}})),\hfill & \end{array}$$  (9.12) 
where ${Q}_{\mathrm{\ell}}$ is defined in (9.7).
9.2 Computational details
We solve problem (9.9) combined with the optimal control (9.12) on the localized domain $(x,b)\in [{x}_{\mathrm{min}},{x}_{\mathrm{max}}]\times [0,{b}_{\mathrm{max}}]$. We discretize in the $x$direction using an equally spaced grid with ${N}_{x}$ nodes, and in the $B$direction using an unequally spaced grid with ${N}_{b}$ nodes. Set ${B}_{\mathrm{max}}={\mathrm{e}}^{{x}_{\mathrm{max}}}$ and denote the discrete solution at $({x}_{m},{b}_{j},{\tau}_{n}^{+})$ by
$$\begin{array}{cc}\hfill {({v}_{m,j}^{n})}^{+}& =v({x}_{m},{b}_{j},{\tau}_{n}^{+}),\hfill \\ \hfill {({v}^{n})}^{+}& ={\{{({v}_{m,j}^{n})}^{+}\}}_{m={N}_{x}/2,\mathrm{\dots},{N}_{x}/21;j=1,\mathrm{\dots},{N}_{b}},\hfill \\ \hfill {({v}_{j}^{n})}^{+}& =[{({v}_{{N}_{x}/2,j}^{n})}^{+},\mathrm{\dots},{({v}_{{N}_{x}/21,j}^{n})}^{+}].\hfill \end{array}\}$$ 
Let ${\mathcal{I}}_{\mathrm{\Delta}x,\mathrm{\Delta}b}(x,b){({v}^{n})}^{}$ be a twodimensional linear interpolation operator acting on the discrete solution values ${({v}^{n})}^{}$. Given the solution at ${\tau}_{n}^{+}$, we use Algorithm 5.4 to advance the solution to ${\tau}_{n+1}^{}$. For the mean–variance problem, we extend this algorithm to approximate (9.11), as described in Algorithm 9.1.
Algorithm 9.1 (Advance time ${({v}^{n})}^{+}\to {({v}^{n+1})}^{}$).
Require: ${({v}^{n})}^{+}$, $\stackrel{~}{G}=\{\stackrel{~}{G}({\omega}_{m},\mathrm{\Delta}\tau ,\alpha )\}$, $m=\frac{1}{2}{N}_{x},\mathrm{\dots},\frac{1}{2}{N}_{x}1$ (from Algorithm 5.1)
 (1)
for $j=1,\mathrm{\dots},{N}_{b}$ do $\{$advance time loop$\}$
 (2)
${v}_{m,j}^{\text{int}}={\mathcal{I}}_{\mathrm{\Delta}x,\mathrm{\Delta}b}({x}_{m},{b}_{j}{\mathrm{e}}^{r\mathrm{\Delta}\tau}){({v}^{n})}^{+}$, $m=\frac{1}{2}{N}_{x},\mathrm{\dots},\frac{1}{2}{N}_{x}1$,
 (3)
$\stackrel{~}{V}=\text{FFT}[{v}_{j}^{\text{int}}]$
 (4)
${({v}_{j}^{n+1})}^{}=\text{iFFT}[\stackrel{~}{V}\circ \stackrel{~}{G}]$ $\{$iFFT(Hadamard product)$\}$
 (5)
end for $\{$end advance time loop$\}$
In order to advance the solution from ${\tau}_{n+1}^{}$ to ${\tau}_{n+1}^{+}$, we approximate the solution to the optimal control problem (9.12). The optimal control is approximated by discretizing the candidate control ${b}^{*}$ using the discretized $b$ grid and exhaustive search:
$$v({x}_{m},{b}_{j},{\tau}_{n}^{+})=\underset{({b}^{*},{c}^{*})}{\mathrm{min}}{\mathcal{I}}_{\mathrm{\Delta}x,\mathrm{\Delta}b}({x}^{*},{b}^{*}){({v}^{n+1})}^{}$$  
$$\text{subject to}\{\begin{array}{cc}{c}^{*}=\mathrm{max}({\mathrm{e}}^{{x}_{m}}+{b}_{j}+{q}_{Mn}{Q}_{Mn},0),\hfill & \\ {W}^{\prime}={\mathrm{e}}^{{x}_{m}}+{b}_{j}+{q}_{Mn}{c}^{*},\hfill & \\ {b}^{*}\in \{{b}_{1},\mathrm{\dots},\mathrm{min}({b}_{\mathrm{max}},{W}^{\prime})\},\hfill & \\ {x}^{*}=\mathrm{log}(\mathrm{max}({W}^{\prime}{b}^{*},{\mathrm{e}}^{{x}_{\mathrm{min}}})).\hfill & \end{array}$$  (9.13) 
This algorithm converges to the solution of the original control problem as ${N}_{x},{N}_{b}\to \mathrm{\infty}$. This can be proved using similar steps to the finitedifference case (Dang and Forsyth 2014). For brevity, we omit the proof.
Using the control determined from solving (9.9), we can determine $E[{W}_{T}]$ and standard deviation $\text{SD}[{W}_{T}]$ by solving an additional linear PIDE (for details, see Dang and Forsyth (2014)).
Remark 9.2.
(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, $\stackrel{~}{G}({\omega}_{k},\mathrm{\Delta}\tau ,\alpha )$ and $\stackrel{~}{V}$ need to be computed and stored only for ${\omega}_{k}\ge 0$. 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 valueweighted (capitalization weighted) total return index (“vwretd”), which includes all distributions for all domestic stocks trading on major US exchanges, and the monthly ninetyday 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 $T$  30 years 

Initial wealth  0 
Rebalancing frequency  Yearly 
Cash injection ${\{{q}_{i}\}}_{i=\text{0},\mathrm{\dots},\text{29}}$  10 
Real interest rate $r$  0.00827 
Volatility $\sigma $  0.14777 
$\mu $  0.08885 
$\lambda $  0.3222 
${\eta}_{1}$  4.4273 
${\eta}_{\text{2}}$  5.262 
${p}_{u}$  0.2758 
${x}_{\mathrm{max}}$  $\mathrm{log}(\text{100})+\text{5}$ 
${x}_{\mathrm{min}}$  $\mathrm{log}(\text{100})\text{10}$ 
${\epsilon}_{\text{1}},{\epsilon}_{\text{2}}$  ${\text{10}}^{\text{6}}$ 
Asymptotic form $E[{({W}_{T}{W}^{*})}^{\text{2}}]$, $x\to \mathrm{\infty}$  $A(x)=\text{0}$ 
As a first test, we fix ${W}^{*}=1022$, and then increase the number of nodes in the $x$direction (${N}_{x}$) and in the $b$direction (${N}_{b}$). We use the monotone scheme, with linear basis functions. In Table 7, we show the value function $v(0,0,T)$ and the mean $E[{W}_{T}]$ and standard deviation $\text{SD}[{W}_{T}]$ of the final wealth, which are of practical importance. The value function shows smooth secondorder convergence, which is to be expected. Even though the optimal control is correct only to order $\mathrm{\Delta}b$ (since we optimize by discretizing the controls and using exhaustive search), the value function is correct to $O{(\mathrm{\Delta}b)}^{2}$ (since it is an extreme point).
We expect that the derived quantities $E[{W}_{T}]$, $\text{SD}[{W}_{T}]$, which are based on the controls computed as a byproduct of computing the value function, should show a lowerorder 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.
${?}_{?}$  ${?}_{?}$  Value function  Ratio  $?\mathbf{[}{?}_{?}\mathbf{]}$  Ratio  $\text{??}\mathbf{[}{?}_{?}\mathbf{]}$  Ratio 

512  305  97 148.899100  N/A  824.02599269  N/A  240.73884508  N/A 
1024  609  97 042.740997  N/A  824.07104985  N/A  240.55534019  N/A 
2048  1217  97 014.471301  3.8  824.09034690  2.3  240.51245396  4.3 
4096  2433  97 007.286530  3.9  824.08961667  $$26.0  240.49691620  2.7 
8192  4865  97 005.451814  3.9  824.09295889  $$0.22  240.49585213  14.6 
${?}_{\text{???}}$  $?\mathbf{[}{?}_{?}\mathbf{]}$  $\text{??}\mathbf{[}{?}_{?}\mathbf{]}$ 

$\text{1.6}\times {\text{10}}^{\text{5}}$  824.3425 (1.55)  240.2263 
$\text{6.4}\times {\text{10}}^{\text{5}}$  823.6719 (0.78)  240.7278 
$\text{2.56}\times {\text{10}}^{\text{6}}$  824.0077 (0.39)  240.4336 
$\text{1.024}\times {\text{10}}^{\text{7}}$  824.1043 (0.19)  240.5217 
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 ${W}^{*}$ to vary, but fixing the expected value so that $E[{W}_{T}]$ 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 ${W}^{*}$ 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 $704$ for the constant proportion strategy. A heat map of the optimal strategy is shown in Figure 3.
$?\mathbf{[}{?}_{?}\mathbf{]}$  $\text{??}\mathbf{[}{?}_{?}\mathbf{]}$  $\text{??????}\mathbf{[}{?}_{?}\mathbf{]}$ 

824.10047  511.8482  704 
${?}_{?}$  ${?}_{?}$  $?\mathbf{[}{?}_{?}\mathbf{]}$  $\text{??}\mathbf{[}{?}_{?}\mathbf{]}$  Ratio 

512  305  824.10047  240.79440842  N/A 
1024  609  824.10047  240.57925928  N/A 
2048  1217  824.10047  240.52022512  3.6 
4096  2433  824.10047  240.50571976  4.1 
8192  4865  824.10047  240.50220544  4.1 
10 Conclusions
Many problems in finance give rise to discretely monitored complex control problems. In many cases, the optimal controls are not of a simple bangbang 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 highorder 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 secondorder 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 timestep 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.
Acknowledgements
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 thirdparty suppliers.
References
 AlonsoGarcia, 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 quasivariational 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/s0018601404892).
 Bokanowski, O., Picarelli, A., and Reisinger, C. (2018). Highorder 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/s002110080152z).
 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). Multiperiod mean–variance portfolio optimization based on MonteCarlo 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.14679965.2010.00461.x).
 Cui, X., Gao, J., Li, X., and Li, D. (2014). Optimal multiperiod mean–variance policy under noshorting 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.14679965.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 precommitment mean–variance portfolio allocation strategies: a semiselffinancing 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 equityindexed 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 longterm targetbased 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 IntegroDifferential 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.15396975.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 timestepping algorithm for valuing guaranteed minimum withdrawal benefits in variable annuities under regimeswitching 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 timestepping 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/14679965.00100).
 Lippa, J. (2013). A Fourier space timestepping 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 earlyexercise 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 timestepping method. MMath Essay, University of Waterloo.
 Surkov, A. (2010). Option pricing using Fourier space timestepping framework. PhD Thesis, University of Toronto.
 Tanskanen, A., and Lukkarinen, J. (2003). Fair valuation of pathdependent 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 earlyexercise 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). Continuoustime 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.
To access these options, along with all other subscription benefits, please contact [email protected] or view our subscription options here: http://subscriptions.risk.net/subscribe
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.infoproinsight.com/termsconditions/insightsubscriptions/
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.infoproinsight.com/termsconditions/insightsubscriptions/
If you would like to purchase additional rights please email [email protected]