Skip to content

Commit

Permalink
Fix chapter 5 notebooks
Browse files Browse the repository at this point in the history
  • Loading branch information
@alessandrozocca
alessandrozocca committed Oct 27, 2023
1 parent 7c3d79b commit 3f02141
Show file tree
Hide file tree
Showing 4 changed files with 308 additions and 319 deletions.
175 changes: 71 additions & 104 deletions notebooks/05/markowitz_portfolio.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,17 +24,13 @@
"source": [
"## Preamble: Install Pyomo and a solver\n",
"\n",
"This cell selects and verifies a global SOLVER for the notebook.\n",
"\n",
"If run on Google Colab, the cell installs Pyomo and ipopt, then sets SOLVER to \n",
"use the ipopt solver. If run elsewhere, it assumes Pyomo and the Mosek solver\n",
"have been previously installed and sets SOLVER to use the Mosek solver via the Pyomo \n",
"SolverFactory. It then verifies that SOLVER is available."
"This cell selects and verifies a global SOLVER for the notebook. If run on Google Colab, the cell installs Pyomo and ipopt, then sets SOLVER to \n",
"use the ipopt solver. If run elsewhere, it assumes Pyomo and the Mosek solver have been previously installed and sets SOLVER to use the Mosek solver via the Pyomo SolverFactory. In both cases it then verifies that SOLVER is available."
]
},
{
"cell_type": "code",
"execution_count": 14,
"execution_count": 1,
"id": "5da22c67-5c34-4c3a-90a4-61222899e855",
"metadata": {
"tags": []
Expand All @@ -47,26 +43,24 @@
" !pip install idaes-pse --pre >/dev/null 2>/dev/null\n",
" !idaes get-extensions --to ./bin \n",
" os.environ['PATH'] += ':bin'\n",
" SOLVER = \"ipopt\"\n",
" \n",
" solver = \"ipopt\"\n",
"else:\n",
" SOLVER = \"mosek_direct\"\n",
" solver = \"mosek_direct\"\n",
"\n",
"import pyomo.environ as pyo\n",
"if not pyo.SolverFactory(SOLVER).available():\n",
" print(f\"Solver {SOLVER} is not available\")"
"SOLVER = pyo.SolverFactory(solver)\n",
"assert SOLVER.available(), f\"Solver {solver} is not available.\""
]
},
{
"cell_type": "code",
"execution_count": 15,
"execution_count": 2,
"id": "b13edf26",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"from IPython.display import Markdown, HTML\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt"
]
Expand All @@ -80,13 +74,13 @@
"source": [
"## Problem description and model formulation\n",
"\n",
"A canonical stochastic optimization problem is the so-called portfolio selection problem, also known as **Markowitz portfolio optimization**. Assume an investor has an initial capital $C$ that she wants to invest in $n$ possible risky assets, each of them with an unknown return rate $r_i$, $i=1,\\dots,n$, or in another risk-free asset with guaranteed return rate $R$. Let $x$ be the vector whose $i$-th component $x_i$ describes the amount invested in asset $i$ and $\\tilde{x}$ the amount invested in the risk-free asset. We consider a stochastic model where the return of the $n$ risky assets is then a random vector $r$ with known expected values $\\mu = \\mathbb E r $ and covariance\n",
"A canonical stochastic optimization problem is the so-called portfolio selection problem, also known as **Markowitz portfolio optimization**. Assume that an investor has an initial capital $C$ that she wants to invest in $n$ possible risky assets, each of them with an unknown return rate $r_i$, $i=1,\\dots,n$, or in another risk-free asset with a guaranteed return rate $R$. Let $x$ be the vector whose $i$-th component $x_i$ describes the amount invested in asset $i$ and $\\tilde{x}$ the amount invested in the risk-free asset. We consider a stochastic model in which the return of the $n$ risky assets is then a random vector $r$ with known expected values $\\mu = \\mathbb E r $ and covariance\n",
"\n",
"$$\n",
" \\Sigma = \\mathbb{E} [ (r-\\mu)(r-\\mu)^\\top].\n",
"$$\n",
"\n",
"The return of the investment $y = R \\tilde{x} + r^\\top x$ then also becomes a random variable with mean\n",
"The investment return $y = R \\tilde{x} + r^\\top x$ then becomes also a random variable with mean\n",
"\n",
"$$\n",
" \\mathbb{E} y = R \\tilde{x} + \\mathbb{E} r^\\top x = R \\tilde{x} + \\mu^\\top x\n",
Expand All @@ -98,9 +92,9 @@
" \\mathrm{Var}(y) = \\mathbb{E}(y-\\mathbb{E}y)^2 = x^\\top \\Sigma x.\n",
"$$\n",
"\n",
"The variance of the return of the investment is one possible way to quantify the risk of the investment $x$.\n",
"The variance of the investment return is a possible way to quantify the risk of the investment $x$.\n",
"\n",
"The problem the investor is facing is how to select a portfolio that achieves a good compromise between risk and expected return. More specifically, one could try to maximize the expected return $\\mathbb{E} y$ subject to an upper bound on the tolerable risk, obtaining the following optimization problem:\n",
"The investor needs to select a portfolio that achieves a good compromise between risk and expected return. One could try to maximize the expected return $\\mathbb{E} y$ subject to an upper bound on the tolerable risk, obtaining the following optimization problem:\n",
"\n",
"$$\n",
"\\begin{align*}\n",
Expand All @@ -113,54 +107,64 @@
"\\end{align*}\n",
"$$\n",
"\n",
"The first constraint describes the fact that the total amount invested must be equal to the initial capital. The second constraint ensures that the variance of the chosen portfolio is upper bounded by a parameter $\\gamma$, which captures the risk the investor is willing to undertake. The last nonnegativity constraint excludes the possibility of short-selling.\n",
"The first constraint describes the fact that the total amount invested must be equal to the initial capital. The second constraint ensures that the variance of the chosen portfolio is upper bounded by a parameter $\\gamma$, which captures the risk the investor is willing to take. The last non-negativity constraint excludes the possibility of short-selling.\n",
"\n",
"One can easily show that the quadratic constraint $x^\\top \\Sigma x \\leq \\gamma$ is convex in $x$ due to the fact that $\\Sigma$ is positive semidefinite, being a covariance matrix. Therefore, the optimization problem is convex. Let us implement it in Pyomo and solve it."
"One can easily show that the quadratic constraint $x^\\top \\Sigma x \\leq \\gamma$ is convex in $x$ due to the fact that $\\Sigma$ is positive semi-definite, being a covariance matrix. Therefore, the optimization problem is convex. Let us implement it in Pyomo and solve it."
]
},
{
"cell_type": "code",
"execution_count": 18,
"execution_count": 3,
"id": "0e194e8d",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"def markowitz(gamma, mu, Sigma):\n",
" m = pyo.ConcreteModel(\"Markowitz portfolio optimization\")\n",
"\n",
" m.xtilde = pyo.Var(domain=pyo.NonNegativeReals)\n",
" m.x = pyo.Var(range(n), domain=pyo.NonNegativeReals)\n",
"\n",
" @m.Objective(sense=pyo.maximize)\n",
" def objective(m):\n",
" return mu @ m.x + R * m.xtilde\n",
"\n",
" @m.Constraint()\n",
" def bounded_variance(m):\n",
" return (m.x @ (Sigma @ m.x)) <= gamma\n",
"\n",
" @m.Constraint()\n",
" def total_assets(m):\n",
" return sum(m.x[i] for i in range(n)) + m.xtilde == C\n",
"\n",
" SOLVER.solve(m)\n",
"\n",
" return m"
]
},
{
"cell_type": "markdown",
"id": "5f6de49c",
"metadata": {},
"source": [
"We now calculate the optimal investment strategy for a portfolio with $n=3$ assets, with expected returns $\\mu = (1.2, 1.1, 1.3)$ and a predefined covariance matrix. We set the risk parameter $\\gamma$ to $q$ and the risk-free return rate $R$ to $1.01$."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "fc935ce0",
"metadata": {},
"outputs": [
{
"data": {
"text/markdown": [
"**Solver status:** *ok, optimal*"
],
"text/plain": [
"<IPython.core.display.Markdown object>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/markdown": [
"**Solution:** $\\tilde x = 0.158$, $x_1 = 0.561$, $x_2 = 0.142$, $x_3 = 0.139$"
],
"text/plain": [
"<IPython.core.display.Markdown object>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/markdown": [
"**Maximizes objective value to:** $1.17$"
],
"text/plain": [
"<IPython.core.display.Markdown object>"
]
},
"metadata": {},
"output_type": "display_data"
"name": "stdout",
"output_type": "stream",
"text": [
"Optimal solution: xtilde = 0.158, x_1 = 0.561, x_2 = 0.142, x_3 = 0.139\n",
"Optimal investment return: 1.17\n"
]
}
],
"source": [
Expand All @@ -177,50 +181,21 @@
"# Check that Sigma is semi-definite positive\n",
"assert np.all(np.linalg.eigvals(Sigma) >= 0)\n",
"\n",
"# If you want to change the covariance matrix Sigma, ensure you input a semi-definite positive one.\n",
"# The easiest way to generate a random covariance matrix is first generating a random m x m matrix A\n",
"# and then taking the matrix A^T A (which is always semi-definite positive)\n",
"# When changing the matrix Sigma, make sure to input a semi-definite positive matrix.\n",
"# The easiest way to generate a random covariance matrix is as follows:\n",
"# first generate a random m x m matrix A and than take the matrix A A^T\n",
"# which is always semi-definite positive by construction.\n",
"#\n",
"# m = 3\n",
"# A = np.random.rand(m, m)\n",
"# Sigma = A.T @ A\n",
"\n",
"m = markowitz(gamma, mu, Sigma)\n",
"\n",
"def markowitz(gamma, mu, Sigma):\n",
" model = pyo.ConcreteModel(\"Markowitz portfolio optimization\")\n",
"\n",
" model.xtilde = pyo.Var(domain=pyo.NonNegativeReals)\n",
" model.x = pyo.Var(range(n), domain=pyo.NonNegativeReals)\n",
"\n",
" @model.Objective(sense=pyo.maximize)\n",
" def objective(m):\n",
" return mu @ m.x + R * m.xtilde\n",
"\n",
" @model.Constraint()\n",
" def bounded_variance(m):\n",
" return (m.x @ (Sigma @ m.x)) <= gamma\n",
"\n",
" @model.Constraint()\n",
" def total_assets(m):\n",
" return sum(m.x[i] for i in range(n)) + m.xtilde == C\n",
"\n",
" result = pyo.SolverFactory(SOLVER).solve(model)\n",
"\n",
" return result, model\n",
"\n",
"\n",
"result, model = markowitz(gamma, mu, Sigma)\n",
"\n",
"display(\n",
" Markdown(\n",
" f\"**Solver status:** *{result.solver.status}, {result.solver.termination_condition}*\"\n",
" )\n",
"print(\n",
" f\"Optimal solution: xtilde = {m.xtilde.value:.3f}, x_1 = {m.x[0].value:.3f}, x_2 = {m.x[1].value:.3f}, x_3 = {m.x[2].value:.3f}\"\n",
")\n",
"display(\n",
" Markdown(\n",
" f\"**Solution:** $\\\\tilde x = {model.xtilde.value:.3f}$, $x_1 = {model.x[0].value:.3f}$, $x_2 = {model.x[1].value:.3f}$, $x_3 = {model.x[2].value:.3f}$\"\n",
" )\n",
")\n",
"display(Markdown(f\"**Maximizes objective value to:** ${model.objective():.2f}$\"))"
"print(f\"Optimal investment return: {m.objective():.2f}\")"
]
},
{
Expand All @@ -233,7 +208,7 @@
},
{
"cell_type": "code",
"execution_count": 31,
"execution_count": 5,
"id": "10e5ab7d",
"metadata": {},
"outputs": [
Expand Down Expand Up @@ -296,23 +271,15 @@
"\n",
"plt.rcParams.update({\"font.size\": 14})\n",
"for gamma in gamma_values:\n",
" _, model = markowitz(gamma, mu, Sigma)\n",
" objective.append(round(model.objective(), 3))\n",
" m = markowitz(gamma, mu, Sigma)\n",
" objective.append(round(m.objective(), 3))\n",
"\n",
"plt.plot(gamma_values, objective, color=plt.cm.tab20c(0))\n",
"plt.xlabel(r\"Risk threshold $\\gamma$\")\n",
"plt.ylabel(\"Optimal objective value\")\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ad37924d",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
Expand Down
15 changes: 9 additions & 6 deletions notebooks/05/milk-pooling.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,20 +36,23 @@
"outputs": [],
"source": [
"import sys, os\n",
" \n",
"\n",
"if 'google.colab' in sys.modules:\n",
" !pip install idaes-pse --pre >/dev/null 2>/dev/null\n",
" !pip install highspy >/dev/null 2>/dev/null\n",
" !idaes get-extensions --to ./bin \n",
" os.environ['PATH'] += ':bin'\n",
" solver_NLO = \"ipopt\"\n",
"else:\n",
" solver_NLO = \"mosek\"\n",
"\n",
"import pyomo.environ as pyo\n",
"import pyomo.environ as pyo \n",
"\n",
"solver_LO = 'appsi_highs'\n",
"solver_LO = \"appsi_highs\"\n",
"SOLVER_LO = pyo.SolverFactory(solver_LO)\n",
"assert SOLVER_LO.available(), f\"Solver {solver_LO} is not available.\"\n",
"\n",
"solver_NLO = \"ipopt\"\n",
"SOLVER_NLO = pyo.SolverFactory(solver_NLO)\n",
"\n",
"assert SOLVER_LO.available(), f\"Solver {solver_LO} is not available.\"\n",
"assert SOLVER_NLO.available(), f\"Solver {solver_NLO} is not available.\""
]
},
Expand Down
170 changes: 92 additions & 78 deletions notebooks/05/ols-regression.ipynb
View file

Large diffs are not rendered by default.

Loading

0 comments on commit 3f02141

Please sign in to comment.