{'fraction_non_alphanumeric': 0.04570056566149839, 'fraction_numerical': 0.01283929516135596, 'mean_word_length': 4.285760378576038, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A probabilistic numerical method for optimal multiple switching problem and application to investments in electricity generation

In this paper, we present a probabilistic numerical algorithm combining dynamic programming, Monte Carlo simulations and local basis regressions to solve non-stationary optimal multiple switching problems in infinite horizon. We provide the rate of convergence of the method in terms of the time step used to discretize the problem, of the size of the local hypercubes involved in the regressions, and of the truncating time horizon. To make the method viable for problems in high dimension and long time horizon, we extend a memory reduction method to the general Euler scheme, so that, when performing the numerical resolution, the storage of the Monte Carlo simulation paths is not needed. Then, we apply this algorithm to a model of optimal investment in power plants. This model takes into account electricity demand, cointegrated fuel prices, carbon price and random outages of power plants. It computes the optimal level of investment in each generation technology, considered as a whole, w.r.t. the electricity spot price. This electricity price is itself built according to a new extended structural model. In particular, it is a function of several factors, among which the installed capacities. The evolution of the optimal generation mix is illustrated on a realistic numerical problem in dimension eight, i.e. with two different technologies and six random factors.


Introduction
This paper presents a probabilistic numerical method for multiple switching problem with an application to a new stylized long-term investment model for electricity generation. Since electricity cannot be stored and building new plants takes several years, investment in new capacities must be decided a long time in advance if a country wishes to be able to satisfy its demand 1 . Before the worldwide liberalization of the electricity sector, electric utilities were monopolies whose objective was to plan the construction of power plants in order to satisfy demand at the minimum cost under a given constraint on the loss of load probability or on the level of energy non-served. This investment process was called generation expansion planning (GEP). Its output was mostly a given set of power plants to build for the next ten or twenty years (see [36] for a comprehensive description of the GEP methodology and related difficulties). Despite thirty years of liberalization of the electricity sector, of the recognition that GEP methods were inadequate within a market context ( [34,24]) and of an important set of alternative methods (see [28] and [22] for recent surveys on generation investment models and softwares), power utilities still heavily rely on GEP methods (see [3]). However, real option methods, which should have been the natural alternative valuation method for firms converted to a value maximizing objective, did not emerge as the method of choice. Despite the important body of literature that followed [45]'s seminal paper, [23]'s monography and implementations for electricity generation investment (see for instance [11]), real options still remain a marginal way of assessing investment decisions both in the electric sector and in the industry in general (see for instance, amongst the recurrent surveys on capital budgeting methods, [8]). Nevertheless, as shown in [44], firms tend to reproduce with heuristic constraints (such as hurdle rate or profitability index) the decision criteria given by real option methodology.
The main reason for this situation lies in the considerable mathematical difficulties involved in the conception of a tractable yet realistic real option model for electricity generation. This difficulty reflects in the literature where the main trend consists in designing a small dimensional (1 or 2) real option model to assess investment behaviour with respect to some specific variables (see for instance [1,9] for models in dimension 2 analysing the effects of uncertainty and time to build). It is still possible to find investment model in dimension 3 based on dynamic programming which are numerically tractable (see for instance [46,11] and Section 5 for comments). But, in higher dimension, because of the curse of dimensionality, investment models mainly rely on decision trees to represent random factors (see [2] for a recent typical implementation of this approach). The resulting tractability is however obtained at the expense of a crude simplification of the statistical properties of the factors.
Our approach in the present paper takes advantage of the considerable progress made in the last ten years by numerical methods for high-dimensional American options valuation problems to propose a probabilistic way to look at future electricity generation mixes. For an up-to-date state of the art on this subject, the reader is referred to the recent book [15].
In this paper, we first adapt the resolution of American option problems by Monte-Carlo methods ( [43,57]) to the more general class of optimal switching problems. The crucial choice of regression basis is done here in the light of the work of [13], so as to obtain a stable algorithm suited to highdimensional problems, aiming at the best possible numerical complexity. The memory complexity, often acknowledged as the major drawback of such a Monte Carlo approach (see [16]), is drastically slashed by generalizing the memory reduction method from [18,19,20] to any stochastic differential equation. We provide a rigorous and comprehensive analysis of the rate of convergence of our algorithm, taking advantage of the works of, most notably, [12], [55] and [29]. Note that such features as infinite horizon and non-stationarity are encompassed here. Finally, we build a long-term structural model for the spot price of electricity, extending the work of [5] and [4] in several directions (cointegrated fuels and CO 2 prices, stochastic availability rate of production capacities, etc.). This model is itself incorporated into an optimal control problem corresponding to the search for the optimal investments in electricity generation. The resolution of this problem using our algorithm is illustrated on a simple numerical example with two different technologies, leading to an eight-dimensional problem (demand, CO 2 price, and, for each technology, fuel price, random outages and the controlled installed capacity). The time evolution of the distribution of power prices and of the generation mix is illustrated on a forty-year time horizon.
To sum up, the contribution of the paper is threefold. Firstly, it provides, for a suitably chosen regression basis, a comprehensive analysis of convergence of a regression-based Monte-Carlo algorithm for a class of infinite horizon optimal multiple switching problems, large enough to handle realistic short term profit functions and investment cost structures with possible seasonality patterns. Secondly, it adapts and generalizes a memory reduction method in order to slash the amount of memory required by the Monte Carlo algorithm. Thirdly, a new stylised investment model for electricity generation is proposed, taking into account electricity demand, cointegrated fuel prices, carbon price and random outages of power plants, used as building blocks of a new structural model for the electricity spot price. A numerical resolution of this investment problem with our algorithm is illustrated on a specific example, providing, among many other outputs, an electricity spot price dynamics consistent with the investment decision process in power generation.
The outline of the paper is the following. Section 2 describes the class of optimal switching problems studied here, including the detailed list of assumptions considered. Section 3 describes the resolution algorithm and analyzes its rate of convergence, in terms of the discretization step, of the size of the local hypercubes from the regression basis, and of the truncating time horizon. Section 4 details the computational complexity of the algorithm, as well as its memory complexity, along with the construction of the memory reduction method. Finally, Section 5 introduces the extended structural model of power spot price, the investment problem, as well as an illustrated numerical resolution. Section 6 concludes the paper.

Notation
Here are some notation that will be used throughout the paper: • The notation 1 {.} stands for the indicator function.
• Throughout the paper, C > 0 denotes a generic constant whose value may differ from line to line, but which does not depend on any parameter of our scheme.
• For any stochastic process X = (X s ) s≥0 taking values in a given set X , and any (t, x) ∈ R + × X , we denote as X t,x = (X t,x s ) s≥t the stochastic process with the same dynamics as X, but starting from x at time t: X t,x t = x.
• ∀p ≥ 1, the norms . p and . Lp denote respectively the p−norm and the L p -norm: ∀x ∈ R n and any R-valued random variable X such that E [|X| p ] < ∞: x p 2 Optimal switching problem

Formulation
Fix a filtered probability space Ω, F, F = (F t ) t≥0 , P , where F satisfies the usual conditions of rightcontinuity and P-completeness. We consider the following general class of (non-stationary) optimal switching problems: where: • I α = (I α s ) s≥0 is a càd-làg, R d -valued, F-adapted piecewise constant process. It is controlled by a strategy α, described below. We suppose it can only take values into a fixed finite set which means that equation (2.1) corresponds to an optimal switching problem.
• An impulse control strategy α corresponds to a sequence (τ n , ι n ) n∈N of increasing stopping times τ n ≥ 0, and F τn -measurable random variables ι n valued in I q . Using this sequence, I α = (I α s ) s≥0 is defined as follows: Alternatively, α can be described by the sequence (τ n , ζ n ) n∈N , where ζ n := ι n − ι n−1 (and ζ 0 := 0). Using this alternative sequence, I α can be written as follows: A is the set of admissible strategies: a strategy α belongs to A if τ n → +∞ a.s. as n → ∞.
• For any (t, i) ∈ R + × I q , the set A t,i ⊂ A is defined as the subset of admissible strategies α such that I α t = i.
• f and k are R-valued measurable functions.

Assumptions
We complete the above formulation with the following relevant assumptions.

Assumption 1. [Diffusion]
The R d -valued uncontrolled process X is a diffusion process, governed by the dynamics where W is a d-dimensional Brownian motion, and b and σ are respectively R d -valued and R d×dvalued functions.
(use Burkholder-Davis-Gundy inequality and Gronwall's Lemma, see for instance [37] Theorem 4.5.4 for the even power case).

Assumption 3.
[Lipschitz&Discount] The functions f and k decrease exponentially in time: where the functionsf andk are Lipschitz continuous with linear growth: Moreover, we assume in the following that ρ > ρ 1 where ρ 1 is defined in equation (2.3).

Assumption 4. [Fixed costs] The cost function
• (triangular inequality) ∀t ∈ R + , ∀ (i, j, k) ∈ (I q ) 3 with i = j and j = k: Remark 2.3. The economic interpretations of Assumption 4 are the following: 1. There is no cost for not switching, but any switch incurs at least a positive fixed cost.
2. At any given date, it is always cheaper to switch directly from i to k than to switch first from i to j and then from j to k.
Remark 2.4. Under those standard assumptions, the value function v is well-defined and finite. Indeed, using equation In particular, the costs being positive (Assumption 4), and recalling (2.1), it holds that:

Outline of the solution
From a theoretical point of view, the value functions v i := v (., ., i), i ∈ I q from equation (2.1) are known to satisfy (under suitable conditions on f i (., .) := f (., ., i) and k, see for instance [53] in a much more general setting) the following Hamilton-Jacobi-Bellman Quasi-Variational Inequalities (HJBQVI): together with suitable limit condition.
Alternatively, the process v (t, X t , i), t ≥ 0 can be characterized as the solution of a particular Reflected Backward Stochastic Differential Equation ( [33,25]).
Moreover, the value function (2.1) satisfies the well-known dynamic programming principle, i.e., for any stopping time τ ≥ t: From a practical point of view, apart from a few simple examples in low-dimension, finding directly the solution of the HJBQVI (2.6) is usually infeasible, and the numerical PDE tools become cumbersome and inefficient in the multi-dimensional setting. Instead, probabilistic methods based on (2.7), in the spirit of [16], are usually more practical and versatile.
Indeed, as the diffusion X is not controlled, this optimal switching problem can be seen as an extended American option problem. This suggests that, up to some adjustments, the probabilistic numerical tools developed in this context (see [13] for instance) may be adapted to solve (2.1).
To be more specific, consider a variantv of (2.1) such that the switching decisions can only take place on a finite time grid Π = {t 0 = 0 < t 1 < . . . < t m = T } for a fixed T > 0. Then ∀i ∈ I q , ∀x ∈ R d , and ∀t k ∈ Π, the dynamic programming principle (2.7) becomes: where: and where the notation X t,x := (X t,x s ) s≥t refers to the process X conditioned on the initial value X t = x.
If, moreover, the cost function k is such that at most one switch can occur on a given date t k (triangular condition), then equation (2.8) can be simplified into: which is explicit in the sense thatv . (t k , .) directly depends onv . (t k+1 , .).
In practice, apart from the potential approximation of the stochastic process X and of the final values (2.9), the difficulty lies in the efficient computation of the conditional expectations (2.10).
In the American option literature, various approaches have been developed to solve (2.11) efficiently. Notable examples are the least-squares' approach ( [43,57]), the quantization approach and the Malliavin calculus based formulation (see [13] for a thorough comparison and improvements of these techniques). In the spirit of [17], one may also consider non-parametric regression (see [38] and [56]) combined with speeding up techniques like Kd-trees ( [32,40]) or the Fast Gauss Transform ( [61,47,50,54,51]) in the case of kernel regression.
Here, we intend to solve (2.1) on numerical applications which bears the particularity of handling stochastic processes in high dimension (dim (X) = d 3, with however dim (I) = d ≈ 3, see Section 5). For such problems, the most adequate technique so far seems to be the local regression method developed in [13]. We are thus going to make use of this specific method to solve (2.11) in practice.
In the following, we provide a detailed analysis of the above suggested computational method.

Numerical approximation and convergence analysis
This section is devoted to the precise description of the resolution of (2.1), along the lines of the discussions from Subsection 2.3. Moreover, the convergence rate of the proposed algorithm will be precisely assessed.

Approximations
We are going to consider the following sequence of approximations: The time horizon will be truncated to a finite horizon T .
• [Time discretization] The continuous state process X and investment process I will be discretized with a time step h.
• [Space localization] The R d -valued process X will be projected into a bounded domain D ε , parameterized by ε.
• [Conditional expectation approximation] The conditional expectation involved in the dynamic programming equation will be replaced by an empirical least-squares regression, computed on a bundle of M Monte Carlo trajectories, on a finite basis of local hypercubes with edges of size δ.
The rate of convergence of the algorithm will then be provided, as a function of these five numerical parameters: T , h, ε, M and δ.

Finite time horizon
The first step is to reduce the set of strategies to a finite horizon: where 0 ≤ t ≤ T < +∞, and A T t,i ⊂ A t,i is the subset of strategies without switches strictly after time T . Hence the final value g f corresponds to the remaining gain after T .
Alternatively, one may choose, for convenience, another final value g instead of g f , as long as it is Lipschitz-continuous and satisfies a suitable condition (cf. equation (3.21)). The set of such functions will be denoted as Θ g f . The difference between the two value functions is quantified in Proposition 3.1.
This freedom on the final values will be used in practice to avoid a computation on an infinite interval [T, ∞[ as in the definition of g f . From now on, we choose and fix one such g ∈ Θ g f .

Time discretization
Then, we discretize the time segment [0, T ]. Introduce a time grid Π = {t 0 = 0 < t 1 < . . . < t N = T } with constant mesh h. Consider the following approximation: i is the subset of strategies such that switches can only occur at dates τ n ∈ Π ∩ [t, T ]. Now, with a slight abuse of notation, we can safely switch from the notation α = (τ n , ζ n ) n≥0 to the notation α = (τ n , ι n ) n≥0 (remember Subsection 2.1), replacing the quantity t≤τn≤T k (τ n , ζ n ) by t≤τn≤T k τ n , I α τn−h , I α τn or by t≤τn≤T k (τ n , ι n−1 , ι n ), where k (t, i, j) = k(t, j − i). The error between v T and v Π is quantified in Proposition 3.2.
Next we also approximate the stochastic process X by its Euler schemeX = X s 0≤s≤T , with dynamics: The error between v Π andv Π is computed in Proposition 3.3.

Space localization
Denote this projection asX ε t :X In other words,X ε t is equal toX t most of the time (i.e. whenX t ∈ D ε t ), except whenX t is outside D ε t , in which caseX ε t corresponds to the projection ofX t onto D ε t . Definev ε Π as the value functionv Π from equation (3.6) withX replaced byX ε . The error between those two value functions is computed in Proposition 3.4.

Example 3.1.
To clarify this construction of space localization, we explicit it on the very simple example of a d-dimensional standard brownian motion Choose D ε t to be a centered, symmetric hypercube: where W 1 is a one-dimensional Brownian motion. Hence, finding a value for C (t, ε) such that (3.7) holds boils down to inverting Bachelier's option pricing formula in order to get the strike as a function of the price of the call option. This is done in [6], see [52], but under the form of a series expansion for small moneyness, which is unsuitable for our purpose (because C (t, ε) → ∞ when ε → 0). Thus, we are here only going to look for a simply invertible upper bound for (3.8). Denoting as N the cumulative distribution function of a standard Gaussian random variable, and using the standard Inverting this last upper bound, the inequality (3.7) is satisfied with

Conditional expectation approximation
From now on, in order to prevent the notation from becoming too cumbersome and clumsy, we are going to drop the ε index in the following, i.e.X t will stand forX ε t , andv Π forv ε Π . For the fully discretized problem (3.6), the dynamic programming principle (2.11) becomes: The last step is to approximate the conditional expectation appearing in equation (3.9). As discussed in Subsection 2.3, we choose to approximate it using the following regression procedure. Consider basis functions (e k (x)) 1≤k≤K , As truncating the approximated conditional expectations is a necessity in theory as well as in practice (see [12,30,55]), suppose that there exist known bounds Γ tn,x (ϕ) and Γ tn, Then, ∀i ∈ I q the quantity E ϕ t n+1 ,X tn,x tn+1 , i is approximated by: which is used to define the next approximationṽ Π of the value function: Interesting discussions on the choice of function basis can be found in [13]. In particular they advocate bases of local polynomials, which is numerically efficient and well-suited to tackle large-dimensional problems (see Subsection 4.1). However, for the sake of simplicity, we will restrict our study in this section to a basis of indicator functions on local hypercubes (as in [55] and [30]) (which is the simplest example of local polynomials) as defined below: For every t n ∈ Π, consider a partition of the domain D ε tn into hypercubes B k tn k=1,...,Kε , i.e., ∪ k=1,...,Kε B k tn = D ε tn and B i tn ∩ B j tn = ∅ ∀i = j. It may be deterministic, or F t -measurable. We only assume that there exists (δ, δ) ∈ R 2 + with δ ≤ δ such that the lengths of the edges of the hypercubes, in each dimension, belong to [δ, δ] (in particular, the volume of each hypercube B k tn belongs to δ d , δ d ). This liberty over the definition of the partition enables to encompass the kind of adaptative partition described in [13]. Then, the basis functions considered here are defined by e k tn (x) : With this choice of function basis, the error betweenv Π andṽ Π is computed in Proposition 3.5.
Finally, let X m tn 1≤m≤M 1≤n≤N be a finite sample of size M of paths of the processX. The final step is to replace the regression (3.10) by a regression on this sample: (3.14) Then ∀i ∈ I q the quantity E ϕ t n+1 ,X tn,x tn+1 , i is approximated by: leading to the final, computable approximationv Π of the value function: The error betweenṽ Π andv Π with the same choice of function basis is given in Proposition 3.6. This proposition will make use of the following quantity: which is strictly positive, as the domains D ε t , t ∈ [0, T ] are bounded. More precisely, only lower bounds of these quantities will be required.

Example 3.2.
Carrying on with Example 3.1 of a d-dimensional Brownian motion, we explicit a lower bound for p (t n , δ, ε) in this simple case. First, Remark however that this lower bound is very crude, and that it can be very far below p (t n , δ, ε) for large δ.
Combining all these results, we obtain a rate of convergence ofv Π towards v: Remark 3.2. The adaptative local basis can be such that each hypercube contains approximately the same number of Monte Carlo trajectories (see [13]). This means that is the number of functions in the regression basis. With this remark in mind, the leading error term in Theorem 3.1 behaves like This is close to the corresponding statistical error term in [42] hM ) in the context of BSDEs. The advantage of their approach is that they can handle any (orthonormal) regression basis, while our approach (in the context of optimal switching) provides a bound on the L p error for every p ≥ 1.

Example 3.3.
In the case of a d-dimensional Brownian motion, the rate of convergence of Theorem 3.1 can be explicited further, using the upper bound on C (T, ε) from Example 3.1 and the lower bound on p (T, δ, ε) from Example 3.2. Moreover, one can express the rate of convergence as a function of only one parameter, choosing the five numerical parameters T , h, ε, δ and M accordingly. For instance, assuming δ = δ, and minimizing over δ, h, ε and T , one can get a convergence rate upper bounded by C p (1 + |x|) . This is admittedly highly demanding in terms of sample size M , but remember that this expression suffers from the crude lower bound on p (T, δ, ε) we chose previously.

Convergence analysis
From now on, we suppose that all the assumptions from Subsection 2.2 are in force.

Finite time horizon
Proof. First, we introduce the following notations: for any admissible strategy α ∈ A t,i . In particular: which provides the first inequality. Consider now the second inequality. Choose In order not to mix up the variables τ n and ζ n from different strategies, we add the name of the strategy in index when needed. Then: as k (s, 0) = 0 and k ≥ 0 (Assumption 4). Hence, using Jensen's inequality and equation (2.4), ∃C > 0 such that Since this is true for any ε > 0, and that C, ρ and ρ 1 do not depend on ε, the proposition is proved. Now, we focus on the final boundary g f . For the time being, denote the value function (3.2) as v g f T to emphasize the dependence of v on the terminal condition. As a consequence of equation (2.4), as v g T the value function defined as in equation (3.2) with g instead of g f . We are going to show that the precise approximation error due to the choice of final value g does not matter much as long as g is chosen in this class Θ g f .
To shorten the proof, we assume that v g f T (resp. v g T ) admits an optimal strategy α * f ∈ A T t,i (resp. α * ∈ A T t,i ) (this assumption can then be relaxed using ε-optimal strategies as in the proof of Proposition 3.1) 1 . Therefore, recalling the notations H (equation (3.18)) and J (equation (3.19)) introduced in the proof of Lemma 3. 1: Symmetrically, the same inequality holds for v g x, i), ending the proof.
From now on, we choose and keep one final value function g ∈ Θ g f , and remove the index g from the notation of v and its subsequent approximations. 1 Note that under the assumptions from Subsection 2.2, one may use Theorem 3.1 from [35] to get the existence of a unique optimal strategy α * for the value function (3.2), satisfying E 0≤τ α * n ≤T k τ α * n , ζ α * n 2 < ∞

Time Discretization
Proposition 3.2. There exists a positive constant C such that for any (t, Proof. Under the assumptions from Subsection 2.2, one can apply Theorem 3.1 in [29] to prove (3.24), noticing that the cost function k does not depend on the state variable x.
Use the discounting factor in the definition of f to factor the e −ρt term and to get that C does not depend on T .
Remark 3.3. Another alternative to get this rate of h 1 2 is to work with the reflected BSDE representation of v Π , as in [16] (adapting [12]) or [21]. , as stated in [29] (making use of results from [27]).

Proposition 3.3. There exists
Proof. T and g being fixed, we can define, in the spirit of equations (3.18) and (3.19), the following quantities: for any admissible strategy α ∈ A Π t,i . For these discretized problems, the existence of optimal controls α * andᾱ * is granted. Hence: using the strong convergence speed of the Euler scheme on [t, T ]. Symmetrically, the same inequality holds forv Π (t, x, i) − v Π (t, x, i), ending the proof.

Conditional expectation approximation
From now on the domains D ε t , t ∈ [0, T ] are fixed once and for all, and, with a slight abuse of notation, we will drop ε from the subsequent notations.

Lemma 3.6. Consider a measurable function
Proof. Define the following centered random variables: Then: Now, for any p ≥ 1: and: using Markov's inequality. Now, the following lemma will provide upper bounds for E ε tn,X Proof. Using Marcinkiewicz-Zygmund's inequality, there exists C p > 0 such that: Now, suppose that ∃φ tn ∈ R + s.t. ϕ t n+1 ,X tn+1 , j ≤φ tn a.s. . Then, using Lemma 3.7, ∃C p > 0 such that: (3.50) where, for the second inequality, the term m = l in the sum was treated separately. Then: In a similar manner: We now apply Lemma 3.6 tov Π in the following Corollary: Proof. First, recall from equation (3.36) and (3.37) that there exists C > 0 such that for every (t n , j) ∈ Π × I q : Hence one can apply Lemma 3.6 toṽ Π with these upper bounds. The final step is to recall that the minimum probability p (T, δ, ε) defined in equation (3.17) is a lower bound on P X tn ∈ B tn X l tn for any (t n , l) Using this result, we can now assess the error betweenṽ Π andv Π .
where I tn q is the set of F tn -measurable random variables taking values in I q .
Proof. For each t n ∈ Π, we look for an upper bound E n , independent of l, such that:

Computational complexity
The number of operations required by the algorithm described below is in O q 2 .N.M , where we recall that q is the number of possible switches, N is the number of time steps and M is the number of Monte Carlo trajectories.
• The q 2 term stems from the fact that for every i ∈ I q , one has to compute a maximum on j ∈ I q (see equation (3.16)). However, this q 2 can be reduced to q as soon as the two following conditions are satisfied: 1. (Irreversibility) The controlled variable can only be increased (or, symmetrically, can only be decreased)

(Cost Separability) There exists two functions
For instance, this is true of affine costs.
Indeed, under those two conditions, equation (3.16) becomes: These maxima can be computed in O(q) instead of O q 2 by starting from the biggest element i = i q down to the smallest element i = i 1 (in lexicographical order) and keeping track of the partial maxima.
Note that these two conditions hold for the numerical application from Section 5, providing the improved complexity O(q.N.M ).
• The N term comes from the backward time induction.
• The M term corresponds to the cost of a regression, which, in the case of a local basis, can be brought down to O (M ) (cf. [13]).

Memory complexity
The memory size required for solving optimal switching problems (as well as the simpler American option problems and the more complex BSDE problems) by Monte Carlo methods is often said to be in O(N.M ), because, as the Euler scheme is a forward scheme and the dynamic programming principle is a backward scheme, the storage of the Monte Carlo trajectories seems inescapable. This fact is the major limitation of such methods, as acknowledged in [16] for instance.
Since such a complexity would be unbearable in high dimension, we describe below a general memory reduction method to obtain a much more amenable O(N + M ) complexity (or, more precisely, of O(m.N + q.M ) with m M ). This improvement really opens the door to the use of Monte Carlo methods for American options, optimal switching and BSDEs on high-dimensional practical applications. Note that this tool can be combined with all the existing Monte Carlo backward methods which (seem to) require the storage of all the trajectories.
A drawback of this tool is that it is limited to Markovian processes. However, one can usually circumvent this restriction by increasing the dimension of the state variable.

Description
The memory reduction method for Monte Carlo pricing of American options was pioneered by [18] for the geometric Brownian motion, and was subsequently extended to multi-dimensional geometric Brownian motions ( [19]) as well as exponential Lévy processes ( [20]). These papers take advantage of the additivity property of the processes considered. However, as briefly hinted in [59], the memory reduction trick can be extended to more general processes. In particular, it can be combined with any discretization scheme, for instance the Euler scheme or Milstein scheme, as long as the value of the stochastic process at one time step can be expressed via its value at the subsequent time step.
From a practical point of view, the production of "random" sequences usually involves wisely chosen deterministic sequences, with statistical properties as close as possible to true randomness (cf. [39] for instance for an overview). These sequences can usually be set using a seed, i.e. a (possibly multidimensional) fixed value aimed at initializing the algorithm which produces the sequence: where the rand() operation consists in going to the next element of the sequence. Now two useful aspects can be stressed. The first is that one can usually recover the current seed at any stage of the sequence. The second is that, if the seed is set later to, say, once again the seed s from equation (4.1), then the following elements of the sequence will be once again ε 1 , ε 2 , . . . In other words, one can recover any previously produced subsequence of the sequence (ε n ) n≥1 , provided one stored beforehand the seed at the beginning of the subsequence. This feature is at the core of the memory reduction method, which we are going to discuss below in a general setting.
Consider a Markovian stochastic process (X t ) t≥0 , for instance the solution of the stochastic differential equation (2.2), recalled below: The application of the Euler scheme to this equation can be denoted as follows: where ∀i ∈ [0, N − 1] and ∀j ∈ [1, M ], ε j i ∈ R d is drawn from a d-dimensional Gaussian random variable. Suppose that for any ε ∈ R d , the function x → f (x, ε) is invertible (call f inv its inverse). Then, starting from the final value x j t N of the sequence (4.2), one can recover the whole trajectory of X: as long as one can recover the previous draws ε j N −1 , . . ., ε j 0 . The following pseudo-code describes an easy way to do it. In other words, at the price of doubling the computation time, one can bring down the required memory storage by the factor min (M, N ), which is a very significant saving. Moreover, the theoretical additional computation time can be insignificant in practice, as the availability of much more physical memory makes the resort to the slower virtual memory much less likely.
Remark 4.1. Even though the storage of the seeds does take O (N ) in memory size, the constant may be much greater than 1. For instance, on Matlab ® , a seed from the Combined Multiple Recursive algorithm (refer for instance to [39] for a description of several random generators) is made of 12 uint32 (32-bit unsigned integer), a seed from the Multiplicative Lagged Fibonacci algorithm is made of 130 uint64, and a seed from the popular Mersenne Twister algorithm is made of 625 uint32.
In order to relieve the storage of the seeds, we now provide a finer memory reduction algorithm (Algorithm 2). Although Algorithm 2 requires three main loops, it enables to perform the last loop without fiddling the seed of the random generator, and without any vector of seeds locked in memory, which will thus be fully dedicated to the regressions and other resolution operations. Moreover, the first two main loops can be performed beforehand once and for all, storing only the last values of the vector X as well as the first seed S [0]. Finally, if the random generator is able to leapfrop a given number of steps, the first loop can be drastically reduced.

Numerical stability
Theoretically, the trajectories produced by the Euler scheme (4.2) and the inverse Euler scheme (4.4) are exactly the same. In practice however, a discrepancy may appear, the cause of which is discussed below.
On a computer, not all real numbers can be reproduced. Indeed, they must be stored on a finite number of bits, using a predefined format (usually the IEEE Standard for Floating-Point Arithmetic (IEEE 754)). In particular, there exists an incompressible distance ε > 0 between two different numbers stored. This causes rounding errors when performing operations on real numbers.
For instance, consider x ∈ R and an invertible function f : R → R. Compute y = f (x) and then computex = f inv (y). One would expect thatx = x, but in practice, because of rounding effects, one may getx = x + z for a small > 0, where z is a discrete variable, which can be deemed random, taking values around zero. This phenomenon is illustrated on  x j ti+1 = f x j ti , ε j i Now, instead of equation (4.4), the inverse Euler scheme will provide something like: for a small > 0, where z j i , i = 0, . . . , N , j = 1, . . . , M , can be deemed realizations of a discrete random variable Z, independent of W . The distribution of Z is unknown, but data suggests it may be innocuously assumed centered, symmetric, and with finite moments.
We are now interested in studying the compound rounding error y ti − x ti as a function of . Of course, its behaviour depends on the choice of f (equation (4.3)). Below, we explicit this error on two simple examples: an arithmetic Brownian motion and an Ornstein-Uhlenbeck process. These two examples illustrate how the compound rounding error can vary dramatically w.r.t. f .

First example: arithmetic Brownian motion
Consider first the case of an arithmetic Brownian motion with drift parameter µ and volatility parameter σ. Here f and its inverse are given by: Hence, using equation (4.5), for every j = 1, . . . , M : In other words, the compound rounding error behaves as a random walk, multiplied by the small parameter . Hence, as long as h (which is always the case as real numbers smaller than cannot be handled properly on a computer), this numerical error is harmless.
Remark that a similar numerical error arises from the algorithms proposed in [18] , [19] and [20], but, fortunately, as discussed above, this error is utterly negligible.
Second example: Ornstein-Uhlenbeck process Now, consider the case of an Ornstein-Uhlenbeck process with mean reversion α > 0, long-term mean µ and volatility σ. Here: Using equation (4.5), for every j = 1, . . . , M the compound error is given by: when h → 0, one can see that, as soon as T > − ln( ) α , this error may become overwhelming. This phenomenon is illustrated on Figure 4.2a on a sample of 100 trajectories.
In order to mitigate this effect, we propose to modify the Algorithm 2 as follows: in its second loop (usual Euler scheme), instead of saving only the last values x j T , one may define a small subsetΠ ⊂ Π and save the intermediate values x j ti , t i ∈Π. Then, in the last loop (inverse Euler scheme), every time that t i ∈Π, the current value of the set x j ti may be recovered from this previous storage. The drawback of this modification, of course, is that it multiplies the required storage space by the factor #Π. However, this remains much smaller than the O (M × N ) required by the naive full storage algorithm.

Application to investment in electricity generation
This section is devoted to an application of the resolution method studied in Section 2. We choose to apply it to an investment problem in electricity generation on a single geographical zone. We intend to show that it is possible to provide a probabilistic outlook of future electricity generation mixes instead of a deterministic outlook provided by planification methods. Nevertheless, the problem presents so many difficulties that addressing all of them in the same model is unresonable. Some aspects have thus to be left aside. Our goal here is to show that the algorithm described in Section 3 can handle high-dimensional investment problems. We focus on the influence of investment decisions on the spot price, consistently with the fundamentals of the electricity spot price formation mechanism.
Although the strategic aspect of investment is an important driver of utilities' decisions, this aspect is beyond the scope of our modeling approach. There exist models limited to a two-stage decision making (see for instance [48]), but in the case of continuously repeated multiplayer game models, defining what is a closed-loop strategy is already a difficulty (see Sec. 2 in [7]).
We did not consider time-to-build in this implementation either. Relying on the fact that it is possible to transform an investment model with time-to-build into a model without time-to-build by replacing capacities with committed capacities (see [9,1] for implementations in dimension one, and [26] in dimension two), we left this aspect for future work.
Finally, we did not consider the dynamic constraints of power generation. Their effect on spot prices is well-known: they tend to increase spot prices during peak hours and to decrease them during off-peak hours (see [41]). However, we assume here that this effect is negligible compared to the effect induced by a lack or an excess of capacity.
Thus, we focused on the following key factors of electricity spot prices: demand, capacities (including random outages) and fuel prices. Our model is based on [5,4], where the electricity spot price is defined as a linear combination of fuel prices multiplied by a scarcity factor. This model exhibits the main feature wanted here, which is that the spot price, being determined both by the fuel prices and the residual capacity, is directly affected by the evolution of the installed capacity. When the residual capacity tends to decrease, spot prices will tend to increase, making investment valuable. Thus, in this model, investments are undertaken not on the specific purpose of satisfying the demand but as soon as they are profitable. In our example, new capacities are invested according to the criterion of value maximization. Energy non-served and loss of load probability may still be adjusted through the price cap on the spot market.
In this section, we first detail the chosen modelling and objective function (which will be shown to be encompassed in the general optimal multiple switching problem (2.1)), and then solve it numerically using the general algorithm developed in the previous sections.

Modelling
The key variable in order to describe our electricity generation investment problem is the price of electricity. More precisely, the key quantities are the spreads between the prices of electricity and other energies. To model these spreads accurately, it may be worth considering a structural model for electricity (cf. the survey [14]). Here we choose such a model, mainly inspired by those introduced in [5] and [4], albeit amended and customized for a long-term time horizon. All the variables involved are detailed below.

Electricity demand
and its main simulation method (the truncated Euler scheme) disables our memory reduction method described in Subsection 4.2. Hence we look for a simpler model.
In [60], a detailed structural model for electricity is developed, which includes renewable energies like wind and solar. In particular, wind power infeed efficiency (which belongs to [0, 1]) is modelled as a logit transform of an Ornstein-Uhlenbeck process with seasonality. Adapting this idea, we model

Electricity price
where κ f + Consider the case of new plants (ζ j > 0). Assuming that the global availability rate(5.6) of technology j applies to the new plants, they can then produce up to ζ j A j s , s ≥ t, or, more precisely, according to the stack order principle: assuming that, in the stack order, the new plants are called before the older plants I t− of the same technology (as they can be expected to have an at least slightly better efficiency rate compared to the older plants of the same technology, a phenomenon that can be seen as partly captured by the function (5.9)).
At time s ≥ t, this production is sold at price P s , but costsS s to produce (if P s <S s , then of course the producer chooses not to produce). In addition, regardless of the output level, there may exist a fixed maintenance cost κ j . Summing up all these gains, discounted to time t using a constant interest rate ρ > 0, the new plant yield a revenue of: . This was the cost-benefit analysis for one quantity ζ j of new plants. Now, consider as a whole the full fleet of the geographical zone considered. Maximizing the expected gains along the potential new plants yields the following value function: where the strategies α affect the installed capacities (equations (5.4)), hence also the available capacities (equation (5.5)) as well as the power price (equation (5.10)), and where the cash flows are purposely discounted up to time 0, the time of interest.
Remark 5.2. Replacing P in (5.11) by its definition (5.10), it is patent that this objective function fits into the mould studied thoroughly in Section 3. In Subsection 5.2 below, our algorithm will be applied to this specific objective function.
Remark 5.3. Remark that under this modelling, the demand is satisfied as long as it does not exceed the total available capacity. Indeed, the effective output of the plant ζ j is equal to

Numerical results
Finally, we solve the control problem described in Subsection 5.1 on a numerical example, using the algorithm detailed in Subsection 3 combined with the general memory reduction method described in Subsection 4.2.
Our purpose here is not to perform a full study of investments in electricity markets, but a more modest attempt at illustrating the practical feasibility of our approach, with some possible outputs that the algorithm can provide.
We consider a numerical example including two cointegrated fuels (in addition to the price of CO 2 ): one "base fuel" and one "peak fuel", starting respectively from 40e/MWh and 80e/MWh. Hence, using the notations from Subsection 5.1, d = 2 (two technologies) and d = 6 ( electricity demand, CO 2 price, two fuel prices and two availability rates). The main choices of parameters for this application (initial fuel prices and volatilities, initial fleet and proportional costs of new power plants) are summed up in Table 5.1. Moreover, the demand process starts from D 0 = 70GW and does not integrate any linear trend.
κ p+ i 1 40e/MWh 5% 67GW 0.24 10 9 e/GW 2 80e/MWh 15% 33GW 2.00 10 9 e/GW In order to take into account the minimum size of one power plant we restrict the values of the installed capacity process(5.4) to a (bi-dimensional) fixed grid Λ d , with a mesh of 1GW. We make the simplifying assumptions that investments are irreversible, and that no dismantling can occur (recall from Subsection 4.1 the computational gain provided by this assumption).
Remark 5.4. If such a grid is indeed manageable in dimension d = 2, it may less be the case if additional technologies were considered. However, as discussed in [55] equation (3.2), instead of performing one regression for each i ∈ Λ d , one can solve equation (3.16) at time t i by only one (d + d )-dimensional regression, by choosing an a priori law for the randomized control ζ ti . The error analysis from Section 2 can be easily generalized to such regressions in higher dimension.
Finally, we consider the following numerical parameters. We choose a time horizon T = 40 years and a time step h = 1 730 (i.e. two time steps per day, allowing for some intraday pattern in the demand process) but allow for only one investment decision per year. For the regression, we consider a basis of b = 2 d = 64 adaptative local functions, chosen piecewise linear on each hypercube (which is a bit more refined than the piecewise constant basis studied in Section 3) on a sample of M = 5000 trajectories.  Figure 5.1a displays the time evolution of the average as well as the variability of the optimal fleet (only the new plants are shown). One can distinguish a first short phase characterised by the construction of several GW of peak load assets, followed by a much slower second phase involving the construction of both base load and peak load assets. Moreover, the variability of the optimal fleet increases over time. The detailed histogram of the optimal strategy at time T = 40 years is displayed on Figure 5.1b, where it is combined with the price of fuel. One can see that the more the peak fuel is expensive (and hence both fuels are expensive on average, as they are cointegrated), the more constructions of base load plants occur.  The fact that the average fleet seem to converge is related to the fact that this numerical example does not consider any growth trend in the electricity demand (see equation (5.2)). Otherwise, more investments would occur, indeed, to keep the pace with consumption.
Then Figure 5.2 provides information on the price of electricity. Figure 5.2a displays the time evolution of the electricity spot price density. For better readability, each density covers one whole year. One can see how the density moves away from the initial bimodal density (with prices clustering around the initial prices of the two fuels) towards a more diffuse density. Moreover, the downward effect of investments on prices can be noticed. This downward effect is even more visible on Figure  5.2b. It compares the effect on electricity prices of three different strategies: the optimal strategy, the optimal deterministic strategy (computed as the average of the optimal strategy), and the donothing strategy. For each strategy, the joint time-evolution of the yearly median price and the yearly interquartile range are drawn. As expected, prices tend to be higher and more scattered without any new plant. Nevertheless, on this specific example, the price distribution under the optimal deterministic strategy is close to that under the optimal strategy (only slightly more scattered).  These few pictures illustrate the kind on information that can be be extracted from the resolution of this control problem. Of course, as a by-product of the resolution, much more can be extracted and analyzed (distribution of income, CO 2 emissions, optimal exercise frontiers, etc) if needed.

Conclusion
In this paper, we presented a probabilistic method to solve optimal multiple switching problems. We showed on a realistic investment model for electricity generation that it can efficiently provide insight into the distribution of future generation mixes and electricity spot prices. We intend to develop this work in several directions in the future. First, we wish to take into account more generation technologies, most notably wind farms, nuclear production, as well as solar distributed production. These additions would raise the dimension of the problem from eight to fifteen. Yet another range of innovations in numerical methods will be necessary to overcome this increase in dimension. Second, we wish to take time-to-build into account. And last but not least, we wish to adapt the problem to a continuous-time multiplayer game and contribute to the quest for an efficient algorithm to solve it.