Homework 7
Projection and portfolio Optimization
1 Projection onto a halfspace
A nonzero n-vector \(a\) and a scalar \(\beta\) define the halfspace
\[H^{-}_{a,β}=\{\,x∈\mathbb R^n\mid ⟨a,x⟩≤β\,\}.\]
1.1 Derive the expression for the projection
Derive an explicit expression for the projection of an arbitrary n-vector b onto \(H^-_{a,β}\):
1.2 Implement the projection
Implement your solution above by completing the following function.
1.3 Implement a unit test
Verify numerically with a simple example that your projection code is correct. No need to formally verify correctness.
2 Alternating projection method
The function projsplx, defined here, computes the projection of any n-vector onto the probability simplex \(\{x\in\mathbb R^n\mid \sum_j^nx_j=1,\ x≥0\}\).
Implementation of projsplx
"""
Projection onto the unit simplex.
projsplx(b) -> x
Variant of `projsplx`.
"""
function proj_simplex(b::Vector; τ=1.0)
x = copy(b)
projsplx!(x, τ)
return x
end;
"""
Projection onto the unit simplex {x | sum(x) = τ, x ≥ 0}.
projsplx!(b, τ)
In-place variant of `projsplx`.
"""
function projsplx!(b::Vector{T}, τ::T) where T
n = length(b)
bget = false
idx = sortperm(b, rev=true)
tsum = zero(T)
@inbounds for i = 1:n-1
tsum += b[idx[i]]
tmax = (tsum - τ)/i
if tmax ≥ b[idx[i+1]]
bget = true
break
end
end
if !bget
tmax = (tsum + b[idx[n]] - τ) / n
end
@inbounds for i = 1:n
b[i] = max(b[i] - tmax, 0)
end
end;The cell below shows the result of the projection of a vector b onto the simplex.
Now suppose that we wish to compute the projection onto the intersection of the probability simplex and an arbitrary halfspace. Can you do this using the projection function proj_halfspace, which you implemented above, and proj_simplex?
Yes, you can, using the Projection Onto Convex Sets algorithm, which projects onto the intersection of sets \(C_1\bigcap C_2\) each convex set by itereratively projecting onto \(C_1\) and \(C_2\). The algorithm is simple, and involves the following iteration:
\[x_{k+1}=P_1(P_2(x_k)),\]
where \(P_1\) and \(P_2\) are projection operators, respectively, onto the sets \(C_1\) and \(C_2\). Implement this algorithm for generic projection operators \(P_1\) and \(P_2\). Exit the loop when \[\|x_{k+1}-x_k\|≤ϵ\|x_0\|\] for some positive tolerance \(\epsilon\).
Note that this method can in some situations take many iterations. But here we’re only dealing with a small problem, so it should be fine.
2.1 Implement the intersection projection
Complete the following function.
"""
Project a vector `x` onto the intersection of two convex sets. Takes as input two functions `p1` and `p2` that each take a vector and return the projection of `b0` onto the corresponding sets.
"""
function proj_intersection(x::Vector, p1::Function, p2::Function; maxits=100, ϵ=1e-7)
err = []
x0 = copy(x)
ϵ *= norm(x0, Inf)
for i ∈ 1:maxits
# your code here
push!(err, norm(x-x0, Inf))
err[end] ≤ ϵ && break
x0 = copy(x)
end
return x
end;2.2 Implement a unit test
Test your implementation of proj_intersection on a small problem that aims to find the projection of a vector onto the intersection of the unit probability simplex and a halfspace. Verify that the answer is correct.
Make sure that your test problem is feasible – i.e., that the intersection isn’t empty.
3 Financial portfolio optimization
Modern portfolio optimization theory is based on the Markovitz model for determining a portfolio of stocks with a desired expected rate of return that has the smallest amount of risk, as measured by the variance of the portfolio’s return. The main idea is that by diversifying (i.e., investing in a mixture of different stocks), one can guard against large amounts of variance in the rates of return of the individual stocks in the portfolio.
Rate of return and risk
Suppose that \(p_1,\ldots,p_m\) are the historical prices of a stock over some period of time. We define the rPPate of return at time \(t\), relative to the intial price \(p_1\) by \[ r_t := (p_t-p_1)/p_1, \quad t=1,\ldots,m. \tag{1}\]
The expected rate of return is the mean \(\mu\) of the rates of return. The risk of the portfolio is measured as the standard deviation \(\sigma\) of the rates of return: \[ \mu := \frac{1}{m}\sum_{t=1}^m r_t, \quad\text{and}\quad \sigma:=\sqrt{\frac{1}{m}\sum_{t=1}^m (r_t-\mu)^2}. \]
Given a collection of \(n\) stocks, let \(r_t^i\) be the rate of return of stock \(i\) at time \(t\). Let \(r\) be the \(n\)-vector of the expected rates of return of the \(n\) stocks. In addition, let \(\Sigma\) be the \(n\times n\) covariance matrix of the rates of return of the \(n\) stocks. Thus, \(r_i\) is the mean of the rates of return of stock \(i\), the \(i\)th diagonal element \(\Sigma_{ii}\) is its variance, and \(\Sigma_{ij}\) is the covariance of the rates of return of stocks \(i\) and \(j\), i.e., \[ r_i:=\frac{1}{m}\sum_{t=1}^m r_t^i \quad\text{and}\quad \Sigma_{ij}:=\frac{1}{m}\sum_{t=1}^m (r_t^i-r_i)(r_t^j-r_j). \]
Portfolio expected return and risk
We let \(x_i\) be the fraction of our investment money we put into stock \(i\), for \(i=1,\ldots,n\). For the sake of this study, we assume there is no short selling (i.e., holding a stock in negative quantity.) Thus, \(x=(x_1,\ldots,x_n)\) is a vector of length \(n\) that has nonnegative entries that sum to one (i.e, \(x\ge0\) and \(\sum_{i=1}^n x_i=1\)). The vector \(x\) represents our portfolio of investments. The expected rate of return and standard deviation of a portfolio \(x\) are then given by \[ \mu:= r^T x \quad\text{and}\quad \sigma:=\sqrt{x^T\Sigma x}. \]
The code below downloads one year (2022-2023) of historical prices of a collection of stocks and implements meancov, which computes the expected rate of return and covariance matrix of the rates of return.
Download stock data and implement meancov
stocks =
["AAPL", "META", "GOOG", "MSFT", "NVDA", # technology
"AMZN", "COST", "EBAY", "TGT", "WMT", # services
"BMO", "BNS", "HBCP", "RY", "TD", # finance
"BP", "CVX", "IMO", "TTE", "XOM" # energy
]
tmpdir = mktempdir()
for s in stocks
url = "https://query1.finance.yahoo.com/v7/finance/download/$(s)?period1=1647815267&period2=1679351267&interval=1d&events=history&includeAdjustedClose=true"
Base.download(url, joinpath(tmpdir,"$(s).csv"))
end
df = map(stocks) do stock
f = joinpath(tmpdir,"$(stock).csv")
dfs = DataFrame(CSV.File(f, select=[1,6])) # keep only date and closing price
rename!(dfs, "Adj Close" => stock)
select!(dfs, :Date, "$stock" .=> p->(p.-p[1])/p[1])
rename!(dfs, "$(stock)_function" => "$stock")
end |> z->innerjoin(z..., on=:Date)
function meancov(df)
r = mean.(eachcol(df[:,2:end]))
Σ = cov(Matrix(df[:,2:end]))
return r, Σ
end;Here are the expected rates of return for each stock.
3.1 Implement a projected gradient minimizer
The Markovitz model that chooses a portfolio allocation among the given \(n\) stocks that minimizes the risk given a desired expected rate of return \(\mu\). The optimization problem is \[ \min{x}\ \set{\tfrac 12x^T \Sigma x \mid r^T x \ge \mu, \ x\ge0, \ \sum_{i=1}^n x_i = 1}, \] which minimizes a positive semidefinite quadratic objective subject to the intersection of a halfspace (\(r^T x \ge \mu\)) and the probability simplex.
The function below is meant to minimize a quadratic function \(\frac12 x^T Qx\) (with no or constant terms) over a convex set \(C\). Complete the implementation.
3.2 Compute a minimum-risk portfolio
Use the function pgrad to obtain the minimum-risk portfolio from a set of assets with returns r and covariances Σ. Use proj_intersection as the input projection function for pgrad.
"""
Given a set of `n` assets with expected returns `r` and covariance matrix `Σ`, compute the minimum-risk portfolio with expected return `μ`.
"""
function efficient_portfolio(r, Σ, μ)
n = length(r)
p1 = # your code here
p2 = # your code here
p = b -> proj_intersection(b, p1, p2, maxits=10000, ϵ=1e-8)[1]
x = pgrad(Σ, p, maxits=1000)
end;Here we use the function efficient_portfolio to compute an optimal portfolio with expected return of 10%. Show in a pie chart the allocation of the portfolio among the stocks.
3.3 Compute the efficient frontier (optional)
The efficient frontier is the set of optimal portfolios that offer the minimum risk for a given expected rate of return. Use the function efficient_portfolio to compute the efficient frontier for expected rates of return between 0 and 10%. Plot the objective value (risk) for 10 values of expected return between 0 and 10%.
# Define target returns
target_returns = range(0, 0.1, length=10)
# Initialize arrays to store results
portfolio_risks = Float64[]
portfolio_returns = Float64[]
# Loop over target returns
for μ in target_returns
# your code here
endPlot the efficient frontier.