Forest of Stochastic Trees: A Method for Valuing Multiple Exercise Options

: We present the Forest of Stochastic Trees (FOST) method for pricing multiple exercise options by simulation. The proposed method uses stochastic trees in place of binomial trees in the Forest of Trees algorithm originally proposed to value swing options, hence extending that method to allow for a multi-dimensional underlying process. The FOST can also be viewed as extending the stochastic tree method for valuing (single exercise) American-style options to multiple exercise options. The proposed valuation method results in positively- and negatively-biased estimators for the true option value. We prove the sign of the estimator bias and show that these estimators are consistent for the true option value. This method is of particular use in cases where there is potentially a large number of assets underlying the contract and/or the underlying price process depends on multiple risk factors. Numerical results are presented to illustrate the method. computational times comprises 2.2 GHz Opteron processors with 4 GB per core and Inﬁni-Band interconnections. Run times are normalized to the run time of a single


Introduction
The broad class of stochastic optimal control problems includes many important applications in management sciences and quantitative finance such as development of natural resources, project initiation or abandonment, maintenance scheduling, land use decisions, and valuation and hedging of complex contracts. Very few problems have closed form solutions. For example the Black-Scholes-Merton formula for pricing European-style options does not have an American-style option analog. Approaches such as binomial lattice methods, partial differential equation (PDE) methods, variational inequalities and integral equations have been adopted for pricing these types of derivatives. However, all of these methods mentioned are limited in the number of sources of uncertainty and the dimensionality of the underlying asset that can be practically incorporated, leaving Monte Carlo (MC) as a general method of solution.
Multiple exercise options (MEOs) are generalizations of American-style options as they provide the holder more than one exercise right and sometimes control over one or more other variables, such as the amount exercised. Similar to pricing American-style options, MEO valuation is a problem in stochastic optimal control. For American-style options, the solution provides both a value and optimal exercise rule-a stopping time. For MEOs the solution gives both a value and optimal exercise policy. In cases in which the holder controls only the exercise times, the exercise policy is a sequence of stopping times. For MEOs in which the holder controls the exercise times and amounts, the exercise policy is a paired sequence of stopping times and exercise amounts. The policy generalizes to other control variables. Multiple exercise option valuation algorithms are generalizations of those used for pricing American-style options. Most of the work in the literature has focused on swing options in which the holder controls the exercise times and amounts. Hence, we discuss MEO valuation in the context of swing options and note that the FOST method and the FOST price estimators' properties apply to general MEOs.
The Forest of (Binomial/Trinomial) Trees method of Lari et al. (2001) and Jaillet et al. (2004) extends the binomial method of Cox et al. (1979) for pricing American style options to value swing options. As with most tree-based methods, the Forest of Trees is not able to handle high-dimensional underlying processes. Tangential to the Forest of Trees extension, the Stochastic Tree method of Broadie and Glasserman (1997) extends the binomial method for pricing American style options to allow for a high-dimensional underlying asset.
In this paper, we replace the binomial and trinomial trees of Lari et al. (2001) and Jaillet et al. (2004) with the stochastic trees of Broadie and Glasserman (1997), hence creating the Forest of Stochastic Trees (FOST) method for valuing multiple exercise options. The FOST can be thought of as generalizations of two different methodologies; specifically, it extends

•
The Forest of Trees method of Lari et al. (2001) and Jaillet et al. (2004) for valuing multiple exercise options on a single asset to allow for high-dimensional underlying assets and processes; and • The Stochastic Tree method of Broadie and Glasserman (1997) for valuing high-dimensional American-style options (single exercise right) to options having multiple exercise rights and additional controls (e.g., volume control).
Visually these generalizations complete the final entry in the two-by-two table shown in Table 1.

1-Dimensional Asset High-Dimensional Asset
Single exercise American option Binomial Tree Stochastic Tree Cox et al. (1979) Broadie and Glasserman (1997) Multiple exercise option Forest of Trees Forest of Stochastic Trees Lari et al. (2001), Jaillet et al. (2004) Marshall and Reesor We construct high and low biased FOST estimators that are analogous to those from the Stochastic Tree. The main theoretical results, presented in Sections 2.2 and 2.3 and proven Appendix B, are listed below.
The low FOST estimator has negative bias. 3.
On any given realization the high FOST estimator is at least as big as the low FOST estimator. 4.
The high and low FOST estimators are asymptotically unbiased.
The remainder of this section provides a literature review. Section 2 describes the Forest of Stochastic Trees estimators in detail; specifies results giving estimator properties (e.g., biasedness and convergence); provides numerical results showing the method is effective at pricing; and supplies discussion on improving computational performance through parallel processing. Section 3 summarizes and concludes the paper. A list of the notation used is given in Appendix A and proofs of the paper's main theoretical contributions on estimator properties are given in Appendix B.

Literature Review
Multiple exercise options arise in many different areas and the structure of these contracts is typically tailored to particular clients/needs, in contrast to standardized derivatives such as interest rate swaps and exchange-traded commodity futures. A non-exhaustive list of examples of MEOs include (i) tolling agreements used in the steel Kim et al. (2019) and electricity Deng and Oren (2006) sectors; (ii) chooser flexible caps which are exotic interest rate derivatives Meinshausen and Hambly (2004); (iii) valuation and control of energy production and storage facilities Chen and Forsyth (2007); Ludkovski and Carmona (2010); Thompson et al. (2009); and (iv) swing options Calvo-Garrido et al. (2017); Jaillet et al. (2004);Lari et al. (2001); Wilhelm and Winter (2008).
Valuation methods for MEOs are extensions of those used for American-style options. There are continuous-time solutions to both the American-style and multiple exercise option valuation problems; these are computed by solving a system of Hamilton-Jacobi-Bellman quasi-variational inequalities Korn et al. (2005). These methods give more accurate and stable price and sensitivity estimates than those computed using simpler tools (e.g., trees). However, these methods are quite complex mathematically and break down in higher dimensions.
In this article, we focus on the mathematically simpler time-discretized version of the valuation problem. Discrete-time tree-based methods for valuing American-style options Cox et al. (1979) have been extended to MEOs via the Forest of Trees Jaillet et al. (2004);Lari et al. (2001). Techniques for pricing American-style options using solutions of PDEs have been modified to MEOs Calvo-Garrido et al. (2017); Chen and Forsyth (2007); Thompson et al. (2009); Wilhelm and Winter (2008). These methods for MEOs inherit properties similar to the corresponding methods for single-exercise options. One crucial property is that these methods fail as the dimensionality of the problem increases.
Monte Carlo is the obvious tool to overcome the curse of dimensionality, as the rate of convergence of Monte Carlo estimators is independent of the dimension. Tilley (1993) was the first to show that the forward-in-time Monte Carlo approach could be used to solve the backward-in-time dynamic programming problem arising from valuation of an American-style option. Since this seminal paper, numerous other methods for the Monte Carlo valuation of American style options have appeared. These include methods that attempt to parameterize the exercise region Barraquand and Martineau (1995) and those that discretize the state space Bally et al. (2005). Methods that parameterize the early-exercise region have been extended to value multiple exercise options by parameterizing the set of exercise level curves Ibánez (1996). Similarly state space aggregation methods have been used for multiple exercise option valuation Ben Latifa et al. (2016). These approaches, however, also suffer from the curse of dimensionality and do not easily generalize to arbitrary payoffs and underlying price processes.
Monte Carlo methods that do not break down with the dimensionality and that accommodate general payoff and price processes include those that solve the optimal stopping-time problem through estimation of the hold or continuation value. These include the stochastic tree and mesh techniques of Glasserman (1997 2004) and the regression-based approach first appearing in Carriere (1996) and then subsequently generalized in Longstaff and Schwartz (2001). For each of these valuation techniques, high-and low-biased estimators are easily generated, along with a hybrid interleaving estimator that has properties of both. Duality-based methods solve the optimal control problem in the dual space by approximating an optimal martingale, typically by regression Andersen and Broadie (2004); Haugh and Kogan (2004).
Least-squares Monte Carlo has been modified for the pricing of swing options in Barrera et al. (2006); Meinshausen and Hambly (2004), respectively. Although increased dimensionality does not decrease the performance of these methods, they suffer from other drawbacks. In least-squares Monte Carlo methods one must select a set of basis functions on which to run regressions to estimate continuation values. In general only a complete (infinite) set of basis functions results in continuation value estimators that are consistent for the true option value. In practice, of course, a finite set of basis functions is used and introduces an approximation error. Continuation value estimators are consistent for the true approximation value and not the true option value Clement et al. (2001);Stentoft (2004).
Duality methods have been extended to MEOs Bender (2011); Chandramouli and Haugh (2012);Gyurko et al. (2015); Meinshausen and Hambly (2004). Duality methods rely on having a sub-optimal exercise policy that produces a low-biased estimate from which the solution to the dual problem can be approximated to yield a high-biased estimate. Typically regression-based methods are used to estimate the sub-optimal exercise policy Chandramouli and Haugh (2012); Gyurko et al. (2015); Meinshausen and Hambly (2004) implying the above noted issues of least-squares Monte Carlo persist when pricing MEOs. Policy iteration methods such as Bender (2011), yield approximations of the time-0 value at each iteration of the dynamic program. As with the pricing of American-style options this method is advantageous because it removes the requirement to calculate nested conditional expectations prior to the time-0 value being approximated.
The stochastic mesh of Broadie and Glasserman (2004) has been extended to MEOs via the Forest of Stochastic Meshes (FOSM) Marshall (2012); . High and low biased FOSM estimators are derived  and their properties shown Marshall (2012), similar to the work presented here for the Forest of Stochastic Trees estimators.

Results
We consider the valuation of multiple exercise options as a stochastic optimal control problem with three relevant state variables-the underlying variable (S), number of exercise rights remaining (N ), and usage level (U) assuming some volume control. At each exercise opportunity and given (S, N , U), the current values of the state variables, the holder must choose between

•
Exercising u units plus continuing with an option having N − 1 remaining exercise rights and usage level U + u; and • Continuing with an option having N exercise rights and usage level U (i.e., no exercise).
Note that with volume control the payoff from exercising u units changes with u (as does the continuation value of the option). Thus, the holder chooses the value-maximizing u when deciding to exercise. Also note that with N = 1 and u constrained to be 1, this is an American-style option.
We work with the time-discretized problem and use dynamic programming to solve for the optimal exercise policy and the corresponding optimal value. In all variables, let the subscript i denote time-t i and let U i be the time-t i set of admissible volume choices which includes the zero volume choice (i.e., hold). The recursive equations for the dynamic program are (1) with the terminal conditions where H i (S, N , U) and B i (S, N , U) are the time-t i , state-Z i continuation and option values, respectively, h i (S, N , U, u) is the payoff from exercising u units with h i (S, N , U, 0) = 0, Z i is the time-t i information set generated by the paths of (S, N , U), I is an indicator function andφ (·) is a cumulative usage penalty term. Estimator properties and their proofs are given for this multiple exercise option setup. However, the dynamic program and estimator properties can be stated and proven for alternative specifications provided there is a finite number of exercise rights and usage levels. For example, a swing option contract may specify a certain number of up and down swing rights, N u and N d . An up swing right allows the holder to take more than the baseline amount of the underlying asset while a down swing right allows the holder to take less. Another variation is to allow for multiple rights to be exercised at each opportunity where each right corresponds to a fixed volume amount Bender and Schoenmakers (2006); Meinshausen and Hambly (2004).

Forest of Stochastic Trees
The FOST generalizes the stochastic tree method for valuing American-style options to the valuation of multiple exercise options and extends the Forest of Trees method to handle a high-dimensional underlying asset. This is done by replacing the binomial/trinomial trees with stochastic trees in the framework of Lari et al. (2001) and Jaillet et al. (2004) hence giving the FOST. The stochastic tree is constructed identically as described in Broadie and Glasserman (1997) and the tree is replicated multiple times, with one replication corresponding to each possible (N , U) combination. This is analogous to the Forest of Trees in which the same underlying binomial/trinomial tree is replicated for each possible (N , U) combination.
The dynamic program is approximately solved by replacing the continuation values in Equations (1) and (3) with stochastic tree-type estimators. As with the original stochastic tree technique, high-and low-biased option value estimators are constructed by using the analogous high-and low-biased estimators, respectively, on each stochastic tree in the forest. The recursive equations for the high estimator arê with the terminal conditionŝ whereĤ i (S, N , U) andV i (S, N , U) are the time-t i , state-Z i continuation and option value estimators, respectively, h i (S, N , U, u) (with h i (S, N , U, 0) = 0) is the time-t i , state-Z i payoff from exercising u units, b is the branching factor, I is an indicator function andφ (U m + u) is a global usage penalty term. The superscript j = {j 0 , j 1 , . . . , j i } indicates the specific node within a given stochastic tree and k = {j, k}. Figure 1 is a diagram of a section of a Forest of Stochastic Trees with two volume choices, u 1 and u 2 . It illustrates the nodes in the forest which need to be considered when making an exercise decision given state (N , U). The three choices are no exercise, exercise u 1 units, and exercise u 2 units. The low estimator is similarly defined using the low estimator on each stochastic tree via the dynamic program, whereĤ il (X j i , N i , U i ) is the l−th leave-one-out hold value estimator andû * is the estimated optimal exercise amount which depends on i and l. The terminal conditions associated with this dynamic programming scheme are, whereφ (U m + u) is a cumulative usage penalty term.

Estimator Bias
In order to justify using the high-and low-biased estimators to construct upper and lower option price bounds, respectively, we prove that the high estimator is always positively biased and that the low estimator is always negatively biased. In addition we include a comparison of the estimators which orders their values on any realization of the simulated forest.
The theorems that follow are direct extensions of those in Broadie and Glasserman (1997). Below, the branching factor, b, appears as an argument in the estimators. For example,V 0 (b, S 0 , N 0 , U 0 ) refers to the time-0, state-Z 0 high-biased estimator with a stochastic tree branching factor of b.
This argument has been suppressed to this point for convenience. We begin with the theorem regarding the bias of the high estimator.

Theorem 1. (High estimator bias)
The high estimator is biased high, i.e., Similarly, the result stating the bias of the low estimator follows.
Theorem 2. (Low estimator bias) The low estimator is biased low, i.e., for all b.
Finally, an ordering result for the high and low estimators is stated in Theorem 3.

Theorem 3. (Comparison of Estimators)
On every realization of the forest the low estimator is less than or equal to the high estimator. That is,v with probability one for all j, i.
Theorems 1-3 are proven in Appendix B.

Estimator Convergence
An advantage of the stochastic tree method over some other MC valuation methods is that its estimators are consistent to the true option value. This property also holds for the FOST estimators. In this section we state two theorems-one for the consistency of the high estimator, and the other for the consistency of the low estimator. Here convergence is in probability to the true option value and as above the argument b that appears with the estimators refers to an arbitrary branching factor size of b with convergence being shown as b → ∞. Before stating the result, defineV 0 (b, S 0 , N 0 , U 0 ) as the average of R independent replications ofV 0 (b, S 0 , N 0 , U 0 ).

Theorem 4. (High estimator convergence) Suppose
This holds for an arbitrary number of repeated valuations of the forest, R. In particularV 0 (b, S 0 , N 0 , U 0 ) converges to B 0 (S 0 , N 0 , U 0 ) in probability and is thus a consistent estimator of the option value.
This results implies that as b → ∞. Hence the estimator is asymptotically unbiased.
Theorem 5. (Low estimator convergence) Suppose that, for u 1 , u 2 ∈ U i , u 1 = u 2 and all i. Then Theorem 4 also holds for the low estimator.
The additional condition imposed in Theorem 5 is analogous to that used in Theorem 3 of Broadie and Glasserman (1997). This condition says that, with probability one, the optimal exercise policy is never indifferent between the choices of volumes to exercise (including u = 0). As in Broadie and Glasserman (1997) imposing this condition simplifies the analysis of the estimator. Theorems 4 and 5 are proven in Appendix B.

Numerical Results
Pricing swing contracts involves a large number of parameters and in this section we provide some results which illustrate the validity of our method across a variety of specifications. We assume that the underlying assets follow a risk neutral stochastic process, there are no transaction costs and other than penalties, there are no other constraints considered. We also assume a constant risk free rate of interest and the volatilities of all assets are known constant functions of time.
The option swing rights may be exercised at discrete times up to and including expiry and the volume choices given are in discrete amounts. That is, anytime the holder chooses to exercise a right, they must choose from a finite list of possible volume amounts. The rationale behind allowing all time steps to be exercise opportunities is the exponential growth in computational time caused by adding intermediate non-exercise times. However, the method can easily be modified to incorporate these extra time steps. As previously mentioned penalties can be implemented globally and are based on the net volume exercised during the contract.

Single Dimension
Beginning with the one dimensional case, we have based our simulations on an underlying asset with a risk neutralized price process that satisfies the following stochastic differential equation, In this equation, r is the riskless interest rate, Z i is a standard Brownian motion process, σ is a constant volatility parameter and the underlying asset itself pays a continuous dividend yield δ.
The parameter values for the underlying asset are specified as r = 0.05, δ = 0.1, and σ = 0.2. The swing options considered have both up and down swing rights and the payoff upon exercise is where u is the volume exercised, S is the price of the underlying asset at the exercise time, and K u and K d are the up and down swing strike prices, respectively. For the examples considered here, we set K u = K d = K which simplifies the payoff function to For all examples in Section 2.4.1 the option expiry is 3.0 years and the options have both up and down swing rights with strike prices K u = K d = 40.0, respectively. In examples where the holder controls the amount exercised, a list of volume choices is given. The volume choices are consecutive integer multiples of a base amount and the up swing and down swing volumes have the same magnitude. For comparison purposes the results in this subsection include a binomial value which is calculated using the Forest of Trees Lari et al. (2001).
All simulations in this subsection were completed on the SHARCNet cluster Whale. Whale is located at the University of Waterloo and consists of Opteron 2.2 GHz processors (four per node) with a Gigabit Ethernet interconnect. Timing results listed below are given in total cpu time accumulated which is approximately equal to (program runtime) × (number of processors used).
Example 1. (Illustration of Bias and Convergence) The swing option in this example has one up and one down swing right, three exercise opportunities, exercise volume of 60 units of the underlying and there is no penalty. The initial price is USD 40. Here we illustrate the effect of the branching factor b on the value estimates. Specifically, we perform R repeated valuations of the FOST with a branching factor of b and hold the total sample size fixed using the relation R = 32,000 10 b . Figure 2 plots the FOST estimates versus branching factor (the estimates are the averages of the R repeated valuations). Taking the Binomial estimate as the true option value, it is clear that the high estimator overestimates the true price while the low estimator underestimates the true price. Furthermore, as the branching factor increases, the high estimator decreases and the low estimator increases to the true option price, clearly illustrating estimator convergence. Estimator standard errors are approximately 0.07% of estimator value. Example 2. (Effect of Usage Penalty and Initial Asset Price) In this example, the option has two up and two down swing rights and, upon exercise there are three volume choices-20, 40 and 60 units of the asset. Should the final net volume exercised exceed 90 units or be below −90 units a penalty is imposed. The penalty is calculated by multiplying the terminal asset price by ten times the excess usage above 90 units or below −90 units. To see the effect of the penalty on option value, we also turn the penalty term off and value the corresponding option with no penalty. The initial asset value ranges from USD 20 to USD 60 in steps of USD 10. There are m = 5 exercise opportunities and we use a branching factor of b = 20 with R = 4000 repeated valuations.
The pricing results are presented in Table 2. In each row of the table, we see that the high and low estimators bound the true price. Unsurprisingly, imposing a penalty on the cumulative volume reduces the option value. Furthermore, as the initial stock price increases, the increase in an up swing right's value is more than the decrease in a down swing right's value. The opposite is true as the initial stock price decreases. The end result is that as the initial stock price moves away from being at-the-money, the option value increases.
The average computing time per row (not including the binomial forest valuation) was 5.6 h for the cases with usage penalty and was 1.1 h without usage penalty. The reduction in runtime for the case with no penalty can be described as follows. If there are no constraints (e.g., penalty, storage) on the option holder then upon exercise it is always optimal to choose the maximum amount. Therefore with no penalty this option is equivalent to that of an otherwise identical swing option with no volume choices and an exercise volume of 60 units. The latter has fewer trees in its forest and is therefore quicker to evaluate. In Table 2 we have chosen to exploit this as a convenient way to save computational time. For the binomial method run times were on the order of a few seconds.  Example 3. (Effect of Number of Exercise Rights) In this example we illustrate that the option value increases with the number of exercise rights and compare the swing option value with that of a corresponding basket of American options. The option has m = 5 exercise opportunities, an exercise volume of 60 units, and there is no usage penalty. Additionally, we set the initial price to S 0 = 40 and use a branching factor of b = 20 with R = 4000 repeated valuations. We consider options having an equal number of up and down swing rights. Table 3 gives the option price estimates for N u = N d = 1, 3, 5 along with prices computed using the Forest of Trees. First notice that the high and low estimates bound the true option from the binomial model. Next note that with N u = N d = 5 exercise rights, the high and low estimates are exactly the same. In this case the number of exercise opportunities is equal to the numbers of up and down swing rights and since the up and down swing strikes are equal, exactly one of these rights will be exercised at each opportunity (see Equation (19)). This makes both the exercise payoff and the continuation value estimates exactly the same at all times and along all branches for both the low and high estimators, yielding identical prices. Second, the option value increases with the number of swing rights. However, the price increases by a factor that is less than the increase in the number of swing rights. For example, when the number rights increases by a factor of 3 (e.g., going from N u = N d = 1 to N u = N d = 3) the option value increases by a factor of 2.5 and when the number of rights increases by a factor of 5 3 (e.g., going from N u = N d = 3 to N u = N d = 5), the option value increases by a factor of 1.2. This result matches with the intuition that a swing option with a given number of rights is less valuable than a basket of American put and call options with otherwise identical parameters and K d ≤ K u . For the case of one up and one down swing right the basket of American options would contain a single call corresponding to the up swing right and a single put corresponding to the down swing right, with equal strike prices for the call and put. Changing the number of up and down swing rights invokes a change in the corresponding number of call and put options in the basket. Figure 3 shows a comparison between the values of a basket of American options and a swing option with a comparable number of exercise rights. The value of the basket is linear in the number of exercise rights and the swing option value is below the value of the corresponding basket when the number of rights is greater than one. This follows from the restriction that only one swing right may be exercised at any exercise opportunity whereas all American options of a particular type could be exercised at a given time. In the case of one up and one down right the two are equal since it would never be optimal to exercise both the put and call style rights at the same time. As the number of rights increases, the difference in values increases due to the swing option restriction allowing only a single right to be exercises at each opportunity. The low-biased estimator is used for both the basket and swing option values in Figure 3.

Calibrated Forward Curve
In this example, we use the trinomial-tree model described in Section 4 of Jaillet et al. (2004) from which price movements are simulated. This model is a 1-factor model with mean reversion that is seasonally adjusted and calibrated to the forward curve. The option we value is Example (a) of Section 4.2 in Jaillet et al. (2004). This option a two up right swing option with each right allowing the holder to take delivery of either 1 or 2 MMBTus of natural gas. We simplify to have four exercise opportunities and 4 months until expiry. Upon exercise the holder gets where U i is the volume chosen, A i is the seasonality factor and S i is the deseasonalized spot price. Figure 4 plots the option price estimates, including 95% confidence intervals, against branching factor. The number of repeated valuations used to generate prices with a branching factor of b is R = 160,000 b . We see that, with a branching factor of only 8, the confidence intervals for the highand low-biased estimators begin to overlap and quickly become almost indistinguishable for higher branching factors, numerical illustration of both estimators' consistency. Additionally, the high and low estimators converge to the true price computed using the trinomial model. We note that the serial computational times (for a single FOST valuation) for branching factors of 8 and 32 were approximately 4.5 and 110 s, respectively using a 2.1 Ghz Core 2 Duo processor. The pricing results shown in Figure 4 are consistent with the results in Jaillet et al. (2004) but we note that the valuation method in that publication breaks down in higher dimensions and in cases where the inclusion of more risk factors is desirable.

Five Dimensions
Due to the computationally intensive nature of this method it only becomes truly useful in cases where PDE or tree based methods fail. In this subsection we present high-dimensional versions of the examples presented in Section 2.4.1. The asset prices are assumed to follow a risk neutralized correlated geometric Brownian motion described by the stochastic differential equations, where Z k i is a standard Brownian motion process and the instantaneous correlation between Z k and Z s is ρ ks . Here the parameter values are specified as r = 0.05, δ k = δ = 0.1, σ k = σ = 0.2 for all k and that ρ ks = 0 for all k = s. In addition all assets have the same initial value, S 0 , and we take the number of assets to be d = 5.

The swing options considered have both up and down swing rights and the payoff upon exercise is
where S k is the price of the kth underlying asset at the exercise time. This payoff is an extension of the example given in Broadie and Glasserman (1997) and Broadie and Glasserman (2004) for single-exercise American-style options. For the examples considered here, we set K u = K d = K which simplifies the payoff function to As in Section 2.4.1, the option expiry is 3.0 years and the options have both up and down swing rights with strike prices K u = K d = 40.0, respectively. In examples where the holder controls the amount exercised, a list of volume choices is given. Note that we present results from our FOST methodology without comparisons to other methods as there is no generally accepted benchmark for the examples considered here.

Example 4. (Illustration of Bias and Convergence) This 5-dimensional example corresponds with Example 1.
The swing option has one up and one down swing right, three exercise opportunities, exercise volume of 60 units of the underlying and there is no penalty. The initial price is USD 40. We perform R repeated valuations of the FOST with a branching factor of b and hold the total sample size fixed using the relation R = 32000 10 b . This results in standard errors ≈ 0.09% of option value. Figure 5 plots the FOST estimates versus branching, showing that the high estimator overestimates the option price while the low estimator underestimates the price. Furthermore, as the branching factor increases, the high estimator decreases and the low estimator increases and they appear to be converging to the same value, illustrating estimator convergence. These findings are consistent with those in Example 1.  Table 4. The effects of a usage penalty and initial stock price are qualitatively the same compared with the 1-dimensional results. We note that the option value estimates in this example are higher than those in Example 2 due to the payoffs depending on the maximum of the five asset prices. The computing times for the 5-dimensional case are similar to those for the 1-dimensional asset.
Example 6. (Effect of Number of Exercise Rights) This 5-dimensional example corresponds with Example 3, with the option specifications the same as presented there, modulo the adjustment to the payoff function to five dimensions. The results are given in Table 5 and Figure 6. The results, intuition, and interpretation are qualitatively the same as the 1-dimensional results.

Algorithmic Enhancement via Parallel Processing
One method for enhancing the computational efficiency of this algorithm is by taking advantage of multi-processor computing techniques. The simplest and most obvious implementation would be to parallelize across repeated valuations of the forest resulting in serial farming of the repeated valuations. Since each repeated valuation results in an iid random value for the option estimate, the generation of all the results may be completed independently of one another, removing the need for communication between processors. This method is simple and effective. However we state here without numerical evidence that it results in a near perfect speed up without the need for expensive interconnections. With this method the minimum run time that can be produced is determined by the number of processors available, the number of repeated valuations necessary for the desired accuracy and the run time of a single forest.
A variation on the aforementioned parallel implementation is to parallelize the FOST computations internally within the forest. In the results shown in Figure 7 the FOST algorithm has been modified so that the computation of the individual trees within the forest is done using multiple processors. Here we have begun the parallelization after the first time step by dividing up the computation of the remaining subtrees across different processors. Upon completion, the results are gathered and the option value at the initial time step is determined. In Figure 7 we see that this method results in a near perfect speed up due to the small ratio of communication time versus computational time. This implementation may be combined with serial farming resulting in further computational time efficiency. This is discussed more fully in .
In Figure 7 the swing option is identical one in Example 4 and pricing is done with a branching factor b = 160. The computational times were generated using the SHARCNET cluster Hound which comprises 2.2 GHz Opteron processors with 4 GB per core and Infini-Band interconnections. Run times are normalized to the run time of a single processor.

Discussion and Conclusion
The FOST can be thought of as generalizations of two existing pricing methodologies. First, it generalizes the Forest of Trees method to a high-dimensional underlying. Second, it generalizes the Stochastic Tree pricing method for single-exercise right options to multiple exercise rights. We construct high and low FOST estimators analogous to those defined for the Stochastic Tree. We prove properties regarding FOST estimator bias, ordering, and convergence and present numerical results as illustrations.
In related work, we have replaced the binomial/trinomial trees in the Forest of Trees method with Stochastic Meshes Broadie and Glasserman (2004), creating the Forest of Stochastic Meshes . This avoids the exponential growth in computing time with the number of exercise opportunities experienced by the FOST. Another avenue of future work involves algorithmic enhancement. In Section 2.5 we discussed the use of parallel processing to reduce computing time. Two alternatives to this are variance reduction and bias reduction. There are some standard variance reduction methods (e.g., antithetic variates, control variates) that could be used to produce more efficient estimators. The bias reduction technique given in Whitehead et al. (2012) for Stochastic Tree estimators successfully reduces the branching factor required to obtain a desired accuracy for an American option value. This technique can be extended to correct the bias in FOST estimators and we have preliminary evidence of its effectiveness Marshall (2012). Combinations of variance reduction, bias reduction, and parallel processing can be investigated to further improve the algorithm's performance.

Acknowledgments:
The authors thank SHARCNet for computational resources and technical support, particular acknowledgement goes to both Tyson Whitehead and Baolai Ge. We also thank graduate students and faculty in the Financial Mathematics group at Western University, in particular, Matt Davison, Adam Metzler, Rogemar Mamon, and Lars Stentoft.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. The views represented herein are the authors' own views and do not necessarily represent the views of Bank of Montreal or its affiliates and are not a product of its research.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Nomenclature
In this appendix we describe the notation used in this paper. If X is a random variable, we write ||X|| for the p-norm (E|X| p ) 1/p of X. The conditional p-norm of X on Z i , (E|X| p |Z i ) 1/p , is denoted ||X|| Z i . Here we include a summary of notation used for the proofs contained in this appendix as well as all subsequent appendices.

•
Time is indexed by i for t i , i = 0, 1, . . . , m. • R is the number of repeated valuations of the forest. • b is the branching factor.
• S j i is the spot price vector at time t i for branch j = {j 0 , j 1 , . . . , j i }. For convenience we may suppress the bold superscript if there is no ambiguity in doing so, in these cases S j i+1 refers to the time-t i+1 price along the branch path j = {j 0 , j 1 , . . . , j i , j}. • Z i represents the time-t i history of the set of state variables (S j i , N i , U i ), where we suppress the branching history index.
is the time-t i , state-Z i leave one out low biased estimator which does not include node l at time-t i+1 .
is the time-t i , state-Z i leave one out hold value estimator for exercising u units which does not include node l at time-t i+1 , where u 0 = 0. • u is the time-t i volume exercised. Here u ∈ U i . • D i+1 is the discount factor from t i+1 to t i .
where I {A} is the indicator function for set A.

Appendix B. Proofs of Main Results and Lemmas
Appendix B.1.

Proofs of Main Results
In this section we prove the main results of the paper, Theorems 1-5, presented in Section 2.
Proof of Theorem 1. Here we prove the more general statement that E . . , m. The proof proceeds by backward induction. At expiry the inequality holds trivially sinceV , and proceed to the time-t i case. We have, where the first equality comes from the definition of the high estimator, the first inequality comes from the conditional Jensen's inequality and note that N i+1 = N i − I {u * =0} and U i+1 = u i + u * where u * is the value-maximizing volume choice, the second equality uses the tower law and the fact that D i+1 is Z i -measurable, and the second inequality invokes the inductive hypothesis.

Proof of Theorem 2.
As with the proof of the bias of the high estimator we prove the more general In what follows we defineû * l ∈ U i to be the volume choice which maximizes a particularv il . That is, given Z i and subsequentlyû * l is also independent ofv i+1,l given Z i since it is a function ofĝ il . Now, where in the second equality we have used the conditional independence ofĝ il andv i+1,l .
Here p 0 = P(û * l = 0|Z i ) and p j = P(û * l = u j |Z i ) for 1 ≤ j ≤ z and p 0 + . . . + p z = 1. Thus, using the tower law, we have, where the first inequality follows from the inductive hypothesis and the remaining steps follow from the definitions for B i and H i .
Using theĝ il as defined above we first consider the case where for a given tree, is the same for all l (i.e.,û * l =û * , for all l).
where the first inequality comes from the inductive hypothesis and the remaining relations come from the parameter definitions. Next consider the case where the low estimator gives two different estimated optimal exercise amounts,û 1 ,û 2 , across all l branches whereû 1 =û 2 . That isû * l =û 1 orû * l =û 2 for all l = 1, . . . , b. As above we takeû * l to be the optimal exercise amount determined by the l-th leave one out estimator, then, Without loss of generality, suppose thatû * l =û 1 for l = 1, . . . , k andû * l =û 2 for l = k + 1, . . . , b. Then the above ratios become respectively. Now for any i * ≤ k < j * ≤ b we havê which from the definition ofĝ il , implies that This implies that Equation (A3) and similarly for Equation (A4) where the second inequality comes from the inductive hypothesis, the third inequality is an application of Jensen's inequality, the fourth inequality comes from maximizing over a larger set, and the final equality is the definition of the high-biased estimator. For the cases where the low estimator gives z * distinct estimated optimal exercise amounts, u 1 , . . . ,û z * , across all z branches, z * = 3, . . . , z, arguments similar to those given above (for 2 distinct estimated optimal exercise amounts) show that, Since we restrict the number of volume choices to be finite, the theorem is proven.
Prior to proving Theorems 4 and 5 we first state the following preliminary result. Lemma A1. If ||h i (S i , N i , U i , u)|| < ∞ for all t i , for some p ≥ 1, then the following are true for all 0 ≤ t i ≤ t k ≤ t m : The proof of this lemma can be found in the appendix. A second preliminary result is Lemma A2. Let a 1 , . . . , a n , b 1 , . . . , b n , c 1 , . . . , c n be real numbers. Then, A n ≡ |max (a 1 + b 1 , . . . , a n + b n ) − max (a 1 + c 1 , . . . , a n + c n )| ≤ 2 Its proof can also be found in the appendix. We now prove Theorems 4 and 5.
Proof of Theorem 4. Here we take R = 1 and state that if the convergence holds for a single realization of the forest then it will hold for the mean of any number of realizations due to the independence of each repeated valuation. Here we prove by backward induction the more where the first equality comes from the definitions of the estimator and the true value. The third step comes as a result of Lemma A2, the fourth step comes from a generalization of the triangle inequality.
In the final step we rewrite the expression for convience in what follows. First we deal with the C k 's. Given Z i we have that D i+1 B i+1 (S j i+1 , N i − I {u k =0} , U i + u k ) for j = 1, . . . , b and k = 0, . . . , z are iid with means of H i (S i , N i − I {u k =0} , U i + u k ) and finite p-norms. Then by Theorem I.4.1 of Gut (1988) we have that all C k 's in the above expression go to zero.
Next we consider the E k 's. Here we have, by the properties of p-norms and the fact that the terms being averaged are iid, that since E k is bounded by the p-norm of any one of the terms being averaged. By the inductive hypothesis Also by a standard condition for uniform integrability (see Gut (1988) p. 178) we have that for some . From Lemma A1 we know that sup b and we note that by (A7). Hence we have proven the consistency of the low-biased estimator.

Appendix B.2. Lemma Proofs
Proof of Lemma A1. If every h i (S i , N i , U i , u) has finite p-th moment, then each ||h i (S k , N k , U k , u)|| Z i is finite. Since the max, discounting, and conditional expectation operators preserve finiteness of moments then it follows that ||B k (S k , N k , U k )|| Z i and also ||H k (S k , N k , U k )|| Z i must also be finite. Proceeding to (A6), fix t i and proceed by backward induction on t k from t m to t i . At expiry (A6) follows from (A5). Then for t k < t m , where h k (S k , N k , U k , 0) = 0. This is the sum of a finite number of terms, each of which is finite. For (iii) the proof is similar to that of (ii).
2. a 2 + b 2 > a 1 + b 1 (a) a 2 + c 2 > a 1 + c 1 Then a 2 + c 2 < a 1 + c 1 Note that conditions (ii) and (b) imply that b 1 − b 2 < a 2 − a 1 < c 1 − c 2 (A12) and we have where again the first inequality comes from the triangle inequality, the second comes from Inequality (A12) and the third inequality comes from another application of the triangle inequality.
Therefore A 2 ≤ B 2 . Now assume that the inductive hypothesis A n ≤ B n is true. We need to show that A n+1 ≤ B n+1 . First define i n and j n such that a i n + b i n = max(a 1 + b 1 , . . . , a n + b n ) and a i n + c i n = max(a 1 + c 1 , . . . , a n + c n ) respectively. Now, A n+1 = |max(a 1 + b 1 , . . . , a n + b n , a n+1 + b n+1 ) − max(a 1 + c 1 , . . . , a n + c n , a n+1 + c n+1 )| = max(a i n + b i n , a n+1 + b n+1 ) − max(a j n + c j n , a n+1 + c n+1 ) Consider the following, 1. a i n + b i n > a n+1 + b n+1 (a) a j n + c j n > a n+1 + c n+1 where the first inequality comes from the inductive hypothesis. (b) a j n + c j n < a n+1 + c n+1 By the definitions of i n and j n and (b) we have a i n + c i n ≤ a j n + c j n < a n+1 + c n+1 .
This combined with (i) gives b n+1 − b i n < a i n − a n+1 < c n+1 − c i n (A13) where again the first inequality comes from the triangle inequality, the second comes from Inequality (A13) and the third inequality comes from another application of the triangle inequality.
2. a i n + b i n < a n+1 + b n+1 (a) a j n + c j n < a n+1 + c n+1 (b) a j n + c j n > a n+1 + c n+1 By the definitions of i n and j n and (ii) we have a j n + b j n ≤ a i n + b i n < a n+1 + b n+1 .
This combined with (b) gives b j n − b n+1 < a n+1 − a j n < c j n − c n+1 (A14) Then where again the first inequality comes from the triangle inequality, the second comes from Inequality (A14) and the third inequality comes from another application of the triangle inequality.
Therefore A n+1 ≤ B n+1 and the Lemma is proven.