import numpy as np
import cvxpy as cvx
import matplotlib.pyplot as plt
plt.style.use('seaborn')
# load data
R = np.load("R.npy")
n = R.shape[0]
theta = np.ones((n, 1)) / n # assume uniform budget
C = R.T@R
r = np.ones((1, n)) # expected return
mu = 10
def find_portfolio(R, C, theta, n, r, mu):
x = cvx.Variable((n, 1))
# TODO: your implementation here
constraints = [
cvx.SOC(x[i] + (C@x)[i],
cvx.vstack([2*(theta[i]**0.5)*(R@x),
cvx.reshape((x[i] - (C@x)[i]), (1,1))])) for i in range(n)
] + [x >= 0, r@x >= mu]
obj = cvx.Minimize(cvx.quad_form(x, C))
prob = cvx.Problem(obj, constraints)
prob.solve()
return x.value
x = find_portfolio(R, C, theta, n, r, mu)
risk = np.diag(C@x@x.T)
total_risk = np.sum(risk)
for i in range(n):
assert(np.isclose(x[i]*C[i]@x, theta[i]*x.T@C@x))
# Plot elements of the optimal portoflio and partial risk as a function of the index
fig, (ax1, ax2) = plt.subplots(2, 1, constrained_layout=True, sharex=True)
ax1.stem(x)
ax1.set_title("Portfolio Positiions")
ax1.set_ylabel("Optimal Portfolio")
# We see that the optimal partial risk is equal to uniform budgeting
ax2.stem(risk/total_risk, label="partial risk")
ax2.plot(theta, 'r*', markersize=5, label="theta")
ax2.legend()
ax2.set_ylim(top=.015)
ax2.set_title("Risk component divided by total risk")
ax2.set_xlabel("Index")
ax2.set_ylabel("Partial risk");