Homework 7

Projection and portfolio Optimization

Modified

March 18, 2025

Import packages

using LinearAlgebra, DataFrames, Statistics, StatsPlots, YFinance, Dates, CSV

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 projecting an arbitrary n-vector b onto the halfspace \(H^-_{a,β}\).

1.2 Implement the projection

Implement your derived solution by completing the following function.

"""
Project the n-vector `b` onto the halfspace {x ∣ ⟨a,x⟩ ≤ β} defined by the nonzero vector `a` and the scalar `β`.
"""
function proj_halfspace(b, a, β)
  # your code here
end;

1.3 Implement a unit test

Verify numerically with a simple example that your projection code is correct. Formal verification is not required.

function test_proj_halfspace()
    a = ...
    β = ...
    b = ...
    expected = ...
    result = proj_halfspace(b, a, β)
    result == expected
end
test_proj_halfspace()

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 demonstrates the projection of a vector b onto the simplex.

let
    b = [12.0, -23.0, 16.0]
    x = proj_simplex(b)
end

3-element Vector{Float64}:
 0.0
 0.0
 1.0

Now, suppose we want to compute the projection onto the intersection of the probability simplex and an arbitrary halfspace. Can we achieve this using the projection functions proj_halfspace and proj_simplex?

Yes, we can use the Projection Onto Convex Sets algorithm, which projects onto the intersection of convex sets \(C_1\bigcap C_2\) by itereratively projecting onto \(C_1\) and \(C_2\). The algorithm is straightforward and involves the following iteration:

\[x_{k+1}=P_1(P_2(x_k)),\]

where \(P_1\) and \(P_2\) are projection operators onto the sets \(C_1\) and \(C_2\), respectively. 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 require many iterations. However, since we’re dealing with a small problem, this should not be an issue.

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.

Warning

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
]

# Calculate dates for one year of data
end_date = today()
start_date = end_date - Year(1)

# Create cache directory if it doesn't exist
cache_dir = joinpath(@__DIR__, "cache")
mkpath(cache_dir)

"""
    fetch_and_process_stock(symbol::String, start_date::Date, end_date::Date, cache_dir::String) -> DataFrame

Fetch stock data for the given symbol and process it to calculate relative returns from the first day.
If cached data exists and is recent enough, use that instead of downloading.

Returns a DataFrame with columns `Date` and `relative returns`.
"""
function fetch_and_process_stock(symbol::String, start_date::Date, end_date::Date, cache_dir::String)
    cache_file = joinpath(cache_dir, "$(symbol)_$(start_date)_$(end_date).csv")
    
    # Use cache if it exists and is recent (less than 1 day old)
    if isfile(cache_file) && (time() - mtime(cache_file)) < 24*60*60
        df = DataFrame(CSV.File(cache_file))
        return df
    end
    
    # Get daily stock data
    stock_data = get_prices(
        symbol,
        startdt=start_date,
        enddt=end_date,
        interval="1d"
    )
    
    # Convert to DataFrame and process
    df = DataFrame(stock_data)
    
    # Calculate relative returns
    initial_price = first(df.adjclose)
    df[!, symbol] = @. (df.adjclose - initial_price) / initial_price
    
    # Keep only date and relative returns
    select!(df, :timestamp => :Date, symbol)
    
    # Save to cache
    CSV.write(cache_file, df)
    
    return df
end

# Fetch and process all stocks
df = mapreduce(
    s -> fetch_and_process_stock(s, start_date, end_date, cache_dir),
    (x, y) -> innerjoin(x, y, on=:Date),
    stocks
)

"""
    meancov(df) -> Tuple{Vector{Float64}, Matrix{Float64}}

Compute the expected rate of return and covariance matrix from a DataFrame of stock returns.

Returns:
- r: Vector of expected returns for each stock
- Σ: Covariance matrix of returns
"""
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.

Plot expected rates of return and risk

r, Σ = meancov(df)
df_plot = DataFrame(stock=stocks, r=r, risk=diag(Σ))
sort!(df_plot, :r)

# Create primary y-axis for returns
p = @df df_plot bar(:stock, :r, legend=false, ylabel="Return", title="Returns and Risk of Stocks")
xlabel!("Stock")

# Add secondary y-axis for risk
p2 = twinx(p)
@df df_plot scatter!(p2, :stock, :risk, color=:red, markersize=4, 
                   marker=:circle, label="Risk", legend=:topleft)
ylabel!(p2, "Risk (Variance)")

In-class: Stock Analysis - NVDA

Let’s examine NVIDIA (NVDA) stock in more detail — it’s been one of the most talked-about technology stocks in recent years because of AI advancements.

# Extract NVDA data
nvda_idx = findfirst(s -> s == "NVDA", stocks)
nvda_returns = df[:, nvda_idx+1] # +1 because first column is Date

# Calculate rate of return and risk
nvda_return = mean(nvda_returns)
nvda_risk = var(nvda_returns)

println("NVDA Rate of Return: ...")
println("NVDA Risk (Variance): ...")

Clicker Question 1

What is the approximate rate of return for NVDA over the past year?

3.1 Markowitz model

The Markowitz portfolio optimization model finds the allocation of investments across \(n\) stocks that minimizes risk while achieving a target expected return \(\mu\). The model is formulated as: \[ \min_{x}\ \set{\tfrac{1}{2} x^T \Sigma x \mid r^T x \ge \mu,\ x \ge 0,\ \sum_{i=1}^n x_i = 1} \] Here, \(\Sigma\) is the covariance matrix representing risk, \(r\) contains the expected returns of each stock, and \(x\) represents the portfolio weights. The constraints ensure the portfolio achieves the target return \(\mu\), maintains non-negative weights (no short-selling), and fully allocates the investment (weights sum to 1). The objective function minimizes portfolio variance subject to these constraints, balancing risk and return.

In-class: Solving the Markowitz Model with Convex.jl

In this exercise, we’ll implement the Markowitz portfolio optimization model using Convex.jl, which provides a domain-specific language for expressing convex optimization problems, and COSMO.jl as our solver.

First, let’s install and load the necessary packages:

# Install packages if needed
# using Pkg
# Pkg.add(["Convex", "COSMO"])

using Convex, COSMO

Now, let’s formulate the Markowitz portfolio optimization problem using Convex.jl. We’ll create a function that takes in the expected returns r, covariance matrix Σ, and target return μ:

"""
Solve the Markowitz portfolio optimization model using Convex.jl.

Given a set of assets with expected returns `r` and covariance matrix `Σ`,
find the portfolio with minimum risk (variance)
that achieves a target return of `μ`.

Returns the optimal portfolio weights.
"""
function markowitz_portfolio(r, Σ, μ)
    n = length(r)
    
    # Define variables
    x = Variable(n)
    
    # Define objective: minimize risk (variance)
    # Hint: Use quadform() function to represent the quadratic term x'Σx
    objective = # Your code here
    
    # Define constraints
    constraints = [
        # 1. Expected return constraint: portfolio return ≥ μ
        # 2. Budget constraint: sum of weights = 1
        # 3. No short selling: all weights ≥ 0
        # Your code here
    ]
    
    # Create problem
    problem = minimize(objective, constraints)
    
    # Solve the problem with a silent optimizer
    solve!(problem, () -> COSMO.Optimizer(verbose=false), silent=true)
    
    # Return optimal portfolio weights
    return evaluate(x)
end;

Use your solution to the previous question to solve the problem for a target return of 10% using our stock data:

# Solve for target return of 10%
target_return = 0.10
optimal_weights = markowitz_portfolio(r, Σ, target_return)

# Calculate the achieved return and risk
achieved_return = dot(r, optimal_weights)
achieved_risk = sqrt(optimal_weights' * Σ * optimal_weights)

println("Optimal portfolio:")
println("Target return: $(round(target_return * 100, digits=2))%")
println("Achieved return: $(round(achieved_return * 100, digits=2))%")
println("Portfolio risk (std dev): $(round(achieved_risk, digits=4))")

Optimal portfolio:
Target return: 10.0%
Achieved return: 10.0%
Portfolio risk (std dev): 0.0312

Clicker Question 2

Which stock has the highest weight in the optimal portfolio?

NVDA
META
GOOG
TTE

3.2 Implement a projected gradient minimizer

Complete the implementation of the function below, which minimizes a quadratic function \(\frac12 x^T Qx\) (with no or constant terms) over a convex set \(C\).

"""
Projected gradient method for the quadratic optimization problem

    minimize_{x}  1/2 x'Qx  subj to  x ∈ C.

The function `proj(b)` compute the projection of the vector `b` onto the convex set 𝒞.
"""
function pgrad(Q, proj::Function; maxits=100, ϵ=1e-5)
    # your code here
    return x
end;

3.3 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)
    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.

xp1 = efficient_portfolio(r, Σ, .1)
pie(stocks, xp1, title="Portfolio Allocation")

3.4 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
end

Plot the efficient frontier.

plot(portfolio_returns, portfolio_risks, label="Efficient Frontier", 
    xlabel="Return", ylabel="Risk", title="Portfolio Risk vs Return")