# Compute a probability using stochastic simulation

My problem is to solve the following optimisation problem using GA (Genetic Algorithm)and stochastic simulation. The goal is to solve the maximisation problem : \begin{equation*} \begin{aligned} \text{Max}_{x_{1},\cdots,x_{n}} \qquad &\mathbb{E}[\xi_{1}x_{1}+\xi_{2}x_{2}+\cdots+\xi_{n}x_{n}]\\ \text{s.t.} \qquad &\left\{ \begin{array}{rcr} & \mathbb{P}r\{(0.4-(\xi_{1}x_{1}+\xi_{2}x_{2}+\cdots+\xi_{n}x_{n}))\geqslant r\}\leqslant \alpha(r)\\ & \forall r\geqslant 0\\ &x_{1}+x_{2}+\cdots+x_{n}=1\\ &x_{1} \qquad \forall i \in \{1,2,\cdots,n\}.\\ \end{array} \right. \end{aligned} \end{equation*} where $n=5$, $\ \xi_{i}\sim N(\mu_{i},\sigma_{i})$ with $\mu_{i}\in\{0.1, 0.2, 0.3,0.4,0.5\}$, $\sigma_{i}=1$ and $\ \alpha(r)=\dfrac{1}{(r+1.1)^{4}}$, $\ r\geqslant 0$.

• Is the "$r$" after the "$\mathbb{P}$" a reference to $r$ or did you intend to write $\Pr$ ("probability of") instead? And presumably $b$ and the $x_i$ are constants, right? Finally, your question appears to have nothing to do with stochastic simulation but that is mentioned in the title: have you left out some important details? – whuber Dec 6 '13 at 16:39
• The main goal is to solve the following portfolio problem : \begin{equation*} \begin{aligned} \text{Max}_{x_{1},\cdots,x_{n}} \qquad &\mathbb{E}[\xi_{1}x_{1}+\xi_{2}x_{2}+\cdots+\xi_{n}x_{n}]\\ \text{s.t.} \qquad &\left\{ \begin{array}{rcr} & \mathbb{P}r\{(b-(\xi_{1}x_{1}+\xi_{2}x_{2}+\cdots+\xi_{n}x_{n}))\geqslant r\}\leqslant \alpha(r)\\ & \forall r\geqslant 0\\ &x_{1}+x_{2}+\cdots+x_{n}=1\\ &x_{1} \qquad \forall i \in \{1,2,\cdots,n\}.\\ \end{array} \right. \end{aligned} \end{equation*} – Lea Dec 6 '13 at 18:56
• Because $\mathbb{E}[\xi_i]=0$ for all $i$, the $x_i$ are constants, and expectation is linear, $\mathbb{E}[\xi_{1}x_{1}+\xi_{2}x_{2}+\cdots+\xi_{n}x_{n}]$ = $x_1\mathbb{E}[\xi_1]+\cdots+x_n\mathbb{E}[\xi_n]$ = $x_1\times 0+\cdots+x_n\times 0=0$ no matter what. That's not much of a portfolio :-). – whuber Dec 6 '13 at 19:32

If the $\xi_i \sim N(\mu_i, 1)$ are independent, then $\sum \xi_i x_i \sim N(\sum \mu_i x_i, \sum x_i^2)$. Thus your probability question can be computed as

\begin{align*} \Pr\left(b - \sum \xi_i x_i \ge r\right) &= \Pr\left(\sum \xi_i x_i \le b - r\right) \\ &= \Pr\left(\frac{\sum \xi_i x_i - \sum \mu_i x_i}{\sqrt{\sum x_i^2}} \le \frac{b - r - \sum \mu_i x_i}{\sqrt{\sum x_i^2}}\right) \\ &= \Phi\left( \frac{b - r - \sum \mu_i x_i}{\sqrt{\sum x_i^2}} \right) \end{align*}

where $\Phi$ is the cdf of the standard normal.

Denoting $\mu$ as the vector of means and $x$ the vector of $x_i$s, your objective is $\mathbb{E}[\sum x_i \xi_i] = \mu^T x$.

The first constraint is then that $x$ must satisfy $\Phi\left( \frac{b - r - \mu^T x}{\lVert x \rVert} \right) \le \alpha(r) \quad \forall r > 0$. The others are simple: $1^T x = 1$, and presumably your last constraint is meant to be $x_i \ge 0 \quad \forall i$.

Note that this constraint is a function of only $\mu^T x$ and $\lVert x \rVert$. For a fixed $\lVert x \rVert$, and assuming an $\alpha(r)$ generally like the one you specified on the other question of $(1.1 - r)^{-4}$, I'm pretty sure there is a certain value of $\mu^T x$ for which all lower values do not satisfy the constraint and all higher values do. Cauchy-Schwarz gives us $\lvert \mu^T x \rvert \le \lVert \mu \rVert \lVert x \rVert$, so for a given $\lVert x \rVert$ we can binary-search for the minimum value of $\mu^T x$ such that the constraint is satisfied.

I don't have time to work this out right now, but I think we can use this to figure out the shape of the constraint set and/or get a projection onto it to solve your optimization problem.

• This is correct, but what if $n=5$ and $\xi_{i}$ are $N(0.1,1)$ ,$N(0.2,1)$,$N(0.4,1)$, $N(0.6,1)$,$N(0.8,1)$. – Lea Dec 6 '13 at 21:12
• That changes the question profoundly, Lea. Please modify your question to show that's what you really intended to ask. – whuber Dec 6 '13 at 21:33
• Dougal, this is exactly the way to do it, but is it possible to have some package or some code in R to do it – Lea Dec 7 '13 at 18:54