# Primordial Black Holes Are True Vacuum Nurseries

Louis Hamaide<EMAIL_ADDRESS>Theoretical Particle Physics and Cosmology,
King’s College London,
Strand, London WC2R 2LS, United Kingdom  AMOPP, University College London,
Gower Place, London WC1E 6BS, United Kingdom  Lucien Heurtier
<EMAIL_ADDRESS>Theoretical Particle Physics and Cosmology,
King’s College London,
Strand, London WC2R 2LS, United Kingdom  IPPP, Durham University, Durham, DH1
3LE, United Kingdom  Shi-Qian Hu<EMAIL_ADDRESS>Theoretical Particle
Physics and Cosmology, King’s College London,
Strand, London WC2R 2LS, United Kingdom  Andrew Cheek<EMAIL_ADDRESS>Astrocent, Nicolaus Copernicus Astronomical Center of the Polish Academy of
Sciences, ul.Rektorska 4, 00-614 Warsaw, Poland

###### Abstract

The Hawking evaporation of primordial black holes (PBH) reheats the Universe
locally, forming hot spots that survive throughout their lifetime. We propose
to use the temperature profile of such hot spots to calculate the decay rate
of metastable vacua in cosmology, avoiding inconsistencies inherent to the
Hartle-Hawking or Unruh vacuum. We apply our formalism to the case of the
electroweak vacuum stability and find that a PBH energy fraction
$\beta>7\times 10^{-80}(M/\mathrm{g})^{3/2}$ is ruled out for black holes with
masses $0.8\,{\rm g}<M<10^{15}\,{\rm g}$.

††preprint: IPPP/23/65

## Introduction. —

Many models of early-universe cosmology predict the formation of light
primordial black holes that have evaporated long before the onset of Big Bang
Nucleosynthesis (BBN). Although such evanescent PBHs constitute unique relics
of primordial cosmic history, confirming their existence with observations can
be challenging, if not impossible, due to their small density fraction in
cosmology.

Hawking’s famous result states that light black holes can be extremely hot,
$T_{H}=(8\pi GM)^{-1}$, which would have substantial impact on the local
environment. By using the so-called Hartle-Hawking vacuum — in which PBHs are
assumed to be in thermal equilibrium with their surrounding plasma Hartle and
Hawking (1976) — the authors of Refs. Gregory  _et al._ (2014); Burda  _et
al._ (2015a, b, 2016) proposed that PBHs could trigger first-order phase
transitions in the early Universe and endanger the meta-stability of the
electroweak vacuum. Indeed, unlike sphaleron processes which pump energy
directly from a plasma with infinite heat capacity to form a true vacuum
bubble, PBHs have the ability to release part of their mass energy to support
the formation of true-vacuum bubbles around them. If true, such a claim would
severely constrain the formation of PBHs with masses smaller than
$\mathcal{O}(10)\mathrm{g}$ in the early Universe or impress the importance of
stabilising the Higg’s vacuum via new physics Hiller  _et al._ (2022).

Since then, this claim has been subject to controversy: PBHs are rarely
surrounded by a plasma of temperature $T=T_{H}$ and several authors pointed
out that the Unruh vacuum Unruh (1976) — in which PBHs radiate energy in an
empty Universe — should be used instead, largely mitigating the aforementioned
result Kohri and Matsui (2018); Hayashi  _et al._ (2020); Shkerin and
Sibiryakov (2021, 2022); Strumia (2022); Briaud  _et al._ (2022). Moreover,
thermal quantum corrections to the Higgs potential were argued to be relevant
for PBHs lighter than $\mathcal{O}(10^{3})\mathrm{g}$, increasing the energy
of true vacuum bubbles, and effectively rescuing the stability of the
electroweak vacuum Strumia (2022).

In this letter, we point out that evaporating PBHs neither live in the vacuum
nor are in thermal equilibrium with their surrounding. Instead, by depositing
energy locally into the thermal plasma, Hawking radiation supports the
formation of hot spots around PBHs Das and Hook (2021); He  _et al._ (2023).
Building on the results derived in Gregory  _et al._ (2014); Burda  _et al._
(2015a, b, 2016); Dai  _et al._ (2020), we provide a new avenue to calculate
false-vacuum decay (FVD) rates around PBHs in cosmology. In particular, we
will show that the presence of such PBH hot spots can indeed seed the
formation of true electroweak vacuum bubbles in the early Universe,
circumventing the criticisms listed above.

## Bubble Action at $T<T_{H}$. —

Let us for the time being consider an homogeneous scalar field configuration
$\phi$ living in a metastable minimum of its potential $V(\phi)$. The decay
rate of this configuration is determined in curved spacetime by a saddle point
“bounce” solution of the Euclidean action. Around a Schwarzschild black hole
living in a background plasma with an arbitrary constant temperature $T$, one
can write the Euclidean action of a time-independent scalar field bubble
configuration as in Burda  _et al._ (2016); Dai  _et al._ (2020)

$I_{\rm b}[T]=\beta\int dx^{3}\,\sqrt{-h}\left(-\frac{R}{16\pi
G}+\frac{1}{2}h^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi+V(\phi)\right)\,.$
(1)

in which $h_{\mu\nu}$ is the spatial part of the metric — assumed to be time-
independent as well — and $\beta=1/T$ denotes the Euclidean periodicity. At
first sight, because the periodicity $\beta$ appears in this expression as an
overall factor, it is tempting to conclude that a low plasma temperature would
lead to a large bubble action, hence suppressing the false-vacuum decay rate
exponentially:

$T\ll T_{H}\ \Rightarrow\ \Gamma_{\rm FVD}\propto\exp\left(-I_{\rm
b}[T]\right)\ll\exp\left(-I_{\rm b}[T_{H}]\right)\,.$ (2)

However, it is important to note that there are two pieces missing in Eq. (1)
Gregory  _et al._ (2014). First, the gravitational background action should be
subtracted to extract the actual energy of the bubble. Second, if the
temperature entering the Euclidean action’s periodicity differs from the
Hawking temperature at the black hole horizon, the system features a conical
deficit at the horizon that needs to be accounted for in the calculation. Ref.
Gregory  _et al._ (2014) calculated both contributions in the thin-wall
approximation. In particular, the contribution from the conical deficit was
evaluated using a mathematical regularisation procedure: using a temperature
profile interpolating over a scale $\epsilon$ between a constant temperature
$T(r)=T$ for $r-r_{H}\gg\epsilon$, and $T(r_{H})=T_{H}$ at the horizon, one
can show that the Ricci contribution to the action is of order

$I_{\rm
b}[T]\supset-\frac{\beta_{H}-\beta}{\beta}\frac{\mathcal{A}}{4G}+\mathcal{O}(\epsilon)\,,$
(3)

in the $\epsilon\to 0$ limit. In this equation $\mathcal{A}$ stands for the
area of the black hole, and $\beta_{H}=1/T_{H}$ the corresponding inverse
Hawking temperature.

Remarkably, after summing up all the contributions, the authors of Ref.
Gregory  _et al._ (2014) proved, in the $\epsilon\to 0$ limit, that the
$\beta$-dependency of the total Euclidean action exactly cancels, leading to
the seminal result

$I_{\rm b}[T]=\frac{\mathcal{A}_{+}}{4G}-\frac{\mathcal{A}_{-}}{4G}=I_{\rm
b}[T_{H}]\,,$ (4)

where $\mathcal{A}_{+}$ ($\mathcal{A}_{-}$) denotes the area of the black hole
horizon before (after) the formation of the bubble. This property was then
generalised beyond the thin-wall approximation, and in the presence of matter
in Ref. Burda  _et al._ (2015b). Despite its remarkable generality, to our
knowledge, this result was only used in the Hartle-Hawking vacuum where
$T=T_{H}$ to derive constraints on PBHs using the electroweak vacuum
stability. In what follows, we will show that in a realistic cosmological set-
up, PBHs do live in a plasma with a temperature lower than their Hawking
temperature, but also that the mathematical smoothing used to derive Eq. (3)
does actually correspond to a physical situation that we will describe.

## Thermal Profile Around PBHs. —

Throughout cosmic history, PBHs are surrounded by two main sources of energy:
$(i)$ their Hawking radiation, made of relativistic particles with average
energies $\langle E\rangle\sim T_{H}$, and $(ii)$ the energy of the ambient
plasma surrounding the black hole, with particles with average energy $\langle
E\rangle\sim T$. Studies considering the Hartle-Hawking vacuum assume that
$T=T_{H}$ and both energy sources are in perfect equilibrium. Instead, studies
working with the Unruh vacuum neglect the presence of any ambient plasma and
only consider the Hawking radiation, effectively working in the limit where
$T\to 0$.

As a matter of fact, both of these pictures are incomplete. When PBHs form,
the Universe’s temperature may be much larger than their Hawking temperature.
However, the end of PBH evaporation typically takes place much later, when the
Universe’s temperature is much colder than its original Hawking temperature.
Nevertheless, a big piece is missing in this discussion: In the intermediate
period, the Standard Model (SM) radiation that forms Hawking radiation
unavoidably interacts with the surrounding plasma. This energy deposition is
localised in space, and instead of reheating the Universe in a homogeneous
way, PBHs are expected to heat up their surrounding environments. This energy
deposition leads to an inhomogeneous temperature profile forming around the
PBHs that persists long after they have completely evaporated. As a
consequence, the plasma temperature around PBHs can differ by orders of
magnitude from the average temperature in the Universe.

The temperature profile evolution around PBHs was explored in Ref. He  _et
al._ (2023); Das and Hook (2021). Assuming a universal coupling constant
$\alpha$ to encode particle interactions between Hawking radiation and the
plasma, the authors of Ref. He  _et al._ (2023) showed that for PBHs with
masses

$\displaystyle M$ $\displaystyle\gtrsim$ $\displaystyle M_{\star}\,,$
$\displaystyle M_{\star}$ $\displaystyle\equiv$ $\displaystyle
0.8\mathrm{g}\left(\frac{\alpha}{0.1}\right)^{-\frac{11}{3}}\,,$ (5)

the formation of a hot spot is quicker than the evaporation of the black hole.
PBHs are thus rapidly surrounded by an initial hot spot featuring a constant
temperature

$\displaystyle T_{\rm plateau}$ $\displaystyle\approx$ $\displaystyle 2\times
10^{-4}\left(\frac{\alpha}{0.1}\right)^{\frac{8}{3}}T_{H}$ (6)
$\displaystyle\times$ $\displaystyle\left(\frac{g_{\star}(T_{\rm
plateau})}{106.75}\right)^{-\frac{2}{3}}\left(\frac{g_{\star}(T_{H})}{106.75}\right)^{\frac{2}{3}}\,,$

over a distance

$\displaystyle r_{\rm plateau}$ $\displaystyle\approx$ $\displaystyle 7\times
10^{8}\left(\frac{\alpha}{0.1}\right)^{-6}r_{H}$ (7) $\displaystyle\times$
$\displaystyle\left(\frac{g_{\star}(T_{\rm
plateau})}{106.75}\right)\left(\frac{g_{\star}(T_{H})}{106.75}\right)^{-1}\,,$

where $g_{\star}(T)$ denotes the number of degrees of freedom that are
relativistic at temperature $T$. At a time smaller than the evaporation time
$t_{\rm ev}\equiv\Gamma_{\rm ev}^{-1}$, where $\Gamma_{\rm
ev}\equiv|\dot{M}/M|\propto M^{-3}$ can be calculated from e.g. Refs.
MacGibbon and Webber (1990); MacGibbon (1991), the value of this plateau
temperature is constant. On larger distances, diffusion dominates over energy
deposition and the temperature decreases, smoothly interpolating between the
plateau temperature and the Universe’s temperature.

Later on, the PBH mass starts decreasing, its Hawking temperature increases,
and so does the plateau temperature, while its radius decreases. After the
black hole mass reaches $M=M_{\star}$, the plateau temperature saturates at a
maximum value He  _et al._ (2023)

$\displaystyle T_{\rm max}$ $\displaystyle\approx$ $\displaystyle 2\times
10^{9}\mathrm{GeV}\left(\frac{\alpha}{0.1}\right)^{\frac{19}{3}}$ (8)
$\displaystyle\times$ $\displaystyle\left(\frac{g_{\star}(T_{\rm
max})}{106.75}\right)^{-\frac{4}{3}}\left(\frac{g_{\star}(T_{H})}{106.75}\right)^{\frac{5}{6}}\,,$

over a radius

$r_{\rm max}=\left.r_{\rm
plateau}\text{{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\Large
j}}\right|_{T_{H}=T_{\rm max}}\,.$ (9)

At this point, the PBH, in the last stage of its evaporation, is not able to
provide enough energy to reheat the hot spot further.

To obtain these results, the authors of Ref. He  _et al._ (2023) assumed that
Hawking radiation particles all deposited their energy at an averaged time
$t_{\rm dep.}\equiv\Gamma(T)^{-1}$, where the first-splitting rate was
estimated to be

$\Gamma(T)\sim\alpha^{2}T\sqrt{\frac{T}{T_{H}}}\,.$ (10)

We argue that this approximation cannot remain valid near the horizon. Indeed,
even if the temperature within the hot spot was perfectly homogeneous, the
energy deposition probability of Hawking radiation

$dP\sim\Gamma(T)e^{-(r-r_{H})\Gamma(T)}dr\,,$ (11)

is a decreasing function of the radius. In addition, the Hawking radiation
energy density also decreases like $\propto r^{-2}$. Energy deposition is thus
more efficient close to the horizon, leading to a local increase of the
temperature. Because an increase of the local temperature also corresponds to
an enhancement of the first-splitting rate of Eq. (10), the plasma temperature
must interpolate between the Hawking temperature at the horizon and the
temperature of the hot spot calculated on larger distances in Eq. (6) and (8).
In Fig. 1, we present the qualitative temperature profile that must be present
around the hot spot at $t\ll\Gamma_{\rm ev}^{-1}$, when the hot spot core
temperature was $T_{\rm plateau}$ (upper panel), and at $t\lesssim\Gamma_{\rm
ev}^{-1}$, when the plateau temperature has saturated at $T_{\rm max}$ (lower
panel).

Figure 1: Sketch of the temperature profile around PBHs throughout the
evaporation process.

Note that by enforcing the smooth behaviour of the plasma $T\rightarrow T_{H}$
as $t\rightarrow r_{H}$, we have introduced a physical realisation of the
regularisation procedure of Ref. Gregory  _et al._ (2014) that was necessary
to evaluate the contribution of the conical deficit to the action in Eq. (3).
In full generality, one would need to solve the heat equation in the vicinity
of PBHs to obtain the exact temperature profile and, in particular, to
evaluate the smoothing size $\varepsilon$ that is sent to zero to obtain Eq.
4. For simplicity, we will assume here that the typical size of this smoothing
is negligible as compared to the typical size of a true vacuum bubble and
leave a precise study of the profile geometry contribution for future work.

## Bounce Solution and Nucleation Rate. —

Let us now consider the case of a homogeneous scalar-field configuration
$\phi_{+}$ living in a metastable minimum of its potential where
$V(\phi_{+})=0$ before tunneling to the true vacuum, located at $\phi_{-}$,
for which $V(\phi_{-})<0$. To evaluate the nucleation rate of true-vacuum
bubbles, one needs to search for a saddle-point ‘bounce’ solution of the
action in Eq. (1). In a finite temperature background, one also needs to
account for thermal quantum corrections to the scalar potential, such that
$V(\phi)=V(\phi,T)$. Following Ref. Dai  _et al._ (2020), we parametrise the
metric as

$ds^{2}=f(r)e^{2\delta(r)}d\tau^{2}+\frac{dr^{2}}{f(r)}+r^{2}d\Omega_{2}^{2}\,,$
(12)

where $f(r)\equiv 1-2G\mu(r)/r$ with $\mu(r)$ representing the local black
hole mass. We then search for solutions of the classical equations of motion
satisfying the boundary conditions111See Ref. Dai  _et al._ (2020) for more
details on the resolution.

$\lim_{r\rightarrow\infty}\phi(r)\rightarrow\phi_{+}\,\,,\,\,M_{+}\equiv\lim_{r\rightarrow\infty}\mu(r)\,.$
(13)

Assuming that the energy of the bubble is entirely provided by the mass
variation of the black hole inside the bubble, its mass at infinity is
unchanged, and $M_{+}$ corresponds to the PBH mass before the transition. In
the near horizon limit, due to the presence of a negative cosmological
constant $\Lambda_{-}\equiv V(\phi_{-})$, the mass parameter at the horizon
$\mu_{-}\equiv\lim_{r\rightarrow r_{H}}\mu(r)$ slightly differs from the
physical (ADM) mass of the black hole inside the bubble, $M_{-}$. The
Schwarzschild radius can be related to these masses by the relation

$r_{H}=2G\mu_{-}=2GM_{-}+\frac{\Lambda_{-}r_{H}^{3}}{3}\,.$ (14)

With such a solution in hand, one can calculate the bounce action of Eq. (1)
numerically. Assuming a local thermal equilibrium around PBHs with a
temperature profile $T(r)$ such that $\lim_{r\rightarrow r_{H}}T(r)=T_{H}$,
the metric does not feature any conical singularity at the horizon, but the
contribution of the effective conical deficit far from the horizon is entirely
contained in its contribution to the Ricci scalar. For simplicity, we will
assume that the interpolation distance $\varepsilon$ is small enough that the
regularisation procedure described in Gregory  _et al._ (2014); Burda  _et
al._ (2015b) is sufficient, and we will use equivalently the result of Eq. (4)
to calculate the action.

Despite the powerful generality of Eq. (4), the calculation of the
corresponding nucleation rate was always restricted to the case of the Hartle-
Hawking vacuum Burda  _et al._ (2015a, 2016, b); Gregory  _et al._ (2014); Dai
_et al._ (2020). Using dimensional analysis to obtain the pre-factor, the
nucleation rate was obtained in Ref. Burda  _et al._ (2015a, 2016, b) by
simply writing

$\Gamma_{\rm FVD}^{\rm HH}\equiv(GM_{+})^{-1}\left(\frac{I_{\rm
b}[T_{H}]}{2\pi}\right)^{1/2}\exp\left(-I_{\rm b}[T_{H}]\right)\,.$ (15)

By noticing that $(GM_{+})^{-1}\sim T_{H}$, one can interpret this result as
the FVD rate in a thermal plasma derived in Ref. Linde (1981). Generalizing
this analogy to the case of a thermal hot spot with temperature $T<T_{H}$, we
thus propose to calculate the nucleation rate using

$\displaystyle\Gamma_{\rm FVD}(T)$ $\displaystyle\approx$ $\displaystyle
T\left(\frac{I_{\rm b}[T]}{2\pi}\right)^{1/2}\exp\left(-I_{\rm
b}[T]\right)\,,$ (16) $\displaystyle\approx$ $\displaystyle
T\left(\frac{I_{\rm b}[T_{H}]}{2\pi}\right)^{1/2}\exp\left(-I_{\rm
b}[T_{H}]\right)\,,$

where we used Eq. (4) to go from the first line to the second. Note that in
this equation, thanks to Eq. (4), the only suppression is due to the fact that
$T<T_{H}$ appears in the pre-factor, as the four-dimensional action does not
depend on the hot-spot temperature value.

Before we derive bounds on PBHs using this result, one needs to discuss what
the actual value of the temperature $T$ should be when calculating the FVD
rate in Eq. 16. In the previous section, we have seen that PBHs with masses
$M\gtrsim 0.8\mathrm{g}$ are surrounded for a long time by a hot spot with
temperature $T_{\rm plateau}\approx 2\times 10^{-4}T_{H}$ before they
evaporate and the hot spot reaches a maximum plateau temperature $T_{\rm
max}\approx 2\times 10^{9}\mathrm{GeV}$ that is independent on the initial PBH
mass. Once the hot spot has formed at $T_{\rm plateau}$, its temperature
remains constant throughout most of the PBH lifetime, enabling the use of the
Euclidean formalism in this regime. Therefore, one can safely use $T_{\rm
plateau}$ to calculate the nucleation rate. However, such a choice is
extremely conservative. Indeed, every PBH heavier than $M_{\star}$ is expected
to evaporate and reach a point where $M=M_{\star}$ and the hot spot
temperature reaches $T_{\rm max}$. Because $T_{\rm max}\geqslant T_{\rm
plateau}$,

$\left.I_{b}[T_{H}]\right|_{M=M_{\star}}\ll\left.I_{b}[T_{H}]\right|_{M>M_{\star}}\,,$
(17)

leading to a much larger FVD rate independent on the PBH mass. A legitimate
concern is whether the Euclidean formalism is reliable given that the PBH and
its environment is dynamical. To ensure this, we restrict ourselves to cases
where the characteristic timescales are much larger than FVD, i.e.
$\Gamma_{\rm FVD}/\Gamma_{\rm ev}\gg 1$.

## Case of the Electroweak Vacuum. —

Let us now apply our results to the case of the electroweak vacuum. In order
to precisely calculate $\Gamma_{\rm FVD}$, we parameterise the Higgs
potential,

$V(\phi)=\lambda_{\text{eff}}(\phi)\phi^{4}/4\,,$ (18)

where

$\lambda_{\text{eff}}=\lambda_{*}+b\left(\ln\frac{\phi}{M_{p}}\right)^{2}+c\left(\ln\frac{\phi}{M_{p}}\right)^{4}\,,$
(19)

in which we determined the parameters $(\lambda_{*},b,c)=(-3.2\times
10^{-3},\,-1.497\times 10^{-6},\,5.42\times 10^{-8})$ by taking the updated
measurements of the top quark mass from the CMS collaboration CMS (2023) and
using the numerical package PyR@te 3 Poole and Thomsen (2019); Sartore and
Schienbein (2021) to solve the renormalization group running at the three-loop
level. To this zero-temperature potential, one needs to add one-loop thermal
corrections:

$\Delta V(\phi,T)=\kappa^{2}T_{\text{plasma}}^{2}\phi^{2}\,,$ (20)

where $\kappa=0.35$ following (See Ref. Strumia (2022) and ref therein). Given
these parameters, the electroweak vacuum appears to be metastable, although
its lifetime is much longer than the age of the Universe in the absence of
evaporating PBHs Degrassi (2014); Alekhin  _et al._ (2012); Espinosa  _et al._
(2015); Buttazzo  _et al._ (2013). In Refs. Gregory  _et al._ (2014); Burda
_et al._ (2015a, 2016); Hayashi  _et al._ (2020); Burda  _et al._ (2015b), the
effect of PBHs on the FVD was used to derive constraints on the amount of PBHs
formed in the early Universe, using the Hartle-Hawking vacuum, and ignoring
thermal corrections to the potential. Here, we propose to use the result of
Eq. (16) to obtain more realistic results.

Figure 2: Evolution of the ratio $\Gamma_{\rm FVD}/\Gamma_{\rm ev}$ for
different values of the PBH mass, using the initial plateau temperature
$T_{\rm plateau}$ to calculate $\Gamma_{\rm FVD}$ (red curve). The gray-shaded
area depicts the region where $M<M_{\star}$, and the red-shaded area the
region where the ratio is larger than unity.

## PBH bounds and Discussion. —

In Fig. 2, we depict the evolution of the ratio $\Gamma_{\rm FVD}/\Gamma_{\rm
ev}$ as a function of the PBH mass. As one can see, $\Gamma_{\rm
FVD}/\Gamma_{\rm ev}$ is larger than one for $M_{\star}<M\lesssim
2\mathrm{g}$. However, $\Gamma_{\rm FVD}$ is exponentially suppressed for
larger PBH masses. Remarkably, the value of this ratio at $M=M_{\star}$ is of
about $5\times 10^{5}\gg 1$. PBHs with large initial masses — for which
$\Gamma_{\rm FVD}(T_{\rm plateau})/\Gamma_{\rm ev}\ll 1$ — thus feature a much
larger bubble nucleation rate once they evaporate and reach the point where
$M=M_{\star}$. The FVD is then a much faster process than the PBH evaporation
at that point. For PBHs with initial masses in the region where $M<M_{\star}$,
the hot spot does not have enough time to form in the first place He  _et al._
(2023). However, PBHs with initial masses larger than $M_{\star}$, will
eventually enter this region when they evaporate. After they pass the
threshold $M=M_{\star}$, such PBHs do not have enough energy to reheat the hot
spot. Therefore the hot spot temperature stays constant at $T=T_{\rm max}$
throughout the end of the evaporation. In this region, we therefore fixed the
temperature at $T_{\rm max}$ in the calculation.

In Ref. Strumia (2022), the effects of thermal corrections were claimed to
become sizeable for $M\lesssim 10^{3}\mathrm{g}$. However, such effects only
become relevant for $M\lesssim 0.1\mathrm{g}<M_{\star}$ in our case and do not
affect any of the conclusions drawn in this letter. Indeed, because we now
consider the hot spot temperature to be $T_{\rm plateau}\approx
10^{-4}T_{\text{H}}$, the effect of thermal corrections is shifted to lower
masses by a factor $\mathcal{O}(10^{-4})$, in agreement with our results.
Another important feature of the bubble profiles we obtain is that they always
have radii much smaller than $r_{\rm plateau}$ and $r_{\rm max}$. This
validates the assumption of homogeneous plateau temperature that we have used
throughout this work. We also checked that the value of the action calculated
using Eq. (1) and substracting the background action and conical deficit
matches exactly with the expression in Eq. (4).

To convert the FVD rate around a black hole into a probability of FVD
throughout cosmic history, one needs to evaluate time $\Delta t$ over which
the FVD rate calculation remains valid during the PBH lifetime. The
corresponding total decay probability can then be evaluated as

$P_{\rm FVD}\equiv 1-e^{-\Gamma_{\rm FVD}\Delta t}\,.$ (21)

First, we have seen that the hot spot temperature stays constant at $T=T_{\rm
plateau}$ during most of the PBH lifetime $t_{\rm ev}$. Therefore, the false
vacuum has about $\Delta t\sim\Gamma_{\rm ev}^{-1}$ to tunnel. In that case,

$P_{\rm FVD}(M)=1-e^{-\Gamma_{\rm FVD}(T_{\rm plateau})/\Gamma_{\rm ev}}\,.$
(22)

Depending on the ratio $\Gamma_{\rm FVD}/\Gamma_{\rm ev}$ given in Fig. 2,
this probability can vary over orders of magnitude. Once the black hole starts
evaporating and eventually reaches $M=M_{\star}$, then we have seen that
$\Gamma_{\rm FVD}/\Gamma_{\rm ev}\approx 5\times 10^{5}$. This ensures that
$P_{\rm FVD}\approx 1$ as long as one considers time scales $\Delta t\lesssim
10^{-6}\times\Gamma_{\rm ev}^{-1}$. Over such a short time, the PBH mass
variation is negligible, and the Euclidean formalism still holds. Finally,
with these decay probabilities, we may now look at the ensuing bounds on the
share of the earlier universe energy which collapsed into PBHs in cosmic
history $\beta_{\rm PBH}=\rho_{\text{PBH}}/\rho_{\rm tot}$. At the earliest,
PBHs formed when an overdensity of size comparable to that of a Hubble patch
collapsed gravitationally. Assuming the Universe to be radiation-dominated
during that time such that $\rho_{\rm tot}=\rho_{\rm rad}$ one can calculate
the Universe’s temperature at formation to be

$T_{f}=\left(\frac{\gamma}{4\pi}\sqrt{\frac{45}{\pi
g_{\star}(T_{f})}}\frac{M_{p}^{3}}{M}\right)^{1/2}\,,$ (23)

and obtain the value of the density fraction at formation

$\beta_{\rm PBH}=\frac{4}{3}\frac{MN_{\rm PBH}H_{0}^{3}}{s_{0}T_{f}}\approx
2\times 10^{-80}N_{\rm PBH}\left(\frac{M}{M_{\star}}\right)^{3/2}\,,$ (24)

where $\gamma=(1/\sqrt{3})^{3}$ is a numerical factor related to the
gravitational collapse Carr (1975); Carr  _et al._ (2010), $N_{\rm PBH}$ is
the total number of PBHs that, if stable, would be contained in a Hubble patch
of size $H_{0}^{-1}$, where $s_{0}\approx 2\times 10^{-38}\mathrm{GeV^{3}}$ is
the entropy density of the Universe today, and $H_{0}\approx 70\
\mathrm{km^{-1}s^{-1}Mpc^{-1}}$ is the Hubble constant today.

Figure 3: Constraints on the PBH energy fraction at formation, $\beta_{\rm
PBH}^{\prime}\equiv\gamma^{1/2}(g_{\star}(T_{f})/106.75)^{-1/4}\beta_{\rm
PBH}$, as a function of its mass. The orange-shaded region represents the most
competitive constraint derived using $T_{\rm max}$ for the hot spot
temperature, whereas the red-shaded region denotes the most conservative bound
using $T_{\rm plateau}$.

Using a Poisson distribution, demanding that the electroweak vacuum has never
decayed in a Hubble patch until today at $95\%$ confidence level is equivalent
to requesting that

$N_{\text{PBH}}\,P_{d}<2.7\,,$ (25)

leading, in the case where $P_{d}\approx 1$, to the constraint

$\beta_{\rm PBH}\lesssim 5\times
10^{-80}\left(\frac{M}{M_{\star}}\right)^{3/2}\,.$ (26)

In Fig. 3, we draw the corresponding constraints both when considering $T_{\rm
plateau}$ and $T_{\rm max}$ for the hot spot temperature. Note that for masses
larger than $M_{\star}$ but smaller than $88$g ($4$g), the temperature $T_{\rm
plateau}$ ($T_{\rm max}$) is smaller than the plasma temperature at the time
of evaporation. In this case, we slightly modified the result of Eq. (16), by
using the plasma temperature at evaporation instead of the hot spot
temperature:

$\Gamma_{\rm FVD}(T)\approx\max\\{T,T_{\rm ev}\\}\left(\frac{I_{\rm
b}[T_{H}]}{2\pi}\right)^{1/2}\exp\left(-I_{\rm b}[T_{H}]\right)\,,$ (27)

where

$T_{\rm ev}\equiv\left(\frac{90}{8\pi^{3}g_{\star}(t_{\rm
ev})}\right)^{1/4}\\!\\!\sqrt{\Gamma_{\rm ev}M_{p}}\,,$ (28)

is the plasma temperature at evaporation. In the regime where such constraints
are relevant, we checked that the ratio $\Gamma_{\rm FVD}/\Gamma_{\rm ev}\gg
1$ guaranteeing the validity of the Euclidean formalism. Note also that
thermal corrections in the action, arising now at $T_{\rm ev}$ rather than
$T_{\rm plateau}$ or $T_{\rm max}$, remain subdominant in this region of
parameter space. For comparison, we indicate other constraints from inflation,
Big Bang Nucleosynthesis, CMB distortion and $\gamma$-rays Carr  _et al._
(2021).

The conclusion one can draw from the results above is twofold: If the
electroweak vacuum is confirmed to be metastable with future SM measurements,
such limits exclude many scenarios predicting a large abundance of evaporating
PBHs in cosmology Inomata  _et al._ (2020); Barman  _et al._ (2022); Bhaumik
_et al._ (2022); Morrison  _et al._ (2019); Gondolo  _et al._ (2020); Bernal
and Zapata (2021a); Jyoti Das  _et al._ (2021); Hamada and Iso (2017); Hooper
and Krnjaic (2021); Chaudhuri and Dolgov (2020); Kitabayashi (2021); Cheek
_et al._ (2022a); Martin  _et al._ (2020a, b); Auclair and Vennin (2021);
Cheek  _et al._ (2022b); Hooper  _et al._ (2019); Cheek  _et al._ (2022c);
Perez-Gonzalez and Turner (2020); Bernal and Zapata (2021b); Datta  _et al._
(2020); Masina (2021); Bernal and Zapata (2021c); Baldes  _et al._ (2020);
Allahverdi  _et al._ (2018); Masina (2020); Lunardini and Perez-Gonzalez
(2020); Dienes  _et al._ (2022); Papanikolaou (2022); Saito and Yokoyama
(2009). In particular, they would exclude entirely the possibility that PBHs
dominate the Universe before evaporating (above the dotted line in Fig. 3).
Alternatively, if such scenarios are confirmed using cosmological data, our
results may be used as a strong indication that new physics is required to
stabilise the electroweak vacuum throughout cosmic history. Finally, we
emphasize that the results presented in this letter can be applied to any
first-order phase transition taking place in cosmology.

## Acknowledgments. —

The authors would like to thank W.-Y. Ai for extremely useful discussions
throughout the realisation of this work, as well as R. Gregory for providing
us with helpful feedback on the manuscript. L. Hamaide is supported by the
King’s Cromwell Scholarship. S. Hu is supported by the Chinese Scholarship
Council and King’s College London. The work of L. Heurtier is supported in
part by the U.K. Science and Technology Facilities Council (STFC) under Grant
ST/P001246/1. A. Cheek is supported by the grant “AstroCeNT: Particle
Astrophysics Science and Technology Centre” carried out within the
International Research Agendas program of the Foundation for Polish Science
financed by the European Union under the European Regional Development Fund.
L. Heurtier would like to thank his recently born daughter, Lilie Ming-Xuan
Heurtier, for providing him with the extreme happiness required to writing
this letter despite the lack of sleep.

## References

  * Hartle and Hawking (1976) J. B. Hartle and S. W. Hawking, Phys. Rev. D 13, 2188 (1976).
  * Gregory  _et al._ (2014) R. Gregory, I. G. Moss, and B. Withers, JHEP 03, 081 (2014), arXiv:1401.0017 [hep-th] .
  * Burda  _et al._ (2015a) P. Burda, R. Gregory, and I. Moss, Phys. Rev. Lett. 115, 071303 (2015a), arXiv:1501.04937 [hep-th] .
  * Burda  _et al._ (2015b) P. Burda, R. Gregory, and I. Moss, JHEP 08, 114 (2015b), arXiv:1503.07331 [hep-th] .
  * Burda  _et al._ (2016) P. Burda, R. Gregory, and I. Moss, JHEP 06, 025 (2016), arXiv:1601.02152 [hep-th] .
  * Hiller  _et al._ (2022) G. Hiller, T. Höhne, D. F. Litim, and T. Steudtner, Phys. Rev. D 106, 115004 (2022), arXiv:2207.07737 [hep-ph] .
  * Unruh (1976) W. G. Unruh, Phys. Rev. D 14, 870 (1976).
  * Kohri and Matsui (2018) K. Kohri and H. Matsui, Phys. Rev. D 98, 123509 (2018), arXiv:1708.02138 [hep-ph] .
  * Hayashi  _et al._ (2020) T. Hayashi, K. Kamada, N. Oshita, and J. Yokoyama, JHEP 08, 088 (2020), arXiv:2005.12808 [hep-th] .
  * Shkerin and Sibiryakov (2021) A. Shkerin and S. Sibiryakov, JHEP 11, 197 (2021), arXiv:2105.09331 [hep-th] .
  * Shkerin and Sibiryakov (2022) A. Shkerin and S. Sibiryakov, JHEP 08, 161 (2022), arXiv:2111.08017 [hep-th] .
  * Strumia (2022) A. Strumia, (2022), arXiv:2209.05504 [hep-ph] .
  * Briaud  _et al._ (2022) V. Briaud, A. Shkerin, and S. Sibiryakov, Phys. Rev. D 106, 125001 (2022), arXiv:2210.08028 [gr-qc] .
  * Das and Hook (2021) S. Das and A. Hook, JHEP 12, 145 (2021), arXiv:2109.00039 [hep-ph] .
  * He  _et al._ (2023) M. He, K. Kohri, K. Mukaida, and M. Yamada, JCAP 01, 027 (2023), arXiv:2210.06238 [hep-ph] .
  * Dai  _et al._ (2020) D.-C. Dai, R. Gregory, and D. Stojkovic, Phys. Rev. D 101, 125012 (2020), arXiv:1909.00773 [hep-ph] .
  * MacGibbon and Webber (1990) J. H. MacGibbon and B. R. Webber, Phys. Rev. D 41, 3052 (1990).
  * MacGibbon (1991) J. H. MacGibbon, Phys. Rev. D 44, 376 (1991).
  * Linde (1981) A. D. Linde, Phys. Lett. B 100, 37 (1981).
  * CMS (2023) (2023), arXiv:2302.01967 [hep-ex] .
  * Poole and Thomsen (2019) C. Poole and A. E. Thomsen, JHEP 09, 055 (2019), arXiv:1906.04625 [hep-th] .
  * Sartore and Schienbein (2021) L. Sartore and I. Schienbein, Comput. Phys. Commun. 261, 107819 (2021), arXiv:2007.12700 [hep-ph] .
  * Degrassi (2014) G. Degrassi, Nuovo Cim. C 037, 47 (2014), arXiv:1405.6852 [hep-ph] .
  * Alekhin  _et al._ (2012) S. Alekhin, A. Djouadi, and S. Moch, Phys. Lett. B 716, 214 (2012), arXiv:1207.0980 [hep-ph] .
  * Espinosa  _et al._ (2015) J. R. Espinosa, G. F. Giudice, E. Morgante, A. Riotto, L. Senatore, A. Strumia, and N. Tetradis, JHEP 09, 174 (2015), arXiv:1505.04825 [hep-ph] .
  * Buttazzo  _et al._ (2013) D. Buttazzo, G. Degrassi, P. P. Giardino, G. F. Giudice, F. Sala, A. Salvio, and A. Strumia, JHEP 12, 089 (2013), arXiv:1307.3536 [hep-ph] .
  * Carr (1975) B. J. Carr, Astrophys. J. 201, 1 (1975).
  * Carr  _et al._ (2010) B. J. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Phys. Rev. D 81, 104019 (2010), arXiv:0912.5297 [astro-ph.CO] .
  * Carr  _et al._ (2021) B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Rept. Prog. Phys. 84, 116902 (2021), arXiv:2002.12778 [astro-ph.CO] .
  * Inomata  _et al._ (2020) K. Inomata, M. Kawasaki, K. Mukaida, T. Terada, and T. T. Yanagida, Phys. Rev. D 101, 123533 (2020), arXiv:2003.10455 [astro-ph.CO] .
  * Barman  _et al._ (2022) B. Barman, D. Borah, S. Jyoti Das, and R. Roshan, (2022), arXiv:2212.00052 [hep-ph] .
  * Bhaumik  _et al._ (2022) N. Bhaumik, A. Ghoshal, R. K. Jain, and M. Lewicki, (2022), arXiv:2212.00775 [astro-ph.CO] .
  * Morrison  _et al._ (2019) L. Morrison, S. Profumo, and Y. Yu, JCAP 1905, 005 (2019), arXiv:1812.10606 [astro-ph.CO] .
  * Gondolo  _et al._ (2020) P. Gondolo, P. Sandick, and B. Shams Es Haghi, Phys. Rev. D 102, 095018 (2020), arXiv:2009.02424 [hep-ph] .
  * Bernal and Zapata (2021a) N. Bernal and O. Zapata, Phys. Lett. B 815, 136129 (2021a), arXiv:2011.02510 [hep-ph] .
  * Jyoti Das  _et al._ (2021) S. Jyoti Das, D. Mahanta, and D. Borah, (2021), arXiv:2104.14496 [hep-ph] .
  * Hamada and Iso (2017) Y. Hamada and S. Iso, PTEP 2017, 033B02 (2017), arXiv:1610.02586 [hep-ph] .
  * Hooper and Krnjaic (2021) D. Hooper and G. Krnjaic, Phys. Rev. D 103, 043504 (2021), arXiv:2010.01134 [hep-ph] .
  * Chaudhuri and Dolgov (2020) A. Chaudhuri and A. Dolgov, (2020), arXiv:2001.11219 [astro-ph.CO] .
  * Kitabayashi (2021) T. Kitabayashi, (2021), arXiv:2101.01921 [hep-ph] .
  * Cheek  _et al._ (2022a) A. Cheek, L. Heurtier, Y. F. Perez-Gonzalez, and J. Turner, Phys. Rev. D 105, 015022 (2022a), arXiv:2107.00013 [hep-ph] .
  * Martin  _et al._ (2020a) J. Martin, T. Papanikolaou, and V. Vennin, JCAP 01, 024 (2020a), arXiv:1907.04236 [astro-ph.CO] .
  * Martin  _et al._ (2020b) J. Martin, T. Papanikolaou, L. Pinol, and V. Vennin, JCAP 05, 003 (2020b), arXiv:2002.01820 [astro-ph.CO] .
  * Auclair and Vennin (2021) P. Auclair and V. Vennin, JCAP 02, 038 (2021), arXiv:2011.05633 [astro-ph.CO] .
  * Cheek  _et al._ (2022b) A. Cheek, L. Heurtier, Y. F. Perez-Gonzalez, and J. Turner, (2022b), arXiv:2207.09462 [astro-ph.CO] .
  * Hooper  _et al._ (2019) D. Hooper, G. Krnjaic, and S. D. McDermott, JHEP 08, 001 (2019), arXiv:1905.01301 [hep-ph] .
  * Cheek  _et al._ (2022c) A. Cheek, L. Heurtier, Y. F. Perez-Gonzalez, and J. Turner, Phys. Rev. D 105, 015023 (2022c), arXiv:2107.00016 [hep-ph] .
  * Perez-Gonzalez and Turner (2020) Y. F. Perez-Gonzalez and J. Turner, (2020), arXiv:2010.03565 [hep-ph] .
  * Bernal and Zapata (2021b) N. Bernal and O. Zapata, JCAP 03, 007 (2021b), arXiv:2010.09725 [hep-ph] .
  * Datta  _et al._ (2020) S. Datta, A. Ghosal, and R. Samanta, (2020), arXiv:2012.14981 [hep-ph] .
  * Masina (2021) I. Masina, (2021), arXiv:2103.13825 [gr-qc] .
  * Bernal and Zapata (2021c) N. Bernal and O. Zapata, JCAP 03, 015 (2021c), arXiv:2011.12306 [astro-ph.CO] .
  * Baldes  _et al._ (2020) I. Baldes, Q. Decant, D. C. Hooper, and L. Lopez-Honorez, JCAP 08, 045 (2020), arXiv:2004.14773 [astro-ph.CO] .
  * Allahverdi  _et al._ (2018) R. Allahverdi, J. Dent, and J. Osinski, Phys. Rev. D97, 055013 (2018), arXiv:1711.10511 [astro-ph.CO] .
  * Masina (2020) I. Masina, Eur. Phys. J. Plus 135, 552 (2020), arXiv:2004.04740 [hep-ph] .
  * Lunardini and Perez-Gonzalez (2020) C. Lunardini and Y. F. Perez-Gonzalez, JCAP 08, 014 (2020), arXiv:1910.07864 [hep-ph] .
  * Dienes  _et al._ (2022) K. R. Dienes, L. Heurtier, F. Huang, D. Kim, T. M. P. Tait, and B. Thomas, (2022), arXiv:2212.01369 [astro-ph.CO] .
  * Papanikolaou (2022) T. Papanikolaou, JCAP 10, 089 (2022), arXiv:2207.11041 [astro-ph.CO] .
  * Saito and Yokoyama (2009) R. Saito and J. Yokoyama, Phys. Rev. Lett. 102, 161101 (2009), [Erratum: Phys.Rev.Lett. 107, 069901 (2011)], arXiv:0812.4339 [astro-ph] .


∎

11institutetext: Pulak Swain 22institutetext: School of Basic Sciences, Indian
Institute of Technology Bhubaneswar
Tel.: +917978772457
22email<EMAIL_ADDRESS>33institutetext: Akshay Kumar Ojha 44institutetext:
School of Basic Sciences, Indian Institute of Technology Bhubaneswar

# Robust Portfolio Optimization under Ambiguous Chance Constraints

Pulak Swain  Akshay Kumar Ojha  Pulak Swain  Akshay Kumar Ojha

(Received: date / Accepted: date)

###### Abstract

In this paper, we discuss the ambiguous chance constrained based portfolio
optimization problems, in which the perturbations associated with the input
parameters are stochastic in nature, but their distributions are not known
precisely. We consider two different families of perturbation distributions–
one is when only upper $\&$ lower bounds of mean values are known to us, and
the second one is along with the mean bounds, we also have the knowledge of
the standard deviations of the perturbations. We derive the safe convex
approximations of such chance constrained portfolio problems by using some
suitable generating functions such as a piecewise linear function, an
exponential function, and a piecewise quadratic function. These safe
approximations are the robust counterparts to our ambiguous chance constrained
problem and they are computationally tractable due to the convex nature of
these approximations.

###### Keywords:

Markowitz Model Robust Counterpart Chance Constraint Safe Convex Approximation
Piecewise Linear Approximation Bernstein Approximation Piecewise Quadratic
Approximation

##  1 Introduction

The problem of optimal asset allocation is an important aspect of investment
studies and it has been truly identified by Harry Markowitz in 1952 Markowitz
J 1952 . Markowitz’s studies show the importance of diversification of
investment through the optimization models Markowitz 1959 . In such problems,
the prime objective of the investor is to find the right balance between two
conflicting components- risk and return. To visualize this opposing behavior
of the above two components, a trade-off curve can be drawn which is known as
the efficient frontier. Markowitz’s mean-variance portfolio model is defined
by using the expected returns and covariance of returns of the individual
assets as the input parameters. For a mathematical overview, let’s assume that
the investor wants to invest his total wealth among $n$ assets where the
return of any arbitrary $i^{th}$ asset over a time period $t$ is given by
$r_{it}\ (t=1,2,\dots,T)$. Let $\mu_{i}$ and $\sigma_{ij}$ are respectively
the expected return of $i^{th}$ asset and the covariance of return between
$i$-$j$ pair of assets, which are defined as,

$\displaystyle\begin{aligned}
&\mu_{i}=\dfrac{1}{T}\sum_{t=1}^{T}r_{it},\quad\sigma_{ij}=\dfrac{1}{T}\sum_{t=1}^{T}(r_{it}-\mu_{i})(r_{jt}-\mu_{j})\end{aligned}$

Let $x_{i}$ be the weight allocated for $i^{th}$ asset in the portfolio. Then
the expected portfolio return ($\mu_{p}$) and the variance of portfolio return
($\sigma_{p}^{2}$) can be calculated as,

$\displaystyle\begin{aligned}
&\mu_{p}=\sum_{i=1}^{n}\mu_{i}x_{i},\quad\sigma_{p}^{2}=\sum_{i=1}^{n}\sum_{j=1}^{n}\sigma_{ij}x_{i}x_{j}\end{aligned}$

Markowitz portfolio model minimizes the variance of portfolio return by
constraining the expected return to become more than or equal to a target
return ($\tau$). Moreover, the model assumes the weights to be non-negative
and the sum of the weights to be $1$. Mathematically the Markowitz portfolio
model is given by,

$\displaystyle\begin{aligned}
&\min_{x_{i}}&&\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\sigma_{ij}x_{i}x_{j}\\\
&\textrm{s.t.:}&&\sum_{j=1}^{n}\mu_{j}x_{j}\geq
1,\quad\sum_{j=1}^{n}{x_{j}}=1,\quad{x_{j}}\geq 0,\
j=1,2,\dots,n.\end{aligned}$ (1)

One of the biggest challenges in portfolio optimization is to get the exact
input parameters such as expected asset returns and the covariance of asset
returns. The main reason for this is the several economic factors involved in
the market which can affect these input parameters. This leads to the presence
of some uncertainty factor in these parameters, which if ignored can
potentially affect the optimal asset allocation. To address this issue,
several approaches such as sensitivity analysis, dynamic programming, fuzzy
programming, etc. are used in the field of optimization. However, in the last
two decades, robust optimization approach has been widely used by researchers
due to its ability to produce ”completely immunized against uncertainty”
solutions. This approach was first studied by El Ghaoui and Lebret Ghaoui 1997
, in which they worked on the uncertain least square problems. Then El Ghaoui
et al. Ghaoui 1998  used the robust optimization approach for uncertain
semidefinite programs and provided sufficient conditions for the robust
solutions. Later on Ben-Tal and Nemirovski Ben-Tal 1999  studied the robust
optimization for linear programming problems for ellipsoidal uncertainty and
showed the tractability of the robust counterpart. Another work of Ben-Tal and
Nemirovski Ben-Tal 2000  was based on some uncertain linear programming
problems from NETLIB collection and their results showed that the robust
solutions nearly lose nothing in optimality. In the literature, there are some
more studies on the geometrical shape based uncertainties such as box,
ellipsoidal, polyhedral, etc. The robust optimization approach was first used
in the uncertain portfolio problems when Goldfarb and Iyenger Goldfarb 2003
introduced some uncertain structures for the market parameters and
reformulated the robust portfolio problems into second order cone programming.
Tütüncü and Koenig Tutuncu 2004  formulated the robust counterparts when the
expected return vector and the covariance matrix of asset returns are defined
by lower and upper bounds and they also showed the robust efficient frontiers.
Fliege and Werner Fliege 2014  studied the application of robust
multiobjective optimization on uncertain mean-variance portfolio problems.
Sehgal and Mehra Sehgal 2019  proposed robust portfolio optimization models
for reward–risk ratios by using Omega, semi-mean absolute deviation ratio, and
weighted stable tail adjusted return ratio and also showed the efficiency of
these models.

On the other hand, the chance constrained based optimization problems Biswal
1998  is another type of uncertain optimization problem in which the uncertain
parameters are stochastic in nature. Moreover, chance constraints basically
indicate the probability of uncertain constraint(s) to be satisfied, where the
perturbations act as random variables. Calafiore and El Ghaoui Calafiore 2006
studied distributional robust chance constrained linear problems with radially
distributed perturbations. They showed that the linear chance constraints can
be converted explicitly into convex second-order cone constraints so that they
can be efficiently solvable. Erdoğan and Iyengar Erdogan 2006  studied the
ambiguous chance constrained problem which is approximated by a robust sampled
problem. Nemirovski and Shapiro Nemirovski 2007  studied the convex
approximations for the ambiguous chance constraints. In particular, they
derived the Bernstein approximation and their approximation was convex and
efficiently solvable. Yanıkoğlu and den Hertog Yanıkoğlu 2013  used the
available historical data of uncertain parameters to get the safe
approximations of ambiguous chance constraints and their approach has the
advantage that it can also easily handle joint chance constraints. Bertsimas
et al. Bertsimas 2018  derived safe approximations to nonlinear robust
individual chance constraints by using statistical hypothesis tests. Zhang et
al. Zhang 2018  solved distributionally robust chance-constrained binary
programs by equivalently reformulating them as $0-1$ second order cone
programs.

The main contribution of the paper is as follows: Sometimes the uncertainty
present in the portfolio optimization problems is stochastic in nature and we
don’t know the exact distribution of the uncertain parameters. Such problems
are in general difficult to solve, as in those problems we not only need to
handle the chance constraints, but also we need to take care of the ambiguous
nature of probability distributions present in the chance constraints. So the
challenge is to convert such problems into deterministic as well as
computationally tractable optimization problems. In this paper we use some
generating functions to get the safe convex approximations to such problems.
We consider two different families of perturbation distributions– (i) when we
know the bounds of mean perturbation values and (ii) when we know the mean and
standard deviations of perturbation values. The safe convex approximations
which we get for the family-i are quadratic programming problems and the
approximations we get for the family-ii are quadratically constrained
quadratic programming problems, which are computationally tractable. Moreover,
we implemented the results into a real life stock market problem to find out
the robust optimal allocation. Also, the risk-return efficient frontiers are
drawn for all the safe approximation cases.

The remainder of the paper is organized as follows: Section 2 presents some
definitions and basic concepts regarding robust optimization and the chance
constrained based portfolio optimization problems. Section 3 constitutes the
discussion on generating function based safe convex approximation and then the
safe approximations for several families of perturbation distributions are
derived in 4. Further, a numerical example consisting of an Indian stock
market problem is given in Section 5. Finally, some concluding remarks are
provided in Section 6.

##  2 Preliminaries

###  2.1 Uncertainty in Optimization Problems and the Robust Optimization
Approach

When the constant coefficients in an optimization problem are not fixed
values; but they perturb around some nominal values, such problems are called
uncertain optimization problems. Mathematically this can be represented as
Bertsimas 2011 ,

$\displaystyle\begin{aligned} &\min_{\bm{x}}&&f(\bm{x},\bm{u})\\\
&\textrm{s.t.:}&&c(\bm{x},\bm{u})\geq
0,\quad\forall\bm{u}=\bm{u}^{(0)}+\displaystyle{\sum_{j=1}^{n}}\zeta_{j}\bm{u}^{(j)}\in\mathscr{U}\end{aligned}$
(2)

where $\bm{x}$ is the vector of decision variables and $\bm{u}$ is the vector
of uncertain variables, $\bm{u}^{(0)}$ is the vector of nominal values of
uncertain parameters and $\bm{u}^{(j)}$ (for $j=1,2,\dots,n$) are the basic
shifts.
Now we will go through some definitions related to robust optimization.

######  Definition 1 (Robust Feasible Solution)

Ben-Tal 2009  A vector $\bm{x}$ is said to be the robust feasible solution of
the uncertain problem (2) if it satisfies the uncertain constraints for all
realizations of the uncertain set $\mathscr{U}$, that is, if $\bm{x}$
satisfies

$\displaystyle\begin{aligned} c(\bm{x},\bm{u})\leq
0,\quad\forall\bm{u}\in\mathscr{U}\end{aligned}$

######  Definition 2 (Robust Value)

Ben-Tal 2009  For a given candidate solution $\bm{x}$, the robust value
$\widehat{f}(\bm{x})$ of the objective in problem (2) is the largest value of
$f(\bm{x},\bm{u})$ over all realizations of the data from the uncertain set,
that is

$\displaystyle\begin{aligned}
\widehat{f}(\bm{x})=\sup_{\bm{u}\in\mathscr{U}}f(\bm{x},\bm{u})\end{aligned}$

It is to be noted that if the problem (2) were a maximization problem, then
the robust value would have been
$\widehat{f}(\bm{x})=\inf_{\bm{u}\in\mathscr{U}}f(\bm{x},\bm{u})$.

######  Definition 3 (Robust Counterpart)

Ben-Tal 2009  The robust counterpart (RC) of the uncertain problem (2) is the
optimization problem

$\displaystyle\begin{aligned}
&\min_{\bm{x}}&&\left\\{\sup_{\bm{u}\in\mathscr{U}}f(\bm{x},\bm{u})\right\\}\\\
&\text{s.t.:}&&c(\bm{x},\bm{u})\leq\bm{0},\quad\forall\bm{u}\in\mathscr{U}\end{aligned}$

of minimizing the robust value of the objective over all the robust feasible
solutions to the uncertain problem.

######  Definition 4 (Robust Optimal Solution)

Ben-Tal 2009  The solution of the RC problem is said to be the robust optimal
solution of the uncertain problem (2).

###  2.2 Ambiguous Chance Constraint

Consider the uncertain constraint given in the problem (2),

$\displaystyle\begin{aligned} c(\bm{x},\bm{u})\geq
0,\quad\bm{u}=\bm{u}^{(0)}+\displaystyle{\sum_{j=1}^{n}}\zeta_{j}\bm{u}^{(j)}\end{aligned}$
(3)

In the ideal scenario, we would like to find a solution $\bm{x}$ which
satisfies the constraint (3) for all realizations of the uncertain set. But
when the data perturbation is stochastic in nature, it is not always possible
to find such a solution. So we look for a candidate solution that satisfies
(3) for ”nearly all” realizations of $\bm{\zeta}$. Thus we want our constraint
to be satisfied with a predefined reliability level. Then the probability
chance constraint to (3) is introduced which is given by,

$\displaystyle\begin{aligned} {Prob}_{\bm{\zeta}\sim
P}\left\\{\bm{\zeta}:c(\bm{x},\bm{u})\geq
0,\quad\bm{u}=\bm{u}^{(0)}+\displaystyle{\sum_{j=1}^{n}}\zeta_{j}\bm{u}^{(j)}\right\\}\geq\beta\end{aligned}$
(4)

Eq. (4) indicates that the constraint (3) would be satisfied with a
probability of at least $\beta$, where $\beta$ may be taken as close to 1.
However, sometimes the probability distribution of $\bm{\zeta}$ is not known
precisely; all we know is that the distribution belongs to a given family
$\mathcal{P}$ of probability distributions on $\mathbb{R}^{n}$. In such cases,
we need to convert the problem (4) into an ambiguously chance constrained
version which is given by, One of the biggest drawbacks of these ambiguous
chance constrained problems is the tractability issues of the robust
counterparts. The problem becomes NP-hard even for the cases when P is simple.

###  2.3 Uncertain Portfolio Problem with Chance Constraint

Consider the portfolio optimization problem (LABEL:eq*) in vector form given
by,

$\displaystyle\begin{aligned}
&\min_{\bm{x}}&&\frac{1}{2}\bm{x}^{\top}\bm{\Sigma}\bm{x}\\\
&\textrm{s.t.:}&&\bm{\mu}^{\top}\bm{x}\geq\tau,\quad\bm{e}^{\top}\bm{x}=1,\quad\bm{x}\geq\bm{0}\end{aligned}$
(5)

where $\bm{\mu}=\begin{bmatrix}\mu_{1}\\\ \mu_{2}\\\ \vdots\\\
\mu_{n}\end{bmatrix}$,
$\bm{\Sigma}=\begin{bmatrix}\sigma_{11}&\sigma_{12}&\dots&\sigma_{1n}\\\
\sigma_{21}&\sigma_{22}&\dots&\sigma_{2n}\\\ \vdots&\vdots&\ddots&\vdots\\\
\sigma_{n1}&\sigma_{n2}&\dots&\sigma_{nn}\\\ \end{bmatrix}$,
$\bm{x}=\begin{bmatrix}x_{1}\\\ x_{2}\\\ \vdots\\\ x_{n}\end{bmatrix}$,
$\bm{e}=\begin{bmatrix}1\\\ 1\\\ \vdots\\\ 1\end{bmatrix}$,
$\bm{0}=\begin{bmatrix}0\\\ 0\\\ \vdots\\\ 0\end{bmatrix}$. The uncertainty in
portfolio optimization problems in general occurs due to the input parameters
– expected returns and covariance of returns of the assets. However, the
studies show that the uncertainty in the covariance of asset returns does not
affect the optimal solution so much Pulak 2021 . So in this study, we assume
that the uncertainty occurs only in the expected returns of assets. Then the
uncertain portfolio model can be written as,

$\displaystyle\begin{aligned}
&\min_{\bm{x}}&&\frac{1}{2}\bm{x}^{\top}\bm{\Sigma}\bm{x}\\\
&\textrm{s.t.:}&&\bm{\mu}^{\top}\bm{x}\geq\tau,\quad\bm{\mu}\in\mathscr{U}_{\bm{\mu}},\\\
&&&\bm{e}^{\top}\bm{x}=1,\quad\bm{x}\geq\bm{0}\end{aligned}$ (6)

where $\mathscr{U}_{\bm{\mu}}$ is the uncertainty set associated with
$\bm{\mu}$ given by,

$\displaystyle\begin{aligned}
\mathscr{U}_{\mu}=\left\\{\bm{\mu}:\bm{\mu}\equiv\bm{\mu}(\bm{\zeta})=\bm{\mu}^{(0)}+\sum_{j=1}^{n}\zeta_{j}\bm{\mu}^{(j)}\right\\}\end{aligned}$

Now the model (LABEL:eq3a) can be written with ambiguous chance constraint as,

$\displaystyle\begin{aligned}
&\min_{\bm{x}}&&\frac{1}{2}\bm{x}^{\top}\bm{\Sigma}\bm{x}\\\
&\textrm{s.t.:}&&\forall P\in\mathcal{P}:{Prob}_{\bm{\zeta}\sim
P}\left\\{[\bm{\mu}(\bm{\zeta})]^{\top}\bm{x}\geq\tau\right\\}\geq\beta,\\\
&&&\bm{e}^{\top}\bm{x}=1,\quad\bm{x}\geq\bm{0}\end{aligned}$ (7)

The chance constraint in problem (LABEL:eq4) can also be written as,

$\displaystyle\begin{aligned} &{Prob}_{\bm{\zeta}\sim
P}\left\\{[\bm{\mu}(\bm{\zeta})]^{\top}\bm{x}<\tau\right\\}\leq
1-\beta\end{aligned}$ (8)

##  3 Generating-Function-Based Safe Convex Approximation for Ambiguous Chance
Constrained Portfolio Optimization Model

The major drawbacks of chance constrained problems are:

  * •

It is difficult to evaluate with high accuracy the probability on the left
side of a chance constraint

* •

In most of the cases the feasible set of the chance constraint is non-convex,
due to which the optimization problem becomes very difficult to solve.

One way to overcome these issues is to replace the chance constraint in the
model (LABEL:eq4) with its computationally tractable safe approximation and
then solve that approximation to get the robust solution.

######  Definition 5 (Safe Convex Approximation)

Ben-Tal 2009  Let $\\{\bm{\mu}^{(j)}\\}_{j=0}^{n}$, $Prob$, $\beta$ be the
data of chance constraint (LABEL:eq5) and let $S$ be a system of convex
constraints on $\bm{x}$ and additional variables $\bm{v}$. Then $S$ is said to
be a safe convex approximation of chance constraint (LABEL:eq5) if the
$\bm{x}$ component of every feasible solution $(\bm{x};\bm{v})$ of $S$ is
feasible for the chance constraint.

Moreover, a safe convex approximation $S$ of (LABEL:eq5) is said to be
computationally tractable, if the convex constraints forming $S$ are
efficiently computable. Nemirovski and Shapiro Nemirovski 2007  suggested a
generating function based approach to obtain the computationally tractable
safe convex approximation, where the generating function $\gamma(\cdot)$ must
satisfy the following properties: (i) $\gamma(\cdot):\mathbb{R}\to\mathbb{R}$
is a non-negative non-decreasing function. (ii) $\gamma(\cdot)$ is convex.
(iii) $\gamma(0)\geq 1$ and $\gamma(t)\to 0$ as $t\to-\infty$. Let’s denote
the left hand side of the chance constraint (LABEL:eq5) by
$p(\bm{\mu},\bm{x})$ and the constraint is given by,

$\displaystyle\begin{aligned} &p(\bm{\mu},\bm{x})\equiv{Prob}_{\bm{\zeta}\sim
P}\left\\{[\bm{\mu}(\bm{\zeta})]^{\top}\bm{x}<\tau\right\\}\leq
1-\beta\end{aligned}$ (9)

The constraint can also be written as,

$\displaystyle\begin{aligned}
&E\left[\chi(\tau-\bm{\mu}(\bm{\zeta})^{\top}\bm{x})\right]\leq
1-\beta,\end{aligned}$ (10)

where $\chi(t)$ is the characteristic function given by,

$\displaystyle\begin{aligned} &\chi(t)=\left\\{\begin{array}[]{ll}0&\mbox{ for
}t\leq 0\\\ 1&\mbox{ for }t>0\end{array}\right.\end{aligned}$

and $E[.]$ is the expectation of the function defined as,

$\displaystyle\begin{aligned}
&E\left[\chi(\tau-\bm{\mu}(\bm{\zeta})^{\top}\bm{x})\right]=\int\chi(\tau-\bm{\mu}(\bm{\zeta})^{\top}\bm{x})\
dP(\bm{\zeta}).\end{aligned}$

Since $\chi(\cdot)$ is not a convex function, it follows that we can not
convert (LABEL:eq6a) into a convex constraint of $\bm{\mu}$ by using the
function $\chi(\cdot)$. To find out a convex approximation for the chance
constraint (LABEL:eq6), let $\gamma(t)$ be a convex function on the axis which
is everywhere $\geq$ $\chi(t)$ and it satisfies the properties (i)-(iii)
mentioned before. Then we have,

$\displaystyle\begin{aligned}
&p(\bm{\mu},\bm{x})=E\left[\chi(\tau-\bm{\mu}(\bm{\zeta})^{\top}\bm{x})\right]\leq
E\left[\gamma(\tau-\bm{\mu}(\bm{\zeta})^{\top}\bm{x})\right]=\psi(\bm{\mu},\bm{x})\text{
(say)}\end{aligned}$ (11)

Now let $\psi^{+}(\bm{\mu},\bm{x})$ is some upper bound of
$\psi(\bm{\mu},\bm{x})$, that is,
$\psi(\bm{\mu},\bm{x})\leq\psi^{+}(\bm{\mu},\bm{x})~{}\forall\bm{\mu}$. Then
the constraint

$\displaystyle\begin{aligned} &\psi^{+}(\bm{\mu},\bm{x})\leq
1-\beta\end{aligned}$ (12)

can be a safe convex approximation of the chance constraint (LABEL:eq5). When
the function $\psi^{+}(\cdot)$ is efficiently computable, the safe
approximation (LABEL:eq7a) is tractable. This observation can be converted
into the following result:

###### Theorem 3.1

For the ambiguous chance constraint (LABEL:eq5) of an uncertain portfolio
optimization problem (LABEL:eq4), a generating function $\gamma(\cdot)$ can be
used to write the safe convex approximation given by (LABEL:eq7a), where
$\psi^{+}(\bm{\mu},\bm{x})\geq
E\left[\gamma(\tau-\bm{\mu}(\bm{\zeta})^{\top}\bm{x})\right]$. Moreover, this
approximation is tractable, provided that $\psi^{+}(\cdot)$ is efficiently
computable.

It is to be noted that, to reduce the conservatism of the above approximation,
we can also do scaling as follows: From the definition of $\chi(t)$, it can be
observed that, $\chi(t)=\chi(\alpha^{-1}t)$, $\forall\alpha>0$, which implies
$\gamma(\alpha^{-1}t)\geq\chi(\alpha^{-1}t)=\chi(t)$, $\forall t$. So from eq.
(LABEL:eq5) we get,

$\displaystyle\begin{aligned}
&p(\bm{\mu},\bm{x})\geq\psi^{+}(\alpha^{-1}\bm{\mu},\bm{x})\end{aligned}$ (13)

and the constraint

$\displaystyle\begin{aligned} &\psi^{+}(\alpha^{-1}\bm{\mu},\bm{x})\leq
1-\beta\end{aligned}$ (14)

is also a safe convex approximation for all $\alpha>0$. This implies
$\alpha\psi^{+}(\alpha^{-1}\bm{\mu},\bm{x})-\alpha(1-\beta)\leq 0$ is a safe
approximation for any $\alpha>0$. So the infimum of the left hand side of the
function can also be used to write the safe convex approximation, which is
given by,

$\displaystyle\begin{aligned}
\inf_{\alpha>0}\left[\alpha\psi^{+}(\alpha^{-1}\bm{\mu},\bm{x})-\alpha(1-\beta)\right]\leq
0\end{aligned}$ (15)

When the function $\psi^{+}(\cdot)$ is efficiently computable, the left hand
side of (15) is also efficiently computable. Then the approximation (15) is
tractable. In the next section, we obtain the robust counterparts of the
problem (LABEL:eq4) for several families $\mathcal{P}$ of perturbation
distributions.

##  4 Safe Convex Approximations for Several Families of Perturbation
Distributions

We consider two different families of perturbation distributions– (i) the
family whose upper $\&$ lower mean bounds are only known to us, and (ii) the
family whose upper $\&$ lower mean bounds, and the standard deviations are
known to us. To derive the safe approximations we use some generating
functions $\gamma(t)$ which are greater or equal to $\chi(t)$ everywhere and
satisfy the properties given earlier.

###  4.1 Family of Perturbation Distributions with Known Upper $\&$ Lower Mean
Bounds

In this case we consider the family of
$\bm{\zeta}=[\zeta_{1},\zeta_{2},\dots,\zeta_{n}]^{\top}$ with the upper $\&$
lower bounds of each $\zeta_{j}$, $(j=1,2,\dots,n)$ are known to be
$m_{j}^{U}$ $\&$ $m_{j}^{L}$ respectively. For deriving the safe convex
approximation we use a piecewise linear function and the exponential function
as the generating functions.

####  4.1.1 Piecewise Linear Generating Function Based Approximation

The function

$\displaystyle\begin{aligned}
&\gamma(t)=[1+t]_{+}=\left\\{\begin{array}[]{ll}0&\mbox{ for }t\leq-1\\\
1+t&\mbox{ for }t>-1\end{array}\right.\end{aligned}$ (16)

is everywhere greater than or equal to $\chi(t)$. Further, it is a non-
negative and non-decreasing convex function satisfying $\gamma(0)\geq 1$ and
$\lim_{t\to-\infty}\gamma(t)=0$ and hence can be considered as a generating
function.

Figure 1: Piecewise Linear Generating Function

Now since $E\left[\chi(\tau-\bm{\mu}(\bm{\zeta})^{\top}\bm{x})\right]\leq
E\left[\gamma(\tau-\bm{\mu}(\bm{\zeta})^{\top}\bm{x})\right]$ by eq.
(LABEL:eq7), it follows that we can replace the chance constraint (LABEL:eq6a)
by the constraint given by,

$\displaystyle\begin{aligned} &\qquad E_{\bm{\zeta\sim
P}}\left[\gamma(\tau-\bm{\mu}(\bm{\zeta})^{\top}x)\right]\leq
1-\beta\end{aligned}$ (17)

By using the definition of $\gamma(.)$ given by (LABEL:eq10a), the constraint
(LABEL:eq10b) can be written as,

$\displaystyle\begin{aligned} &\qquad E_{\bm{\zeta\sim
P}}\left[\max\left\\{0,1+\tau-\bm{\mu}(\bm{\zeta})^{\top}x\right\\}\right]\leq
1-\beta\end{aligned}$ (18)

The constraint (LABEL:eq10c) can be replaced with two equivalent constraints,

$\displaystyle E_{\bm{\zeta\sim P}}\left[0\right]\leq 1-\beta$ (19a)
$\displaystyle E_{\bm{\zeta\sim
P}}\left[1+\tau-\bm{\mu}(\bm{\zeta})^{\top}x\right]\leq 1-\beta$ (19b)

Here the constraint (19a) is trivial, as by the assumption of chance
constraint problems the value of confidence level $\beta$ lies between $0$ and
$1$, and hence we have $1-\beta\geq 0$. Therefore now we only need to proceed
with the constraint (19b). The left hand side of (19b) can be simplified in
component form as,

$\displaystyle\begin{aligned} E_{\bm{\zeta\sim
P}}\left[1+\tau-\bm{\mu}(\bm{\zeta})^{\top}x\right]&=E_{{\zeta_{j}\sim
P_{j}}}\left[1+\tau-\sum_{j=1}^{n}{\mu_{j}^{(0)}x_{j}}-\sum_{j=1}^{n}{(\mu_{j}^{(j)}x_{j})\zeta_{j}}\right]\\\
&=1+\tau-\sum_{j=1}^{n}{\mu_{j}^{(0)}x_{j}}-E_{{\zeta_{j}\sim
P_{j}}}\left[\sum_{j=1}^{n}{(\mu_{j}^{(j)}x_{j})\zeta_{j}}\right]\\\
&=1+\tau-\sum_{j=1}^{n}{\mu_{j}^{(0)}x_{j}}-\sum_{j=1}^{n}{(\mu_{j}^{(j)}x_{j})E_{{\zeta_{j}\sim
P_{j}}}\left[\zeta_{j}\right]}\end{aligned}$ (20)

It is to be noted that $\zeta_{j}$’s are the random variables in eq. (20), so
the terms without $\zeta_{j}$ are treated as constants. Now our objective is
to find out an upper bound of the right hand side of eq. (20), which would be
our $\psi^{+}(\bm{\mu},\bm{x})$ given in Theorem 3.1. Since the upper $\&$
lower mean bounds each perturbation $\zeta_{j}$ in the family are known to be
$m_{j}^{U}$ $\&$ $m_{j}^{L}$ respectively, that is, $m_{j}^{L}\leq
E_{{\zeta_{j}\sim P_{j}}}\left[\zeta_{j}\right]\leq m_{j}^{U}$, it follows
from eq. (20) that,

$\displaystyle\begin{aligned} E_{\bm{\zeta\sim
P}}\left[1+\tau-\bm{\mu}(\bm{\zeta})^{\top}x\right]\leq
1+\tau-\sum_{j=1}^{n}{\mu_{j}^{(0)}x_{j}}-\sum_{j=1}^{n}{(\mu_{j}^{(j)}x_{j})m_{j}^{L}}\end{aligned}$
(21)

Thus by taking
$\psi^{+}(\bm{\mu},\bm{x})=1+\tau-\sum_{j=1}^{n}{\mu_{j}^{(0)}x_{j}}-\sum_{j=1}^{n}{(\mu_{j}^{(j)}x_{j})m_{j}^{L}}$,
the condition in Theorem 3.1 is satisfied and the constraint

$\displaystyle\begin{aligned}
&1+\tau-\sum_{j=1}^{n}{\mu_{j}^{(0)}x_{j}}-\sum_{j=1}^{n}{(\mu_{j}^{(j)}x_{j})m_{j}^{L}}\leq
1-\beta\end{aligned}$

can be replaced with the chance constraint. Therefore a safe approximation of
the robust counterpart to the portfolio model (LABEL:eq4) is given by,

$\displaystyle\begin{aligned}
&\min_{x_{i}}&&\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\sigma_{ij}x_{i}x_{j}\\\
&\textrm{s.t.:}&&1+\tau-\sum_{j=1}^{n}{\mu_{j}^{(0)}x_{j}}-\sum_{j=1}^{n}{(\mu_{j}^{(j)}x_{j})m_{j}^{L}}\leq
1-\beta,\\\ &&&\sum_{j=1}^{n}{x_{j}}=1,\quad{x_{j}}\geq 0,\
j=1,2,\dots,n.\end{aligned}$ (22)

The robust counterpart model (LABEL:eq10g) is a quadratic programming problem,
so it is computationally tractable. The above result can be represented as the
following theorem:

###### Theorem 4.1

Consider the ambiguous chance constrained portfolio optimization model
(LABEL:eq4), in which the probability distribution of the perturbations
$\zeta_{j}$, $(j=1,2,\dots,n)$ are not known precisely, and the only
information regarding each $\zeta_{j}$ is that it belongs to a family of
perturbations having known upper $\&$ lower mean bounds $m_{j}^{U}$ $\&$
$m_{j}^{L}$ respectively. Then a safe tractable convex approximation to this
model can be a quadratic programming problem given by,

$\displaystyle\begin{aligned}
&\min_{x_{i}}&&\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\sigma_{ij}x_{i}x_{j}\\\
&\textrm{s.t.:}&&1+\tau-\sum_{j=1}^{n}{\mu_{j}^{(0)}x_{j}}-\sum_{j=1}^{n}{(\mu_{j}^{(j)}x_{j})m_{j}^{L}}\leq
1-\beta,\\\ &&&\sum_{j=1}^{n}{x_{j}}=1,\quad{x_{j}}\geq 0,\
j=1,2,\dots,n.\end{aligned}$

####  4.1.2 Bernstein Approximation

One of the first studies of the Bernstein approximation scheme was done by
Nemirovski and Shapiro Nemirovski 2007 . The scheme has two assumptions: (i)
$\mathcal{P}$ can be written as,

$\displaystyle\mathcal{P}=\left\\{P_{1}\times P_{2}\times\dots\times
P_{n}:P_{j}\in\mathcal{P}_{j},\ j=1,2,\dots,n\right\\}$

where for $j=1,2,\dots,n$, the distributions ${P}_{j}$ are independent and
belong to the families $\mathcal{P}_{j}$ of distributions. (ii) The moment
generating functions $M_{j}(t)=E[e^{t\zeta_{j}}]=\int e^{tz}dP_{j}(z)$
$(j=1,2,\dots,n)$ are finite valued for all $t\in\mathbb{R}$ and are
efficiently computable. The generating function for this approximation is
given by,

$\displaystyle\begin{aligned} \gamma(t)=e^{t}\end{aligned}$ (23)

which is everywhere greater than or equal to $\chi(t)$. Further, it can be
easily verified that $e^{t}$ satisfies the properties of generating function.

Figure 2: Exponential Generating Function

By using the definition of $\gamma(.)$ given by (23), the constraint
(LABEL:eq10b) can be written as,

$\displaystyle\begin{aligned} &E_{{\bm{\zeta}\sim
P}}\left[e^{\tau-\bm{\mu}(\bm{\zeta})^{\top}\bm{x}}\right]\leq
1-\beta\end{aligned}$ (24)

Taking the logarithm in both sides of the constraint (LABEL:eq10i) we get,

$\displaystyle\begin{aligned} &E_{{\bm{\zeta}\sim
P}}\left[\tau-\bm{\mu}(\bm{\zeta})^{\top}\bm{x}\right]\leq\log(1-\beta)\end{aligned}$
(25)

Now the left hand side of (LABEL:eq10j) can be simplified in component form
as,

$\displaystyle\begin{aligned} E_{\bm{\zeta}\sim
P}\left[{\tau-\bm{\mu}(\bm{\zeta})^{\top}\bm{x}}\right]&=E_{\bm{\zeta_{j}}\sim
P_{j}}\left[{\tau-\sum_{j=1}^{n}\mu_{j}^{(0)}x_{j}-\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})\zeta_{j}}\right]\\\
&=\tau-\sum_{j=1}^{n}\mu_{j}^{(0)}x_{j}-\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})E_{\bm{\zeta_{j}}\sim
P_{j}}\left[\zeta_{j}\right]\end{aligned}$ (26)

Now our objective is to find out an upper bound of the right hand side of eq.
(26), which would be our $\psi^{+}(\bm{\mu},\bm{x})$ given in Theorem 3.1.
Since we have, $m_{j}^{L}\leq E_{{\zeta_{j}\sim
P_{j}}}\left[\zeta_{j}\right]\leq m_{j}^{U}$, it follows from eq. (26) that,

$\displaystyle\begin{aligned} E_{\bm{\zeta}\sim
P}\left[{\tau-\bm{\mu}(\bm{\zeta})^{\top}\bm{x}}\right]&\leq\tau-\sum_{j=1}^{n}\mu_{j}^{(0)}x_{j}-\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})m_{j}^{L}\end{aligned}$
(27)

Then we get
$\psi^{+}(\bm{\mu},\bm{x})=\tau-\sum_{j=1}^{n}\mu_{j}^{(0)}x_{j}-\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})m_{j}^{L}$
and thus the safe approximation of the chance constraint becomes,

$\displaystyle\tau-\sum_{j=1}^{n}\mu_{j}^{(0)}x_{j}-\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})m_{j}^{L}\leq\log(1-\beta)$

Therefore a safe approximation of the robust counterpart to the portfolio
model (LABEL:eq4) is given by,

$\displaystyle\begin{aligned}
&\min_{x_{i}}&&\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\sigma_{ij}x_{i}x_{j}\\\
&\textrm{s.t.:}&&\tau-\sum_{j=1}^{n}{\mu_{j}^{(0)}x_{j}}-\sum_{j=1}^{n}{(\mu_{j}^{(j)}x_{j})m_{j}^{L}}\leq\log(1-\beta),\\\
&&&\sum_{j=1}^{n}{x_{j}}=1,\quad{x_{j}}\geq 0,\ j=1,2,\dots,n.\end{aligned}$
(28)

The robust counterpart model (LABEL:eq12) is a quadratic programming problem,
so it is computationally tractable. The above result can be represented as the
following theorem:

###### Theorem 4.2

Consider the ambiguous chance constrained portfolio optimization model
(LABEL:eq4), in which the probability distribution of the perturbations
$\zeta_{j}$, $(j=1,2,\dots,n)$ are not known precisely, and the only
information regarding each $\zeta_{j}$ is that it belongs to a family of
perturbations having known upper $\&$ lower mean bounds $m_{j}^{U}$ $\&$
$m_{j}^{L}$ respectively. Then a safe tractable convex approximation to this
model can be a quadratic programming problem given by,

$\displaystyle\begin{aligned}
&\min_{x_{i}}&&\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\sigma_{ij}x_{i}x_{j}\\\
&\textrm{s.t.:}&&\tau-\sum_{j=1}^{n}{\mu_{j}^{(0)}x_{j}}-\sum_{j=1}^{n}{(\mu_{j}^{(j)}x_{j})m_{j}^{L}}\leq\log(1-\beta),\\\
&&&\sum_{j=1}^{n}{x_{j}}=1,\quad{x_{j}}\geq 0,\ j=1,2,\dots,n.\end{aligned}$

###  4.2 Family of Perturbation Distributions with Known Upper $\&$ Lower Mean
Bounds and Standard Deviations

In this case we consider the family of
$\bm{\zeta}=[\zeta_{1},\zeta_{2},\dots,\zeta_{n}]^{\top}$ with the upper $\&$
lower bounds of each $\zeta_{j}$, $(j=1,2,\dots,n)$ are known to be
$m_{j}^{U}$ $\&$ $m_{j}^{L}$ respectively and the standard deviation is known
to be $s_{j}$. For deriving the safe convex approximation we use a piecewise
quadratic generating function.

####  4.2.1 Piecewise Quadratic Generating Function Based Approximation

The function

$\displaystyle\begin{aligned} \gamma(t)=\left\\{\begin{array}[]{ll}0&\mbox{
for }t\leq-1\\\ (1+t)^{2}&\mbox{ for }t>-1\end{array}\right.\end{aligned}$
(29)

is everywhere greater than or equal to $\chi(t)$. Further, it is a non-
negative and non-decreasing convex function satisfying $\gamma(0)\geq 1$ and
$\lim_{t\to-\infty}\gamma(t)=0$ and hence can be considered as a generating
function.

Figure 3: Piecewise Quadratic Generating Function

We can also rewrite $\gamma(t)$ given in (29) as,

$\displaystyle\begin{aligned}
&\gamma(t)=(1+t)\cdot[1+t]_{+}=(1+t)\cdot\max\left\\{0,1+t\right\\}\end{aligned}$
(30)

Similar to the previous two approximations, let’s begin with the constraint

$\displaystyle\begin{aligned} &E_{\bm{\zeta\sim
P}}\left[\gamma(\tau-\bm{\mu}(\bm{\zeta})^{\top}x)\right]\leq
1-\beta\end{aligned}$ (31)

By the use of (LABEL:eq12b), the constraint (LABEL:eq12c) becomes,

$\displaystyle\begin{aligned} &E_{\bm{\zeta\sim
P}}\left[(1+\tau-\bm{\mu}(\bm{\zeta})^{\top}x)\cdot\max\left\\{0,(1+\tau-\bm{\mu}(\bm{\zeta})^{\top}x)\right\\}\right]\leq
1-\beta\end{aligned}$ (32)

Now here arises two cases for
$\max\left\\{0,1+\tau-\bm{\mu}(\bm{\zeta})^{\top}x\right\\}$. Case I: When
$\max\left\\{0,1+\tau-\bm{\mu}(\bm{\zeta})^{\top}x\right\\}=0$ Then the
constraint (LABEL:eq12d) becomes,

$\displaystyle\begin{aligned} E_{\bm{\zeta\sim P}}\left[0\right]\leq
1-\beta\end{aligned}$

Case II: When
$\max\left\\{0,1+\tau-\bm{\mu}(\bm{\zeta})^{\top}x\right\\}=1+\tau-\bm{\mu}(\bm{\zeta})^{\top}x$
Then the constraint (LABEL:eq12d) becomes,

$\displaystyle\begin{aligned} E_{\bm{\zeta\sim
P}}\left[1+\tau-\bm{\mu}(\bm{\zeta})^{\top}x\right]^{2}\leq
1-\beta\end{aligned}$

So the constraint (LABEL:eq12d) can be replaced by two equivalent constraints,

$\displaystyle E_{\bm{\zeta\sim P}}\left[0\right]\leq 1-\beta$ (33a)
$\displaystyle E_{\bm{\zeta\sim
P}}\left[1+\tau-\bm{\mu}(\bm{\zeta})^{\top}x\right]^{2}\leq 1-\beta$ (33b)

Here the constraint (33a) is trivial, as by the assumption of chance
constraint problems we have $1-\beta\geq 0$. Therefore now we only need to
proceed with the constraint (33b). Now the left hand side of (33b) can be
simplified in component form as,

$\displaystyle\begin{aligned} E_{\bm{\zeta}\sim
P}\left[1+\tau-\bm{\mu}(\bm{\zeta})^{\top}\bm{x}\right]^{2}&=E_{\zeta_{j}\sim
P_{j}}\left[1+\tau-\left(\sum_{j=1}^{n}\mu_{j}^{(0)}x_{j}+\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})\zeta_{j}\right)\right]^{2}\\\
&=E_{\zeta_{j}\sim
P_{j}}\left[(1+\tau-\sum_{j=1}^{n}\mu_{j}^{(0)}x_{j})-\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})\zeta_{j}\right]^{2}\\\
&=E_{\zeta_{j}\sim
P_{j}}\left[\left(1+\tau-\sum_{j=1}^{n}\mu_{j}^{(0)}x_{j}\right)^{2}+\left(\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})\zeta_{j}\right)^{2}\right.\\\
&\qquad\qquad\qquad\qquad\qquad\left.-2\left(1+\tau-\sum_{j=1}^{n}\mu_{j}^{(0)}x_{j}\right)\cdot\left(\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})\zeta_{j}\right)\right]\end{aligned}$
(34)

Since in our problem the perturbations $\zeta_{j}$’s are the random variables,
the terms without $\zeta_{j}$ can be treated as constants. So the eq. (34) can
be written as,

$\displaystyle\begin{aligned} &E_{\bm{\zeta}\sim
P}\left[1+\tau-\bm{\mu}(\bm{\zeta})^{\top}\bm{x}\right]^{2}=&&(1+\tau-\sum_{j=1}^{n}\mu_{j}^{(0)}x_{j})^{2}+E_{\zeta_{j}\sim
P_{j}}\left[\left(\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})\zeta_{j}\right)\right]^{2}\\\
&&&-2\left(1+\tau-\sum_{j=1}^{n}\mu_{j}^{(0)}x_{j}\right)\cdot
E_{\zeta_{j}\sim
P_{j}}\left[\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})\zeta_{j}\right]\end{aligned}$
(35)

Now our objective is to find out an upper bound of the right hand side of eq.
(LABEL:eq18), which would be our $\psi^{+}(\bm{\mu},\bm{x})$ given in Theorem
3.1. For this, we use the following lemma.

###### Lemma 1

Let $\zeta_{j},\ j=1,2,\dots,n$ be independently distributed random variables,
each with its upper $\&$ lower mean bounds are respectively $m_{j}^{U}$ $\&$
$m_{j}^{L}$ and its standard deviation is $s_{j}$. Then we have,

$\displaystyle\begin{aligned}
&(i)&&\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})m_{j}^{L}\leq E_{\zeta_{j}\sim
P_{j}}\left[\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})\zeta_{j}\right]\leq\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})m_{j}^{U},\\\
&&&\text{ and, }\\\
&(ii)&&\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})^{2}s_{j}^{2}+\left[\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})m_{j}^{L}\right]^{2}\leq
E_{\zeta_{j}\sim
P_{j}}\left[\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})\zeta_{j}\right]^{2}\leq\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})^{2}s_{j}^{2}+\left[\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})m_{j}^{U}\right]^{2}\end{aligned}$

###### Proof

We first denote the random variable
$\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})\zeta_{j}$ by $Y$. So our objective is to
find out the bounds of $E\left[Y\right]$ and $E\left[Y^{2}\right]$

Now the mean of $Y$ is given by,

$\displaystyle\begin{aligned}
E\left[Y\right]&=\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})\cdot E_{\zeta_{j}\sim
P_{j}}\left[\zeta_{j}\right]\end{aligned}$ (36)

Using the fact that $m_{j}^{L}\leq E_{\zeta_{j}\sim
P_{j}}\left[\zeta_{j}\right]\leq m_{j}^{U}$ in (36), we get,

$\displaystyle\begin{aligned} \sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})m_{j}^{L}\leq
E\left[Y\right]\leq\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})m_{j}^{U},\end{aligned}$
(37)

which proves (i). Now the variance of $Y$ is given by,

$\displaystyle\begin{aligned}
Var\left[Y\right]&=Var\left[\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})\zeta_{j}\right]\\\
&=\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})^{2}\cdot
Var\left[\zeta_{j}\right]+2\cdot\displaystyle{\sum_{\begin{subarray}{c}i,j\\\
i\neq j\end{subarray}}}(\mu_{i}^{(i)}x_{i})(\mu_{j}^{(j)}x_{j})\cdot
Cov(\zeta_{i},\zeta_{j})\end{aligned}$ (38)

Since $\zeta_{j}$’s are independently distributed random variables, we have
$Cov(\zeta_{i},\zeta_{j})=0$, which reduces (38) to,

$\displaystyle\begin{aligned}
Var\left[Y\right]&=\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})^{2}s_{j}^{2}\end{aligned}$
(39)

Again from the variance formula we have,

$\displaystyle\begin{aligned}
E\left[Y^{2}\right]=Var\left[Y\right]+\left(E\left[Y\right]\right)^{2}\end{aligned}$
(40)

Using (37) and (39) in (40), we obtain

$\displaystyle\begin{aligned}
&\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})^{2}s_{j}^{2}+\left[\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})m_{j}^{L}\right]^{2}\leq
E\left[Y\right]^{2}\leq\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})^{2}s_{j}^{2}+\left[\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})m_{j}^{U}\right]^{2},\end{aligned}$
(41)

and this proves (ii). ∎ Now applying the results (i) and (ii) of Lemma 1 in
eq. (LABEL:eq18) we get,

$\displaystyle\begin{aligned} E_{\bm{\zeta}\sim
P}\left[1+\tau-\bm{\mu}(\bm{\zeta})^{\top}\bm{x}\right]^{2}\leq&\left(1+\tau-\sum_{j=1}^{n}\mu_{j}^{(0)}x_{j}\right)^{2}+\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})^{2}s_{j}^{2}+\left(\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})m_{j}^{U}\right)^{2}\\\
&\qquad-2(1+\tau-\sum_{j=1}^{n}\mu_{j}^{(0)}x_{j})\cdot\left(\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})m_{j}^{L}\right)\end{aligned}$
(42)

Then we get
$\psi^{+}(\bm{\mu},\bm{x})=\left(1+\tau-\sum_{j=1}^{n}\mu_{j}^{(0)}x_{j}\right)^{2}+\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})^{2}s_{j}^{2}+\left(\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})m_{j}^{U}\right)^{2}-2(1+\tau-\sum_{j=1}^{n}\mu_{j}^{(0)}x_{j})\cdot\left(\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})m_{j}^{L}\right)$
and thus the safe approximation of the chance constraint becomes,

$\displaystyle\begin{aligned}
&\left(1+\tau-\sum_{j=1}^{n}\mu_{j}^{(0)}x_{j}\right)^{2}+\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})^{2}s_{j}^{2}+\left(\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})m_{j}^{U}\right)^{2}-2(1+\tau-\sum_{j=1}^{n}\mu_{j}^{(0)}x_{j})\cdot\left(\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})m_{j}^{L}\right)\\\
&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\leq
1-\beta\end{aligned}$

Therefore a safe approximation of the robust counterpart to the portfolio
model (LABEL:eq4) is given by,

$\displaystyle\begin{aligned}
&\min_{x_{i}}&&\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\sigma_{ij}x_{i}x_{j}\\\
&\textrm{s.t.:}&&\left(1+\tau-\sum_{j=1}^{n}\mu_{j}^{(0)}x_{j}\right)^{2}+\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})^{2}s_{j}^{2}+\left(\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})m_{j}^{U}\right)^{2}\\\
&&&\qquad\qquad\qquad-2(1+\tau-\sum_{j=1}^{n}\mu_{j}^{(0)}x_{j})\cdot\left(\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})m_{j}^{L}\right)\leq
1-\beta,\\\ &&&\sum_{j=1}^{n}{x_{j}}=1,\quad{x_{j}}\geq 0,\
j=1,2,\dots,n.\end{aligned}$ (43)

The robust counterpart model (LABEL:eq22a) is a quadratically constrained
quadratic programming (QCQP) problem, so it is computationally tractable. The
above result can be represented as the following theorem:

###### Theorem 4.3

Consider the ambiguous chance constrained portfolio optimization model
(LABEL:eq4), in which the probability distribution of the perturbations
$\zeta_{j}$, $(j=1,2,\dots,n)$ are not known precisely, and the only
information regarding each $\zeta_{j}$ is that it belongs to a family of
perturbations having known upper $\&$ lower mean bounds, and the standard
deviation $m_{j}^{U}$ $\&$ $m_{j}^{L}$, and $s_{j}$ respectively. Then a safe
tractable convex approximation to this model can be a quadratic programming
problem given by,

$\displaystyle\begin{aligned}
&\min_{x_{i}}&&\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\sigma_{ij}x_{i}x_{j}\\\
&\textrm{s.t.:}&&\left(1+\tau-\sum_{j=1}^{n}\mu_{j}^{(0)}x_{j}\right)^{2}+\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})^{2}s_{j}^{2}+\left(\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})m_{j}^{U}\right)^{2}\\\
&&&\qquad\qquad\qquad-2(1+\tau-\sum_{j=1}^{n}\mu_{j}^{(0)}x_{j})\cdot\left(\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})m_{j}^{L}\right)\leq
1-\beta,\\\ &&&\sum_{j=1}^{n}{x_{j}}=1,\quad{x_{j}}\geq 0,\
j=1,2,\dots,n.\end{aligned}$

###  4.3 Extension to Higher Orders

In Section 4.1, the considered perturbation family has the information
regarding first-moment statistics of the distribution, that is, the mean. On
the other hand, the perturbation family in Section 4.2 has information
regarding the first two moments of the distribution, that is, the mean and the
variance (or the standard deviation). A piecewise linear function
(LABEL:eq10a) is used to derive the safe convex approximation for the former
family, whereas a piecewise quadratic function (29) is used for the latter
family. Inductively, if we have the information regarding the first $n$
moments of the perturbation distribution, then we can use the piecewise
generating function of degree $n$ given by,

$\displaystyle\begin{aligned} \gamma(t)=\left\\{\begin{array}[]{ll}0&\mbox{
for }t\leq-1\\\ (1+t)^{n}&\mbox{ for }t>-1\end{array}\right.\end{aligned}$
(44)

to derive the safe convex approximation. This is because, while deriving the
left hand side of the chance constraint componentwise by using the above
generating function, we obtain

$\displaystyle\begin{aligned} &E_{\bm{\zeta}\sim
P}\left[1+\tau-\bm{\mu}(\bm{\zeta})^{\top}\bm{x}\right]^{n}=&&E_{\zeta_{j}\sim
P_{j}}\left[(1+\tau-\sum_{j=1}^{n}\mu_{j}^{(0)}x_{j})-\sum_{j=1}^{n}(\mu_{j}^{(j)}x_{j})\zeta_{j}\right]^{n}\end{aligned}$
(45)

After applying the binomial expansion to the right hand side of eq.
(LABEL:eq22aaa), we require the information regarding the values of
$E_{\zeta_{j}\sim P_{j}}[\zeta_{j}]$, $E_{\zeta_{j}\sim
P_{j}}[\zeta_{j}^{2}]$, $\dots$, $E_{\zeta_{j}\sim P_{j}}[\zeta_{j}^{n}]$ to
derive the safe convex approximation of the uncertain model. Therefore the
generating function (44) can be applied for the family of perturbations, whose
first $n$ moments’ information is known.

##  5 Numerical Example

In this section, we consider an Indian stock market problem for our study. The
objective is to find out the optimal allocation by using the safe
approximation models which we derive in the last section. Finally, we compare
the obtained results.

###  5.1 Data Description

The stock price data of three major sectors of India such as Nifty Bank, Nifty
Infra, and Nifty IT are taken into consideration for our numerical example.
The data are collected for the period June 2017 to May 2022 from
https://finance.yahoo.com. The quarterly returns of the sectors are derived
from their stock price data. Using these we calculate the input parameters for
the portfolio model namely, the expected returns of the sectors and the
covariance of returns between each pair of sectors. The calculated input
parameters are the nominal values, which are given in Table 1.

Table 1: Nominal Values of the Input Parameters Sectors | Expected Returns (in $\%$) | Covariance of Returns (in $\%$)
---|---|---
|  | Nifty Bank | Nifty Infra | Nifty IT
Nifty Bank | $2.609$ | $24.126$ | $-1.460$ | $11.032$
Nifty Infra | $-1.430$ | $-1.460$ | $8.237$ | $0.461$
Nifty IT | $6.329$ | $11.032$ | $0.461$ | $18.034$

###  5.2 Problem Formulation

First, we construct the nominal problem with the input parameters given in
Table 1. The problem is given by,

$\displaystyle\begin{aligned}
&\min_{x_{i}}&&\frac{1}{2}\left[24.126x_{1}^{2}+8.237x_{2}^{2}+18.034x_{3}^{2}+2\cdot(-1.460)x_{1}x_{2}+2\cdot
11.032x_{1}x_{3}+2\cdot 0.461x_{2}x_{3}\right]\\\
&s.t.:&&2.609x_{1}-1.430x_{2}+6.329x_{3}\geq\tau,\\\
&&&x_{1}+x_{2}+x_{3}=1,\quad x_{1},x_{2},x_{3}\geq 0\end{aligned}$ (51)

Taking the uncertainty in expected returns into account, we aim to formulate
the problem with ambiguous chance constraint. Let $\bm{\zeta}=[\zeta_{1}\
\zeta_{2}\ \zeta_{3}]^{\top}$ be the perturbation vector associated with the
expected returns. Here the nominal expected return vector is given by,

$\displaystyle\bm{{\bm{\mu}^{(0)}}}=\begin{bmatrix}\mu_{1}^{(0)}\\\
\mu_{2}^{(0)}\\\ \mu_{3}^{(0)}\\\ \end{bmatrix}=\begin{bmatrix}2.609\\\
-1.430\\\ 6.329\\\ \end{bmatrix}$

Let the basic shifts of the perturbations are assumed to be,

$\displaystyle\bm{{\bm{\mu}^{(1)}}}=\begin{bmatrix}\mu_{1}^{(1)}\\\
\mu_{2}^{(1)}\\\ \mu_{3}^{(1)}\\\ \end{bmatrix}=\begin{bmatrix}0.2\\\ 0\\\
0\\\ \end{bmatrix},\quad\bm{{\bm{\mu}^{(2)}}}=\begin{bmatrix}\mu_{1}^{(2)}\\\
\mu_{2}^{(2)}\\\ \mu_{3}^{(2)}\\\ \end{bmatrix}=\begin{bmatrix}0\\\ 0.1\\\
0\\\ \end{bmatrix},\quad\bm{{\bm{\mu}^{(3)}}}=\begin{bmatrix}\mu_{1}^{(3)}\\\
\mu_{2}^{(3)}\\\ \mu_{3}^{(3)}\\\ \end{bmatrix}=\begin{bmatrix}0\\\ 0\\\
0.3\\\ \end{bmatrix}.$

Then the uncertain expected return of the portfolio is given by,

$\displaystyle\begin{aligned}
&[\bm{\mu}(\bm{\zeta})]^{\top}\bm{x}&&=\left[\bm{\mu}^{(0)}+\sum_{j=1}^{3}\zeta_{j}\bm{\mu}^{(j)}\right]^{\top}\bm{x}\\\
&&&=\begin{bmatrix}2.609+0.2\zeta_{1}\\\ -1.430+0.1\zeta_{2}\\\
6.329+0.3\zeta_{3}\\\ \end{bmatrix}^{\top}\begin{bmatrix}x_{1}\\\ x_{2}\\\
x_{3}\\\ \end{bmatrix}\\\
&&&=2.609x_{1}-1.430x_{2}+6.329x_{3}+0.2x_{1}\zeta_{1}+0.1x_{2}\zeta_{2}+0.3x_{3}\zeta_{3}\end{aligned}$

Therefore the ambiguous chance constrained problem for $\beta=0.95$ is given
by,

$\displaystyle\begin{aligned}
&\min_{x_{i}}&&\frac{1}{2}\left[24.126x_{1}^{2}+8.237x_{2}^{2}+18.034x_{3}^{2}+2\cdot(-1.460)x_{1}x_{2}+2\cdot
11.032x_{1}x_{3}+2\cdot 0.461x_{2}x_{3}\right]\\\ &s.t.:&&\forall
P\in\mathcal{P}:{Prob}_{\bm{\zeta}\sim
P}\left\\{2.609x_{1}-1.430x_{2}+6.329x_{3}+0.2x_{1}\zeta_{1}+0.1x_{2}\zeta_{2}+0.3x_{3}\zeta_{3}\geq\tau\right\\}\geq
0.95\\\ &&&x_{1}+x_{2}+x_{3}=1,\quad x_{1},x_{2},x_{3}\geq 0\end{aligned}$
(52)

We aim to solve the above ambiguous chance constrained problem for the
families of perturbation distributions by using the safe approximations
discussed in Section 4.

###  5.3 Optimal Solutions for Several Families of Perturbation Distributions

####  5.3.1 Family of Perturbation Distributions with Known Upper $\&$ Lower
Mean Bounds

Consider the family $\mathcal{P}$ of perturbation distributions, whose upper
$\&$ lower mean bounds are respectively known as,

$\displaystyle\begin{aligned} \bm{{m}}^{U}=\begin{bmatrix}m_{1}^{U}\\\
m_{2}^{U}\\\ m_{3}^{U}\\\ \end{bmatrix}=\begin{bmatrix}0.3\\\ 0.2\\\ 0.1\\\
\end{bmatrix},\quad\bm{{m}}^{L}=\begin{bmatrix}m_{1}^{L}\\\ m_{2}^{L}\\\
m_{3}^{L}\\\ \end{bmatrix}=\begin{bmatrix}-0.3\\\ -0.2\\\ -0.1\\\
\end{bmatrix}\end{aligned}$ (53)

a. Piecewise Linear Generating Function Based Approximation Using Theorem 4.1
we get the safe convex approximation with the piecewise linear generating
function as,

$\displaystyle\begin{aligned}
&\min_{x_{i}}&&\frac{1}{2}\left[24.126x_{1}^{2}+8.237x_{2}^{2}+18.034x_{3}^{2}+2\cdot(-1.460)x_{1}x_{2}+2\cdot
11.032x_{1}x_{3}+2\cdot 0.461x_{2}x_{3}\right]\\\
&s.t.:&&1+\tau-(2.609x_{1}-1.430x_{2}+6.329x_{3})-(0.2x_{1}m_{1}^{L}+0.1x_{2}m_{2}^{L}+0.3x_{3}m_{3}^{L})\leq
0.05,\\\ &&&x_{1}+x_{2}+x_{3}=1,\quad x_{1},x_{2},x_{3}\geq 0,\end{aligned}$
(54)

where the values of $m_{1}^{L}$, $m_{2}^{L}$, $m_{3}^{L}$ are given in eq.
(53). We solve this problem for different values of $\tau$ varying from $1.5$
to $3.5$ and calculate the optimal portfolio risks. Such choice of $\tau$ is
because of the fact that the portfolio is a combination of all the individual
assets and thus there is a high possibility of the portfolio return to lie
somewhere near the average of the asset returns. The obtained optimal results
are provided in Table 2.

Table 2: Optimal Results for Piecewise Linear Generating Function Based Approximation Target Return($\tau$) | Optimal Allocation | Optimal Portfolio Risk
---|---|---
| Nifty Bank | Nifty Infra | Nifty IT |
$1.5$ | $0.0979$ | $0.4493$ | $0.4528$ | $3.3142$
$1.7$ | $0.0891$ | $0.4278$ | $0.4831$ | $3.4685$
$1.9$ | $0.0803$ | $0.4062$ | $0.5134$ | $3.6382$
$2.1$ | $0.0716$ | $0.3847$ | $0.5438$ | $3.8231$
$2.3$ | $0.0628$ | $0.3631$ | $0.5741$ | $4.0232$
$2.5$ | $0.0540$ | $0.3415$ | $0.6045$ | $4.2386$
$2.7$ | $0.0452$ | $0.3200$ | $0.6348$ | $4.4693$
$2.9$ | $0.0364$ | $0.2984$ | $0.6652$ | $4.7152$
$3.1$ | $0.0276$ | $0.2769$ | $0.6955$ | $4.9763$
$3.3$ | $0.0189$ | $0.2553$ | $0.7259$ | $5.2528$
$3.5$ | $0.0101$ | $0.2337$ | $0.7562$ | $5.5444$

Using these results, we plot the optimal portfolio allocation graph in Fig. 4,
where the black, red, and blue respectively indicate the weights of Nifty
Bank, Nifty Infra, and Nifty IT in the optimal portfolio, when the target
return varies from 1.5 to 3.5. Moreover, the tradeoff between target returns
and their corresponding optimal risks is illustrated by the efficient
frontier, which is given in Fig. 5.

Figure 4: Portfolio Allocation for the Piecewise Linear Generating Function
Based Approximation Figure 5: Efficient Frontier of the Portfolio for the
Piecewise Linear Generating Function Based Approximation

b. Bernstein Approximation Using Theorem 4.2 the safe convex approximation for
Bernstein approximation is given by,

$\displaystyle\begin{aligned}
&\min_{x_{i}}&&\frac{1}{2}\left[24.126x_{1}^{2}+8.237x_{2}^{2}+18.034x_{3}^{2}+2\cdot(-1.460)x_{1}x_{2}+2\cdot
11.032x_{1}x_{3}+2\cdot 0.461x_{2}x_{3}\right]\\\
&s.t.:&&\tau-(2.609x_{1}-1.430x_{2}+6.329x_{3})-(0.2x_{1}m_{1}^{L}+0.1x_{2}m_{2}^{L}+0.3x_{3}m_{3}^{L})\leq\log(0.05),\\\
&&&x_{1}+x_{2}+x_{3}=1,\quad x_{1},x_{2},x_{3}\geq 0,\end{aligned}$ (68)

where the values of $m_{1}^{L}$, $m_{2}^{L}$, $m_{3}^{L}$ are given in eq.
(53). We solve this problem for the values of $\tau$ varying from $1.5$ to
$3.5$ and calculate the optimal portfolio risks. The obtained optimal results
are provided in Table 3.

Table 3: Optimal Results for the Bernstein Approximation Target Return($\tau$) | Optimal Allocation | Optimal Portfolio Risk
---|---|---
| Nifty Bank | Nifty Infra | Nifty IT |
$1.5$ | $0.0081$ | $0.2288$ | $0.7631$ | $5.6133$
$1.7$ | $0.0000$ | $0.2069$ | $0.7931$ | $5.9237$
$1.9$ | $0.0000$ | $0.1811$ | $0.8189$ | $6.2503$
$2.1$ | $0.0000$ | $0.1553$ | $0.8447$ | $6.5939$
$2.3$ | $0.0000$ | $0.1295$ | $0.8705$ | $6.9543$
$2.5$ | $0.0000$ | $0.1037$ | $0.8963$ | $7.3316$
$2.7$ | $0.0000$ | $0.0779$ | $0.9221$ | $7.7257$
$2.9$ | $0.0000$ | $0.0520$ | $0.9480$ | $8.1368$
$3.1$ | $0.0000$ | $0.0262$ | $0.9738$ | $8.5648$
$3.3$ | $0.0000$ | $0.0004$ | $0.9996$ | $9.0096$
$3.5$ | $0.0000$ | $0.0000$ | $1.0000$ | $9.5839$

Using these results, we plot the optimal allocation graph and the efficient
frontier, which are illustrated in Fig. 6 and Fig. 7 respectively.

Figure 6: Portfolio Allocation for the Bernstein Approximation Figure 7:
Efficient Frontier of the Portfolio for the Bernstein Approximation

####  5.3.2 Family of Perturbation Distributions with Known Upper $\&$ Lower
Mean Bounds and Standard Deviations

In this case, we consider the family $\mathcal{P}$ of perturbation
distributions, whose upper $\&$ lower mean bounds are given by eq. (53) and
the standard deviations are given by,

$\displaystyle\begin{aligned} \bm{s}=\begin{bmatrix}s_{1}\\\ s_{2}\\\ s_{3}\\\
\end{bmatrix}=\begin{bmatrix}0.1\\\ 0.1\\\ 0.1\\\ \end{bmatrix}\end{aligned}$
(82)

a. Piecewise Quadratic Approximation Using Theorem 4.3 we get the safe convex
approximation with the piecewise quadratic generating function as,

$\displaystyle\begin{aligned}
&\min_{x_{i}}&&\frac{1}{2}\left[24.126x_{1}^{2}+8.237x_{2}^{2}+18.034x_{3}^{2}+2\cdot(-1.460)x_{1}x_{2}+2\cdot
11.032x_{1}x_{3}+2\cdot 0.461x_{2}x_{3}\right]\\\
&s.t.:&&[(1+\tau-(2.609x_{1}-1.430x_{2}+6.329x_{3})]^{2}+[(0.2x_{1})^{2}s_{1}^{2}+(0.1x_{2})^{2}s_{2}^{2}+(0.3x_{3})^{2}s_{3}^{2}]\\\
&&&+[(0.2x_{1})m_{1}^{U}+(0.1x_{2})m_{2}^{U}+(0.3x_{3})m_{3}^{U}]^{2}\\\
&&&-2\left[1+\tau-(2.609x_{1}-1.430x_{2}+6.329x_{3})\right]\cdot\left[0.2x_{1}(m_{1}^{L})+0.1x_{2}(m_{2}^{L})+0.3x_{3}(m_{3}^{L})\right]\leq
0.05,\\\ &&&x_{1}+x_{2}+x_{3}=1,\quad x_{1},x_{2},x_{3}\geq 0,\end{aligned}$
(83)

where the values of $m_{j}^{L}$ $\&$ $m_{j}^{U}$, ($j=1,2,3$) are given in eq.
(53), and the values of $s_{j}$, ($j=1,2,3$) are given in eq. (82). We solve
this problem for the values of $\tau$ varying from $1.5$ to $3.5$ and
calculate the optimal portfolio risks. The obtained optimal results are
provided in Table 4.

Table 4: Optimal Results for Piecewise Quadratic Function Based Approximation Target Return($\tau$) | Optimal Allocation | Optimal Portfolio Risk
---|---|---
| Nifty Bank | Nifty Infra | Nifty IT |
$1.5$ | $0.1060$ | $0.4584$ | $0.4356$ | $3.2423$
$1.7$ | $0.0974$ | $0.4362$ | $0.4664$ | $3.3924$
$1.9$ | $0.0888$ | $0.4140$ | $0.4972$ | $3.5585$
$2.1$ | $0.0803$ | $0.3917$ | $0.5280$ | $3.7405$
$2.3$ | $0.0717$ | $0.3695$ | $0.5588$ | $3.9384$
$2.5$ | $0.0631$ | $0.3473$ | $0.5896$ | $4.1523$
$2.7$ | $0.0546$ | $0.3250$ | $0.6204$ | $4.3821$
$2.9$ | $0.0460$ | $0.3028$ | $0.6512$ | $4.6279$
$3.1$ | $0.0374$ | $0.2806$ | $0.6820$ | $4.8896$
$3.3$ | $0.0289$ | $0.2583$ | $0.7128$ | $5.1673$
$3.5$ | $0.0203$ | $0.2361$ | $0.7436$ | $5.4609$

Using these results, we plot the optimal portfolio allocation graph and
efficient frontier, which are illustrated in Fig. 8 and Fig. 9 respectively.

Figure 8: Portfolio Allocation for the Piecewise Quadratic Function Based
Approximation Figure 9: Efficient Frontier of the Portfolio for the Piecewise
Quadratic Generating Function Based Approximation

##  6 Conclusion

In this paper, we derive safe convex approximations for the ambiguous chance
constraint based uncertain portfolio optimization problems. Our study
comprises two families of perturbation distributions– (i) family with known
upper $\&$ lower mean bounds, and (ii) family with known upper $\&$ lower mean
bounds and standard deviations. We use a piecewise linear function and an
exponential function as the generator to derive the safe approximation for the
former case, whereas a piecewise quadratic function is used for the latter
case. The approximations which we get from family-i are quadratic programming
problems and the approximation for family-ii is a quadratically constrained
quadratic programming problem. Thus these approximations are computationally
tractable and by solving them we get the robust solutions to our original
ambiguous chance constrained problems. Moreover, we solve a numerical problem
by the use of these approximations and we plot the optimal allocations and
efficient frontiers for each case. From the two approximations of family-ii,
it is observed that the optimal portfolio risks for the Bernstein
approximation are higher than that of the piecewise linear generator based
approximations. Furthermore, our study can be used for any optimization
problem with the linear chance constraint(s). In this study, we use piecewise
linear and quadratic functions as the generators for deriving the safe
approximations. For the derivation of piecewise linear generator based
approximations, we require the information regarding the mean bounds of the
perturbations, whereas for piecewise quadratic generator based approximations
we need the standard deviation (or variance) information of the perturbations
along with their mean bounds. Similarly, if we have the information regarding
the first $n$ moment statistics of the perturbations, then we can extend our
study to use piecewise generating functions of any degree $n$, where $n$ is a
positive integer. Data Availability: The Nifty price data of the three
sectors– Nifty Bank, Nifty Infra and Nifty IT used in the study are collected
for the time period June 2017 to May 2022 and the data is available in the
website https://finance.yahoo.com. We also attach the data in the file named
Nifty price data Conflict of Interest Statement: The authors declare that they
have no conflicts of interest.

## References

  * (1) Markowitz, H., Portfolio Selection, Journal of Finance, Vol. 7, pp. 77-91 (1952).
  * (2) Markowitz, H., Portfolio Selection: efficient diversification of investments, Basil Blackwell, New York (1959).
  * (3) Estrada, J. Systematic risk in emerging markets: the D-CAPM. Emerging Markets Review, 3(4), pp. 365-379 (2002).
  * (4) Estrada, J. Mean-semivariance behavior: Downside risk and capital asset pricing. International Review of Economics & Finance, 16(2), pp. 169-185 (2007).
  * (5) El Ghaoui, L., and Lebret, H. Robust solutions to least-squares problems with uncertain data. SIAM Journal on matrix analysis and applications, 18(4), pp. 1035-1064 (1997).
  * (6) El Ghaoui, L., Oustry, F., and Lebret, H. Robust solutions to uncertain semidefinite programs. SIAM Journal on Optimization, 9(1), pp. 33-52 (1998).
  * (7) Ben-Tal, A., and Nemirovski, A. Robust solutions of uncertain linear programs. Operations research letters, 25(1), pp. 1-13 (1999).
  * (8) Ben-Tal, A., and Nemirovski, A. Robust solutions of linear programming problems contaminated with uncertain data. Mathematical programming, 88(3), pp. 411-424 (2000).
  * (9) Goldfarb, D., and Iyengar, G., Robust portfolio selection problems. Mathematics of operations research, 28(1), pp. 1-38 (2003).
  * (10) Tütüncü, R. H., and Koenig, M., Robust asset allocation. Annals of Operations Research, 132(1-4), pp. 157-187 (2004).
  * (11) Fliege, J., and Werner, R. Robust multiobjective optimization & applications in portfolio optimization. European Journal of Operational Research, 234(2), pp. 422-433 (2014).
  * (12) Sehgal, R. and Mehra, A. Robust reward-risk ratio portfolio optimization. International Transactions in Operational Research (2019).
  * (13) Biswal, M. P., Biswal, N. P., and Li, D. (1998). Probabilistic linear programming problems with exponential random variables: A technical note. European Journal of Operational Research, 111(3), 589-597.
  * (14) Calafiore, Giuseppe Carlo, and L. El Ghaoui, ”On distributionally robust chance-constrained linear programs.” Journal of Optimization Theory and Applications 130, no. 1, pp. 1-22 (2006).
  * (15) Erdoğan, E., and Iyengar, G., Ambiguous chance constrained problems and robust optimization. Mathematical Programming, 107(1), pp. 37-61 (2006).
  * (16) Nemirovski, A., and Shapiro, A., Convex approximations of chance constrained programs. SIAM Journal on Optimization, 17(4), pp. 969-996 (2007).
  * (17) Yanıkoğlu, İ., and den Hertog, D. Safe approximations of ambiguous chance constraints using historical data. INFORMS Journal on Computing, 25(4), pp. 666-681 (2013).
  * (18) Bertsimas, D., Gupta, V., and Kallus, N., Data-driven robust optimization. Mathematical Programming, 167, pp. 235-292 (2018).
  * (19) Zhang, Y., Jiang, R., and Shen, S. Ambiguous chance-constrained binary programs under mean-covariance information. SIAM Journal on Optimization, 28(4), pp. 2922-2944 (2018).
  * (20) Bertsimas, D., Brown, D. B., and Caramanis, C. Theory and applications of robust optimization. SIAM review, 53(3), pp. 464-501 (2011).
  * (21) Swain, P., & Ojha, A. K., Bi-level optimization approach for robust mean-variance problems. RAIRO-Operations Research, 55(5), pp. 2941-2961 (2021).
  * (22) Ben-Tal, A., El Ghaoui, L., and Nemirovski, A. Robust optimization (Vol. 28). Princeton University Press (2009).