{'fraction_non_alphanumeric': 0.05204928337943173, 'fraction_numerical': 0.011549744363423016, 'mean_word_length': 3.907939119703826, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 29, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Bayesian Risk Approach to Data-driven Stochastic Optimization: Formulations and Asymptotics

A large class of stochastic programs involve optimizing an expectation taken with respect to an underlying distribution that is unknown in practice. One popular approach to addressing the distributional uncertainty, known as the distributionally robust optimization (DRO), is to hedge against the worst case over an uncertainty set of candidate distributions. However, it has been observed that inappropriate construction of the uncertainty set can sometimes result in over-conservative solutions. To explore the middle ground between optimistically ignoring the distributional uncertainty and pessimistically fixating on the worst-case scenario, we propose a Bayesian risk optimization (BRO) framework for parametric underlying distributions, which is to optimize a risk functional applied to the posterior distribution of an unknown distribution parameter. Of our particular interest are four risk functionals: mean, mean-variance, value-at-risk, and conditional value-at-risk. To unravel the implication of BRO, we establish the consistency of objective functions and optimal solutions, as well as the asymptotic normality of objective functions and optimal values. More importantly, our analysis reveals a hidden interpretation: the objectives of BRO can be approximately viewed as a weighted sum of posterior mean objective and the (squared) half-width of the true objective's confidence interval.


Introduction
Stochastic simulation is a widely used tool for decision making under uncertainty. When we build a simulation model for a real-world system, the randomness in the system is typically modeled by random variables, which leads to what we call the "stochastic uncertainty" in the model. The distributions of these random variables, viewed as an input fed into the model, are referred to as input distributions. For example, if we model a queueing system as an M/M/1 queue, then the customer interarrival times and the service times are modeled as random variables representing stochastic uncertainty, while the exact distributions of the interarrival times and the service times are the input distributions that must be specified to drive the simulation process.
In practice, input distributions are often estimated from a finite amount of input data, and the resulting model is subject to the so-called "input uncertainty" caused by finitesample error. Therefore, the error of simulation output is a mixed outcome of both stochastic uncertainty and input uncertainty. However, in contrast to stochastic uncertainty, which can be reduced by increasing the number of simulation replications, input uncertainty is often uncontrollable due to the difficulty in acquiring additional input data. As a result, if we optimize a model built on input distributions estimated from a small number of data samples, the obtained optimal solution may have a poor performance under the true input distribution. To see the risk of ignoring input uncertainty, let us consider a newsvendor problem with lost sales, where the goal is to maximize the expected profit under stochastic demand (details can be found in section 6.3). Suppose the demand is exponentially distributed with rate θ c = 1/20, and the unit price and unit cost are 2 and 1, respectively. We draw 1,000 sets of independent and identically distributed (i.i.d.) demand data, each dataset of size 20. Then, for each dataset we use maximum likelihood to get a point estimate of θ c , and we plug the estimate into the model to solve a maximization problem. The histograms of optimal solutions and their corresponding true performances (under θ c ) are given in figure  1. Although the average optimal solution is close to the true optimal solution 13.8629, the extreme values could reach 5 and 25. Similarly, the extreme performance of solutions can be as low as 3 compared with the true optimal performance 6.1371. In this example, a decision maker may suffer a great loss if s/he blindly assumes that the estimated input distribution is accurate.  There has been an extensive study on how to quantify the impact of input uncertainty on simulation output in the presence of stochastic uncertainty. The methods proposed so far in the literature can be roughly grouped into two major categories. One is the frequentist methods that allow nonparametric input distributions and use direct or bootstrap resampling techniques to assess input uncertainty (e.g., [1,2,3]). The other is the Bayesian model averaging methods (e.g., [4,5,6]) that assume a parametric model and use the posterior distribution of unknown parameters as the sampling distribution during simulation process. In addition to these methods, [3] uses delta method to decompose the variance of simulation output into two parts corresponding to stochastic uncertainty and input uncertainty. [7] proposes a method for quickly assessing the relative contribution of each input distribution to the overall effect of input uncertainty. [8] quantifies the risk of input uncertainty via their proposed efficient estimators based on nested simulation. Recent advances in stochastic kriging [9] also gives rise to the application of meta-model assisted method to quantify input uncertainty (see, e.g., [10,11,12]) Simulation optimization under input uncertainty, on the other hand, has not been studied as extensively. Since we are concerned with the extreme performance caused by input uncertainty, perhaps the most related study is distributionally robust optimization (DRO). The idea of DRO is to optimize the worst case over a family of possible input distributions, where the crux is to construct a reasonable and tractable ambiguity set. An abundant literature is available on this topic, and we refer the reader to [13,14,15,16,17] for some reviews and recent development. In simulation literature, a few works have studied the case of a finite solution space. For example, in a situation where solutions have independent but different input distributions, [18] develops a procedure for selecting a subset of solutions whose performances are within a user-specified distance from the true optimal performance. [19] considers a case of shared input distributions across solutions and explores the impact of input uncertainty on the indifference-zone (IZ) ranking and selection method. They show that a direct application of IZ selection may provide invalid correct-selection guarantee. These works provide general guidelines for finding solutions with good average performances, but they do not consider the extreme cases of input uncertainty.
However, in many applications, finding solutions with an acceptable level of potential risk is at least as important as achieving good average performance. For example, in the aforementioned newsvendor problem with input uncertainty, a risk-averse manager would base the order quantity on both the expected profit and the risk of large loss. In other cases, such as the operation of hospital emergency rooms and large-scale electrical power systems, there is even a greater need for controlling the risk of extreme performance because the failure of these systems, no matter how unlikely it may seem, can result in disastrous consequences. Thus, a more comprehensive framework is needed to incorporate the decision maker's risk preference while accounting for the impact of input uncertainty.
In financial engineering, one way to achieve a similar goal is to optimize risk measures, such as mean-variance [20], value-at-risk (VaR) [21], and conditional value-at-risk (CVaR) [22]. In simulation literature, [23] borrows this idea and proposes risk formulations for simulation optimization under input uncertainty. They consider a continuous problem where all solutions share the same parametric input distributions with unknown parameters. A Bayesian risk formulation is proposed to minimize the risk measure of mean output with respect to (w.r.t.) the posterior distribution of the unknown parameters. The purpose of risk formulations is to adaptively model input uncertainty and provide a range of choices between optimizing the average performance (risk-neutral) and the worst case performance (risk-averse). Four risk formulations have been studied in detail: the mean, mean-variance, VaR and CVaR formulations. They show that as the size of input data goes to infinity, the objective functions of the four risk formulations converge pointwise to the true objective function in probability.
Despite the reported robustness of risk formulations in the numerical results of [23], several fundamental questions on the effectiveness of risk formulations remain to be addressed : (i) does the set of optimal solutions converge to the true optimal solutions as the data size goes to infinity? (ii) can we characterize the asymptotic distributions of the objective functions and optimal values? (iii) what is the essential meaning of optimizing risk formulations? To answer the first two questions, we borrow well-established results from Bayesian asymptotics to show the consistency of risk formulations' objective functions and optimal solutions, as well as the asymptotic normality of their objective functions and optimal values. Our proof also uses a general paradigm proposed in sample average approximation (SAA) literature [24], but our detailed proof techniques are very different since SAA is under a frequentist setting.
The third question deserves the most attention, because even though the intent of risk formulations is to improve the robustness of solutions, it is not obvious at all in what sense the solutions are more robust. In this paper, we analyze the asymptotic normality of risk formulations and shed light on a rather interesting fact: the risk formulations are approximately equivalent to a weighted sum of posterior mean performance and the (squared) half-width of the true performance's confidence interval (CI), leading to an explicit tradeoff between average performance and the risk of actual performance. For example, if we apply risk formulations to the aforementioned newsvendor problem and minimize the mean and the CVaR risk measures, then the average solution of the mean formulation (13.9763) is closer to the true optimal solution (13.8629) than the CVaR formulation's average solution (9.0794). Consequently, the mean formulation outperforms the CVaR formulation in terms of the average expected profit of solutions (5.9175 vs. 5.3980). However, if we construct CIs for the true performance of the mean formulation's solutions, then the average width of the CIs is 5.5116, which is larger than the CVaR formulation's 2.7476. The reason is that although 13.8629 is the true optimal solution, 0 is the solution least affected by input uncertainty since 0 order quantity always results in 0 profit. Therefore, the CVaR formulation attempts to strike a balance between the expected profit and the risk of large loss, and is able to come up with a solution between the solution with the highest expected profit and the solution whose actual performance can be estimated with the most accuracy.
We summarize the contributions in this paper as follows.
1. We establish the consistency and asymptotic normality of risk formulations.
2. We show that when the size of input data is small, risk formulations are explicitly seeking a tradeoff between average performance and the risk of actual performance.
3. We use numerical examples to verify our results and intuition, and provide some insight into the price of robustness by comparing different risk formulations' behavior.
The rest of the paper is organized as follows. In section 2, we formally introduce risk formulations, and then lay a measure-theoretic foundation for Bayesian risk framework. In section 3 and 4, we establish the consistency and asymptotic normality of risk formulations. We reveal the underlying meaning of optimizing risk formulations in section 5. The numerical results are provided in section 6, followed by concluding remarks in section 7.

Risk Formulations under Bayesian Framework
In this section, we specify the simulation optimization problem and the type of input uncertainty we consider. Then, we formally introduce risk formulations under Bayesian framework. Since our main goal is to rigorously prove the consistency and asymptotic normality of risk formulations, it is necessary to formalize Bayesian statistics in a measure-theoretic context and explicitly construct an underlying probability space. As we shall see later, our construction of the probability space clarifies under which measure almost sure (a.s.) convergence occurs, which is crucial to characterizing a typical mode of convergence under Bayesian framework.

Risk formulations and highlights
In simulation optimization, we often want to optimize the expected performance of a simulation model. Mathematically put, the following problem is of our interest.
where x is a decision variable that takes value in X ⊂ R d , ξ is a random variable (r.v.) that takes value in R m , and h : R d × R m → R is a deterministic function capturing the system dynamics. Here ξ represents the stochastic uncertainty in the system, and sometimes its distribution may depend on x. For example, the distribution of a product's market demand may depend on the pricing of the product. We consider a less general case where all x ∈ X share the same stochastic uncertainty ξ. For instance, if we optimize the energy consumption of a building, the weather condition, modeled as a r.v., is hardly affected by the decisions we make. In simulation literature, the case of a finite X is the most often studied, where the optimization problem is called a ranking and selection problem. In this paper, we allow X to be a compact set and h may have certain structures. We also assume that we know ξ's cumulative distribution function (c.d.f.) has a parametric form F (· | θ c ), but the parameter θ c is unknown can only be estimated from a finite stream of input data ψ n = (ξ 1 , ξ 2 , . . . , ξ n ), where ξ 1 , . . . , ξ n are i.i.d. samples from F (· | θ c ). Let us define a bridging function Since the expectation is estimated by Monte-Carlo simulation, the strong law of large numbers (SLLN) guarantees that the estimate of H at a given x and θ converges a.s. to the true value as we increase simulation effort. Asymptotically speaking, it is easy to account for stochastic uncertainty via another layer of limit, but it only adds more complication than insight if our focus is on input uncertainty. For this reason, we assume that there is no stochastic uncertainty and H can be evaluated exactly. As we mentioned in section 1, it can be dangerous to blindly assume the estimate of θ c is accurate, and in many applications there is a need for incorporating a decision maker's risk preference into the formulation. In DRO, this is done by solving a min-max problem where D is an ambiguity set representing input uncertainty. In practice, D is constructed from ψ n through some statistical procedures, and the "size" of D (e.g., the diameter of D) determines the extent of conservativeness. However, constructing a D based on one's risk preference is a nontrivial task, because it is often unclear how the optimal solutions to (2.2) change as we modify D. Instead, we take a Bayesian perspective and use the posterior distribution of θ c to model input uncertainty. We minimize a risk measure ρ on H w.r.t. the posterior distribution, i.e., min where θ n is a belief r.v. following the posterior distribution of θ c constructed from ψ n . This Bayesian risk formulation is motivated by two aspects: (i) the posterior distribution adaptively models input uncertainty since it gradually concentrates around θ c as n increases; (ii) optimizing a risk measure allows some flexibility between optimizing the average performance and optimizing the worst-case scenario (e.g., VaR 100% can be viewed as the worst case). Of our particular interest are the mean formulation, the mean-variance formulation with a weight constant c > 0, the VaR formulation with risk level α ∈ (0, 1), and the CVaR formulation with risk level α ∈ (0, 1), In particular, VaR and CVaR are two commonly used risk measures in financial industry for controlling large loss. For a r.v. X, VaR α (X) is defined as the α-quantile of X, i.e., VaR α (X) := inf{t : P(X ≤ t) ≥ α}, and CVaR α (X) is defined as the following conditional tail expectation.
We remark here that CVaR is a coherent measure and it enjoys some nice properties such as convexity. The mean-variance risk measure and VaR, however, are not coherent risk measures. For more details on coherent risk measures, we refer the reader to [25]. How to solve the risk formulations is still an ongoing study. One way is optimize an sample average approximation (SAA) of the original objectives [26]. But since estimating VaR and CVaR involves rare-event simulation, applying SAA directly may not be very efficient. If H is a smooth function, an alternative is to apply stochastic gradient methods, such as Robbins-Monro algorithm [27]. The gradient of VaR and CVaR has been studied in [28] and [29]. Since part of our study is to analyze the behavior of optimal solutions under the the change of parameters such as c and α, we assume that the risk formulations can be solved to optimality without bothering the choice of algorithm.
To provide a sense of the main contribution of this paper, we highlight the following results on the asymptotic normality of risk formulations. We should point out that the convergence results presented here are in the sense of weak convergence of random functions rather than just pointwise convergence.
(i) The mean and mean-variance formulation In (i)-(iii), N stands for a normal distribution, φ and Φ are the probability density function (p.d.f.) and c.d.f. of a standard normal distribution, and σ x is a term that depends on the sensitivity of H to θ and the Fisher information ξ carries about the true parameter θ c . An immediate consequence of (i)-(iii) that we can construct CIs for the true performance of solutions. More importantly, as we will show in section 5, these results imply that the risk formulations are approximately equal to a weighted sum of posterior mean performance and the half-width of the true performance's CI.

Construction of Probability Space
To study the asymptotic properties of risk formulations, we will borrow machineries from Bayesian asymptotics. This requires us to put Bayesian statistics in a measure-theoretic context. Since we describe the uncertainty of θ c via a belief r.v. θ, it is necessary to construct a common underlying probability space (Ω, F , µ) for θ and ξ. Our construction is very similar to the standard construction in Bayesian literature (see, e.g., [30]). Let (Θ, B Θ , π) be the space in which θ takes value, where Θ ⊂ R l , B Θ is the Borel σ-algebra on Θ, and π is the prior distribution. Also let (Ξ, B Ξ , P 1 θ ) be the space in which ξ takes value, where Ξ ⊂ R m , B Ξ is the Borel σ-algebra on Ξ and P 1 θ is the distribution of ξ conditional on θ. Then the n i.i.d. copies of ξ take values in the product space (Ξ n , B n Ξ , P n θ ), where B n Ξ is the product σ-algebra B Ξ ⊗ · · · ⊗ B Ξ and P n θ is the product measure P 1 θ × · · ·× P 1 θ . Next, use Kolmogorov's extension theorem [31] to extend (Ξ n , B n Ξ , P n θ ) to a sequence space (Ξ N , B N Ξ , P θ ), where Ξ N is the space of all infinite sequences in Ξ, and B N Ξ is the σ-algebra generated by all cylinder sets of the form ξ ∈ Ξ N : (ξ 1 ,ξ 2 , . . . ,ξ n ) ∈ H , H ∈ B n Ξ , whereξ i is the ith entry in the sequenceξ. Correspondingly, P θ is the product measure that coincides with P n θ on B n Ξ , i.e., P θ ξ ∈ Ξ N : (ξ 1 ,ξ 2 , . . . ,ξ n ) ∈ H = P n θ (H), ∀H ∈ B n Ξ . We now arrive at the final step of constructing the probability space (Ω, F , µ). Let Ω = Θ × Ξ N be the sample space equipped with a σ-algebra F = B Θ ⊗ B N Ξ . We make the following assumption throughout the paper.
Under assumption 2.1, we define a measure µ on the collection of all rectangle sets in F (a semialgebra 1 ) with the following property. 3) The integrals in (2.3) are well-defined by assumption 2.1. Moreover, there exists a unique extension of µ to F by Carathéodory's extension theorem [31].
Proposition 2.1 characterizes µ for all sets in F . Loosely speaking, µ is the joint distribution of θ and an infinite sequence of i.i.d. ξ's. It is important to understand what it implies for an event to occur µ-a.s. : if C ∈ F has µ(C) = 1, then P θ (C θ ) = 1 for almost all θ ∈ Θ possibly except for a subset of measure 0 w.r.t. π. This particular mode of convergence is due to 1 = µ(C) = Θ P θ (C θ )π(dθ).

Consistency of Risk Formulations
Let P n (·) := P(θ ∈ · | ψ n ) = P(θ ∈ · | F n ) denote the posterior distribution of θ conditional on input data ψ n . For ease of notation, we also let θ n denote a belief r.v. with distribution P n . Then two of the most fundamental questions on risk formulations are (i) does the objective function ρ[H(x, θ n )] converge to the true objective function H(x, θ c ) pointwise as n → ∞? (ii) does the minimizers of ρ[H(x, θ n )] converge to the minimizers of H(x, θ c ) as n → ∞? Clearly such convergence can only be in a probabilistic sense. To answer these questions, we prove the pointwise strong consistency of ρ[H(x, θ n )] w.r.t. P θ c , and the a.s. convergence of optimal solutions w.r.t. µ. Let us begin by making the following assumption.
(ii) For all n, P n θ has a density p n θ that is B Θ × B Ξ -measurable.
(iii) For any neighborhood V ∈ B Θ of θ c , there is a sequence of uniformly consistent tests 2 of the hypothesis P θ = P θ c against the alternative P ∈ {P θ , θ ∈ V c }.
(iv) For any ǫ > 0 and any neighborhood V ∈ B Θ of θ c , V contains a subset W such that π(W ) > 0 and D KL (p θ c p θ ) < ǫ on W .
Assumption 3.1 provides sufficient conditions for the strong consistency of posterior distribution. We remark here that (ii)-(iv) are mild assumptions that are usually satisfied in practice, but (i) is often too strong because many conjugate priors, such as Gamma distribution and normal distribution, are not supported on a compact set. Therefore, our plan is to assume (i) holds throughout the derivation of theoretical results, and then test general cases in the numerical experiment. We will need the following theorem on the strong consistency of posterior distribution. Lemma 3.1 (Theorem 6.1, [30]). Suppose assumption 3.1 holds. Then for any neighborhood As is noted in [32], except for some special priors, the method used to prove lemma 3.1 is "practically the only general method for establishing consistency". The merit of having strong consistency of P n is that we can find a subsetΩ of Ω on which P n converges in distribution to the point mass on θ c . This kind of result cannot be obtained if we only have weak consistency, i.e., P n (A) → 1 in probability (P θ c ) for any neighborhood A ∈ B Θ of θ c .

Consistency of objective functions
The pointwise convergence of ρ[H(x, θ n )] to H(x, θ c ) is important, because when there is little input uncertainty, we do not want to be solving a problem that is very different from the problem we actually want to solve. It has been shown in [23] that for a fixed x, ρ[H(x, θ n )] converges to H(x, θ c ) in probability (P θ c ) as n → ∞, where the proof is based on the weak consistency of posterior distribution. We strengthen this result by proving the pointwise strong consistency of risk formulations. Our proof technique is different from the proof in [23]. Moreover, the result we obtain helps us establish the consistency of optimal solutions. Lemma 3.2. Suppose assumption 3.1 holds. Then P n ⇒ δ θ c a.s. (P θ c ), where "⇒" means convergence in distribution and δ θ c is a point mass on θ c .

Proof. By lemma 3.1 we have for any
, an open ball centered at θ c with radius 1/m. It follows that for every m, there exists Ω m ∈ B N Ξ such that P θ c (Ω m ) = 1 and P n (A m ) → 1 on Ω m . Definẽ Ω := ∩ ∞ m=1 Ω m , then P θ c (Ω) = 1 and it suffices to show P n ⇒ δ θ c onΩ. Let θ n be a belief r.v. with distribution P n . Note that for any η > 0, we can choose k > 1/η such that onΩ we have P( θ n − θ c ≤ η) ≥ P(θ n ∈ A k ) = P n (A k ) → 1, so θ n → θ c in probability and hence in distribution P θ c -a.s.. Proof. The continuity of H in θ for all x allows us to suppress x and write H(x, θ) as H(θ). We only need to prove convergence on the eventΩ constructed in the proof of lemma 3.2. Let us begin with the mean formulation. Since Θ is compact, H is bounded on Θ. Recall that θ n denotes a belief r.v. with distribution P n . It follows from θ n ⇒ θ c that E[H(θ n )] → H(θ c ).
For the mean-variance formulation, it suffices to show Var[H(θ n )] → 0. Noting that For the VaR formulation, notice that H(θ n ) ⇒ H(θ c ) by continuous mapping theorem. Take any ǫ > 0, since H(θ c ) − ǫ and H(θ c ) + ǫ are continuity points of the c.d.f. of δ H(θ c ) , for any α ∈ (0, 1) we have It is left to prove the case of CVaR. Take α ∈ (0, 1), by the translation invariance and monotonicity of CVaR, The proof is complete by noting |H(θ) − H(θ c )| is bounded and continuous on Θ.

Consistency of optimal solutions
Our next task is to show the consistency of optimal solutions. Let S n be the set of optimal solutions of a risk formulation when we have n i.i.d. input data samples, and let S be the set of optimal solutions to 2.1. If S n and S are not singletons, the convergence of solutions will be in the mode of set convergence, which can be defined in a variety of ways. In particular, we look at the following "distance" between S n and S. The deviation defined above is not a proper distance metric due to its lack of symmetry, but it is helps characterize the convergence of solutions. We are interested in finding out when D(S n , S) → 0 as n → ∞. In the rest of the paper, we will assume the solution space X is compact. This is not a strong assumption because very often the minimizers of a function are contained in a compact set. Let us make the following mild assumption.
By Doob's consistency theorem [33], assumption 3.2 implies that for any neighborhood V ∈ B Θ of θ c , P n (V ) → 1 a.s. (µ), which together with the compactness of Θ implies that ρ[H(x, θ n )] converges pointwise to H(x, θ c ) a.s. (µ) by an argument similar to the proof of theorem 3.1. Compared with assumption 3.1, which guarantees the consistency of posterior distribution for a specific θ c , assumption 3.2 is significantly weaker and only implies consistency for almost all θ c possibly except for a subset of Θ with measure 0 w.r.t. π (recall the meaning of µ-a.s. in section 2.2). The reason for making this weaker assumption will be specified later. Let us first look at the following theorem on the convergence of S n . Theorem 3.2. Suppose assumption 3.2 holds, Θ and X are compact, and H is continuous in θ for every x ∈ X . If there exists a measurable function κ : Θ → R + such that (i) |H(x, θ) − H(y, θ)| ≤ κ(θ) x − y , ∀x, y ∈ X ; (ii) Θ κ(θ)π(dθ) < ∞, then for the mean and the CVaR formulations, we have D(S n , S) → 0 a.s. (P θ c ) for almost all θ c ∈ Θ possibly except for a subset of measure 0 w.r.t. π. If in addition H is jointly continuous on X × Θ, then the same convergence holds for the mean-variance formulation. If we further have sup θ∈Θ κ(θ) < ∞, then the VaR formulation has the same convergence.
The proof of theorem 3.2 uses the following lemmas.  Lemma 3.5. Let X, Y be two r.v.s defined on some probability space (Ω, F , P) and ρ is a risk measure satisfying monotonicity and translation invariance.
Let us briefly explain why Lemma 3.3-3.5 are useful in proving theorem 3.2. It is not hard to see that the µ-a.s. pointwise convergence of ρ[H(x, θ n )] by Doob's consistency theorem is not strong enough to imply the convergence of S n . Instead, a sufficient condition given by lemma 3.3 is that ρ[H(x, θ n )] converges uniformly to H(x, θ c ). In the asymptotic analysis of SAA, such uniform convergence is guaranteed by uniform SLLN. In the case of risk formulations, however, we establish uniform convergence using lemma 3.4, which says that pointwise convergence together with a common Lipschitz constant implies uniform convergence. The final step is to find a common Lipschitz constant for each risk formulation. For the mean and CVaR formulations, this is done through a Lipschitz-type condition on H, where κ is an integrable function on Θ. This condition is useful because we can write posterior mean as a conditional expectation and apply martingale convergence theorem. It is noted here that the martingale converges µ-a.s., thus we introduce assumption 3.2 to ensure the µ-a.s. pointwise convergence of ρ[H(x, θ n )]. The proof for the cases of mean-variance and VaR relies on stronger assumptions on H: the continuity of H on X ×Θ is assumed to ensure the boundedness of H because joint continuity does not follow from separate continuity in each variable; the uniform Lipschitz assumption for the VaR formulation is made because VaR is not a coherent risk measure and lacks convexity. Under the assumption for the VaR formulation, all four risk formulations are trivially uniformly convergent. An immediate corollary of theorem 3.2 on the convergence of optimal values is as follows. Proof. Since theorem 3.2 guarantees the uniform convergence of the risk formulations, the conclusion follows from the second part of lemma 3.3.
As a special case, convex functions have the following nice property regarding uniform convergence: if a sequence of convex functions converges pointwise on an open set O ⊂ R n , then it also converges uniformly on any compact subset of O (theorem 3.1.4, [35]). This leads to the following corollary of theorem 3.2.
Corollary 3.2. Suppose assumption 3.1 holds, and H is continuous in θ for every x ∈ X . Let H : R d × Θ → R be convex in x for all θ ∈ Θ, then for the mean and the CVaR formulations, we have D(S n , S) → 0 a.s. (P θ c ).
Proof. By the uniform convergence property of convex functions, we only need to show that the two risk formulations are convex in x. The mean formulation is convex in x due to linearity. The convexity of the CVaR formulation follows from the properties of coherent risk measures, i.e., for x, y ∈ R n and λ ∈ [0, 1], In other words, CVaR preserves the convexity of the function it is taken w.r.t..

Asymptotic Normality of Risk Formulations
Another important aspect we would like to investigate is the asymptotic normality of risk formulations. Specifically, for a risk measure ρ among the four choices and a fixed x in X , we want to know if √ n{ρ[H(x, θ n )] − H(x, θ c )} converges in distribution to some r.v. depending on x. Moreover, we want to find out whether √ n{ρ[H(x, θ n )] − H(x, θ c )} converges weakly to a random function of x. In this section, we will show that if H is a well-behaved function, the asymptotic normality of risk formulations is inherited from the asymptotic normality of posterior distribution.
be the total variation metric between two probability measures µ and ν defined on (Ω, F ). We sometimes use P X to denote the distribution of a r.v. X. The following regularity assumption is needed for the asymptotic normality of posterior distribution.
(i) The experiment (P θ : θ ∈ Θ) is differentiable in quadratic mean at θ c with nonsingular Fisher information matrix I(θ c ).
(ii) The prior measure π is absolutely continuous in a neighborhood of θ c with a continuous density at θ c .
where N is the normal distribution and For more details on the conditions of lemma 4.1, we refer the reader to [36]. We note here that the (iii) of assumption 3.1 is a separability condition on θ c . Since we assume Θ is compact, the existence of uniformly consistent tests is guaranteed so long as P θ = P θ ′ for every pair θ = θ ′ and the maps θ → P θ are continuous for the total variation metric (lemma 10.6, [36]). The reader may rest assured that the conditions of lemma 4.1 are mild and therefore will not be of our concern.
We provide a little heuristics on how lemma 4.1 might be useful in our proof. The lemma says that the total variation distance between the distribution of the scaled residual √ n(θ n −θ c ) and N (∆ n,θ c , [I(θ c )] −1 ), a normal distribution with random mean ∆ n,θ c , goes to 0 in probability. It is important to notice that N (∆ n,θ c , [I(θ c )] −1 ) has a fixed variance [I(θ c )] −1 , and the mean ∆ n,θ c converges in distribution to N (0, [I(θ c )] −1 ) as n → ∞ by the central limit theorem (CLT). Suppose for illustration that Θ ⊂ R, then first we may try to argue that for some risk measure ρ, the difference between ρ[ √ n(θ n − θ c )] and (with a slight abuse In addition to lemma 4.1, the following lemma on the asymptotic normality of posterior mean (theorem 10.8, [36]) will be useful in characterizing the asymptotic distribution of the first-order term in Taylor expansion.

Asymptotic normality at a fixed x
The asymptotic distribution of risk formulations at a fixed x ∈ X is important in its own right. Such results allow us to construct CIs for the true performance of solutions.
Assumption 4.2 requires the (1+a)th moment of √ n(θ n −θ c ) to be bounded in probability (P θ c ). This can be viewed as an "in probability" version of uniform integrability (UI) because on the event {E [ √ n(θ n − θ c ) 1+a ] ≤ M ǫ }, we have for any K > 0, Thus, for a sufficiently large K, the tail expectation of √ n(θ n − θ c ) can be bounded by an arbitrarily small number with probability (w.p.) greater than 1 − ǫ w.r.t. P θ c for all n.
Another implication is the tightness of P √ n(θn−θ c ) since As we will see in the proof, assumption 4.2 plays an important role in bounding the remainder term in Taylor expansion. We do not attempt to give the exact conditions for assumption 4.2 to hold. Instead, we verify assumption 4.2 for some commonly used conjugate priors.
We know the first part of the RHS converges in distribution by lemma 4.2, so we only need to check if the second part is bounded in probability (P θ c ).
The proof of theorem 4.1 is a direct application of lemma 4.2, where lemma 4.2 is in fact built upon lemma 4.1. In the upcoming proof of VaR and CVaR's normality, we will follow the heuristics mentioned in the beginning of section 4 and work with lemma 4.1 explicitly. For the VaR formulation, we will prove the following theorem. Proof. See appendix A.3.
The following lemmas will be useful in proving theorem 4.2.
Lemma 4.3. Let X and Y be r.v.s taking values on the same probability space, then for any measurable function h, we have Lemma 4.4. Let X and Y be two r.v.s defined on the same probability space, where X and X + Y both have positive densities. Also let P(|Y | > δ) < ǫ for some δ > 0 and ǫ > 0. Then, for any α ∈ (0, 1) and ǫ ∈ (0, min{α, 1 − α}).
The proofs of lemmas 4.3, 4.4 and 4.5 can be found in appendix A.2. We provide some intuition on why theorem 4.2 holds and how lemma 4.3-4.5 can be used to prove theorem 4.2. The first thing to notice is that we can approximate √ n[H(θ n ) −H(θ c )] by its first-order Taylor approximation ∇H(θ c ) ⊺ [ √ n(θ n − θ c )]. So intuitively we would expect Then, by lemma 4.1 we know that √ n(θ n − θ c ) and N (∆ n,θ c , [I(θ c )] −1 ) are very close in total variation distance as n gets large. This motivates us to guess that Finally, a normal distribution N (µ, σ 2 ) is known to have VaR α [N (µ, σ 2 )] = µ + σΦ −1 (α). It follows that VaR α {∇H(θ c ) ⊺ N (∆ n,θ c , [I(θ c )] −1 )} = ∇H(θ c ) ⊺ ∆ n,θ c +σΦ −1 (α), which converges in distribution to N (σΦ −1 (α), σ 2 ). This is exactly the limiting distribution in theorem 4.2. Therefore, the result in theorem 4.2 makes sense, and the normality of the VaR formulation is inherited from the normality of posterior distribution. Nevertheless, a rigorous proof does not easily follow from the preceding intuition. There are two problems to take care of: (i) the gap between total variation distance and the difference of VaRs; (ii) the effect of the remainder term in Taylor expansion. The first problem is not hard. We will prove lemma 4.3, which guarantees that the total variation distance between ∇H(θ c ) ⊺ [ √ n(θ n − θ c )] and ∇H(θ c ) ⊺ N (∆ n,θ c , [I(θ c )] −1 ) converges to 0 in probability.
Although in general P Xn − P Yn TV → 0 does not imply VaR α (X n ) − VaR α (Y n ) → 0, we have a special structure in our favor, that is, {∇H(θ c ) ⊺ N (∆ n,θ c , [I(θ c )] −1 )} is a sequence of normal distributions with equal variance. This structure leads to lemma 4.5, which helps us close the gap between total variation distance and the difference of VaRs.
The main difficulty lies in the second problem, namely, how to deal with the remainder term in Taylor expansion. From the consistency of posterior distribution and assumption 4.2, we know that the remainder term converges to 0 in distribution a.s. (P θ c ). This motivates us to prove lemma 4.4, which basically says that if two r.v.s have a difference term that is very small with a large probability, then the difference between their VaRs, albeit at different risk levels, can be bounded. In other words, we can bound the difference between where ǫ is a number that can be made arbitrarily small. We already know from lemma 4.5 that we can bound the difference between where the latter is the VaR of a sequence of normal distributions with equal variance. From the structure of a normal distribution's VaR, we may choose ǫ small enough to make as small as we want for all n. To this point, everything is pieced together. The remainder term in the Taylor expansion of √ n[H(θ n ) − H(θ c )] is the main trouble in the proof of VaR's asymptotic normality. The reason is that VaR lacks structural properties, and we cannot move the remainder outside VaR. This problem, however, is not a concern for CVaR. Rather, the challenge in proving the case of CVaR is to bridge the gap between total variation distance and the difference of expectations. We will prove the following theorem on the asymptotic normality of the CVaR formulation.
where φ is the p.d.f. of N (0, 1) and Proof. See appendix A.3.
The proof of theorem 4.3 relies on the following lemmas.
Lemma 4.7. Suppose X n ∼ N (µ n , σ 2 ) and there is a constant C > 0 such that |µ n | < C for all n. Then for a fixed α ∈ (0, 1), for any ǫ > 0, there exists δ C > 0 such that The intuition behind theorem 4.3 follows the same pattern as theorem 4.2. It is not surprising that we will prove the normality of CVaR by bounding the difference Again, the idea is to take the Taylor expansion of √ n[H(θ n ) − H(θ c )] inside CVaR. But since CVaR is a conditional expectation, it is not difficult to bound The problem is how to bound the difference We know the distributions of ∇ θ H(θ c ) ⊺ [ √ n(θ n − θ c )] and ∇H(θ c ) ⊺ N (∆ n,θ c , [I(θ c )] −1 ) will be very close in total variation, but the convergence of their expected difference requires much more than that. In particular, we need ∇H(θ c ) ⊺ N (∆ n,θ c , [I(θ c )] −1 ) to be UI in probability (see assumption 4.2). Fortunately, we can further exploit the structure of normal r.v.s with equal variance. Note that the means of which are bounded in probability due to the weak convergence of ∆ n,θ c . We will prove lemma 4.6 and lemma 4.7 (see appendix A.2), which state that for a sequence of normal r.v.s with equal variance and uniformly bounded means, we can bound their tail expectations and their restricted expectations around VaR α uniformly. On the other hand, we have by assumption Putting these together, we are able to proceed from the convergence of total variation to the convergence of expectations.

Asymptotic normality of optimal values
Apart from the asymptotic normality of objective functions at a fixed x, it is also of practical interest to study the asymptotic distribution of the optimal values where ρ is a risk measure among the four choices. The reasons for studying the asymptotics of optimal values are twofold: (i) we can use the result to construct CIs for the true optimal objective value; (ii) its derivation requires establishing the weak convergence of risk formulations as random functions, which is an indispensable result for uncovering the hidden meaning of optimizing risk formulations. Let C(X ) denote the Banach space of all continuous functions on a compact set X equipped with the sup-norm. We will resort to a framework in [37], which says that if √ n(f n −f ) ⇒ Y x , where f n ,f and Y x are random elements of C(X ), then where S := arg min x∈Xf . Define To apply the framework (4.4), we need to carry out two main steps: (i) show that ρ{H(·, θ n )} and H(·, θ c ) are continuous functions on X ; (ii) show that g n (x) converges weakly to some random element of C(X ). Different from the study of SAA's asymptotics, where the functional CLT can be directly applied to step (ii), we have to start from scratch and decompose step (ii) into two parts. One is to show the convergence of g n 's finite dimensional distributions (f.d.d.s), and the other is to prove the tightness of g n 's distributions.
Step 1. For ρ{H(·, θ n )} and H(·, θ c ) to be continuous on X , it suffices to assume that H is continuous on X × Θ. The reason is as follows. Since X × Θ is compact, H is uniformly continuous on X × Θ by Heine-Cantor theorem, so for any ǫ > 0, there exists δ > 0 such that |H(x, θ) − H(x ′ , θ ′ )| < ǫ as long as (x, θ) − (x ′ , θ ′ ) < δ. Thus, for the mean, VaR and CVaR formulations, if (x, θ n ) − (x ′ , θ n ) = x − x ′ < δ, then by lemma 3.5, So ρ{H(x, θ n )} is uniformly continuous on X for mean, VaR and CVaR. The mean-variance formulation is also continuous because Step 2. The next step is to show the weak convergence of f.d.d.s. Since a finitedimensional random vector converges in probability to 0 if and only if each dimension converges in probability to 0, using the proofs of theorems 4.1, 4.2 and 4.3, it is straightforward to verify that each formulation's f.d.d. converges weakly. We shall defer the details of this step to appendix A.4.
Step 3. To show the tightness of g n 's distributions, we only need to check the following condition provided by theorem 7.8 in [38]. lim δ→0 lim sup n→∞ P θ c (ω(g n , δ) ≥ ǫ) = 0, ∀ǫ > 0, (4.5) where ω(f, δ) is the modulus of continuity of f ∈ C(X ) and is defined as The condition in (4.5) is also known as stochastic equicontinuity (s.e.). We now check the s.e. of all four risk formulations.
(i) Mean formulation: By Taylor expansion, so it suffices to show that the two terms of the RHS are both s.e.. Since E[ √ n(θ n − θ c )] converges in distribution, ∀η > 0, ∃M η > 0 such that Assuming that ∇ θ H(·, θ c ) is continuous on X (and hence uniformly continuous on X ), for any ǫ > 0, there exists δ η > 0 such that

It follows that on the event
So the first term is s.e.. For the second term, we only need to show that it suffices to show that sup x∈X |e(x, θ n )| ⇒ 0 a.s. (P θ c ). However, the fact that e is continuous on X × Θ implies that sup x∈X |e(x, ·)| is continuous on Θ. If we set e(x, θ c ) = 0, ∀x ∈ X as before to avoid technicalities, we would have and the rest follows from the proof of theorem 4.1.
(ii) Mean-variance formulation: Since ∇ θ H(x, θ c ) is bounded on X and e is bounded on X × Θ, it follows from the proof of theorem 4.1 that (iii) VaR formulation: Recall that in the proof of theorem 4.2 we try to bound where the bound of ( †) depends on x via e(x, θ n ), and the bound of ( † †) does not depend on x due to lemma 4.3. It follows from sup x∈X |e(x, θ n )| ⇒ 0 a.s. (P θ c ) that in probability (P θ c ). But we know is s.e. since ∇ θ H(·, θ c ) and σ 2 x are uniformly continuous on X and ∆ n,θ c converges in distribution. So the case of VaR is proved.
(iv) CVaR formulation: Since we have and the rest follows from the proof of the s.e. of the mean formulation.
The preceding derivation proves the following theorem. Further suppose X is a compact set, and H is continuous on X × Θ and partially differentiable at θ c for all x ∈ X , where ∇ θ H(·, θ c ) is continuous on X . Then, where S := arg min x∈X H(x, θ c ) and

Interpretation of Optimizing Risk Formulations
Although the intent of minimizing ρ[H(x, θ n )] is to improve the robustness of solutions against input uncertainty, it is not clear yet in what sense the solutions are more robust. Suppose ρ[H(x, θ n )] has a unique minimizer x * n , then x * n is a r.v. (assuming measurability) whose distribution depends on the distribution of input data ψ n . One may expect ρ[H(x * n , θ c )], i.e., the actual risk of x * n , to be controlled by solving risk formulations. However, our numerical results suggest otherwise. For example, the minimizer x * n of CVaR α [H(x, θ n )] does not necessarily have a smaller CVaR α [H(x * n , θ c )] than other formulations. We now interpret optimizing risk formulations based on the asymptotic normality results established in section 4.2. Take the VaR formulation as an example. From the derivation of theorem 4.4 we have the following weak convergence of random functions.
This can be rewritten as where o p (1/ √ n) stands for a term whose product with √ n converges to 0 in probability (P θ c ) uniformly in x. The LHS of (5.1) is the risk formulation we are minimizing, and the RHS can be viewed as the sum of the true objective H(x, θ c ) and some error terms. Compared with the mean formulation, which has we see that (5.1) has an extra deterministic bias term Φ −1 (α)σ x / √ n that vanishes as n → ∞. Combine (5.1) and (5.2) and we have So the VaR formulation is approximately equivalent to a weighted sum of posterior mean and a bias σ x / √ n, where the weight of the bias is Φ −1 (α). Even though the bias diminishes as n → ∞, it has an undeniable impact on the VaR formulation when n is small. In particular, if n is not too large (e.g. 20) and α is large (e.g. 99%), it is possible for the bias term to dominate E[H(x, θ c )] and we are close to solving min x∈X σ x / √ n. Similarly, the CVaR formulation can be rewritten as, For the mean-variance formulation, by making appropriate assumptions of UI it can be shown that the variance satisfies Therefore, the mean-variance, VaR and CVaR formulations are approximately equivalent to a weighted sum of the mean formulation and σ x / √ n (or (σ x / √ n) 2 ), where the weight of σ x / √ n is controlled by α in the VaR and CVaR formulations, and the weight of (σ x / √ n) 2 is controlled by the constant c in the mean-variance formulation.
But what does it mean to minimize σ x / √ n? For a given solution x, if we construct a CI for its true performance H(x, θ c ) using theorem 4.1 -4.3, then σ x / √ n is exactly proportional to the half-width of the CI. Thus, a narrower CI implies higher accuracy of estimating the true performance H(x, θ c ), while a wider CI indicates more risk because we are less confident about how a solution actually performs. In other words, the risk formulations are essentially seeking a tradeoff between average performance and the robustness of actual performance. It is also interesting to notice that σ x depends on ∇ θ H(x, θ c ) and I(θ c ), where ∇ θ H(x, θ c ) is the sensitivity of the function H to the perturbation of the parameter θ c (i.e. the input uncertainty), while I(θ c ) captures how much information a data sample carries about the true parameter.

Numerical Examples
We use three examples to verify the theoretical results in section 3, 4, and 5. All three examples assume that the stochastic uncertainty ξ follows exponential distribution with rate θ c . This is no loss of generality because the asymptotics of posterior distribution holds for most choices of prior in practice. We choose a (conjugate) Gamma prior on ξ so that the posterior distribution is easy to compute. The following aspects are examined: (i) the convergence of x * n to the set of true optimal solutions as n increases; (ii) the coverage probabilities of CIs constructed from theorems 4.1, 4.2 and 4.3; (iii) how x * n changes as α and c change; (iv) the actual risk of solutions, i.e., ρ[H(x * n , θ c )]. Let z 1−β denote the (1 − β)quantile of a standard normal distribution. For a fixed x ∈ X , the asymptotically valid 100(1 − β)% CIs of H(x, θ c ) are listed as follows (in the form of center ± half-width).
(i) The mean and mean-variance formulation (ii) The VaR formulation (iii) The CVaR formulation Since the risk formulations do not have simple closed-forms, we use SAA with 10,000 samples to optimize the risk formulations, namely, we draw 10,000 θ's independently from the Gamma posterior distribution, and optimize a sample average version of the risk formulations. Although the solution space X and the parameter space Θ in these examples are not compact, our numerical results suggest that the consistency and asymptotic normality of risk formulations remains valid.

Quadratic objective
In the first example, h is a quadratic function (x − ξ) 2 . We use this simple example because it has a clear structure for illustrative purpose. Since we assume ξ ∼ Expo(θ c ), we have and σ x can be computed as Set θ c = 1 so that the true objective function H(x, θ c ) has a unique minimizer 1 on X = R, and σ x has a unique minimizer 2. We shall verify our theoretical results on this example via simulation experiment. To begin with, we would like to see if x * n converges to 1 as n increases. For a given n, we draw 1,000 sets of i.i.d. input data from Expo(1), where each dataset is of size n. Then, for each dataset, we optimize the risk formulations using SAA. The histograms of x * n for n = 20, 50, 100 and 200 are given in Figure 2. It is clear that the optimal solutions of all four risk formulations converge to 1 as n increases.
Next, we verify the validity of the CIs constructed from theorem 4.1, 4.2 and 4.3. For x = 2, α = 0.95, and c = 1, we use 1,000 i.i.d. input datasets to estimate the probabilities with which the 95% CIs cover the true objective value H(2, 1). When we compute the CIs, the unknown θ c is replaced with the Bayes estimator E[θ n ]. As is shown in figure 3, all of the CIs attain the target value when n gets large, even though the CIs of the mean-variance, VaR and CVaR formulations have relatively low coverage probabiltiies when n is small.
We now examine how x * n respond to the changes of α and c when n is fixed at 20. As before, we draw 1,000 i.i.d. datasets to estimate the distribution of x * n , and the risk formulations are solved using SAA. The histograms of x * n for the cases of mean-variance, VaR and CVaR are given in figure 4. Observe that as α and c increase, the average value of x * n shifts from 1 towards 2, i.e., the minimizer of σ x . This is consistent with our insights that the mean-variance, VaR and CVaR formulations put more weight on σ x / √ n when α and c are larger. Consequently, the solutions of the mean-variance, VaR and CVaR formulations have narrower CIs of true performance compared with the mean formulation as indicated by table 6.1, demonstrating the robustness gained from optimizing risk formulations. In particular, the average x * n of the CVaR formulation at α = 0.99 is 1.9717, which almost coincides with 2. We remark here that σ x = 0 when x = 2, but it does not make sense for a CI to have 0 width. A better way to interpret this is that the first-order term in the Taylor expansion disappears, so the convergence in theorem 4.1, 4.2 and 4.3 will be at a faster rate, e.g., o p (1/n), leading to a narrower CI and hence lower risk.  formulations. This is because the solutions of CVaR are closer to 2 and thus have larger objective values. In fact, the mean-variance, VaR and CVaR all have worse ρ[H(x * n , θ c )] than the mean formulation for every choice of ρ. One can view this as the price of robustness.

Multimodal objective
In the second example, we want to test the convergence of optimal solutions when the true objective function has multiple minimizers. We directly specify H instead of h as follows.
where ǫ is a positive number we can tune to modify the shape of H. The graph of H(x, θ c ) under different choices of ǫ is plotted in figure 5. One can see that for any ǫ < 2, H(x, θ c ) is a multimodal function with two global minimizers on X = R. If ǫ decreases to 0, then the value of H(x, θ c ) at x = 1 increases to ∞ and eventually becomes a point of discontinuity. Suppose θ c = 1, then H(x, θ c ) has two minimizers 0 and 2. Furthermore, we have which has a unique minimizer 1 on R. Similar to section 6.1, we examine the convergence of optimal solutions, the coverage probabilities of CIs, the response of x * n to the change of c and α, and the average width of true performances' CIs. We can see from figure 6 that as n increases, the optimal solutions of risk formulations tend to concentrate around the true optimal solutions 0 and 2. Also, figure 7 shows the effectiveness of CIs of H(2, 1) as n gets large. The histograms in figure 8 indicates the trend of x * n moving towards 1 with the increase of c and α, since we are closer to minimizing σ x / √ n. When n = 20, the average width of the mean formulation's CIs is 4.5796, which is much larger than the other three formulations as shown in table 6.3. Thus, when the true objective function has multiple minimizers, the risk formulations still produce more robust solutions.
We make an extra remark on the price of robustness through this example. If a decision maker is very conservative and chooses a large c for the mean-variance formulation, then s/he tends to get solutions that are close to 1. Thus, the true performance of x = 1 determines how bad the solutions can be. In theory, we can make ǫ arbitrarily close to 0 to make H(1, θ c ) as large as possible. As a result, the solutions of the mean-variance formulation may have much worse performances than the true optimal solutions. Therefore, the choice of c heavily depends on the structure of the true objective function. However, this weight choosing problem is less of a concern for the VaR and CVaR formulations, because the normal distribution is light-tailed. Even for α = 0.99, the weights of the VaR and CVaR formulations are only 2.3263 and 2.6652, respectively.

Newsvendor problem
The last example is a newsvendor problem with lost sales. The setting is as follows. A decision maker needs to decide the amount of newspaper to order in an effort to maximize the expected profit. Suppose the market demand ξ ∼ Expo(θ c ), where θ c is unknown and must be estimated from finite historical data. If the inventory exceeds the demand, then there are newspapers left unsold. If the demand is higher than the inventory, then the unfulfilled demand is simply lost. Let x be the order quantity, p be the unit price and c be the unit cost. To apply risk formulations, we consider minimizing the expected loss, which is equivalent to maximizing the expected profit. The expected loss is given by which has a unique minimizer 1 θ c ln p c on X = R + . We also have which has a unique minimizer 0 on X = R + . This is because there is no uncertainty at x = 0: the expected loss is always 0 if we do not order anything, regardless of the demand distribution. Assuming p = 2, c = 1 and θ c = 1/20, we verify our interpretation of optimizing risk formulations in the same way as in the first two examples. The response of x * n to the change of c and α is shown in figure 9. We can still see a general trend of x * n shifting towards 0, but the magnitude of shift is much less significant compared with previous examples. The reason is that σ x / √ n is dominated by E[H(x, θ n )], and we have to make the weight very large to see x * n approaching 0. As we mentioned in section 6.2, the VaR and CVaR formulations have a limited ability of increasing the weight of σ x / √ n due to the light-tailedness of normal distribution. This turns out to be desirable in this example, because a solution close to 0 is too conservative: the main purpose of selling newspapers is to make profit rather than control loss. Therefore, if the goal is to achieve a balance between expected performance and the risk of actual performance, then more effort should be devoted to exploring the

Concluding Remarks
This paper addresses some fundamental questions on the consistency and asymptotic normality of risk formulations. The most important insight is that risk formulations explicitly seek a tradeoff between posterior mean performance and the half-width of the actual performance's confidence interval. A question of practical interest is how to choose the weight c in the mean-variance formulation and the risk level α in the VaR and CVaR formulations, since these two parameters control the balance between average performance and robustness. In addition, our proof assumes the compactness of the parameter space Θ, but it is worth studying general cases since many priors are not supported on compact sets. It is also interesting to consider a nonparametric setting with Dirichlet prior, although the asymptotics may be much more complicated. We hope our work will inspire more study on the application of risk measures in simulation optimization under input uncertainty.

A.2 Proof of lemmas
Proof of lemma 3.5. Note that X ≤ Y + X − Y ∞ , so ρ(X) ≤ ρ (Y + X − Y ∞ ) by monotonicity, and ρ(X) − ρ(Y ) ≤ X − Y ∞ by translation invariance. Now reverse X and Y to get the first result. For the second result, by subadditivity ρ(Y ) + ρ(X − Y ) ≥ ρ(X), and by monotonicity ρ( . Reverse X and Y and the proof is complete.
where the last term is P X − P Y TV . The result follows from taking supremum over B ∈F on both sides.
It suffices to prove the following two inequalities.
So the proof is complete.
Proof of lemma 4.6. Let Z ∼ N (C, σ 2 ), then ∃M > 0 such that E[Z½ {Z>M } ] < ǫ/2, and we claim that this M responds to the ǫ challenge in the lemma. Indeed, since X n − µ n + C ∼ N (C, σ 2 ), we have for all n, By symmetry, we also have −E[X n ½ {Xn<−M } ] < ǫ/2 for all n.
The rest follows from the proof of the mean formulation. ≤ 2ME[κ(θ) | F n ] x − y , thus the first term has a Lipschitz constant independent of θ n . The second term can be bounded in a similar way. Finally, for the VaR formulation we assume sup θ∈Θ κ(θ) < ∞, which means there exists L > 0 such that H(x, θ) − H(y, θ) ≤ L x − y for all θ ∈ Θ and x, y ∈ X . By the second part of lemma 3.5, Expand H around θ c to obtain VaR α ∇H(θ c ) ⊺ [ √ n(θ n − θ c )] + e(θ n ) √ n(θ n − θ c ) .
Proof of theorem 4.3. To simplify notations, let X n := ∇H(θ c ) ⊺ [ √ n(θ n − θ c )] and let Y n be a r.v. with distribution ∇H(θ c ) ⊺ N (∆ n,θ c , [I(θ c )] −1 ). Also let v n x := VaR α (X n ) and v n y := VaR α (Y n ). Notice that which converges to 0 in distribution from the proof of theorem 4.1. So it suffices to show that CVaR α (X n ) − CVaR α (Y n ) → 0 in probability (P θ c ). First of all, since E[Y n ] = ∆ n,θ c converges in distribution, it is bounded in probability (P θ c ). Thus for any ǫ > 0, there exists In addition, by lemma 4.6, on the event E 1,n we can find M 2 > 0 such that Since the proof of theorem 4.2 implies that |v n x − v n y | → 0 in probability (P θ c ), we can find N 1 > 0 such that the event E 2,n := {|v n x − v n y | < δ M 1 } occurs w.p. (P θ c ) greater than 1 − ǫ/4 for all n > N 1 . Furthermore, since E[ √ n(θ n − θ c ) ] is bounded in probability (P θ c ), there exists M ǫ > 0 such that P θ c (E[ √ n(θ n − θ c ) ] > M ǫ ) < ǫ/4, ∀n.