|
Article Excerpt
Adaptive Monte Carlo methods are simulation efficiency improvement techniques designed to adaptively tune simulation estimators. Most of the work on adaptive Monte Carlo methods has been devoted to adaptively tuning importance sampling schemes. We instead focus on adaptive methods based on control variate schemes. We introduce two adaptive control variate methods, and develop their asymptotic properties. The first method uses stochastic approximation to adaptively tune control variate estimators. It is easy to implement, but it requires some nontrivial tuning of parameters. The second method is based on sample average approximation. Tuning is no longer required, but it can be computationally expensive. Numerical results for the pricing of barrier options are presented to demonstrate the methods.
Key words: variance reduction; adaptive Monte Carlo; control variates
1. Introduction. Suppose that we wish to estimate EX, where X is a real-valued random variable. Suppose also that {Y([theta]): [theta] [member of] [THETA]} is a parametric collection of random variables such that EY([theta]) = for any [theta] in the parameter set [THETA]. Then X - Y([theta]) is an unbiased estimator for EX, where Y([theta]) serves as a control variate, and the parameter [theta] can be selected so as to minimize the variance of X - Y([theta]). When the parameterization is linear, we can appeal to the standard theory of (linear) control variates, e.g., Law and Kelton [14]. Identifying the [theta] that minimizes the variance is straightforward in the linear case because the variance is a convex quadratic in [theta]. When the parameterization is nonlinear, the problem is not so straightforward. We propose two adaptive procedures that tune the parameter [theta] while estimating EX, and study the large sample properties of the procedures.
It was shown in Glynn and Whitt [4] that, from an asymptotic perspective (asymptotic in the simulation runlength), there is no advantage gained by nonlinear parameterizations over linear parameterizations. However, the form of the parameterization considered in Glynn and Whitt [4] is different from the form considered in this paper. Glynn and Whitt [4] analyze the asymptotic performance of control variate estimators that have the form f([[bar.X].sub.n], [[bar.Y].sub.n]) where n is the simulation runlength; ([[bar.X].sub.n] and [[bar.Y].sub.n] are sample averages of X and the control Y, respectively; and f is a function that combines the two sample averages. They show that, from the asymptotic point of view, we may restrict the choice of f to linear functions. Indeed, in our parametrization, X and the control variate Y([theta]) are linearly related, so the negative result does not apply.
Our interest in this problem stems from several application areas. An extended example in this paper (see [section][section]2 and 6) shows that one can apply adaptive control variates in pricing certain financial derivatives. A second example arises in the simulation analysis of multiclass processing networks. When these networks are heavily loaded, simulation estimators can suffer from large variance, so some form of variance reduction is needed. The simulation estimators developed in Henderson and Meyn [6, 7] give large variance reductions, but the asymptotic rates of growth in the variance are the same as for the naive estimator; see Meyn [ 18]. One approach to potentially improving these estimators is to develop parameterized estimators. A third class of examples arises in the problem of estimating the expected cost to absorption in a Markov chain. This problem has received a great deal of attention because of its applications in radiation transport problems; see, e.g., Kollman et al. [12]. The common thread underlying these applications is that they involve the simulation of a Markov process. This allows us to construct a parameterized family of control variates using the approximating martingales developed in Henderson and Glynn [5]. Once we have a parameterized class of control variates at hand we need a procedure for selecting a control from within the class.
The first of our procedures is based on a stochastic approximation (SA) scheme. At iteration k, several independent replications of X - Y([[theta].sub.k-1]) are generated, conditional on the parameter choice [[theta].sub.k-1] from the previous iteration. The sample mean and the gradient (with respect to [theta]) of the sample variance are then computed, and the parameter [[theta].sub.k-1] is updated to [[theta].sub.k] in an SA step. This procedure is easily implemented and exhibits good performance with appropriately chosen step sizes. But the selection of the step size is nontrivial, and the finite-time performance of the algorithm is strongly affected by the choice of step sizes.
The second procedure is based on the theory of sample average approximation (SAA). In an initial stage, a random sample is generated, and a sample variance function is defined with the generated sample. The sample variance function is deterministic in terms of the parameter [theta], and the optimal value of [theta] that minimizes this sample variance function is determined using a nonlinear optimization solver. Then one makes a production run using the value of [theta] returned in the first stage. The initial optimization can be computationally expensive when compared with one step of the SA procedure. However, SAA does not require tuning parameters beyond the choice of run length, and for very long simulation runs, a vanishingly small fraction of the effort is required in the initial optimization.
Henderson and Simon [8] also develop an adaptive control variate method for finite-horizon simulations. They give conditions under which adaptive control variate estimators converge at an exponential rate. One of the key assumptions there is the existence of a "perfect" control variate, i.e., a parameter value [[theta].sup.*] such that var(X - Y([[theta].sup.*])) = 0. For the applications we have in mind, this assumption is unlikely to hold. Bolia and Juneja [2] use the martingale control variates developed in Henderson and Glynn [5], as we do, but in the case of linearly parameterized controls. Maire [17] expresses the estimation problem as an integration problem over the unit hypercube and uses the expansion of the integrand for an approximate orthonormal basis as a control variate. An iterative procedure estimates the coefficients of the expansion so that the variance of each estimated coefficient has a polynomial decay. The residual terms are not estimated iteratively so, in general, the convergence rate of the procedure cannot exceed the canonical rate. Henderson et al. [9] develop adaptive control variate schemes for Markov chains in the steady-state setting. They use an SA procedure for tuning control variate estimators developed in Tadie and Meyn [22] and provide conditions for minimization of an approximation of the steady-state variance.
The main contribution of this paper is to develop adaptive control variate methods for finite-horizon simulation when nonlinearly parameterized control variables are available and provide conditions under which the adaptive estimators are consistent. In general, it will be the case that the optimal variance is positive; in this paper we focus our attention on this case. Our proposed estimators then satisfy central limit theorems and, as a consequence, converge at the canonical [n.sup.-1/2] rate, where n is proportional to the computational effort. In addition, to the best of our knowledge, this is the first application of SA and SAA methods in this context. Finally, we discuss implementation issues regarding the practical use of these adaptive methods.
This paper is organized as follows. Section 2 sketches some of the main ideas in pricing barrier options using adaptive control variates, as motivation for the remainder of the paper. Section 3 outlines the general problem and discusses estimation of the gradient of the variance function. Section 4 explores the use of SA and [section]5 explores SAA. In [section] 6, we describe the results of some limited experiments with the example of [section]2. Finally, [section]7 contains some concluding remarks.
Unless otherwise stated, all vectors are column vectors and all norms are Euclidean. Suffixes can either indicate different instances of a random vector or components of a single vector, with the context clarifying what is intended.
2. A motivating example. In this section, we describe the problem of pricing barrier options, and explain how parameterized controls can be found. Our goal in this section and in [section] 6 is not to develop the most efficient known estimators for pricing barrier options, but rather to demonstrate the adaptive control variate methodology in a familiar, but nontrivial, setting, and bring out some of the practical issues involved in applications.
A barrier option is a derivative security that is either activated (knocked-in) or extinguished (knocked-out) when the price of the underlying asset reaches a certain level (barrier) at any time during the lifetime of the option; see, e.g., Glasserman [3].
The price of the underlying stock at time t is denoted by S(t), for t [greater than or equal to] 0. Suppose that the underlying stock price is monitored at discrete times [t.sub.i] = i[DELTA]t, i = 0, 1, 2, ..., l, where T is the (deterministic) expiration date of the option and [DELTA]t = T/l is the time between consecutive monitoring dates. For notational convenience, let [S.sub.i] denote the underlying stock price at the ith monitoring point (i.e., S([t.sub.i])). Assume that the initial stock price [S.sub.0] takes a value in an interval H and the barrier is the boundary of H. When the stock price crosses the barrier, the option is knocked out and the payoff is zero. If the option has not been knocked out by time T, then the payoff at time T is [([S.sub.l] - K).sup.+], where K > is the strike price. Hence, the option payoff depends on the complete path {[S.sub.i], i = 0, ..., l}. Define
[tau] = inf{n [greater than or equal to] 0: [S.sub.n][not member of] H} and [A.sub.i] = [1.sub{[tau]>i}, i = 0, ..., l.
Then [A.sub.i] is the indicator that determines whether or not the option is alive at time [t.sub.i]. We assume that the market is arbitrage free. Then the price of a knock-out call option is given by
[e-sup.rT] E[[A.sub.l][([S.sub.l] - K).sup.+],
where r is the (assumed constant) risk-free interest rate and the expectation is taken under the risk-neutral measure. Since the discount factor [e.sup.-rT] is constant, pricing the option reduces to estimating the expected payoff with the initial stock price x, i.e., estimating
E[[A.sub.l][([S.sub.l] - K).sup.+]|[S.sub.0] = x].
Assume that the underlying stock price process {S(t): t [greater than or equal to] 0} is a (time homogeneous) Markov process. Then {[S.sub.n]: n = 0, 1, 2, ...}, where [S.sub.n] is the stock price at time [t.sub.n] = n[DELTA]t, is a discrete time Markov chain on the state space [0, [infinity]). For i = 0, 1, ..., define
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],
so that [U.sup.*](x, i) is the expected payoff of the option with the initial stock price x and maturity [t.sub.i]. Our goal is to estimate [U.sup.*](x, l).
We now describe the martingale that serves as a control variate, drawing from the general results of Henderson...
|
