University of Southern California Dissertations and Theses

Stochastic models: simulation and heavy traffic analysis

(USC Thesis Other)

Stochastic models: simulation and heavy traffic analysis

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content STOCHASTIC MODELS: SIMULATION AND HEAVY TRAFFIC ANALYSIS by Samim Ghamami A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulllment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (INDUSTRIAL AND SYSTEMS ENGINEERING) December 2009 Copyright 2009 Samim Ghamami Dedication To Pouneh and my parents ii Acknowledgements Like many I was rst attracted to probability by Sheldon Ross' \A First Course In Probability" during my undergraduate studies in Tehran. A few years later, I was truly fortunate to become Professor Ross' PhD student. Now, at the end of my PhD, I have acquired the ability to see how these valuable references have been created. These years, what Professor Ross walked me through was well beyond probability. I observed an extraordinary optimism, a profound happiness, and a lucid way of thinking in everyday life. Our meetings and interactions have been memorable and extremely in uential for me. I would like to thank Professor Ross for his continuous support, patience, and guidance. After her PhD course in spring 2007, Professor Amy Ward generously provided a research assistantship opportunity so that I can learn from her and work with her on Heavy Trac analysis. Amy has been a source of encouragement and a valuable friend from whom I continue to learn. Our interactions have been very constructive and grati- fying for me. Our technically demanding problem brought a set of unique consequences including enjoyable trips to San Diego where I beneted from Ruth Williams' insight and advice. I would like to thank Professor Ward for her guidance, patience, continuous support, and friendship. During these years at USC, I have beneted greatly from Professor Jim Moore's support. In spite of his busy schedule, Professor Moore has been a supportive member iii of both my qualifying examination and dissertation committee. I would like to thank Professor Moore for his guidance and support. I am grateful to Professor Ruth Williams for her enlightening advice on weak conver- gence argument of chapter 6. She was extremely generous with her time. I would like to thank Professor Fernando Ordonez and Professor Detlof von Winter- feldt for participating in my qualifying examination. I met with Pouneh, my wife and best friend, early in my life, and since then, almost everyday I discover new exceptional inspiring aspects of her personality. She has been extremely enthusiastic about my PhD studies. None of this would be possible without her endless support. I would like to thank Pouneh for her sel essness and love. I grew up enjoying the loving companionship of my parents. Shohreh, my mother, was more anxious than I while I prepared for the screening and qualifying exams. Majid, my father, still listens to me on the phone, thousands of miles away in Tehran, and helps me view myself and the world more critically. The physical distance is something I confront everyday. I would like to thank my parents for their many sacrices, loving support, and friendship. iv Table of Contents Dedication ii Acknowledgements iii List of Tables viii List of Figures ix Abstract x Chapter 1 : Ecient Monte Carlo Barrier Option Pricing When The Underlying Se- curity Price Follows a Jump-Diusion Process 1 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Estimating the Value of a Knock-Out Option Using M&A's Raw Simula- tion Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2.1 Sequential Simulation of the Jump-Diusion Sample Path based on Jump Times . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2.2 Maximum of a Brownian Bridge . . . . . . . . . . . . . . . . . . . 5 1.2.3 Metwally and Atiya's approach to estimate UOC . . . . . . . . . 6 1.3 Variance Reduction in Estimating the Value of an Up-and-Out Barrier Option . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.3.1 Stratied Sampling . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.3.2 The Proposed Simulation Estimator of E[RjN(T ) =m] . . . . . . 8 1.4 Variance Reduction in Estimating the Value of an Up-and-In Barrier Op- tion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.4.1 The Proposed Simulation Estimator of E[RjN(T ) =m] . . . . . . 14 1.5 Measuring the Eciency of the Proposed Estimators with Numerical Ex- amples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 1.5.1 Monte Carlo Estimation of UOC . . . . . . . . . . . . . . . . . . . 20 1.5.2 Monte Carlo Estimation of UIC . . . . . . . . . . . . . . . . . . . . 24 Chapter 2 : Ecient Simulation of a Random Knockout Tournament 26 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.2 Method 1: Observed Survivals and the Product Estimator . . . . . . . . . 28 2.3 Method 2: Conditional Expectation and Post Stratication . . . . . . . . 30 v 2.4 Classical Knockout Tournament . . . . . . . . . . . . . . . . . . . . . . . . 34 2.5 Numerical examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 Chapter 3 : Ecient Monte Carlo Estimation of System Reliability 39 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.1.1 Improved Stratied Sampling Approach . . . . . . . . . . . . . . . 42 3.2 Binary Dependent Components . . . . . . . . . . . . . . . . . . . . . . . . 43 3.2.1 Random Scan Gibbs Sampler . . . . . . . . . . . . . . . . . . . . . 44 3.2.2 Case 1: Exchangeable Components . . . . . . . . . . . . . . . . . . 45 3.2.3 Case 2: Non-exchangeable Components with a Small Number of Conditional Dependency Functions . . . . . . . . . . . . . . . . . . 47 3.2.4 Case 3: Non-exchangeable Components with a Large Number of Conditional Dependency Functions . . . . . . . . . . . . . . . . . . 50 3.3 3-State Independent Components . . . . . . . . . . . . . . . . . . . . . . . 52 3.3.1 3-State i.i.d Components and the Estimator of E[H] . . . . . . . . 52 3.3.2 The 3-Stage Simulation Procedure . . . . . . . . . . . . . . . . . . 54 Chapter 4 : Dynamic Scheduling of an N-System with Reneging 57 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.1.1 Notation and Terminology . . . . . . . . . . . . . . . . . . . . . . . 63 4.2 The N System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 4.2.1 Stochastic Primitives . . . . . . . . . . . . . . . . . . . . . . . . . . 65 4.2.2 Scheduling Control . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.3 Heavy Trac Assumptions, Scaling, and the Cost Function . . . . . . . . 69 Chapter 5 : Approximating Brownian Control Problem and the Proposed Policy 76 5.1 The Brownian Control Problem Solution . . . . . . . . . . . . . . . . . . . 81 5.1.1 The Bellman Equation . . . . . . . . . . . . . . . . . . . . . . . . . 81 5.1.2 The Case that a Threshold Control Policy Solves the Brownian Control Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 5.1.3 The Case that a Static Policy Solves the Brownian Control Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 5.2 The Proposed Policy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 Chapter 6 : Weak Convergence and Asymptotic Optimality 96 6.1 Weak Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 6.1.1 Proof of Proposition 6.1 . . . . . . . . . . . . . . . . . . . . . . . . 98 6.1.2 Proof of Proposition 6.2 . . . . . . . . . . . . . . . . . . . . . . . 105 6.1.3 Proof of Theorem 6.1 . . . . . . . . . . . . . . . . . . . . . . . . . 113 6.2 Asymptotic Optimality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 Chapter 7 : Proof of the Technical Lemmas 124 7.1 An Oscillation Inequality for Solutions of a One-dimensional, Linearly Gen- eralized, Perturbed Skorokhod Problem . . . . . . . . . . . . . . . . . . . 124 7.2 The Details Omitted from the Proof of Proposition 6.1 . . . . . . . . . . . 127 7.3 Proofs of Lemmas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 vi Bibliography 164 vii List of Tables Table 1.1 Optimal Stratication . . . . . . . . . . . . . . . . . . . . . . . . . . 23 Table 1.2 Example 1. 10 5 simulation runs, B=45, b C = 3:4774. \Method 1" is the most ecient method. . . . . . . . . . . . . . . . . . . . . . . . . . . 23 Table 1.3 Example 2. 10 5 simulation runs, B=60, b C = 3:4822. \Method 1" is the most ecient method. . . . . . . . . . . . . . . . . . . . . . . . . . . 23 Table 1.4 Example 3. 10 5 simulation runs, B=90, b C = 3:4435. d UOC is very close to b C. \Method 2" combined with our estimator of C is the most ecient method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 Table 1.5 Example 1. 10 5 simulation runs, B=75, b C = 7:1299. \Method 1" combined with our estimator of C is the most ecient method. . . . . . . 25 Table 1.6 Example 2. 10 5 simulation runs, B=100, b C = 7:1316. \Method 1" combined with our estimator of C is the most ecient method. . . . . . . 25 Table 1.7 Example 3. 10 5 simulation runs, B=150, b C = 7:0964. d UOC is very close to b C. \Method 2" is the most ecient method. . . . . . . . . . . . . 25 Table 2.1 Estimates of P i 's . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 Table 2.2 Variance of the estimators based on 10 5 simulation runs (10 9 scale) 37 Table 2.3 Elapsed time (in second) . . . . . . . . . . . . . . . . . . . . . . . . 37 Table 2.4 Estimates of P i 's . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 Table 2.5 Variance of the estimators based on 10 5 simulation runs (10 9 scale) 38 Table 2.6 Elapsed time (in second) . . . . . . . . . . . . . . . . . . . . . . . . 38 Table 5.1 Our proposed policy. . . . . . . . . . . . . . . . . . . . . . . . . . . 93 viii List of Figures Figure 4.1 The N-system with reneging. . . . . . . . . . . . . . . . . . . . . . 59 ix Abstract This dissertation consists of two main parts: Monte Carlo Simulation and Heavy Trac analysis. Much of the simulation part, chapter 1-3, is devoted to ways of improving sim- ulation estimators in 3 dierent problems: Barrier Option Pricing, Random Knockout Tournaments, and System Reliability. Our eciency criteria for comparing alternative estimators are variance and computing time. To achieve variance reduction, in addition to using combinations of dierent classical variance reduction techniques, innovative sub- stantial sources of variance reduction are introduced by exploiting specic features of the problems under consideration. Second part of this dissertation, chapter 4-7, focuses on analyzing a stochastic control problem associated with a parallel server queueing sys- tem by employing Heavy Trac approach. The stochastic control problem can not be solved exactly. We use heavy trac limit theorems to derive an approximating control problem, referred to as Brownian control problem. Based on the Brownian control prob- lem solution, we propose a control policy and then prove that our proposed policy is asymptotically optimal for the original control problem. Part 1: Ecient Monte Carlo Simulation It is well-known that continuous-sample- path and constant volatility assumptions of (geometric) Brownian motion used in the classical Black-Scholes framework are inadequate assumptions to characterize nancial markets. Market data reveal that logarithm of a price process has a high peak and heavy tails. Variations of jump-diusion (Levy) processes and stochastic volatility models have x been extensively used to capture these empirical observations (see Cont and Tankov [10]). In chapter 1, we consider the problem of valuating barrier call options when the underlying security price follows a geometric Brownian motion with mutiplicative jumps (the jump part is a compound poisson process). Analytical formulas are very dicult to derive in this setting. Metwally and Atiya [44] introduce a Monte Carlo approach for pricing knock-out barrier options in the same setting. No variance reduction, however, has been introduced in M&A's raw simulation approach. Our proposed estimators for both types of barrier options, i.e., knock-out and knock-in barrier call options, achieve substantial variance reduction. Recall that a knock-in option becomes alive only if a barrier value is crossed by the underlying security price during the life of the option, and a knock-out option remains alive only if the barrier is never crossed during the life of the option. A signicant source of variance reduction in deriving our estimator of knock-out barrier options follows from conditioning the expected return of the knock-out option on a given value of security price at expiration time which lies between the strike price and the barrier value. Then, we illustrate how to generate the jump-diusion process sequentially based on jump times and conditional on this given sampled value. To derive our estimator of knock-in options, we rst condition the expected return on jump times. Suppose there aren jumps and jump times are indexed by i. Then, we illustrate that given jump times the expected return of the knock-in option can be expressed as the summation of expected returns of n knock-in options where the i th option becomes alive only if the barrier is crossed in the i th inter-jump interval. Next, a sequential conditional expectation is used to estimate the aforementioned identity. Finally, we use the parity between vanilla and barrier call options, both of our proposed estimators, and numerical examples to carry xi out a general rule for pricing barrier options. This general rule identies the condition under which each of our proposed simulation methods are most ecient in terms of both variance and computing time. In chapter 2, we consider the problem of using simulation to eciently estimate the win probabilities for participants in a general random knockout tournament. In a random knockout tournament players are paired at random at the outset of each round and play head-to-head in matches; with losers being eliminated from the tournament. Let P ij = 1P ji denote the probability that playeri defeats playerj if they are paired in a round. The assumption isP ij 's are given, and the outcome of a match between 2 players depends only on the pair. We develop two very ecient Monte Carlo estimators to estimate win probabilities in a general random knockout tournament, one based on the notion of \observed survivals" and the other based on an innovative use of conditional expectation and post-stratication (referred to as method 1 and method 2, respectively). In the special case of a 2 n -player classical random knockout tournament a second conditional expectation-based estimator is also introduced. In [53], Ross and Schechner, use the notion of \observed hazards" to introduce very ecient Monte Carlo estimators for the problem of rst passage distribution of a Markov chain to a given set of states. Our rst ecient estimator is derived based on the notion of \observed hazards" in [53]. Our numerical examples show that method 2 leads to more variance reduction. However, in terms of both variance reduction and computing time, we conclude that method 1, i.e., the method of observed survivals is the most ecient Monte Carlo method to estimate win probabilities. xii In chapter 3, we consider the problem of Monte Carlo estimation of system reliabil- ity. Specifying the probability that a component system functions, i.e., system reliability (Barlow and Proschan [4]), is often a computationally complex problem. Monte Carlo simulation has been extensively used for the estimation of system reliability. In [48], Ross introduces a very ecient stratied-sampling based simulation approach to esti- mate the expected values of functions of Bernoulli random variables under certain types of dependencies. Ross' improved stratied sampling approach (ISS) plays a key role in deriving our ecient Monte Carlo estimators. First, we consider a system of binary de- pendent components. We model the dependency among system components by assuming that whether each component functions or not dependents on the remaining number of working components. Then, we modify ISS to estimate system reliability under this prob- abilistic structure. To utilize ISS we use the notion of Gibbs Sampler (see chapter 6 of Liu [39]) to analytically compute stratication probabilities, i.e., the probability mass function of our stratication variable. Next, in a dierent setting, we consider a system with 3-state independent components. There, we introduce a 3-stage simulation proce- dure which is in fact the generalization of Ross' ISS from binary to 3-state independent components. Part 2: Heavy Trac analysis: Dynamic Scheduling of an N-System with Reneging N-system is a simple and well-known class of a more general parallel-server queueing systems which arise in a variety of applications. In an N-system, there are two independent renewal arrival streams of two classes of jobs to two innite capacity buers. Jobs will abandon the system independently if not being served within an exponential xiii amount of time. Server 1 can only serve class 1, and server 2 can serve both classes. We dene class dependent holding and abandonment costs. The objective is dynamic scheduling of waiting jobs on to available servers to minimize the expected innite horizon discounted holding and reneging cost. Bell and Williams in [5] address the same problem using heavy trac approach but in the absence of reneging. Our scheduling control problem associated with the N-system can not be solved ex- actly, and so using heavy trac approach we solve it approximately. We rst, in chapter 4, dene the heavy trac and complete resource pooling assumptions for the N-system and derive the diusion and uid scaled processes associated with the N-system which are indexed by r > 0. In words, we would like to analyze the normalized (diusion-scaled) N-system at long intervals of time when it is heavily loaded and balanced. To do so, in chapter 5, we derive our approximating multidimensional Brownian control problem by formally passingr to innity which achieves the N-system in heavy trac limit. A useful feature of Brownian control problems associated with parallel server systems is that they can be converted to an eectively one-dimensional equivalent control problem referred to as workload control problem (see [29] and [28]). In our setting workload process is a specic linear combination of queue-length 1 (number of class 1 jobs present in buer 1 at a specic time) and queue-length 2 processes. In contrast to Bell and Williams' setting where the workload control problem admits a straightforward path-wise solution, the in- clusion of reneging leads us to solve a non-trivial control problem. In a specic parameter regime where the c- rule of scheduling theory suggests buer 1 is more expensive than buer 2 (and so server 2 should give priority to class 1 over class 2 jobs), and also reneging rate of class 1 jobs is larger than reneging rate of class 2 jobs according to (4.2) in chapter xiv 4, BCP solution suggests placing a threshold on the workload process. In this parameter regime, our proposed policy suggests monitoring the N-system by placing one threshold on workload and a second smaller-sized threshold on buer 1. Server 2 gives priority to class 1 only if workload is below the workload threshold and the number of class 1 jobs in buer 1 exceeds the safety stock buer threshold. Otherwise, server 2 gives priority to class 2 jobs (Our proposed policy in chapter 5 is exhaustive in that it cover the whole parameter space based on BCP solution). Comparing our BCP, BCP solution, and pro- posed policy to that of Bell and Williams in a comparable parameter regime, we conclude that the inclusion of reneging fundamentally changes the asymptotically optimal policy derived from heavy trac analysis of the control problem associated with the N-system. Chapter 6 is devoted to verifying that our proposed policy derived from the BCP solution is an asymptotically optimal policy for the original control problem. We rst in Theorem (6.1) prove that for an N-system operating under our proposed policy, diusion scaled queue-length and idle-time processes weakly converge to the BCP solution as r! 1. Next, using our weak convergence result, we show that when the N-system operates under any arbitrary scheduling control policy (our proposed control policy), \lim inf" (lim) of a sequence of rindexed diusion-scaled cost functions is larger than (equal to) the optimal value of the cost function achieved from the BCP solution, Theorem (6.2). Chapter 7 is devoted to the proof of Lemmas in chapter 4-6. xv Chapter 1 Ecient Monte Carlo Barrier Option Pricing When The Underlying Security Price Follows a Jump-Diusion Process 1.1. Introduction Let S(t) denote the price of a security at time t, and suppose that S(t) =S (t)J(t) ; t 0 (1.1) whereJ(t) = Q N(t) i=1 J i ,N(t) is a Poisson process with meant,J i 's are independent ran- dom variables having a common general probability distribution, and Q N(t) i=1 J i is dened to equal 1 when N(t) = 0. S (t) is a geometric Brownian motion with volatility param- eter and drift parameter . It is assumed that J i 's, N(t), and S (t) are independent. Throughout this chapter, we occasionally refer toS(t) as the jump-diusion process. Un- der the preceding jump-diusion framework, let C, UIC, and UOC represent the expected discounted payo of a European call option, a European up-and-in barrier call option, 1 and a European up-and-out barrier call option with strike priceK, expiration timeT , and initial price s, respectively. Let B denote the barrier value, and recall that an up-and-in barrier option becomes alive only if the underlying security's price crosses the barrier during the life of the option, and an up-and-out barrier option remains alive only if the barrier is never crossed during the life of the option. That is, C = e rT E[(S(T )K) + ] and, UIC =e rT E[(S(T )K) + 1f B Tg] UOC =e rT E[(S(T )K) + 1f B >Tg] (1.2) where s K B, 1fg denotes the indicator of the event in the braces, B = infft > 0 : S(t) Bg, and r is the continuously compounded interest rate. For the purpose of risk-neutral pricing, drift of the preceding geometric Brownian motion satises = r 2 =2+E[J] by the usual risk neutral assumption, i.e. E[S(t)] =S(0)e rt , (see [51] or [22]). It is well known that deriving analytical closed form formulas for the valuation of vanilla and specially exotic options under the aforementioned jump-diusion framework is very dicult. Monte Carlo methods are often used for pricing these derivatives. In our study, we are mainly interested in developing ecient simulation methods to estimate the value of UOC and UIC dened above. Note that down-and-in and down-and-out call options are dened and their value is estimated analogously, and so our ecient simulation procedures cover the two general types of barrier options, i.e., knock-out and knock-in barrier options. In [44], Metwally and Atiya introduce a simulation approach for pricing knock-out barrier options under the preceding jump-diusion framework. No variance reduction, however, has been introduced in M&A' raw simulation approach. Our 2 criteria for eciency of simulation estimators in the sequel are variance and computing time. We signicantly improve upon the proposed Monte Carlo estimators in [44] both in terms of variance and computing time . In section 1.2, we summarize the simulation method proposed by Metwally and Atiya to estimate the value of an up-and-out call option under the preceding jump-diusion framework, and we refer to it as the M&A's raw simulation approach. In section 1.3 and 1.4 we illustrate how to achieve signicant variance reduction when estimating the price of knock-out and knock-in barrier options, respectively. Stratied sampling with N(T ) as the \stratication variable" is used in both section 1.2 and 1.3. The innovative source of variance reduction in deriving the estimator of the value of a knock-out barrier option, section 1.3, originates from conditioning the expected dis- counted payo of the knock-out option on a given value of the underlying security at expiration time which lies between the strike price and the barrier value. Then, we show how to sequentially generate the sample path of the jump-diusion process at jump times conditional on this sampled value of S(T ). In section 1.4, we derive our estimator of the value of a knock-in barrier option in 2 stages. A conditional expectation on jump times is used in stage 1. There, we use the simple fact that given jump times, the event that the jump-diusion process crosses the barrier in (0;T ] is equivalent to the union of the disjoint events that the jump-diusion process crosses the barrier in any of the inter-jump intervals. Next, in stage 2, we estimate the identity derived in stage 1 by a sequential conditional expectation. We also use control variate method in section 1.4. 3 In section 1.5, we use numerical examples to compare the eciency of our estimators to Metwally and Atiya's raw estimator. Also, based on our estimators introduced in section 1.3, 1.4 and the parity between vanilla call and barrier options, C=UOC+UIC, we propose a general rule for eciently pricing barrier options using our simulation procedures in the jump-diusion framework considered in this study. Numerical examples show that following this rule signicant amount of variance reduction and reduction in computing time is achieved which indicates that the analytical work and additional programming eort on this problem has been worthwhile. 1.2. Estimating the Value of a Knock-Out Option Using M&A's Raw Simulation Approach In section 1.2.3 we summarize Metwally and Atiya's raw simulation estimator of an up- and-out barrier call option. In [44] and also in our proposed simulation methods sample paths of the diusion process are generated sequentially based on jump times. This procedure is described in 1.2.1. Moreover, as will be seen in section 1.2.2, the probability that the maximum of a Brownian motion conditioned on its endpoints, i.e. a Brownian bridge, stays below a certain level (or crosses a certain level) is computable and used in both M&A's raw and our ecient simulation estimation of barrier options. 4 1.2.1 Sequential Simulation of the Jump-Diusion Sample Path based on Jump Times Suppose N(T ) = m > 0 and jump time are given. We would like to simulate S(t) sequentially at jump timesT 1 ;:::;T m . Note that from one jump to the other,S(t) evolves like a geometric Brownian motion with volatility and drift, so conditional onT 1 ;:::;T m S(T i ) =S(T i1 )e X i S(T i ) =S(T i )J i ; i = 1;:::;m (1.3) where T 0 = 0, S(T i ) denotes the price of the security right before the i th jump, X i 's are independent and X i N((T i T i1 ); 2 (T i T i1 )), where N(a;b 2 ) denotes the normal distribution with mean a and variance b 2 . So given jump times, sample path of the jump-diusion process can be generated at jump times sequentially by simulating fX i ;J i ; i = 1;:::;mg. 1.2.2 Maximum of a Brownian Bridge The probability that the maximum of a Brownian bridge remains below the barrier (or crosses the barrier) can be computed explicitly (see [44]). Here, we just introduce the notation and state the result. Let B(t) be a Brownian motion in the interval [s;u] with drift parameter , volatility , and B(s) = a and B(u) = b. The probability that the maximum of B(t) stays below a certain level V in the interval [s;u] is 5 P ( sup stu B(t)<VjB(s) =a;B(u) =b) = 8 > > < > > : 1 exp 2 (Va)(Vb) (us) 2 if a;b<V 0 otherwise Given T i1 and T i , the preceding probability is used to compute i = 1 i P ( sup T i1 uT i S(u)<BjS(T i1 ) =a;S(T i ) =b) in our setting. In words, i 's are barrier-crossing probabilities by the jump-diusion process between two consecutive jumps. 1.2.3 Metwally and Atiya's approach to estimate UOC The raw estimator of UOC in a single simulation run is computed based on the following algorithm. Suppose N(T ) =m 0. Let i = 1, T 0 = 0, and T m+1 =T , 1. SimulateT i by generating inter-jump times, T i T i1 , which are exponential. Also simulate J i . 2. Simulate S(T i ) and S(T i ) as described in section 1.2.1. If either of them exceeds B, return 0 as the estimate of UOC. Otherwise, go to step 3. 3. Compute i and generate a random number,U i . IfU i > i , return 0 as the estimate of UOC. Otherwise, go to step 4. 4. Ifi<m + 2, seti =i + 1 and go to step 1. Otherwise, return e rT (S(T )K) + as the estimate of UOC. 6 1.3. Variance Reduction in Estimating the Value of an Up- and-Out Barrier Option LetRe rT (S(T )K) + 1f B >Tg. In section 1.3.2, we introduce our proposed ecient simulation estimator of UOC = E[R]. In all of our proposed simulation methods, part of the variance reduction gained compared to the raw simulation estimator is due to stratied sampling whereN(T ) is the stratication variable. In section 1.3.1, we describe the poisson stratication employed in our study. 1.3.1 Stratied Sampling Let P m P (N(T ) = m) and E m [R] E[RjN(T ) = m], the following basic identity is used in our proposed ecient simulation methods. E[R] = l X m=0 E m [R]P m +E[RjN(T )>l]P >l (1.4) whereP >l = 1 P l i=0 P i andl is chosen such that P (N(T )>l) is small. Let R m be the single-run simulation estimator ofE m [R], andR >l be the single-run simulation estimator of E[RjN(T )>l]. To estimate E[RjN(T )>l] in a simulation run, we rst generate the value ofN(T ) conditional on it exceedingl using the inverse transform method. Suppose the generated value is l 0 , then R >l =R l 0. In section 1.3.2, we illustrate how to compute R m and R >l . Let b R be the estimator of E[R] based on n simulation runs. Let n m be the number of simulation runs allocated to estimating R m , and n >l be the number of 7 simulation runs allocated to estimating R >l , where P l m=0 n m +n >l =n. The stratied sampling estimator of E[R] based on n simulation runs is b R = l m=1 R m P m + R >l P >l (1.5) where R m and R >l are averages of R m and R >l based on n m and n >l simulation runs, respectively. While proportional sampling, i.e. setting n m = nP m and n >l = nP >l guaranties variance reduction over the raw simulation estimator of E[R], to specify the optimal value of n >l and n m 's, we rst perform a small simulation study to estimate Var(R m ) and Var(R >l ), then we set up and solve a simple optimization problem based on those estimates, (see [50] chapter 8 for details). 1.3.2 The Proposed Simulation Estimator of E[RjN(T ) =m] We use stratication on N(T ) as described in section 1.3.1 to estimate UOC = E[R]. Consider the identity (1.4) and supposeN(T ) =m. In what follows, we illustrate how to estimate E m [R]E[RjN(T ) =m] with R m . We rst consider S(T ) as a stratication- type variable and take the strata to be the disjoint intervals (K;B) and its complement. That is, E m [R] =E m [RjK TjT) = m+1 Y i=1 P ( B >T i j B >T i1 ) where T 0 0, T m+1 T , and for i = 1 the probability inside Q is simply P ( B > T 1 ). Conditional on jump times ands T which denotes a given value ofS(T ) that lies in (K;B), the expectation on the right of (1.6) can be written as E m [RjK T i j B >T i1 ;s T ) (1.7) where fori = 1 the probability inside Q is simplyP ( B >T 1 js T ). Let S i (S(0);S(T 1 ); :::;S(T i )). Now, consider thei th probability inside Q above and and note that conditional on S i , P ( B >T i jf B >T i1 ;s T g; S i ) = 8 > > < > > : i if S(T i ) and S(T i ) > < > > : e rT (s T K) Q m+1 i=1 i if I = 0 0 otherwise ThatR m;s T is an unbiased estimator follows from noting that two conditional expectations have been used to derive R m;s T . We rst condition E m [RjK T i j B >T i1 ;s T ) we use (1.8) by sequentially conditioning on S i . Note that it remains to show how to sequentially generate sample path of the jump- diusion process conditional on s T , i.e., a given value of S(T ). Recall (1.3) and let X m+1 X i=1 X i N(T; 2 T ) whereX i N((T i T i1 ); 2 (T i T i1 )); T 0 = 0; T m+1 =T . Note thatS(T ) =J(T )e X and based on (1.3), it suces to show how to sequentially simulate X 1 conditional on X =x, X 2 conditional on X 1 =x 1 and X =x, and in general ~ X i X i jX 1 =x 1 ;:::;X i1 =x i1 ;X =x (1.9) 10 In other words, we can write, ~ S(T i ) = ~ S(T i1 )e ~ X i ~ S(T i ) = ~ S(T i )J i ; i = 1;:::;m (1.10) where ~ S(T i ) denotes thatS(T i ) is to be simulated conditional on a given value ofS(T ). To illustrate how to sample ~ X i 's, we use lemma 1 below to conclude that ~ X i 's are normally distributed as follows ~ X i N 0 @ (x i1 X j=1 x j ) T i T i1 TT i1 ; 2 (T i T i1 )(TT i ) TT i1 1 A (1.11) Lemma 1 LetY i N( i ; 2 i ); i = 1;:::;n. Assume thatY i 's are independent. Also, let S = P n i=1 Y i . Then, ~ Y i (Y i jY 1 =y 1 ;:::;Y i1 =y i1 ;S =s)N i + 2 i (a P n j=i j ) P n j=i 2 j ; 2 i ( P n j=i+1 2 j ) P n j=i 2 j ! (1.12) where ~ Y i denotes Y i conditional on Y 1 =y 1 ;:::;Y i1 =y i1 ;S =s and a =s P i1 j=1 y j . Proof Let Y = P n j=i+1 Y j and note that the event (Y 1 = y 1 ;:::;Y i1 = y i1 ;S = s) is equal to the event (Y i +Y = s P i1 j=1 y j = a). So, ~ Y i (Y i jY i +Y = a). That is, to specify the distribution of ~ Y i dened in (1.12), we derive the distribution ofY i conditional onY i +Y =a, whereY i andY are independent normal random variables with given mean and variance. One simple way to specify the distribution of ~ Y i = (Y i jY i +Y = a) and obtain (1.12) is by considering the normal random vector (Y i ;Y i +Y ) with given mean and 11 covariance matrix and using conditional formula for the multivariate normal distribution, (see, for example, [22] page 65). So far we have shown how to estimate E m [RjK < S(T ) < B] with R m;s T derived in the previous page. Consider (1.6) again, it remains to show how to estimate P (K < S(T ) < BjN(T ) = m). Let m be the estimator of P (K < S(T ) < BjN(T ) = m), and note that m =P (K l P >l where, as mentioned in section 1.3.1, R m and R >l are the averages of the n m and n >l simulation runs allocated to estimating E m [R] and E[RjN(T ) > l], respectively, and P l m=0 n m +n >l =n. Our proposed simulation procedure is summarized algorithmically as follows. Recall that when N(T )>l we begin the simulation by generating the value of N(T ) conditional on it exceeding l. Let i = 1 and suppose N(T ) =m, 1. Simulate jump sizes, J. Simulate S(T )jK < S(T ) < B. Let s T be that simulated value. 12 2. Simulate jump times, T, by generating and then sorting m [0;T ] uniform random variables. 3. Simulate ~ S(T i ) and ~ S(T i ), i = 1;:::;m, sequentially using (1.10) and (1.11) based on lemma 1. If either of them exceedsB, stop and return 0 as the estimate of UOC. Otherwise, compute i ; i = 1;:::;m + 1. 4. Compute m based on (13) and return e rT (s T K)( Q m+1 i=1 i ) m as the estimator of UOC. Remark To simulate Xju 1 < X c where Z is a standard normal random variable and c is a positive constant. This procedure is based on the acceptance-rejection technique that uses the density function of an exponential random variable. 1.4. Variance Reduction in Estimating the Value of an Up- and-In Barrier Option Let R e rT (S(T )K) + 1f B Tg. In this section, we illustrate how to eciently estimate the expected discounted payo of an up-and-in barrier call option, UIC =E[R]. Stratied sampling withN(T ) as the stratication variable is used as mentioned in section 1.3.1. That is, similar to our approach in section 1.3.1, we use the following identity E[R] = l X m=0 E m [R]P m +E[RjN(T )>l]P >l (1.14) 13 to estimateE[R]. Moreover, sinceR is an increasing function of Q N(T) i=1 J i , we useJ(m) Q m i=1 J i as the control variate for the simulation estimator ofE m [R]. Let J(m) =E[J] m , we estimate E[R] by ^ R = l m=0 R m + ^ c m ( J(m) J(m) ) P m + R >l + ^ c >l ( J(>l) J(>l) ) P >l (1.15) where J(m), J(>l), R m , and R >l are averages of J(m), J(>l), R m , and R >l , based on n m and n >l simulation runs. Also, R m and R >l are the single-run unbiased estimators of E m [R] and E[RjN(T ) > l], respectively. Recall that to estimate E[RjN(T ) > l] in a simulation run, we rst generate the value ofN(T ) conditional on it exceedingl . Suppose the generated value is l 0 , then R >l = R l 0, J(> l) = J(l 0 ), and ^ c >l = ^ c l 0. Also, ^ c m is an approximate of c m =Cov(R m ;J(m))=Var(J(m)), (see chapter 8 of [50] on the use of control variates and on how to compute ^ c m , for instance, using simple linear regression models). 1.4.1 The Proposed Simulation Estimator of E[RjN(T ) =m] Suppose N(T ) = m 0, we derive our estimator of E m [R] E[RjN(T ) = m] in two stages where in stage one we consider the expected discounted payo of an up-and-in barrier call option conditional on jump times, T (T 1 ;:::;T m ), and derive (1.16). In stage two we use a sequential conditional expectation to estimate E m [RjT] based on (1.17)-(1.19). 14 Recall Re rT (S(T )K) + 1f B Tg. Suppose that T is given and let R i e rT (S(T )K) + 1f B 2 (T i1 ;T i ]g wherei = 1;:::;m + 1; andT 0 = 0;T m+1 =T . In words, given jump times, one can think ofR i as the return of an up-and-in barrier option which becomes alive only if the barrier is crossed in the i th inter-jump interval. Given T,f B 2 (T i1 ;T i ]g are disjoint events, and so we have E[1f B TgjT] =E[1f m+1 [ i=1 B 2 (T i1 ;T i ]gjT] = m+1 X i=1 E[1f B 2 (T i1 ;T i ]gjT] Using the preceding identity we can write E m [RjT] = m+1 X i=1 E m [R i jT] (1.16) Now, we illustrate how to estimate E m [RjT] using the preceding identity by sequentially estimating E m [R i jT]. Let S i (S(0);S(T 1 );::;S(T i )) and J = (J 1 ;:::;J m ). Consider E m [R i jT] conditional on jump sizes and S i , i.e., E m [R i jT; J; S i ]. We consider B as a stratication-type variable for the preceding conditional expectation with strata being (T i1 ;T i ] and its complement. That is, E m (R i jT; J; S i ) =E m (R i jT; J; S i ; B 2 (T i1 ;T i ])P ( B 2 (T i1 ;T i ]jT; J; S i ) (1.17) 15 Now, letI = minf1i m + 1 :S(T i ) orS(T i )Bg if no suchI exists, setI = 0. In words, given the jump times and jump sizes, I is the rst time that the jump-diusion process crosses the barrier right before or at a jump time.Observe that the probability on the right above can be computed as follows P ( B 2 (T i1 ;T i ]jT; J; S i ) = 8 > > > > > > > > > > > > > > < > > > > > > > > > > > > > > : ( Q i1 j=1 j ) i if I >i or I = 0 Q i1 j=1 j if I =i 0 if I > > > > > < > > > > > > : P I1 k=1 ( Q k1 l=1 l ) k B k + ( Q I1 l=1 l )B I if I6= 0 P m+1 k=1 ( Q k1 l=1 l ) k B k otherwise where Q k1 l=1 l is dened to equal 1 when k = 1, and R m = B 1 when I = 1. Recall the poisson stratication and control variate mentioned at the beginning of section 1.4, our proposed n-simulation-run estimator of UIC is ^ R dened in (1.15). The following algo- rithm summarizes the preceding procedure to estimateE m [R] based on a single simulation run. 1. Simulate T and J as mentioned in the previous section. 2. Sequentially simulate S(T i ) and S(T i ) up to I. That is, stop the sequential simu- lation when the barrier is crossed. Otherwise, simulate the whole sample path up to S(T ). 3. Compute inter-jump barrier-crossing probabilities. That is, i ; i = 1;:::;I1 when I6= 0 and 1 ;::: m+1 , otherwise. 17 4. R m as dened above is the single-run simulation estimator of UIC. Remark Instead of control variate J(m) in (1.15), one could use Antithetic method when estimating E m [R]. To implement the idea of Antithetic method in our set up, one would use m independent random numbers, U 1 ;:::;U m and 1U 1 ;:::1U m at two consecutive simulation runs to simulate (J i 1 ;:::;J i m ) and (J i+1 1 ;:::;J i+1 m ) whereJ i refers to theJ generated ati th simulation run (as opposed to simulating 2m independent random numbers in two consecutive runs). Since R is a monotone function of the jump process, J(T ), variance reduction is guaranteed by the corollary of theorem 8.10 in [50]. However, based on various numerical examples, we found that more variance reduction is gained in our set up by using the control variate method, i.e., using J(m) as the control variable. 1.5. Measuring the Eciency of the Proposed Estimators with Numerical Examples In this section \Method 1" and \Method 2" refer to our proposed ecient simulation procedures to estimate UOC and UIC introduced in section 1.3 and 1.4, respectively. In the following remark we introduce our ecient estimator of the value of the vanilla call option, C, in our set up. Given the Monte Carlo estimator of the vanilla call op- tion, the parity between vanilla call and barrier options, C=UOC+UIC, implies that in Monte Carlo estimation of barrier call options both Method 1 and Method 2 can be em- ployed. That is, one could use Method 1's (Method 2's) estimator \directly" to estimate UOC (UIC). Or, another alternative is to use our proposed estimator of C and Method 2 18 (Method 1) \indirectly" to estimate CUIC=UOC (C-UOC=UIC). In addition to com- paring the eciency of Method 1 developed in section 1.3 to that of Metwally and Atiya's raw estimator, we are also interested in devising a general rule for pricing barrier options by comparing the \direct" and \indirect" use of Method 1 and Method 2 in estimating both UOC and UIC. Monte Carlo estimation of UOC and UIC are considered in section 1.5.1 and 1.5.2, respectively. In 1.5.1 and 1.5.2 we present 6 of various numerical examples used to compare the eciency of our proposed estimators. As stated in section 1.1, e- ciency criteria considered here are variance of the estimators and the time to complete the simulation experiments (CPU time in MATLAB). The numerical examples consistently lead to the following conclusion. In Monte Carlo estimation of knock-out (knock-in) options we suggest using Method 2 indirectly (directly), when the price of the knock-out barrier option is very close to the price of vanilla call option. Otherwise, we suggest using \Method 1" directly (indirectly) . Remark: Ecient Estimator of C LetR =e rT E[(S(T )K) + ], to estimateE[R] we use the poisson stratication and control variate as we did in section 1.4. That is our n-simulation run estimator of E[R] is ^ R = l m=0 R m + ^ c m ( J(m) J(m) ) P m + R >l + ^ c >l ( J(>l) J(>l) ) P >l 19 Now, it remains to introduce the estimator of E[RjN(T ) = m]. As we mentioned in section 1.4, it is well-known that given the value of J(T ), E[e rT (S(T )K) + jJ(T )] is analytically computable by a formula similar to Black-Scholes formula for risk neutral pricing of a call option under geometric Brownian motion. That is, the preceding con- ditional expectation is in fact equivalent to the expected discounted payo of a vanilla call option with initial price S(0)J(T ), expiration time T , and strike price K, where the underlying security price follows a geometric Brownian motion with volatility and drift . So, to estimate C using Monte Carlo, we employ 3 well-known variance reduction techniques: stratication, control variate, and conditional expectation. Similar to the remark of section 1.4.1, instead of control variate, one could use Antithetic method when estimating C. However, in our set up, using numerical examples, we found that control variate leads to more variance reduction. 1.5.1 Monte Carlo Estimation of UOC In this subsection we compare three methods to estimate UOC in terms of variance of the estimators and computing time. As mentioned before we refer to Metwally and Atiya's approach introduced in [44] and summarized in 1.2.3 as raw simulation since no variance reduction has been introduced there. Method 1 is our proposed method of estimating UOC introduced in section 1.3. Finally, we have also used Method 2, introduced in section 1.4, combined with our proposed method of estimating C to estimate C-UIC=UOC. 20 Remark 1 Let R e rT (S(T )K) + and R i e rT (S(T )K) + 1f B Tg. We suggest beneting from the idea of common random numbers (see chapter 8 of [50]) when using our ecient estimator of C and Method 2 to estimate CUIC =UOC. We do so by using the same simulated jump sizes in each simulation run to estimate C and UIC as opposed to simulating independent jump sizes. Note that both R and R i are increasing functions ofJ(t), and so using the same stream of generatedJ's in each run to estimate C and UIC induces positive correlation betweenR andR i which leads to variance reduction. To see this suppose N(T ) =m. LetR m andR i m denote single-run simulation estimators of E[e rT (S(T )K) + jN(T ) = m] and E[e rT (S(T )K) + 1f B TgjN(T ) = m], respectively. Using the same simulatedfJ 1 ;:::;J m g for both R m and R i m , we have Var(R m R i m ) = Var(R m ) + Var(R i m ) Cov(R m ;R i m ) which implies variance reduction as opposed to using two independent streams of J's for R m and R i m . Note that using the aforementioned idea of common ranom numbers, in addition to variance reduction, computing time will also be reduced since in a simulation run with m jumps, m jump sizes are generated instead of 2m. Numerical Examples Suppose risk-neutral pricing is under consideration, and so the drift of the geometric Brownian motion S (t) satises = r 2 =2 +E[J]. Also, suppose that jump sizes are uniformly distributed in [:5; 1:5]. Let r = :08, T = 1 year, K = 34, S(0) = 30, = 1 and = :2. The following numerical examples are based on 10 5 simulation runs. Variance of the estimators stated in the second column of each table are based on 10 5 simulations runs. By b C, [ UOC, and d UIC we mean 10 5 -simulation run 21 estimates of C,UOC, and UIC, respectively. Also, to estimate C, we use the simulation procedure described in the rst remark of section 1.5. The CPU time in the following examples is in second. Also, CPU time for the main simulation studies to estimate barrier options and the initial simulation studies, which are required when using poisson strati- cation and control variate method (see the following remark), are reported separately. In the the following examples we refer to the computing time required for main simulation studies and initial simulation studies as \Main CPU Time" and \Initial CPU Time", respectively. Remark 2: Initializations Recall that poisson stratication has been used in all our proposed simulation procedures. When we explained the poisson stratication used in our study, section 1.3.1, it was mentioned that to outperform proportional sampling, we specify the optimal number of n m and n >l . Recall that n m and n >l are number of simulation runs allocated to estimate E m [R] and E[RjN(T ) > l], respectively, where P l m=0 n m +n >l =n. To specify optimal values ofn m andn >l , we perform a small initial simulation study which consists of specifying l, estimating Var(R m ) and Var(R >l ), and solving a simple optimization problem. In our numerical examples estimates of Var(R m ) and Var(R >l ) are based on 1000 simulation runs. Consider, for instance, Method 1 in Example 1 below. In the following table the rst row represents the number of jumps, the second row, number of simulation runs allocated to estimate E m [R]. The third row is simply P (N(T ) =m). Note that l = 8 and P >8 = 10 6 . Recall that in developing Method 2's estimator, section 1.4, and our estimator of C, in addition to poisson stratication we also useJ(m) = Q m i=1 J i as the control variable in the 22 m 0 1 2 3 4 5 6 7 8 > 8 n m 29487 43492 19664 5693 1377 244 32 6 3 1 P m .3679 .3679 .1839 .0613 .0153 .0031 .0005 7:3 10 5 9 10 6 10 6 Table 1.1: Optimal Stratication stratum whereN(T ) =m (see (1.15)). The application of control variate method requires an initial simulation study to estimate c m =Cov(R m ;J(m))=Var(J(m)). In our study, the link between regression and control variates is used to estimate c m 's. Our estimates of c m 's are based on 1000 simulation runs. So, when using Method 1's estimator, the initial simulation study species optimal values of n m and n >l . When using Method 2's estimator and our estimator of C, in addition to optimal values of n m andn > , the initial simulation study species c m 's. [ UOC Var( [ UOC) Main CPU Time Initial CPU Time M&A approach .7367 4:27 10 5 318:7969 s - Method 1 .7355 5.210 6 119:1719 s 16:3594 s Method 2 (C-UIC) .7269 1:36 10 4 198:0812 s 19:625 s Table 1.2: Example 1. 10 5 simulation runs, B=45, b C = 3:4774. \Method 1" is the most ecient method. [ UOC Var( [ UOC) Main CPU Time Initial CPU Time M&A approach 2.3461 2:6303 10 4 365:7969 s - Method 1 2.3369 8:7939 10 5 196:5781 s 24:5469 s Method 2 (C-UIC) 2.3210 1:8541 10 4 258:2500 s 30:0156 s Table 1.3: Example 2. 10 5 simulation runs, B=60, b C = 3:4822. \Method 1" is the most ecient method. Remark 3 In addition to signicant variance reduction gained by using our proposed estimators, it is very interesting to note that computing time of our ecient estimators is considerably lower than computing time of the raw simulation estimator. This signicant 23 [ UOC Var( [ UOC) Main CPU Time Initial CPU Time M&A approach 3.3302 4:3425 10 4 466:3594 s - Method 1 3.3269 2:3838 10 4 213:5156 s 26:049 s Method 2 (C-UIC) 3.3188 1:1243 10 4 296:3125 s 32:5938 s Table 1.4: Example 3. 10 5 simulation runs, B=90, b C = 3:4435. [ UOC is very close to b C. \Method 2" combined with our estimator of C is the most ecient method. amount of \time reduction" is due to poisson stratication. Recall that (see 1.2.3) in each single run of the raw simulation approach, one is to generate exponential random variables and stop when their sum exceeds the expiration time, T . When using poisson stratication (see 1.3.1), however, there is no need to generate exponentials in each single run as a small simulation study pre-determines the number of runs conditional on specied values ofN(T ) (except for the runs allocated to estimatingE[RjN(T )>l] where the value of N(T ) conditional on it exceeding l is generated using inverse transform method). 1.5.2 Monte Carlo Estimation of UIC In the previous subsection it was shown that our proposed estimators of UOC lead to sig- nicant variance reduction and also reduction in computing time compared to Metwally and Atiya's raw estimator. Here, we present 2 numerical examples to compare the e- ciency of Method 2 which \directly" estimates UIC with that of Method 1 combined with our proposed estimator of C to estimate CUOC=UIC. As stated before, our simulation experiments suggest that Method 2 is the most ecient method when UOC is very close to C. Otherwise, we suggest using Method 1 \indirectly" to estimate UIC. Remark Let R o e rT (S(T )K) + 1f B > Tg. When using ecient estimator of C and Method 1's estimator to estimate CUOC=UIC, the idea of common random 24 numbers as mentioned in remark 1 in section 1.5.1 would not lead to variance reduction. This is so because R o is not a monotone function of J(t), and a positive correlation can not be induced between R and R o by using the same stream of jump sizes in each simulation run. Numerical examples We use the same numerical example used in section 1.5.1 except that we chose = 2, =:5 in the following 3 examples. d UIC Var( d UIC) Main CPU Time Initial CPU Time Method 2 4.4482 19 10 4 89:6875 s 21:4219 s Method 1 (C-UOC) 4.4374 3:24 10 4 119:9531 s 22:1406 s Table 1.5: Example 1. 10 5 simulation runs, B=75, b C = 7:1299. \Method 1" combined with our estimator of C is the most ecient method. d UIC Var( d UIC) Main CPU Time Initial CPU Time Method 2 2.5614 16 10 4 125:3438 s 25:7031 s Method 1 (C-UOC) 2.5478 7:293 10 4 151:0946 s 23:8125 s Table 1.6: Example 2. 10 5 simulation runs, B=100, b C = 7:1316. \Method 1" combined with our estimator of C is the most ecient method. d UIC Var( d UIC) Main CPU Time Initial CPU Time Method 2 .7574 6:61 10 4 174:5313 s 37:0469 s Method 1 (C-UOC) .7623 12 10 4 161:9375 s 26:7813 s Table 1.7: Example 3. 10 5 simulation runs, B=150, b C = 7:0964. [ UOC is very close to b C. \Method 2" is the most ecient method. 25 Chapter 2 Ecient Simulation of a Random Knockout Tournament 2.1. Introduction Selecting a single winner from a set of many contestants is the common theme in many sporting events. Tournaments of varying structures are used for choosing a winner in situations where players can be compared pairwise. In a knockout tournament players compete head-to-head in matches; with the losers being eliminated from the tournament. In a random knockout tournament at the outset of each round players are paired at random. A general knockout tournament involving N players is specied by positive integer parameters n;m 1 ;:::;m n1 ;m n = 1; satisfying 2m i N i1 X j=1 m j ; i = 1;:::;n 1 n X j=1 m j = N 1 26 The interpretation of these parameters is that the tournament is to consist of n rounds, with round i involving m i matches. Let r k = N P k1 j=1 m j be the number of players remaining at the beginning of round k. We say that the knockout tournament is a random knockout tournament if the players involved in roundi are randomly chosen, and then randomly paired, from theN P i1 j=1 m j players who have not been eliminated in any of the preceding rounds. We assume that the outcome of a match between two players in any round depends only on the pair and not on the round or previous rounds, and we let P ij = 1P ji denotes the probability that player i defeats player j if they are paired. The player who wins the match in round n is called the winner of the tournament. LettingP i be the probability that playeri is the winner of the tournament i = 1;:::;N; we are interested in using simulation to eciently estimate these P i : In section 2.2 we propose an ecient simulation procedure based on the notion of \observed survivals". In section 2.3 we combine conditional expectation and poststrat- ication for another simulation approach. In 2.4, we consider the \classical knockout tournament" where N = 2 n and round k consists of 2 nk matches, and present a second conditional expectation estimator. Two numerical examples are considered in section 2.5. In these examples, a considerable variance reduction over the raw simulation estimator is gained when using any of our suggested methods. The time to complete the simula- tion experiment is also considered for each method. Considering the variance reduction criterion in parallel with the time to obtain the simulation results, we conclude that the observed survival method is best. 27 2.2. Method 1: Observed Survivals and the Product Esti- mator We can model the evolution of the random knockout tournament from each round to the next one by the transitions of a Markov chain whose states are the sets of players alive at the outset of each round. Let N i be the total number of rounds that player i survives (that is, for 0j <n,N i =j if playeri loses in roundj + 1:) Our problem of estimating P i =P (N i >n1) is in fact a special case of a more general problem of using simulation to estimate the distribution of the number of transitions it takes a Markov chain to enter a particular set of states. Specically,N i +1 is the number of transitions that it takes the Markov chain whose state is the set of surviving players at the beginning of each round to enter a state that does not contain i. In [53] simulation estimators based on observed hazards were developed to eciently estimate rst passage distributions to a given set of states of a Markov chain. Our method 1 uses the estimator developed in [53]. Each simulation run involves simulating the tournament throughn1 rounds. (There is no need to simulate round n.) Observe that P i =P (N i >n 1) = Q n k=1 P (N i >k 1jN i k 1) Note that s k (i) P (N i > k 1jN i k 1) is the discrete survival rate value representing the probability that player i survives roundk given that he has survived up to roundk. Using simulation, we illustrate how to estimate P (N i >n 1) = Q n k=1 s k (i). 28 Remark s 1 (i) = 1 2m 1 N P l6=i P li =(N 1) Let N j i be the number of rounds that i survives in the j'th simulation run. We dene \observed survivals" at the j'th simulation run as point estimates of s k (i) as follows. Let j k (i) =P (N j i >k 1j j k1 ) where j k is the simulated set of players present at the beginning of round k. Thus, ifN j i is the number of rounds that i survives in the j th simulation run, then j k (i) = 8 > > < > > : 1 2m k r k (r k 1) P l2 j k ;l6=i P li ifN j i k 1 0 ifN j i <k 1 Note that j k (i) is an unbiased estimator of s k (i) provided that N j i k 1: Suppose that we make a total of r simulation runs. Let k (i) = P r j=1 j k (i)=f]j :N j i k 1g Then k (i) is an estimate of the discrete survival rate values k (i) =P (N i >k1jN i k 1), and Q n k=1 k (i) is a consistent estimate of P (N i > n 1) = Q n k=1 s k as r, the number of simulation runs, goes to innity. Remark Suppose performingr simulation runs, we have not been able to obtain point estimates of s m (i);:::;s n (i). In other words, consider the situation where N j i <m 1< n1 for allj = 1;::;r. To estimateP (N i >n1) based on the aforementioned estimator it is useful to derive analytic bounds for s k (i). 29 Consider the set of probabilities for player i,fP ij :j6=i ;j = 1; 2;:::;Ng, and order them so that P iI 1 P iI 2 :::P iI N1 , lower and upper bounds for s k (i) are as follows 1 2m k r k (r k 1) P r k 1 j=1 P I j i s k (i) 1 2m k r k (r k 1) P N1 j=Nr k +1 P I j i Remark In [53] a two-step method to approximate Var( Q n k=1 k ) is developed. First Var( k ) ; k = 1;:::;n is approximated, and then using the Delta method and combining these approximations, an approximation of Var( Q n k=1 k ) is obtained. See [53] for details. 2.3. Method 2: Conditional Expectation and Post Strati- cation Until player i loses a match, player i is designated as the surrogate for i. If a player beats i in some round, then at the beginning of the next round that player becomes the surrogate for i. Any player beating a surrogate for i during a round takes over as the surrogate for i at the beginning of the next round. In this manner, at the beginning of each round exactly one of the remaining players is the surrogate for i. Let T i be the number of matches played by players while they are the surrogate of i. Suppose that T i =t. Let F i;1 denote the rst player to play i in a match; let W i;1 be the winner of that match. Let F i;2 denote the next player to play W i;1 in a match; let W i;2 be the winner of that match. LetF i;3 denote the next player to playW i;2 in a match; letW i;3 be the winner of that match, and so on up to F i;t : Consequently, in order to win the tournament, player i would have to successively beat all the players F i;j ;j = 1;:::;t: 30 Now, let p j = 2m j r j ; j = 1;:::;n denote the proportion of the players that have survived the rst j 1 rounds who play in round j. Also, let I j ;j = 1;:::;n; be independent Bernoulli random variables such that P (I j = 1) =p j ; j = 1;:::;n It is easy to see that P (T i =k) =P ( n X j=1 I j =k) To compute the preceding probabilities, let P s (k) =P ( s X j=1 I j =k); sn; ks Conditioning on I s gives that P s (k) =p s P s1 (k 1) + (1p s )P s1 (k) Starting with P 1 (1) =p 1 ; P 1 (0) = 1p 1 ; we can recursively solve the preceding until we have the quantities P n (k); k = 0;:::;n: Let I i be the indicator of the event that i wins the tournament. Because E[I i jT i =t;F i;1 ;:::;F i;t ] = t Y j=1 P i;F i;j 31 we could use Q t j=1 P i;F i;j as a conditional expectation estimator of P i . However, be- cause there would appear to be a strong negative correlation between I i and T i , a post stratication estimator should have an even smaller variance. That is, write P i = n X k=1 E[I i jT i =k]P n (k) and then, after completingr simulation runs, estimateE[I i jT i =k] by the average of the quantities Q T i j=1 P i;F i;j obtained in those runs havingT i =k. Using Delta method, it can be shown that the asymptotic variance of the post-stratied estimator (i.e. the variance as the number of simulation runs goes to innity) is equal to that of stratied estimator, see e.g. [22] page 235. In the large sample limit, we can also estimate the variance of the estimator of E[I i jT i =k] with the sample variance of the observations falling in the k'th stratum divided by P n (k). Remark: For a xed value of i, there may be values of k for which there are no runs for which T i =k. If the value of P n (k) is suciently small, we can use 0 as the estimate of E[I i jT i = k]P n (k): However, if there are enough such values of k; the preceding post stratication will not work for estimating P i . In this case, we suggest using the simu- lated data to estimate P i by using the conditional expectation estimator Q T i j=1 P i;F i;j in conjunction with using T i as a control variable with known mean E[T i ] = P n j=1 p j : (Of course, we would still use the post stratication procedure to estimate the values P j for which the preceding diculty did not arise.) 32 Example To get an idea of the amount of variance reduction the preceding methods yield, suppose we want to estimate P i for a specied i when P ij = ; j6= i: Now, with I j ;j = 1;:::;n; as previously dened P i = E[ T i ] = E[ P n j=1 I j ] = n Y j=1 E[ I j ] = n Y j=1 (p j + 1p j ) = n Y j=1 (1 (1)2m j =r j ) The straight conditional expectation estimator reduces, in this case, to the estimator T i . Its variance is Var( T i ) = E[( T i ) 2 ]P 2 i = E[( 2 ) T i ]P 2 i = n Y j=1 (1 (1 2 )2m j =r j )P 2 i For instance, when N = 5; m j 1 ; = 1=2, the preceding gives that Var((1=2) T i ) = 4 Y j=1 (1 3 2(6j) ) (1=5) 2 :015 whereas the variance of the raw simulation estimator is P i (1P i ) =:16: 33 Both the conditional expectation stratication estimator and the previous section's estimator based on survival probabilities have 0 variance when P ij : 2.4. Classical Knockout Tournament In a classical knockout tournament there are 2 n players,N = 2 n ,n rounds, and in thek'th round of the tournament, there are 2 nk matches,m k = 2 nk ,r k = 2 nk+1 . In a classical random knockout tournament at the outset of each round players are paired at random. In performing a simulation run of the classical 2 n -player random knockout tournament we can determine pairings in each round by generating a single random permutation R 1 ;:::;R 2 n of 1;:::; 2 n that would determine how to assign players to starting positions that are numbered from 1 to 2 n . Matches in the rst round are specied by pairing players assigned to adjacent starting positions. In the k'th round, adjacent winners of the previous round are paired,k = 2;:::;n. Consider, for instance, an eight-player random knockout tournament with, [1; 2; 3; 4; 5; 6; 7; 8], as the generated random permutation. In the rst round, player 1 plays against player 2, player 3 plays against player 4, player 5 plays against player 6, and player 7 plays against player 8. In the second round, the winner of player 1 and 2 plays against the winner of player 3 and 4, the winner of player 5 and 6 plays against the winner of player 7 and 8. In the last round survivors of the second round play against each other. A recursive procedure has been introduced in [15] to compute win probabilities for a classical random knockout tournament when players' starting positions are known and from that point on pairings are determined as described 34 above. Thus, given the random permutation R, we can compute the conditional win probabilities, P i ; i = 1;::; 2 n . Remark When compared to method 2, the second conditional expectation estimator is derived by conditioning on less information and so it is preferable to method 2 in terms of the variance reduction criterion. 2.5. Numerical examples Consider a knockout tournament with 8 players. It is common to represent players's probabilities, P ij ; i = 1; 2;:::; 8 ; j = 1; 2;:::; 8, in a matrix which is usually called the \preference matrix". Consider the following preference matrix for an 8-player random knockout tournament: 0 B B B B B B B B B B B B B B B B B B B B B B B B B B B @ :2 :3 :5 :4 :7 :8 :6 :8 :4 :6 :7 :8 :5 :5 :7 :6 :2 :3 :9 :5 :3 :5 :4 :8 :7 :6 :1 :6 :6 :3 :7 :3 :8 :4 :5 :3 :2 :1 :4 :2 :6 :6 :2 :5 :5 :9 :6 :4 :3 :4 :5 :7 :4 :5 :4 :7 1 C C C C C C C C C C C C C C C C C C C C C C C C C C C A We rst consider a one-match-per-round random knockout tournament,N = 8,n = 7, m 1 =::: =m 7 = 1. A classical 8-player random knockout tournament,n = 3,m k = 2 nk , 35 is then considered for our empirical study. Using 10 5 simulation runs, we would like to estimate P i 's for both types of tournaments. The result The following rst three tables refer to the simulation result of the one- match-per-round random knockout tournament, and the last three tables refer to the classical random knockout tournament. (Note that in tables given below, numbers below each Pi refer to simulation estimates related to player i). Recall that method 1 is the method of observed survivals, and method 2 refers to the conditional expectation com- bined with post-stratication. In the simulation of the classical knockout tournament method 3 refers to the second conditional expectation introduced in section 2.4. Based on our numerical example, method 2 in the one-match-per-round tournament and both method 2 and 3 in the classical knockout tournament are preferable to the method of observed survivals in terms of variance reduction. However, in addition to the variance reduction, the factor of time has also been taken in to account by measuring the amount of time MATLAB takes to complete our simulation experiments. It is interesting to ob- serve that our estimator based on the notion of observed survivals takes considerably less time compared to the conditional expectation based estimators in both types of tour- naments. This can be explained by noting that when computing simulation estimators based on the notion of conditional expectation introduced in this study, there is more than one multiplication operation involved in a single simulation run. In contrast, in method 1, point estimates of s k (i) are rst collected for k = 1; 2;::;n over 10 5 simulation runs, and then averaged. P i is then estimated by the product of theses averages at the 36 end of the simulation experiment. Based on the values obtained in our numerical exam- ple, method 1 is the most ecient simulation procedure for estimating P i 's in random knockout tournaments. One Match Per Round Random Knockout Tournament P1 P2 P3 P4 P5 P6 P7 P8 Raw Simulation .1097 .2129 .1132 .1388 .1195 .0481 .1208 .1369 Method 1 .1112 .2112 .1141 .1394 .1184 .0483 .1207 .1365 Method 2 .1113 .2098 .1131 .1394 .1183 .0483 .1194 .137 Table 2.1: Estimates of P i 's P1 P2 P3 P4 P5 P6 P7 P8 Raw Simulation 977 1676 1004 1195 1052 458 1062 1182 Method 1 448.45 122.2 182.1 208.43 104 65.67 130.08 52.31 Method 2 59 111 63 84 55 21 68 62 Table 2.2: Variance of the estimators based on 10 5 simulation runs (10 9 scale) Raw Simulation Method 1 Method 2 101.926 184.9 7362 Table 2.3: Elapsed time (in second) 37 Classical Random Knockout Tournament P1 P2 P3 P4 P5 P6 P7 P8 Raw Simulation .1083 .2249 .1107 .1422 .1183 .0367 .1209 .138 Method 1 .1078 .2245 .1123 .1439 .1174 .0363 .1188 .1384 Method 2 .1082 .2244 .1124 .1437 .118 .0366 .1192 .1385 Method 3 .108 .2241 .1122 .1437 .1176 .0365 .1193 .1387 Table 2.4: Estimates of P i 's P1 P2 P3 P4 P5 P6 P7 P8 Raw Simulation 966 1743 985 1220 1043 353 1063 1189 Method 1 107 81 114 130 78 38 117 42 Method 2 51 60 68 92 40 10 51 21 Method 3 14 23 22 31 13 3 16 9 Table 2.5: Variance of the estimators based on 10 5 simulation runs (10 9 scale) Raw Simulation Method 1 Method 2 Method 3 9.1791 13.2288 2580.4 3475.3 Table 2.6: Elapsed time (in second) Remark To estimate the variance of Method 1's estimator in our empirical study, we have used the approximation result derived in [53] and mentioned in the last re- mark of section 2.2. We call it the internal variance estimator, Var(int). To study it's accuracy , we compared it with an external variance estimator Var(ext) based on a meta-experiment of 1000 independent experiments. From the aforementioned numerical example,Var(int) = 120 10 9 andVar(ext) = 118 10 9 , and so we conclude that the Var(int) is accurate enough to be used instead of Var(ext) in our empirical study. The external variance estimator is Var(ext) = (1=999) P n i=1 (X i X) 2 where X = Q n k=1 k 38 Chapter 3 Ecient Monte Carlo Estimation of System Reliability 3.1. Introduction Specifying the probability that a component system functions, i.e., system reliability (see Barlow and Proschan [4]) is often a computationally complex problem. It is well-known that Monte Carlo simulation is an alternative eective approach for estimating system reliability. We are interested in developing ecient simulation procedures to estimate system reliability. In our study, we have addressed two dierent cases: a system of binary components under a certain type of dependency, section 3.2, and a system of 3-state independent components, section 3.3. In [48], Ross introduces a new very ecient simulation approach to estimate expected values of functions of Bernoulli random variables under certain types of dependencies. The simulation approach introduced in [48] is based on an innovative use of stratied sampling. In our study, in both cases mentioned above, Ross' improved stratied sampling approach (ISS) plays a key role in deriving our ecient estimators of system reliability. For future references ISS is summarized in section 3.1.1. 39 Consider an n-component system. Let X i , a discrete-valued random variable, be the state of component i, i = 1; 2;:::;n. The random vector X = (X 1 ;:::;X n ) represents the state of an n-component system. Recall that a random vector (X 1 ;:::;X n ) is said to be exchangeable if for each permutation R 1 ;:::;R n off1; 2;:::;ng, (X 1 ;:::;X n ) and (X R 1 ;:::;X Rn ) have the same distribution. In the context of system reliability, a structure function, H = h(X), is a binary nondecreasing function dened on the n-dimensional vector X, and it determines whether the component system functions. As mentioned earlier, the main objective of our study is developing ecient simulation procedures to estimate the probability that an n-component system functions, P (h(X) = 1) = E[H]. Binary Dependent Components Very few system reliability models allow the com- ponents to be dependent. In practice, however, it is not unlikely that system components be dependent. In section 3.2, rst, a certain type of dependency among system compo- nents is dened. We then develop ecient simulation procedures to estimate the reliability of a system whose components are dependent according to our denition. Consider a sys- tem of binary components where the state of each component is dependent on the state of the rest of the system through the number of working components as stated below: P (X i = 1jX [i] = x [i] )P (X i = 1j X k6=i X k =j) = i (j) (3.1) where X [i] = (X 1 ;:::;X i1 ;X i+1 ;:::;X n ) and components of x [i] sum to j, and i = 1;:::;n. 40 We refer to i as the \conditional dependency function" of component i. It should be noted that binary components of a system with one conditional dependency function, = i ; i = 1; 2;:::;n , are exchangeable while capturing the dependency among sys- tem components with more than one conditional dependency function leaves the system components non-exchangeable. Three cases are studied in section 3.2. In case 1 binary dependent components are exchangeable, i = ; i = 1; 2;:::;n, section 3.2.2. In cases 2 and 3 binary dependent components are non-exchangeable, section 3.2.3 and 3.2.4, respectively. In cases 1 and 2, we illustrate how the notion of Markov Chain Monte Carlo sampling in general and that of random scan Gibbs sampler in particular will enable us to utilize ISS and develop an ecient simulation procedure to estimate the system reliability. In case 3 the probabilistic structure of the system model is dened so that improved stratied sampling loses its practicality. In this case, we suggest using Gibbs sampler combined with Antithetic method to estimate system reliability. 3-State Independent Components In section 3.3, we consider a system whose com- ponents are independent and can be in state 1, 2, or 3. Components' probability mass functions are known, P (X i = j) P i (j) ;i = 1; 2;::;n ; j = 1; 2; 3. To estimate the reliability of a system of this type we introduce a 3-stage simulation procedure. ISS as developed in [48], could only be employed when dealing with Bernoulli random variables or binary components in the context of system reliability. Our proposed 3-stage simula- tion procedure is in fact the generalization of the improved stratied sampling method from Bernoulli to 3-state independent random variables. 41 3.1.1 Improved Stratied Sampling Approach As mentioned earlier, system components (random variables) in [48] are considered to be binary (Bernoulli). When components are exchangeable the following stratied sampling- based algorithm is used to estimate system reliability, and when dealing with independent non-exchangeable components, the proposed simulation procedure in [48] is a 2-stage simulation method which is also summarized below. Exchangeable case Let X = (X 1 ;X 2 ;:::;X n ) be an exchangeable Bernoulli random vector and S = P n i=1 X i . Suppose that the probability mass function P j P (S = j), j = 0; 1;:::;n is computable. Leth be a non-decreasing function dened onn-dimensional binary vectors, h(X) = H. To use simulation to estimate E[H], S is used as the \stratication variable". That is, E[H] = n j=0 E[HjS =j]P (S =j) (3.2) Each simulation run of the improved stratied sampling approach performed as follows: 1. Generate R 1 ;:::;R n , a random permutation of 1; 2;:::;n. 2. For each j, j = 0; 1;:::;n, compute H j = h(x 1 ;x 2 ;:::;x n ) where x R i = 1 for i = 1; 2;::;j and x R i = 0 for i =j + 1;:::;n 3. Return the estimate ^ = j H j P j In [48], it is proved that ^ is an unbiased estimator of and has a smaller variance than does the raw simulation estimator. 42 Non-exchangeable case Consider binary independent non-exchangeable components. Let P (X i = 1) P i ; i = 1; 2;:::;n and choose P to be close to most of the P i 's. For instance it can be arithmetic mean. Let L be the set of component indices i such that P i P and let H be the remaining component indices. The simulation is performed in two stages, with the second stage using ISS approach to estimate the reliability of a system whose components all work with probability P . In stage 1; generate random numbers U 1 ;:::;U n ; and set J i = 8 > > > > > > < > > > > > > : 0 ifi2L andU i >P i =P 0 ifi2H andU i > (1P i )=(1P ) 1 otherwise If J i = 0 then set X i = 0 if i2 L; or set X i = 1 if i2 H: Now, consider the system conditional on these values of X i ; under the supposition that all the other components (i.e., those components i for which J i = 1) have reliability P . Now, use the ISS method when all P i =P to estimate the reliability of the resulting system. When most of the P i are close to P this method should work very well. 3.2. Binary Dependent Components The notion of random scan Gibbs sampler proves to be useful in all the following three cases where components are binary and the state of each component depends on the state of the remaining ones through the number of working components as dened in section 3.1: 43 P (X i = 1jX [i] = x [i] )P (X i = 1j P k6=i X k =j) = i (j) where X [i] = (X 1 ;:::;X i1 ;X i+1 ;:::;X n ) and components of x [i] sum to j, and i = 1;:::;n. To facilitate our analysis in sections 3.2.2 - 3.2.4 and introduce the notation, we rst provide a brief non-technical introduction to Gibbs sampler in section 3.2.1, (see, for example, chapter 6 of [39]). 3.2.1 Random Scan Gibbs Sampler While in traditional Markov Chain analysis given a transition rule the main interest lies in knowing what the stationary distribution is, in Markov Chain Monte Carlo simulation the stationary distribution is known and the objective is reaching this distribution by prescribing an ecient transition rule. A Markov process is evolved to achieve the sam- pling of the desirable target distribution. In Gibbs Sampler, a popular scheme of MCMC algorithms, the underlying Markov Chain is constructed by composing a sequence of con- ditional distributions along a set of directions. Let X t = (X t 1 ;:::;X t n ) be then-dimensional random vector simulated at iteration t. In Gibbs sampler, one randomly, random scan Gibbs sampler, or systematically, systematic-scan Gibbs sampler, chooses a coordinate, sayX t i , and updates it with a new sample X t+1 i , at iteration t + 1, drawn from the con- ditional distribution, P (X i jX [i] ) where X [i] refers to (X 1 ;::;X i1 ;X i+1 ;:::;X n ). The remaining components then will be leaved unchanged, X t+1 [i] = X t [i] . 44 3.2.2 Case 1: Exchangeable Components Consider a system of n binary components with one conditional dependency function, i = for all i = 1; 2;:::;n, P (X i = 1jX [i] = x [i] )P (X i = 1j k6=i X k =j) = j (3.3) where components of x [i] sum to j. Equation (3.3) implies components are exchangeable. Recall thath(X) is the structure function of the system under study, and we would like to estimate the system reliability, E[h(X)]. Based on the probabilistic structure of the system components in case 1, the improved stratied sampling approach, as introduced in section 3.1.1, can be used if the probability mass function of the number of working components can be computed analytically. In other words, given , we would like to specify p.m.f of S. That is, we would like to analytically compute: P ( n i=1 X i =j) =P (S =j) =P j ; j = 0; 1;:::;n (3.4) We now illustrate how the Markov Chain representation of S and notion of a Random Scan Gibbs sampler would enable us to specify the distribution of S. Consider a discrete time Markov Chain,fS t ; t = 0; 1; 2;:::g with the state spacef0; 1;:::;ng. Given , we specify the transition rule, P (S t+1 =jjS t =i) =P ij , of this Markov Chain as follows: P 00 = 1 0 ; P ii = i i1 + (ni)(1 i ) n ; P ii1 = i(1 i1 ) n ; P ii+1 = (ni) i n (3.5) 45 The preceding transition probabilities are obtained by evolving the underlying n-dimen- sional discrete time Markov Chain X t = (X t 1 ;X t 2 ;:::;X t n ) based on iterations of a random scan Gibbs sampler where knowingP (X i = 1j k6=i X k =j) = j makes the transition from iteration t to iteration t + 1 a standard random scan Gibbs sampler step for simulating the random vector X = (X 1 ;X 2 ;:::;X n ). Note thatfS t ; t = 0; 1; 2;:::g is an irreducible ergodic Markov Chain, so the limiting probabilities, j = lim t!1 P t ij ; j = 0; 1;:::;n, exist and can be computed according to the limiting probabilities equations. Based on (5) limiting probabilities equations are as follows, j = 1;:::;n 1. 0 = 0 (1 0 ) + 1 1 0 n j = j1 (nj + 1) j1 n + j j j1 + (nj)(1 j ) n + j+1 (j + 1)(1 j ) n n = n1 n1 n + n n1 n X j=0 j = 1 46 and so limiting probabilities are: P (S = 0) =P 0 = 0 = " 1 + n i=1 n i i Y k=1 k1 1 k1 # 1 (3.6) P (S =j) =P j = j = " n j j Y k=1 k1 1 k1 # P 0 (3.7) So, to estimate the reliability of a system whose binary components are exchangeable and dependent according to a given \conditional dependency function" (see (3.3)) we act as follows. We rst analytically compute P j 's based on (3.6) and (3.7). Next, the improved stratied sampling, see section 3.1.1, is employed to estimate the probability that the system functions. 3.2.3 Case 2: Non-exchangeable Components with a Small Number of Conditional Dependency Functions Suppose that dependency across system components is taken in to account by a small number of conditional dependency functions. As will be seen below, this assumption that the number of \conditional dependency functions" is small would make analytical computations, necessary to develop a modied version of improved stratied sampling, practical. Suppose that there are m conditional dependency functions, 1 ;:::; m and m < n. Consider m distinct sets of components, A 1 ;::;A m where we say X i 2 A j if j is the conditional dependency function of component i, i = 1; 2;:::;n ; j = 1; 2;:::;m. Let n j be the number of components whose conditional dependencies are dened by j , P m j=1 n j = n. Note that since the components are not exchangeable in this setting, the 47 improved stratied sampling approach as introduced in [48] and summarized in section 1:1 can not be used. We now illustrate how to modify ISS to make it applicable in this case. The idea is to consider subsets of exchangeable components. That is, considering subsets where each subset consists of components with common \conditional dependency function". Let the random variableN j be the number of working components in setA j . Consider the random vector, N = (N 1 ;:::;N m ), whereN j = 0;:::;n j . We use N as the stratication variable to estimate E[H]. That is, =E[H] = k E[HjN=k]P k (3.8) where h(X) =H, k = (k 1 ;:::;k m ), and P k =P (N=k). Components that belong to each set, A j , are exchangeable as their common conditional dependency function is j . Now, given 1 ;:::; m , we would like to specify the p.m.f of the mdimensional random vector N. That is, given 's, we would to analytically derive, P (N=k) =P ((N 1 ;:::;N m ) = (k 1 ;:::;k m )) (3.9) for all possible values of k. Similar to our approach in case 1, section 3.2.2, we dene and consider an m-dimensional Markov Chain , fN t ;t = 0; 1;:::g, with state space fN = (N j ;j = 1; 2;:::;m) ; N j = 0; 1;:::;n j g where the transition rule P (N t+1 = jjN t = i) =P ij (3.10) 48 would be completely specied by considering the underlyingn-dimensional Markov chain X = (X 1 ;X 2 ;:::;X n ), and evolving it according to iterations of a random scan Gibbs sampler where conditional dependency functions 1 ;:::; m are given. In other words, from each iteration to the next one, X is updated based on j if the randomly chosen coordinate belongs toA j ,j = 1; 2;:::;m. So the limiting probabilities, P k =P (N=k), of the \main" Markov Chain,fN t ;t = 0; 1;:::g, are analytically computable in a practical manner as long as Q m j=1 (1 +n j ) is a small number. After analytically computing P k as stated above, our method, which is a modied version of ISS, will, on each simulation run, estimate the quantities E[HjN = k] for all possible values of k = (k 1 ;:::;k m ). This is done by letting each simulation run be as follows: 1. Generate m dierent random permutations, R 1 ;:::; R m , where R j = (R 1 ;:::;R n j ) is a random permutation of 1; 2;:::;n j 2. For all possible values of k = (k 1 ;:::;k m ) compute H k = h(x 1 ;:::;x n ) by utilizing each R j in the following way. When X i 2 A j consider R j and set x R i = 1 for i = 1; 2;::;k j . Set x R i = 0 for i =k j + 1;:::;n j 3. Return the estimate ^ = k H k P k Note that when N = (N 1 ;:::;N m ) = (k 1 ;:::;k m ) each of the n j k j sets of k j components in set A j is equally likely to be the set of ones that work, j = 1; 2;:::;m. So H k 's are unbiased estimators of E[HjN], and consequently ^ is an unbiased estimator of . Now, let H =h(X) be the raw simulation estimator of . 49 Proposition Var(^ ) Var(H) Proof We show that ^ can be written as a conditional expectation of the raw simulation estimator. Suppose that X = (X 1 ;X 2 ;:::;X n ) is generated by rst generating the value of N, then generating m random permutations R 1 ;:::; R m , where R j = (R 1 ;:::;R n j ) is a random permutation of 1; 2;:::;n j , and setting x R i = 1 for i = 1; 2;::;k j and x R i = 0 for i =k j + 1;:::;n j for X i 2A j . The result follows because ^ =E[HjR 1 ;:::; R m ]. 3.2.4 Case 3: Non-exchangeable Components with a Large Number of Conditional Dependency Functions When the number of system components and the number of conditional dependency func- tions 1 ;:::; m is so that Q m j=1 (1 +n j ) becomes a large number, our proposed simulation procedure in case 2 loses its practicality. In this case, based on the way we modeled the dependency among system components, the natural way to estimate system reliability is Gibbs Sampler. We suggest using Gibbs sampler combined with Antithetic method in this case. Let us consider case 3 at its extreme where each component has its own distinguished conditional dependency function, i ;i = 1;:::n. We use the Batch Means method (see chapter 10 of [50]) to estimate the mean square error of the simulation estimator and to make \raw" Gibbs sampler and Gibbs sampler combined with Antithetic method comparable in terms of variance reduction. Suppose length of each simulation run, number of Gibbs sampler iterations, is m , and let r be 50 the number of runs. Based on r simulation runs, the \raw" Gibbs sampler estimator of =E[h(X)] would be ^ ^ RG where ^ ^ RG = r i=1 ^ i r ; ^ i = m t=1 h(X t ) m (3.11) Gibbs Sampler Combined with Antithetic Method Suppose U 1 ;:::;U m is the sequence of random numbers used in computing ^ i at the i th simulation run. Instead of generating a new sequence of random numbers to estimate ^ i+1 at the i + 1 th run, we implement the notion of Antithetic method by using (1U 1 );:::; (1U m ) to compute ^ i+1 . Proposition Var( ^ ^ AG ) Var( ^ ^ RG ) Proof Supposeh(X 1 ;:::;X n ) is a monotone function of each of its coordinates, then for a set of independent random numbers, U 1 ;:::;U n , it can be proved that Cov [h(U 1 ;:::;U n );h(1U 1 ;:::; 1U n )] 0 for the proof see the corollary of theorem 8.10 on page 209 of [50]. In Raw Gibbs sampler, if one choosesm large enough, ^ i 's, can be considered asr independent random variables. Now consider Antithetic Gibbs sampler, and note that ^ i can be considered as a monotone function ofU 1 ;:::;U m . The result follows since i+1 is a monotone function of 1U 1 ;:::; 1 U m , and we can use the above mentioned corollary of [50] to conclude Cov(^ i ; ^ i+1 ) 0. 51 Remark In Gibbs Sampler since the early iterations of the Markov Chain can be strongly in uenced by the initial state chosen, it is common in practice to disregard the rst m 0 iterations (a burn-in period of m 0 iterations). 3.3. 3-State Independent Components In this section system components are independent and can be in state 1, 2, or 3. Compo- nents' probability mass functions, P i (j)P (X i =j), are known, i = 1;:::;n, j = 1; 2; 3. We are interested in using simulation to estimate system reliability, E[H]. To introduce our 3-stage simulation procedure, we rst illustrate how to modify improved stratied sampling approach when 3-state components are i.i.d. 3.3.1 3-State i.i.d Components and the Estimator of E[H] Let P (X i = j) P j , i = 1;:::;n, j = 1; 2; 3, and N = (N 1 ;N 2 ;N 3 ) where N j represents the number of components that are at state j. Note that N is a multinomial random vector, N 1 +N 2 +N 3 =n, and P k P (N = k) = n k 1 ;k 2 ;k 3 P k 1 1 P k 2 2 P k 3 3 (3.12) is easily computable, k = (k 1 ;k 2 ;k 3 ). Also, recall that we set H = h(X) which is a structure function in the context of system reliability. To estimate E[H], we use the multinomial random vector N as the \stratication variable". That is, the following identity is used, =E[H] = k E[HjN = k]P k (3.13) 52 To eciently estimate based on (3.13), on each simulation run, we estimate the quan- tities E[HjN = k] for all possible values of k = (k 1 ;k 2 ;k 3 ). This is done by letting each simulation run be as follows: 1. Generate R 1 ;:::;R n , a random permutation of 1; 2;:::;n. 2. For all possible values of k, compute H k = h(x 1 ;x 2 ;:::;x n ) where x R i = 1 for i = 1; 2;::;k 1 ,x R i = 2 fori =k 1 +1;:::;k 1 +k 2 , andx R i = 3 fori =k 1 +k 2 +1;:::;n 3. Return the estimate ^ = k H k P k That ^ is an unbiased estimator of E[H] follows from noting that the exchangeable assumption implies that H k 's are unbiased estimators of E[HjN]. This is so because given that N = k each of the n k 1 ;k 2 ;k 3 possible sets is equally likely to be any set where k j component are at state j ; j = 1; 2; 3. Let H =h(X) be the raw simulation estimator of . Proposition Var(^ ) Var(H) Proof We show that ^ can be written as a conditional expectation of the raw simulation estimator. Suppose that X = (X 1 ;X 2 ;:::;X n ) is generated by rst generating the value of N = (N 1 ;N 2 ;N 3 ), then generating a random permutation R = (R 1 ;R 2 ;:::;R n ) and settingX R i = 1 fori = 1; 2;::;k 1 ,X R i = 2 fori =k 1 +1; 2;::;k 1 +k 2 , andX R i = 3 fori = k 1 +k 2 + 1; 2;::;n. The result follows because ^ =E[HjR]. 53 3.3.2 The 3-Stage Simulation Procedure Recall that components' probability mass functions, P i (j) P (X i = j), are given, i = 1; 2;:::;n, j = 1; 2; 3. Let P (j) P n i=1 P i (j) n . For xed 0 i 1 dene 0 ~ P i 1 as a function of j, where j = 1; 2; 3, such that ~ P i = P i i P 1 i (3.14) Note that (3.14) implies ~ P i (1) + ~ P i (2) + ~ P i (3) = 1. i 's are chosen to satisfy (3.14) and to be as large as possible. From (3.14) we can conclude that, i = min min 1j3 ( P i (j) P (j) ); min 1j3 ( 1P i (j) 1 P (j) ) (3.15) Note that rearranging (3.14) would give P i = i P + (1 i ) ~ P i (3.16) for each j = 1; 2; 3. Intuitively, based on equation (3.16), we would like to express the probability mass function of each component as a mixture of two probability mass functions, P and ~ P i . To employ the modied ISS introduced in the previous section, we would like to draw most of our samples from the components' common probability mass function, P . That is, our objective is computing the components' new probability mass functions ~ P i while obtaining the maximum of i for each i = 1; 2;:::;n. Our proposed 3-stage simulation procedure follows: 54 Stage 1: Preliminary computations 1. Compute P , components' common probability mass function. 2. Compute i for all i = 1;::;n using equation (3.15). 3. Compute components' \new" probability mass functions, ~ P i 's. Now a simulation run is performed as follows : Stage 2 1. To simulate state of each component, generate a random number u i . If u i i , go to stage 3 where X i is simulated based on components' common p.m.f, P , and by employing modied ISS . 2. If u i > i then simulate X i based on the component i's \new" probability mass function, ~ P i . Set X i =j if the simulated value is j. Stage 3: Employing Modied ISS approach Consider the system conditional on the values of X i simulated at stage 2, and suppose that all the other components have reliability P . Now, use the modied improved stratied sampling introduced in the previous section to estimate the reliability of the resulting system. Remark Recall the dependency among binary components considered in section 3.2, similarly we can dene the following type of dependency among 3-state components where 55 state of each component depends on the state of the rest of the system only through the number of components that are at state 1,2, or 3, P (X i =jjX[i] = x[i])P (X i =jjN[i] = k [i] ) = i j (k [i] ) (3.17) where X [i] = (X 1 ;:::;X i1 ;X i+1 ;:::;X n ), N = (N 1 ;N 2 ;N 3 ), and N [i] j denote number of components, excluding component i, that are at state j, j = 1; 2; 3. (Note that in the preceding denition we have assumed that x[i] is so that N[i] = k [i] ). Simulation estimation of the network reliability for the 3-state components in the presence of the preceding component dependencies can be dealt with similar to our approach presented in section 3.2. When i = for all i = 1; 2;:::;n, ISS is employed after analytically computing the p.m.f of the random vector N using the Markov chain representation of N by evolving X based on iterations of a Gibbs Sampler. In the non-exchangeable case, we suggest using Gibbs sampler combined with Antithetic method. 56 Chapter 4 Dynamic Scheduling of an N-System with Reneging 4.1. Introduction Parallel server systems arise in a variety of applications contexts, including computer systems, manufacturing, and service systems settings. One specic model that is well- studied is the N-system; see, for example, Harrison [26], Bell and Williams [5], Osogami et al [45], Tezcan [56], and Down and Lewis [13]. This is because the N-system is one of the simplest parallel server system models that retains much of the complexity inherent in more general models. The basic N-system described in [26] and [5] consists of two job classes and two servers. Server 1 is dedicated, and serves only one job class. Server 2 can serve either job class. Each time server 2 nishes processing a job, and there are jobs in both buers, there is a real-time control decision concerning which class server 2 should serve next. In the case that there are only jobs in one buer, server 2 may either serve that class or idle. In the case that the system is empty, server 2 idles, and there is no control decision. 57 The basic N-system can be made more realistic by incorporating reneging; i.e., by accounting for the fact that some jobs may abandon the system without receiving ser- vice. Our objective is to show that the inclusion of reneging fundamentally changes the asymptotically optimal control policy for the N-system given in [5]. This then implies that the inclusion of reneging in a more general parallel server system, such as the model considered in the later Bell and Williams paper [6], fundamentally changes their asymp- totically optimal control policy in that more general setting as well. This conclusion is not surprising given that the results in Ward and Glynn [57] [58] imply that the inclu- sion of reneging in a single server model can dramatically impact performance measure approximation, even when only a small percentage of customers abandon. Specically, we study the system shown in Figure 4.1. Jobs arrive to buerk,k = 1; 2; in accordance with a renewal process having rate k . We call a job residing in buer k a class k job. Within a buer, jobs are served in the order of their arrival times. Server 1 can serve class 1 jobs, whereas server 2 can serve both job classes. The mean time for server 1 to serve a class 1 job is 1= 1 , the mean time for server 2 to serve a class 1 job is 1= 2 , and the mean time for server 2 to serve a class 2 job is 1= 3 . We incorporate reneging into the model by assuming that each class k job independently abandons his buer if his service has not begun within an exponentially distributed amount of time with mean 1= k (meaning that a job in service will not renege). The cost of holding a class k job ish k > 0 per unit time, and the cost of a classk job reneging isc k . The objective is to minimize the expected innite horizon discounted cost by dynamically deciding which job class server 2 should serve. 58 Figure 4.1: The N-system with reneging. This problem is too dicult to solve exactly, and so we solve it asymptotically. We focus on the heavy-trac parameter regime in [5]. Specically, 1 1 1 2 2 3 ; so that the long-run fraction of time server 2 has left over after helping server 1 serve class 1 jobs is approximately the amount that the class 2 jobs require. Furthermore, 1 > 1 so that server 2 must help server 1 in order that the number of class 1 jobs in the buer, and the number of jobs reneging, do not become very large. This is known as a complete resource pooling condition (cf. [26], Harrison and Lopez [28], Laws [36], and Kelly and Laws [35]), because in the heavy-trac limit, under an asymptotically optimal control policy, the two servers combine to form a single \super-server". 59 The main results of our study are as follows. First, assume that (h 2 +c 2 2 ) 3 < (h 1 +c 1 1 ) 2 ; (4.1) so that class 2 is the cheapest class under the c-rule. We then focus on the case that 1 is large enough so that 1 + 2 + > 2 (h 1 +c 1 1 ) 3 (h 2 +c 2 2 ) ; (4.2) where is the discount rate. Then, an asymptotically optimal policy for an expected innite horizon discounted cost setting is a dynamic priority threshold control policy that is based on two thresholds. There is one threshold L on queue 1 that serves the same purpose as in [5], to prevent excessive idling of server 1 caused by server 2 helping \too much". The other threshold M is on what is known as the workload process in the framework of [28]. More specically, for a specied threshold level M (see Denition 5.2 in Section 5.2), if Q 1 (t) andQ 2 (t) are the number of class 1 and class 2 jobs respectively at time t, and Q 1 (t)>L, then server 2 should give priority to class 1 over class 2 when Q 1 (t) + 2 3 Q 2 (t)M; (4.3) and server 2 should give priority to class 2 over class 1 otherwise. If Q 1 (t) < L, then server 2 does not serve class 1, even if there are no class 2 jobs present. Then, as shown in Chapter 6, when (4.3) holds the number of class 1 jobs is very small, and, when (4.3) does not hold, the number of class 2 jobs is very small. When (4.1) holds but not (4.2), 60 an asymptotically optimal policy is the policy proposed in Denition 5.1 in [5]. In the case that (4.1) does not hold, the condition 2 + 1 + > 3 (h 2 +c 2 2 ) 2 (h 1 +c 1 1 ) (4.4) determines whether or not a workload-based dynamic priority threshold policy is optimal. The implications of our results are that in heavy trac the system manager is able to rapidly obtain a completely dierent system conguration by switching the priority of server 2. In particular, the system manager may move from a state in which all jobs are held in one buer to a state in which all jobs are held in the other buer in a very short amount of time. This is because demand is large and service is fast, so that switching the priority of server 2 enables the jobs in one buer to be served very quickly, while the new arrivals backlog in the other buer. This is important because this insight should hold for more general parallel server system congurations. The workload-based dynamic priority threshold policy we propose is very similar to the policy proposed in Ata and Rubino [2] for a general parallel server system (that need not be anN-system) that models a make-to-order production system with job cancelation. For the N-system, the only dierence between the Brownian control problem they solve and the one we solve is that they use an average cost criterion, and we use a discounted cost criterion. However, they prove neither weak convergence of their proposed policy nor asymptotic optimality. Instead, they verify its eectiveness via simulation. This is because the proof is very challenging. This study provides the rigorous proof that 61 the solution to the Brownian control problem for the N-system does indeed motivate an asymptotically optimal policy for the original system, under the discounted cost criterion. We rst provide the notation and terminology to be used throughout our study. In Section 4.2, we complete the description of theN-system considered here. This includes a description of the primitive stochastic processes in our model, allowed scheduling control policies, and a specication of dynamic equations satised by the queue length processes. In Section 4.3, we describe the asymptotic regime in which we seek to analyze the per- formance of control policies for our system. In particular, we specify assumptions on the stochastic primitives that imply our system is asymptotically in heavy trac (cf. Harri- son [27]), and satises the complete resource pooling condition of [28]. Also, we describe the normalization of the queue length process via diusive scaling and we specify the as- sociated cost function. Next, following the general scheme proposed by Harrison [25], in Chapter 5, we state the formal Brownian control problem associated with our N-system, and, in Section 5.1, we solve the Brownian control problem. In contrast to [5], this is a non-trivial control problem to solve, and is an important part of this work. In Section 5.2, we translate the solution to the Brownian control problem into a policy for the original system. We then focus on the parameter regime in which conditions (4.1) and (4.2) hold, because this is the part of the parameter regime in which an asymptotically optimal pol- icy diers from [5]. In Chapter 6, we prove the weak convergence of the diusion-scaled system processes under our proposed policy to the Brownian control problem solution. This establishes that the diusion-scaled number of class 1 jobs weakly converges to zero when (4.3) holds, and that the diusion-scaled number of class 2 jobs weakly converges to zero when (4.3) does not hold, so that there is a state-space collapse. Finally, in 62 Section 6.2, we show that no other policy can achieve a lower limiting cost in the heavy trac limit, which establishes that our proposed policy is asymptotically optimal. Chapter 7 provides the following supporting material: the statement and proof of an oscillation inequality for the solution of a one-dimensional, linearly generalized, perturbed Skorokhod problem, a technical argument omitted from the proof of the proposition that shows weak convergence (Proposition 6.2), and the proofs of all the Lemmas stated in this study. 4.1.1 Notation and Terminology The m-dimensional (m 1) Euclidean space will be denoted by< m and< + will denote [0;1). Letjj denote the norm on< m given byjxj := P m i=1 jx i j for x2< m . Vectors in < m should be treated as column vectors unless indicated otherwise, inequalities between vectors should be interpreted componentwise, the transpose of a vector a will be denoted by a 0 , the diagonal matrix with the entries of a vector a on its diagonal will be denoted by diag(a), and the dot product of two vectors a and b will be denoted by ab. The identity function is denoted by e, so that e(t) =t for all t 0. For x2<, let [x] denote the integer part of x. For each positive integer m, let D m be the set of all functions ! :< + !< m that are right continuous on< + and have nite left limits on (0;1). The member of D m that stays at the origin in< m for all time will be denoted by 0. For !2 D m and t 0, let k!k t := sup s2[0;t] j!(s)j; 63 and for 0t 1 t 2 <1, let Osc(w; [t 1 ;t 2 ]) := sup t 1 stt 2 jw(t)w(s)j: Consider D m to be endowed with the usual Skorokhod J 1 -topology (cf. Billingsley [7] or Ethier and Kurtz [19]). LetM m denote the Borel -algebra on D m associated with the J 1 -topology. This is the same -algebra as the one generated by the coordinate maps; i.e.,M m =f!(s) : 0s<1g. All of the continuous-time processes in this study will be assumed to have sample paths in D m for some m 1. SupposefW n g 1 n=1 is a sequence of processes with sample paths in D m for some m 1. Then we say thatfW n g 1 n=1 is tight if and only if the probability measures induced by the W n 's on (D m ;M m ) form a tight sequence; i.e., they form a weakly relatively compact sequence in the space of probability measures on (D m ;M m ). The notation \W n ) W ", where W is a process with sample paths in D m , will mean that the probability measures induced by the W n 's on (D m ;M m ) converge weakly to the probability measure on (D m ;M m ) induced by W . If for each n, W n and W are dened on the same probability space, we say that W n converges to W uniformly on compact time intervals in probability if P(kW n Wk t )! 0 as n!1 for each > 0 and all t 0. The abbreviations \u.o.c.", \i.p.", and \a.s." refer respectively to \uniformly on compact intervals", \in probability", and \almost surely". 64 4.2. The N System We follow the notation and model assumptions in [5], except that we add additional nota- tion necessary for the model to include reneging. Subsection 4.2.1 sets up the stochastic primitives in our model formulation, and Subsection 4.2.2 details the scheduling control and system evolution equations. We delay the precise denition of the cost function until Section 4.3, since it is formulated in terms of normalized queue lengths where the nor- malization is on diusion scale, commensurate with the heavy trac limiting regime in which we consider our model. 4.2.1 Stochastic Primitives All random variables and stochastic processes in our model are assumed to be dened on a complete probability space ( ;F; P). The expectation operator under P will be denoted by E and P(A;B) will mean P(A\B). We assume that the system is initially empty. For k = 1; 2; we take as given a sequence of strictly positive i.i.d. random variables fu k (i);i = 1; 2;:::g with mean 1 k 2 (0;1) and squared coecient of variation 2 k 2 [0;1). Fori = 1; 2;:::, we interpretu k (i) as the time between the arrival of the (i 1)st and the ith arrival to class k, where the \0th arrival" occurs at time zero. Setting k (0) = 0 and k (n) := n X i=1 u k (i); for n = 1; 2;:::; we dene A k (t) := supfn 0 : k (n)tg for all t 0: 65 Then A k is a renewal process, A k (t) counts the number of arrivals to class k that have occurred in [0;t], and k is the long run arrival rate to class k. Forj = 1; 2; 3, we take as given a sequence of strictly positive, i.i.d. random variables fv j (i);i = 1; 2;:::g with mean 1 j 2 (0;1) and squared coecient of variation 2 j 2 [0;1). For each j, we interpret v j (i);i 1; as the amount of service time required by the ith job to be processed by activity j, where activity 1 = processing of class 1 jobs by server 1 activity 2 = processing of class 1 jobs by server 2 activity 3 = processing of class 2 jobs by server 2 Note that j is the long run rate at which activity j can process its associated class if the associated server works continuously and exclusively on this class. For j = 1; 2; 3, let j (0) := 0, j (n) := n X i=1 v j (i); for n = 1; 2;:::; and S j (t) := supfn 0 : j (n)tg for all t 0: Then, S 1 , S 2 , S 3 are renewal processes and S 1 (t) represents the number of class 1 jobs that server 1 could complete if that server worked continuously in [0;t], and for j = 2; 3, S j (t) is the number of class j 1 jobs that server 2 could complete if that server worked continuously and exclusively on class j 1 jobs in [0;t]. Also, letN 1 andN 2 be two standard, independent Poisson processes. We will use N 1 and N 2 to track the cumulative number of jobs that have reneged from each buer. 66 We assume that the interarrival time sequencesfu k (i);i = 1; 2;:::g, k = 1; 2, service time sequencesfv j (i);i = 1; 2;:::g,j = 1; 2; 3, and Poisson processesN k ,k = 1; 2, are all mutually independent. 4.2.2 Scheduling Control Scheduling control of the system is exerted through a three-dimensional service time allocation process T (t) := (T 1 (t);T 2 (t);T 3 (t)) 0 ;t 0: Forj = 1; 2; 3,T j (t) is the cumulative amount of service time devoted to activity j in the time interval [0;t]. Then I 1 (t) :=tT 1 (t) (4.5) is the cumulative idletime of server 1 up to time t, I 2 (t) :=tT 2 (t)T 3 (t) (4.6) is the cumulative idletime of server 2 up to time t,S 1 (T 1 (t)) is the number of jobs served by server 1 up to timet, andS j (T j (t)) is the number of classj1 jobs served by server 2, forj = 2; 3, up to timet. Let _ T k (t) := d dt T k (t) fort> 0 represent whether or not activity j is being undertaken at time t. The number of class k, k = 1; 2, jobs at time t is Q 1 (t) := A 1 (t)S 1 (T 1 (t))S 2 (T 2 (t))R 1 (t) (4.7) Q 2 (t) := A 2 (t)S 3 (T 3 (t))R 2 (t); (4.8) 67 where R 1 (t) := N 1 Z t 0 1 Q 1 (s) 1f _ T 1 (s)> 0g 1f _ T 2 (s)> 0g ds (4.9) R 2 (t) := N 2 Z t 0 2 Q 2 (s) 1f _ T 3 (s)> 0g ds (4.10) are the cumulative number of class 1 and 2 jobs that have reneged by time t, and 1 and 2 are the respective class 1 and class 2 reneging rates. Note that a job in service does not renege. The scheduling control T must satisfy the following properties. For j = 1; 2; 3 and k = 1; 2, and I(t) := (I 1 (t);I 2 (t)) 0 and Q(t) := (Q 1 (t);Q 2 (t)) 0 for all t 0 that satisfy (4.5)-(4.8), T j (t)2F for each t 0; (4.11) T j () is continuous and non-decreasing with T j (0) = 0; (4.12) I k () is continuous and non-decreasing with I k (0) = 0; (4.13) Q k (t) 0 for all t 0: (4.14) Note that conditions (4.12)-(4.13) imply that for j = 1; 2; 3, T j is uniformly Lipschitz continuous with a Lipschitz constant less than or equal to one. 68 4.3. Heavy Trac Assumptions, Scaling, and the Cost Function The problem of nding a control policy that minimizes costs associated with holding jobs and reneging is not solvable via an exact analysis. (However, it may be possible to obtain bounds, as Chan and Farias [8] do in a discrete-time parallel server system model.) Therefore, we perform an asymptotic analysis. In particular, we consider the asymptotic regime associated with heavy trac limit theorems in which the queue length process is normalized with diusive scaling. This corresponds to viewing the system over long intervals of time of order r 2 (where r will tend to innity in the asymptotic limit) and regarding a single job as only having a small contribution to the overall cost, where this is quantied to be of order 1=r. Formally, we consider a sequence of N-systems indexed by r, where r tends to innity through a sequence of values in [1;1). These systems all have the same basic structure as that described in the last section; however, the arrival, service, and reneging rates, scheduling control and cost function (dened below) may vary with r. We shall indicate the dependence of relevant parameters and processes on r by appending a superscript to them. We assume that the interarrival and service times are given for each r, k = 1; 2, j = 1; 2; 3, i = 1; 2;:::, by u r k (i) := 1 r k u k (i);v r j (i) := 1 r j v j (i); 69 where the u k (i); v j (i) do not depend on r, have mean one and squared coecients of variation 2 k , 2 j respectively. The Poisson processes N 1 and N 2 do not vary with r, but the reneging rates r 1 and r 2 do. The following assumption is Assumption 3.1 in [5]. Assumption 4.1. There are strictly positive nite constants 1 ; 2 ; 1 ; 2 ; 3 such that (i) 1 > 1 , (ii) 1 1 1 2 = 2 3 , and, as r!1, (iii) r k ! k ;k = 1; 2, (iv) r j ! j ;j = 1; 2; 3, (v) r r 2 r 1 r 1 r 2 1 1 2 ! 0, (vi) r r 3 r 2 r 3 2 3 ! 0. We also require an assumption on the reneging rates. Assumption 4.2. There are strictly positive nite constants 1 ; 2 such that as r!1 r 2 r 1 ! 1 and r 2 r 2 ! 2 : (This structure is a convenient means of allowing the sequence of systems to approach heavy trac by simply changing arrival, service, and reneging rates while keeping the underlying sources of variability u k (i); v j (i);N 1 , and N 2 xed. For a rst reading, the reader may like to simply choose r k = k , r j = j , and r k = k =r 2 for k = 1; 2 and j = 1; 2; 3, and for all r.) 70 Note that we assume the convergences in (v) and (vi) are to zero for convenience in solving the diusion control problem in Section 5.1; however, we expect the same results structure to hold when the convergences in (v) and (vi) are to any real numbers. For each xedr and control policyT r with associated queue lengths Q r and idletime processes I r in the rth system, dene for each t 0 T r (t) := r 2 T r (r 2 t); Q r (t) := r 2 Q r (r 2 t); ^ A r k (t) := r 1 A r k (r 2 t) r k r 2 t ; k = 1; 2; ^ S r j (t) := r 1 S r (r 2 t) r j r 2 t ; j = 1; 2; 3; ^ Q r k (t) := r 1 Q r k (r 2 t); k = 1; 2 ^ I r k (t) := r 1 I r k (r 2 t); k = 1; 2 ^ N r k (t) := r 1 N k (r 2 t)r 2 t : Then, equations (4.7)-(4.8) yield the following expressions for the normalized (diusion scaled) queue length processes: ^ Q r 1 (t) := ^ A r 1 (t) ^ S r 1 (T r 1 (t)) ^ S r 2 (T r 2 (t)) (4.15) ^ N r 1 Z t 0 r 2 r 1 Q r 1 (s) 1 r 2 1f _ T r 1 (r 2 s)> 0g 1 r 2 1f _ T r 2 (r 2 s)> 0g ds +r r 1 t r 1 T r 1 (t) r 2 T r 2 (t) Z t 0 r 2 r 1 ^ Q r 1 (s) 1 r 1f _ T r 1 (r 2 s)> 0g 1 r 1f _ T r 2 (r 2 s)> 0g ds ^ Q r 2 (t) := ^ A r 2 (t) ^ S r 3 T r 3 (t) ^ N r 2 Z t 0 r 2 r 2 Q r 2 (s) 1 r 2 1f _ T r 3 (r 2 s)> 0g ds (4.16) +r r 2 t r 3 T r 3 (t) Z t 0 r 2 r 2 ^ Q r 2 (s) 1 r 1f _ T r 3 (r 2 s)> 0g ds: 71 On combining Assumption 4.1 with the nite variance and mutual independence of the stochastic primitivesf u k (i);i = 1; 2;:::g;f v j (i);i = 1; 2;:::g;N k for k = 1; 2 and j = 1; 2; 3, we may deduce from renewal process functional central limit theorems (see, for example Iglehart and Whitt [33]) that ^ A r ; ^ S r ; ^ N r ) ~ A; ~ S; ~ N as r!1; where ~ A; ~ S; ~ N are mutually independent, ~ A is a two-dimensional driftless Brownian mo- tion that starts from the origin and has diagonal covariance matrix diag( 1 2 1 ; 2 2 2 ), ~ S is a three-dimensional driftless Brownian motion that starts from the origin and has diagonal covariance matrix diag( 1 2 1 ; 2 2 2 ; 3 2 3 ), and ~ N is a two-dimensional driftless Brownian motion that starts from the origin and has diagonal covariance matrix that is the identity matrix. For the rth system, we consider a mean cumulative discounted cost for the use of a control T r having associated normalized queue length process ^ Q r : ^ J r :=E Z 1 0 e t v ^ Q r (t)dt ; (4.17) where > 0 is a constant discount rate and v := (h 1 +c 1 1 ;h 2 +c 2 2 ) 0 72 is the vector of holding costs that includes both the cost of physically having a job in the buer and reneging costs. Then, the control problem we would like to solve is as follows. Denition 4.1. (The control problem for the rth system) minimize ^ J r (T r ) by choosing a scheduling control (T r 1 ;T r 2 ;T r 3 ) that satises (4.1)-(4.8). Note that the k fork = 1; 2 in the above is the limit in Assumption 4.2 asr becomes large, whereas the h k andc k fork = 1; 2 are taken to represent costs associated with the diusion-scaled queue-length. Remark 4.1. For simplicity in notation in this remark, assume that a job in service can renege, so that R k (t) =N k Z t 0 k Q k (s)ds : The heavy-trac limit will be the same under either this denition or the one given in (4.9)-(4.10). Instead of directly dening the mean cumulative discounted cost in (4.17) in terms of the normalized queue-length processes, we could have rst dened the mean cumulative discounted cost for the original system J(T ) :=E Z 1 0 e t (h1Q1(t) +h2Q2(t))dt + c1dN1 Z t 0 1Q1(s)ds +c2dN2 Z t 0 2Q2(s)ds ; 73 and then allowed the parameters h 1 ;h 2 ; ;c 1 , and c 2 to also depend on r so that J r (T ) = E 2 6 6 4 Z 1 0 e r t 2 6 6 4 0 B B @ h r 1 Q r 1 (t) +h r 2 Q r 2 (t)dt +c r 1 dN 1 R t 0 r 1 Q r 1 (s)ds +c r 2 dN 2 R t 0 r 2 Q r 2 (s)ds 1 C C A 3 7 7 5 3 7 7 5 = E 2 6 6 4 Z 1 0 e r r 2 u 2 6 6 4 r 3 h r 1 ^ Q r 1 (u) +h r 2 ^ Q r 2 (u) du + c r 1 dN 1 R r 2 u 0 r 1 Q r 1 (s)ds +c r 2 dN 2 R r 2 u 0 r 2 Q r 2 (s)ds 3 7 7 5 3 7 7 5 : The strong law of large numbers for Poisson processes implies that N(r 2 t)!r 2 t a.s. as r!1. Therefore, we have the following, dN k Z r 2 u 0 r k Q r k (s)ds ! r k r 2 Q r k (r 2 u)du; and so J r (T ) =E Z 1 0 e r r 2 u h r 3 h r 1 ^ Q r 1 (u) +h r 2 ^ Q r 2 (u) du +r c r 1 r 2 r 1 ^ Q r 1 (u) +c r 2 r 2 r 2 ^ Q r 2 (u) i : Sincer 2 r k ! k , k = 1; 2, by Assumption 4.2, this implies that we would need to assume r 3 h r k !h k and rc r k !c k ; k = 1; 2; as r!1: In other words, holding costs in the original system should be very small (of order r 3 ) compared to the cost of reneging (of order r 1 ). Intuitively, this makes sense given that the parameter regime is one in which the reneging rate is becoming very small and the number of jobs in each buer is becoming large (because the system is in heavy trac). 74 In addition to Assumptions 4.1 and 4.2, we make the following exponential moment assumptions which ensure that certain large deviations estimates hold for the renewal processes A r k ;S r j associated with the interarrival and service times. This is Assumption 3.3 in [5]. Assumption 4.3. For k = 1; 2;j = 1; 2; 3, and all i 1, let (cf. [33]): u k (i) = 1 k u k (i); v j (i) = 1 j v j (i): Assume that there is a non-empty open neighborhood,O, of 02< such that for all l2O and all i 1, a k (l) := log E e lu k (i) <1; for k = 1; 2; and s j (l) := log E e lv j (i) <1; for j = 1; 2; 3: Note that we have not made the Assumption 3.2 in [5], that h 1 2 h 2 3 . This is because the appropriate division of the parameter space, that was outlined in the sixth paragraph of the Introduction, will become apparent in Section 5.1, where we solve the (formal) Brownian control problem associated with our sequence of N-systems. 75 Chapter 5 Approximating Brownian Control Problem and the Proposed Policy We formally show how to derive the approximating Brownian control problem that arises in the heavy trac limit when passing to the limit as r!1 in the control problem in Denition 4.1. This requires the assumption that in uid or law of large numbers scale, the scheduling control achieves the long run rates required in Assumption 4.1 for a balanced system in the heavy trac limit; i.e., that T r )T ? ; as r!1; for T ? (t) := t; 1 1 2 t; 2 3 t : The preceding assumption implies that the uid scaled queue length process converges to 0 in the limit. That is, Q r 1 ; Q r 1 ) (0; 0); as r!1: 76 The rst step is to note that the normalized queue length processes in (4.15) and (4.16) can be written in terms of the diusion-scaled processes that we expect to converge to a two-dimensional Brownian motion ^ X r 1 (t) := ^ A r 1 (t) ^ S r 1 (T r 1 (t)) ^ S r 2 (T r 2 (t)) ^ N r 1 ( r 1 (t)) (5.1) +r r 2 r 1 r 1 r 2 1 1 2 t ^ X r 2 (t) := ^ A r 2 (t) ^ S r 3 (T r 3 (t)) ^ N r 2 ( r 2 (t)) +r r 3 r 2 r 3 2 3 t; (5.2) for r 1 (t) := Z t 0 r 2 r 1 Q r 1 (s) 1 r 2 1f _ T r 1 (r 2 s)> 0g 1 r 2 1f _ T r 2 (r 2 s)> 0g ds r 2 (t) := Z t 0 r 2 r 2 Q r 2 (s) 1 r 2 1f _ T r 3 (r 2 s)> 0g ds; and the normalized deviation process ^ Y r (t) :=r T ? T r ; (5.3) as follows ^ Q r 1 (t) = ^ X r 1 (t) Z t 0 r 2 r 1 ^ Q r 1 (s) 1 r 1f _ T r 1 (r 2 s)> 0g 1 r 1f _ T r 2 (r 2 s)> 0g ds (5.4) + r 1 ^ Y r 1 (t) + r 2 ^ Y r 2 (t) ^ Q r 2 (t) = ^ X r 2 (t) Z t 0 r 2 r 2 ^ Q r 2 (s) 1 r 1f _ T r 3 (r 2 s)> 0g ds + r 3 ^ Y r 3 (t): (5.5) 77 The second step is to observe that there is a connection between the normalized deviation process ^ Y r and the diusion-scaled idletime processes ^ I r 1 and ^ I r 2 . In particular, from (4.5), ^ Y r 1 (t) = 1 r (r 2 tT r (r 2 t)) = ^ I r 1 (t); and from (4.6), ^ Y r 2 (t) + ^ Y r 3 (t) = 1 r (r 2 tT r 2 (r 2 t)T r 3 (r 2 t)) = ^ I r 2 (t): Finally letting r!1 and using Assumptions 4.1 and 4.2, we arrive at the following Brownian control problem. Denition 5.1. (Brownian Control Problem) minimize E Z 1 0 e t v ~ Q(t)dt using a three-dimensional control process ~ Y = ( ~ Y 1 ; ~ Y 2 ; ~ Y 3 ) 0 such that ~ Q 1 (t) := ~ X 1 (t) Z t 0 1 ~ Q 1 (s)ds + 1 ~ Y 1 (t) + 2 ~ Y 2 (t) 0 for all t 0; ~ Q 2 (t) := ~ X 2 (t) Z t 0 2 ~ Q 2 (s)ds + 3 ~ Y 3 (t) 0 for all t 0; ~ I 1 := ~ Y 1 is non-decreasing, ~ I 1 (0) = 0; ~ I 2 := ~ Y 2 + ~ Y 3 is non-decreasing, ~ I 2 (0) = 0; ~ X is a two-dimensional Brownian motion with drift zero that starts from the origin and has diagonal covariance matrix diag( 1 2 1 + 1 2 1 + 2 2 ( 1 1 ); 2 ( 2 2 + 2 3 )). 78 The key to solving the Brownian control problem in Denition 5.1 is to write the \equivalent workload formulation", which follows the general procedure outlined in [29], and is very close to the specic procedure in [2]. For ~ Q and ~ I as in Denition 5.1, let ~ W (t) := y ~ Q(t) = ~ (t) Z t 0 1 ~ Q 1 (s) + 2 3 2 ~ Q 2 (s) ds + ~ L(t); where y := (1; 2 = 3 ), ~ :=y ~ X, and ~ L(t) = 1 ~ I 1 (t) + 2 ~ I 2 (t) for all t 0: (5.6) The process ~ has zero drift and variance 2 := 1 2 1 + 1 2 1 + 2 2 ( 1 1 ) + 2 3 2 2 ( 2 2 + 2 3 ) : Also let A(w) :=fq2< 2 + :q 1 + 2 3 q 2 =wg; and dene the functions f and g having domain< 2 + and range< + as follows f(q) := 1 q 1 + 2 3 2 q 2 g(q) := vq = (h 1 +c 1 1 )q 1 + (h 2 +c 2 2 )q 2 : 79 The equivalent workload formulation assumes that the value of the process ( ~ Q 1 ; ~ Q 2 ) at time t> 0 can be chosen as function of ~ W (t), and is as follows, minimize q2A(w) E R 1 0 e t g q( ~ W (t)) dt subject to: ~ W (t) = ~ (t) R t 0 f(q( ~ W (s)))ds + ~ L(t) 0 ~ L(0) = 0 and is non-descreasing (5.7) We end this section with a Proposition that establishes the equivalence between the equivalent workload formulation (5.7) and the Brownian control problem in Denition 5.1. The following proposition has a proof very similar to that of Propositions 2 and 3 in [2], and so is omitted. Proposition 5.1. Every admissible policy ~ Y for the Brownian control problem in De- nition 5.1 yields an admissible policy ( ~ L;q) for the workload control problem formulation in (5.7), and these two policies have the same cost. Similarly, for any admissible policy ( ~ L;q) for the workload control problem formulation in (5.7), there exists an admissible policy ~ Y for the Brownian control problem in Denition 5.1, and its cost is less than or equal to that of the policy ( ~ L;q) for the workload control problem formulation in (5.7). Given a pair (q; ~ L) that satises the constraint in (5.7), ~ W (q; ~ L) (t) = ~ (t) Z t 0 f(q( ~ W ( ~ ; ~ L) (s)))ds + ~ L(t) 0 is the associated workload at time t 0. 80 5.1. The Brownian Control Problem Solution It follows from Proposition 5.1 that to solve the Brownian control problem, we can rst solve the equivalent workload formulation (5.7). In Section 5.1.1, we provide the Bellman equation and verication Lemma for (5.7). We then separate the parameter space into the part in which the solution to the Bellman equation denes a dynamic threshold control policy, and the part in which the solution to the Bellman equation denes a static control policy in which either q 1 is always 0 or q 2 is always 0. Next, in Section 5.1.2, we solve the Bellman equation in the part of the parameter space in which a threshold control policy is optimal. This requires showing that the solutions to two ordinary dierential equations (ode's) can be pasted together in such a way that the resulting function is twice continuously dierentiable. Finally, in Section 5.1.3, we verify that the Bellman equation solution implies that a static control policy is optimal in the complementary part of the parameter space. 5.1.1 The Bellman Equation The Bellman equation is minimize q2A(w) 2 2 V 00 (w)f(q)V 0 (w) V (w) +g(q) = 0 for all w 0 subject to: V is twice-continuously dierentiable. V 0 is bounded, increasing, and has V 0 (0) = 0: (5.8) Letq ? denote a solution to (5.8). The following verication Lemma (whose proof utilizes the generalized Ito formula, for which Section 4.7 in Harrison [24] is a convenient reference) 81 provides a solution to the equivalent workload formulation (5.7) in terms of q ? and the process ~ L ? that satises Z t 0 ~ W ? (s)d ~ L ? (s) = 0; ~ L ? (0) = 0; ~ L ? is continuous and non-decreasing (5.9) for ~ W ? (t) := ~ W (q ? ; ~ L ? ) (t); t 0: Lemma 5.1. The policy (q ? ; ~ L ? ) solves the equivalent workload formulation (5.7). In particular, for any other pair (~ q; ~ L) that satisfy the constraints in (5.7), and ~ W (t) := ~ W (~ q; ~ L) (t);t 0, E Z 1 0 e t g(q( ~ W (t)))dtE Z 1 0 e t g(q ? ( ~ W ? (t)))dt: Furthermore, for V that satises (5.8) V (0) =E Z 1 0 e t g(q ? ( ~ W ? (t)))dt: The solution to the Bellman equation (5.8) depends on the values of the parameters h 1 ,c 1 , 1 ,h 2 ,c 2 , 2 , 2 , 3 , and . We next show the conditions under which the solution to (5.8) will be a threshold control policy, and the conditions under which the solution to (5.8) will be a static control policy. We prove the conditions we provide in this subsection are correct in Sections 5.1.2 and 5.1.3. 82 First, we re-write (5.8) in a form that is more insightful for analysis. Substituting for f and g, and performing some simple algebraic manipulations, we nd that min q2A(w) 2 2 V 00 (w)f(q)V 0 (w) V (w) +g(q) = 2 2 V 00 (w) + (h 1 +c 1 1 1 V 0 (w))w V (w) +U w (V 0 (w)) for U w (x) := min q 2 2 h 0; 3 2 w i (h 2 +c 2 2 ) 2 3 (h 1 +c 1 1 ) + 2 3 ( 1 2 )x q 2 Hence (5.8) is equivalent to 2 2 V 00 (w) + (h 1 +c 1 1 1 V 0 (w))w V (w) +U w (V 0 (w)) = 0: (5.10) It then follows that q ? (w) = 8 > > < > > : (w; 0) if (h 2 +c 2 2 ) 2 3 (h 1 +c 1 1 ) + 2 3 ( 1 2 )V 0 (w)> 0 0; 3 2 w if (h 2 +c 2 2 ) 2 3 (h 1 +c 1 1 ) + 2 3 ( 1 2 )V 0 (w) 0 : We explicitly solve (5.10) in the case that (4.1) holds, which assumes that (h 1 +c 1 1 ) 2 > (h 2 +c 2 2 ) 3 : This is exactly the case in [5] when 1 = 2 = 0. We show in Section 5.1.2 that when also (4.2) holds, which assumes that 83 1 + 2 + > 2 (h 1 +c 1 1 ) 3 (h 2 +c 2 2 ) ; there exists a unique b ? such that q ? (w) = 8 > > < > > : (0; 3 2 w) if wb ? (w; 0) if w>b ? (5.11) Otherwise, when (4.2) does not hold, we show in Section 5.1.3 that q ? (w) = (0; 3 2 w) for all w 0: (5.12) In summary, a threshold control policy solves the Brownian control problem in Deni- tion 5.1 when (4.2) holds, and a static policy consistent with [5] is optimal otherwise. The intuition is that in general jobs should be held in the cheaper buer 2. However, if the reneging rate 1 is large compared to the other system parameters, it may lower the long-run average cost to take advantage of the higher reneging rate from buer 1 when there are many jobs in the system. This is a non-greedy, forward-looking policy. For the complementary parameter regime; i.e., when condition (4.1) does not hold, so that (h 2 +c 2 2 ) 3 > (h 1 +c 1 1 ) 2 ; the solution is symmetric. Specically, when also 2 + 1 + > (h 2 +c 2 2 ) 3 (h 1 +c 1 1 ) 2 ; 84 there exists a unique b ? 2 such that q ? (w) = 8 > > < > > : (w; 0) if wb ? 2 (0; 3 2 w) if w>b ? 2 ; and otherwise q ? (w) = (w; 0) for all w 0: We omit the proofs for this case because the arguments are essentially identical to the case that condition (4.1) holds. 5.1.2 The Case that a Threshold Control Policy Solves the Brownian Control Problem We prove the following Theorem. Theorem 5.1. Assume conditions (4.1) and (4.2) hold. Then, there exists ab ? such that q ? in (5.11) solves (5.10). This b ? is specied so that (V;b ? ) satisfy V (w) = 8 > > < > > : V 2 (w) when w2 [0;b ? ] V 1 (w) when w>b ? for V 2 and V 1 that solve 85 2 2 V 00 2 (w) 2 wV 0 2 (w) V 2 (w) + 3 2 (h 2 +c 2 2 )w = 0 for all w2 [0;b ? ] 2 2 V 00 1 (w) 1 wV 0 1 (w) V 1 (w) + (h 1 +c 1 1 )w = 0 for all wb ? subject to: V 0 2 is increasing on [0;b ? ] V 0 1 is increasing on (b ? ;1) and bounded V 0 2 (0) = 0 V 0 1 (b ? ) =V 0 2 (b ? ) = (h 1 +c 1 1 ) 3 2 (h 2 +c 2 2 ) ( 1 2 ) V 1 (b ? ) =V 2 (b ? ) and V 00 1 (b ? ) =V 00 2 (b ? ): (5.13) Proof: It is sucient to solve (5.13). To do this, we rst consider an arbitrary b > 0 and solve 2 2 V 00 2 (w) 2 wV 0 2 (w) V 2 (w) + 3 2 (h 2 +c 2 2 )w = 0 for all w2 [0;b] subject to: V 0 2 (0) = 0 and V 0 2 (b) = (h 1 +c 1 1 ) 3 2 (h 2 +c 2 2 ) ( 1 2 ) (5.14) and 2 2 V 00 1 (w) 1 wV 0 1 (w) V 1 (w) + (h 1 +c 1 1 )w = 0 for all wb subject to: V 0 1 (b) = (h 1 +c 1 1 ) 3 2 (h 2 +c 2 2 ) ( 1 2 ) V 0 1 is bounded and increasing on [b;1) (5.15) 86 We then show there exists b ? such that V 1 (b ? ) =V 2 (b ? ) and V 00 1 (b ? ) =V 00 2 (b ? ) and V 0 2 is increasing on [0;b ? ] (5.16) Together (5.14)-(5.16) provide a solution to (5.13). The solutions to (5.14) and (5.15) are expressed in terms of Kummer's function M(a;b;z) := 1 + az b + (a) 2 (b) 2 z 2 2! +::: + (a) n (b) n z n n! + ; where (a) 0 = 1 and (a) n =a(a + 1)(a + 2):::(a +n 1); n = 1; 2;:::: Specically, V 2 (w) := d 1 (b)M 2 2 ; 1 2 ; 2 2 w 2 +d 2 M 2 2 ; 1 2 ; 2 2 w 2 Z w 0 exp( 2 2 y 2 ) M( 2 2 ; 1 2 ; 2 2 y 2 ) 2 dy + 3 2 (h 2 +c 2 2 ) 2 + w; 87 where d 1 (b) := (h 1 +c 1 1 ) 3 2 (h 2 +c 2 2 ) 1 2 3 2 (h 2 +c 2 2 ) 2 + 2 2 bM( 2 2 + 1; 3 2 ; 2 2 b 2 ) + 2 2 1 b 3 2 (h 2 +c 2 2 ) 2 + exp( 2 2 b 2 ) M( 2 2 ; 1 2 ; 2 2 b 2 )M( 2 2 + 1; 3 2 ; 2 2 b 2 ) + 3 2 (h 2 +c 2 2 ) 2 + Z b 0 exp( 2 2 y 2 ) M( 2 2 ; 1 2 ; 2 2 y 2 ) 2 dy d 2 := 3 2 (h 2 +c 2 2 ) 2 + solves (5.14), and V 1 (w) := (h 1 +c 1 1 ) 3 2 (h 2 +c 2 2 ) 1 2 h 1 +c 1 1 1 + M( 2 1 ; 1 2 ; 1 2 w 2 ) R 1 w exp( 1 2 y 2 ) M( 2 1 ; 1 2 ; 1 2 y 2 ) 2 dy 2 2 bM( 2 1 + 1; 3 2 ; 1 2 b 2 ) R 1 b exp( 1 2 y 2 ) M( 2 1 ; 1 2 ; 1 2 y 2 ) 2 dy exp( 1 2 b 2 ) M( 2 1 ; 1 2 ; 1 2 b 2 ) + (h 1 +c 1 1 ) 1 + w solves (5.15). To see this, as in Ward and Kumar [59], use equation #103 on page 143 of [46] to nd one solution of the homogeneous equation having the form 2 2 V 00 (w)wV 0 (w) V (w) = 0. Then, apply formula (2) on page 129 of Polyanin and Zaitsev [46] to obtain the general solution to the homogeneous equation. Next, use the method of undetermined coecients 88 to nd the general solution to the non-homogeneous equations in (5.25) and (5.26). Then, use relation 13.4.9 in Slater [55] M 0 (a;b;z) = a b M(a + 1;b + 1;z) (5.17) to nd V 0 2 (w) := d 1 (b)2 2 wM 2 2 + 1; 3 2 ; 2 2 w 2 +d 2 2 2 wM 2 2 + 1; 3 2 ; 2 2 w 2 Z w 0 exp( 2 2 y 2 ) M( 22 ; 1 2 ; 2 2 y 2 ) 2 dy + exp( 2 2 w 2 ) M( 22 ; 1 2 ; 2 2 w 2 ) ! + 3 2 (h 2 +c 2 2 ) 2 + and V 0 1 (w) := (h 1 +c 1 1 ) 3 2 (h 2 +c 2 2 ) 1 2 (h 1 +c 1 1 ) 1 + ! 0 @ 2 2 wM( 2 1 +1; 3 2 ; 1 2 w 2 ) R 1 w exp( 1 2 y 2 ) M( 2 1 ; 1 2 ; 1 2 y 2 ) 2 dy exp( 1 2 w 2 ) M( 2 1 ; 1 2 ; 1 2 w 2 ) 1 A 2 2 bM( 2 1 +1; 3 2 ; 1 2 b 2 ) R 1 b exp( 1 2 y 2 ) M( 2 1 ; 1 2 ; 1 2 y 2 ) 2 dy exp( 1 2 b 2 ) M( 2 1 ; 1 2 ; 1 2 b 2 ) + (h 1 +c 1 1 ) 1 + which shows that the conditions V 0 2 (0) = 0 and V 0 1 (b) =V 0 2 (b) = (h 1 +c 1 1 ) 3 2 (h 2 +c 2 2 ) ( 1 2 ) are satised. The relation 13.1.4 in [55] M(a;b;z) = (b) (a) e z z ab (1 +o(jzj 1 )); (5.18) 89 asjzj!1 can be used to show that V 0 1 is bounded. Finally, the following Lemma, which capitalizes on Lemma 3.2 in Weerasinghe and Mandelbaum [60], shows that V 0 1 is increasing. Lemma 5.2. For all wb, V 00 1 (w)> 0. In summary, we have providedV 1 andV 2 that solve (5.14)-(5.15). In order to establish that this V 1 and V 2 satisfy (5.13), and so solve the Bellman equation (5.10), it remains to show that (5.16) holds. This is done in the following two Lemmas. Lemma 5.3. For anyb under whichV 00 1 (b) =V 00 2 (b), we haveV 00 2 (w)> 0 for allw2 [0;b]. Lemma 5.4. There exists b ? under which V 1 (b ? ) = V 2 (b ? ) and V 00 1 (b ? ) = V 00 2 (b ? ). Fur- thermore, for such a b ? , V 0 2 is increasing on [0;b ? ]. 5.1.3 The Case that a Static Policy Solves the Brownian Control Problem The following Theorem establishes that a static priority policy is optimal when (4.1) holds, but (4.2) does not hold, either because 2 > 1 (5.19) or 1 > 2 but 1 + 2 + < 2 (h 1 +c 1 1 ) 3 (h 2 +c 2 2 ) : (5.20) 90 Theorem 5.2. Assume condition (4.1) holds, but (4.2) does not hold. Then, q ? in (5.12) solves (5.10). Proof: It is sucient to show that there exists a solution to 2 2 V 00 (w) 2 wV 0 (w) V (w) + 3 2 (h 2 +c 2 2 )w = 0 for all w 0 subject to: V is twice continuously dierentiable V 0 is bounded, has V 0 (0) = 0; and ( 1 2 )V 0 (w) (h 1 +c 1 1 ) 3 2 (h 2 +c 2 2 ) for all w 0: (5.21) This is done in the following Lemma. Lemma 5.5. There exists a unique solution to (5.21) when (4.1) holds and either (5.19) or (5.20) hold. 5.2. The Proposed Policy We rst describe our candidate for an asymptotically optimal policy. The policy we pro- pose is motivated by the solution of the associated Brownian control problem described in the previous section, and the proposed policy diers depending on the system parameters in accordance with the Brownian control problem solution. In particular, the proposed policy depends rst on which class is cheaper under the c rule; i.e., on whether or not 91 (4.1) is satised. Then, if that class is class 2, the proposed policy also depends on whether or not (4.2) that states 1 + 2 + > (h 1 +c 1 1 ) 2 (h 2 +c 2 2 ) 3 is satised, and if that class is class 1 the proposed policy also depends on whether or not condition (4.4) that states 2 + 1 + > (h 2 +c 2 2 ) 3 (h 1 +c 1 1 ) 2 is satised. The proposed policy is summarized in Table 5.1 and described in detail in Deni- tions 5.2 and 5.3 below. In this paragraph, we give a brief overview of the proposed policy. The proposed policy is exactly the queue 1 based threshold policy stated in [5] when condition (4.1) holds and condition (4.2) does not hold. In the case that both conditions (4.1) and (4.2) hold, the proposed policy is a dynamic priority control policy determined by two thresholds, one on the workload and one on queue 1, which we state in Denition 5.2. In the case that condition (4.1) does not hold and also condition (4.4) does not hold, the proposed policy is the static priority policy suggested by the c rule. Finally, in the case that condition (4.1) does not hold and condition (4.4) does hold, the proposed policy is a dynamic priority control policy determined by one threshold on the workload process, which we state in Denition 5.3. Note that in all cases the threshold placed on the workload level is of order r, while the threshold placed on queue 1 is of 92 order logr. The fact that the threshold placed on queue 1 is of order logr is consistent with the policy in [5]. (h 2 +c 2 2 ) 3 < (h 1 +c 1 1 ) 2 (h 2 +c 2 2 ) 3 (h 1 +c 1 1 ) 2 Condition (4.2) holds. Condition (4.4) holds. Workload and queue 1 based (M r ;L r ) M r threshold policy given in Denition 5.3 threshold policy, given in Denition 5.2 M r is of order r. M r is of order r; L r is of order logr. Condition (4.2) does not hold. Condition (4.4) does not hold. Queue 1 based L r threshold policy given in [5] Static priority policy suggested by L r is of order logr. the c-rule, in which server 2 gives priority to class 2 over class 1. Table 5.1: Our proposed policy. Denition 5.2. (M r ;L r ) threshold policy Let b ? be as dened in Theorem 5.1. For each r 1, let M r = [rb ? ], and let L r = [c logr] for c a positive, nite constant. In the rth system, the dynamic priority threshold control policy is described as follows: server 1 operates whenever possible, or, equivalently, server 1 is never idle when there are jobs in buer 1 or at server 1. When the workload process W r exceeds the threshold level M r , server 2 gives preemptive resume priority to class 2 jobs over class 1 jobs. In particular, when the workload process reaches M r + 1 from below, server 2 immediately suspends work on class 1 jobs and begins to serve class 2 jobs (if there is a suspended class 2 job, it resumes work on this; otherwise, the server starts work on the next job in buer 2). Otherwise, when the workload process reaches M r or goes below the level M r , and the class 1 queue length exceeds the level L r , server 2 suspends its service of class 2 jobs, and either resumes service of a suspended class 1 job or starts service of a new class 1 job. Finally, when the workload process is at or below the level M r , and the class 1 queue length is at or below the level L r , server 2 does not 93 serve class 1 jobs, and either serves a class 2 job if one is available or idles if there are no class 2 jobs. Denition 5.3. M r threshold policy Letb ? 2 solve (5.13), except withV 1 dened on [0;b ? 2 ] (instead of [b ? 2 ;1)) andV 2 dened on [b ? 2 ;1) (instead of [0;b ? 2 ]). The conditions remain the same except that V 0 1 is required to be increasing on [0;b ? 2 ] and V 0 2 is required to be increasing on [b ? 2 ;1). For each r 1, let M r = [rb ? 2 ]. In the rth system, the dynamic priority threshold control policy is described as follows: server 1 operates whenever possible, or, equivalently, server 1 is never idle when there are jobs in buer 1 or at server 1. When the workload process W r exceeds the threshold level M r , server 2 gives preemptive resume priority to class 1 jobs over class 2 jobs. In particular, when the workload process reaches M r + 1 from below, server 2 immediately suspends work on class 2 jobs and begins to serve class 1 jobs (if there is a suspended class 1 job, it resumes work on this; otherwise, the server starts work on the next job in buer 1). Otherwise, when the workload process reaches M r or goes below the level M r , server 2 suspends its service of any class 1 jobs, and either resumes service of a suspended class 2 job, or starts service of a new class 2 job. This is provided there is such a job in class 2; if not, server 2 serves a class 1 job if there is one and idles if there is no class 1 job. Remark 5.1. The purpose of the policies in Denitions 5.2 and 5.3 is to ensure that, in the heavy trac limit, all of the jobs in the system reside in the cheaper buer when the workload process is below the threshold and in the other buer otherwise. The purpose of having a threshold on queue 1 in the policy in [5] and in Denition 5.2 is to ensure that 94 server 2 does not help server 1 so much that server 1 ends up incurring too much idle time. We show weak convergence and asymptotic optimality of our proposed policy in the case that conditions (4.1) and (4.2) hold; that is, in the case that we propose the system operate in accordance with the policy described in Denition 5.2. We let T r;? denote the allocation processes associated with the use of the policy in Denition 5.2 in the rth system. In the case that condition (4.1) holds and condition (4.2) does not hold, the proofs in [5] can be adapted to the setting in which jobs renege. When condition (4.1) does not hold, the two subcases should be a little more simple, because the proposed policy either only has a single threshold, or is a static priority policy. 95 Chapter 6 Weak Convergence and Asymptotic Optimality 6.1. Weak Convergence We devote this Chapter to showing that under the threshold policy proposed in Deni- tion 5.2, when conditions (4.1) and (4.2) hold, the diusion-scaled queue-length processes in (5.4) and (5.5) weakly converge to the solution of the Brownian control problem in Denition 5.1. Specically, we prove the following Theorem. Theorem 6.1. Consider the sequence of N-systems indexed by r, where the rth system operates under the threshold policy T r;? described in Denition 5.2. Then, ^ Q r 1 ; ^ Q r 2 ; r 1 ^ I r 1 + r 2 ^ I r 2 ) q ? ~ W ? ; ~ L ? ; as r!1: The proof of Theorem 6.1 is done in several steps. First, in Section 6.1.1, we prove the following state-space collapse result. 96 Proposition 6.1. Consider the sequence ofN-systems indexed byr, where therth system operates under the threshold policy T r;? in Denition 5.2. Then, for any > 0 ^ Q r 1 1f ^ W r <bg) 0 and ^ Q r 2 1f ^ W r >b +g) 0 as r!1: Next, in Subsection 6.1.2, we prove a weak convergence result on the workload process ^ W r (t) := ^ Q r 1 (t) + r 2 r 3 ^ Q r 2 (t): In order to state this result, dene ^ X r (t) := ^ X r 1 (t) + r 2 r 3 ^ X r 2 (t) ^ V r 1 (t) := Z t 0 r 2 r 1 ^ Q r 1 (s)1f ^ W r (s)bgds ^ V r 2 (t) := Z t 0 r 2 r 3 r 2 r 2 ^ Q r 2 (s)1f ^ W r (s)>bgds ^ r (t) := Z t 0 1 ^ W r (s)1f ^ W r (s)>bgds Z t 0 r 2 r 1 ^ Q r 1 (s)1f ^ W r (s)>bgds + Z t 0 2 ^ W r (s)1f ^ W r (s)bgds Z t 0 r 2 r 3 r 2 r 2 ^ Q r 2 (s)1f ^ W r (s)bgds + 1 r Z t 0 r 2 r 1 1 n _ T r 1 (r 2 s)> 0 o +1 n _ T2(r 2 s)> 0 o + r 2 r 3 r 2 r 2 1 n _ T r 3 (r 2 s)> 0 o ds ^ V r (t) := Z t 0 1 ^ W r (s)1f ^ W r (s)>bgds + Z t 0 2 ^ W r (s)1f ^ W r (s)bgds; so that ^ W r (t) = ^ X r (t) ^ V r 1 (t) ^ V r 2 (t) + ^ r (t) ^ V r (t) + r 1 ^ I r 1 (t) + r 2 ^ I r 2 (t): The weak convergence result we prove is the following. 97 Proposition 6.2. Consider the sequence ofN-systems indexed byr, where therth system operates under the threshold policy T r;? in Denition 5.2. Let b > 0, and suppose the function q in the constraint in (5.7) is dened as q(w) = 8 > > < > > : 0; 3 2 w wb (w; 0) w>b ; so that ~ W (t) = ~ (t) Z t 0 1 ~ W (s)1f ~ W (s)>bg + 2 ~ W (s)1f ~ W (s)bgds + ~ L(t): Then, ^ W r ; ^ X r ; ^ V r 1 ; ^ V r 2 ; ^ V; ^ I r 1 ; ^ I r 2 ) ~ W; ~ ; 0; 0; 0; ~ V; 0; ~ L as r!1; where ~ V (t) = Z t 0 1 ~ W (s)1f ~ W (s)>bg + 2 ~ W (s)1f ~ W (s)bgds: Finally, in Subsection 6.1.3, we prove Theorem 6.1. 6.1.1 Proof of Proposition 6.1 We establish that ^ Q r 2 1f ^ W r >b +g) 0 as r!1: (6.1) The argument to establish ^ Q r 1 1f ^ W r <bg) 0 as r!1 has many similarities, and so we have put the details in the appendix. 98 We begin with some useful denitions r 0 := 0 r n := inffs r n1 :W r (s)rb + 1g r;k n := inffs r n :Q r k (s) = 0g;k2f1; 2g r n := inffs r n :W r (s)r(b +)g r n := inffs r n :W r (s)rbg; for n2f0; 1; 2;:::g. We make use of the shifted processes A r;n k (s) := A r k ( r n +s)A r k ( r n ); k = 1; 2; S r;n j (s) := S r j T r j ( r n +s) S r j T r j ( r n ) ; j = 1; 2; 3; R r;n k (s) := R r k ( r n +s)R r k ( r n ); k = 1; 2: Next, we dene the \good set", in which the primitive processes stay close to their mean rates. For ~ > 0 arbitrarily small, dene s r := r 3 2 b + 2 3 2 + 1 r 3 r 2 2~ ; 99 and assume ~ is small enough so that r 3 r 2 2~ > 0 for all large enough r. Let r;n := 8 > > > > > > < > > > > > > : A r;n k () r k () s r < ~ s r ; k2f1; 2g; S r;n j T r j () r j T r j () s r < ~ s r ; j2f1; 2; 3g; R r;n k () s r ~ s r 9 > > > > > > = > > > > > > ; : The proof of the following Lemma, which relies on Lemma 9 in Ata and Kumar [1], can be found in the Appendix. Lemma 6.1. P \ m r n=1 r;n C ! 0; as r!1: We are now ready to establish (6.1). Fix t> 0. Let > 0 be arbitrarily small. It is sucient to show P Q r 2 (s)1fW r (s)>r(b +)gr for some s2 [0;r 2 t] ! 0; as r!1: (6.2) Let m r := r 1 + r 2 r 3 r 2 + 2~ r 2 t: Observe that P Q r 2 (s)1fW r (s)>r(b +)gr for some s2 [0;r 2 t] P r m rr 2 t +P r m r >r 2 t;Q r 2 (s)r for some s2[ m r n=1 [ r n ; r n ]\ [0;r 2 t] P r m rr 2 t +P \ m r n=1 r;n C +P r m r >r 2 t; \ m r n=1 r;n ;Q r 2 (s)r for some s2[ m r n=1 [ r n ; r n ]\ [0;r 2 t] : 100 We now argue that the rst term on the right-hand side in the above expression converges to 0. Since every \up" excursion from rb torb + 1 requires at least one arrival P r m rr 2 t P A r 1 r 2 t + r 2 r 3 A r 2 r 2 t m r : Furthermore, P A r 1 (r 2 t) + r 2 r 3 A r 2 (r 2 t)m r P A r 1 (r 2 t) ( r 1 + ~ )r 2 t +P A r 2 (r 2 t) ( r 2 + ~ )r 2 t ; and, assuming ~ is also small enough so that ~ < min 2 1 1 ; 2 1 3 , by (190) in [5], P A r k (r 2 t)> ( r k + ~ )r 2 t exp k r 2 t 1 a;? k 1 k 1 1 + ~ =(3 k ) ; k2f1; 2g; where a;? k is the Legendre-Fenchel transform of a k ; i.e., a;? k (x) := sup l2< (lx a k (l)); x2<; k2f1; 2g: Hence P r m rr 2 t ! 0; as r!1: Lemma 6.1 establishes that P \ m r n=1 r;n C ; as r!1: 101 Hence to show (6.2), it is sucient to establish that P r m r >r 2 t; \ m r n=1 r;n ;Q r 2 (s)r for some s2[ m r n=1 [ r n ; r n ]\ [0;r 2 t] ! 0; as r!1: For this, rst note that P r m r >r 2 t; \ m r n=1 r;n ;Q r 2 (s)r for some s2[ m r n=1 [ r n ; r n ]\ [0;r 2 t] m r X n=1 P r m r >r 2 t; r;n ;Q r 2 (s)r for some s2 [ r n ; r n ]\ [0;r 2 t] : Assume we can show that P r m r >r 2 t; r;n ;Q r 2 (s)r for some s2 [ r n ; r n ]\ [0;r 2 t] (6.3) P 0 B B @ r m r >r 2 t; r;n ; r;2 n < r n ; r;2 n r n +s r ; r;2 n r n ; Q r 2 (s)r for some s2 h r;2 n ;r 2 t^ r n i 1 C C A : Note that Q r 2 (s) =Q r 2 r;2 n + A r 2 (s)A r 2 r;2 n S r 3 (T r 3 (s))S r 3 (T r 3 ( r;2 n )) (R r 2 (s)R r 2 ( r n )) for s r;2 n : This follows exactly the evolution dynamics of an underloaded GI/GI/1 queue. It is well- known that in heavy-trac an underloaded GI/GI/1 queue converges to 0 exponentially fast. We conclude that the probability on the right-hand side of (6.3) converges to 0 exponentially fast, which implies m r X n=1 P r m r >r 2 t; \ m r n=1 r;n ;Q r 2 (s)r for some s2 [ r n ; r n ]\ [0;r 2 t] ! 0; as r!1; 102 and completes the proof. Hence it remains only to show that (6.3) is valid to complete the argument that ^ Q r 2 1f ^ W r >b +g) 0 as r!1. The argument that (6.3) holds. Assume that ~ has also been chosen small enough so that 1 + 2 3 4~ s r < 2 r; ~ 3 + 2 2 3 s r < 2 ; and r 1 r 1 r 1 r 2 ~ > 0; for large enough r. We rst show that on r;n \f r n <1g, r;2 n r n +s r . To see this, note that Q r 2 ( r n +s) =Q r 2 ( r n ) +A r;n 2 (s)S r;n 3 (s)R r;n 2 (s) for r n +s r;2 n : Since W r ( r n )rb + 1 + (1 + 2 = 3 ), it follows that Q r 2 ( r n ) ( 3 = 2 )(rb + 2 + 2 = 3 ). If r n +s r > r;2 n , then Q r 2 ( r n +s r ) 3 2 rb + 2 + 2 3 s r ( r 3 r 2 2~ ): This contradicts the denition of r;2 n , because the denition of s r implies that 3 2 rb + 2 + 2 3 s r ( r 3 r 2 2~ ) = 0: 103 We next show that on r;n \f r n <1g, r;2 n r n . For this, we rst show that if queue 1 is small, the workload cannot become large. Dene r n := inffu 0 :Q r 1 ( r n +u) 3~ s r g; and note that if r n > 0, for 0s r n ^s r , W r ( r n +s) = Q r 1 ( r n +s) + 2 3 Q r 2 ( r n +s) 3~ s r + 2 3 (Q r 2 ( r n ) +A r;n 2 (s)S r;n 3 (s)R r;n 2 (s)) 3~ s r + 2 3 3 2 rb + 2 + 2 3 + r 2 s r 3 s + 2~ s r ~ 3 + 2 2 3 s r +rb + 2 + 2 3 < r b + 2 + 1 r 2 + 2 3 : Also, once queue 1 becomes large, it stays large. In particular, for ss r , on r;n , Q r 1 ( r n +s) = Q r 1 ( r n + r n ) + (A r;n 1 (s)A r;n 1 ( r n )) (S r 1 (T r 1 (s))S r 1 (T r 1 ( r n ))) (R r;n 1 (s)R r;n 1 ( r n )) 3~ s r + r 1 (s r n ) r 1 (s r n ) 3~ s r > 0: Finally, for r n + s r;2 n ^ r;1 n , s s r , since Q r 1 (s) > 0 and Q r 2 (s) > 0, so that T r 1 ( r n +s) =T r 1 ( r n ) +s and T r 3 ( r n +s) =T r 3 ( r n ) +s, 104 W r ( r n +s) = W r ( r n ) +A r;n 1 (s)A r;n 1 ( r n ) (S r;n 1 (T r 1 ( r n +s))S r;n 1 (T r 1 ( r n ))) (R r;n 1 (s)R r;n 1 ( r n )) + 2 3 0 B B B B B B @ A r;n 2 (s)A r;n 2 ( r n ) (S r;n 3 (T r 3 ( r n +s))S r;n 3 (T r 3 ( r n ))) (R r;n 2 (s)R r;n 2 ( r n )) 1 C C C C C C A r b + 2 + 1 r 2 + 2 3 + r 1 (s r n ) r 1 (s r n ) + 2 3 ( r 2 (s r n ) r 3 (s r n )): Assume r is large enough so that r k k + ~ , k2f1; 2g, and r j j ~ , j2f1; 2; 3g. Then, noting that 1 1 + ( 2 = 3 ) 2 2 = 0, it follows that W r ( r n +s)r b + 2 + 1 r 2 + 2 3 4~ s r <r(b +) for large enough r. This implies that on r;n , r;2 n < r n . Finally, note that r;2 n < n is not an additional constraint. 6.1.2 Proof of Proposition 6.2 The proof of Proposition 6.2 requires the following three Lemmas, whose proofs can all be found in the appendix. 105 Lemma 6.2. Consider the sequence of N-systems indexed by r, where the rth system operates under the threshold policy T r;? in Denition 5.2. Then, ^ I r 1 ) 0 as r!1: Lemma 6.3. Consider the sequence of N-systems indexed by r, where the rth system operates under the threshold policy T r;? in Denition 5.2. Then, Q r 1 ;Q r 2 ;T r;? ; 1 r ^ I r 1 ; 1 r ^ I r 2 ! 0; 0;T ? ; 0; 0 u.o.c., a.s., as r!1: Furthermore, ^ X r ) ~ ; as r!1; where ~ is a Brownian motion with zero drift and variance 2 . Denition 6.1. In the following, a sequence of processes with paths in D m for some m 1 is called C-tight if it is tight in D m and any weak limit point of the sequence (obtained as a weak limit along a subsequence) has continuous paths almost surely. Lemma 6.4. Consider the sequence of N-systems indexed by r, where the rth system operates under the threshold policy T r;? in Denition 5.2. Then, n ^ W r ; ^ X r ; ^ r ; ^ V r 1 ; ^ V r 2 ; ^ V r ; ^ I r 1 ; ^ I r 2 o is C-tight. 106 It is also useful to dene the function m(w) := 8 > > < > > : 2 wb 1 w>b ; and note that it is equivalent to write ^ V r as ^ V r (t) = Z t 0 m ^ W r (s) ^ W r (s)ds: We are now ready to prove Proposition 6.2. It follows from Lemma 6.4 that n ^ W r ; ^ X r ; ^ r ; ^ V r 1 ; ^ V r 2 ; ^ V r ; ^ I r 1 ; ^ I r 2 o is C-tight. Lemma 6.2 establishes that ^ I r 1 ) 0 as r!1, and Lemma 6.3 establishes that ^ X r ) ~ as r!1, where ~ is a Brownian motion with zero drift and variance 2 . Let ( ~ W; ~ ; ~ E; ~ V 1 ; ~ V 2 ; ~ V; 0; ~ L) be any weak limit point of the above sequence. Then, ~ W (t) = ~ (t) + ~ E(t) + ~ V 1 (t) + ~ V 2 (t) + ~ V (t) + 2 ~ L(t): Furthermore, ~ W; 2 ~ L is a solution of the Skorokhod problem for ~ + ~ E ~ V 1 ~ V 2 ~ V , so that ~ W; 2 ~ L = (; ) ~ + ~ E ~ V 1 ~ V 2 ~ V : 107 To see this, it is necessary to rst establish that ^ W r ; r 2 ^ I r 2 is a solution to the perturbed Skorokhod problem dened in Theorem 5.1 in Williams [61] for ^ X r ^ V r 1 ^ V r 2 + ^ r ^ V r + r 1 ^ I r 1 with the in that Theorem equalling L r =r. Conditions (i), (ii), (iii)a, and (iii)b in Theorem 5.1 in [61] are immediate. For condition (iii)c, realize that the denition of the (M r ;L r ) threshold policy in Denition 5.2 implies that for any 0t 1 <t 2 , Z t 2 t 1 1 ^ Q r 1 (t)> L r r or ^ Q r 2 (t)> 0 d ^ I r 2 (t) = 0: Since W r = ^ Q r 1 + ( r 2 = r 3 ) ^ Q r 2 , for any t> 0, ^ Q r 1 (t) L r r and ^ Q r 2 (t) = 0 ^ W r (t) L r r ; from which it follows that Z t 2 t 1 1 ^ W r (t)> L r r d ^ I r 2 (t) = 0; which establishes that condition (iii)c holds. 108 LetfF t g be the ltration generated by ~ W; ~ ; ~ E; ~ V 1 ; ~ V 2 ; ~ V; 2 ~ L . Note that ^ X r is uniformly integrable, because an argument similar to that used to show (170) in [5] shows that E sup 0st ^ A r k (s) 2 ; k2f1; 2g and E sup 0st ^ S r j (s) 2 ; j2f1; 2; 3g are all uniformly bounded in r, and so uniformly integrable. Furthermore, Lemma 6.3 establishes that Q r 1 ) 0 and Q r 2 ) 0 as r!1, and so ^ N r k ( r k )! 0 in probability for k2f1; 2g, which implies uniform integrability by Theorem 5.2 in Durrett [14]. Hence ~ is a martingale relative to fF t g. Next, ~ E + ~ V 1 + ~ V 2 + ~ V is a sum of continuous, monotone,fF t g-adapted processes, and so is of nite variation. Hence ~ W is a continuous fF t g-semimartingale with quadratic variation process h ~ W i t = h ~ X i t =t;t 0. It then follows that ~ W has a continuous local time process L a () for eacha2 [0;1), and, by the occupation time formula (see, for example Corollary 1.6 in Revuz and Yor [47]) Z t 0 1f ~ W (s) =bgd h ~ W i (s) = Z 1 0 1fb = 0gL a (t)da = 0 for all t 0: We conclude that Z 1 0 1 n ~ W (s) =b o ds = 0;a:s:: (6.4) 109 We next show that ~ V 1 = ~ V 2 = 0 a.s.. Without loss of generality, by the Skorokhod representation theorem, we may replace the weak convergence along a subsequence to ( ~ W; ~ ; ~ ; ~ V 1 ; ~ V 2 ; ~ V; 0; ~ L) to a.s. u.o.c. convergence. Let > 0 and dene ^ ^ r; (t) := Z t 0 ^ Q r 2 (s)1f ^ W r (s)>b +gds ^ r; (t) := Z t 0 ^ Q r 2 (s)1fb< ^ W r (s)b +gds: Then, for t 0, a.s., lim sup r!1 ^ V r 2 (t) lim sup r!1 ^ ^ r; (t) + lim sup r!1 ^ r; (t): It follows from Proposition 6.1 that lim sup r!1 ^ ^ r; (t) = 0: Also, lim sup r!1 ^ r; (t) lim sup r!1 3 2 Z t 0 ^ W r (s)1fb< ^ W r (s)b +gds 3 2 Z t 0 ~ W (s)1f ~ W (s)2 [b;b +]gds: Hence lim sup r!1 ^ V r 2 (t) 3 2 Z t 0 ~ W (s)1 n ~ W (s)2 [b;b +] o ds: Since the left-hand side of the above inequality does not depend on , we can let# 0 to conclude that lim sup r!1 ^ V r 2 (t) 3 2 Z t 0 ~ W (s)1f ~ W (s) =bgds: 110 It now follows from (6.4) that lim sup r!1 ^ V r 2 (t) = 0; and so ~ V 2 = 0 a.s.. We omit the argument that ~ V 1 = 0 a.s. because it is similar. We now claim that also ~ E = 0 a.s.. To see this, again replace the weak convergence along a subsequence to ( ~ W; ~ ; ~ E; ~ V 1 ; ~ V 2 ; ~ V; ~ L) to a.s. u.o.c. convergence by the Skorokhod representation theorem. Then, since ^ V r 1 ! 0 and ^ V r 2 ! 0 a.s. u.o.c. as r!1; it follows that ^ Q r 1 (t)1f ^ W r (t)bg! 0 and ^ Q r 2 (t)1f ^ W r (t)>bg! 0 a.s. u.o.c. as r!1: This then implies that ^ Q r 2 1f ^ W r (t)bg! 3 2 ~ W (t)1f ^ W r (t)bg and ^ Q r 1 (t)1f ^ W r (t)>bg! ~ W (t)1f ^ W r (t)>bg a.s. u.o.c., as r!1. We can now use the triangle inequality and the denition of ^ r to conclude that ~ E = 0 a.s.. Finally, to complete the proof, x an ! such that (6.4) holds for !, and ^ W r (;!)! ~ W (;!) u.o.c. as r!1: 111 The function m is continuous except at b, and so m ^ W r (;!) !m ~ W (;!) as r!1 for alls 0 at which ~ W (s;!)6=b. From (6.4), this is true for Lebesgue-almost everywhere s2 [0;1). Hence for such ! m ^ W r (;!) !m ~ W (;!) almost everywhere on [0;1); as r!1. We can also multiply by ^ W r (;!) to nd m ^ W r (;!) ^ W r (;!)!m ~ W (;!) ~ W (;!); as r!1, almost everywhere on [0;1): It follows by bounded convergence that Z t 0 m ^ W r (s;!) ^ W r (s;!)ds! Z t 0 m ~ W (s;!) ~ W (s;!)ds u.o.c., as r!1: Hence ^ V r ! ~ V u.o.c. a.s. as r!1: To nish the proof, note that we have shown ~ W (t) = ~ (t) + Z t 0 m ~ W (s) ~ W (s)ds + 2 ~ L(t) for all t 0: (6.5) 112 It follows from a Girsanov argument that there is uniqueness in law for a solution ~ W of (6.5). Hence the entire sequencef ^ W r g converges in distribution to ~ W , and, furthermore, ^ W r ; ^ X r ; ^ V r 1 ; ^ V r 2 ; ^ V r ; ^ I r 1 ; ^ I r 2 ) ~ W; ~ ; 0; 0; 0; ~ V; 0; ~ L as r!1: We note that another strategy for obtaining the weak convergence of the diusion- scaled workload process is to modify Theorem 2.11 in Atar, Budhiraja, and Ramanan [3] to account for a discontinuity in the linear portion of the drift (in addition to the assumed discontinuity in the constant portion of the drift). 6.1.3 Proof of Theorem 6.1 Note that everywhere we apply Proposition 6.2, we let b = b ? that is dened in Theo- rem 5.1. From Proposition 6.2, ^ V r 1 ; ^ V r 2 ) (0; 0) as r!1; which implies ^ Q r 1 ()1f ^ W r ()<bg) 0 and ^ Q r 2 ()1f ^ W r ()>bg) 0 as r!1: 113 Since ^ W r (t) = ^ Q r 1 (t)1f ^ W r (t)<bg + ^ Q r 1 (t)1f ^ W r (t)bg + r 2 r 3 ^ Q r 2 (t)1f ^ W r (t)<bg + ^ Q r 2 (t)1f ^ W r (t)bg ; and also by Proposition 6.2 ^ W r ) ~ W ? as r!1; we conclude that ^ Q r 1 ()1f ^ W r ()bg ) ~ W ? ()1f ~ W ? ()bg ^ Q r 2 ()1f ^ W r ()bg ) 3 2 ~ W ? ()1f ~ W ? ()1f ~ W ? ()<bg; as r!1; by the continuous mapping theorem, because P ~ W ? (t) =b for any t> 0 = 0. Hence ^ Q r 1 ; ^ Q r 2 )q ? ~ W ? as r!1: The joint convergence ^ Q r 1 ; ^ Q r 2 ; ^ L r ) q ? ~ W ? ; ~ L ? as r!1 follows from the weak convergence of ^ W r and ^ L r in Proposition 6.2, and the fact that q ? is a deterministic function. 114 6.2. Asymptotic Optimality We prove the following Theorem. Theorem 6.2. Suppose thatfT r g is any sequence of scheduling control policies (one for each member of the sequence of N-systems). Then, when conditions (4.1) and (4.2) both hold, lim inf r!1 ^ J r (T r )E Z 1 0 e t g q ? ~ W ? (t) dt ; (6.6) whereq ? is as given in (5.11), forb ? that solves (5.13), and ~ W ? that satises the constraint in (5.7) withq =q ? . Furthermore, when therth system operates under the threshold policy in Denition 5.2, lim r!1 ^ J r (T r;? ) =E Z 1 0 e t g q ? ~ W ? (t) dt : (6.7) Theorem 6.2 shows that E Z 1 0 e t g q ? ~ W ? (t) dt is the minimum cost that can be achieved asymptotically, and that this asymptotically minimal cost is achieved by the sequence of dynamic threshold policiesfT r;? g. Thus we conclude that the sequence of threshold policiesfT r;? g is asymptotically optimal. It is useful for the proof of Theorem 6.2 to rst state two Lemmas, which are the equivalent of Lemmas 9.2 and 9.3 in [5] in our setting. In order to do this, in addition to the uid-scaled processesQ r andT r dened directly before (4.15) and (4.16), also dene 115 A r (t) := 1 r 2 A r (r 2 t) S r (t) := 1 r 2 S r (r 2 t) N r (t) := 1 r 2 N r (r 2 t) I r 1 (t) := 1 r 2 I r 1 (r 2 t) I r 2 (t) := 1 r 2 I r 2 (r 2 t): Then, it follows from equations (4.5)-(4.8) for the queue-length and idletime processes that Q r 1 (t) = A r 1 (t)S r 1 T r 1 (t) S r 2 T r 2 (t) N r 1 ( r 1 (t)) (6.8) Q r 2 (t) = A r 2 (t)S r 3 T r 3 (t) N r 2 ( r 2 (t)) (6.9) I r 1 (t) = tT r 1 (t) (6.10) I r 2 (t) = tT r 2 (t)T r 3 (t): (6.11) Lemma 6.5. LetfT r g be any sequence of scheduling control policies (one for each mem- ber of the sequence of N systems). Then, Q r ;A r ;S r ;N r ; r ;T r ;L r is C-tight. 116 Lemma 6.6. LetfT r g be a sequence of scheduling controls such that lim inf r!1 ^ J r (T r )<1: Consider a subsequencefT r 0 g offT r g along which the lim inf is achieved; i.e., lim r 0 !1 ^ J r 0 T r 0 = lim inf r!1 ^ J r (T r ): Then, as r 0 !1, Q r 0 ;A r 0 ;S r 0 ;N r 0 ; r 0 ;T r 0 ;I r 0 ) 0;e;e;e; 0;T ? ; 0 : The proofs of Lemmas 6.5 and 6.6 can be found in the appendix. Proof of Theorem 6.2: We rst show (6.6) and then show (6.7). The argument to establish (6.6). If lim inf r!1 ^ J r (T r ) =1, then (6.6) holds trivially. Hence assume lim inf r!1 ^ J r (T r )< 1. It follows from Lemma 6.5 that on any subsequencefT r 0 g offT r g along which the lim inf is achieved, Q r 0 ;A r 0 ;S r 0 ;N r 0 ; r 0 ;T r 0 ;I r 0 ) 0;e;e;e; 0;T ? ; 0 ; 117 as r 0 !1. Hence exactly as in the proof of Lemma 6.3 ^ A r 0 ; ^ S r 0 ;T r 0 ; ^ N r 0 ; r 0 ) ~ A; ~ S;T ? ; ~ N; 0 as r 0 !1. Then, ^ X r 0 1 ; ^ X r 0 2 ) ~ X as r 0 !1; where ~ X is a two-dimensional Brownian motion with drift zero that starts from the ori- gin and has diagonal covariance matrix diag 1 2 1 + 1 2 1 + 2 2 ( 1 1 ); 2 2 2 + 2 3 . Furthermore, without loss of generality, by appealing to the Skorokhod representation theorem, we may choose an equivalent distributional representation (for which we use the same symbols) such that ^ A r 0 ; ^ S r 0 ;T r 0 ; ^ N r 0 ; r 0 ! ~ A; ~ S;T ? ; ~ N; 0 u.o.c., a.s., as r 0 !1; and ^ X r 0 1 ; ^ X r 0 2 ! ~ X u.o.c., a.s., as r 0 !1: We now argue that on the subsequencer 0 (for which lim inf ^ J(T r ) = lim r 0 !1 ^ J(T r 0 )), lim r 0 !1 ^ J(T r 0 )E Z 1 0 e t g q ? ~ W ? (t) dt; (6.12) which is sucient to complete the proof. Since lim r 0 !1 ^ J(T r 0 ) <1, it follows that lim inf r 0 !1 Q r 0 k (t)<1 for all t 0 and k2f1; 2g. Also, from (5.4) and (5.5), 118 lim inf r 0 !1 r 0 1 ^ Y r 0 1 (t) + r 0 2 ^ Y r 0 2 (t) = lim inf r 0 !1 ^ Q r 0 1 (t) + Z t 0 1 lim inf r 0 !1 ^ Q r 0 1 (s) ds ~ X 1 (t)<1 lim inf r 0 !1 r 0 3 ^ Y r 0 3 (t) = lim inf r 0 !1 ^ Q r 0 2 (t) + Z t 0 2 lim inf r 0 !1 ^ Q r 0 2 (s) ds ~ X 2 (t)<1; so that ~ Q k (t) := lim inf r 0 !1 ^ Q r 0 k (t); k2f1; 2g ~ Y j (t) := lim inf r 0 !1 ^ Y r 0 j (t); j2f1; 2; 3g are nite for all t 0. Next, taking the lim inf as r 0 !1 in (5.4) and (5.5) shows that ~ Q 1 (t) = ~ X 1 (t) Z t 0 1 ~ Q 1 (s)ds + 1 ~ Y 1 (t) + 2 ~ Y 2 (t) ~ Q 2 (t) = ~ X 2 (t) Z t 0 2 ~ Q 2 (s)ds + 3 ~ Y 3 (t): Furthermore, since ^ Y r 1 = ^ I r 1 and ^ Y r 2 + ^ Y r 3 = ^ I r 2 , and ^ I r 1 and ^ I r 2 are non-decreasing pro- cesses, the limiting processes ~ Y 1 and ~ Y 2 + ~ Y 3 are non-decreasing. Hence ( ~ Y 1 ; ~ Y 2 ; ~ Y 3 ) is an admissible control for the Brownian control problem in Denition 5.1. Therefore, by Proposition 5.1, there exists an equivalent workload control problem formulation, and so Lemma 5.1 and Theorem 5.1 allow us to conclude that E Z 1 0 e t h ~ Q(t)E Z 1 0 e t g q ? ~ W ? (t) dt: 119 Finally, to establish (6.12), note that by Fatou's lemma, lim r 0 !1 ^ J(T r 0 )E Z 1 0 e t lim inf r 0 !1 h ^ Q r 0 (t)dt ; and from the denition of ~ Q, E Z 1 0 e t lim inf r 0 !1 h ^ Q r 0 (t)dt =E Z 1 0 e t h ~ Q(t)dt : The argument to establish (6.7). It follows from Theorem 6.1 that ^ Q r 1 ; ^ Q r 2 ; r 1 ^ I r 1 + r 2 ^ I r 2 ) q ? ~ W ? ; ~ L ? ; as r!1: Hence to establish (6.7), it is sucient to show that the families n R 1 0 e t ^ Q r 1 (t)dt o and n R 1 0 e t ^ Q r 2 (t)dt o are uniformly integrable. For this, since 0 ^ Q r 1 (t) + r 2 r 3 ^ Q r 2 (t) = ^ W r (t) for all t 0, it is enough to show that the familyf R 1 0 e t ^ W r (t)dtg is uniformly integrable. We rst establish an upper bound on sup 0tT ^ W r (t). Note that ^ W r (t) = ^ X r (t) + ^ r T Z t 0 r 2 r 1 ^ Q r 1 (s) + r 2 r 3 r 2 r 2 ^ Q r 2 (s) ds + r 1 ^ I r 1 (t) + r 2 ^ I r 2 (t); for ^ r T (t) := 1 r Z t 0 r 2 r 1 1 n _ T r 1 (r 2 s)> 0 o + 1 n _ T r 2 (r 2 s)> 0 o + r 2 r 3 (r 2 r 2 )1 n _ T r 3 (r 2 s)> 0 o 120 Then, w = ^ W r x = ^ X r + r 1 ^ I r 1 + ^ r T z = r 2 r 1 ^ Q r 1 + r 2 r 3 r 2 r 2 ^ Q r 2 y = r 2 ^ I r 2 satises the conditions of Theorem 7.1 with L = L r =r. (Note that ^ W r (0) = 0. Also, the fact that condition (iii)c is satised follows as in the rst paragraph of the proof of Proposition 6.2 for the perturbed, but not linearly generalized, Skorokhod problem.) Hence there exists a constant c> 0 that does not depend on r such that for any T > 0, Osc ^ W r ; [0;T ] c Osc ^ X r + ^ r T + r 1 ^ I r 1 ; [0;T ] + L r r : Since ^ W r (0) = 0 and ^ W r (t) 0 for all t 0, sup 0tT ^ W r (t) = Osc ^ W r ; [0;T ] : Also Osc ^ X r + ^ r T + r 1 ^ I r 1 ; [0;T ] 2 sup 0tT ^ X r (t) + 2 sup 0tT j^ r T (t)j + 2 r 1 sup 0tT ^ I r 1 (t): 121 Hence it follows that sup 0tT ^ W r (t) 2 sup 0tT ^ X r (t) + 2 sup 0tT j^ r T (t)j + 2 r 1 sup 0tT ^ I r 1 (t): We conclude that to establish the family n R 1 0 e t ^ W r (t)dt o is uniformly integrable, it suces to show that the families Z 1 0 e t sup 0st ^ X r (s) dt ; Z 1 0 e t sup 0st j^ r T (s)jdt ; and Z 1 0 e t sup 0st ^ I r 1 (s) dt are uniformly integrable. The family R 1 0 e t sup 0st j^ r T (s)jdt is uniformly inte- grable because ^ r T (t) t r 1 + 2 3 2 + 1 for all t 0. The fact that the family n R 1 0 e t sup 0st ^ X r (s) dt o is uniformly integrable follows because the same argu- ments in the proof of Proposition 6.2 that show ^ X r is uniformly integrable also apply here. Next, to establish that the family n R 1 0 e t sup 0st ^ I r 1 (s) dt o is uniformly inte- grable, it suces to show that Z 1 0 e 2 t E h ^ I r 1 (t) 2 i dt<1: 122 This is because E " Z 1 0 e t sup 0st ^ I r 1 (s)dt 2 # E " Z 1 0 e 2 t sup 0st ^ I r 1 (s) 2 dt # = Z 1 0 e 2 t E h ^ I r 1 (t) 2 i dt: It follows that E ^ I r 1 (t) 2 = Z r 2 t 0 P I r 1 (r 2 t)>r p s ds: Next, for r 0 := inffs 0 :Q r 1 (s)L r g, P I r 1 (r 2 t)>r p s = P I r 1 (r 2 t)>r p sjI r 1 ( r 0 )>r p s P I r 1 ( r 0 )>r p s +P I r 1 (r 2 t)>r p sjI r 1 ( r 0 )r p s P I r 1 ( r 0 )r p s P I r 1 ( r 0 )>r p s +P I r 1 (r 2 t)I r 1 ( r 0 )> 0 : Hence for t r := 8L r =( 1 1 ), E h ^ I r 1 (t) 2 i t r r 2 + Z r 2 t ( t r r ) 2 P (I r 1 ( r 0 )>t r )ds +r 2 t 2 P I r 1 (r 2 t)I r 1 ( r 0 )> 0 : The fact that the last two terms are bounded on (0;1) follows because we have shown that r 2 P (I r 1 ( r 0 )>t r ) and r 2 P I r 1 (r 2 t)I r 1 ( r 0 )> 0 both converge to 0 as r!1 in the proof of Lemma 6.2 in the Appendix. 123 Chapter 7 Proof of the Technical Lemmas This chapter is divided into three sections. Section 7.1 proves an oscillation inequality for solutions of a one-dimensional, linearly generalized, perturbed Skorokhod problem. Section 7.2 provides the details omitted from the proof of Proposition 6.1 that establish ^ Q r 1 1 n ^ W r Lgdy(t) = 0. Then, there is a constant c> 0 such that Osc(w; [0;T ]) c (Osc(x; [0;T ]) +L): Proof of Theorem 7.1: Let w;y2 D be such that w(0) =w(0) and 1. w(t) =x(t) +y(t) for all t2 [0;T ], 2. w(t) 0 for all t2 [0;T ], 3. (a) y(0) 0, (b) y is non-decreasing, (c) R [0;T] 1fw(t)>Lgdy(t) = 0. Since w and w are non-negative processes, and w(0) =w(0) = 0, Osc(w; [0;T ]) = sup 0tT w(t) Osc(w; [0;T ]) = sup 0tT w(t): Furthermore, Theorem 5.1 in [61] establishes that there exists c 0 > 0 such that Osc(w; [0;T ])c 0 (Osc(x; [0;T ]) +L): 125 Therefore, if we can show that w(t)w(t) +L for all t 0; (7.1) then it follows that Osc(w; [0;T ]) Osc (w; [0;T ]) +L c 0 (Osc(x; [0;T ]) +L) +L c 0 1 + 1 c 0 (Osc(x; [0;T ]) +L): This completes the proof, by letting the constant c =c 0 1 + 1 c 0 . The proof that (7.1) holds is by contradiction. Suppose that there exists t > 0 for which w(t) > w(t) +L. Since w(0) = w(0), and w;w2 D, there exists s2 [0;t) such that w(s)w(s) +L and w(u)>w(u) +L for s<ut. Then, w(t)w(t)L = w(s)w(s) Z t s z(v)dv + (y(t)y(s)) y(t)y(s) L y(t)y(s); sincew(s)w(s)L,z is non-negative, andy is non-decreasing. Butw(u)>w(u)+L L for s<ut implies y(t)y(s) = 0. We conclude that w(t)w(t)L 0, which is a contradiction. 126 7.2. The Details Omitted from the Proof of Proposition 6.1 In Section 6.1.1, we established that ^ Q r 2 1f ^ W r >b +g) 0 as r!1. To complete the proof of Proposition 6.1, it remains to establish that ^ Q r 1 1f ^ W r <bg) 0 as r!1: (7.2) The general outline for this argument is the same as the one that shows ^ Q r 2 1f ^ W r > b +g) 0 as r!1; however, there are changes in the details. Specically, we dene \down" excursion intervals for the workload process (rather than \up" excursion inter- vals), and then dene shifted processes based on the times at which \down" excursions begin. Let r 0 := 0 r n := inffs r n1 :W r (s)rb 1g r;1 n := inffs r n :Q r 1 (s)L r g; r n := inffs r n :W r (s)r(b)g r n := inffs r n :W r (s)rbg s r := rb 1 r 1 + r 2 r 1 3~ ; where ~ is small enough so that r 1 + r 2 r 1 3~ > 0 for all large enough r. The shifted processes A r;n k , S r;n k , and R r;n k , k2f1; 2g, j2f1; 2; 3g, and n2f1; 2; 3;:::g, and the \good set" r;n ;n2f1; 2; 3;:::g, are dened as in the rst paragraph of Section 6.1.1 (that proves Proposition 6.1). Note that in order to emphasize the similarity in proof 127 concept, we have abused notation, and re-dened r n ; r n ; r;k n ; r n ;s r ;A r;n k ;S r;n j , and R r;n k , k2f1; 2g, j2f1; 2; 3g, and n2f1; 2; 3;:::g. Fix t> 0. Let > 0 be arbitrarily small. To show (7.2), it is sucient to show that P Q r 1 (s)1fW r (s)<r(b)gr for some s2 [0;r 2 t] ! 0; as r!1: (7.3) For m r dened as in the third paragraph of Section 6.1.1, it follows that P Q r 1 (s)1fW r (s)<r(b)gr for some s2 [0;r 2 t] P r m rr 2 t +P r m r >r 2 t;Q r 1 (s)r for some s2[ m r n=1 [ r n ; r n ]\ [0;r 2 t] P r m rr 2 t +P \ m r n=1 r;n C +P r m r >r 2 t; \ m r n=1 r;n ;Q r 1 (s)r for some s2[ m r n=1 [ r n ; r n ]\ [0;r 2 t] : The argument that the rst term on the right-hand side in the above expression converges to 0 is exactly as in the fourth paragraph of Section 6.1.1. Also, Lemma 6.1 establishes that P \ m r n=1 r;n C ! 0; as r!1: Hence to show (7.3), it is sucient to establish that P r m r >r 2 t; \ m r n=1 r;n ;Q r 1 (s)r for some s2[ m r n=1 [ r n ; r n ]\ [0;r 2 t] ! 0; as r!1: 128 For this, rst note that P r m r >r 2 t; \ m r n=1 r;n ;Q r 1 (s)r for some s2[ m r n=1 [ r n ; r n ]\ [0;r 2 t] m r X n=1 P r m r >r 2 t; r;n ;Q r 1 (s)r for some s2 [ r n ; r n ]\ [0;r 2 t] : Assume we can show that P r m r >r 2 t; r;n ;Q r 1 (s)r for some s2 [ r n ; r n ]\ [0;r 2 t] (7.4) P 0 B B @ r m r >r 2 t; r;n ; r;1 n < r n ; r;1 n r n +s r ; r;1 n r n ; Q r 1 (s)r for some s2 h r;1 n ;r 2 t^ r n i 1 C C A : Then, for r large enough so that, r> 2L r , it suces to show that as r!1 m r X n=1 P Q r 1 (s) 2L r for some s2 r;1 n ;r 2 t^ r n ! 0 as r!1: (7.5) In words, it is enough to show that in the n th down-excursion interval for the workload process once queue 1 reaches L r , it stays close to it. For this, as in [5], we need to construct up-excursion intervals for buer 1 with respect to theL r threshold, and consider the behavior of queue 1 in these up-excursions intervals. Let r 0;n := inffs r;1 n :Q r 1 (s)L r g r 1;n := inffs r 0;n :Q r 1 (s) =L r + 1g r 2;n := inffs r 1;n :Q r 1 (s)L r g; 129 and dene recursively r 2m1;n := inffs r 2m2;n :Q r 1 (s) =L r + 1g; r 2m;n := inffs r 2m1;n :Q r 1 (s)L r g; for m2f2; 3;:::g. Then, [ r 2m1;n ; r 2m;n ) is the m th up-excursion interval for buer 1 (in the n th workload down-excursion interval). The arguments in [5] that establish (52) in the statement of Theorem 7.2 in their paper also establish (7.5). It remains to establish (7.4) to complete the proof. The argument to establish (7.4). First, we argue by contradiction that on r;n \f r n <1g, r;1 n r n +s r . Suppose not; that is, suppose that on r;n , r;1 n > r n +s r . Then, since server 2 is processing jobs from buer 1 during the interval [ r n ; r n +s r ], for any s2 [ r n ; r n +s r ], Q r 1 ( r n +s r )Q r 1 ( r n ) +A r;n 1 (s)S r;n 1 (s)S r;n 2 (s)<rb 1 + ( r 1 r 1 r 2 + 3~ )s r = 0; where the last equality follows directly from the denition of s r . This is a contradiction. Hence we conclude that on r;n , r;1 n r n +s r . Next, we show that on r;n \f r n <1g, r;1 n r n . For this, it suces to show that on r;n , for anys2 [0;s r ] such that r n +s r;1 n , it is also true thatW ( r n +s)rbr. In words, in the n th workload down-excursion interval and on the good set, when queue 130 1 reaches the buer threshold,L r , the workload is still above the levelrbr. Let R be small enough so that 2 + r 2 r 3 + 1 + 3 + r 2 r 3 ~ + 1 + r 2 r 3 R s r <r for all large enough r. For any s2 [0;s r ], W r ( r n +s) rb 2 r 2 r 3 + (A r;n 1 (s)S r;n 1 (s)S r;n 2 (s)R r;n 1 (s)) + r 2 r 3 (A r;n 2 (s)R r;n 2 (s)) rb 2 r 2 r 3 j r 1 r 1 r 2 + r 2 r 3 r 2 js r ((3 + r 2 r 3 )~ + (1 + r 2 r 3 ) R )s r : Since r 1 r 1 r 2 + r 2 r 3 r 2 = r 2 r 1 r 1 r 2 1 1 2 + r 2 r 3 2 3 ; it follows from Assumption 4.1 that for large enough r s r r 2 r 1 r 1 r 2 1 1 2 + r 2 r 3 2 3 < 1: 131 Hence for all s2 [0;s r ] W r ( r n +s) rb 2 + r 2 r 3 + 1 3 + r 2 r 3 ~ + 1 + r 2 r 3 R s r > rbr; and so (7.4) holds. 7.3. Proofs of Lemmas Proof of Lemma 5.1: Let4 ~ L(t) := ~ L(t) ~ L(t ) and4V ( ~ W (t)) := V ( ~ W (t))V ( ~ W (t )). Also let ~ L C (t) := ~ L(t) P 0st 4 ~ L(s) denote the continuous part of ~ L. It follows from the generalized Ito formula (see, for example Section 4.7 in [24]) that V ( ~ W (t)) = V (0) + Z t 0 2 2 V 00 ( ~ W (s))f(q( ~ W (s)))V 0 ( ~ W (s)) ds + Z t 0 V 0 ( ~ W (s))d ~ (s) + Z t 0 V 0 ( ~ W (s))d ~ L(s) + X 0st 4V ( ~ W (s)): Also, integration by parts (see, for example, Section 4.8 in [24]) shows that e t V ( ~ W (t)) = V (0) + Z t 0 e s 2 2 V 00 ( ~ W (s))f(q( ~ W (s)))V 0 ( ~ W (s)) V ( ~ W (s)) ds + Z t 0 e s V 0 ( ~ W (s))d ~ (s) + Z t 0 e s V 0 ( ~ W (s))d ~ L(s) + X 0st e s 4V ( ~ W (s)): 132 Since V 0 is bounded, the stochastic integral is a martingale, and so taking expectations shows that Ee t V ( ~ W (t)) = Z t 0 e s 2 2 V 00 ( ~ W (s))f(q( ~ W (s)))V 0 ( ~ W (s)) V ( ~ W (s)) ds (7.6) +V (0) + Z t 0 e s V 0 ( ~ W (s))d ~ L(s) + X 0st e s 4V ( ~ W (s)): We now argue that it follows from (7.6) that E Z 1 0 e t g(q( ~ W (t)))dtV (0): Since (5.8) holds, (7.6) implies that Ee t V ( ~ W (t)) +E Z t 0 e s g(q( ~ W (s)))ds V (0) +E Z t 0 e s V 0 ( ~ W (s))d ~ L(s) +E X 0st e s 4V ( ~ W (s)): Since V is non-increasing, V 0 (x) 0 for all x 0, and so Ee t V ( ~ W (t)) +E Z t 0 e s g(q( ~ W (s)))dsV (0) +E X 0st e s 4V ( ~ W (s)) 133 Next, since ~ W (s ) = ~ W (s)4 ~ L(s) and V 0 (x) 0 for all x 0, 4V ( ~ W (s)) = V ( ~ W (s))V ( ~ W (s)4 ~ L(s)) = Z ~ W(s) ~ W(s)4 ~ L(s) V 0 (x)dx 0; and so Ee t V ( ~ W (t)) +E Z t 0 e s g(q( ~ W (s)))dsV (0): Finally, the fact that V 0 is bounded implies that V cannot grow faster than linearly, and so Ee t V ( ~ W (t))! 0 as t!1. Hence, letting t!1 on both sides of the above inequality and applying dominated convergence shows E Z 1 0 e s g(q( ~ W (s)))dsV (0): To complete the proof, it is sucient to show that V (0) =E Z 1 0 e t g(q ? ( ~ W ? (t)))dt: (7.7) For this, note that when ~ L = ~ L ? , ~ W = ~ W ? , and q = q ? , it follows from (5.8) and (5.9) that (7.6) becomes Ee t V ( ~ W ? (t)) +E Z t 0 e s g(q( ~ W ? (s)))ds =V (0): 134 As in the previous paragraph, Ee t V ( ~ W ? (t))! 0 as t!1, and so letting t!1 on both sides of the above equality establishes (7.7). Proof of Lemma 5.2: The key idea of the proof is to note that the unique solution to (5.15) can also be written in terms of Q 1 that satises 2 2 Q 00 1 (w) 1 wQ 0 1 (w) ( 1 + )Q 1 (w) = 0 for all w< 0 Q 1 (0) = 1 and lim w!1 Q 1 (w) = 0: Let U 1 (w) = V 0 1 (w) for all w 0. Then, to solve (5.15), it is enough to nd a V 1 that satises 2 2 U 00 1 (w) 1 wU 0 1 (w) ( 1 + )U 1 (w) + (h 1 +c 1 1 ) = 0 for all wb subject to :U 1 (b) = (h 1 +c 1 1 ) 3 2 (h 2 +c 2 2 ) 1 2 U 1 is bounded and non-decreasing. Such a U 1 can be written as, U 1 (w) = h 1 +c 1 1 1 + +Q 0 (w); where 135 Q 0 (w) = Q 1 (w) Q 1 (b) " (h 1 +c 1 1 ) 3 2 (h 2 +c 2 2 ) 1 2 h 1 +c 1 1 1 + # solves the homogeneous ode 2 2 Q 00 0 (w) 1 wQ 0 0 (w) ( 1 + )Q 0 (w) = 0 subject to :Q 0 (b) = (h 1 +c 1 1 ) 3 2 (h 2 +c 2 2 ) 1 2 h 1 +c 1 1 1 + Then, Q 0 0 (w) =Q 0 1 (w)> 0 for all wb because Q 0 1 (w)> 0 by Lemma 3.2 in [60] and (h 1 +c 1 1 ) 3 2 (h 2 +c 2 2 ) 1 2 h 1 +c 1 1 1 + = 2 (h 1 +c 1 1 ) 1 3 2 (h 2 +c 2 2 ) + [(h 1 +c 1 1 ) 3 2 (h 2 +c 2 2 )] ( 1 2 )( 1 + ) < 0 by (4.2). Proof of Lemma 5.3: Suppose we can establish that d 1 (b)> 0. Then, since V 2 satises (7.25), V 00 2 (0) = 2 2 V 2 (0) = 2 2 d 1 (b)> 0: Also, by Lemma 1, V 00 1 (b) = V 00 2 (b) > 0. Now, let F (w) = V 00 2 (w) and note that F must satisfy 2 2 F 00 (w) 2 wF 0 (w) ( + 2 2 )F (w) = 0 for all w2 [0;b]: (7.8) 136 The argument that F (w) > 0 for all w 2 [0;b] is by contradiction. Let F (c) := min w2[0;b] F (w), and suppose that F (c) 0. Since F (0) > 0 and F (b) > 0, it must be that c2 (0;b). Then, F 0 (c) = 0 and F 00 (c) > 0 because c is a minimum. However (7.8) implies that 2 2 F 00 (c) = ( + 2 2 )F (c) 0; which is a contradiction. Hence, F (c) 0, and we conclude V 00 2 (w) 0 for all w2 [0;b]. To complete the proof, it remains to show d 1 (b)> 0. For this, set h(b) := exp( 2 2 b 2 ) M( 2 2 ; 1 2 ; 2 2 b 2 ) + 2 2 bM 2 2 + 1; 3 2 ; 2 2 b 2 Z b 0 exp( 2 2 y 2 ) M( 2 2 ; 1 2 ; 2 2 y 2 ) 2 dy: Then, d 1 (b) := (h 1 +c 1 1 ) 3 2 (h 2 +c 2 2 ) 1 2 3 2 (h 2 +c 2 2 ) 2 + + 3 2 (h 2 +c 2 2 ) 2 + h(b) 2 2 bM( 2 2 + 1; 3 2 ; 2 2 b 2 ) Note that h(0) = 1 and by (4.1), (h1 +c11) 3 2 (h2 +c22) 12 3 2 (h2 +c22) 2 + + 3 2 (h2 +c22) 2 + h(0) = (h1 +c11) 3 2 (h2 +c22) 12 > 0 Furthermore, from (5.17) h 0 (b) = 2 2 2 b exp( 2 2 b 2 ) M( 22 ; 1 2 ; 2 2 b 2 ) + 2 2 (1 + 2 3 b + 2 2 2 )M( 2 2 + 1; 3 2 ; 2 2 b 2 ) Z b 0 exp( 2 2 y 2 ) M( 22 ; 1 2 ; 2 2 y 2 ) 2 dy > 0; 137 which implies d 1 (b)> 0 for all b> 0. Proof of Lemma 5.4: First note that by Lemma 5.3, it is sucient to show there existsb ? under whichV 1 (b ? ) = V 2 (b ? ) and V 00 1 (b ? ) =V 00 2 (b ? ). We next argue that it is sucient to show there existsb ? under whichV 1 (b ? ) =V 2 (b ? ). For this, note that sinceV 1 andV 2 solve (5.14) and (5.15), for anyb under whichV 1 (b) = V 2 (b), 2 2 V 00 2 (b) = b 2 V 0 2 (b) 3 2 (h 2 +c 2 2 ) V 2 (b) = b 2 (h 1 +c 1 1 ) 3 2 (h 2 +c 2 2 ) 1 2 3 2 (h 2 +c 2 2 ) ! V 2 (b) = b 1 (h 1 +c 1 1 ) 3 2 (h 2 +c 2 2 ) 1 2 (h 1 +c 1 1 ) ! V 1 (b) = b 1 V 0 1 (b) (h 1 +c 1 1 ) V 1 (b) = 2 2 V 00 1 (b): 138 Finally, to complete the proof, we show that there existsb under whichV 1 (b) =V 2 (b). It follows by evaluating the explicit expressions for V 1 and V 2 at b that V 2 (b) = 2 2 1 b M( 2 2 ; 1 2 ; 2 2 b 2 ) M( 2 2 + 1; 3 2 ; 2 2 b 2 ) 2 6 6 6 6 6 6 6 6 6 6 6 4 3 2 (h 2 +c 2 2 ) 2 + exp( 2 2 b 2 ) M( 2 2 ; 1 2 ; 2 2 b 2 ) 0 B B B B B B B B B B B @ 1 3 2 (h 2 +c 2 2 ) 2 (h 1 +c 1 1 ) + [ 3 2 (h 2 +c 2 2 ) (h 1 +c 1 1 )] ( 1 2 )( 2 + ) 1 C C C C C C C C C C C A 3 7 7 7 7 7 7 7 7 7 7 7 5 V 1 (b) = 2 2 1 b 1 3 2 (h 2 +c 2 2 ) 2 (h 1 +c 1 1 )+ [ 3 2 (h 2 +c 2 2 )(h 1 +c 1 1 )] ( 1 2 )( 1 + ) M( 2 1 +1; 3 2 ; 1 2 b 2 ) M( 2 1 ; 1 2 ; 1 2 b 2 ) exp( 1 2 b 2 ) bM( 2 1 ; 1 2 ; 1 2 b 2 ) 2 R 1 b exp( 1 2 y 2 ) M( 2 1 ; 1 2 ; 1 2 y 2 ) 2 + h 1 +c 1 1 1 + b: Since 3 2 (h 2 +c 2 2 ) 2 + 1 3 2 (h 2 +c 2 2 ) 2 (h 1 +c 1 1 ) + [ 3 2 (h 2 +c 2 2 ) (h 1 +c 1 1 )] ( 1 2 )( 2 + ) = (h 1 +c 1 1 ) 3 2 (h 2 +c 2 2 ) 1 2 < 0; it follows that V 2 (b)!1 as b# 0. Since also lim b#0 V 1 (b) is nite, it follows that for small b, V 1 (b)>V 2 (b): 139 Next, it follows from (5.18) that V 1 (b) b ! h 1 +c 1 1 1 + and V 2 (b) b ! 3 2 (h 2 +c 2 2 ) 2 + : Since from (4.2) 3 2 (h2 +c22) 2 + h1 +c11 1 + = 1 3 2 (h2 +c22)2(h1 +c11) + [ 3 2 (h2 +c22) (h1 +c11)] (12)(2 + ) > 0 for large b, V 1 (b)<V 2 (b): Therefore, sinceV 1 andV 2 are continuous, there must existsb ? such thatV 1 (b ? ) =V 2 (b ? ). Proof of Lemma 5.5: We construct a solution to (5.21) in terms of Q 1 that satises 2 2 Q 00 1 (w) 2 wQ 0 1 (w) ( 2 + )Q 1 (w) = 0 for all w< 0 Q 1 (0) = 1 and lim w!1 Q 1 (w) = 0: Lemma 3.2 in [60] guarantees that Q 1 exists and is unique. We can also conclude from Lemma 3.2 in [60] that 0Q 1 (w) 1 for all w< 0. Let U(w) :=V 0 (w) for all w 0. Then, to solve (5.21) it is enough to nd a U that satises 140 2 2 U 00 (w) 2 wU 0 (w) ( 2 + )U(w) + 3 2 (h 2 +c 2 2 ) = 0 for all w 0 subject to: U is twice continuously dierentiable U is bounded, has U(0) = 0; and U(w) (h 1 +c 1 1 ) 3 2 (h 2 +c 2 2 ) 1 2 for all w 0: Such a U can be written as U(w) = 3 2 (h 2 +c 2 2 ) 2 + +Q 0 (w); where Q 0 is a solution to the homogeneous ode, 2 2 Q 00 0 (w) 2 wQ 0 0 (w) ( 2 + )Q 0 (w) = 0; (7.9) and has Q 0 (0) = 3 2 (h 2 +c 2 2 ) 2 + : (7.10) Note that Q 0 (w) =Q 1 (w) 3 2 (h 2 +c 2 2 ) 2 + is a bounded solution to (7.9)-(7.10). Then, U is bounded and has U(0) = 0. Also, Q 0 0 (w) =Q 0 1 (w) 3 2 (h 2 +c 2 2 ) 2 + > 0 141 by Lemma 3.2 in [60], so thatQ 0 (w) is increasing inw. In the case that (5.19) holds, the proof is complete. Otherwise, in the case that (5.20) holds, it is enough to show lim w!1 U(w)< (h 1 +c 1 1 ) 3 2 (h 2 +c 2 2 ) 1 2 : This is true because lim w!1 Q 0 (w) = 0, which implies lim w!1 U(w) = 3 2 (h 2 +c 2 2 ) 2 + ; and under (5.20) (h 1 +c 1 1 ) 3 2 (h 2 +c 2 2 ) 1 2 3 2 (h 2 +c 2 2 ) 2 + = 2 (h 1 +c 1 1 ) 1 3 2 (h 2 +c 2 2 ) + [(h 1 +c 1 1 ) 3 2 (h 2 +c 2 2 )] ( 1 2 )( 2 + ) > 0: Proof of Lemma 6.1: It follows from DeMorgan's Laws that P \ m r n=1 C =P [ m r n=1 ( r;n ) C : Next, note that P [ m r n=1 ( r;n ) C m r X n=1 P ( r;n ) C ; 142 and that P ( r;n ) C 2 X k=1 P kA r;n k () r k ()k s r > ~ s r + 3 X j=1 P kS r;n j T r;n j () r j T r;n j ()k s r > ~ s r + 2 X k=1 P kR r;n k ()k s r > ~ s r : Sincem r =O(r 2 ), it is enough to show that there exists a nite constant b> 0 such that for each n2f1; 2;:::g, k2f1; 2g, and j2f1; 2; 3g, P kA r;n k () r k ()k s r > ~ s r b r 3 (7.11) P kS r;n k T r;n k () r k T r;n k ()k s r > ~ s r b r 3 (7.12) P kR r;n k ()k s r > ~ s r b r 3 : (7.13) We rst establish (7.11). Dene ~ u k =( r k ) to be the residual inter-arrival time at r n . Let v> 0 and note that if ~ u k =( r k )v, then jA r k ( r n +v)A r k ( r n ) r k vj st 1 +jA r k (v) r k vj: Since kA r;n k () r k ()k s r = sup 0vs r jA r k ( r n +v)A r k ( r n ) r k vj; 143 it follows that P kA r;n k () r k ()k s r > ~ s r P kA r;n k () r k ()k s r > ~ s r \ 1 r k ~ u k s r +P 1 r k ~ u k >s r P 1 +kA r;n k () r k ()k s r > ~ s r +P 1 r k ~ u k >s r : Also Assumption 4.3 implies that the inter-arrival and service times have exponential moments, so that we can let 1 = 2 in Lemma 9 in [1]. Then, letting c 2 (~ 1) be the constant dened in Lemma 9 in [1], which does not depend on r, it follows from Lemma 9 in [1] and the observation that ~ s r 1> (~ 1)s r for large enough r that P 1 +kA r;n k () r k ()k s r > ~ s r P kA r;n k () r k ()k s r > (~ 1)s r c 2 (~ 1) (s r ) 3 : Next, (113) in [1] and Markov's inequality shows that P 1 r k ~ u k >s r P 1 r k u k (1)>s r E u k (1) 3 r k 3 (s r ) 3 : We conclude that (7.11) holds. We next establish (7.12). Since for any s> 0, 0T r;n k (s)s, it is enough to show that P kS r;n k () r k ()k s r > ~ s r ! 0; as r!1: This follows exactly as in the preceding paragraph. 144 Finally, we show (7.13), which completes the proof. First observe that since a job in service cannot renege, P kR r;n k ()k s r > ~ s r P N k Z s r 0 r k Q r k (s)ds > ~ s r : Next note that P N k Z s r 0 r k Q r k (s)ds ! > ~ s r ! P (kA r;n k () r k ()ks r > ~ s r ) +P N k Z s r 0 r k Q r k (s)ds ! > ~ s r jkA r;n k () r k ()ks r ~ s r ! : It follows from (7.11) that P kA r;n k () r k ()k s r > ~ s r b=2 r 3 : Hence to show (7.13), we must show that P N k Z s r 0 r k Q r k (s)ds > ~ s r jkA r;n k () r k ()k s r ~ s r b=2 r 3 : 145 For this, we have that P N k Z s r 0 r k Q r k (s)ds > ~ s r jkA r;n k () r k ()k s r ~ s r P N k (s r r k A r k (s r ))> ~ s r jkA r;n k () r k ()k s r ~ s r P N k s r r 2 r 2 r k ( r k + ~ ) ! > ~ s r ! P N k s r r 2 r 2 r k ( r k + ~ ) ! > ~ s r jkN r k () ()k~ 2 s r ~ 2 s r ! +P kN r k () ()k~ 2 s r > ~ 2 s r : An argument similar to the one that established (7.11) shows that P kN r k () ()k~ 2 s r > ~ 2 s r b r 3 : Furthermore, wheneverkN r k () ()k~ 2 s r ~ 2 s r , it is also true that N r k (s) ~ s r for all s2 0; ~ 2 s r . Since s r r 2 < ~ 2 s r for all large r, P N k s r r 2 r 2 r k ( r k + ~ ) ! > ~ s r jkN r k () ()k~ 2 s r ~ 2 s r ! P (~ s r > ~ s r ) = 0: Proof of Lemma 6.2: The argument is similar to the proof of Theorem 7.1 in [5], except that we must also handle the reneging. 146 Fix t> 0. It is sucient to show that for any > 0, P (kI r 1 k r 2 t >r)! 0 as r!1: Let r 0 := inffs 0 :Q r 1 (s)L r g: As in the last paragraph of the proof of Theorem 6.2, the following inequality holds P (kI r 1 k r 2 t >r)P (I r 1 ( r 0 )r) +P I r 1 (r 2 t)I r 1 ( r 0 )> 0 We establish the stronger result that r 2 P (I r 1 ( r 0 )r)! 0 as r!1; (7.14) and r 2 P I r 1 (r 2 t)I r 1 ( r 0 )> 0 ! 0 as r!1; (7.15) because this is needed in the proof of Theorem 6.2. The argument to establish (7.14). Let ~ > 0 be such that ~ < min(( 1 1 )=12; 1 =2; 1 =2), and let r ~ 1 be such that for all rr ~ , r 1 r 1 > ( 1 1 )=2,j 1 r 1 j< ~ , andj 1 r 1 j< ~ . Dene t r := 8L r 1 1 V r := fA r 1 (t r ) ( r 1 ~ )t r ;S r 1 (t r ) ( r 1 + ~ )t r ;R r 1 (t r ) ~ t r g: 147 Then, for r large enough so that r>t r , since I r 1 (t)t for all t> 0, P (I r 1 ( r 0 )r)P (I r 1 ( r 0 )t r )P ( r 0 >t r ): Also, P ( r 0 >t r )P (A r 1 (t r )S r 1 (t r )R r 1 (t r )L r ): On the set V r , A r 1 (t r )S r 1 (t r )R r 1 (t r ) ( r 1 r 1 3~ )t r 1 1 4 t r = 2L r ; which implies P (A r 1 (t r )S r 1 (t r )R r 1 (t r )L r ;V r ) = 0: Therefore, noting that similar argument to Lemma 6.1 establishes r 2 P (V r ) C ! 0 as r!1; and that P (A r 1 (t r )S r 1 (t r )R r 1 (t r )L r )P (V r ) C ; it follows that r 2 P (A r 1 (t r )S r 1 (t r )R r 1 (t r )L r )! 0 as r!1; which establishes (7.14). 148 The argument to establish (7.15). We rst note that if server 1 incurs some idletime in [ r 0 ;r 2 t], so thatI r 1 (r 2 t)I r 1 ( r 0 )> 0, then it must be that Q r 1 (s) 1 for some s2 [ r 0 ;r 2 t]. Hence P I r 1 (r 2 t)I r 1 ( r 0 )> 0 =P Q r 1 (s) 1 for some s2 [ r 0 ;r 2 t] : We next make some useful denitions. Let r 1 := infft> r 0 :Q r 1 (s)<L r g r 2 := infft> r 1 :Q r 1 (s)L r g; and, in general, r 2n1 := infft> r 2n2 :Q r 1 (s)<L r g r 2n := infft> r 2n1 :Q r 1 (s)L r g: Dene the shifted processes A r;;n 1 (s) := A r 1 r 2n1 +s A r 1 r 2n1 S r;;n 1 (s) := S r 1 T r 1 ( r 2n1 +s) S r 1 (T r 1 ( r 2n1 )) R r;;n 1 (s) := R r 1 r 2n1 +s R r 1 r 2n1 ; 149 and the set V r;n := kA r;;n 1 () r 1 ()k 1 rL r < L r 6 ;kS r;;n 1 () r 1 ()k 1 rL r < L r 6 ;kR r;;n 1 ()k 1 rL r < L r 6 : It follows that for m r := [( r 1 + ~ )r 2 t] + 1, P Q r 1 (s) 1 for some s2 [ r 0 ;r 2 t] P 2n r 1 r 2 t +P \ m r n=1 V r;n C + P Q r 1 (s) 1 for some s2 [ r 0 ;r 2 t]; 2n r 1 >r 2 t;\ m r n=1 V r;n P 2n r 1 r 2 t +P \ m r n=1 V r;n C + m r X n=1 P Q r 1 (s) 1 for some s2 [ r 2n1 ; r 2n ]; r 2n1 r 2 t;V r;n : The bound in (74) in [5] implies that r 2 P 2n r 1 r 2 t ! 0 as r!1: Similar argument to the proof of Lemma 6.1 shows that P \ m r n=1 V r;n C ! 0 as r!1: Hence it is sucient to show as r!1 m r X n=1 P Q r 1 (s) 1 for some s2 [ r 2n1 ; r 2n ]; r 2n1 r 2 t;V r;n ! 0 150 Note that P Q r 1 (s) 1 for some s2 [ r 2n1 ; r 2n ]; r 2n1 r 2 t;V r;n P Q r 1 r 2n1 +A r;;n 1 (s)S r;;n 1 (s)R r;;n 1 (s) 1 for some s2 [0; r 2n r 2n1 ];V r;n : On the setV r;n , for any interval [ r 2n1 ; r 2n ] such thatQ r 1 (s) 1 for somes2 [ r 2n1 ; r 2n ], it is also true that r 2n r 2n1 infft 0 :A r;;n 1 (t)L r g infft 0 : r 1 t + 1 6 L r L r g = 1 r 5 6 L r : (The rst inequality follows because in the worst case scenario the queue goes all the way down to 0, and we must have L r arrivals.) Hence it is sucient to show that as r!1 m r X n=1 P Q r 1 r 2n1 +A r;;n 1 (s)S r;;n 1 (s)R r;;n 1 (s) 1 for some s2 0; 1 r 5 6 L r ;V r;n ! 0: On the set V r;n , Q r 1 r 2n1 +A r;;n 1 (s)S r;;n 1 (s)R r;;n 1 (s)L r 1 + ( r 1 r 1 )s 1 2 L r 1 2 L r 1: 151 Since 1 2 L r 1> 0 for large enough r, we conclude that, for large enough r, P Q r 1 r 2n1 +A r;;n 1 (s)S r;;n 1 (s)R r;;n 1 (s) 1 for some s2 0; 1 r 5 6 L r ;V r;n = 0 for each n2f1;:::;m r g, which completes the proof. Proof of Lemma 6.3: We rst note that iffZ r g is a sequence of processes and Z is a continuous deterministic process, thenZ r )Z asr!1 is equivalent toZ r !Z u.o.c. in probability as r!1. This is implicitly used several times in the proof below to combine statements involving convergence in distribution to deterministic processes. It follows directly from Lemma 6.2 and the fact that for all t 0,T r 1 (t) + 1 r ^ I r 1 (t) =t, that T r 1 )e and 1 r ^ I r 1 ) 0 as r!1: We next show that 1 r ^ W r ) 0 as r!1: (7.16) From (5.1)-(5.5), 1 r ^ W r (t) = 1 r ^ X r 1 (t) + r 2 r 3 1 r ^ X r 2 (t) + r 1 1 r ^ I r 1 (t) (7.17) Z t 0 r 2 r 1 Q r 1 (s) 1 r 2 1 n _ T r 1 (r 2 s)> 0 o 1 r 2 1 n _ T r 2 (r 2 s)> 0 o ds Z t 0 r 2 r 3 r 2 r 2 Q r 2 (s) 1 r 2 1 n _ T r 3 (r 2 s)> 0 o ds + r 2 1 r ^ I r 2 (t); for all t 0: 152 Then, w = 1 r ^ W r x = 1 r ^ X r 1 + r 2 r 3 1 r ^ X r 2 + r 1 1 r ^ I r 1 + 1 r 2 Z 0 r 2 r 1 2 X j=1 1 n _ T r j (r 2 u)> 0 o + r 2 r 3 r 2 r 2 1 n _ T r 3 (r 2 u)> 0 o du z = r 2 r 1 Q r 1 + r 2 r 3 r 2 r 2 Q r 2 y = r 2 1 r ^ I r 2 satises the conditions of Theorem 7.1 with L = L r =r 2 . (Note that ^ W r (0) = 0. Also, the fact that condition (iii)c is satised follows as in the rst paragraph of the proof of Proposition 6.2 for the perturbed, but not linearly generalized, Skorokhod problem.) Hence there exists a constant c> 0 that does not depend on r such that for any T > 0 Osc 1 r ^ W r ; [0;T ] c 0 B B B B B B @ Osc 0 B B B B B B @ 1 r ^ X r 1 + r 2 r 3 1 r ^ X r 2 + r 1 1 r ^ I r 1 + 1 r 2 R 0 r 2 r 1 P 2 j=1 1 n _ T r j (r 2 u)> 0 o du + R 0 r 2 r 3 r 2 r 2 1 n _ T r 3 (r 2 u)> 0 o du ; [0;T ] 1 C C C C C C A + L r r 2 1 C C C C C C A : (7.18) Since ^ W r (0) = 0 and ^ W r (t) 0 for all t 0, sup 0tT ^ W r (t) = Osc 1 r ^ W r ; [0;T ] ; and so to show 1 r ^ W r ) 0 as r!1; it is enough to show that the right-hand side of (7.18) weakly converges to 0 as r!1. 153 We know from Lemma 6.2 that r 1 1 r ^ I r 1 ) 0 as r!1, and it is obvious that 1 r 2 Z 0 r 2 r 1 2 X j=1 1 n _ T r j (r 2 u)> 0 o + r 2 r 3 r 2 r 2 1 n _ T r 3 (r 2 u)> 0 o du) 0 as r!1: Hence it remains to show that 1 r ^ X r 1 + r 2 r 3 1 r ^ X r 2 ) 0 as r!1: For this, note that by the functional strong law of large numbers and the fact that T r j (t)t, j2f1; 2; 3g, for all t 0, 1 r ^ A r 1 ; 1 r ^ A r 2 ; 1 r ^ S r 1 T r 1 ; 1 r ^ S r 2 T r 2 ; 1 r ^ S r 3 T r 3 ) (0; 0; 0; 0; 0) as r!1: Also, since Q r k (t) 1 r 2 A r k (r 2 t);k2f1; 2g, for all t 0, and P sup 0st 1 r 2 A r k (r 2 s) k t > 1 ! 0 for any t> 0;k2f1; 2g; it follows that for any t> 0 P Z t 0 r 2 r 1 Q r 1 (s) 1 r 2 1 n _ T r 1 (r 2 s)> 0 o 1 r 2 1 n _ T r 2 (r 2 s)> 0 o ds> 1 ( 1 t + 1)t ! 0 P Z t 0 r 2 r 2 Q r 2 (s) 1 r 2 1 n _ T r 3 (r 2 s)> 0 o ds> 2 ( 2 t + 1)t ! 0; asr!1. Hence we can also use the functional strong law of large numbers to conclude that 154 1 r ^ N r 1 Z 0 r 2 r 1 Q r 1 (s) 1 r 2 1 n _ T r 1 (r 2 s)> 0 o 1 r 2 1 n _ T r 2 (r 2 s)> 0 o ds ) 0 1 r ^ N r 2 Z 0 r 2 r 2 Q r 2 (s) 1 r 2 1 n _ T r 3 (r 2 s)> 0 o ds ) 0; as r!1. Also using Assumption 4.1, we conclude that 1 r ^ X r 1 (t) + r 2 r 3 ^ X r 2 ) 0 as r!1: Hence we have established (7.16). It is a direct consequence of (7.16) that Q r 1 ;Q r 2 ) 0 as r!1: Hence we can also conclude from the representation (7.17) that 1 r ^ I r 2 ) 0 as r!1. Therefore, in order to complete the argument that Q r 1 ;Q r 2 ;T r;? ; 1 r ^ I r 1 ; 1 r ^ I r 2 ) 0; 0;T ? ; 0; 0 as r!1; (7.19) it remains to show T r 2 ) 1 1 2 e and T r;? 3 ) 0 as r!1: 155 This follows because from (4.15) Q r 1 (t) + Z t 0 r 2 r 1 Q r 1 (s) 1 r 2 1 n _ T r 1 (r 2 s)> 0 o 1 r 2 1 n _ T r 2 (r 2 s)> 0 o ds = 1 r ^ A r 1 (t) 1 r ^ S r 1 T r 1 (t) 1 r ^ S r 2 T r 2 (t) + r 1 t r 1 T r 1 (t) r 2 T r 2 (t) 1 r ^ N r 1 Z t 0 r 2 r 1 Q r 1 (s) 1 r 2 1 n _ T r 1 (r 2 s)> 0 o 1 r 2 1 n _ T r 2 (r 2 s)> 0 o ds : Since Q r 1 ) 0 as r!1, it follows that Z 0 r 2 r 1 Q r 1 (s) 1 r 2 1 n _ T r 1 (r 2 s)> 0 o 1 r 2 1 n _ T r 2 (r 2 s)> 0 o ds) 0 as r!1; and 1 r ^ N r 1 Z t 0 r 2 r 1 Q r 1 (s) 1 r 2 1 n _ T r 1 (r 2 s)> 0 o 1 r 2 1 n _ T r 2 (r 2 s)> 0 o ds ) 0 as r!1: Hence r 1 e r 1 T r 1 () r 2 T r 2 ()) 0: Since T r 1 )e as r!1, we conclude that T r 2 ) 1 1 2 e: Since T r 2 +T r 3 + 1 r ^ I r 2 =e, and we have already shown 1 r ^ I r 2 ) 0 as r!1, it also follows that T r 3 ) 1 1 1 2 e as r!1: 156 Since 1 ( 1 1 )= 2 = 2 = 3 from Assumption 4.1, this completes the argument to establish (7.19). Finally, to complete the proof, it remains to show ^ X r ) as r!1. For this, rst note that the functional central limit theorem, Lemma 6.3, and the random time change theorem establish that ^ A r 1 (t) + r 2 r 3 ^ A r 2 ^ S r 1 T r 1 r 2 r 3 ^ S r 2 T r 2 r 2 r 3 ^ S r 3 T r 3 ) ~ ; as r!1: From Lemma 6.3, Q r 1 ;Q r 2 ! (0; 0) u.o.c., a.s., as r!1; and so the functional central limit theorem and random time change theorem also estab- lish that ^ N r 1 Z 0 r 2 r 1 Q r 1 (s)ds + ^ N r 2 Z 0 r 2 r 2 Q r 2 (s)ds ) 0 as r!1: Assumption 4.1 completes the proof. Proof of Lemma 6.4: Let t > 0. We rst observe that for any > 0, there exists M > 0 such that for all r large enough P sup 0st ^ W r (s)M <: (7.20) 157 To see (7.20), note that it follows from the fact that ^ W r = ^ Q r 1 + r 2 r 3 ^ Q r 2 and the represen- tations for ^ Q r 1 and ^ Q r 2 in (5.4) and (5.5) that for any s 0 ^ W r (s) = ^ X r 1 (s) + r 2 r 3 ^ X r 2 (s) + r 1 ^ I r 1 (s) + ^ r T (s) Z s 0 r 2 r 1 ^ Q r 1 (u) + r 2 r 3 r 2 r 2 ^ Q r 2 (u)du + r 2 ^ I r 2 (t); where ^ r T is as dened in the proof of Theorem 6.2. Then, w = ^ W r x = ^ X r 1 + r 2 r 3 ^ X r 2 + r 1 ^ I r 1 + ^ r T z = r 2 r 1 ^ Q r 1 + r 2 r 3 r 2 r 2 ^ Q r 2 y = r 2 ^ I r 2 satises the conditions of Theorem 7.1 with L = L r =r. (Note that ^ W r (0) = 0. Also, the fact that condition (iii)c is satised follows as in the rst paragraph of the proof of Proposition 6.2 for the perturbed, but not linearly generalized, Skorokhod problem.) Hence there exists a constant c> 0 that does not depend on r such that Osc ^ W r ; [0;T ] c Osc ^ X r 1 + r 2 r 3 ^ X r 2 + ^ r T + r 1 ^ I r 1 ; [0;t] + L r r The process ^ X r 1 + r 2 r 3 ^ X r 2 + r 1 ^ I r 1 is C-tight because Lemmas 6.2 and 6.3 establish weak convergence, and weak convergence establishes C-tightness. It is obvious that ^ r T ) 0 as r!1. Hence there exists M > 0 such that for all r large enough P c Osc ^ X r 1 + r 2 r 3 ^ X r 2 + ^ r T + r 1 ^ I r 1 ; [0;t] + L r r >M <: 158 Since ^ W r (0) = 0, sup 0st ^ W r (s) = Osc ^ W r ; [0;T ] ; and so it follows that P sup 0st ^ W r (s)M =P Osc ^ W r ; [0;T ] M <; which establishes (7.20). We next argue that the processes ^ r , ^ V r 1 , ^ V r 2 , and ^ V r are all C-tight. For this, it follows from Theorem 15.1 in [7] that it is sucient to verify the following two conditions. (i) For each > 0, there exists M such that P sup 0st max ^ r (s); ^ V r 1 (s); ^ V r 2 (s); ^ V r (s) >M ; for all large enough r. (ii) For each > 0 and > 0, there exists and 2 (0; 1) such that P sup 0t1 max w (^ r ; [t;t +]);w ^ V r 1 ; [t;t +] ;w ^ V r 2 ; [t;t +] ;w ^ V r ; [t;t +] ! for all large enough r, where for any set S [0; 1] and x2 D, w(x;S) := sup u;v2S jx(u)x(v)j: 159 Condition (i) follows in a straightforward manner from (7.20). For condition (ii), note that for any t2 [0; 1], since for large enough r, r 2 1 + 1 and r 2 r 3 r 2 r 2 < 2 3 2 + 1, max w (^ r ; [t;t +]);w ^ V r 1 ; [t;t +] ;w ^ V r 2 ; [t;t +] ;w ^ V r ; [t;t +] 2 1 + 2 3 2 + 2 sup 0st ^ W r (s) : Therefore, given > 0, for any M > 0, if we choose > 0 small enough so that 2 1 + 2 3 2 + 2 M <; then P 0 B B @ sup 0t1 max w (^ r ; [t;t +]);w ^ V r 1 ; [t;t +] ;w ^ V r 2 ; [t;t +] ;w ^ V r ; [t;t +] \ sup 0st ^ W r (s) M 1 C C A = 0: Hence for M large enough so that (7.20) holds, P sup 0t1 max w (^ r ; [t;t +]);w ^ V r 1 ; [t;t +] ;w ^ V r 2 ; [t;t +] ;w ^ V r ; [t;t +] ! 0 +P sup 0st ^ W r (s) M <; and we can conclude that condition (ii) holds. The process ^ X r is C-tight because it follows from Lemma 6.3 that ^ X r weakly converges to Brownian motion as r!1. The process r 1 ^ I r 1 is C-tight because it follows from Lemma 6.2 that ^ I r 1 ) 0 as r!1. 160 Finally, to complete the proof, it is sucient to show that sequencef ^ W r ; r 2 ^ I r 2 g is C-tight. For this, recall from the proof of Proposition 6.2 that ^ W r ; r 2 ^ I r 2 ; ^ X r ^ V r 1 ^ V r 2 + ^ r ^ V r + r 1 ^ I r 1 is a solution to the perturbed Skorokhod problem dened in Theorem 5.1 in [61], with the in that Theorem equallingL r =r. Therefore, by Theorem 5.1 in [61], for any 0t 1 <t 2 , there exists a constant c that does not depend on r such that w ^ W r ; [t 1 ;t 2 ] c w ^ X r ^ V r 1 ^ V r 2 + ^ r ^ V r + r 1 ^ I r 1 ; [t 1 ;t 2 ] + L r r w r 2 ^ I r 2 ; [t 1 ;t 2 ] c w ^ X r ^ V r 1 ^ V r 2 + ^ r ^ V r + r 1 ^ I r 1 ; [t 1 ;t 2 ] + L r r : Since the sum of C-tight sequences is again C-tight, it follows from the preceding para- graphs that the sequence n ^ X r ^ V r 1 ^ V r 2 + ^ r ^ V r + r 1 ^ I r 1 o is C-tight. Then, noting that L r =r! 0 as r!1, it is straightforward to use (7.21) to show that the conditions in Theorem 15.1 in [7] are satised to establish that the sequencef ^ W r ; r 2 ^ I r 2 g is C-tight. Proof of Lemma 6.5: It follows from the functional strong law of large numbers that A r ;S r ;N r ) (e;e;e) as r!1: 161 Also, since they correspond to cumulative allocations of time, each of the three compo- nents ofT r is uniformly Lipschitz continuous with a Lipschitz constant less than or equal to one, and this property is preserved by the uid-scaled processes T r . It now follows thatfA r ;S r ;N r ;T r g is C-tight. We now argue thatfA r ;S r ;N r ;T r ; r g is C-tight, for which it is enough to establish thatf r g is C-tight. For this, rst observe that for any 0st 1, r k (t) r k (s) r 2 r k (ts)A r k (1): The above inequality and the functional strong law of large numbers then implies that for any > 0 and k = (2 k k ) 1 lim sup r!1 P sup 0t1 j r k (t)j 2 k k lim sup r!1 P r 2 r k A r k (1) 2 k k = 0 lim sup r!1 P sup jstj k j r k (t) r k (s)j ! lim sup r!1 P r 2 r k k A r k (1) = 0: We conclude from Theorem 15.1 in [7] thatf g is C-tight. Finally, to complete the proof, note that equations (6.8)-(6.11) combined with the C-tightness established above and a random time change theorem shows that A r ;S r ;N r ; r ;T r ;L r is C-tight. 162 Proof of Lemma 6.6: The proof is very similar to the proof of Lemma 9.3 in [5]. From Lemma 6.5, it follows that n Q r 0 ;A r 0 ;S r 0 ;N r 0 ; r 0 ;T r 0 ;L r 0 o (7.21) is C-tight. Thus it suces to show that all weak limit points of this sequence are given by (0;e;e;e; 0;T ? ; 0). For this, suppose that (Q;A;S;N;;T;L) is obtained as a weak limit of (7.21) along a subsequence indexed byr 00 . Without loss of generality, by appealing to the Skorokhod representation theorem (see, for example, [19], Theorem 3.1.8), we may choose an equivalent distributional representation (for which we use the same symbols) such that all of the random processes in (7.21) indexed by r 00 , as well as the limit, are dened on the same probability space and the convergence in distribution is replaced by almost sure convergence on compact time intervals, so that Q r 00 ;A r 00 ;S r 00 ;N r 00 ; r 00 ;T r 00 ;L r 00 ! Q;A;S;N;;T;L ; u.o.c., a.s., as r 00 !1: From the rst sentence in the proof of Lemma 6.5, a.s., A =e;S =e, and N =e. It follows exactly as in (132) in [5] that Q = 0 a.s.. Hence also = 0 a.s.. Then, by letting r 00 !1 in (6.8)-(6.11), we have a.s., for each t 0, 163 0 = 1 t 1 T 1 (t) 2 T 2 (t) 0 = 2 t 3 T 3 (t) I 1 (t) = tT 1 (t) I 2 (t) = tT 2 (t)T 3 (t): This is exactly (133)-(136) in [5], and so the conclusion that T = T ? and L = 0 follows exactly as in the last three sentences of the proof of Lemma 9.3 in their paper, on page 25. 164 Bibliography [1] B. Ata and S. Kumar. Heavy trac analysis of open processing networks with complete resource pooling: Asymptotic optimality of discrete review policies. Annals of Applied Probability, 15:331{391, 2005. [2] B. Ata and M. Rubino. Dynamic control of a make-to-order parallel server system with cancellations, 2008. Forthcoming in Operations Research. [3] R. Atar, A. Budhiraja, and K. Ramanan. Deterministic and stochastic dieren- tial inclusions with multiple surfaces of discontinuity. Probab. Theory Relat. Fields, 142:249{283, 2008. [4] Barlow, R., and F. Proschan, Statistical Theory Of Reliability and Life Testing: Probability Models, Holt-Rhinehart-Winston, (1975) [5] S. L. Bell and R. J. Williams. Dynamic scheduling of a system with two parallel servers in heavy trac with resource pooling: Asymptotic optimality of a threshold policy. Annals of Applied Probability, 11:608{649, 2001. [6] S. L. Bell and R. J. Williams. Dynamic scheduling of a parallel server system in heavy trac with complete resource pooling: Asymptotic optimality of a threshold policy. Electronic Journal of Probability, 10:1044{1115, 2005. [7] P. Billingsley. Convergence of Probability Measures. John Wiley & Sons, Inc., New York, 1999. Second Edition. [8] C. W. Chan and V. F. Farias. Stochastic depletion problems: Eective myopic policies for a class of cynamic optimization problems. Mathematics of Operations Research, 35:333{350, 2009. [9] Chung, F. R. D. and Hwang, F. K., \Do stronger players win more knockout tour- naments?", J. Am. Statist. Ass., 73, 593-596, (1978) [10] Cont, R. and Tankov, P., Financial Modelling with Jump Processes, Chapman & Hall (2004) [11] David, H. A., \Tournaments and Paired comparisons", Biometrika, 46, 139-149, (1959) [12] Dobson, I., Carreras, B.A., and D.E. Newman, \A Loading Dependent Model of Probabilistic Cascading Failure, Probability in the Engineering and Informational Sciences, 19, 15-32, (2005) 165 [13] D. G. Down and M. E. Lewis. A call center model with upgrades, 2008. Working Paper. [14] R. Durrett. Probability: Theory and Examples. Wadsworth Publishing Company, Pacic Grove, CA, 1995. Second Edition. [15] Edwards, C. T., \Double-elimination tournaments: Counting and calculating.", Amer. Statis., 50, 27-33, (1996) [16] El Khadiri, M., and H. Cancela, \An Improvement to the Total Hazard Method for System Reliability Simulation", Probability in the Engineering and Informational Sciences, 10, 2, (1996) [17] Elperin, T., Gertsbakh, I., and Lomonsov, M., `Estimation of Network Reliability using Graph Evolution Models", IEEE Transactions in Reliability, 40, 5, 572-581, (1991) [18] Elperin, T., Gertsbakh, I., and Lomonsov, M., \An Evolution Model for Monte Carlo Estimation of Equilibrium Network Renewal Parameters", Probability in the Engineering and Informational Sciences, 6, 457-469, (1992) [19] S. N. Ethier and T. G. Kurtz. Markov Procdesses: Characterization and Conver- gence. Wiley, New York, 1986. [20] Fishman, G., A Monte Carlo Sampling Plan for Estimating Network Reliability, Operations Research, 34, 581-594, (1986) [21] Gertsbakh, I., and Y. Shpungin, \Combinatorial approaches to Monte Carlo Estima- tion of Network Lifetime Distribution", Applied Stochastic Models in Business and Industry, 20, 49-57, (2004) [22] Glasserman, P., \Monte Carlo Methods in Financial Engineering", Springer, 2004 [23] Glenn, W. A. , \A comparison of the eectiveness of tournaments.", Biometrika, 47, 253-262, (1960)a [24] J. M. Harrison. Brownian Motion and Stochastic Flow Systems. Krieger, Malabar, Florida, 1985. [25] J. M. Harrison. The BIGSTEP approach to ow management in stochastic process- ing networks. In F. P. Kelly, S. Zachary, and I. Ziedins, editors, Stochastic Networks: Theory and Applications, pages 57{90. Oxford University Press, 1996. [26] J. M. Harrison. Heavy trac analysis of a system with parallel servers: Asymptotic optimality of discrete-review policies. Annals of Applied Probability, 8:822{848, 1998. [27] J. M. Harrison. Brownian models of open processing networks: Canonical represen- tation of workload. Annals of Applied Probability, 10:75{103, 2000. Correction Vol. 13 (2003) 390-393. 166 [28] J. M. Harrison and M. J. Lopez. Heavy trac resource pooling in parallel-server systems. Queueing Systems, 33:339{368, 1999. [29] J. M. Harrison and J. A. Van Mieghem. Dynamic control of Brownian networks: State space collapse and equivalent workload formulations. Annals of Applied Prob- ability, 7:747{771, 1997. [30] Horen, J. , Riezman R., \Comparing draws for single elimination tournaments." Operations Research, 33,249-262, (1985) [31] Hull, J. \Options, Futures, and Other Derivatives", 3rd ed. Englewood Clis, NJ: Prentice-Hall, 1997 [32] Hwang, F. K. , \New concepts in seeding knockout tournaments.", Am. Math. Monthly, 107, 140-150, (1982) [33] D. L. Iglehart and W. Whitt. The equivalence of functional central limit theorems for counting processes and associated partial sums. Annals of Mathematical Statistics, 42:1372{1378, 1971. [34] Jun, C. H., and S. M. Ross, \System Reliability by Simulation: Random Hazards versus Importance Sampling", Probability in the Engineering and Informational Sci- ences, 6, 119-126, (1992) [35] F. P. Kelly and C. N. Laws. Dynamic routing in open queueing networks: Brownian models, cut constraints and resource pooling. Queueing Systems, 13:47{86, 1993. [36] C. N. Laws. Resource pooling in queueing networks with dynamic routing. Adv. Appl. Prob., 24:699{726, 1992. [37] Lefevre, C., \On the Outcome of a Cascading Failure Model, Probability in the Engineering and Informational Sciences, 20, 413-427, (2006) [38] Lieber, D., Rubinstein, R., and D. Elmakis, \Quick Estimation of Rare Events in Stochastic Networks," IEEE Transactions on Reliability, 46, 2, 254-265, (1996) [39] Liu, J. S., Monte Carlo Strategies in Scientic Computing, Springer, 2001 [40] Lomonsov, M., \On Monte Carlo Estimates in Network Reliability", Probability in the Engineering and Informational Sciences, 8, 245-264, (1994) [41] Lomonsov, M., and Y. Shpungin, \Combinatorics of Reliability Monte Carlo", Ran- dom Structures and Algorithms 14, 4, 329-343, (1999) [42] Marchand, E., \On the comparison between standard and knockout tournaments", The Statistician, 51, 2, 169-178, (2002) [43] Mazumdar, N., \Importance Sampling in Reliability Estimation", in R. Barlow, J. Fussell, and N. Singpurwalla eds., Reliability and Fault Tree Analysis, Philadelphia, SIAM, 153-165, (1975) 167 [44] Metwally S. and Atiya A. , \Using Brownian Bridge for Fast Simulation of Jump- Diusion Processes and Barrier Options,", The Journal of Derivatives, Vol. 10, no. 1, pp. 43-54, (2002) [45] Takayuki Osogami, Mor Harchol-Balter, and Alan Scheller-Wolf. Analysis of cycle stealing with switching times and thresholds. Performance Evaluation, 61:347{369, 2005. [46] A. D. Polyanin and V. F. Zaitsev. Handbook of Exact Solutions for Ordinary Dier- ential Equations. CRC Press, Boca Raton, Florida, 1995. [47] D. Revuz and M. Yor. Continuous Martingales and Brownian Motion. Springer, Berlin, 1999. Third Edition. [48] Ross, S. M., \A New Simulation Approach to Estimating Expected Values of Func- tions of Bernoulli Random Variables under Certain Types of Dependencies", IIE Transactions (2009) [49] Ross, S. M., Introduction to Probability Models, ninth edition, Academic Press, 2007 [50] Ross, S. M., SIMULATION, fourth ed., Academic Press, 2006 [51] Ross, S. M., An Elementray Introduction to Mathematical Finance, second ed., Cambridge University Press, 2003 [52] Ross, S. M., \A New Simulation Estimator of System Reliability," Journal of Applied Mathematics and Stochastic Analysis, 7, 331-336, (1994) [53] Ross, S. M., and Z. Schechner , \Using Simulation to Estimate First Passage Distri- bution", Management Science, 31, 2, (1985) [54] Searls, D. T., \On the probability of winning with dierent tournament procedures.", J. Amer. Stat. Assoc.,58, 1064-1081, (1963) [55] L. J. Slater. Con uent hypergeometric functions. In M. Abramowitz and I. A. Stegun, editors, Handbook of Mathematical Functions, pages 503{536. Dover, New York, 1972. [56] T. Tezcan and J. Dai. Dynamic control of N-systems with many servers: Asymptotic optimality of a static priority policy in heavy trac, 2008. Forthcoming in Operations Research. [57] A. R. Ward and P. W. Glynn. A diusion approximation for a Markovian queue with reneging. Queueing Systems, 43:103{128, 2003. [58] A. R. Ward and P. W. Glynn. Properties of the re ected Ornstein-Uhlenbeck process. Queueing Systems, 44:109{123, 2003. [59] A. R. Ward and S. Kumar. Asymptotically optimal admission control of a queue with impatient customers. Mathematics of Operations Research, 33:167{202, 2008. 168 [60] A. Weerasinghe and A. Mandelbaum. A many server controlled queueing system with impatient customers, 2008. Working Paper. [61] R. J. Williams. An invariance principle for semimartingale re ecting brownian mo- tions in an orthant. Queueing Systems, 30:5{25, 1998. 169

Abstract (if available)

Abstract This dissertation consists of two main parts: Monte Carlo Simulation and Heavy Traffic analysis. Much of the simulation part, chapter 1-3, is devoted to ways of improving simulation estimators in 3 different problems: Barrier Option Pricing, Random Knockout Tournaments, and System Reliability. Our efficiency criteria for comparing alternative estimators are variance and computing time. To achieve variance reduction, in addition to using combinations of different classical variance reduction techniques, innovative substantial sources of variance reduction are introduced by exploiting specific features of the problems under consideration. Second part of this dissertation, chapter 4-7, focuses on analyzing a stochastic control problem associated with a parallel server queueing system by employing Heavy Traffic approach. The stochastic control problem can not be solved exactly. We use heavy traffic limit theorems to derive an approximating control problem, referred to as Brownian control problem. Based on the Brownian control problem solution, we propose a control policy and then prove that our proposed policy isasymptotically optimal for the original control problem.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Computational validation of stochastic programming models and applications

PDF

Novel queueing frameworks for performance analysis of urban traffic systems

PDF

Strategic and transitory models of queueing systems

PDF

Bayesian analysis of stochastic volatility models with Levy jumps

PDF

Discounted robust stochastic games with applications to homeland security and flow control

PDF

Stochastic and multiscale models for urban and natural ecology

PDF

Model, identification & analysis of complex stochastic systems: applications in stochastic partial differential equations and multiscale mechanics

PDF

Novel techniques for analysis and control of traffic flow in urban traffic networks

PDF

Learning and control in decentralized stochastic systems

PDF

On the interplay between stochastic programming, non-parametric statistics, and nonconvex optimization

PDF

Queueing loss system with heterogeneous servers and discriminating arrivals

PDF

On practical network optimization: convergence, finite buffers, and load balancing

PDF

The fusion of predictive and prescriptive analytics via stochastic programming

PDF

Provable reinforcement learning for constrained and multi-agent control systems

PDF

The next generation of power-system operations: modeling and optimization innovations to mitigate renewable uncertainty

PDF

Active state tracking in heterogeneous sensor networks

PDF

Learning and decision making in networked systems

PDF

Integrated control of traffic flow

PDF

Modern nonconvex optimization: parametric extensions and stochastic programming

PDF

Some bandit problems

Asset Metadata

Creator Ghamami, Samim (author)

Core Title Stochastic models: simulation and heavy traffic analysis

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Industrial and Systems Engineering

Publication Date 08/17/2009

Defense Date 07/02/2009

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag barrier options,Brownian motion,efficient simulation,Gibbs sampler,heavy traffic analysis,jump-diffusion,Markov chain Monte Carlo,network reliability,OAI-PMH Harvest,queueing theory,stochastic control,stochastic processes limits,variance reduction,weak convergence

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Ross, Sheldon M. (committee chair), Moore, Jame E. (committee member), Ward, Amy R. (committee member)

Creator Email ghamami@usc.edu,samim_ghamami@yahoo.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-m2569

Unique identifier UC1197293

Identifier etd-Ghamami-3182 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-255496 (legacy record id),usctheses-m2569 (legacy record id)

Legacy Identifier etd-Ghamami-3182.pdf

Dmrecord 255496

Document Type Dissertation

Rights Ghamami, Samim

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Repository Name Libraries, University of Southern California

Repository Location Los Angeles, California

Repository Email cisadmin@lib.usc.edu