Class 7 OutlineWhy Variance Reduction?Why Variance Reduction?The WallCommon Random NumbersFormal CRN ArgumentAntithetic Variables (AV)AV ImplementationControl VariatesVariance Reduction: ResultsSimulation-Based OptimizationNewsvendor ModelNewsvendor Model ParametersModel DerivationNewsvendor FormulaModule Wrap-UpFollow-up ClassesModule Wrap-Up2002 - Jérémie GallienClass 7 OutlineAdvanced Topics in Simulation:1. Variance ReductionA. Common Random NumbersB. Antithetic VariablesC. Control Variates2. Simulation-Based Optimization2002 - Jérémie GallienWhy Variance Reduction?• Example: Rare Event Probability (Insurance, Risk Analysis, Reliability, Public Health, etc.)• Simulation Model for Catastrophic Event C with p < 1/1M (exact value unknown)How many trials to estimate p with 1% accuracy?• Remember:• Here (UB – LB) / E[C] < 1% Æ n > 2.6 x 1011!!!Yn  z1−/2S2nn2002 - Jérémie GallienWhy Variance Reduction?• Cost estimate within 1% for policy RCNC2 in Ontario Gateway Æ 200,000 iterations !!!• Instead of brute-force sample size approach, let’s try to reduce the estimator variance…Yn  z1−/2S2nnIdea:Reducethis term!2002 - Jérémie GallienThe WallRelative Estimation Accuracy0.00%5.00%10.00%15.00%20.00%25.00%30.00%35.00%40.00%45.00%50.00%3100320033003400350036003700380039003Sample SizeCI Width / Estimate2002 - Jérémie GallienCommon Random Numbers• This technique (CRN) is used when comparing alternatives• Intuition? • Implication: Use the same (synchronized) seeds of random numbers during simulation runs intended to compare alternatives• Is this always beneficial?2002 - Jérémie GallienFormal CRN Argument• Want to estimate E[ g(X)-h(X) ]• Generate X1, X2, X3, …, Xn, Zi= g(Xi)-h(Xi)• Estimator Z(n) = (Z1+Z2+Z3+ …+Zn) / n Var(Z(n)) = Var( g(X)-h(X) ) / n= ( Var g(X) + Var h(X) - 2Cov[g(X),h(X)] ) / n• CRN is only a good idea when g(X) and h(X) are positively correlated!2002 - Jérémie GallienAntithetic Variables (AV)• Idea: Take the average of negatively correlated (unbiased) estimators!• What’s the intuition?2002 - Jérémie GallienAV Implementation• Want to estimate E[ h(X) ]. Let U1, U2, … , Unthe Uniform[ 0,1 ] numbers used to generate X1, … , Xn• Idea: Compute X’iusing 1 – Ui!!! Why do X’iand Xihave the same distribution?• Estimators: Z(n) = (Z1+Z2+Z3+ …+Zn) / nZ’(n) = (Z’1+Z’2+Z’3+ …+Z’n) / n with Zi= h(Xi) and Z’i= h(X’i)• Take W(n) = (Z(n) + Z’(n)) / 2 !When does it work best / worst?2002 - Jérémie GallienControl Variates• Let Y be the raw output variableLet X be some variable correlated to Y• Definition: Z = Y + c( X-E[X] )• Intuition?• Var Z = Var Y + c2 Var X + 2c Cov(X,Y), soVar Z is minimized when c = - Cov(X,Y) / Var X2002 - Jérémie GallienVariance Reduction: Results• For the reliability example, the results obtained where:EstimatorRaw Antithetic Control VariateStandard Absolute 0.144 0.100 0.023Deviation Relative 100% 69% 16%2002 - Jérémie GallienSimulation-Based Optimization• Remember Monte-Carlo framework:• The associated optimization problem is:Estimate E[ h(X) ]where X = {X1,…,Xm} is a random vector in RmMax dE[ h(X,q) ]s.t. q ∈ Fwhere X is a random vector in Rmq a vector of decision variables in RwF is a subset of Rw(feasible region)2002 - Jérémie GallienNewsvendor Model• One time decision under uncertainty• Trade-off: Ordering or producing– too much (waste, salvage value < cost) versus– too little (excess demand is lost)• Examples:– Restaurant;– Fashion;– High Tech;– Capacity and Inventory decisions…2002 - Jérémie GallienNewsvendor Model Parameters• q = Order Quantity decision• c = Unit Cost • r = Unit Revenue parameters• s = Unit Salvage Value (r > c > s)• d = Demand (unknown) random variable2002 - Jérémie GallienModel Derivation∆ Profit:• IF d > q • IF d < qProfit:()crq−⋅()crd−⋅r -c s -cEAP:()()()()s-cqdPr-cqdP⋅≤+⋅>Incremental Analysis:q → q +1:()()csdq−⋅−+(demand > order qty) (demand < order qty)As long as the Expected Additional Profit[EAP] is positive, it is lucrative to increase q to q + 1 !!!2002 - Jérémie GallienNewsvendor Formula• Set EAP = 0 to find:()()() stocking-over ofcost stocking-under ofcost ouukkksccrcrsrcrqdP+=−+−−=−−=<qsrcr−−Demand Distribution2002 - Jérémie GallienModule Wrap-Up• Class 1 Definitions, The Simulation Process• Class 2 Monte-Carlo Theory and Applications, Crystal Ball• Class 3 Ontario Gateway: Monte-Carlo Case• Class 4 Discrete-Event Theory and Applications, Simul8• Class 5 Human Genome Project: Discrete-Event Case• Class 6 Design of Simulation Experiments• Class 7 Advanced Topics: Variance Reduction, Simulation-Based Optimization2002 - Jérémie GallienFollow-up Classes• ESD.76J / 1.019 Systems SimulationSpring, 12 credits. Comparable footprint.• 1.021J Introduction to Modeling and SimulationSpring, 12 credits. Continuous models, engineering applications.• 2.141 Modeling and Simulation of Dynamics SystemsFall, 12 credits. Advanced engineering simulation class.2002 - Jérémie GallienModule Wrap-UpIndustrial decision-making is interdisciplinary:ContextualknowledgeOperationsFinanceAccountingMethodologicalskillsProbabilityOptimizationSoftwareImplementationSolverCrystal BallSimul8…Output AnalysisSensitivity