Risk Measures and Portfolio Optimization

Prepared by: Madelynn Watson (mwatson4@nd.edu, 2023)

# Imports
import sys

if "google.colab" in sys.modules:
    !wget "https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/helper.py"
    import helper

    helper.easy_install()
else:
    sys.path.insert(0, "../")
    import helper
helper.set_plotting_style()

Introduction¶

In engineering, a common objective for optimization problems is maximizing overall profit; however, typically, static market prices are considered for different commodities which do not reflect actual market behavior. Failing to account for this uncertainty can lead to risk-inclined investments. Markowitz initially developed portfolio optimization in 1959 to de-risk financial asset portfolio returns in volatile market situations. In this framework proposed by Markowitz, a multi-objective optimization problem is considered where the expected value of the portfolio return is maximized, and the financial risk is minimized by selecting “weights” or quantities of each asset [1]. One challenge in this framework is how to quantify the financial risk. Throughout this notebook, we will examine five different risk measures (mean-variance (MV), mean-absolute deviation (MAD), Minimax (MM), value-at-risk (VaR), and conditional value-at-risk (CVaR)) and demonstrate their use in a motivating example for product portfolio optimization.

Common Risk Measures for Portfolio Optimization¶

Mean-Variance (MV) [1,2]¶

MV was the original quantifier of risk used by Markowitz in 1959. Here the risk is defined as the variance for a given expected return. The mathematical model is detailed in the following equations where $x_i$ and $x_j$ represent the quantity of assets i and j, and $\sigma_{i,j}$ represents the covariance of asset i and j.

\begin{align*} MV = \sum_{i=1}^{n} \sum_{j=1}^{n} \sigma_{i,j} x_i x_j \end{align*}

(1)

Although MV is easy to implement for simple portfolio problems, the computation of the covariance matrix makes it challenging to implement on large-scale problems. Additionally, if data is not normally distributed (a common challenge with real data) and there lies asymmetry in the probability distribution leading to a need for different risk measures to better quantify the risk-return tradeoffs.

Literature Examples Using MV¶

Xu, D., Bai, Z., Jin, X., Yang, X., Chen, S., & Zhou, M. (2022). A mean-variance portfolio optimization approach for high-renewable energy hub. Applied Energy, 325. https://doi.org/10.1016/j.apenergy.2022.119888

Chen, W., Zhang, H., Mehlawat, M. K., & Jia, L. (2021). Mean–variance portfolio optimization using machine learning-based stock price prediction. Applied Soft Computing, 100. https://doi.org/10.1016/j.asoc.2020.106943

Mean-Absolute Deviation (MAD) [1,3]¶

MAD was proposed by Konno and Yamazaki in 1991. The MAD model is a variant of the MV model in which the measure of risk is replaced by the absolute deviation. The mathematical model is detailed in the following equations where $x_j$ represents the quantity of asset j and $R_j$ represents the return of asset j.

\begin{align*} MAD = E \left[ \mid \sum_{j=1}^{n} R_j x_j -E \left[ \sum_{j=1}^{n} R_j x_j \right] \mid \right] \end{align*}

(2)

MAD can be employed for large-scale and highly diversified portfolio selection problems. One drawback of this risk measure is that it penalizes both positive and negative deviations.

Graphical Representation [3]¶

Literature Examples Using MAD¶

Hosseini-Nodeh, Z., Khanjani-Shiraz, R., & Pardalos, P. M. (2023). Portfolio optimization using robust mean absolute deviation model: Wasserstein metric approach. Finance Research Letters. https://doi.org/10.1016/j.frl.2023.103735

Ferreira, L., Borenstein, D., Righi, M. B., & de Almeida Filho, A. T. (2018). A fuzzy hybrid integrated framework for portfolio optimization in private banking. Expert Systems with Applications, 92, 350–362. https://doi.org/10.1016/j.eswa.2017.09.055

Minimax (MM) [1]¶

MM uses the minimum return as a measure of risk. In scenarios where the asset returns are multivariate and normally distributed, both MM and MV lead to the same result. This risk measure is simply shown below, where $M_p$ represents the minimum portfolio return.

\begin{align*} MM = M_p \end{align*}

(3)

MM has certain advantages when the returns are not normally distributed. Additionally, MM is fast due to its property of linear programming and can be employed for more complex models and constraints. One of the disadvantages of MM is its sensitivity to outliers. Hence, it cannot be used when the historical data is missing.

Literature Examples Using MM¶

Li, B., Sun, Y., Aw, G., & Teo, K. L. (2019). Uncertain portfolio optimization problem under a minimax risk measure. Applied Mathematical Modelling, 76, 274–281. https://doi.org/10.1016/j.apm.2019.06.019

Polak, G. G., Rogers, D. F., & Sweeney, D. J. (2010). Risk management strategies via minimax portfolio optimization. European Journal of Operational Research, 207(1), 409–419. https://doi.org/10.1016/j.ejor.2010.04.025

Value-at-Risk (VaR) [1,4]¶

VaR is a measure of how the market value of an asset is likely to reduce over a period of time under certain market conditions. VaR requires the determination of three parameters (i) time horizon, (ii) confidence level, and (iii) unit of VaR. The VaR at a given confidence level α corresponds to the 1−α percentile of the profit distribution, namely the lowest yearly profit after excluding all worse profits whose joint probability is at most 1−α. The formulation from Mutran et al. (2020) is shown below. Here the maximum V represents VaR, $R_j$ is the return of asset j, and 1[·] stands for the Heaviside step function, such that 1[x] = 1 for x $\geq$ 0 and 1[x]=0 otherwise.

\begin{align*} VaR_{\alpha} = \max V \\ s.t. \sum_{j=1}^{N} 1 \left[ V - R_j \right] \leq (1-\alpha)N \end{align*}

(4)

VaR can be misleading as it does not consider the worst-case loss. Additionally, VaR is discrete in nature and difficult to implement mathematically.

Literature Examples Using VaR¶

al Janabi, M. A. M., Ferrer, R., & Shahzad, S. J. H. (2019). Liquidity-adjusted value-at-risk optimization of a multi-asset portfolio using a vine copula approach. Physica A: Statistical Mechanics and Its Applications, 536. https://doi.org/10.1016/j.physa.2019.122579

Conditional Value-at-Risk (CVaR) [1,4]¶

CVaR is a scenario-based approach where the condition associated with VaR is the weighted average of VaR and losses exceeding VaR. This quantity represents the expected value of the 100(1-α) % worst scenarios at a given confidence level α. CVaR is an alternative to VaR that is more sensitive to the shape of the tail of the scenario distribution. Additionally, CVaR is a coherent and convex measure of risk.

\begin{align*} CVaR_{\alpha} = VaR_{\alpha} - \frac{1}{(1-\alpha)N} \sum_{j=1}^{N} max \left[ 0, VaR_{\alpha} - R_j \right] \end{align*}

(5)

One disadvantage of CVaR is that it cannot indicate the maximum loss that can be incurred.

The lecture notes write CVaR the other way around: they minimize the CVaR of the loss, while this notebook maximizes the CVaR of the profit. Substituting loss = -profit turns one into the other, so the two are the same measure at the same confidence level $\alpha$ . Birge and Louveaux, Introduction to Stochastic Programming, 2nd ed., §2.9, p. 85, give the minimizing form.

Graphical Representation [4]¶

In this figure EP represents the expected value of the profit distribution (average profit) where profit is representative of the portfolio return

Literature Examples Using CVaR¶

Mutran, V. M., Ribeiro, C. O., Nascimento, C. A. O., & Chachuat, B. (2020). Risk-conscious optimization model to support bioenergy investments in the Brazilian sugarcane industry. Applied Energy, 258. https://doi.org/10.1016/j.apenergy.2019.113978

Xuan, A., Shen, X., Guo, Q., & Sun, H. (2021). A conditional value-at-risk based planning model for integrated energy system with energy storage and renewables. Applied Energy, 294. https://doi.org/10.1016/j.apenergy.2021.116971

Motivating Example: Product Portfolio Optimization in the Brazilian Sugarcane Industry¶

Brazil is the world leader in sugar production and exportation. In 2015, sugarcane GDP reached US$ 28.5 billion, with a production scale that achieved more than 30 million tons of sugar and 21 million cubic meters of ethanol [5]. The price of ethanol in Brazil is limited by the domestic price of oil and international fuel market conditions, and sugar prices are defined at the international commodity market. Weather conditions, international sugar prices, and instability of governmental fuel policies and subsidies create huge uncertainties in these markets [6]. Therefore, it is essential to quantify the financial risk from these volatile markets. Here we use portfolio optimization and three different risk measures to mitigate the impacts of financial risk in this industry by changing the sugarcane mill product portfolio.

Problem Definition¶

For this problem, we model a sugarcane mill that can produce sugar, ethanol, and electricity with the ability to sell electricity to free and regulated markets. The overall capacity of the sugarcane mill is 3,000,000 tonnes of sugarcane, and conversion and cost data are provided.

Process Flow Diagram¶

Process Model¶

Sets

\begin{align*} p \in P: \text{Products} \\ s \in S: Saleable\ Products \\ c \in C: Commodities \\ u \in U: Process\ Units \\ q \in Q: Historical\ Price\ Scenarios \end{align*}

(6)

Parameters

\begin{align*} Ca: Yearly\ production\ of\ sugarcane \\ \tau_{u}: Maximum\ capacity\ of\ process\ unit\ u \\ PC_{s}: Production\ cost\ of\ saleable\ product\ s \\ \theta_{f,u,p}: Conversion\ of\ resource\ f\ in\ process\ unit\ u\ to\ product\ p \\ \gamma_{u,p}: Generation\ of\ product\ p\ in\ unit\ u \\ HP_{q,c}: Historical\ Prices\ at\ price\ scenario\ q\\ \pi_{reg}: Price \ of \ electricity \ at \ the \ regulated \ market \end{align*}

(7)

Variables

\begin{align*} x_{p}: Amount\ of\ product\ p \\ profit_{q}: Profit\ for\ each\ historical\ price\ scenario\ q\\ Eprofit: Expected\ value\ of\ the\ profit \end{align*}

(8)

Constraints

Mass Balances

\begin{align*} x_{jui} = Ca \cdot \theta_{cane,mill,jui} \\ x_{bag} = Ca \cdot \theta_{cane,mill,bag} \\ x_{jui} = x_{jui-fact} + x_{jui-dist} \\ x_{sug} = x_{jui-fact} \cdot \theta_{jui,fact,sug} \\ x_{mol} = x_{sug} \cdot \gamma_{mol,fact} \\ x_{eth} = x_{jui-dist} \cdot \theta_{jui,dist,eth} + x_{mol} \cdot \theta_{mol,dist,eth} \\ x_{el-r} = 0.053 \cdot Ca \\ x_{reg} + x_{free} = x_{el-r} \end{align*}

(9)

Profit Calculations

\begin{align*} profit_{q} = \sum_{c \in C} x_{c} \cdot HP_{q,c} + x_{reg} \cdot \pi_{reg} - \sum_{s \in S} x_{s} \cdot PC_{s} \ \ \forall \ \ q \in Q\\ Eprofit = \frac{1}{Q} \sum_{q \in Q} profit_{q} \\ \end{align*}

(10)

Degree of Freedom (DOF) Analysis¶

Number of Variables

\begin{align*} x_{p} \cdot 10 \ products \\ profit_{Q} \cdot Q \\ Eprofit \cdot 1 \\ Total\ Variables = 10 + Q + 1 \end{align*}

(11)

Number of Constraints

\begin{align*} 8 \ mass \ balances \\ profit \ scenarios \cdot Q \\ Eprofit \cdot 1 \\ Total\ Constraints = 8 + Q + 1 \end{align*}

(12)

DOF

\begin{align*} DOF = Variables - Constraints \\ DOF = 10 + Q + 1 - 8 - Q - 1 \\ DOF = 2 \end{align*}

(13)

Visualize Input Data¶

import logging

import pandas as pd
import numpy as np
import pyomo.environ as pyo
import matplotlib.pyplot as plt
from pyomo.environ import units as u
from pyomo.util.check_units import assert_units_consistent
from pyomo.core.base.units_container import InconsistentUnitsError

# Pyomo's unit library has tonnes, cubic meters, and MWh but no money, so declare
# a currency.
u.load_definitions_from_strings(["USD = [currency]"])

# Shorthands for the compound units this problem uses
USD_PER_TONNE = u.USD / u.metric_ton
USD_PER_M3 = u.USD / u.m**3
USD_PER_MWH = u.USD / u.MWh
M3_PER_TONNE = u.m**3 / u.metric_ton
MWH_PER_TONNE = u.MWh / u.metric_ton

# Load Data From CSV files stored on Github
path_cap = "https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/data/riskmeasures_capacity.csv"
path_opex = "https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/data/riskmeasures_opex.csv"
path_conv = "https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/data/riskmeasures_conversions.csv"
path_gen = "https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/data/riskmeasures_generation.csv"
path_hp = "https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/data/riskmeasures_historicalprices.csv"
df_maxcap = pd.read_csv(path_cap)
df_prodcost = pd.read_csv(path_opex)
df_conv = pd.read_csv(path_conv)
df_gen = pd.read_csv(path_gen)
df_hp = pd.read_csv(path_hp)

## Uncomment to Display Additional Data in Tables
# display(df_maxcap)
# display(df_prodcost)
# display(df_conv)
# display(df_gen)

# Plot historical prices
x = np.arange(1, len(df_hp["eth"]) + 1)
fig, ax = plt.subplots(figsize=(6.4, 4))
plt.scatter(x, df_hp["eth"], 1, label="Ethanol")
plt.xlabel("Week (2013 - 2022)", fontsize=16, fontweight="bold")
plt.ylabel("USD/m$^3$", fontsize=16, fontweight="bold")
plt.legend(fontsize=14)
ax.xaxis.set_tick_params(labelsize=15)
ax.yaxis.set_tick_params(labelsize=15)
ax.tick_params(direction="in")
plt.show()

fig, ax = plt.subplots(figsize=(6.4, 4))
plt.scatter(x, df_hp["sug"], 1, label="Sugar")
plt.xlabel("Week (2013 - 2022)", fontsize=16, fontweight="bold")
plt.ylabel("USD/tonne", fontsize=16, fontweight="bold")
plt.legend(fontsize=14)
ax.xaxis.set_tick_params(labelsize=15)
ax.yaxis.set_tick_params(labelsize=15)
ax.tick_params(direction="in")
plt.show()

fig, ax = plt.subplots(figsize=(6.4, 4))
plt.scatter(x, df_hp["fre"], 1, label="Electricity to Free Market")
plt.xlabel("Week (2013 - 2022)", fontsize=16, fontweight="bold")
plt.ylabel("USD/MWh", fontsize=16, fontweight="bold")
plt.legend(fontsize=14)
ax.xaxis.set_tick_params(labelsize=15)
ax.yaxis.set_tick_params(labelsize=15)
ax.tick_params(direction="in")
plt.show()

Define Base Model in Pyomo¶

Units. This mill sells three commodities measured in three different physical dimensions: sugar in tonnes, ethanol in cubic meters, and electricity in MWh. Pyomo attaches units to a component, not to an index, so one Var indexed over all nine resources cannot carry them. The builder below declares one variable per dimension and stitches them back together with pyo.Reference, so the rest of the notebook can still write m.x["eth"]. Every model then ends with assert_units_consistent, which raises if any constraint or the objective does not balance. A units declaration nobody checks only looks verified.

# Create Base Model In Pyomo
def create_model():
    """
    This function builds a superstructure model in Pyomo for a sugarcane mill that
    can produce sugar, ethanol, and electricity.

    Returns: Pyomo model m
    """
    m = pyo.ConcreteModel()

    # SETS
    resources = ["sug", "eth", "fre", "reg", "cane", "jui", "mol", "bag", "el-r"]
    saleable_products = ["sug", "eth", "fre", "reg"]
    processes = ["mill", "fact", "dist", "cogen"]
    commodity = ["fre", "sug", "eth"]

    # The nine resources are measured in three different physical dimensions.
    mass = ["sug", "cane", "jui", "mol", "bag"]  # tonne
    volume = ["eth"]  # m3
    energy = ["fre", "reg", "el-r"]  # MWh

    # Price of each saleable product carries the matching units
    price_units = {
        "sug": USD_PER_TONNE,
        "eth": USD_PER_M3,
        "fre": USD_PER_MWH,
        "reg": USD_PER_MWH,
    }

    # PARAMETERS
    # Scalars
    Ca = 3000000 * u.metric_ton  # Annual Sugarcane Capacity
    price_reg = 72.5 * USD_PER_MWH  # Price of electricity sold to the regulated market

    # Fill in Dictionaries with Excel Data
    max_cap = {}
    prodcost = {}
    conv = {}
    gen = {}
    hp = {}

    for i in range(len(df_maxcap["process"])):
        max_cap[df_maxcap["process"][i]] = df_maxcap["capacity"][i] * u.metric_ton

    for i in range(len(df_prodcost["saleable_product"])):
        s = df_prodcost["saleable_product"][i]
        prodcost[s] = df_prodcost["cost"][i] * price_units[s]

    for k in resources:
        for i in range(len(df_conv[k])):
            conv[(df_conv["resource"][i], df_conv["process"][i], k)] = df_conv[k][i]

    for k in resources:
        for i in range(len(df_gen[k])):
            gen[(df_gen["process"][i], k)] = df_gen[k][i]

    for k in commodity:
        for i in range(len(df_hp[k])):
            hp[(df_hp["q"][i], k)] = df_hp[k][i] * price_units[k]

    N = len(df_hp["q"])
    q = []
    for i in range(1, N + 1):
        q.append("t%d" % (i))

    # PYOMO SETS
    m.resources = pyo.Set(initialize=resources)
    m.saleable_products = pyo.Set(initialize=saleable_products)
    m.processes = pyo.Set(initialize=processes)
    m.price_time_obs = pyo.Set(initialize=q)
    m.commodities = pyo.Set(initialize=commodity)

    m.N = pyo.Param(initialize=N)

    # VARIABLES
    # Positive Variables. Pyomo attaches units to a component, not to an index, so
    # the amount of each resource needs one variable per physical dimension.
    m.x_mass = pyo.Var(mass, domain=pyo.NonNegativeReals, units=u.metric_ton)
    m.x_volume = pyo.Var(volume, domain=pyo.NonNegativeReals, units=u.m**3)
    m.x_energy = pyo.Var(energy, domain=pyo.NonNegativeReals, units=u.MWh)

    # ... and a Reference stitches the three back into a single indexed view, so
    # the rest of the notebook can write m.x["eth"] regardless of dimension.
    m.x = pyo.Reference(
        {
            **{i: m.x_mass[i] for i in mass},
            **{i: m.x_volume[i] for i in volume},
            **{i: m.x_energy[i] for i in energy},
        },
        ctype=pyo.Var,
    )

    # Juice split between the factory and the distillery
    m.r = pyo.Var(["fact", "dist"], domain=pyo.NonNegativeReals, units=u.metric_ton)
    m.profit = pyo.Var(m.price_time_obs, units=u.USD)
    m.EP = pyo.Var(units=u.USD)

    # CONSTRAINTS
    # Superstructure Constraints
    def mill1(m):
        return m.x["jui"] == Ca * conv["cane", "mill", "jui"]

    m.milleq1 = pyo.Constraint(rule=mill1)

    def mill2(m):
        return m.x["bag"] == Ca * conv["cane", "mill", "bag"]

    m.milleq2 = pyo.Constraint(rule=mill2)

    def juice1(m):
        return m.x["jui"] == m.r["fact"] + m.r["dist"]

    m.juiceeq1 = pyo.Constraint(rule=juice1)

    def juice2(m):
        return m.r["fact"] <= max_cap["fact"]

    m.juiceeq2 = pyo.Constraint(rule=juice2)

    def juice3(m):
        return m.r["dist"] <= max_cap["dist"]

    m.juiceeq3 = pyo.Constraint(rule=juice3)

    def sugar1(m):
        return m.x["sug"] == conv["jui", "fact", "sug"] * m.r["fact"]

    m.su1 = pyo.Constraint(rule=sugar1)

    def sugar2(m):
        return m.x["mol"] == gen["fact", "mol"] * m.x["sug"]

    m.su2 = pyo.Constraint(rule=sugar2)

    def ethanol1(m):
        return (
            m.x["eth"]
            == conv["jui", "dist", "eth"] * M3_PER_TONNE * m.r["dist"]
            + conv["mol", "dist", "eth"] * M3_PER_TONNE * m.x["mol"]
        )

    m.et1 = pyo.Constraint(rule=ethanol1)

    def el_prod(m):
        # 53 kWh produced per tonne of sugarcane processed
        return m.x["el-r"] == Ca * 0.053 * MWH_PER_TONNE

    m.el_produced = pyo.Constraint(rule=el_prod)

    def el_sales(m):
        return m.x["fre"] + m.x["reg"] == m.x["el-r"]

    m.electricity = pyo.Constraint(rule=el_sales)

    # Expected Profit
    def profit(m, q):
        return m.profit[q] == m.x["sug"] * hp[q, "sug"] + m.x["eth"] * hp[
            q, "eth"
        ] + m.x["fre"] * hp[q, "fre"] + m.x["reg"] * price_reg - sum(
            m.x[j] * prodcost[j] for j in m.saleable_products
        )

    m.profiteq = pyo.Constraint(m.price_time_obs, rule=profit)

    def eprofit1(m):
        return m.EP == (1 / N) * sum(m.profit[q] for q in m.price_time_obs)

    m.ep1 = pyo.Constraint(rule=eprofit1)

    # Monetary exposure. The mean historical price of each commodity is the
    # conversion factor that puts tonnes, cubic meters, and MWh on one basis:
    # exposure[c] is the revenue exposed to the price of commodity c, in USD.
    p_bar = {c: df_hp[c].mean() * price_units[c] for c in commodity}

    def exposure(m, c):
        return p_bar[c] * m.x[c]

    m.exposure = pyo.Expression(m.commodities, rule=exposure)

    # Product Revenues to be used in analysis
    def sug_prof(m):
        return sum(m.x["sug"] * hp[q, "sug"] for q in m.price_time_obs) * (1 / N)

    m.avg_sug_prof = pyo.Expression(rule=sug_prof)

    def eth_prof(m):
        return sum(m.x["eth"] * hp[q, "eth"] for q in m.price_time_obs) * (1 / N)

    m.avg_eth_prof = pyo.Expression(rule=eth_prof)

    def el_prof(m):
        return sum(m.x["fre"] * hp[q, "fre"] for q in m.price_time_obs) * (1 / N)

    m.avg_el_prof = pyo.Expression(rule=el_prof)

    return m

Maximize Profit Only (No Risk Measure)¶

# Load the base model from the function created above
m = create_model()


# Set the objective to maximize the expected value of the profit
def obj_rule(m):
    return m.EP


m.obj = pyo.Objective(rule=obj_rule, sense=pyo.maximize)

# Every constraint and the objective must balance dimensionally
assert_units_consistent(m)

# Solve the Model
sol = pyo.SolverFactory("ipopt")

results = sol.solve(m)
assert pyo.check_optimal_termination(results), (
    f"Solve failed: status={results.solver.status}, "
    f"termination={results.solver.termination_condition}"
)

print("Results")
print("----------------------------------------------------")
print("Expected Profit", np.round(pyo.value(m.obj), 2), "USD")
print("Ethanol Produced", np.round(pyo.value(m.x["eth"]), 2), "m3")
print("Sugar Produced", np.round(pyo.value(m.x["sug"]), 2), "tonne")
print("Electricity Produced", np.round(pyo.value(m.x["el-r"]), 2), "MWh")
print("Electricity to Free Market", np.round(pyo.value(m.x["fre"]), 2), "MWh")
print("Electricity to Regulated Market", np.round(pyo.value(m.x["reg"]), 2), "MWh")

Results
----------------------------------------------------
Expected Profit 64780608.4 USD
Ethanol Produced 155062.24 m3
Sugar Produced 189129.78 tonne
Electricity Produced 159000.0 MWh
Electricity to Free Market 159000.0 MWh
Electricity to Regulated Market -0.0 MWh

Visualize Results¶

Profit Distribution¶

# Collect the profit distribution
profits = []
for i in m.price_time_obs:
    profits.append(pyo.value(m.profit[i]))

# Collect the minimum profit
min_prof = min(profits)

# Plot the Profit Distribution
fig, ax = plt.subplots(figsize=(6.4, 4))
plt.hist(np.array(profits) / 1e6, alpha=0.5)
plt.vlines(
    ymin=0,
    ymax=160,
    x=pyo.value(m.obj) / 1e6,
    label="Expected Profit",
    color="green",
    linestyles="solid",
    linewidth=3,
)
plt.vlines(
    ymin=0,
    ymax=160,
    x=min_prof / 1e6,
    label=("Minimum Profit"),
    color="red",
    linestyles="dashed",
    linewidth=3,
)
plt.xlabel("profit$_q$ M USD", fontsize=16, fontweight="bold")
plt.ylabel("Frequency", fontsize=16, fontweight="bold")
plt.legend(fontsize=14)
ax.xaxis.set_tick_params(labelsize=15)
ax.yaxis.set_tick_params(labelsize=15)
ax.tick_params(direction="in")
plt.show()

print("Results")
print("-------------------------------------")
print("Expected Profit:", np.round(pyo.value(m.obj) / 1e6, 2), "M USD")
print("Minimum Profit:", np.round(min_prof / 1e6, 2), "M USD")
difference = pyo.value(m.obj) - min_prof
print("Difference:", np.round(difference / 1e6, 2), "M USD")

# Collect Results for Conclusion
final_EP = {}
final_min = {}
final_diff = {}
final_exposure = {}
final_profits = {}

final_EP["No Risk"] = pyo.value(m.obj) / 1e6
final_min["No Risk"] = min_prof / 1e6
final_diff["No Risk"] = difference / 1e6
final_exposure["No Risk"] = {c: pyo.value(m.exposure[c]) for c in m.commodities}
final_profits["No Risk"] = np.array(profits)

Results
-------------------------------------
Expected Profit: 64.78 M USD
Minimum Profit: -39.0 M USD
Difference: 103.78 M USD

Product Distribution¶

fig, ax = plt.subplots(figsize=(6.4, 4))
plt.bar("Sugar", pyo.value(m.avg_sug_prof) / 1e6)
plt.bar("Ethanol", pyo.value(m.avg_eth_prof) / 1e6)
plt.bar("Electricity \n Free \n Market", pyo.value(m.avg_el_prof) / 1e6)
plt.bar("Electricity \n Regulated \n Market", pyo.value(m.x["reg"]) * 72.5 / 1e6)
plt.xlabel("Product", fontsize=16, fontweight="bold")
plt.ylabel("Revenue from \n Each Product (M USD)", fontsize=16, fontweight="bold")
ax.xaxis.set_tick_params(labelsize=15)
ax.yaxis.set_tick_params(labelsize=15)
ax.tick_params(direction="in")
plt.show()

Solve Model with Risk Measures¶

Two Objectives, Not One¶

The introduction described portfolio optimization as a multi-objective problem: the expected profit is maximized and the risk is minimized. Those two goals conflict, so there is no single best portfolio. There is instead a set of portfolios, the Pareto set, in which no portfolio can be improved on one objective without giving up something on the other.

The $\epsilon$ -constraint method finds that set one point at a time. Keep one objective, and move the other into a constraint whose right-hand side $\rho$ we choose:

\begin{align*} minimize \ \ risk \\ s.t. \ \ Eprofit \geq \rho \end{align*}

(14)

Each value of $\rho$ gives one Pareto-optimal portfolio, and sweeping $\rho$ from the lowest return worth considering up to the highest attainable one traces the whole frontier. This is the form the lecture notes use for the Markowitz problem. Recall from the DOF analysis that the base model has 2 degrees of freedom; the required-return constraint spends one of them, leaving one to trade risk against return.

The alternative is a weighted sum: minimize $risk - \lambda \cdot Eprofit$ and sweep the weight $\lambda$ . That is easier to write, but a weighted sum can only ever return portfolios on the convex hull of the frontier. Where the frontier has a non-convex (dented) section, no value of $\lambda$ returns a point inside the dent, and those Pareto-optimal portfolios stay invisible no matter how finely $\lambda$ is swept. The $\epsilon$ -constraint method has no such blind spot, which is why it is used here: all three risk measures below are solved subject to the same required return, so the portfolios they choose can be compared.

# The maximum attainable expected profit. This is the top of every frontier
# below, and it sets the scale for the required return rho.
sol = pyo.SolverFactory("ipopt")

m = create_model()
m.obj = pyo.Objective(expr=m.EP, sense=pyo.maximize)
results = sol.solve(m)
assert pyo.check_optimal_termination(results), (
    f"Solve failed: status={results.solver.status}, "
    f"termination={results.solver.termination_condition}"
)
EP_max = pyo.value(m.EP)

# Required return for the three risk measures below: 95% of the maximum, so we
# are willing to give up 5% of the expected profit to buy a reduction in risk.
rho_required = 0.95 * EP_max

print("Maximum expected profit :", np.round(EP_max / 1e6, 3), "M USD")
print("Required return (rho)   :", np.round(rho_required / 1e6, 3), "M USD")


def add_required_return(m, rho):
    """Attach the epsilon-constraint Eprofit >= rho to a model.

    Arguments:
        m: a Pyomo model from create_model()
        rho: required expected profit [USD]
    """
    m.required_return = pyo.Constraint(expr=m.EP >= rho * u.USD)


def initialize_from_max_profit(m):
    """Solve the expected-profit problem and leave its solution in m.

    Ipopt cannot start the MAD or CVaR models from all zeros, so call this
    BEFORE adding any risk-measure variables and use the result as the
    starting point.

    Arguments:
        m: a Pyomo model from create_model()
    """
    m.init_obj = pyo.Objective(expr=m.EP, sense=pyo.maximize)
    results = sol.solve(m)
    assert pyo.check_optimal_termination(results), (
        f"Initialization solve failed: status={results.solver.status}, "
        f"termination={results.solver.termination_condition}"
    )
    m.init_obj.deactivate()

Maximum expected profit : 64.781 M USD
Required return (rho)   : 61.542 M USD

Formulate the MV Objective¶

Additional Parameters

\begin{align*} \sigma_{i,j}: \ Covariance \ matrix \ between \ returns \ i \ and \ j \in \ commodities \\ \bar{p}_{c}: \ Mean \ historical \ price \ of \ commodity \ c \\ \rho: \ Required \ expected \ profit \end{align*}

(15)

Calculating Returns from Historical Price Data

\begin{align*} r^{lin} & = \frac{HP_{q+\Delta q} - HP_q}{HP_q} \ for \ c \in C \end{align*}

(16)

Monetary Exposure

\begin{align*} w_{c} = \bar{p}_{c} \cdot x_{c} \ \ \forall \ \ c \in C \end{align*}

(17)

Objective

\begin{align*} minimize \ MV = \sum_{i=1}^{n} \sum_{j=1}^{n} \sigma_{i,j} w_i w_j \\ s.t. \ \ Eprofit \geq \rho \end{align*}

(18)

Note

Why $w_c$ and not $x_c$ ? Markowitz’s $x_i$ is the money invested in asset $i$ , so $\sigma_{i,j} x_i x_j$ has units of currency squared. Here $x_c$ is physical production --- tonnes, cubic meters, and MWh --- and the covariance is of dimensionless price returns, so $\sum_{i,j} \sigma_{i,j} x_i x_j$ would add tonne $^2$ to m $^6$ to MWh $^2$ . Multiplying by the mean price $\bar{p}_c$ puts every commodity on a common monetary basis: $w_c$ is the revenue exposed to the price of commodity $c$ , in USD, and $MV$ comes out in USD $^2$ --- the right dimension for the variance of a quantity of money. Its square root is a standard deviation in USD, directly comparable to $MAD$ and $CVaR$ below.

Calculate the Covariance Matrix¶

# Calculate the Covariance Matrix

# Drop the scenario number column
price = df_hp.drop(columns=["q"])

# Calculate Returns
Returns = price.diff() / price.shift(1)

covar = Returns.cov()
print("Covariance Matrix")
print("-------------------")
print(covar)

Covariance Matrix
-------------------
          fre       sug       eth
fre  0.094993 -0.000350 -0.000756
sug -0.000350  0.000922  0.000390
eth -0.000756  0.000390  0.001578

The Units Check Bites¶

# Build MV over the raw production variables -- the dimensionally incoherent
# version -- and Pyomo refuses it before any solver is called.
m_bad = create_model()
m_bad.obj = pyo.Objective(
    expr=sum(
        m_bad.x[i] * covar.loc[i, j] * m_bad.x[j]
        for i in m_bad.commodities
        for j in m_bad.commodities
    ),
    sense=pyo.minimize,
)

# Pyomo also logs the offending expression; the exception below says it all.
logger = logging.getLogger("pyomo.util.check_units")
logger.setLevel(logging.CRITICAL)

try:
    assert_units_consistent(m_bad)
    print("No error -- this should not happen!")
except InconsistentUnitsError as e:
    print("InconsistentUnitsError:", e)

logger.setLevel(logging.WARNING)

InconsistentUnitsError: Error in units found in expression: 0.09499265505563526*x_energy[fre]*x_energy[fre] - 0.0003500722699504732*x_energy[fre]*x_mass[sug] - 0.0007555516171915655*x_energy[fre]*x_volume[eth] - 0.0003500722699504732*x_mass[sug]*x_energy[fre] + 0.0009217779801648573*x_mass[sug]*x_mass[sug] + 0.00039027914171650425*x_mass[sug]*x_volume[eth] - 0.0007555516171915655*x_volume[eth]*x_energy[fre] + 0.00039027914171650425*x_volume[eth]*x_mass[sug] + 0.0015784904740984824*x_volume[eth]*x_volume[eth]: megawatt_hour ** 2 not compatible with megawatt_hour * metric_ton.

# Reload the base Model
m = create_model()

# The epsilon-constraint: this portfolio must earn at least rho in expectation
add_required_return(m, rho_required)


# Define MV as the objective. This is written over the monetary exposures
# w_c = pbar_c * x_c, so every term is in USD^2.
def mean_value(m):
    return sum(
        m.exposure[i] * covar.loc[i, j] * m.exposure[j]
        for i in m.commodities
        for j in m.commodities
    )


m.obj = pyo.Objective(rule=mean_value, sense=pyo.minimize)

assert_units_consistent(m)

# Solve the Model
sol = pyo.SolverFactory("ipopt")

results = sol.solve(m)
assert pyo.check_optimal_termination(results), (
    f"Solve failed: status={results.solver.status}, "
    f"termination={results.solver.termination_condition}"
)

print("Results")
print("--------------------------------------")
print("MV", "%.4g" % pyo.value(m.obj), "USD^2")
print("sqrt(MV)", np.round(np.sqrt(pyo.value(m.obj)) / 1e6, 2), "M USD")
print("Expected Profit", np.round(pyo.value(m.EP), 2), "USD")
print("Required Return", np.round(rho_required, 2), "USD")
print("Ethanol Produced", np.round(pyo.value(m.x["eth"]), 2), "m3")
print("Sugar Produced", np.round(pyo.value(m.x["sug"]), 2), "tonne")
print("Electricity Produced", np.round(pyo.value(m.x["el-r"]), 2), "MWh")
print("Electricity to Free Market", np.round(pyo.value(m.x["fre"]), 2), "MWh")
print("Electricity to Regulated Market", np.round(pyo.value(m.x["reg"]), 2), "MWh")

Results
--------------------------------------
MV 2.904e+13 USD^2
sqrt(MV) 5.39 M USD
Expected Profit 61541577.98 USD
Required Return 61541577.98 USD
Ethanol Produced 163286.44 m3
Sugar Produced 175592.78 tonne
Electricity Produced 159000.0 MWh
Electricity to Free Market 14099.97 MWh
Electricity to Regulated Market 144900.03 MWh

Visualize Results¶

Profit Distribution¶

# Collect the profit distribution
profits = []
for i in m.price_time_obs:
    profits.append(pyo.value(m.profit[i]))

# Collect the minimum profit
min_prof = min(profits)

# Plot the Profit Distribution
fig, ax = plt.subplots(figsize=(6.4, 4))
plt.hist(np.array(profits) / 1e6, alpha=0.5)
plt.vlines(
    ymin=0,
    ymax=220,
    x=pyo.value(m.EP) / 1e6,
    label="Expected Profit",
    color="green",
    linestyles="solid",
)
plt.vlines(
    ymin=0,
    ymax=220,
    x=min_prof / 1e6,
    label=("Minimum Profit"),
    color="red",
    linestyles="dashed",
)
plt.xlabel("profit$_q$ M USD", fontsize=16, fontweight="bold")
plt.ylabel("Frequency", fontsize=16, fontweight="bold")
plt.legend(fontsize=14)
ax.xaxis.set_tick_params(labelsize=15)
ax.yaxis.set_tick_params(labelsize=15)
ax.tick_params(direction="in")
plt.show()

print("Results")
print("-------------------------------------")
print("Expected Profit:", np.round(pyo.value(m.EP) / 1e6, 2), "M USD")
print("Minimum Profit:", np.round(min_prof / 1e6, 2), "M USD")
difference = pyo.value(m.EP) - min_prof
print("Difference:", np.round(difference / 1e6, 2), "M USD")

# Collect results for conclusion
final_EP["MV"] = pyo.value(m.EP) / 1e6
final_min["MV"] = min_prof / 1e6
final_diff["MV"] = difference / 1e6
final_exposure["MV"] = {c: pyo.value(m.exposure[c]) for c in m.commodities}
final_profits["MV"] = np.array(profits)

Results
-------------------------------------
Expected Profit: 61.54 M USD
Minimum Profit: -32.36 M USD
Difference: 93.91 M USD

Product Distribution¶

fig, ax = plt.subplots(figsize=(6.4, 4))
plt.bar("Sugar", pyo.value(m.avg_sug_prof) / 1e6)
plt.bar("Ethanol", pyo.value(m.avg_eth_prof) / 1e6)
plt.bar("Electricity \n Free \n Market", pyo.value(m.avg_el_prof) / 1e6)
plt.bar("Electricity \n Regulated \n Market", pyo.value(m.x["reg"]) * 72.5 / 1e6)
plt.xlabel("Product", fontsize=16, fontweight="bold")
plt.ylabel("Revenue from \n Each Product (M USD)", fontsize=16, fontweight="bold")
ax.xaxis.set_tick_params(labelsize=15)
ax.yaxis.set_tick_params(labelsize=15)
ax.tick_params(direction="in")
plt.show()

Formulate the MAD Objective¶

\begin{align*} minimize \ \ MAD = E \left[ \mid \sum_{j=1}^{n} R_j w_j -E \left[ \sum_{j=1}^{n} R_j w_j \right] \mid \right] \\ s.t. \ \ Eprofit \geq \rho \end{align*}

(19)

# Reload the base Model
m = create_model()

# Ipopt cannot start this model from all zeros either, so solve the
# expected-profit problem first and keep its solution as the starting point.
initialize_from_max_profit(m)

# Create set of return indicies.  Returns.loc[0] is NaN -- the first week has
# no previous week to difference against -- so the deviations run over 1..N-1.
J = np.arange(1, m.N)
n_returns = len(J)

# Mean return of each commodity over those n_returns observations.  Dividing
# the sum by m.N instead of n_returns would use the wrong sample size.
mean_return = {c: Returns.loc[J, c].mean() for c in m.commodities}

# Create new variables. The deviations are monetary, like the exposures.
m.y_aux = pyo.Var(J, domain=pyo.NonNegativeReals, units=u.USD)
m.z_aux = pyo.Var(J, domain=pyo.NonNegativeReals, units=u.USD)


# Constrain y - z to be the abs portion
def aux_con(m, j):
    return m.y_aux[j] - m.z_aux[j] == sum(
        m.exposure[i] * (Returns.loc[j, i] - mean_return[i]) for i in m.commodities
    )


m.aux = pyo.Constraint(J, rule=aux_con)

# Initialize the deviations from the exposures we already have
for j in J:
    dev = sum(
        pyo.value(m.exposure[i]) * (Returns.loc[j, i] - mean_return[i])
        for i in m.commodities
    )
    m.y_aux[j].set_value(max(0.0, dev))
    m.z_aux[j].set_value(max(0.0, -dev))

# The epsilon-constraint: this portfolio must earn at least rho in expectation
add_required_return(m, rho_required)


# Define MAD as the objective
def mad(m):
    return (1 / n_returns) * sum(m.y_aux[j] + m.z_aux[j] for j in J)


m.obj = pyo.Objective(rule=mad, sense=pyo.minimize)

assert_units_consistent(m)

# Solve the Model
sol = pyo.SolverFactory("ipopt")

results = sol.solve(m)
assert pyo.check_optimal_termination(results), (
    f"Solve failed: status={results.solver.status}, "
    f"termination={results.solver.termination_condition}"
)

print("Results")
print("------------------------------")
print("MAD", np.round(pyo.value(m.obj) / 1e6, 2), "M USD")
print("Expected Profit", np.round(pyo.value(m.EP), 2), "USD")
print("Required Return", np.round(rho_required, 2), "USD")
print("Ethanol Produced", np.round(pyo.value(m.x["eth"]), 2), "m3")
print("Sugar Produced", np.round(pyo.value(m.x["sug"]), 2), "tonne")
print("Electricity Produced", np.round(pyo.value(m.x["el-r"]), 2), "MWh")
print("Electricity to Free Market", np.round(pyo.value(m.x["fre"]), 2), "MWh")
print("Electricity to Regulated Market", np.round(pyo.value(m.x["reg"]), 2), "MWh")

Results
------------------------------
MAD 3.61 M USD
Expected Profit 61541577.98 USD
Required Return 61541577.98 USD
Ethanol Produced 163217.88 m3
Sugar Produced 175705.63 tonne
Electricity Produced 159000.0 MWh
Electricity to Free Market 10856.52 MWh
Electricity to Regulated Market 148143.48 MWh

Visualize Results¶

Profit Distribution¶

# Collect the profit distribution
profits = []
for i in m.price_time_obs:
    profits.append(pyo.value(m.profit[i]))

# Collect the minimum profit
min_prof = min(profits)

# Plot the Profit Distribution
fig, ax = plt.subplots(figsize=(6.4, 4))
plt.hist(np.array(profits) / 1e6, alpha=0.5)
plt.vlines(
    ymin=0,
    ymax=220,
    x=pyo.value(m.EP) / 1e6,
    label="Expected Profit",
    color="green",
    linestyles="solid",
)
plt.vlines(
    ymin=0,
    ymax=220,
    x=min_prof / 1e6,
    label=("Minimum Profit"),
    color="red",
    linestyles="dashed",
)
plt.xlabel("profit$_q$ M USD", fontsize=16, fontweight="bold")
plt.ylabel("Frequency", fontsize=16, fontweight="bold")
plt.legend(fontsize=14)
ax.xaxis.set_tick_params(labelsize=15)
ax.yaxis.set_tick_params(labelsize=15)
ax.tick_params(direction="in")
plt.show()

print("Results")
print("-------------------------------------")
print("Expected Profit:", np.round(pyo.value(m.EP) / 1e6, 2), "M USD")
print("Minimum Profit:", np.round(min_prof / 1e6, 2), "M USD")
difference = pyo.value(m.EP) - min_prof
print("Difference:", np.round(difference / 1e6, 2), "M USD")

# Collect results for conclusion
final_EP["MAD"] = pyo.value(m.EP) / 1e6
final_min["MAD"] = min_prof / 1e6
final_diff["MAD"] = difference / 1e6
final_exposure["MAD"] = {c: pyo.value(m.exposure[c]) for c in m.commodities}
final_profits["MAD"] = np.array(profits)

Results
-------------------------------------
Expected Profit: 61.54 M USD
Minimum Profit: -32.17 M USD
Difference: 93.71 M USD

Product Distribution¶

fig, ax = plt.subplots(figsize=(6.4, 4))
plt.bar("Sugar", pyo.value(m.avg_sug_prof) / 1e6)
plt.bar("Ethanol", pyo.value(m.avg_eth_prof) / 1e6)
plt.bar("Electricity \n Free Market", pyo.value(m.avg_el_prof) / 1e6)
plt.bar("Electricity \n Regulated Market", pyo.value(m.x["reg"]) * 72.5 / 1e6)
plt.xlabel("Product", fontsize=16, fontweight="bold")
plt.ylabel("Revenue from \n Each Product (M USD)", fontsize=16, fontweight="bold")
ax.xaxis.set_tick_params(labelsize=15)
ax.yaxis.set_tick_params(labelsize=15)
ax.tick_params(direction="in")
plt.show()

Formulate the CVaR Objective¶

Additional Parameters

\begin{align*} \alpha: Confidence\ Level \end{align*}

(20)

Additional Variables

\begin{align*} \nu: VaR \\ shortfall_{q}: Difference \ between \ VaR \ and \ profit_{q} \ at \ price \ scenario \ q\\ \end{align*}

(21)

Additional Constraints

\begin{align*} shortfall_{q} \ \geq \ 0 \ \ \forall \ \ q \in Q \\ profit_{q} + shortfall_{q} - \nu \geq 0 \ \ \forall \ \ q \in Q \\ CVaR = \nu - \left[ \frac{1}{Q(1-\alpha)} \sum_{q \in Q} shortfall_{q} \right] \\ Eprofit \geq \rho \end{align*}

(22)

CVaR Objective

\begin{align*} \max{CVaR} \end{align*}

(23)

Note

This model needs to be initialized. The CVaR formulation above is a linear program and it is perfectly feasible, but if you hand it to Ipopt from the default all-zero starting point the solver reports

Converged to a locally infeasible point. Problem may be infeasible.

At $\text{profit}_q = \text{shortfall}_q = \nu = 0$ every one of the $2N$ CVaR inequalities is active simultaneously, so the linear independence constraint qualification (LICQ) fails there and Ipopt cannot compute a useful step; it enters its restoration phase immediately and never leaves. The IDAES Diagnostics Toolbox flags exactly this --- report_numerical_issues() refuses to run at all until the variables are given values.

The cell below therefore solves the expected-profit problem first (which Ipopt handles without difficulty) and uses that solution, plus an analytic guess for $\nu$ and the shortfalls, as the starting point. This is a good habit in general: scenario-based risk models almost always need to be initialized.

# Reload the base Model
m = create_model()

# Add parameters
alpha = 0.9  # confidence level

# ---------------------------------------------------------------------------
# STEP 1: INITIALIZE.
#
# Ipopt starts every uninitialized variable at zero, and for the CVaR model
# that point is badly degenerate: with profit_q = shortfall_q = nu = 0,
# *all* 2N CVaR inequalities (cvar2eq and cvar4eq, 1736 of them) are active at
# once, so LICQ cannot hold.  Ipopt drops into its restoration phase on the
# first iteration and never leaves it, reporting "Converged to a locally
# infeasible point."  The model is NOT infeasible -- it just needs a starting
# point that is strictly inside the CVaR inequalities.
#
# The expected-profit problem from the top of this notebook solves without any
# trouble, so we use it as the initialization: solve it first, then set alpha
# and the shortfalls analytically from the resulting profit distribution.
# ---------------------------------------------------------------------------
initialize_from_max_profit(m)

# STEP 2: build the CVaR model on top of that initialized point.

# Add variables
m.shortfall = pyo.Var(m.price_time_obs, units=u.USD)
m.nu = pyo.Var(units=u.USD)
m.CVaR = pyo.Var(units=u.USD)


# Add CVaR Constraints
def CVaR1(m):
    return m.CVaR == m.nu - (
        (1 / (m.N * (1 - alpha))) * sum(m.shortfall[q] for q in m.price_time_obs)
    )


m.cvar1eq = pyo.Constraint(rule=CVaR1)


def CVaR2(m, q):
    return m.shortfall[q] >= 0


m.cvar2eq = pyo.Constraint(m.price_time_obs, rule=CVaR2)


def CVaR4(m, q):
    return m.profit[q] + m.shortfall[q] - m.nu >= 0


m.cvar4eq = pyo.Constraint(m.price_time_obs, rule=CVaR4)

# Initialize the new variables.  At the optimum nu is the value-at-risk, so
# start it at the (1 - alpha) quantile of the profit distribution and set each
# shortfall to max(0, nu - profit_q), which is exactly what the two CVaR
# inequalities require.
init_profits = np.array([pyo.value(m.profit[q]) for q in m.price_time_obs])
nu_init = float(np.percentile(init_profits, 100 * (1 - alpha)))
m.nu.set_value(nu_init)
for q in m.price_time_obs:
    m.shortfall[q].set_value(max(0.0, nu_init - pyo.value(m.profit[q])))
m.CVaR.set_value(
    nu_init
    - sum(pyo.value(m.shortfall[q]) for q in m.price_time_obs)
    / (pyo.value(m.N) * (1 - alpha))
)


# The epsilon-constraint: this portfolio must earn at least rho in expectation
add_required_return(m, rho_required)


# OBJECTIVE
def obj_rule(m):
    return m.CVaR


m.obj = pyo.Objective(rule=obj_rule, sense=pyo.maximize)

assert_units_consistent(m)

# Solve the Model
results = sol.solve(m)
assert pyo.check_optimal_termination(results), (
    f"Solve failed: status={results.solver.status}, "
    f"termination={results.solver.termination_condition}"
)

print("Results")
print("------------------------------")
print("CVaR", np.round(pyo.value(m.CVaR), 2), "USD")
print("VaR (nu)", np.round(pyo.value(m.nu), 2), "USD")
print("Expected Profit", np.round(pyo.value(m.EP), 2), "USD")
print("Required Return", np.round(rho_required, 2), "USD")
print("Ethanol Produced", np.round(pyo.value(m.x["eth"]), 2), "m3")
print("Sugar Produced", np.round(pyo.value(m.x["sug"]), 2), "tonne")
print("Electricity Produced", np.round(pyo.value(m.x["el-r"]), 2), "MWh")
print("Electricity to Free Market", np.round(pyo.value(m.x["fre"]), 2), "MWh")
print("Electricity to Regulated Market", np.round(pyo.value(m.x["reg"]), 2), "MWh")

Results
------------------------------
CVaR 4423136.02 USD
VaR (nu) 15088239.22 USD
Expected Profit 63816117.69 USD
Required Return 61541577.98 USD
Ethanol Produced 155062.24 m3
Sugar Produced 189129.78 tonne
Electricity Produced 159000.0 MWh
Electricity to Free Market -0.0 MWh
Electricity to Regulated Market 159000.0 MWh

Visualize Results¶

Profit Distribution¶

# Collect the profit distribution
profits = []
for i in m.price_time_obs:
    profits.append(pyo.value(m.profit[i]))

# Collect the minimum profit
min_prof = min(profits)

# Plot the Profit Distribution
fig, ax = plt.subplots(figsize=(6.4, 4))
plt.hist(np.array(profits) / 1e6, alpha=0.5)
plt.vlines(
    ymin=0,
    ymax=220,
    x=pyo.value(m.EP) / 1e6,
    label="Expected Profit",
    color="green",
    linestyles="solid",
)
plt.vlines(
    ymin=0,
    ymax=220,
    x=min_prof / 1e6,
    label=("Minimum Profit"),
    color="red",
    linestyles="dashed",
)
plt.xlabel("profit$_q$ M USD", fontsize=16, fontweight="bold")
plt.ylabel("Frequency", fontsize=16, fontweight="bold")
plt.legend(fontsize=14)
ax.xaxis.set_tick_params(labelsize=15)
ax.yaxis.set_tick_params(labelsize=15)
ax.tick_params(direction="in")
plt.show()

print("Results")
print("-------------------------------------")
print("Expected Profit:", np.round(pyo.value(m.EP) / 1e6, 2), "M USD")
print("Minimum Profit:", np.round(min_prof / 1e6, 2), "M USD")
difference = pyo.value(m.EP) - min_prof
print("Difference:", np.round(difference / 1e6, 2), "M USD")

# Collect results for conclusion
final_EP["CVaR"] = pyo.value(m.EP) / 1e6
final_min["CVaR"] = min_prof / 1e6
final_diff["CVaR"] = difference / 1e6
final_exposure["CVaR"] = {c: pyo.value(m.exposure[c]) for c in m.commodities}
final_profits["CVaR"] = np.array(profits)

Results
-------------------------------------
Expected Profit: 63.82 M USD
Minimum Profit: -29.94 M USD
Difference: 93.76 M USD

Product Distribution¶

fig, ax = plt.subplots(figsize=(6.4, 4))
plt.bar("Sugar", pyo.value(m.avg_sug_prof) / 1e6)
plt.bar("Ethanol", pyo.value(m.avg_eth_prof) / 1e6)
plt.bar("Electricity \n Free \n Market", pyo.value(m.avg_el_prof) / 1e6)
plt.bar("Electricity \n Regulated \n Market", pyo.value(m.x["reg"]) * 72.5 / 1e6)
plt.ylabel("Revenue from \n Each Product (M USD)", fontsize=16, fontweight="bold")
plt.xlabel("Product", fontsize=16, fontweight="bold")
ax.xaxis.set_tick_params(labelsize=15)
ax.yaxis.set_tick_params(labelsize=15)
ax.tick_params(direction="in")
plt.show()

Trace the Frontier: An $\epsilon$ -Constraint Sweep¶

One value of $\rho$ gives one Pareto-optimal portfolio. Solving the same three models over a range of $\rho$ traces the whole frontier for each measure.

Two anchors bound the sweep. The top is the maximum expected profit, which is the same for all three measures. The bottom is each measure’s own minimum-risk portfolio, solved with no required-return constraint at all, because no value of $\rho$ can push a measure below its own minimum. Those two anchors are what a payoff table provides in the multi-objective literature, and they are why the three frontiers below start at three different heights.

The three risk measures do not share units, so each frontier is plotted against its own measure in M USD: $\sqrt{MV}$ , which is a standard deviation; $MAD$ ; and the CVaR of the loss, which is $-CVaR$ of the profit. All three are then oriented the same way, smaller being less risk, but they are still three different quantities. Compare the shape of the three curves and the return each one starts at, not their left-to-right positions.

# The three models again, wrapped as functions of rho so they can be called in
# a loop. Nothing here is new: these are the same three models solved above.
def build_mv(rho=None):
    """Build the mean-variance model, optionally with a required return.

    Arguments:
        rho: required expected profit [USD], or None for no constraint

    Returns:
        a Pyomo ConcreteModel, not yet solved
    """
    m = create_model()
    if rho is not None:
        add_required_return(m, rho)
    m.obj = pyo.Objective(
        expr=sum(
            m.exposure[i] * covar.loc[i, j] * m.exposure[j]
            for i in m.commodities
            for j in m.commodities
        ),
        sense=pyo.minimize,
    )
    return m


def build_mad(rho=None):
    """Build the mean-absolute-deviation model, optionally with a required return.

    Arguments:
        rho: required expected profit [USD], or None for no constraint

    Returns:
        a Pyomo ConcreteModel, not yet solved
    """
    m = create_model()
    initialize_from_max_profit(m)

    m.y_aux = pyo.Var(J, domain=pyo.NonNegativeReals, units=u.USD)
    m.z_aux = pyo.Var(J, domain=pyo.NonNegativeReals, units=u.USD)

    def aux_con(b, j):
        return b.y_aux[j] - b.z_aux[j] == sum(
            b.exposure[i] * (Returns.loc[j, i] - mean_return[i]) for i in b.commodities
        )

    m.aux = pyo.Constraint(J, rule=aux_con)

    for j in J:
        dev = sum(
            pyo.value(m.exposure[i]) * (Returns.loc[j, i] - mean_return[i])
            for i in m.commodities
        )
        m.y_aux[j].set_value(max(0.0, dev))
        m.z_aux[j].set_value(max(0.0, -dev))

    if rho is not None:
        add_required_return(m, rho)

    m.obj = pyo.Objective(
        expr=(1 / n_returns) * sum(m.y_aux[j] + m.z_aux[j] for j in J),
        sense=pyo.minimize,
    )
    return m


def build_cvar(rho=None):
    """Build the CVaR model, optionally with a required return.

    Arguments:
        rho: required expected profit [USD], or None for no constraint

    Returns:
        a Pyomo ConcreteModel, not yet solved
    """
    m = create_model()
    initialize_from_max_profit(m)

    m.shortfall = pyo.Var(m.price_time_obs, units=u.USD)
    m.nu = pyo.Var(units=u.USD)
    m.CVaR = pyo.Var(units=u.USD)

    m.cvar1eq = pyo.Constraint(
        expr=m.CVaR
        == m.nu
        - (1 / (m.N * (1 - alpha))) * sum(m.shortfall[q] for q in m.price_time_obs)
    )

    def cvar2(b, q):
        return b.shortfall[q] >= 0

    m.cvar2eq = pyo.Constraint(m.price_time_obs, rule=cvar2)

    def cvar4(b, q):
        return b.profit[q] + b.shortfall[q] - b.nu >= 0

    m.cvar4eq = pyo.Constraint(m.price_time_obs, rule=cvar4)

    # Same analytic initialization as the CVaR cell above
    init_profits = np.array([pyo.value(m.profit[q]) for q in m.price_time_obs])
    nu_init = float(np.percentile(init_profits, 100 * (1 - alpha)))
    m.nu.set_value(nu_init)
    for q in m.price_time_obs:
        m.shortfall[q].set_value(max(0.0, nu_init - pyo.value(m.profit[q])))
    m.CVaR.set_value(
        nu_init
        - sum(pyo.value(m.shortfall[q]) for q in m.price_time_obs)
        / (pyo.value(m.N) * (1 - alpha))
    )

    if rho is not None:
        add_required_return(m, rho)

    m.obj = pyo.Objective(expr=m.CVaR, sense=pyo.maximize)
    return m

# Each measure reports its risk in its own units. Put all three in M USD and
# orient them so that a smaller number always means less risk.
builders = {"MV": build_mv, "MAD": build_mad, "CVaR": build_cvar}
risk_in_musd = {
    "MV": lambda m: np.sqrt(pyo.value(m.obj)) / 1e6,  # sqrt(MV) [M USD]
    "MAD": lambda m: pyo.value(m.obj) / 1e6,  # MAD [M USD]
    "CVaR": lambda m: -pyo.value(m.CVaR) / 1e6,  # CVaR of the loss [M USD]
}

n_points = 9
frontier = {}

for name, build in builders.items():
    # Lower anchor: the minimum-risk portfolio, with no required return imposed
    m_anchor = build()
    results = sol.solve(m_anchor)
    assert pyo.check_optimal_termination(results), (
        f"{name} minimum-risk solve failed: "
        f"termination={results.solver.termination_condition}"
    )
    risks = [risk_in_musd[name](m_anchor)]
    returns = [pyo.value(m_anchor.EP) / 1e6]

    # Sweep rho from that anchor up to the maximum attainable expected profit
    for rho in np.linspace(pyo.value(m_anchor.EP), EP_max, n_points)[1:]:
        m_rho = build(rho)
        results = sol.solve(m_rho)
        assert pyo.check_optimal_termination(results), (
            f"{name} sweep failed at rho = {rho / 1e6:.3f} M USD: "
            f"termination={results.solver.termination_condition}"
        )
        risks.append(risk_in_musd[name](m_rho))
        returns.append(pyo.value(m_rho.EP) / 1e6)

    frontier[name] = (risks, returns)
    print(
        f"{name:5s} risk {risks[0]:8.4f} -> {risks[-1]:8.4f} M USD   "
        f"expected profit {returns[0]:7.3f} -> {returns[-1]:7.3f} M USD"
    )

MV    risk   5.3870 ->   6.4439 M USD   expected profit  60.157 ->  64.781 M USD

MAD   risk   3.5831 ->   4.5157 M USD   expected profit  53.349 ->  64.781 M USD

CVaR  risk  -4.4231 ->   0.1173 M USD   expected profit  63.816 ->  64.781 M USD

# Colour AND linestyle AND marker: the handouts get printed in black and white,
# so no distinction here is carried by colour alone.
fig, ax = plt.subplots(figsize=(6.4, 4.5))
ax.plot(
    frontier["MV"][0],
    frontier["MV"][1],
    color="tab:blue",
    linestyle="-",
    marker="o",
    linewidth=2,
    label="MV: $\\sqrt{MV}$",
)
ax.plot(
    frontier["MAD"][0],
    frontier["MAD"][1],
    color="tab:orange",
    linestyle="--",
    marker="s",
    linewidth=2,
    label="MAD",
)
ax.plot(
    frontier["CVaR"][0],
    frontier["CVaR"][1],
    color="tab:green",
    linestyle=":",
    marker="^",
    linewidth=2,
    label="CVaR of loss",
)
ax.axhline(
    EP_max / 1e6,
    color="black",
    linestyle="-.",
    linewidth=1,
    label="Maximum expected profit",
)
ax.set_xlabel("Risk measure (M USD)", fontsize=16, fontweight="bold")
ax.set_ylabel("Expected profit (M USD)", fontsize=16, fontweight="bold")
ax.legend(fontsize=12)
ax.xaxis.set_tick_params(labelsize=15)
ax.yaxis.set_tick_params(labelsize=15)
ax.tick_params(direction="in")
plt.show()

The MV and MAD frontiers are nearly vertical over most of their length: the mill can be pushed to earn several million USD more per year for around one percent more risk. Both bend sharply only at the very top. That bend is where the factory’s sugar capacity runs out. Up to it, the extra return comes from shifting juice toward sugar, whose price is far more stable than free-market electricity; past it, the only remaining way to earn more is to sell electricity into the volatile free market.

The CVaR frontier is the opposite shape, and it is much shorter. Minimizing the CVaR of the loss already earns 63.8 M USD, only 1.5% below the maximum, so there is very little return left to demand -- but across that narrow range the average of the worst 10% of years deteriorates by 4.5 M USD. CVaR is a one-sided measure: it averages the worst outcomes and is indifferent to how good the good years are, so it does not charge for the upside variance that MV and MAD both penalize, and in exchange it reacts sharply to the last of the return.

All three curves end at the same top-right corner, which is the maximum-profit portfolio. The vertical extent of the three curves is comparable -- that is expected profit in M USD on a shared axis -- but the horizontal extent is not, because each curve is drawn against its own measure.

Conclusion¶

Final Results Comparison¶

k = 0
for i in final_EP.keys():
    if k < 0.5:
        plt.scatter(i, final_EP[i], color="green", marker="o", label="Expected Profit")
        plt.scatter(i, final_min[i], color="red", marker="s", label="Minimum Profit")
        plt.bar(
            i,
            final_diff[i],
            color="blue",
            alpha=0.5,
            bottom=final_min[i],
            hatch="//",
            label="Potential Loss",
        )
        plt.text(k - 0.2, 20, np.round(final_diff[i], 2), fontsize=14)
    else:
        plt.scatter(i, final_EP[i], color="green", marker="o")
        plt.scatter(i, final_min[i], color="red", marker="s")
        plt.bar(
            i, final_diff[i], color="blue", alpha=0.5, bottom=final_min[i], hatch="//"
        )
        plt.text(k - 0.2, 20, np.round(final_diff[i], 2), fontsize=14)
    k = k + 1

plt.axhline(
    rho_required / 1e6,
    color="black",
    linestyle="dotted",
    linewidth=2,
    label="Required return $\\rho$",
)

# `ax` from the previous cell refers to a figure that is already closed
ax = plt.gca()
ax.xaxis.set_tick_params(labelsize=15)
ax.yaxis.set_tick_params(labelsize=15)
ax.tick_params(direction="in")
plt.legend(bbox_to_anchor=(1, 1), fontsize=14)
plt.xlabel("Risk Measure", fontsize=16, fontweight="bold")
plt.ylabel("Profit Measure (M USD)", fontsize=16, fontweight="bold")
plt.show()

Every Portfolio, Measured Every Way¶

The figure above compares the four portfolios on the same two profit statistics. A sharper question is how each portfolio scores under the other measures’ definitions of risk, since a portfolio that is best under MV is not automatically good under CVaR. The table below evaluates all three risk measures on all four portfolios.

def portfolio_risk(exposure, profits, alpha):
    """Evaluate all three risk measures on one portfolio.

    Arguments:
        exposure: dict of monetary exposure w_c by commodity [USD]
        profits: array of profit_q over the price scenarios [USD]
        alpha: CVaR confidence level

    Returns:
        (sqrt(MV), MAD, CVaR of the loss), all in USD
    """
    w = np.array([exposure[c] for c in covar.columns])
    sigma = np.sqrt(w @ covar.to_numpy() @ w)

    # Deviation of the portfolio return from its mean, in USD, per observation
    R = Returns[list(covar.columns)].iloc[1:]
    mad = np.abs((R - R.mean()).to_numpy() @ w).mean()

    # Empirical CVaR of the loss = min over nu of nu + E[(loss - nu)+]/(1 - alpha).
    # The minimizer is always one of the observed losses, so check them all.
    loss = -np.asarray(profits)
    cvar = min(nu + np.maximum(loss - nu, 0.0).mean() / (1 - alpha) for nu in loss)

    return sigma, mad, cvar


rows = []
for name in final_EP:
    sigma, mad, cvar = portfolio_risk(final_exposure[name], final_profits[name], alpha)
    rows.append(
        {
            "Objective": name,
            "E[profit]": final_EP[name],
            "sqrt(MV)": sigma / 1e6,
            "MAD": mad / 1e6,
            "CVaR of loss": cvar / 1e6,
            "Min profit": final_min[name],
        }
    )

summary = pd.DataFrame(rows).set_index("Objective").round(4)
print("All values in M USD. Smaller is less risk in the three risk columns.")
print("---------------------------------------------------------------------")
print(summary.to_string())

All values in M USD. Smaller is less risk in the three risk columns.
---------------------------------------------------------------------
           E[profit]  sqrt(MV)     MAD  CVaR of loss  Min profit
Objective                                                       
No Risk      64.7806    6.4439  4.5157        0.1173    -39.0034
MV           61.5416    5.3888  3.6118       -3.1634    -32.3643
MAD          61.5416    5.3894  3.6085       -3.2252    -32.1660
CVaR         63.8161    5.4108  3.6312       -4.4231    -29.9420

Every portfolio is best under the measure it was optimized for and worse under the others: the diagonal wins in all three risk columns. That is what makes the three measures genuinely different rather than three names for the same thing, and it is only visible because all three were solved subject to the same required return.

Notice how narrow the gap between the MV and MAD portfolios is -- they agree to three decimal places on $\sqrt{MV}$ -- because the same lever dominates both. Free-market electricity carries sixty to a hundred times the price variance of sugar or ethanol, so nearly all of the available risk reduction comes from selling less of it, whichever measure is asked.

The CVaR row is the interesting one. Its required-return constraint is inactive: minimizing the CVaR of the loss already earns more than $\rho$ demands, so it is the one portfolio here that was not pushed down to the return floor. Compared with the no-risk portfolio it gives up about 1.5% of the expected profit and, in exchange, turns the average of the worst 10% of years from a 0.1 M USD loss into a 4.4 M USD profit.

Discussion Questions¶

References¶

[1] Gunjan, A., & Bhattacharyya, S. (2022). A brief review of portfolio optimization techniques. Artificial Intelligence Review. https://doi.org/10.1007/s10462-022-10273-7

[2] Markowitz, H. (1959). Portfolio selection, efficient diversification of investments. Wiley, New York

[3] Konno, H., & Koshizuka, T. (2005). Mean-absolute deviation model. IIE Transactions (Institute of Industrial Engineers), 37(10), 893–900. https://doi.org/10.1080/07408170591007786

[4] Mutran, V. M., Ribeiro, C. O., Nascimento, C. A. O., & Chachuat, B. (2020). Risk-conscious optimization model to support bioenergy investments in the Brazilian sugarcane industry. Applied Energy, 258. https://doi.org/10.1016/j.apenergy.2019.113978

[5] Li, L., Lin, J., Wu, N., Xie, S., Meng, C., Zheng, Y., Wang, X., & Zhao, Y. (2022). Review and outlook on the international renewable energy development. In Energy and Built Environment (Vol. 3, Issue 2, pp. 139–157). KeAi Communications Co. https://doi.org/10.1016/j.enbenv.2020.12.002

[6] Oliveira, S. M. de, Ribeiro, C. de O., & Cicogna, M. P. V. (2018). Uncertainty effects on production mix and on hedging decisions: The case of Brazilian ethanol and sugar. Energy Economics, 70, 516–524. https://doi.org/10.1016/j.eneco.2018.01.025

Introduction¶

Common Risk Measures for Portfolio Optimization¶

Mean-Variance (MV) [1,2]¶

Literature Examples Using MV¶

Mean-Absolute Deviation (MAD) [1,3]¶

Graphical Representation [3]¶

Literature Examples Using MAD¶

Minimax (MM) [1]¶

Literature Examples Using MM¶

Value-at-Risk (VaR) [1,4]¶

Literature Examples Using VaR¶

Conditional Value-at-Risk (CVaR) [1,4]¶

Graphical Representation [4]¶

Literature Examples Using CVaR¶

Motivating Example: Product Portfolio Optimization in the Brazilian Sugarcane Industry¶

Problem Definition¶

Process Flow Diagram¶

Process Model¶

Degree of Freedom (DOF) Analysis¶

Visualize Input Data¶

Define Base Model in Pyomo¶

Maximize Profit Only (No Risk Measure)¶

Visualize Results¶

Profit Distribution¶

Product Distribution¶

Solve Model with Risk Measures¶

Two Objectives, Not One¶

Formulate the MV Objective¶

Calculate the Covariance Matrix¶

The Units Check Bites¶

Visualize Results¶

Profit Distribution¶

Product Distribution¶

Formulate the MAD Objective¶

Visualize Results¶

Profit Distribution¶

Product Distribution¶

Formulate the CVaR Objective¶

Visualize Results¶

Profit Distribution¶

Product Distribution¶

Trace the Frontier: An ϵ\epsilonϵ-Constraint Sweep¶

Conclusion¶

Final Results Comparison¶

Every Portfolio, Measured Every Way¶

Discussion Questions¶

References¶

Trace the Frontier: An $\epsilon$ -Constraint Sweep¶