Python application 3: classical portfolio theory

Consider a financial economy with three risky securities and one risk-free asset. The (percentage terms) expected returns of assets 1, 2, and 3 are $(7, 10, 7)$, respectively. The variance-covariance matrix of these returns (also in percentage terms) is: $$ \Omega = \left[ \begin{array}{c c c} 49 & 10 & -10 \\ 10 & 100 & 10 \\ -10 & 10 & 49 \end{array} \right].$$

An investor cannot hold negative positions in any of those risky assets. The risk-free rate is $r_F=5\%.$ The investor can borrow and lend however much she wants at that rate. Denote by $\alpha_i$ the non-negative share of her wealth she puts in security $i \in \{1,2,3\}$ while $\alpha_F$ is the (possibly negative) share of her wealth she invests in the risk-free asset.

The investor has mean-variance preferences. Specifically, she seeks to maximize $$E(r_P) - \frac{\sigma\left(r_P\right)^2}{6} ,$$ where $$r_P= \alpha_F r_F + \alpha_1 r_1 + \alpha_2 r_2 + \alpha_3 r_3 $$ is her portfolio return, $E(r_P)$ is the expected value of that return, while $\sigma\left(r_P\right)$ is the standard deviation of that return. For future reference and in matrix notation: $$\sigma\left(r_P\right) = \sqrt{ (\alpha_1 , \alpha_2 , \alpha_3) \Omega (\alpha_1 , \alpha_2 , \alpha_3) '}.$$

First we import the modules we will probably need.

In [50]:
import pandas as pd # pandas is excellent for creating and manipulating dataframes, R-style
import numpy as np # great for simulations, we may not use for running regressions here
import scipy.linalg as la # for matrix algebra
import matplotlib.pyplot as plt #graphing module with matlab-like properties
%matplotlib inline 
import requests # to make requests to webpages (like mine)
import statsmodels.api as sm # module full of standard statistical and econometric models, including OLS and time-series stuff
from IPython.display import Latex # to be able to display latex in python command cells

Now we start building the feasible set. This involves building lots of random portfolios and then calculating the expected return and standard deviation of the returns on those portfolios.

In [51]:
from numpy import random # for random numbers below
import math

r=np.array([7,10,7])
Omega=np.array([[49, 10, -10], [10, 100, 10], [-10 ,10 ,49]])

n=100000 # number of random portfolios

ret=np.zeros((n,1)) # this will store expected returns
std=np.zeros((n,1)) # this will store standard deviations
sharpe=np.zeros((n,1)) # this will sharpe ratios
sharpemarket=0 # market portfolio's sharpe ratio
retmarket=0
stdmarket=0

market=[1/3, 1/3, 1/3] # this will store the (approximate) market portfolio

# and off we go

for i in range(n):
    a=np.array([random.rand(),random.rand(),random.rand()]) # can compactify this syntax but more clear this way
    a=1/sum(a)*a # normalizes things so they sum to one
    ret[i,:]=a @ r.T
    std[i,:]=math.sqrt(a @ Omega @ a.T)
    sharpe[i,:]=(ret[i,:]-5)/std[i,:]
    
    if sharpe[i,:]>sharpemarket:
        sharpemarket=sharpe[i,:]
        market=a
        retmarket=ret[i,:]
        stdmarket=std[i,:]
        
fig, ax = plt.subplots()

plt.plot(std, ret, 'o',markersize=1)

plt.plot(stdmarket, retmarket, 'ro')

plt.plot(0, 5, 'go')

plt.title('Feasible set (market portfolio in red, risk free portfolio in green)')

plt.plot(std, sharpemarket*std + 5)

xnow=[0,stdmarket]

plt.plot(xnow, sharpemarket*xnow + 5)

ax.set_ylabel('Expected return')
ax.set_xlabel('Standard deviation')

print('The approximate market portfolio is: ', market)
print('The market return is: ', retmarket)
print('The std of the market return is: ', stdmarket)
print('The sharpe ratio of the market return is: ', sharpemarket)
The approximate market portfolio is:  [0.32943626 0.34046698 0.33009676]
The market return is:  [8.02140093]
The std of the market return is:  [4.95630484]
The sharpe ratio of the market return is:  [0.60960757]

We can also find the market portfolio by solving $$ \max_{(\alpha_1, \alpha_2,\alpha_3) \geq (0,0,0)} \frac{E(r_P)- r_F}{\sigma(r_P)}$$ subject to $$\alpha_1+ \alpha_2 + \alpha_3 =1.$$

In [52]:
import numpy as np
from scipy.optimize import minimize
from scipy.optimize import Bounds
from scipy.optimize import LinearConstraint

bounds=[(0,1),(0,1),(0,1)]
linear_constraint = LinearConstraint([[1, 1, 1]], 1,1) 

def objective(x):
    return -(x @ r.T-5)/math.sqrt(x @ Omega @ x.T)

x0=np.array([1/2,1/4,1/4])

res = minimize(objective, x0, method='trust-constr', 
               constraints=[linear_constraint],
               options={'verbose': 1}, bounds=bounds)

print('The market portfolio is: ',res.x)



print('The sharpe ratio of the market portfolio is', -res.fun)

   
`gtol` termination condition is satisfied.
Number of iterations: 15, function evaluations: 28, CG iterations: 6, optimality: 6.31e-09, constraint violation: 0.00e+00, execution time: 0.065 s.
The market portfolio is:  [0.32967028 0.34065937 0.32967035]
The sharpe ratio of the market portfolio is 0.6096077604670229

Now we can plot the CAPM line somewhat more exactly.

In [53]:
fig, ax = plt.subplots()

marketexact=np.array([res.x[0],res.x[1],res.x[2]])
retmarketexact=marketexact @ r.T
stdmarketexact=math.sqrt(marketexact @ Omega @ marketexact.T)

sharpemarketexact=(retmarketexact-5)/stdmarketexact
                      
plt.plot(std, ret, 'o',markersize=1)


plt.plot(stdmarketexact, retmarketexact, 'ro')

plt.plot(0, 5, 'go')

plt.title('Feasible set (market portfolio in red, risk free portfolio in green)')

plt.plot(std, sharpemarketexact*std + 5)

xtry=np.array([0,stdmarketexact])

plt.plot(xtry, sharpemarketexact*xtry + 5)

ax.set_ylabel('Expected return')
ax.set_xlabel('Standard deviation')
Out[53]:
Text(0.5, 0, 'Standard deviation')

Now we find the optimal portfolio. We could once again use simulations and get as close to the answer as we want, but we will use optimization instead. We need to solve:

solve: $$ \max_{\alpha_F \leq 1} E(r_P) - \frac{\sigma\left(r_P\right)^2}{6}$$ subject to: \begin{eqnarray*} E(r_P) & = & \alpha_F r_F +(1- \alpha_F) E(r_M) \\ \sigma(r_P) & = & (1-\alpha_F) \sigma(r_M) \end{eqnarray*} Traditional variational arguments show the solution to this to be:

$$\alpha_F=1-3\frac{E(r_M)-r_F}{\sigma(r_M)^2}\approx 0.63107.$$

Below we check with Python's optimization routine that it all works out.

In [54]:
bounds=[(-1,1)]


def objective2(x):
    Eret=x*5+(1-x)*retmarketexact
    Stdret=(1-x)*stdmarketexact
    return -Eret+Stdret**2/6

x0=0.75

res2 = minimize(objective2, x0, method='trust-constr', 
               options={'verbose': 1}, bounds=bounds)

print('The optimal risk-free share is: ',res2.x)
print('At the optimal, the expected return is: ',res2.x*5+(1-res2.x)*retmarketexact)
print('At the optimal, the standard deviation of the return is: ',(1-res2.x)*stdmarketexact)
`gtol` termination condition is satisfied.
Number of iterations: 14, function evaluations: 18, CG iterations: 8, optimality: 6.35e-09, constraint violation: 0.00e+00, execution time: 0.082 s.
The optimal risk-free share is:  [0.63107289]
At the optimal, the expected return is:  [6.11488966]
At the optimal, the standard deviation of the return is:  [1.82886396]

Let's now locate this portfolio on our graph and plot the indifference curve that runs through it to get a nice tangency.

In [55]:
fig, ax = plt.subplots()

marketexact=np.array([res.x[0],res.x[1],res.x[2]])
retmarketexact=marketexact @ r.T
stdmarketexact=math.sqrt(marketexact @ Omega @ marketexact.T)

sharpemarketexact=(retmarketexact-5)/stdmarketexact
                      
plt.plot(std, ret, 'o',markersize=1)


plt.plot(stdmarketexact, retmarketexact, 'ro')

plt.plot((1-res2.x)*stdmarketexact, res2.x*5+(1-res2.x)*retmarketexact, 'yo',markersize=10)

plt.plot(0, 5, 'go')

plt.title('Feasible set (optimal portfolio is big yellow dot)')

plt.plot(std, sharpemarketexact*std + 5)

xtry=np.array([0,stdmarketexact])

plt.plot(xtry, sharpemarketexact*xtry + 5)

# now for good measure we add the indifference curve that runs thru the optimal portfolio

Ustar=res2.x*5+(1-res2.x)*retmarketexact-((1-res2.x)*stdmarketexact)**2/6 # the value of the objective at the optimal portfolio
stdopt=np.linspace(0, stdmarketexact, num=500)
plt.plot(stdopt, Ustar + stdopt**2/6)


ax.set_ylabel('Expected return')
ax.set_xlabel('Standard deviation')
Out[55]:
Text(0.5, 0, 'Standard deviation')