**1. Introduction**

Simulation and regression based methods are by now standard methods used for American option pricing. Simulation methods are flexible, work in multiple dimensions, can be used with even the most complicated models and payoff functions, and are straightforward to use. Essentially, if you can simulate it (the underlying dynamics), you can price it (the derivatives contract).<sup>1</sup> Moreover, under standard assumptions, the estimates have nice properties since averages over independent samples converge to expected values. Almost all the simulation methods that have been proposed for American option pricing share the feature that in one way or another, they involve estimating or approximating the early exercise boundary or, more generally, estimating or approximating the conditional expectations that are used to decide whether or not to exercise at each step.<sup>2</sup> In the Least-Squares Monte Carlo (LSM) method of Longstaff and Schwartz (2001), for example, the conditional expectations are approximated with polynomials, and based on these approximations,

<sup>1</sup> See also Boyer and Stentoft (2013) for applications of this method to longevity risk products in insurance.

<sup>2</sup> See Ibanez and Zapatero (2004) for an example of a simulation method that directly approximates the optimal early exercise boundary.

it is then decided at each early exercise point and along each simulated path whether the option should be exercised or not; implicitly determining the early exercise region.

While the asymptotic properties of the LSM method are well established (see for example (Stentoft 2004b)), the quality of the price estimate obtained with this method depends on the order of the polynomial, *M*, used in the cross-sectional regression and the number of paths, *N*, used in the simulation (see, e.g., Moreno and Navas 2003; Stentoft 2004a). With a finite choice of polynomial terms, *M*, and simulated paths, *N*, the early exercise boundary is estimated with some error. Because of this, sub-optimal decisions are made, and everything else equal, low biased price estimates are obtained. Moreover, because sub-optimal decisions are made early on, errors are propagated backwards through the algorithm. This paper presents an innovative bias reduction technique directly applicable to the LSM method, as well as to many other similar methods. The method successfully delivers at each early exercise step, starting from just before maturity, a "better" estimated early exercise boundary. Because of this, less errors are accumulated, and by recursively correcting the estimated optimal early exercise boundary, we end up with essentially unbiased price estimates even for long maturity options.

To motivate our methodology, we first demonstrate that significantly less biased estimates can be obtained if the precision of the early exercise boundary is improved upon by simply averaging it over several independently repeated simulations. Next, we show that by averaging at each step, and therefore obtaining a less noisy estimate of the early exercise strategy, less errors are propagated backwards through the algorithm. Our method is based on the simple observation that if one uses the average of the approximated early exercise strategies instead of, say, the *I* = 100 individually determined strategies, the estimated prices are much closer to the benchmark prices. Our proposed method exploits and leverages this observation to the fullest by averaging at each step in the backwards induction algorithm instead of only after having gone through all the steps of the backward induction algorithm. We refer to our method as a bootstrapping approach because of the similarities it has with how the term structure of interest rates is bootstrapped and the way bootstrapping is used for inference in regression models. As a result of the bootstrapping approach, the final bias of the estimated price is almost completely eliminated.

We compare the results from our bootstrapping approach to the regular method of Longstaff and Schwartz (2001) and show that our proposed method works extremely well across various levels of moneyness, maturity, and volatility. Moreover, compared to the regular LSM method, the bootstrapping method is much less sensitive to the choice of polynomial order used in the cross-sectional regressions and to the number of paths used in the Monte Carlo simulation. When a reasonable order, say *M* = 9, is used in the regressions and when a standard number of paths, say *N* = 100,000, is used in the simulation, our proposed method essentially provides unbiased price estimates. We also demonstrate that the technique generalizes straightforwardly to multi-asset options with various payoffs. To illustrate this, we price arithmetic and geometric average options and options on the maximum and minimum of three assets. For all four payoff functions, the results hold, and the bootstrapped method always performs the best and delivers price estimates with negligible or very small bias for reasonable choices of the number of regressors and paths used in the Monte Carlo simulation.

Our method successfully corrects the bias even in very small samples and allows for using, e.g., *N* = 10,000 paths or less to determine the early exercise boundary. Our detailed results also show that one does not need hundreds of repeated simulations to obtain a good estimate of the early exercise boundary, but a much lower number, e.g., *I* = 10 suffices. Determining the early exercise boundary is the most computationally complex part of the algorithm since this involves performing several cross-sectional regressions, and our method therefore offers significant computational speedups. For example, one can run *I* = 10 initial simulations with *N* = 10,000 paths to estimate the optimal early exercise boundary. One can then use *N* = 1,000,000 paths for pricing, whereby an estimate with very low bias and low variance is obtained. Note that in this setup, obtaining the initial estimate of the optimal early exercise boundary takes no more time than running the standard LSM method once, and since it is straightforward to parallelize our method, it could run in a fraction of this time on multicore processors or modern computer clusters.

It should be noted that several papers have proposed alternative refinements to this type of simulation method. In particular, it has been suggested to use well known variance reduction techniques like antithetic simulation, control variates, importance sampling, as well as initial dispersion together with the LSM method (see, among others: (Juneja and Kalra 2009; Lemieux and La 2005; Rasmussen 2005)). However, there are very few alternative suggestions for how to reduce the bias of the simulated estimates of American option prices in general and of the LSM method in particular. One potentially interesting method, which could be combined with our method, is that of Kan et al. (2009), although their application was to the value function iteration method of, e.g., Carriere (1996) and Tsitsiklis and Van Roy (2001). Our proposed method could also be used together with the inequality constrained least-squares method suggested in Létourneau and Stentoft (2014) or the method that corrects for heteroskedasticity in the cross-sectional regression proposed by Fabozzi et al. (2017).

Finally, our current applications are to standard financial options, but our results have important implications for other applications of simulation based pricing as well. One particularly important and challenging area in which the LSM method has been used is the field of real option valuation. For example, the work in Gamba and Fusari (2009) proposed a general valuation approach for capital budgeting decision involving modularization, and the work in Kang and Létourneau (2016) studied the investment decision and choice between coal and natural gas power plants under political risk with options to turn production on or off, while the work in Power et al. (2015) studied the valuation and timing of complex infrastructure projects. In all these cases, the regular LSM method was used to approximate complex decision problems in settings with multiple risk factors with complicated dynamics and for problems where the payoff function is non-standard and "exotic". We conjecture that our proposed bootstrapping method will allow more efficient determination of the optimal controls and more precise valuation of these assets.

The remainder of the paper is structured as follows: Section 2 describes the discrete time framework used for valuation, introduces the least-squares Monte Carlo method of Longstaff and Schwartz (2001), and demonstrates how our proposed bootstrapping method can be implemented and improves estimation of the early exercise boundary. Section 3 presents an extensive numerical analysis of the performance of our proposed method demonstrating its robustness to changes in the simulation setup, to different option characteristics, and its applicability to higher dimensional problems. Section 4 provides a discussion of the potential benefits of separating the estimation of the early exercise boundary from the pricing and demonstrates that our method can be implemented in a very efficient way that works well in empirically relevant settings. Section 5 concludes.
