import numpy as np
import cvxpy as cp
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 = cp.Variable((n, 1))
# TODO: your implementation here
objective = cp.Minimize(cp.quad_form(x, C))
constraints = [x >= 0]
soc_constraints = [cp.norm2(cp.vstack([2 * R @ x, cp.reshape((x[i] / theta[i]) - (C @ x)[i], (1,1))])) <= (x[i] / theta[i]) + (C @ x)[i] for i in range(n)]
prob = cp.Problem(objective, constraints + soc_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");