Monte carlo within simulated annealing for integral constrained optimizations

Abstract

For years, Value-at-Risk and Expected Shortfall have been well established measures of market risk and the Basel Committee on Banking Supervision recommends their use when controlling risk. But their computations might be intractable if we do not rely on simplifying assumptions, in particular on distributions of returns. One of the difficulties is linked to the need for Integral Constrained Optimizations. In this article, two new stochastic optimization-based Simulated Annealing algorithms are proposed for addressing problems associated with the use of statistical methods that rely on extremizing a non-necessarily differentiable criterion function, therefore facing the problem of the computation of a non-analytically reducible integral constraint. We first provide an illustrative example when maximizing an integral constrained likelihood for the stress-strength reliability that confirms the effectiveness of the algorithms. Our results indicate no clear difference in convergence, but we favor the use of the problem approximation strategy styled algorithm as it is less expensive in terms of computing time. Second, we run a classical financial problem such as portfolio optimization, showing the potential of our proposed methods in financial applications.

Appendices

Appendix A proof to theorem 1: ergodicity of MCWMSA

The proof for Theorem 1 follows standard arguments for the convergence of Inhomogeneous Markov Chains (e.g., Isaacson and Madsen (1976) and Seneta (1981) for related results). More precisely, we first form a lower bound on the probability of reaching any solution from any local minimum, and then show that this bound does not approach zero too quickly. Let G be the generation probability that satisfies 19.

• a. Let:

  $$\begin{aligned} \displaystyle \varDelta _{G} = \min _{\begin{array}{c} (x',\lambda ')\in \mathcal {N}({x,\lambda }),\\ {x,\lambda }\in \mathcal {S} \end{array}} G(x,\lambda ;x',\lambda ') \end{aligned}$$

  and:

  $$\begin{aligned} \displaystyle \varDelta _{L^{M}} = 2\max _{\begin{array}{c} (x',\lambda ')\in \mathcal {N}(x,\lambda ),\\ (x'\lambda ')\in \mathcal {S},y,y'\in \mathcal {Y}^M \end{array}}\{|L^{M}(x',\lambda ') -L^{M}(x,\lambda )|\}, \end{aligned}$$

  be, respectively, the value of the smallest one-step solution generation probabilities and the size of the largest potential jump in the value of the auxiliary penalty function between any pair of trial solutions. $\mathcal {S}=\mathcal {X}\times \mathbb {R}_{+}^{r}$ denotes the search space while $\mathcal {N}(x,\lambda )$ represents the neighborhood of $(x,\lambda )$. If $(x,\lambda )\ne (x',\lambda ')$, then for all $(x,\lambda )\in \mathcal {S}$ and $(x',\lambda ')\in \mathcal {N}(x,\lambda )$, we have:

  $$\begin{aligned} \begin{aligned}&P_{MCWMSA}(x,\lambda ;x',\lambda ') = \int _{\mathcal {Y}^{M}}\int _{\mathcal {Y}^{M}} G(x,\lambda ;x',\lambda ')A^{M}_{T_{k}}(x,\lambda ;x',\lambda ')q^{M}_{x}(y)q^{M}_{x'}(y')dydy' \\&\quad = G(x,\lambda ;x',\lambda ')\int _{\mathcal {Y}^{M}}\int _{\mathcal {Y}^{M}} A^{M}_{T_{k}}(x,\lambda ;x',\lambda ')q^{M}_{x}(y)q^{M}_{x'}(y')dydy' \\&\quad = {\left\{ \begin{array}{ll} G(x,\lambda ;x',\lambda ')\int _{\mathcal {Y}^{M}}\int _{\mathcal {Y}^{M}} \min \left\{ \exp \left( -\frac{(L^{M}(x',\lambda ')-L^{M}(x,\lambda ))}{T}\right) ,1\right\} q^{M}_{x}(y)q^{M}_{x'}(y')dydy', &{} {\text {if }u=1 }\\ G(x,\lambda ;x',\lambda ')\int _{\mathcal {Y}^{M}}\int _{\mathcal {Y}^{M}} \min \left\{ \exp \left( -\frac{(L^{M}(x,\lambda )-L^{M}(x',\lambda '))}{T}\right) ,1\right\} q^{M}_{x}(y)q^{M}_{x'}(y')dydy' , &{} {\text {if }u=0 } \end{array}\right. }\\&\quad \ge G(x,\lambda ;x',\lambda ')\exp \left( -\frac{\varDelta _{L^{M}}}{T_{k}}\right) \int _{\mathcal {Y}^{M}}\int _{\mathcal {Y}^{M}} q^{M}_{x}(y)q^{M}_{x'}(y')dydy'\\ \\&\quad \ge \varDelta _{G}\exp \left( -\frac{\varDelta _{L^{M}}}{T_{k}}\right) . \end{aligned} \end{aligned}$$

  (A.1)

• b. For any $\lambda$, let:

  $$\begin{aligned} \displaystyle \varDelta _{min}^{M} = \min _{\begin{array}{c} (x',\lambda ')\in \mathcal {N}(x,\lambda ),\\ A^{M}_{T_{k}}<1, y,y'\in \mathcal {Y}^{M} \end{array}} \left( L^{M}(x',\lambda ') - L^{M}(x,\lambda ) \right) \ge 0, \end{aligned}$$

  be the size of the smallest positive jump in the value of the objective function between any solution pair. If $(x',\lambda ')=(x,\lambda )$, we have:

  $$\begin{aligned}&P_{MCWMSA}(x,\lambda ;x',\lambda ') \nonumber \\&\quad = 1-\int _{\mathcal {N}(x,\lambda )}\int _{\mathcal {Y}^{M}}\int _{\mathcal {Y}^{M}} A_{T_{k}}^{M}(x,\lambda ;x',\lambda ') q^{M}_{x}(y)q^{M}_{x'}(y')dydy'G(x,\lambda ;dx',d\lambda ') \nonumber \\&\quad = \int _{\begin{array}{c} \mathcal {N}(x,\lambda ),A^{M}_{T_{k}}<1 \end{array}}\int _{\mathcal {Y}^{M}}\int _{\mathcal {Y}^{M}} \left( 1- A_{T_{k}}^{M}(x,\lambda ;x',\lambda ')\right) q^{M}_{x}(y)q^{M}_{x'}(y')dydy'G(x,\lambda ;dx',d\lambda ') \nonumber \\&\quad \ge \int _{\begin{array}{c} \mathcal {N}(x,\lambda ),A^{M}_{T_{k}}<1 \end{array}}\int _{\mathcal {Y}^{M}}\int _{\mathcal {Y}^{M}} \varDelta _{G}(1 - \exp (-\varDelta _{min}^{M}/T_{k}))q^{M}_{x}(y)q^{M}_{x'}(y')dydy'dx'd\lambda '. \end{aligned}$$

  (A.2)

  Since, $T_{k}$ is a decreasing sequence, it is always possible to find a $k_{0}>0$ so that for all $k\ge k_{0}$, $\displaystyle (1-\exp (-\varDelta _{min}^{M}/T_{k}))\ge \exp (-\varDelta _{L^{M}}/T_{k})$.

  Thus:

  $$\begin{aligned} P_{MCWMSA}(x,\lambda ;x,\lambda )\ge \exp (-\varDelta _{L^{M}}/T_{k}). \end{aligned}$$

• c. Based on the result obtained above, for all $(x,\lambda ),(x',\lambda ')\in \mathcal {S}$ and $k\ge k_{0}$, the N-step transition probability from $(x,\lambda )=(x_{0},\lambda _{0})$ to $(x',\lambda ')=(x_{N},\lambda _{N})$ satisfies the following:

  $$\begin{aligned} P_{MCWMSA}^{N}(x,\lambda ;x',\lambda ')= \prod _{i=0}^{N-1}P_{MCWMSA}(x_{i},\lambda _{i};x_{i+1},\lambda _{i+1}) \ge \left( \varDelta _{G}\exp (-\varDelta _{L^{M}}/T_{k}) \right) ^{N}. \end{aligned}$$

  Let $\displaystyle {\tau _{1}\left( P\right) }$ denote the coefficient of ergodicity for the transition probability matrix P, defined by:

  $$\begin{aligned} \begin{aligned}&\tau _{1}\left( P_{MCWMSA}^{N}\right) = 1 - \min _{(x,\lambda ),(x',\lambda ')\in \mathcal {S}}\sum _{(x'',\lambda '')\in \mathcal {S}} d(x,\lambda ,x',\lambda ',x'',\lambda '') \end{aligned} \end{aligned}$$

  where $d(x,\lambda ,x',\lambda ',x'',\lambda '')\!=\!\min \left\{ P_{MCWMSA}^{N}(x,\lambda ; x'',\lambda ''), P_{MCWMSA}^{N}(x',\lambda '; x'',\lambda '' ) \right\}$ An Inhomogeneous Markov Chain is weakly ergodic if and only if there is a strictly increasing sequence of positive numbers so that:

  $$\begin{aligned} \sum _{k=0}\left( 1-\tau _{1}(P_{MCWMSA}^{N})\right) = \infty , \end{aligned}$$

  (A.3)

  using the bound on $P_{MCWMSA}^{N}(x,\lambda ;x',\lambda ')$, we find that:

  $$\begin{aligned}&1- \tau _{1}(P_{MCWMSA}^{N})=\min _{(x,\lambda ),(x',\lambda ')\in \mathcal {S}}\sum _{(x'',\lambda '')\in \mathcal {S}}d(x,\lambda ,x',\lambda ',x'',\lambda '') \nonumber \\&\quad \ge \min _{(x,\lambda ),(x',\lambda ')\in \mathcal {S}}\min _{(x'',\lambda '')\in \mathcal {S}}\left\{ P_{MCWMSA}^{N}(x,\lambda ; x'',\lambda ''), P_{MCWMSA}^{N}(x',\lambda '; x'',\lambda '' ) \right\} \nonumber \\&\quad \ge \left( \varDelta _{G}\exp (-\varDelta _{L^{M}}/T_{k}) \right) ^{N} = \varDelta _{G}^{N}\exp (-N\varDelta _{L^{M}}/T_{k}). \end{aligned}$$

  (A.4)

  Substituting $\displaystyle {T_{k}\ge \frac{N\varDelta _{L^{M}}}{\log {(k+1)}}, ~~k\ge 0 }$ in Eq. A.3 gives:

  $$\begin{aligned} \begin{aligned} \sum _{k=0}\left( 1-\tau _{1}(P_{MCWMSA}^{N})\right)&\ge \sum _{k=k_{0}}\varDelta _{G}^{N}\exp (-N\varDelta _{L^{M}}/T_{k})\\&\ge \varDelta _{G}^{N}\sum _{k=k_{0}}\dfrac{1}{k+1} = \infty . \end{aligned} \end{aligned}$$

  (A.5)

  This implies that the Markov Chain, $(x,\lambda ,x',\lambda ')$, is weakly ergodic.

• d. Furthermore, since the acceptance probability function $A_{T_{K}}^{M}$ is an exponential rational function in $1/T_{k}$ and belongs to the closed class of asymptotically monotone functions (Anily and Federgruen (1987), Proposition 1) then the sequence of $(x_{n},\lambda _{n})$ converges in distribution (Anily and Federgruen (1987), Corollary 1). This implies that the Markov Chain is strongly ergodic.$\square$

Appendix B robustness checks

This section provides the robustness checks of the results in Tables 2 and 3 with respect the main parameters of the algorithms, namely a, N and $T^{*}$ (boldface values), showing similar results in terms of efficiency, computational time and GIMHSA dominance.

Table 7 Parameter Estimates for Dataset#1, comparing CSA of Smith et al. (2015), MCWMSA and GIMHSA Algorithms when solving the Integral Constrained Log-likelihood in Eq. 23 using Dataset#1 as reported in Table 1, with $a=\mathbf {0.80}$, $N=\mathbf {250}$, $T_{0}=10$ and $T^{*}= 10^{-6}$

Table 8 Reported Characteristics of the Methods in terms of Performances of the Algorithms on the Dataset#1 of Smith et al. (2015) with $a=\mathbf {0.80}$, $N=\mathbf {250}$, $T_{0}=10$ and $T^{*}= 10^{-6}$

Table 9 Parameter Estimates for Dataset#1, comparing CSA of Smith et al. (2015), MCWMSA and GIMHSA Algorithms when solving the Integral Constrained Log-likelihood in Eq. 23 using Dataset#1 as reported in Table 1, with $a=0.95$, $N=\mathbf {250}$, $T_{0}=10$ and $T^{*}= \mathbf {10^{-4}}$

Table 10 Reported Characteristics of the Methods in terms of Performances of the Algorithms on the Dataset#1 of Smith et al. (2015) with $a=0.95$, $N=\mathbf {500}$, $T_{0}=10$ and $T^{*}= \mathbf {10^{-4}}$

Table 11 Parameter Estimates for Dataset#1, comparing CSA of Smith et al. (2015), MCWMSA and GIMHSA Algorithms when solving the Integral Constrained Log-likelihood in Eq. 23 using Dataset#1 as reported in Table 1, with $a=0.95$, $N=\mathbf {250}$, $T_{0}=10$ and $T^{*}= \mathbf {10^{-4}}$

Table 12 Reported Characteristics of the Methods in terms of Performances of the Algorithms on the Dataset#1 of Smith et al. (2015) for $a=0.95$, $N=\mathbf {250}$, $T_{0}=10$ and $T^{*}= \mathbf {10^{-4}}$

Appendix C Pseudo-codes