Optimized Portfolio Optimization using EIA Data in Python with AMPL#

# Install dependencies
%pip install amplpy yfinance matplotlib seaborn pandas numpy requests -q

# Google Colab & Kaggle integration
from amplpy import AMPL, ampl_notebook

ampl = ampl_notebook(
    modules=["gurobi"],  # modules to install
    license_uuid="default",  # license to use
)  # instantiate AMPL object and register magics

# Import libraries
import yfinance as yf
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import requests

# Define commodities and fetch data
COMMODITIES = [
    "CL=F",
    "GC=F",
    "NG=F",
    "SI=F",
    "SPY",
]  # Crude Oil, Gold, Natural Gas, Silver, SPY,

START_DATE = "2020-01-01"
END_DATE = "2025-01-01"

def fetch_crude_inventories(api_key):
    """Fetch crude oil inventories data from EIA API."""
    base_url = "https://api.eia.gov/v2/petroleum/stoc/wstk/data/"
    params = {
        "frequency": "weekly",
        "data[0]": "value",
        "facets[series][]": "WTTSTUS1",
        "sort[0][column]": "period",
        "sort[0][direction]": "desc",
        "offset": 0,
        "length": 5000,
        "api_key": api_key,
    }
    response = requests.get(base_url, params=params)
    if response.status_code == 200:
        data = response.json()
        df = pd.DataFrame(data["response"]["data"])
        df["period"] = pd.to_datetime(df["period"])
        df.set_index("period", inplace=True)
        return df[["value"]].rename(columns={"value": "Inventories"})
    else:
        raise Exception(f"Failed to fetch data: {response.status_code}")

def fetch_commodities_data(tickers, start_date, end_date):
    """Fetch historical data for given tickers."""
    data = yf.download(tickers, start=start_date, end=end_date)["Close"]
    return data

def calculate_statistics(data):
    """Calculate daily returns and covariance matrix."""
    returns = data.pct_change(fill_method=None).dropna()
    cov_matrix = returns.cov().values
    mean_returns = returns.mean().values
    return returns, cov_matrix, mean_returns

def plot_pie_chart(weights, title):
    """Plot a pie chart of portfolio weights with improved visuals."""

    # Define a dictionary to map tickers to their respective names
    ticker_to_name = {
        "CL=F": "Crude Oil",
        "GC=F": "Gold",
        "NG=F": "Natural Gas",
        "SI=F": "Silver",
        "SPY": "S&P 500",
    }

    # Map the tickers in the portfolio to their actual names
    labels = [ticker_to_name[ticker] for ticker in weights.index]

    # Filter out the tickers with weight less than or equal to 1e-4
    filtered_weights = weights[weights > 1e-4]
    filtered_labels = [
        label
        for label, ticker in zip(labels, weights.index)
        if ticker in filtered_weights.index
    ]

    # Set color palette
    palette = sns.color_palette("Set3", len(filtered_weights))

    # Plot the pie chart
    plt.figure(figsize=(8, 8))
    wedges, texts, autotexts = plt.pie(
        filtered_weights,
        labels=filtered_labels,
        autopct="%1.1f%%",
        startangle=90,
        colors=palette,
        explode=[0.1]
        * len(filtered_weights),  # Slightly "explode" all slices for emphasis
        wedgeprops={"edgecolor": "black", "linewidth": 1, "linestyle": "solid"},
    )

    # Improve text labels
    for text in texts:
        text.set_fontsize(14)
        text.set_fontweight("bold")
    for autotext in autotexts:
        autotext.set_fontsize(12)
        autotext.set_fontweight("bold")

    # Add title
    plt.title(title, fontsize=18, fontweight="bold")

    # Equal aspect ratio ensures that pie is drawn as a circle.
    plt.axis("equal")

    # Show the plot
    plt.show()

def mean_variance_optimization(tickers, cov_matrix, mean_returns, risk_tolerance):
    """Perform Markowitz mean-variance optimization."""
    ampl = AMPL()
    # Define the AMPL model
    ampl.eval(
        r"""
        set A ordered;  # Set of assets
        param S{A, A};  # Covariance matrix
        param mu{A};    # Expected returns
        param gamma;    # Risk tolerance
        var w{A} >= 0 <= 1;  # Portfolio weights

        maximize expected_return:
            sum {i in A} mu[i] * w[i];  # Objective: maximize expected return

        s.t. portfolio_weights:
            sum {i in A} w[i] = 1;  # Constraint: weights sum to 1

        s.t. risk_constraint:
            sum {i in A, j in A} w[i] * S[i, j] * w[j] <= gamma;  # Risk constraint
    """
    )

    # Set data in AMPL
    ampl.set["A"] = tickers
    ampl.param["S"] = pd.DataFrame(cov_matrix, index=tickers, columns=tickers)
    ampl.param["mu"] = pd.Series(mean_returns, index=tickers)
    ampl.param["gamma"] = risk_tolerance

    # Solve the optimization problem
    ampl.solve(solver="gurobi", mp_options="outlev=1")

    # Extract and return portfolio weights
    return pd.Series(ampl.var["w"].to_dict())

def risk_parity_optimization(tickers, cov_matrix, total_volatility):
    """Perform risk parity optimization."""
    ampl = AMPL()
    # Define the AMPL model
    ampl.eval(
        r"""
        set A ordered;  # Set of assets
        param S{A, A};  # Covariance matrix
        param total_volatility;  # Total portfolio volatility
        var w{A} >= 0 <= 1;  # Portfolio weights

        minimize risk_parity:
            sum {i in A} (w[i] * sum {j in A} S[i, j] * w[j]) / total_volatility;  # Objective: risk parity

        s.t. portfolio_weights:
            sum {i in A} w[i] = 1;  # Constraint: weights sum to 1
        """
    )

    # Set data in AMPL
    ampl.set["A"] = tickers
    ampl.param["S"] = pd.DataFrame(cov_matrix, index=tickers, columns=tickers)
    ampl.param["total_volatility"] = total_volatility

    # Solve the optimization problem
    ampl.solve(solver="gurobi", mp_options="outlev=1")

    # Extract and return portfolio weights
    return pd.Series(ampl.var["w"].to_dict())

def max_sharpe_optimization(tickers, cov_matrix, mean_returns, risk_free_rate=0.02):
    """Perform maximum Sharpe ratio optimization."""
    ampl = AMPL()
    # Define the AMPL model
    ampl.eval(
        r"""
        set A ordered;  # Set of assets
        param S{A, A};  # Covariance matrix
        param mu{A};  # Mean returns
        param risk_free_rate;  # Risk-free rate
        var w{A} >= 0 <= 1;  # Portfolio weights

        maximize sharpe_ratio:
            (sum {i in A} mu[i] * w[i] - risk_free_rate) /
            sqrt(sum {i in A, j in A} w[i] * S[i, j] * w[j]);  # Objective: maximize Sharpe ratio

        s.t. portfolio_weights:
            sum {i in A} w[i] = 1;  # Constraint: weights sum to 1
        """
    )

    # Set data in AMPL
    ampl.set["A"] = tickers
    ampl.param["S"] = pd.DataFrame(cov_matrix, index=tickers, columns=tickers)
    ampl.param["mu"] = pd.Series(mean_returns, index=tickers).to_dict()
    ampl.param["risk_free_rate"] = risk_free_rate

    # Solve the optimization problem
    ampl.solve(solver="gurobi", mp_options="outlev=1")

    # Extract and return portfolio weights
    return pd.Series(ampl.var["w"].to_dict())

crude_inventories = fetch_crude_inventories(os.environ["EIA_API_KEY"])
commodities_data = fetch_commodities_data(COMMODITIES, START_DATE, END_DATE)

returns, cov_matrix, mean_returns = calculate_statistics(commodities_data)

risk_tolerance = 0.01  # Set risk tolerance (gamma)
mean_variance_weights = mean_variance_optimization(
    COMMODITIES, cov_matrix, mean_returns, risk_tolerance
)
plot_pie_chart(mean_variance_weights, "Mean-Variance Optimization Portfolio")

total_volatility = np.sqrt(np.dot(mean_returns, np.dot(cov_matrix, mean_returns)))
risk_parity_weights = risk_parity_optimization(
    COMMODITIES, cov_matrix, total_volatility
)
plot_pie_chart(risk_parity_weights, "Risk Parity Optimization Portfolio")

# Maximum Sharpe Ratio Optimization
max_sharpe_weights = max_sharpe_optimization(COMMODITIES, cov_matrix, mean_returns)
plot_pie_chart(max_sharpe_weights, "Maximum Sharpe Ratio Portfolio")