Robust_Portfolio.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bc31e2d0",
   "metadata": {},
   "source": [
    "# Robust Portfolio Optimization\n",
    "We use the returns of six portfolios formed on size and book-to-market (2$\\times$3) obtained from Kenneth French's website. \n",
    "We choose monthly returns from 1984 January to 2014 January. \n",
    "This gives us a total of $N=360$ return realizations for each of the six portfolios.\n",
    "\n",
    "Let $\\mathbf{a}\\in \\mathbb{R}^6$ be an investment strategy on the six portfolios, subject to the constraints $a_i\\geq 0, \\sum^6_{i=1}a_i=1$. Let $\\mathbf{R}$ be a $360\\times 6$ matrix with $\\mathbf{r}_1,\\ldots,\\mathbf{r}_{360}\\in \\mathbb{R}^6$ the vectors of realisations of returns for all six portfolios. For each realisation $\\mathbf{r}_i$, the wealth level after a single period is defined as $1+\\mathbf{a}^T\\mathbf{r}_i$. We consider the following worst-case risk evaluation minimization problem \n",
    "\\begin{align}\n",
    "\\begin{split}\n",
    "\\min_{\\mathbf{a}\\in \\mathbb{R}^{6},\\mathbf{a}\\geq \\mathbf{0}, c\\in \\mathbb{R}}\\left\\{c~\\middle|~\\sup_{(\\mathbf{q},\\bar{\\mathbf{q}})\\in U_{\\phi,h}(\\mathbf{p})}-\\sum^{360}_{i=1}\\bar{q}_iu(1+\\mathbf{a}^T\\mathbf{r}_i)\\leq c,\\sum^6_{i=1}a_i=1\\right\\},   \n",
    "\\end{split}\n",
    "\\end{align}\n",
    "where we choose the modified chi-squared function $\\phi(x)=(x-1)^2$ with radius $r=\\frac{1}{N}\\chi^2_{359,0.95}$. Furthermore, we choose $h(p)=1-(1-p)^2,~p\\in[0,1]$ as the distortion function and $u(x)=1-e^{-x/\\lambda}$ with $\\lambda=10$ the exponential utility function.\n",
    "The solution of this problem will be compared to the nominal risk evaluation minimization problem, where the nominal probability is the uniform distirbution $p_i=\\frac{1}{360}$.\n",
    "\\begin{align}\n",
    "\\begin{split}\n",
    "\\min_{\\mathbf{a}\\in \\mathbb{R}^6,\\mathbf{a}\\geq \\mathbf{0}, c\\in \\mathbb{R}}\\left\\{c~\\middle|~\\sup_{\\bar{\\mathbf{q}}\\in M_h(\\mathbf{p})}-\\sum^{360}_{i=1}\\bar{q}_iu(1+\\mathbf{a}^T\\mathbf{r}_i)\\leq c, \\sum^6_{i=1}a_i=1,\\right\\}.    \n",
    "\\end{split}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "7cb9fe4b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import cvxpy as cp\n",
    "import mosek\n",
    "import matplotlib.pyplot as plt\n",
    "import datetime as date\n",
    "from datetime import datetime as dt\n",
    "from dateutil.relativedelta import *\n",
    "import distortion_function as hf\n",
    "import phi_divergence as phi\n",
    "import affine_approx as af\n",
    "import Utility_functions as ut\n",
    "import Cutting_plane as ct\n",
    "from time import process_time\n",
    "import importlib\n",
    "import scipy.stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "30a759b2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<module 'Cutting_plane' from 'C:\\\\Users\\\\gjin\\\\Robust\\\\Cutting_plane.py'>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "importlib.reload(ct)    #### reload the module when changes to the module are made"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "6d2c4eb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "############# This code checks if you have installed Mosek with a valid academic license\n",
    "\n",
    "import mosek\n",
    "\n",
    "# Try to initialize MOSEK\n",
    "with mosek.Env() as env:\n",
    "    pass\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "36fae8a8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Date</th>\n",
       "      <th>SMALL LoBM</th>\n",
       "      <th>ME1 BM2</th>\n",
       "      <th>SMALL HiBM</th>\n",
       "      <th>BIG LoBM</th>\n",
       "      <th>ME2 BM2</th>\n",
       "      <th>BIG HiBM</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1926-07-01</td>\n",
       "      <td>1.0874</td>\n",
       "      <td>0.9349</td>\n",
       "      <td>-0.0695</td>\n",
       "      <td>5.7168</td>\n",
       "      <td>1.9620</td>\n",
       "      <td>1.4222</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1926-08-01</td>\n",
       "      <td>0.7030</td>\n",
       "      <td>1.2300</td>\n",
       "      <td>5.3842</td>\n",
       "      <td>2.7154</td>\n",
       "      <td>2.6930</td>\n",
       "      <td>6.3154</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1926-09-01</td>\n",
       "      <td>-2.9117</td>\n",
       "      <td>-0.1303</td>\n",
       "      <td>-0.4374</td>\n",
       "      <td>1.4287</td>\n",
       "      <td>0.0704</td>\n",
       "      <td>-0.7967</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1926-10-01</td>\n",
       "      <td>-3.8196</td>\n",
       "      <td>-4.5860</td>\n",
       "      <td>-2.0112</td>\n",
       "      <td>-3.5898</td>\n",
       "      <td>-2.3398</td>\n",
       "      <td>-4.0970</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1926-11-01</td>\n",
       "      <td>3.1806</td>\n",
       "      <td>3.7233</td>\n",
       "      <td>2.0944</td>\n",
       "      <td>3.1292</td>\n",
       "      <td>2.8952</td>\n",
       "      <td>3.4614</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1139</th>\n",
       "      <td>2021-06-01</td>\n",
       "      <td>5.6058</td>\n",
       "      <td>0.4400</td>\n",
       "      <td>-1.0979</td>\n",
       "      <td>4.8188</td>\n",
       "      <td>-1.2594</td>\n",
       "      <td>-4.0036</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1140</th>\n",
       "      <td>2021-07-01</td>\n",
       "      <td>-5.5593</td>\n",
       "      <td>-1.8623</td>\n",
       "      <td>-3.6521</td>\n",
       "      <td>3.1048</td>\n",
       "      <td>-0.0099</td>\n",
       "      <td>-2.3000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1141</th>\n",
       "      <td>2021-08-01</td>\n",
       "      <td>2.3903</td>\n",
       "      <td>1.5124</td>\n",
       "      <td>2.6680</td>\n",
       "      <td>3.5667</td>\n",
       "      <td>1.4122</td>\n",
       "      <td>3.0371</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1142</th>\n",
       "      <td>2021-09-01</td>\n",
       "      <td>-4.3421</td>\n",
       "      <td>-3.4661</td>\n",
       "      <td>0.6445</td>\n",
       "      <td>-5.4525</td>\n",
       "      <td>-3.8570</td>\n",
       "      <td>-0.2526</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1143</th>\n",
       "      <td>2021-10-01</td>\n",
       "      <td>2.9626</td>\n",
       "      <td>3.7688</td>\n",
       "      <td>3.8617</td>\n",
       "      <td>8.0441</td>\n",
       "      <td>3.1530</td>\n",
       "      <td>6.2389</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1144 rows × 7 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "           Date  SMALL LoBM  ME1 BM2  SMALL HiBM  BIG LoBM  ME2 BM2  BIG HiBM\n",
       "0    1926-07-01      1.0874   0.9349     -0.0695    5.7168   1.9620    1.4222\n",
       "1    1926-08-01      0.7030   1.2300      5.3842    2.7154   2.6930    6.3154\n",
       "2    1926-09-01     -2.9117  -0.1303     -0.4374    1.4287   0.0704   -0.7967\n",
       "3    1926-10-01     -3.8196  -4.5860     -2.0112   -3.5898  -2.3398   -4.0970\n",
       "4    1926-11-01      3.1806   3.7233      2.0944    3.1292   2.8952    3.4614\n",
       "...         ...         ...      ...         ...       ...      ...       ...\n",
       "1139 2021-06-01      5.6058   0.4400     -1.0979    4.8188  -1.2594   -4.0036\n",
       "1140 2021-07-01     -5.5593  -1.8623     -3.6521    3.1048  -0.0099   -2.3000\n",
       "1141 2021-08-01      2.3903   1.5124      2.6680    3.5667   1.4122    3.0371\n",
       "1142 2021-09-01     -4.3421  -3.4661      0.6445   -5.4525  -3.8570   -0.2526\n",
       "1143 2021-10-01      2.9626   3.7688      3.8617    8.0441   3.1530    6.2389\n",
       "\n",
       "[1144 rows x 7 columns]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "### Reading the Kenneth-French return data's\n",
    "df_returns6 = pd.read_csv('6_Portfolios_2x3.csv', skiprows = 15) \n",
    "df_returns = df_returns6[0:1144].copy()\n",
    "df_returns['Date'] = pd.to_datetime(df_returns['Date'], format = '%Y%m')\n",
    "for i in range(1, len(df_returns.columns)):\n",
    "    df_returns[df_returns.columns[i]] = pd.to_numeric(df_returns[df_returns.columns[i]])\n",
    "df_returns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "3c7a79e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Selecting the return values between 1984 and 2014, since everything was in percentage, we divde the return by 100\n",
    "startdate = dt(1984,1,3)\n",
    "enddate = dt(2014,1,3)\n",
    "X = df_returns[np.logical_and(df_returns.Date >= startdate, df_returns.Date <= enddate )][df_returns.columns[1:7]]\n",
    "X = X.reset_index(drop = True)\n",
    "R = X.to_numpy()     \n",
    "R = R/100\n",
    "W0 = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "a823a7ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "phi_func = phi.mod_chi2_cut\n",
    "phi_conj = phi.mod_chi2_conj\n",
    "phi_dot = 2\n",
    "\n",
    "h_conj = hf.h_spw_conj\n",
    "h_func = hf.h_spw_cut\n",
    "h_eva = hf.h_spw_eva\n",
    "### parameter of h function, in this case we have h(p)=1-(1-p)^2, hence we take par = 2\n",
    "h_par = 2 \n",
    "\n",
    "utility = ut.exp_utility\n",
    "utility_eva = ut.exp_utility_eva\n",
    "\n",
    "### Defining the nominal probability vector and denote I the number of assets\n",
    "N=R.shape[0]\n",
    "p = np.zeros(N)+1/N\n",
    "I = R.shape[1]\n",
    "\n",
    "\n",
    "### defining phi-divergence uncertainty set radius\n",
    "alpha_phi_set = 0.95\n",
    "r = phi_dot/(2*N)*scipy.stats.chi2.ppf(alpha_phi_set, N-1)   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "d4926095",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lb: -0.08786117865952753 ub: [-0.09629894] iter: 0\n",
      "lb: -0.08919039635560265 ub: [-0.09044157] iter: 1\n",
      "lb: -0.08940651729343968 ub: [-0.09009738] iter: 2\n",
      "lb: -0.08948381007093242 ub: [-0.08966657] iter: 3\n",
      "lb: -0.08948143999397146 ub: [-0.08952706] iter: 4\n",
      "time 110.0625 seconds\n",
      "('sol:', array([3.54055574e-11, 1.71782966e-10, 1.15974067e-10, 2.50944809e-01,\n",
      "       7.49055190e-01, 3.12074820e-10]), 'lb:', -0.0895270557704481, 'ub:', -0.08948143999397146, 'iter:', 4)\n"
     ]
    }
   ],
   "source": [
    "###### Cutting-plane robust risk minimization solution\n",
    "\n",
    "e_tol = 0.0001\n",
    "tb = process_time()\n",
    "res_rob = ct.cut_rob_pmin(R,p,e_tol,utility,utility_eva,h_func,phi_func,h_eva,r,W0, par=2, par_u=10)  \n",
    "te = process_time()\n",
    "print('time', te-tb, 'seconds')\n",
    "print(res_rob)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "7ab37c66",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lb: -0.0937186092796091 ub: [-0.09629894] iter: 0\n",
      "lb: -0.09393146527231497 ub: [-0.09439286] iter: 1\n",
      "lb: -0.09387991226421565 ub: [-0.09414542] iter: 2\n",
      "lb: -0.09397411716061542 ub: [-0.09408978] iter: 3\n",
      "lb: -0.09396835079972576 ub: [-0.0940359] iter: 4\n",
      "time 61.4375 seconds\n",
      "('sol:', array([6.12265646e-10, 1.40642318e-07, 1.43623410e-01, 3.84557777e-01,\n",
      "       4.71818668e-01, 3.77053986e-09]), 'lb:', -0.09403589599270797, 'ub:', -0.09396835079972576, 'iter:', 4)\n"
     ]
    }
   ],
   "source": [
    "###### Cutting-plane nominal risk minimization solution\n",
    "e_tol = 0.0001\n",
    "tb = process_time()\n",
    "res_nom = ct.cut_nom_pmin(R,p,e_tol,utility,utility_eva,h_eva,W0,par=h_par,par_u=10) \n",
    "te = process_time()\n",
    "print('time', te-tb, 'seconds')\n",
    "print(res_nom)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "8a9909ea",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.08937877108425232\n"
     ]
    }
   ],
   "source": [
    "#### Calculating the worst-case risk of the nominal solution from the cutting-plane\n",
    "w_nom = res_nom[1]       \n",
    "x_nom = ut.exp_utility_eva(R,w_nom,W0,10)\n",
    "print(ct.robustcheck(x_nom,p,h_func,phi_func,r,h_par)[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "e3aa6dfe",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.0939484452866081\n"
     ]
    }
   ],
   "source": [
    "### Calculating the nominal risk of the robust solution from the cutting-plane\n",
    "w_rob = res_rob[1]\n",
    "x_rob = ut.exp_utility_eva(R,w_rob,W0,10)\n",
    "[risk_robnom, qb_nomrb]=ct.nominal_risk(x_rob,p,h_eva,h_par) \n",
    "print(risk_robnom)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "c448fc9e",
   "metadata": {},
   "outputs": [],
   "source": [
    "### The precision of the PL approximation\n",
    "e_pl = 0.001       \n",
    "\n",
    "### This gives the slope and constants that constitutes the PL approximation of h(p)=1-(1-p)^2\n",
    "x_points = af.affine_approx_hspw(h_par,e_pl)     \n",
    "[slope, const] = af.makepoints(af.sing_pw,x_points,par = h_par)  \n",
    "                                                               "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "6f2ceaf2",
   "metadata": {},
   "outputs": [],
   "source": [
    "### The following function solves the robust problem using the picewise-linear approximation method\n",
    "\n",
    "def af_rob_exp_pmin(p,R,r,phi_conj,slope,const,W0,par = 2, par_u = 1):\n",
    "    N = len(p)\n",
    "    I = len(R[0])\n",
    "    K = len(slope)\n",
    "    c = cp.Variable(1)\n",
    "    lbda = cp.Variable((N,K), nonneg = True)\n",
    "    a = cp.Variable(I, nonneg = True)\n",
    "    v = cp.Variable(K, nonneg = True)\n",
    "    alpha = cp.Variable(1)\n",
    "    beta = cp.Variable(1)\n",
    "    gamma = cp.Variable(1,nonneg = True)\n",
    "    t = cp.Variable(N)\n",
    "    s = cp.Variable(N)\n",
    "    w = cp.Variable(N)\n",
    "    constraints = [cp.abs(a)<= W0, cp.sum(a) == 1]\n",
    "    for i in range(N):\n",
    "        arg = -W0*(1+(R @ a)[i])/par_u\n",
    "        constraints.append(-(1-cp.exp(arg)) - cp.sum(lbda[i]) - beta <= 0)\n",
    "        constraints.append(s[i] == -alpha + lbda[i]@slope)\n",
    "        constraints.append(lbda[i] <= v)\n",
    "        constraints = phi_conj(gamma,s[i],t[i],w[i],constraints)\n",
    "    constraints.append(alpha + beta + gamma * r  + v@const + p@t <= c)\n",
    "    obj = cp.Minimize(c)\n",
    "    prob = cp.Problem(obj,constraints)\n",
    "    prob.solve(solver=cp.MOSEK)\n",
    "    return(a.value, prob.value)\n",
    "\n",
    "\n",
    "\n",
    "### The following function solves the nominal problem using the picewise-linear approximation method\n",
    "\n",
    "def af_nom_exp_pmin(p,R,slope,const,W0,par = 2, par_u = 1):\n",
    "    N = len(p)\n",
    "    I = len(R[0])\n",
    "    K = len(slope)\n",
    "    c = cp.Variable(1)\n",
    "    lbda = cp.Variable((N,K), nonneg = True)\n",
    "    a = cp.Variable(I, nonneg = True)\n",
    "    v = cp.Variable(K, nonneg = True)\n",
    "    beta = cp.Variable(1)\n",
    "    constraints = [cp.abs(a)<= W0, cp.sum(a) == 1]\n",
    "    for i in range(N):\n",
    "        arg = -W0*(1+(R @ a)[i])/par_u\n",
    "        constraints.append(-(1-cp.exp(arg)) - cp.sum(lbda[i]) - beta <= 0)\n",
    "        constraints.append(lbda[i] <= v)\n",
    "    constraints.append(beta  + v@const + p@lbda@slope <= c)\n",
    "    obj = cp.Minimize(c)\n",
    "    prob = cp.Problem(obj,constraints)\n",
    "    prob.solve(solver=cp.MOSEK)\n",
    "    return(a.value, prob.value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "eca4bc8c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sol: [7.13420613e-10 3.78209461e-09 2.51930935e-09 3.79717189e-01\n",
      " 6.20282798e-01 6.19859795e-09] lb: -0.08951029319587173\n",
      "time 12.671875 seconds\n",
      "sol: [1.87609197e-09 9.98318381e-09 1.30755000e-08 3.76560089e-01\n",
      " 6.23439871e-01 1.45642906e-08] ub: -0.08948262687815382\n",
      "time 12.46875 seconds\n"
     ]
    }
   ],
   "source": [
    "t1_af_res = process_time()\n",
    "af_res = af_rob_exp_pmin(p,R,r,phi_conj,slope,const,W0,par = h_par, par_u = 10)\n",
    "t2_af_res = process_time()\n",
    "### the solution and objective value of the piecewise-linear approximation method\n",
    "print('sol:',af_res[0], 'lb:', af_res[1])   \n",
    "### Total run time of the algorithm\n",
    "print('time', t2_af_res- t1_af_res, 'seconds')\n",
    "\n",
    "\n",
    "t1_af_res_ub = process_time()\n",
    "af_res_ub = af_rob_exp_pmin(p,R,r,phi_conj,slope,const+e_pl,W0,par = h_par, par_u = 10)\n",
    "t2_af_res_ub = process_time()\n",
    "\n",
    "print('sol:',af_res_ub[0], 'ub:', af_res_ub[1])   \n",
    "### Total run time of the algorithm\n",
    "print('time', t2_af_res_ub- t1_af_res_ub, 'seconds')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "2be58821",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sol: [1.00618210e-09 8.59290405e-02 2.07493483e-01 1.30601768e-01\n",
      " 5.75975698e-01 8.89528303e-09] lb: -0.09400888974788543\n",
      "time 5.984375 seconds\n",
      "sol: [1.50352434e-09 7.84892167e-02 2.07079630e-01 1.29813633e-01\n",
      " 5.84617505e-01 1.41332539e-08] ub: -0.09397821096000436\n",
      "time 6.4375 seconds\n"
     ]
    }
   ],
   "source": [
    "t1_af_nom = process_time()\n",
    "af_nom = af_nom_exp_pmin(p,R,slope,const,1,par = h_par, par_u = 10)\n",
    "t2_af_nom = process_time()\n",
    "### the solution and objective value of the piecewise-linear approximation method\n",
    "print('sol:',af_nom[0], 'lb:', af_nom[1])   \n",
    "### Total run time of the algorithm\n",
    "print('time', t2_af_nom- t1_af_nom, 'seconds')\n",
    "\n",
    "\n",
    "t1_af_nom_ub = process_time()\n",
    "af_nom_ub = af_nom_exp_pmin(p,R,slope,const+e_pl,1,par = h_par, par_u = 10)\n",
    "t2_af_nom_ub = process_time()\n",
    "\n",
    "print('sol:',af_nom_ub[0], 'ub:', af_nom_ub[1])   \n",
    "### Total run time of the algorithm\n",
    "print('time', t2_af_nom_ub- t1_af_nom_ub, 'seconds')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}