:\n",
+ " \n",
+ "Notice that the above optimization problem can be translated to an Ising Hamiltonian with the following change of variables $x_{i}=\\frac{1-\\sigma_{i}}{2}$ as:\n",
+ " $$H(s)=-\\sum_{j}h_{j}\\sigma_{j}-\\sum_{j
"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "f7f83e50-386f-4a6b-a137-fd3cbbb82a47",
+ "metadata": {
+ "id": "f7f83e50-386f-4a6b-a137-fd3cbbb82a47"
+ },
+ "source": [
+ " NOTE 2 (Hint):\n",
+ " \n",
+ "Encoding requires finding the minimum and maximum values for each one of the variables of the problem. You will need to find them.\n",
+ "Notice that $-z \\le r_{\\epsilon} \\le 0$, with $z$ a positive, real number. Are the variables needed to be normalized and why?\n",
+ "\n",
+ "
"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "9cabd342-4061-488b-9b3d-5d93d96344ec",
+ "metadata": {
+ "id": "9cabd342-4061-488b-9b3d-5d93d96344ec"
+ },
+ "source": [
+ " NOTE 3 (Hint):\n",
+ "\n",
+ "Ignore how the choice of encoding might affect the noise level of a quantum device and only consider how it affects the granularity of the problem as well as the number of binary variables needed for the encoding of each variable of the problem. One could also read\n",
+ "[A real world test of Portfolio Optimization with Quantum Annealing](https://arxiv.org/abs/2303.12601).\n",
+ "
"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "897a8e15-2da5-4ae2-9108-31aa005d6057",
+ "metadata": {
+ "id": "897a8e15-2da5-4ae2-9108-31aa005d6057"
+ },
+ "source": [
+ "**Answer:**"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "We use binary encoding due to its efficiency in representing the most number of values with the given bits. Additionally, all of our variables are normalized to $[0,1)$, so the variables represent fractions of constants or inputs (such as $z$).\n",
+ "\n",
+ "Assuming that we are using an p qubit quantum computer, where p is the number of qubits and is sufficiently large, the total number of binary variables we use is dependent on the level of granularity we want. This is because we can map each qubit 1 to 1 to the number of binary variables that our system has to keep track of.\n",
+ "\n",
+ "So, allocating our qubits (with an explanation of each variable afterwards):\n",
+ "\n",
+ "$x_n$: $n * [log_2 (n) + g]$\n",
+ "\n",
+ "- bits derived from investment granularity, where g is a variable denoting how fine we want the granularity, with $log_2 (n)$ being the minimum.\n",
+ "\n",
+ "$y_n$: $T$\n",
+ "\n",
+ "- bits dedicated to keeping track of current time\n",
+ "\n",
+ "$-r_ϵ$: $Φ$\n",
+ "\n",
+ "- bits dedicated to varying VaR\n",
+ "\n",
+ "\n",
+ "$S_{yt}$: $log_2(T)$\n",
+ "\n",
+ "- bits dedicated to y's slack variables\n",
+ "\n",
+ "$S_r: Φ$\n",
+ "\n",
+ "- bits dedicated to $r_ϵ$'s slack variable\n",
+ "\n",
+ "$S_{...}$: $Φ * T$\n",
+ "\n",
+ "- bits dedicated to slack variables for the second equality\n",
+ "\n",
+ "So, assuming p = 80, which is what we will be assuming for the rest of the challenge, so we end up using the binary variable distribution:\n",
+ "\n",
+ "$x_n$: $n * [log_2 (n) + g] = 18$\n",
+ "\n",
+ "\n",
+ "$y_n$: $T = 10$\n",
+ "\n",
+ "\n",
+ "$-r_ϵ$: $Φ = 4$\n",
+ "\n",
+ "\n",
+ "$S_{yt}$: $log_2(T) = 4$\n",
+ "\n",
+ "\n",
+ "$S_r$ : $Φ = 4$\n",
+ "\n",
+ "\n",
+ "$S_{...}$: $Φ * T = 40$\n",
+ "\n",
+ "\n",
+ "\n",
+ "Therefore the total binary variables used is 80.\n",
+ "\n"
+ ],
+ "metadata": {
+ "id": "rNa9yosXoxAV"
+ },
+ "id": "rNa9yosXoxAV"
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "553d8b31",
+ "metadata": {
+ "id": "553d8b31"
+ },
+ "outputs": [],
+ "source": [
+ "# write your code here\n",
+ "from pyqubo import Array, Binary, LogEncInteger, Placeholder, Constraint\n",
+ "\n",
+ "# Number of qubits we have\n",
+ "QUBIT_COUNT = 80\n",
+ "\n",
+ "# Number of binary variables per variable\n",
+ "g = 4\n",
+ "x_size = int(np.ceil(np.log2(number_assets)) + g)\n",
+ "y_size = 1\n",
+ "s_3_size = int(np.ceil(np.log2(T)))\n",
+ "phi = (QUBIT_COUNT - (x_size * number_assets + y_size * T + s_3_size)) // (T+2)\n",
+ "r_size = phi\n",
+ "s_1_size = phi\n",
+ "s_2_size = phi\n",
+ "\n",
+ "# Set up the variables\n",
+ "x = [1/(2**x_size - 1) * LogEncInteger(f\"x{i}\", (0, 2**x_size - 1)) for i in range(number_assets)]\n",
+ "y = Array.create('y', shape = T, vartype = 'BINARY')\n",
+ "# Note r is range of -1 to 0, instead of -z to 0\n",
+ "r = -1/(2**r_size - 1)*LogEncInteger(f\"r\", (0, 2**r_size - 1))\n",
+ "s_1 = 1/(2**s_1_size - 1)*LogEncInteger(f\"s_1\", (0, 2**s_1_size - 1))\n",
+ "s_2 = [1/(2**s_2_size - 1)*LogEncInteger(f\"s_2{i}\", (0, 2**s_2_size - 1)) for i in range(T)]\n",
+ "s_3 = 1/(2**s_3_size - 1)*LogEncInteger(f\"s_3\", (0, 2**s_3_size - 1))\n",
+ "\n",
+ "lambda_1 = Placeholder(\"lmd1\")\n",
+ "lambda_2 = [Placeholder(f\"lmbd2{i}\") for i in range(T)]\n",
+ "lambda_3 = Placeholder(\"lmd3\")\n",
+ "lambda_4 = Placeholder(\"lmd4\")\n",
+ "# Temporary values for lambdas\n",
+ "feed_dict = {\n",
+ " 'lmd1': 1,\n",
+ " 'lmd3': 1,\n",
+ " 'lmd4': 1\n",
+ "}\n",
+ "for i in range(T):\n",
+ " feed_dict[f\"lmbd2{i}\"] = 1\n",
+ "\n",
+ "z = 0.04238\n",
+ "M = -z-np.min(np.min(R))\n",
+ "e = 0.01\n",
+ "\n",
+ "H1 = 0 # The first term\n",
+ "for k in range(number_assets):\n",
+ " for j in range(number_assets):\n",
+ " H1 += x[k] * x[j] * covariance_matrix[k][j]\n",
+ "\n",
+ "H2 = 0 # The second term\n",
+ "for k in range(number_assets):\n",
+ " H2 += -expected_returns[k] * x[k]\n",
+ "\n",
+ "# The constraints\n",
+ "soft_constraint = lambda x: abs(x) < 1e-3\n",
+ "\n",
+ "H3 = lambda_1 * Constraint((-r * z + s_1 - z) ** 2, \"First Constraint\", soft_constraint)\n",
+ "\n",
+ "for t in range(T):\n",
+ " H3 += lambda_2[t] * Constraint((sum([R[i][t] * x[i] for i in range(number_assets)]) + M * y[t] + r + s_2[t] - M) ** 2, f\"Second Constraint{t}\", soft_constraint)\n",
+ "\n",
+ "H3 += lambda_3 * Constraint((-sum(y) + s_3 + T * (1-e)) ** 2, \"Third Constraint\", soft_constraint)\n",
+ "H3 += lambda_4 * Constraint((sum(x) - 1) ** 2, \"Fourth Constraint\", soft_constraint)\n",
+ "\n",
+ "H = H1 + H2 + H3\n",
+ "model = H.compile()\n",
+ "qubo, offset = model.to_qubo(feed_dict = feed_dict)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ed9d5b11-b2f3-4c79-bdf9-000cd7685f85",
+ "metadata": {
+ "id": "ed9d5b11-b2f3-4c79-bdf9-000cd7685f85"
+ },
+ "source": [
+ "\n",
+ " \n",
+ "BONUS EXERCISE:\n",
+ "As you notice the number of binary operators one has to use scales with the number of variables. Introducing the slack variables is therefore very costly. Try to suggest a new way of implementing the inequalities, without a lot of details. You can use the help of reference [Unbalanced penalization: A new approach to encode inequality constraints of combinatorial problems for quantum optimization algorithms](https://arxiv.org/abs/2211.13914).\n",
+ "\n",
+ "
"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "51466cca-e371-4f55-bdee-7b2eea6159d8",
+ "metadata": {
+ "id": "51466cca-e371-4f55-bdee-7b2eea6159d8"
+ },
+ "source": [
+ "**Answer:**\n",
+ "\n",
+ "Our current penalty terms for the constraints looks like\n",
+ "\n",
+ "$$\\sum_{j = 1}^{T+3}\\lambda_j\\left(\\sum_{i = 1}^{n+T+1}(A_{ji} x_i) - b_j+ s_j \\right)^2$$\n",
+ "or if we write out the digits of $s_j$,\n",
+ "$$\\sum_{j = 1}^{T+3}\\lambda_{j} \\left(\\sum_{i=1}^{n+T+1} (A_{ji}x_{i}) - b_j + \\beta_t\\sum_{k=1}^{D} 2^{-k} s_{jk}\\right)^2$$\n",
+ "where $s_{jk}$ are the digits of the slack variables of slack variable $s_j$. This means we have to use $D(T+3)$ extra binary digits for the slack variables.\n",
+ "\n",
+ "Instead, we implement the inequalities as unbalanced penalization functions, specifically an exponential decay. For QUBO to work, we need to approximate the exponential as\n",
+ "$$e^{-h(x)}-1 = -h(x)+\\frac{1}{2} h(x)^2$$\n",
+ "where we are trying to constrain $h(t)\\leq 0$, and the penalty to go to $0$ as $h$ goes to $0$. The constant term can be ignored in the minimization. So our new penalty would be\n",
+ "\n",
+ "$$\\sum_{j = 1}^{T+3}\\left[-\\lambda_{j,1}\\left(\\sum_{i = 1}^{n+T+1}(A_{ji} x_i) - b_j\\right)+\\lambda_{j,2}\\left(\\sum_{i = 1}^{n+T+1}(A_{ji} x_i) - b_j\\right)^2\\right]$$\n",
+ "\n",
+ "We have effectively traded $D(T+3)$ slack variable digits (so less qubits required) for more penalty term coefficients that we need to tune."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "9344d008-0e1a-45fd-9aa3-1399e56b859c",
+ "metadata": {
+ "id": "9344d008-0e1a-45fd-9aa3-1399e56b859c"
+ },
+ "source": [
+ "### Section 3: Simulated Annealing"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "1a9abb2c-9062-451e-9b4b-8009bd455e3c",
+ "metadata": {
+ "id": "1a9abb2c-9062-451e-9b4b-8009bd455e3c"
+ },
+ "source": [
+ "Optimization problems, like the one above, can be solved with the use of simulated annealing. It is one of the most preferred heuristic methods to solve optimization problems, since it avoids local minima."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "49533860-07ce-4bea-a2b8-c1a8cd3d2670",
+ "metadata": {
+ "id": "49533860-07ce-4bea-a2b8-c1a8cd3d2670"
+ },
+ "source": [
+ "- Using simulated annealing and the **data** calculate the ground state and find the best solution of the QUBO problem, $L_{fin}$, found in **Section 2**. Note that the lagrange multipliers in front of the penalty terms/constraints will need to be fine-tuned so that the equalities/inequalities of the problem are all satisfied. The output of the algorithm will give you the optimal values of the binary variables of the encoding of $\\{x_{i}\\}$, $r_{\\epsilon}$, $\\{y_{t}\\}$ and the slack variables so that the return is maximized, the risk is minimized while at the same time the inequality/equality constraints of the problem are all satisfied. In the next section you will learn how to evaluate your solutions."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "9a897de7-1efd-414e-bf16-a3f2353c9f6b",
+ "metadata": {
+ "id": "9a897de7-1efd-414e-bf16-a3f2353c9f6b"
+ },
+ "source": [
+ " NOTE 1 (Hint):\n",
+ "\n",
+ "\n",
+ "Feel free to first read any references provided online for Simulated Annealing like [Optimization Simulated Annealing](https://www.science.org/doi/10.1126/science.220.4598.671). You do not need to write your own code for this Section. Instead try [pyqubo](https://pyqubo.readthedocs.io/en/latest/getting_started.html#) and [D-Wave neal](https://docs.ocean.dwavesys.com/projects/neal/en/latest/).\n",
+ "
"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "54c2a1b7-feda-4461-bd63-b5e183e0b8f7",
+ "metadata": {
+ "id": "54c2a1b7-feda-4461-bd63-b5e183e0b8f7"
+ },
+ "source": [
+ " NOTE 2 (Hint):\n",
+ "\n",
+ " [pyqubo](https://pyqubo.readthedocs.io/en/latest/getting_started.html#) offers an easy way to evaluate whether the constraints are satisfied or not. Use it so that you are able to fine-tune the lagrange multipliers and get your solution to be as close as possible in satisfying the constraints.\n",
+ " \n",
+ " The lagrange multipliers act as relative weights between the different terms of the QUBO problem. Fine-tuning them so that the constraints of the problem are satisfied is **not** trivial. The fine-tuning process goes as following:\n",
+ " \n",
+ " Consider the following constraint optimization problem:\n",
+ " $$\\text{min} f(\\textbf{x})$$\n",
+ " subject to the constraint:\n",
+ " $$c(\\textbf{x})=0.$$\n",
+ "We can find a solution to this problem with the use of the **penalty method** approach, where we define a new objective function that includes the constraint as an added penalty term in the objective fuction. The $k^{th}$ step of the now unconstraint optimization problem is the following:\n",
+ "$$\\text{min}\\Psi_{k}(\\textbf{x})=f(\\textbf{x})+\\mu_{k}c(\\textbf{x})^{2}$$\n",
+ " with $\\mu_{k}$ the lagrange multiplier of this step (or else the penalty term coefficient). In order to fine-tune the langrange multiplier you start with a large-enough value of $\\mu_{1}$ (in comparison to the coefficients of the non-penalty terms of the QUBO problem). In the next steps of the iteration you increase the value of the coefficient until the constraints are satisfied.\n",
+ " \n",
+ " Specifically take a look at [pyqubo--validation of the constraints](https://pyqubo.readthedocs.io/en/latest/getting_started.html#validation-of-constraints). One could check whether the constraints are satisfied by making sure the energy contribution of each penalty term is 0. [pyqubo--validation of the constraints](https://pyqubo.readthedocs.io/en/latest/getting_started.html#validation-of-constraints) outputs the energy contribution of each penalty term. If you do not manage to get the energy contribution of the penalty terms to exactly 0, make sure they are as close as possible to that.\n",
+ "
"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ce724e8f-cc21-49fb-89b6-7e2b1e93296f",
+ "metadata": {
+ "id": "ce724e8f-cc21-49fb-89b6-7e2b1e93296f"
+ },
+ "source": [
+ " NOTE 3 (Hint):\n",
+ "The binary variables used for the encoding of $\\{x_{i}\\}$, $r_{\\epsilon}$, $\\{y_{t}\\}$ and all the slack variables affect the granularity of the problem. Choose the number of binary variables you use for the encoding of each one of the variables $\\{x_{i}\\}$, $r_{\\epsilon}$ and the slack variables. Why do you think that's a reasonable number of binary variables? How many binary variables do you need to use for the variables $\\{y_{t}\\}$?\n",
+ "
"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "a359a291-bbb0-4f4f-97fd-acd5f90d8eb6",
+ "metadata": {
+ "id": "a359a291-bbb0-4f4f-97fd-acd5f90d8eb6",
+ "colab": {
+ "base_uri": "https://localhost:8080/"
+ },
+ "outputId": "bf9e1174-ada9-49cb-a635-cd40b2f248c7"
+ },
+ "outputs": [
+ {
+ "output_type": "stream",
+ "name": "stdout",
+ "text": [
+ "{'Fourth Constraint': (True, -1.8041124150158794e-16), 'First Constraint': (True, 3.246334309293835e-08), 'Second Constraint2': (True, 1.5039738098141906e-07), 'Third Constraint': (True, 0.0008999999999996788), 'Second Constraint0': (True, 7.690880993535025e-09), 'Second Constraint1': (True, 2.8736800918583425e-08)}\n",
+ "[DecodedSolution({x2[5]:0, x2[1]:1, x2[0]:0, x1[4]:0, x1[3]:0, x1[1]:0, x1[0]:1, x0[5]:1, x0[2]:1, s_20[3]:0, s_22[10]:0, s_20[2]:0, s_20[9]:0, s_1[7]:0, s_20[1]:1, s_22[1]:0, x0[0]:0, s_1[6]:1, s_21[6]:1, s_21[8]:0, s_1[4]:0, s_20[0]:1, s_21[2]:1, s_3[1]:0, y[0]:1, s_1[1]:1, s_1[2]:0, y[2]:1, r[3]:1, y[1]:1, s_1[3]:1, s_20[4]:0, x2[3]:1, s_21[0]:0, s_22[5]:0, x2[2]:0, s_1[8]:0, s_22[6]:0, r[10]:0, s_1[5]:0, r[1]:1, r[5]:0, x1[2]:0, s_20[6]:0, s_21[10]:0, s_22[7]:0, r[7]:0, r[9]:0, r[8]:1, s_22[3]:1, r[0]:0, x2[4]:0, s_22[0]:0, r[2]:1, r[6]:0, x0[3]:0, s_21[9]:0, s_22[2]:1, r[4]:1, s_1[0]:1, s_1[10]:0, s_3[0]:0, s_20[7]:0, s_20[8]:1, s_21[1]:1, s_21[3]:1, x1[5]:0, s_21[4]:0, x0[4]:1, s_21[5]:0, s_21[7]:1, s_22[4]:0, s_22[8]:1, s_20[5]:1, s_22[9]:0, s_1[9]:0, s_20[10]:0, x0[1]:0}, energy=0.075229)]\n"
+ ]
+ }
+ ],
+ "source": [
+ "# write your code here\n",
+ "import neal\n",
+ "from itertools import product\n",
+ "\n",
+ "# Ranges for lambda\n",
+ "a, b = 100,102\n",
+ "ranges = {\n",
+ " 'lmd1':list(range(a, b)),\n",
+ " 'lmd3': list(range(a, b)),\n",
+ " 'lmd4': list(range(a, b))\n",
+ "}\n",
+ "for i in range(T):\n",
+ " ranges[f\"lmbd2{i}\"] = list(range(a, b))\n",
+ "\n",
+ "possibilities = list(product(*ranges.values()))\n",
+ "\n",
+ "sampler = neal.SimulatedAnnealingSampler()\n",
+ "feasible_sols =[]\n",
+ "for p in possibilities:\n",
+ " feed_dict = {key: value for key, value in zip(ranges.keys(), p)}\n",
+ " # qubo, offset = model.to_qubo(feed_dict = feed_dict)\n",
+ " bqm = model.to_bqm(feed_dict = feed_dict)\n",
+ " bqm.normalize()\n",
+ " sampleset = sampler.sample(bqm, num_reads = 10)#, sweeps = 1000, beta_schedule_type = 'custom', beta_schedule = [0, 1, 2, 3])\n",
+ " #sampleset = sampler.sample_qubo(qubo)\n",
+ " dec_samples = model.decode_sampleset(sampleset, feed_dict = feed_dict)\n",
+ " best = min(dec_samples, key = lambda x: x.energy)\n",
+ "\n",
+ " print(best.constraints())\n",
+ "\n",
+ " if not best.constraints(only_broken = True):\n",
+ " feasible_sols.append(best)\n",
+ "\n",
+ "print(feasible_sols)\n",
+ "best_feasible = min(feasible_sols,key=lambda x: x.energy)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "5e4786ea-1db6-46a6-b9f8-d1bfcdb88b20",
+ "metadata": {
+ "id": "5e4786ea-1db6-46a6-b9f8-d1bfcdb88b20"
+ },
+ "source": [
+ "### Section 4: Simulated Annealing--Verifying your solutions"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "37e4d030-3484-4fe2-b2de-7c2d374e58af",
+ "metadata": {
+ "id": "37e4d030-3484-4fe2-b2de-7c2d374e58af"
+ },
+ "source": [
+ "In the previous section you solved the **Mean-Variance-VaR** problem with Simulated Annealing, making sure that the penalty terms are imposed properly, so that the equalities/inequalities of the problem are satisfied. In this section, you will further test your solutions.\n",
+ "- From the output of **Section 3** you get the optimal values of the binary variables of the encoding. For example, for the variable $x_{i}$ we have:\n",
+ "$$x_{i}=\\sum_{d=1}^{D}f(d)z_{di}$$\n",
+ "with $\\{ z_{di} \\}$ the output of **Section 3** and $f(d)$ the function you have chosen for the encoding. Use this output in order to calculate the optimal percentages of the total portfolio invested in each one of the assets, $\\{x_{i}\\}$.\n",
+ "- After obtaining $\\{ x_{i} \\}$ use them to calculate the **expected total returns**: $$\\mu_{P}(x):=\\sum_{k=1}^{n}\\mu_{k}x_{k}$$\n",
+ " and **variance**: $$V(x):=\\sum_{k=1}^{n}\\sum_{j=1}^{n}x_{k}x_{j}\\sigma_{kj}$$\n",
+ "- In a similar way, from the output of **Section 3** calculate the $VaR_{\\epsilon}=-r_{\\epsilon}$.\n",
+ "- Run the code of **Section 3** multiple times. Each time calculate the expected returns, the variance and $VaR_{\\epsilon}$. Peak the best result and justify your answer.\n",
+ "- Calculate the historical VaR with the use of the following function:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "95761549",
+ "metadata": {
+ "id": "95761549"
+ },
+ "outputs": [],
+ "source": [
+ "def calculate_historical_VaR(weights, mu_R, confidence_level=0.95):\n",
+ "\n",
+ " portfolio_returns = np.dot(mu_R, np.transpose(np.array(weights)))\n",
+ "\n",
+ " VaR = np.percentile(portfolio_returns, 100 * (1 - confidence_level))\n",
+ "\n",
+ " return -VaR"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "09d0c42b",
+ "metadata": {
+ "id": "09d0c42b"
+ },
+ "source": [
+ "with $\\textit{weights}$ being a list that includes the optimal values (your solution) for the variables $\\{x_{i}\\}$ and $\\textit{mu_R}$ the expected returns of the assets (given by the data). Compare the return value got from this function to $VaR_{\\epsilon}$ calculated above."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "e32da760-2d7c-44cc-b541-0fd576eaedf2",
+ "metadata": {
+ "id": "e32da760-2d7c-44cc-b541-0fd576eaedf2"
+ },
+ "source": [
+ "**Answer:**"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "7c6d5972",
+ "metadata": {
+ "id": "7c6d5972"
+ },
+ "outputs": [],
+ "source": [
+ "x_value = []\n",
+ "for i in range(len(x)):\n",
+ " x_str = ''.join([str(best_feasible.sample[f\"x{i}[{pos}]\"]) for pos in range(x_size - 1, -1, -1)])\n",
+ " x_value.append(1/(2**x_size - 1) * int(x_str, base = 2))\n",
+ "print(x_value)\n",
+ "\n",
+ "mu_P = sum([expected_returns[k]*x_value[k] for k in range(number_assets)])\n",
+ "print(mu_P)\n",
+ "\n",
+ "V = 0\n",
+ "for k in range(number_assets):\n",
+ " for j in range(number_assets):\n",
+ " V += x_value[k] * x_value[j] * covariance_matrix[k][j]\n",
+ "print(V)\n",
+ "\n",
+ "r_value = None\n",
+ "r_str = ''.join([str(best_feasible.sample[f\"r[{pos}]\"]) for pos in range(r_size - 1, -1, -1)])\n",
+ "r_value = -1/(2**r_size - 1) * int(r_str, base = 2) * z\n",
+ "print(r_value)\n",
+ "\n",
+ "print(calculate_historical_VaR(x_value, expected_returns, confidence_level = 1 - e))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "9acc19d6-5b85-4ef4-aef4-7cdfe709b566",
+ "metadata": {
+ "id": "9acc19d6-5b85-4ef4-aef4-7cdfe709b566"
+ },
+ "source": [
+ "### Section 5: Quantum optimization (Quantum Annealing and/or QAOA)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ac6a4eab-e29a-470d-ba29-a3995c8eac10",
+ "metadata": {
+ "id": "ac6a4eab-e29a-470d-ba29-a3995c8eac10"
+ },
+ "source": [
+ "\n",
+ "**Note**: For the exercises below use the $\\hat{H}$ you found in **Section 2**, with the same encoding you used in that exercise and with the use of the **data** provided in the beginning of the challenge.\n",
+ "\n",
+ "**STEP 1:** Study the **Mean-Variance-VaR Simplified Version** problem with QAOA and Quantum Annealing.\n",
+ "\n",
+ "The Quantum Approximate Optimization Algorithm (QAOA) was first introduced in [A Quantum Approximate Optimization Algorithm](https://arxiv.org/abs/1411.4028). QAOA is a popular method to solve combinatorial optimization problems in noisy intermediate-scale quantum (NISQ) devices.\n",
+ "\n",
+ "Useful tutorials on implementing QAOA can be found here:\n",
+ "- [Pulser tutorial QAOA](https://pulser.readthedocs.io/en/latest/tutorials/qubo.html)\n",
+ "- [Qiskit QAOA](https://docs.quantum.ibm.com/api/qiskit/qiskit.algorithms.minimum_eigensolvers.QAOA), [Qiskit tutorial QAOA](https://qiskit.org/documentation/stable/0.40/tutorials/algorithms/05_qaoa.html)\n",
+ "- [pyQAOA tutorial QAOA](https://grove-docs.readthedocs.io/en/latest/qaoa.html)\n",
+ "- [PennyLane tutorial QAOA](https://pennylane.ai/qml/demos/tutorial_qaoa_intro/)\n",
+ "- [OpenQAOA-EntropicaLabs](https://openqaoa.entropicalabs.com/)\n",
+ "\n",
+ "For Quantum Annealing, you may read this reference [D-Wave, Quantum Annealing](https://docs.dwavesys.com/docs/latest/c_gs_2.html#getting-started-qa) and follow [D-Wave Ocean Software documentation](https://docs.ocean.dwavesys.com/en/stable/index.html) for the implementation.\n",
+ "In [Solving the Optimal Trading Trajectory Problem Using a Quantum Annealer](https://arxiv.org/abs/1508.06182), the authors explain how the choice of encoding might differ when considering solving an optimization problem with quantum annealing instead of simulated annealing.\n",
+ " \n",
+ "Do not implement anything at the moment. This step is for introducing you to different quantum or quantum-inspired algorithms one could utilize for this problem.\n",
+ "\n",
+ "**STEP 2:** In **Section 3** you used [D-Wave neal](https://docs.ocean.dwavesys.com/projects/neal/en/latest/) to solve the problem with Simulated Annealing. As also seen in the previous step, [D-Wave Ocean](https://docs.ocean.dwavesys.com/en/stable/getting_started.html) gives you immediate, secure access to D-Wave quantum and hybrid solvers. Study the capabilities of [D-Wave Ocean](https://docs.ocean.dwavesys.com/en/stable/getting_started.html) and answer the following questions:\n",
+ "\n",
+ "- For reasonable granularity (i.e. number of qubits per encoded variable) how many assets $n$ and events $T$ you can study with a quantum device and D-Wave? For that you will need to find the total number of qubits that you can explore with [D-Wave Ocean](https://docs.ocean.dwavesys.com/en/stable/getting_started.html) for the algorithm of your choice. Notice that the resources available to the public are much less than the capabilities of these devices.\n",
+ "\n",
+ "**STEP 3:** After exploring the capabilities of [D-Wave Ocean Software documentation](https://docs.ocean.dwavesys.com/en/stable/index.html), you are encouraged to use a **quantum simulator/backend** to study the **Mean-Variance-VaR Simplified Version** problem with the **data** provided in the beginning of the challenge (feel free to change $T$ and $n$ if needed) and any method of your choice. For this, you can use your qBraid account. Explore the options below:\n",
+ "\n",
+ "- [Qiskit](https://docs.quantum.ibm.com/api/qiskit)\n",
+ "- [Amazon Braket](https://docs.aws.amazon.com/braket/latest/developerguide/what-is-braket.html)\n",
+ "\n",
+ "Or even more, you may also opt to explore the usage of GPUs for simulating quantum algorithms. For example, see a recent [work](https://www.jpmorgan.com/technology/technology-blog/quantum-optimization-research) that studies QAOA scaling performance on a fast GPU-specific QAOA simulator.\n",
+ "\n",
+ "No matter what algorithm you choose, **you are encouraged to run small-scale simulations of the problems on quantum backends and simulators, always keeping in mind the resources that you are utilizing and their cost**. What problems/limitations do you encounter in comparison with simulated annealing?"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "dbc17d65-627d-4805-ac9f-3a3803bfc70d",
+ "metadata": {
+ "id": "dbc17d65-627d-4805-ac9f-3a3803bfc70d"
+ },
+ "source": [
+ " Note 1 (Hint):\n",
+ " \n",
+ " One could read the following references:\n",
+ "\n",
+ "- [Quantum Optimization: Potential, Challenges, and the Path Forward](https://arxiv.org/abs/2312.02279) and references within.\n",
+ " \n",
+ "
"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "8b4a5f08-14de-4d26-8031-8c1b2cd4f003",
+ "metadata": {
+ "id": "8b4a5f08-14de-4d26-8031-8c1b2cd4f003"
+ },
+ "source": [
+ " Note 2 (Hint):\n",
+ " \n",
+ "One can take a look at the [IBM Quantum Development Roadmap to 2033](https://newsroom.ibm.com/2023-12-04-IBM-Debuts-Next-Generation-Quantum-Processor-IBM-Quantum-System-Two,-Extends-Roadmap-to-Advance-Era-of-Quantum-Utility) and [QuEra's Quantum Computing Roadmap](https://www.quera.com/press-releases/quera-computing-releases-a-groundbreaking-roadmap-for-advanced-error-corrected-quantum-computers-pioneering-the-next-frontier-in-quantum-innovation) for an idea about the current and predicted quantum capabilities.\n",
+ " \n",
+ "
"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "990a82c7-6ccc-4552-8b86-97a95fece3cd",
+ "metadata": {
+ "id": "990a82c7-6ccc-4552-8b86-97a95fece3cd"
+ },
+ "source": [
+ "**Answer:**"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "Physically, the D-Wave quantum annealer is implemented as an array of superconducting loops for qubits, with current circulating clockwise or counterclockwise in order to represent +1s and -1s. In order to program a specific Hamiltonian function into the quantum annealer, bias terms can be applied to each qubit using external magnetic fields, and qubits can also be coupled together, which acts as a tuneable correlation weight between them. With the bias and correlation terms, any Hamiltonian of interest can be programmed, which is written as a quadratic function of the spins of each qubit. By making the change of variables x = (s+1)/2, the spins become 0s and 1s and the system has been recast as a quadratic unconstrained binary optimization (QUBO) problem. However, not every qubit can be correlated with every other qubit. There is a very specific network topology, so in order to represent the Hamiltonian as an arbitrary graph and embed it into the qubit network, some qubits need to be chained together, effectively acting as a single variable by highly correlating with each other.\n",
+ "\n",
+ "Quantum annealing with many variables is therefore harder than simulated annealing because of hardware limitations and noise. We have a much smaller set of qubits that we can use to run the simulation, and DWave chains multiple qubits together in order to make a single 'logical' qubit. There is also the issue of how the qubits are physically connected. Since not every qubit has a correlation with every other qubit (it is only connected to the ~10 nearest neighbors), not every possible graph can be represented without some kind of embedding, which could possibly take many more qubits. The upside is that because it is a physical process, it is much faster: even though it is an adiabatic process, which sounds like it should be slow, the timescale depends mainly on the temperature and scales like exp(sqrt(N)) rather than exp(N) because of quantum tunneling.\n",
+ "\n",
+ "Given that D-wave allows us to use 80 qubits, and at the reasonable granularity of 1/64 of total assets, we can simulate 3 assets for 6 weeks. This is an answer backed by data and will be covered briefly in the presentation."
+ ],
+ "metadata": {
+ "id": "ikLFsJ3FmbVW"
+ },
+ "id": "ikLFsJ3FmbVW"
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "The below code is our above code adapted for the dwave quantum annealer. This code accesses the DwaveSampler and performs the previous quantum annealment. We then convert from binary to readable data using the code from Problem 4."
+ ],
+ "metadata": {
+ "id": "VStBngFVkARt"
+ },
+ "id": "VStBngFVkARt"
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "aa60be68",
+ "metadata": {
+ "id": "aa60be68"
+ },
+ "outputs": [],
+ "source": [
+ "# write your code here\n",
+ "import neal\n",
+ "from itertools import product\n",
+ "from dwave.system import DWaveSampler, EmbeddingComposite\n",
+ "\n",
+ "# Ranges for lambda\n",
+ "a, b = 100,102\n",
+ "ranges = {\n",
+ " 'lmd1':list(range(a, b)),\n",
+ " 'lmd3': list(range(a, b)),\n",
+ " 'lmd4': list(range(a, b))\n",
+ "}\n",
+ "for i in range(T):\n",
+ " ranges[f\"lmbd2{i}\"] = list(range(a, b))\n",
+ "\n",
+ "possibilities = list(product(*ranges.values()))\n",
+ "\n",
+ "sampler = EmbeddingComposite(DWaveSampler())\n",
+ "feasible_sols =[]\n",
+ "for p in possibilities:\n",
+ " feed_dict = {key: value for key, value in zip(ranges.keys(), p)}\n",
+ " # qubo, offset = model.to_qubo(feed_dict = feed_dict)\n",
+ " bqm = model.to_bqm(feed_dict = feed_dict)\n",
+ " bqm.normalize()\n",
+ " sampleset = sampler.sample(bqm, num_reads = 10)#, sweeps = 1000, beta_schedule_type = 'custom', beta_schedule = [0, 1, 2, 3])\n",
+ " #sampleset = sampler.sample_qubo(qubo)\n",
+ " dec_samples = model.decode_sampleset(sampleset, feed_dict = feed_dict)\n",
+ " best = min(dec_samples, key = lambda x: x.energy)\n",
+ "\n",
+ " print(best.constraints())\n",
+ "\n",
+ " if not best.constraints(only_broken = True):\n",
+ " feasible_sols.append(best)\n",
+ "\n",
+ " if len(feasible_sols) > 0:\n",
+ " break\n",
+ "print(feasible_sols)\n",
+ "best_feasible = min(feasible_sols,key=lambda x: x.energy)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "95e10967",
+ "metadata": {
+ "id": "95e10967"
+ },
+ "source": [
+ "### Section 6 (BONUS):\n",
+ "\n",
+ "Above you followed specific steps and used simulated annealing as well as quantum or quantum-inspired methods to solve the problem of interest. To do this, you were instructed to incorporate the constraints of the problem (equalities/inequalities) in the objective function and then to find the QUBO formulation of the problem in question. In this exercise you are asked to think of another way to solve the **Mean-Variance-VaR Simplified Version** problem with quantum, quantum-inspired or hybrid methods, without following any of the steps mentioned in this challenge."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d57153d6",
+ "metadata": {
+ "id": "d57153d6"
+ },
+ "source": [
+ "**Answer:**"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "**Hybrid HHL**\n",
+ "\n",
+ "Portfolio optimization can be reformulated as convex quadratic program. This is equivalent to solving a linear system Ax=b using the method of Lagrange multipliers (so A is an augmented covariance matrix). Note that this is NOT encoded using binary variables yet. While we turned this into a QUBO problem for sections 1-5 so that we could use simulated and quantum annealing solvers, we do not necessarily have to do this.\n",
+ "\n",
+ "Instead, we propose an alternative way to solve this problem using a quantum/classical hybrid linear solvers. Quantum state representing the solution can be obtained by solving the corresponding QLSP using the HHL algorithm. This can be done because the covariance matrix is Hermitian. The resulting quantum state allows us to recover the optimal portfolio.\n",
+ "\n",
+ "HHL is a method of solving linear systems of equations using quantum computers, taking advantage of the speedup of the quantum Fourier transform. Noteably, HHL on its own is hard to implement on NISQ hardware because of the need for lots of connections between different qubits compared to the limited topology of the computers. In order to limit the number of swap (and other) gates in the implementation of the HHL algorithm, a hybrid approach can make it more efficient.\n",
+ "\n",
+ "Classical methods are used to determine how to scale the matrix A. In the foundational HHL article, the eigenvalues are restricted to [1/k, 1] to\n",
+ "account for the periodicity of the imaginary exponential and ensure well-conditioning. We need A to have a spectrum in this range, but it is more efficient to estimate λmax first so we don't waste bits for encoding unnecessarily and then scaling A by γ = 1/λmax, so that the maximum eigenvalue of γA is 1. We also use classical methods to guess a scaling parameter and then determine by how much it overestimates.\n",
+ "\n",
+ "Sources:\n",
+ "- \"Hybrid HHL with Dynamic Quantum Circuits on Real Hardware\" (https://arxiv.org/abs/2110.15958)\n",
+ "- Qiskit Documentation (https://docs.quantum.ibm.com/api/qiskit/0.35/qiskit.algorithms.linear_solvers.HHL)"
+ ],
+ "metadata": {
+ "id": "djc0kzxLe-rq"
+ },
+ "id": "djc0kzxLe-rq"
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "**QRNNs**\n",
+ "\n",
+ "Recurrent Neural Network (RNN) are a type of deep learning models in which the subcomponent, NN, contains state data that remembers information over time and is stacked together with other subcomponents, which is ideal for predicting sequences and time series. QRNN are the quantum version of RNNs: the component stacking is staggered to improve coherent time (the staggering removes the need for some qubits to have long term persistent states). It is a fully quantum ML model: input and output requires minor prior and post processing. QRNN performance on classical data has been confirmed on weather data, natural language, and stock data. It shows theoretical better performance than RNNs and other QNN with needing under 10 quantum tunable parameters and even fewer qubits.\n",
+ "\n",
+ "Source:\n",
+ "- Quantum Recurrent Neural Networks for Sequential Learning: https://arxiv.org/pdf/2302.03244.pdf"
+ ],
+ "metadata": {
+ "id": "eKHj_nxGjLfV"
+ },
+ "id": "eKHj_nxGjLfV"
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "067cea0c",
+ "metadata": {
+ "id": "067cea0c"
+ },
+ "outputs": [],
+ "source": [
+ "# write your code here"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "71ae8d5b-93a8-4305-95a1-0cebf1f29cf6",
+ "metadata": {
+ "id": "71ae8d5b-93a8-4305-95a1-0cebf1f29cf6"
+ },
+ "source": [
+ "### Section 7: Pitch your quantum strategy to a client"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ecb9a766-b6df-4b91-b0ff-13ced8fc5b52",
+ "metadata": {
+ "id": "ecb9a766-b6df-4b91-b0ff-13ced8fc5b52"
+ },
+ "source": [
+ "Imagine that you are part of the Quantum team at Moody's. You are meeting the board of stakeholders of an important company in the US. Given what you now know about the current quantum hardware capabilities, in which quantum algorithms should the stakeholders invest and why? Prepare your pitch focusing only on portfolio optimization problems. How do you think quantum could potentially avoid the drawbacks of classical optimization algorithms and why should a company invest in quantum today?"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "a440d8cd-2ed2-45e9-8a89-895847db823c",
+ "metadata": {
+ "id": "a440d8cd-2ed2-45e9-8a89-895847db823c"
+ },
+ "source": [
+ " Hint:\n",
+ " \n",
+ " One could read the following references:\n",
+ "\n",
+ "- [Quantum Optimization: Potential, Challenges, and the Path Forward](https://arxiv.org/abs/2312.02279)\n",
+ "- [QAOA with N⋅p≥200](https://arxiv.org/abs/2303.02064)\n",
+ "- [Evidence of Scaling Advantage for the Quantum Approximate Optimization Algorithm on a Classically Intractable Problem](https://arxiv.org/pdf/2308.02342.pdf)\n",
+ "- [Towards optimization under uncertainty for fundamental models in energy markets using quantum computers](https://arxiv.org/abs/2301.01108)\n",
+ " \n",
+ "
"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "Citations:\n",
+ "\n"
+ ],
+ "metadata": {
+ "id": "-NtU1rdF-92P"
+ },
+ "id": "-NtU1rdF-92P"
+ },
+ {
+ "cell_type": "markdown",
+ "id": "550f2bc1-7726-417e-88b5-8253ac12baeb",
+ "metadata": {
+ "id": "550f2bc1-7726-417e-88b5-8253ac12baeb"
+ },
+ "source": [
+ "**Answer:**\n",
+ "\n",
+ "After carefully reviewing the quantum algorithms suited to solving the problem of portfolio optimization and thoroughly testing the implementations of several of them, we would recommend that stakeholders invest in both Quantum Annealment and Hybrid HHL (Harrow-Hassidim-Lloyd). We recommend Quantum Annealment due to its reliability, existing infrastructure, and highly noise-resistant processes. Furthermore, our personal testing has confirmed the speed and reliability of popular Quantum Annealing processes. Next, Hybrid HHL manages to take advantage of quantum advantage for portfolio optimization while still using classical computing to dramatically increase the speed of several subroutines.\n",
+ "\n",
+ "Quantum manages to avoid many of the drawbacks of classical optimization algorithms by relying on physical processes that scale linearly with size of input, while classical processes scale exponentially. A company should invest in quantum in order to stay ahead in terms of computing power and take advantage of differences in the market.\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "c098599c-7ab0-48b1-9e33-51d22502310d",
+ "metadata": {
+ "id": "c098599c-7ab0-48b1-9e33-51d22502310d"
+ },
+ "source": [
+ "# This is the end of the challenge. Good luck!"
+ ]
+ }
+ ],
+ "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.5"
+ },
+ "colab": {
+ "provenance": []
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
\ No newline at end of file