Portfolio Optimization with Python and R

Libraries and Modules

We need the following Python libraries and modules:

import numpy as np
import pandas as pd
import yfinance as yf
import scipy.optimize as sco
import scipy.interpolate as sci

Data

         For our sample, we select 12 different assets for the analysis with exposure to all of the GICS sectors—– energy, materials, industrials, utilities, healthcare, financials, consumer discretionary, consumer staples, information technology, communication services, and real estate. The start date is 2011-09-13 and the end date is by default today’s date, which is 2023-09-30.

# Create a list of symbols
symbols = [
  "XOM", "SHW", "JPM", "AEP", "UNH", "AMZN", 
  "KO", "BA", "AMT", "DD", "TSN", "SLG"
]
# Create data frame
asset_data = yf.download(
  tickers = ' '.join(symbols), 
  start = '2013-09-13'
)['Adj Close']
# Set column indices

[                       0%%                      ]
[********              17%%                      ]  2 of 12 completed
[************          25%%                      ]  3 of 12 completed
[****************      33%%                      ]  4 of 12 completed
[********************  42%%                      ]  5 of 12 completed
[**********************50%%                      ]  6 of 12 completed
[**********************58%%**                    ]  7 of 12 completed
[**********************67%%******                ]  8 of 12 completed
[**********************75%%**********            ]  9 of 12 completed
[**********************83%%**************        ]  10 of 12 completed
[**********************92%%******************    ]  11 of 12 completed
[*********************100%%**********************]  12 of 12 completed

asset_data.columns = symbols
# Examine the first 5 rows
asset_data.head(n = 5)

                  XOM        SHW      JPM  ...         DD        TSN        SLG
Date                                       ...                                 
2013-09-13  29.889591  60.253777  14.8960  ...  25.158703  63.897133  57.630928
2013-09-16  30.029844  60.711533  14.8030  ...  25.200283  64.446182  57.806953
2013-09-17  30.296333  60.425430  15.2085  ...  25.175337  63.751278  57.982964
2013-09-18  31.425438  61.921352  15.6015  ...  25.308403  62.661716  58.400208
2013-09-19  31.327253  61.388027  15.6030  ...  25.300085  60.774311  58.204639

[5 rows x 12 columns]

Next, we find the simple daily returns for each of the 12 assets using the pct_change() method, since our data object is a Pandas DataFrame. We use simple returns since they have the property of being asset-additive, which is necessary since we need to compute portfolios returns:

# Compute daily simple returns
daily_returns = (
  asset_data.pct_change()
            .dropna(
              # Drop the first row since we have NaN's
              # The first date 2011-09-13 does not have a value since it is our cut-off date
              axis = 0,
              how = 'any',
              inplace = False
              )
)
# Examine the last 5 rows
daily_returns.tail(n = 5)

                 XOM       SHW       JPM  ...        DD       TSN       SLG
Date                                      ...                              
2023-09-25 -0.006189 -0.002948  0.016651  ... -0.009781  0.008259  0.011223
2023-09-26 -0.029741 -0.016956 -0.040299  ...  0.012841 -0.009465  0.001549
2023-09-27 -0.006288 -0.019766  0.000000  ... -0.022625 -0.003403  0.032557
2023-09-28 -0.015687  0.016469  0.000000  ... -0.003592  0.012646 -0.006073
2023-09-29  0.007366  0.013122  0.009049  ...  0.011216 -0.011586 -0.015820

[5 rows x 12 columns]

         The simple daily returns may be visualized using line charts, density plots, and histograms, which are covered in my other post on visualizing asset data. Even though the visualizations in that post use the ggplot2 package in R, the plotnine package, or any other Python graphics librarires, can be employed to produce them in Python. For now, let us annualize the daily returns over the 10-year period from 2011-9-13 to 2023-09-30. We assume the number of trading days in a year is computed as follows:

\[\begin{align*} 365.25 \text{(days on average per year)} \times \frac{5}{7} \text{(proportion work days per week)} \\ - 6 \text{(weekday holidays)} - 3\times\frac{5}{7} \text{(fixed date holidays)} = 252.75 \approx 253 \end{align*}\]

daily_returns.mean() * 253

XOM     0.113546
SHW     0.130445
JPM     0.269720
AEP     0.147516
UNH     0.099439
AMZN    0.165946
KO      0.084762
BA      0.191212
AMT     0.026794
DD      0.111788
TSN     0.239387
SLG     0.109048
dtype: float64

         As can be seen, there are significant differences in the annualized performances between these assets. Amazon dominated the period with an annualized rate of return of \(32.2\%\). Exxon Mobil Corporation, on the other hand, is at the bottom of the rankings, recording an annualized rate of return of only \(4.2\%\) over the period. The annualized variance-covariance matrix of the returns can be computed using built-in pandas method cov:

daily_returns.cov() * 253

           XOM       SHW       JPM  ...        DD       TSN       SLG
XOM   0.042263  0.026698  0.011350  ...  0.014352  0.018433  0.014551
SHW   0.026698  0.059889  0.026605  ...  0.017013  0.024488  0.017036
JPM   0.011350  0.026605  0.110491  ...  0.015029  0.024801  0.018723
AEP   0.021211  0.029408  0.038162  ...  0.036746  0.036390  0.052296
UNH   0.016328  0.023618  0.030571  ...  0.027182  0.029955  0.043092
AMZN  0.014389  0.023525  0.027258  ...  0.026251  0.032357  0.042200
KO    0.020587  0.022202  0.013947  ...  0.017513  0.019160  0.019685
BA    0.017234  0.027274  0.029680  ...  0.023044  0.028145  0.020734
AMT   0.028434  0.037496  0.028200  ...  0.037967  0.031699  0.047775
DD    0.014352  0.017013  0.015029  ...  0.083850  0.019883  0.025544
TSN   0.018433  0.024488  0.024801  ...  0.019883  0.065081  0.025917
SLG   0.014551  0.017036  0.018723  ...  0.025544  0.025917  0.075370

[12 rows x 12 columns]

The variance-covariance matrix of the returns will be needed to compute the variance of the portfolio returns.

The Optimization Problem

         The portfolio optimization problem, therefore, given a universe of assets and their characteristics, deals with a method to spread the capital between them in a way that maximizes the return of the portfolio per unit of risk taken. There is no unique solution for this problem, but a set of solutions, which together define what is called an efficient frontier— the portfolios whose returns cannot be improved without increasing risk, or the portfolios where risk cannot be reduced without reducing returns as well. The Markowitz model for the solution of the portfolio optimization problem has a twin objective of maximizing return and minimizing risk, built on the Mean-Variance framework of asset returns and holding the basic constraints, which reduces to the following:

Minimize Risk given Levels of Return

\[\begin{align*} \min_{\vec{w}} \hspace{5mm} \sqrt{\vec{w}^{T} \hat{\Sigma} \vec{w}} \end{align*}\]

subject to

\[\begin{align*} &\vec{w}^{T} \hat{\mu}=\bar{r}_{P} \\ &\vec{w}^{T} \vec{1} = 1 \hspace{5mm} (\text{Full investment}) \\ &\vec{0} \le \vec{w} \le \vec{1} \hspace{5mm} (\text{Long only}) \end{align*}\]

Maximize Return given Levels of Risk

\[\begin{align*} \max _{\vec{w}} \hspace{5mm} \vec{w}^{T} \hat{\mu} \end{align*}\]

subject to

\[\begin{align*} &\vec{w}^{T} \hat{\Sigma} \vec{w}=\bar{\sigma}_{P} \\ &\vec{w}^{T} \vec{1} = 1 \hspace{5mm} (\text{Full investment}) \\ &\vec{0} \le \vec{w} \le \vec{1} \hspace{5mm} (\text{Long only}) \end{align*}\]

         In absence of other constraints, the above model is loosely referred to as the “unconstrained” portfolio optimization model. Solving the mathematical model yields a set of optimal weights representing a set of optimal portfolios. The solution set to these two problems is a hyperbola that depicts the efficient frontier in the \(\mu-\sigma\) -diagram.

Monte Carlo Simulation

         The first task is to simulate a random set of portfolios to visualize the risk-return profiles of our given set of assets. To carry out the Monte Carlo simulation, we define two functions that both take as inputs a vector of asset weights and output the expected portfolio return and standard deviation:

Returns

# Function for computing portfolio return
def portfolio_returns(weights):
    return (np.sum(daily_returns.mean() * weights)) * 253

Standard Deviation

# Function for computing standard deviation of portfolio returns
def portfolio_sd(weights):
    return np.sqrt(np.transpose(weights) @ (daily_returns.cov() * 253) @ weights)

         Next, we use a for loop to simulate random vectors of asset weights, computing the expected portfolio return and standard deviation for each permutation of random weights. Again, we ensure that each random weight vector is subject to the long-positions-only and full-investment constraints.

Monte Carlo Simulation

         The empty containers we instantiate are lists; they are mutable and so growing them will not be memory inefficient.

# instantiate empty list containers for returns and sd
list_portfolio_returns = []
list_portfolio_sd = []
# For loop to simulate 5000 random weight vectors (numpy array objects)
for p in range(5000):
  # Return random floats in the half-open interval [0.0, 1.0)
  weights = np.random.random(size = len(symbols)) 
  # Normalize to unity
  # The /= operator divides the array by the sum of the array and rebinds "weights" to the new object
  weights /= np.sum(weights) 
  # Lists are mutable so growing will not be memory inefficient
  list_portfolio_returns.append(portfolio_returns(weights))
  list_portfolio_sd.append(portfolio_sd(weights))
  # Convert list to numpy arrays
  port_returns = np.array(object = list_portfolio_returns)
  port_sd = np.array(object = list_portfolio_sd)

Let us examine the simulation results. In particular, the highest and the lowest expected portfolio returns are as follows:

# Max expected return
round(max(port_returns), 4)
# Min expected return

0.1835

round(min(port_returns), 4)

0.102

On the other hand, the highest and lowest volatility measures are recorded as:

# Max sd
round(max(port_sd), 4)
# Min sd

0.2357

round(min(port_sd), 4)

0.1636

         We may also visualize the expected returns and standard deviations on a \(\mu-\sigma\) trade-off diagram. For this task, I will leverage R’s graphics engine and the plotly graphics library. The reticulate package in R allows for relatively seamless transition between Python and R. Fortunately, the NumPy arrays created in Python can be accessed as R vector objects; this makes plotting in R using Python objects simple:

# Plot the sub-optimal portfolios
plotly::plot_ly(
  x = py$port_sd, y = py$port_returns, color = (py$port_returns / py$port_sd),
  mode = "markers", type = "scattergl", showlegend = FALSE,
  marker = list(size = 5, opacity = 0.7)
) %>%
  plotly::layout(
    title = "Mean-Standard Deviation Diagram",
    yaxis = list(title = "Expected Portfolio Return (Annualized)", tickformat = ".2%"),
    xaxis = list(title = "Portoflio Standard Deviation (Annualized)", tickformat = ".2%")
  ) %>% 
  plotly::colorbar(title = "Sharpe Ratio")

         Each point in the diagram above represents a permutation of expected-return-standard-deviation value pair. The points are color coded such that the magnitudes of the Sharpe ratios, defined as \(SR ≡ \frac{\mu_{P} – r_{f}}{\sigma_{P}}\), can be readily observed for each expected-return-standard-deviation pairing. For simplicity, we assume that \(r_{f} ≡ 0\). It could be argued that the assumption here is restrictive, so I explored using a different risk-free rate in my previous post.

The Optimal Portfolios

         Solving the optimization problem defined earlier provides us with a set of optimal portfolios given the characteristics of our assets. There are two important portfolios that we may be interested in constructing— the minimum variance portfolio and the maximal Sharpe ratio portfolio. In the case of the maximal Sharpe ratio portfolio, the objective function we wish to maximize is our user-defined Sharpe ratio function. The constraint is that all weights sum up to 1. We also specify that the weights are bound between 0 and 1. In order to use the minimization function from the SciPy library, we need to transform the maximization problem into one of minimization. In other words, the negative value of the Sharpe ratio is minimized to find the maximum value; the optimal portfolio composition is therefore the array of weights that yields that maximum value of the Sharpe ratio.

# User defined Sharpe ratio function
# Negative sign to compute the negative value of Sharpe ratio
def sharpe_fun(weights):
  return - (portfolio_returns(weights) / portfolio_sd(weights))

We will use dictionaries inside of a tuple to represent the constraints:

# We use an anonymous lambda function
constraints = ({'type': 'eq', 'fun': lambda x: np.sum(x) - 1})

Next, the bound values for the weights:

# This creates 12 tuples of (0, 1), all of which exist within a container tuple
# We essentially create a sequence of (min, max) pairs
bounds = tuple(
  (0, 1) for w in weights
)

We also need to supply a starting list of weights, which essentially functions as an initial guess. For our purposes, this will be an equal weight array:

# Repeat the list with the value (1 / 12) 12 times, and convert list to array
equal_weights = np.array(
  [1 / len(symbols)] * len(symbols)
)

We will use the scipy.optimize.minimize function and the Sequential Least Squares Programming (SLSQP) method for the minimization:

# Minimization results
max_sharpe_results = sco.minimize(
  # Objective function
  fun = sharpe_fun, 
  # Initial guess, which is the equal weight array
  x0 = equal_weights, 
  method = 'SLSQP',
  bounds = bounds, 
  constraints = constraints
)

The class of the optimization results is scipy.optimize.optimize.OptimizeResult, which contains many objects. The object of interest to us is the weight composition array, which we employ to construct the maximal Sharpe ratio portfolio:

# Extract the weight composition array
max_sharpe_results["x"]

array([1.38645164e-01, 0.00000000e+00, 2.67055226e-01, 0.00000000e+00,
       7.00394603e-17, 1.18973033e-02, 3.52636757e-17, 1.51940155e-01,
       1.83257272e-16, 1.35163801e-02, 4.16945772e-01, 0.00000000e+00])

# Expected return
max_sharpe_port_return = portfolio_returns(max_sharpe_results["x"])
round(max_sharpe_port_return, 4)

0.2201

# Standard deviation
max_sharpe_port_sd = portfolio_sd(max_sharpe_results["x"])
round(max_sharpe_port_sd, 4)

0.1947

# Sharpe ratio
max_sharpe_port_sharpe = max_sharpe_port_return / max_sharpe_port_sd
round(max_sharpe_port_sharpe, 4)

1.1304

# Minimize sd
min_sd_results = sco.minimize(
  # Objective function
  fun = portfolio_sd, 
  # Initial guess, which is the equal weight array
  x0 = equal_weights, 
  method = 'SLSQP',
  bounds = bounds, 
  constraints = constraints
)

# Expected return
min_sd_port_return = portfolio_returns(min_sd_results["x"])
round(min_sd_port_return, 4)

0.1288

# Standard deviation
min_sd_port_sd = portfolio_sd(min_sd_results["x"])
round(min_sd_port_sd, 4)

0.154

# Sharpe ratio
min_sd_port_sharpe = min_sd_port_return / min_sd_port_sd
round(min_sd_port_sharpe, 4)

0.8365

Efficient Frontier

         As an investor, one is generally interested in the maximum return given a fixed risk level or the minimum risk given a fixed return expectation. As mentioned in the earlier section, the set of optimal portfolios— whose expected portfolio returns for a defined level of risk cannot be surpassed by any other portfolio— depicts the so-called efficient frontier. The Python implementation is to fix a target return level and, for each such level, minimize the volatility value. For the optimization, we essentially “fit” the twin-objective described earlier into an optimization problem that can be solved using quadratic programming. (The objective function is the portfolio standard deviation formula, which is a quadratic function) Therefore, the two linear constraints are the target return (a linear function) and that the weights must sum to 1 (another linear function). We will again use dictionaries inside of a tuple to represent the constraints:

# We use anonymous lambda functions
# The argument x will be the weights
constraints = (
  {'type': 'eq', 'fun': lambda x: portfolio_returns(x) - target}, 
  {'type': 'eq', 'fun': lambda x: np.sum(x) - 1}
)

         The full-investment and long-positions-only specifications will remain unchanged throughout the optimization process, but the value for target will be different during each iteration. Since dictionaries are mutable, the first constraint will be updated repeatedly during the minimization process. However, because tuples are immutable, the references held by the constraints tuple will always point to the same objects. This nuance makes the implementation Pythonic. We again constrain the weights such that all weights fall within the interval \([0,1]\):

# This creates 12 tuples of (0, 1), all of which exist within a container tuple
# We essentially create a sequence of (min, max) pairs
bounds = tuple(
  (0, 1) for w in weights
)

Here’s the implementation in Python. We will use the scipy.optimize.minimize function again and the Sequential Least Squares Programming (SLSQP) method for the minimization:

# Initialize an array of target returns
target = np.linspace(
  start = 0.15, 
  stop = 0.30,
  num = 100
)
# instantiate empty container for the objective values to be minimized
obj_sd = []
# For loop to minimize objective function
for target in target:
  min_result_object = sco.minimize(
    # Objective function
    fun = portfolio_sd, 
    # Initial guess, which is the equal weight array
    x0 = equal_weights, 
    method = 'SLSQP',
    bounds = bounds, 
    constraints = constraints
    )
  # Extract the objective value and append it to the output container
  obj_sd.append(min_result_object['fun'])
# End of for loop
# Convert list to array
obj_sd = np.array(obj_sd)
# Rebind target to a new array object
target = np.linspace(
  start = 0.15, 
  stop = 0.30,
  num = 100
)

Before we plot the efficient frontier, we may wish to highlight the two portfolios— the maximal Sharpe ratio and the minimum variance portfolios:

# Annotations for maximal Sharpe ratio and minimum variance portfolio
annotation_data <- tibble::tibble(
  x = c(py$max_sharpe_port_sd, py$min_sd_port_sd),
  y = c(py$max_sharpe_port_return, py$min_sd_port_return),
  type = c("Maximal Sharpe Ratio Portfolio", "Minimum Variance Portfolio")
)

Since the optimal expected portfolio returns and standard deviations are both array objects, we can access them via the reticulate package and plot them in R:

plotly::plot_ly(
  x = py$port_sd, y = py$port_returns, color = (py$port_returns / py$port_sd),
  mode = "markers", type = "scattergl", showlegend = FALSE,
  marker = list(size = 5, opacity = 0.7)
) %>%
  # Efficient frontier
  plotly::add_trace(
    data = tibble::tibble(
      Risk = py$obj_sd,
      Return = py$target,
      SharpeRatio = py$target / py$obj_sd
    ),
    x = ~Risk,
    y = ~Return,
    color = ~SharpeRatio,
    marker = list(size = 7)
  ) %>%
  # Maximal Sharpe ratio portfolio
  plotly::add_trace(
    data = tibble::tibble(
      Risk = py$max_sharpe_port_sd,
      Return = py$max_sharpe_port_return,
      SharpeRatio = py$max_sharpe_port_return / py$max_sharpe_port_sd
    ),
    x = ~Risk,
    y = ~Return,
    color = ~SharpeRatio,
    marker = list(size = 7)
  ) %>%
  # Minimum variance portfolio
  plotly::add_trace(
    data = tibble::tibble(
      Risk = py$min_sd_port_sd,
      Return = py$min_sd_port_return,
      SharpeRatio = py$min_sd_port_return / py$min_sd_port_sd
    ),
    x = ~Risk,
    y = ~Return,
    color = ~SharpeRatio,
    marker = list(size = 7)
  ) %>%
  plotly::layout(
    title = "Mean-Standard Deviation Diagram",
    yaxis = list(title = "Expected Portfolio Return (Annualized)", tickformat = ".2%"),
    xaxis = list(title = "Portoflio Standard Deviation (Annualized)", tickformat = ".2%")
  ) %>%
  plotly::add_annotations(
    x = annotation_data[["x"]],
    y = annotation_data[["y"]],
    text = annotation_data[["type"]],
    xref = "x",
    yref = "y",
    showarrow = TRUE,
    arrowhead = 5,
    arrowsize = .5,
    ax = 20,
    ay = -40
  ) %>% 
plotly::colorbar(title = "Sharpe Ratio")

         Evidently, the efficient frontier is comprised of all optimal portfolios with a higher return than the global minimum variance portfolio. These portfolios dominate all other sub-optimal portfolios in terms of expected returns given a certain risk level.

Capital Market Line

         Now that we have determined a set of efficient portfolios of risky assets, we may be interested in adding a risk-free asset to the mix. By adjusting the wealth allocation for the risk-free asset, it is possible for us to achieve any risk-return profile that lies on the straight line (in the \(\mu-\sigma\) diagram) between the risk-free asset and an efficient portfolio. This is called the capital market line, which is drawn from the point of the risk-free asset to the feasible region for risky assets. The capital market line has the equation:

\[\begin{align*} R_{p}=r_{f}+\frac{R_{T}-r_{f}}{\sigma_{T}} \sigma_{p} \end{align*}\]

where

• \(R_{p}=\text{portfolio return}\)

• \(r_{f}=\text{risk free rate}\)

• \(R_{T}=\text{market return}\)

• \(\sigma_{T}=\text{standard deviation of market portfolio}\)

• \(\sigma_{p}=\text{standard deviation of portfolio}\)

         The problem then becomes: which efficient portfolio on the efficient frontier should we use? The answer is the tangency portfolio, at which point on the efficient frontier the slope of the tangent line is equal to that of the capital market line. In order to find the slope of the tangent line, a functional approximation of the efficient frontier and its first derivative must be obtained. Yves Hilpisch’s book recommends the use of cubic splines interpolation to obtain a continuously differentiable functional approximation of the efficient frontier, \(f(x)\). To carry out cubic splines interpolation, we need to create two array objects representing the expected returns and standard deviations of each of the efficient portfolios on the efficient frontier. In other words, we need to exclude from obj_sd and target (our two arrays of optimization results) the portfolios whose standard deviations are greater than that of the minimum variance portfolio but whose expected returns are smaller; they are not optimal and are therefore not on the efficient frontier:

# Efficient portfolio returns
# The function np.argmin() returns the indices of the minimum value(s)
# Take a slice starting from the minimum variance portfolio and exclude previous elements
efficient_sd = obj_sd[np.argmin(a = obj_sd):]
# Take a slice starting from the same index and exclude previous elements
efficient_returns = target[np.argmin(a = obj_sd):]

         We will leverage the the facilities from the scipy.interpolate sub-package, in particular, the splrep() function. Given the set of data points \((x_{i}, y_{i})\), the scipy.interpolate.splrep() function determine a smooth spline approximation of degree \(k\). The function will return a tuple \((t,c,k)\) with three elements— the vector of knots \(t\), the B-spline coefficients \(c\), and the degree of the spline \(k\).

# Cubic spline interpolation
tck = sci.splrep(
  x = np.sort(efficient_sd),
  y = efficient_returns,
  k = 3
)

         Next, we use the scipy.interpolate.splev() function to define the functional approximation and the first derivative of the efficient frontier— \(f(x)\) and \(f^{\prime}(x)\). Given the knots and the B-spline coefficients, the function scipy.interpolate.splev() evaluates the value of the smoothing polynomial and its derivative:

# Define functional approximation of efficient frontier
def f(x):
  return sci.splev(
  x = x,
  tck = tck,
  der = 0
  )
# Define first derivative of the efficient frontier
def df(x):
  return sci.splev(
  x = x,
  tck = tck,
  der = 1
  )

         Let us re-examine the equation for the capital market line. Notice that the line is a linear equation of the form:

\[\begin{align*} f(x)&=a + bx \end{align*}\]

where the parameters are as follows:

• \(a=r_{f}\)

• \(b=\frac{R_{T}-r_{f}}{\sigma_{T}}=f^{\prime}(x) \hspace{5mm} \text{Slope of the capital market line = slope of the tangent line of the efficient frontier}\)

• \(x=\sigma_{p}\)

There are three unknown parameters and therefore we must have three equations to solve for the values of these parameters:

\[\begin{align*} a=r_{f} \hspace{5mm} \Longrightarrow \hspace{5mm} 0&=r_{f}-a \\ b=f^{\prime}(x) \hspace{5mm} \Longrightarrow \hspace{5mm} 0&=b-f^{\prime}(x) \\ f(x)=a + bx \hspace{5mm} \Longrightarrow \hspace{5mm} 0&=a+b\cdot x-f(x) \end{align*}\]

         To solve the system, we will use the scipy.optimize.fsolve() function, which finds the roots of any input function. First, we need to define a function that returns the values of all three equations given a set of parameters \(p = (a, b, x)\). The set of parameters \(p\) can be represented using a python list where \(p[0]=a\), \(p[1]=b\), and \(p[2]=x\):

# System of equations
def sys_eq(p, r_f = 0.01):
  # Equations
  eq1 = r_f - p[0]
  eq2 = p[1] - df(p[2])
  eq3 = r_f + df(p[2]) * p[2] -  f(p[2])
  # Output values
  return eq1, eq2, eq3

The function scipy.optimize.fsolve() returns an ndarray object, which is the solution set to the linear system above:

# Solve the linear system
sol_set = sco.fsolve(
  # Equations to solve for
  func = sys_eq,
  # Initial guess for p
  # This can be determined by trial and error or from the plot above (may take more than one try)
  x0 = [0.05, 1, 0.15]
)

As a sanity check, we an plug our solution set into our user defined function sys_eq to see if the results are indeed zeros:

# Sanity check
np.round(
  sys_eq(
    p = sol_set
  ),
  4
)

array([ 0.,  0., -0.])

         We have confirmed that our solution set for the linear system is as follows:

• \(a=\) 0.01

• \(b=\) 1.0795

• \(x=\) 0.1986

# Constraints
# The target return is now f(x) where x is taken from the solution set above
constraints = (
  {'type': 'eq', 'fun': lambda x: portfolio_returns(x) - f(x = sol_set[2])}, 
  {'type': 'eq', 'fun': lambda x: np.sum(x) - 1}
)
# Optimize
min_result_object_tangency = sco.minimize(
  # Objective function
  fun = portfolio_sd, 
  # Initial guess, which is the equal weight array
  x0 = equal_weights, 
  method = 'SLSQP',
  bounds = bounds, 
  constraints = constraints
)

# Expected return
tangency_port_return = portfolio_returns(min_result_object_tangency["x"])
round(tangency_port_return, 4)

0.2244

# Standard deviation
tangency_port_sd = portfolio_sd(min_result_object_tangency["x"])
round(tangency_port_sd, 4)

0.1986

# Sharpe ratio
tangency_port_sharpe = tangency_port_return / tangency_port_sd
round(tangency_port_sharpe, 4)

1.1299

# Capital market line
def cml(x):
  return sol_set[0] + sol_set[1] * x

# Standard deviations
cml_sd = np.linspace(
  start = 0, 
  stop = 0.25,
  num = 100
)
# Expected returns
cml_exp_returns = cml(x = cml_sd)

# Efficient frontier and Capital market line
plotly::plot_ly(
  x = py$port_sd, y = py$port_returns, color = (py$port_returns / py$port_sd),
  mode = "markers", type = "scattergl", showlegend = FALSE,
  marker = list(size = 5, opacity = 0.7)
) %>%
  # Random portfolios
  plotly::add_trace(
    data = tibble::tibble(
      Risk = py$obj_sd,
      Return = py$target,
      SharpeRatio = py$target / py$obj_sd
    ),
    x = ~Risk,
    y = ~Return,
    color = ~SharpeRatio,
    marker = list(size = 7)
  ) %>%
  # Capital market line
  plotly::add_lines(
    data = tibble::tibble(
      Risk = py$cml_sd,
      Return = py$cml_exp_returns
    ),
    x = ~Risk,
    y = ~Return,
    line = list(color = "#9467bd"),
    # Do no inherit attributes from plot_ly()
    inherit = FALSE
  ) %>%
  plotly::layout(
    title = "Mean-Standard Deviation Diagram",
    yaxis = list(title = "Expected Portfolio Return (Annualized)", tickformat = ".2%"),
    xaxis = list(title = "Portoflio Standard Deviation (Annualized)", tickformat = ".2%"),
    annotations = list(
      x = py$tangency_port_sd,
      y = py$tangency_port_return,
      text = "Tangency Portfolio \n Expected Return: 26.64% \n Volatility: 17.83%",
      xref = "x",
      yref = "y",
      showarrow = TRUE,
      arrowhead = 7,
      ax = 20,
      ay = -40
    )
  ) %>%
  plotly::colorbar(title = "Sharpe Ratio")