Stock Portfolio Optimization with Python

# Retreive Stock index from "Companies by Index" API
import time
import requests
from google.colab import userdata

# Retrieve the API key securely
api_key = userdata.get('SECTORS_API_KEY')

# Define the API URL
url = "https://api.sectors.app/v1/index/sminfra18/"

# Pass the API key in the header
headers = {"Authorization": api_key}

# Make the API request
response_company_index = requests.get(url, headers=headers)

print(response_company_index.text)

[
  {
    "symbol": "ADHI.JK",
    "company_name": "PT Adhi Karya (Persero) Tbk."
  },
  {
    "symbol": "AKRA.JK",
    "company_name": "PT AKR Corporindo Tbk."
  },
  {
    "symbol": "BBNI.JK",
    "company_name": "PT Bank Negara Indonesia (Persero) Tbk"
  },
  ...
]

# Retrieve date and price from "Daily Transaction Data" API
from datetime import datetime, timedelta

# Function to calculate the date 90 days ago from today
def calculate_start_date(days_ago=90):
    return (datetime.now() - timedelta(days=days_ago)).strftime('%Y-%m-%d')

# Calculate the start date 90 days ago
start_date = calculate_start_date()

# Looping API info
history_sminfra18 = []

for i, x in enumerate(response_company_index.json()):

  # Define the URL for the API endpoint
  url = "https://api.sectors.app/v1/daily/" + response_company_index.json()[i]['symbol'] + "/"

  # Define the query string with the calculated start date
  querystring = {"start": start_date}

  headers = {"Authorization": userdata.get('SECTORS_API_KEY')}

  response_daily_transaction_data = requests.request("GET", url, headers=headers, params=querystring)

  # Append the result into target list
  history_sminfra18.append(response_daily_transaction_data.json())

  time.sleep(1)

history_sminfra18

[[{'symbol': 'ADHI.JK',
 'date': '2024-07-10',
 'close': 250,
 'volume': 80145800,
 'market_cap': 2101902442496},
{'symbol': 'ADHI.JK',
 'date': '2024-07-11',
 'close': 248,
 'volume': 37656100,
 'market_cap': 2085087215616},
{'symbol': 'ADHI.JK',
 'date': '2024-07-12',
 'close': 262,
 'volume': 138657900,
 'market_cap': 2202793672704},
 ...
]]

# Import all useful libraries
import numpy as np
import pandas as pd

# Flatten the list of lists into a single list of dictionaries
flattened_data = [item for sublist in history_sminfra18 for item in sublist]

# Convert to a pandas DataFrame
df_history_sminfra18 = pd.DataFrame(flattened_data)

# Ensure 'date' is in datetime format
df_history_sminfra18['date'] = pd.to_datetime(df_history_sminfra18['date'])

# Enforce in case any column needs conversion (e.g., 'close', 'volume', 'market_cap')
df_history_sminfra18['close'] = pd.to_numeric(df_history_sminfra18['close'], errors='coerce')
df_history_sminfra18['volume'] = pd.to_numeric(df_history_sminfra18['volume'], errors='coerce')
df_history_sminfra18['market_cap'] = pd.to_numeric(df_history_sminfra18['market_cap'], errors='coerce')

# Check for the first few rows to ensure everything is correct
print(df_history_sminfra18)

      symbol       date  close     volume      market_cap
0    ADHI.JK 2024-07-10    250   80145800   2101902442496
1    ADHI.JK 2024-07-11    248   37656100   2085087215616
2    ADHI.JK 2024-07-12    262  138657900   2202793672704
3    ADHI.JK 2024-07-15    256   49577700   2152348123136
4    ADHI.JK 2024-07-16    248   29521600   2085087257962
..       ...        ...    ...        ...             ...
409  UNTR.JK 2024-08-05  24425    8787200  88071391412224
410  UNTR.JK 2024-08-06  24150    7129900  87708215017472
411  UNTR.JK 2024-08-07  24725    4726300  89796500258816
412  UNTR.JK 2024-08-08  24625    3126500  89433323864064
413  UNTR.JK 2024-08-09  24700    2528100  91340020908032

[414 rows x 5 columns]

 We aim to have all the stock price evolution in the past 90 days in one dataframe

#Pivot the DataFrame so that the 'symbol' becomes the column names and 'date' becomes the index
df_pivot_history_sminfra18 = df_history_sminfra18.pivot_table(index='date', columns='symbol', values='close')
df_pivot_history_sminfra18.head()

df_sminfra18_company_report = pd.DataFrame(prepared_data)
df_sminfra18_company_report

# Categorize stocks by industries
df_sminfra18_company_report.industry.value_counts()

# Categorize stocks by sectors
df_sminfra18_company_report.sector.value_counts()

# Calculate the correlation of the 18 stock prices
correlation_matrix = df_pivot_history_sminfra18.corr()
correlation_matrix

# Import library for normalization
from sklearn.preprocessing import MinMaxScaler

# Extract useful information for normalization analysis
columns_of_interest = ['symbol', 'eps_growth', 'avg_yield_dividends', 'revenue_growth', 'revenue_estimate']
df_sminfra18_company_report_norm = df_sminfra18_company_report[columns_of_interest]
df_sminfra18_company_report_norm

# Check missing data
df_sminfra18_company_report_norm.isnull().sum()

# Impute missing data with mean - we do a simple imputation here, can also check the reason and correlation with other variables behind missing data
df_sminfra18_company_report_norm.fillna(df_sminfra18_company_report_norm.eps_growth.mean(), inplace=True)
df_sminfra18_company_report_norm.isnull().sum()
# As result, all variable should show 0 missing 

# Normalize the data for the metrics using MinMaxScaler
scaler = MinMaxScaler()
df_sminfra18_company_report_norm[['eps_growth', 'avg_yield_dividends', 'revenue_growth', 'revenue_estimate']] = scaler.fit_transform(df_sminfra18_company_report_norm[['eps_growth', 'avg_yield_dividends', 'revenue_growth', 'revenue_estimate']])

# Calculate the score by summing the normalized scores for all metrics
df_sminfra18_company_report_norm['score'] = df_sminfra18_company_report_norm.iloc[:, 1:5].sum(axis=1)

# Reorder columns to place 'score' as the third column, after 'symbol' and before 'eps_growth'
columns = list(df_sminfra18_company_report_norm.columns)
columns.insert(1, columns.pop(columns.index('score')))
df_sminfra18_company_report_norm = df_sminfra18_company_report_norm[columns]

# Display the DataFrame
df_sminfra18_company_report_norm

# Sort the companies by total score
df_sminfra18_company_report_norm = df_sminfra18_company_report_norm.sort_values(by='score', ascending=False).reset_index(drop=True)
df_sminfra18_company_report_norm

# Pick the top 5 companies
top_5_companies = df_sminfra18_company_report_norm.head(5)
top_5_companies

# Store the top 5 company symbol into a list
top_5_symbols = top_5_companies['symbol'].tolist()
top_5_symbols

# Query the sector & industry information for the top 5 companies
top_5_info = df_sminfra18_company_report[df_sminfra18_company_report['symbol'].isin(top_5_symbols)]
top_5_info

# Query the price history info for selected 3 tickers and check correlation
selected_3_symbols =  ['BBRI.JK', 'TLKM.JK', 'UNTR.JK']
selected_3_prices = df_pivot_history_sminfra18[selected_3_symbols]
correlation_matrix_selected = selected_3_prices.corr()
correlation_matrix_selected

# Normalize prices to start from 100 for easy comparison
price_norm = price_Indra / price_Indra.iloc[0] * 100

# Calculate percentage gains from the initial price
price_gains = ((price_Indra - price_Indra.iloc[0]) / price_Indra.iloc[0]) * 100

# Plotting normalized prices and percentage gains
plt.figure(figsize=(14, 6))

# Plot normalized prices
plt.subplot(1, 2, 1)
price_norm.plot(ax=plt.gca(), title="Normalized Price Performance")
plt.ylabel('Normalized Price (Start = 100)')
plt.xlabel('Date')
plt.legend(title='Stock Symbol')
plt.grid(True)

# Plot percentage gains
plt.subplot(1, 2, 2)
price_gains.plot(ax=plt.gca(), title="Percentage Gains Over Time")
plt.ylabel('Percentage Gain (%)')
plt.xlabel('Date')
plt.legend(title='Stock Symbol')
plt.grid(True)

plt.tight_layout()
plt.show()

# Import all useful libraries
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# Randomly select 4 stocks from sminfra18
sminfra18_tickers = df_pivot_history_sminfra18.columns.tolist()

selected_stocks = np.random.choice(sminfra18_tickers, size=4, replace=False)
selected_stocks

array(['PGAS.JK', 'UNTR.JK', 'TBIG.JK', 'TLKM.JK'], dtype='<U7')

# Query the price info for the selected stocks
prices = df_pivot_history_sminfra18[selected_stocks]
prices.head()

# Calculate Returns
returns = pd.DataFrame()

for stock in prices:
  returns[stock + ' Returns'] = prices[stock].pct_change().dropna()

returns.head()

# Calculate expected returns (mean) and covariance matrix
# Assuming 252 is the trading days in a year = 365 - weekends 104 - holidays 9
# We need to annualize the return and risk since MVO models are long-term orientated and long-term in economics perspectives refers to 12 months +

expected_returns = returns.mean() * 252  # Annualize by multiplying by trading days
covariance_matrix = returns.cov() * 252   # Annualize covariance

# Print the expected returns and covariance matrix
print("Expected Returns:")
print(expected_returns)
print("\nCovariance Matrix:")
print(covariance_matrix)

Expected Returns:
PGAS.JK Returns    0.394443
UNTR.JK Returns    0.582308
TBIG.JK Returns   -0.056850
TLKM.JK Returns   -1.207224
dtype: float64

Covariance Matrix:
                 PGAS.JK Returns  UNTR.JK Returns  TBIG.JK Returns  \
PGAS.JK Returns         0.048921         0.022824         0.002496   
UNTR.JK Returns         0.022824         0.074070         0.008815   
TBIG.JK Returns         0.002496         0.008815         0.067872   
TLKM.JK Returns         0.013206        -0.006651         0.018117   

                 TLKM.JK Returns  
PGAS.JK Returns         0.013206  
UNTR.JK Returns        -0.006651  
TBIG.JK Returns         0.018117  
TLKM.JK Returns         0.110532 

# Here we can see the distribution of each stock returns
returns.hist(figsize=(8,6), bins = 50)
plt.show()

# Define Risk Free rate
# Assume RF as 6.64% as the yield of 10 year Indonesian Government bond
# Source link https://tradingeconomics.com/indonesia/government-bond-yield
rf = 0.0664

# Monte Carlo simulation
#Single Run
np.random.seed(101)
print(selected_stocks)

# Generates an array of random numbers representing initial weights for each asset in the portfolio
weights = np.array(np.random.random(4))
print('\nRandom Weights')
print(weights)

# Normalizing the randomly generated weights to ensure they sum up to 1, representing a fully invested portfolio.
print('\nRebalanced Weights')
weights = weights / np.sum(weights)
print(weights)

# Calculating the portfolio’s expected return using the weighted average of individual asset returns.
# It multiplies the annual mean returns of each asset by its respective weight and aggregates them.
print('\nPortfolio Return')
portfolio_return = np.sum(returns.mean()*252*weights)
print(portfolio_return)

# Expected Volatility (Standard Deviation): Using the formula for portfolio volatility, it calculates the square root of the dot product of weights, the covariance matrix of asset returns (multiplied by 252 for annualization), and weights transpose.
# The square root of covariance matrix is taken because square root of variance is standard deviation or volatility.
print('\nPortfolio Volatility')
portfolio_volatility = np.sqrt(np.dot(weights.T, np.dot(returns.cov()*252, weights)))
print(portfolio_volatility)

# Sharpe Ratio measures the risk-adjusted return by subtracting the risk-free rate from the portfolio return and dividing by its volatility.
print('\nPortfolio Sharpe Ratio')
sharpe_ratio = (portfolio_return - rf) / portfolio_volatility
print(sharpe_ratio)

['PGAS.JK' 'UNTR.JK' 'TBIG.JK' 'TLKM.JK']

Random Weights
[0.51639863 0.57066759 0.02847423 0.17152166]

Rebalanced Weights
[0.40122278 0.44338777 0.02212343 0.13326603]

Portfolio Return
0.25430815938385554

Portfolio Volatility
0.18304080881645118

Portfolio Sharpe Ratio
1.0265916141808857

# Specify the number of simulated portfolios to generate (in our case, 5000)
num_ports = 5000

# A 2D array to store the randomly generated weights for each asset in each portfolio.
all_weights = np.zeros((num_ports, len(selected_stocks)))

# Create arrays to store portfolio returns, volatilities, and Sharpe ratios for each simulated portfolio.
ret_arr = np.zeros(num_ports)
vol_arr = np.zeros(num_ports)
sharp_arr = np.zeros(num_ports)

# Looping over each portfolio
for i in range(num_ports):

  # Generate random weights randomly assigning weights to each asset in the portfolio and rebalancing them to sum up to 1
  weights = np.array(np.random.random(4))
  weights = weights / np.sum(weights)
  all_weights[i, :] = weights

  #Expected return
  ret_arr[i] = np.sum((returns.mean() * 252 * weights))

  #Expected volatility
  vol_arr[i] = np.sqrt(np.dot(weights.T, np.dot(returns.cov() * 252, weights)))

  #Sharpe ratio
  sharp_arr[i] = (ret_arr[i] - rf) / vol_arr[i]

# Retrieve the maximum value from the sharp_arr array
max_sr = sharp_arr.max()

# Find the index of the portfolio with the maximum Sharpe ratio in the sharp_arr array.
max_sr_index = sharp_arr.argmax()

# Retrieve the weights of the assets in the portfolio corresponding to the index 4091 in the all_weights array.
opt_weights = all_weights[max_sr_index,:]

# Retrive the optimal expected returns and volaticity
optimal_return = ret_arr[max_sr_index]
optimal_volatility = vol_arr[max_sr_index]

# Print Result
print('Max Sharpe Ratio: ', max_sr)
print('Optimal Return: ', optimal_return)
print('Optimal Volatility: ', optimal_volatility)
print('Optimal Weights: ', opt_weights)
print('Stock List: ', selected_stocks)

Max Sharpe Ratio:  2.0346018098713943
Optimal Return:  0.4819911385679692
Optimal Volatility:  0.20426165776105273
Optimal Weights:  [0.47205837 0.51425067 0.01118475 0.00250621]
Stock List:  ['PGAS.JK' 'UNTR.JK' 'TBIG.JK' 'TLKM.JK']

plt.figure(figsize = (10,5))
plt.scatter(vol_arr, ret_arr, c = sharp_arr, cmap = 'plasma')
plt.colorbar(label = 'Sharpe Ratio')
plt.xlabel('Volatility', fontweight = 'bold')
plt.ylabel('Return', fontweight = 'bold')

plt.scatter(optimal_volatility, optimal_return, c='red', s=200, edgecolors='black', marker='*')
plt.grid(True, ls=':', lw=1)

def ret_vol_sr(weights):
  weights = np.array(weights)

  # Calculate Annualized Expected Returns
  ret = np.sum(returns.mean() * weights * 252)

  # Calculate Portfolio Volaticity
  vol = np.sqrt(np.dot(weights.T, np.dot(returns.cov() * 252, weights)))

  # Calculate Sharp Ratio, Rist Free rate 6.4%
  sr = (ret - rf) / vol

  return np.array([ret, vol, sr])

from scipy.optimize import minimize

# This function below is designed to negate the Sharpe ratio because in optimization, the objective function is minimized by default.
# We want to maximize the Sharpe ratio, index 2, so multiplying by -1 helps flip the optimization to achieve that.
def sharp_neg(weights):
  return ret_vol_sr(weights)[2] * -1

# constraint function ensures the weight sum is 1
# In optimization, constraints are typically written such that the result equals 0
# If the sum of weights is exactly 1, then np.sum(weights) - 1 will equal 0. This satisfies the constraint that the sum of the weights must be 1.
def sum_check(weights):
  return np.sum(weights) - 1

# We want equity constraint, the function must return a value equals to 0
# Specify the equality constraint that the sum of weights must equal 1.
cons = ({'type': 'eq', 'fun': sum_check})

# Ensure that each weight lies between 0 and 1, meaning no short-selling or leverage.
bounds = ((0, 1), (0, 1), (0, 1), (0, 1))

# Initial guess for the optimization equal weighting of 25% for each stock initially.
init_guess = [0.25, 0.25, 0.25, 0.25]

# SLSQP stands for Sequential Least Squares Quadratic Programming, which is an algorithm that can handle both equality and inequality constraints, making it well-suited for portfolio optimization problems.
# We are minimize the negative sharp ratio - maximize sharp ratio

opt_results = minimize(sharp_neg, init_guess, method = 'SLSQP', bounds = bounds, constraints = cons)

opt_results

message: Optimization terminated successfully
 success: True
  status: 0
     fun: -2.0678015196575585
       x: [ 4.136e-01  5.864e-01  0.000e+00  1.062e-14]
     nit: 6
     jac: [-3.133e-01 -3.133e-01  5.538e-01  5.769e+00]
    nfev: 31
    njev: 6

# opt_results.x gives the optimized portfolio weights.
opt_results.x

array([4.13647417e-01, 5.86352583e-01, 0.00000000e+00, 1.06241346e-14])

#return the expected return, volatility, and Sharpe ratio for the optimized portfolio.
ret_vol_sr(opt_results.x)

array([0.50459828, 0.21191506, 2.06780152])

# Print Result
print('Return: ', ret_vol_sr(opt_results.x)[0])
print('Volatility: ', ret_vol_sr(opt_results.x)[1])
print('Sharpe Ratio: ', ret_vol_sr(opt_results.x)[2])
print (selected_stocks)
print ('Optimal Weights: ', opt_results.x)

Return:  0.5045982793617794
Volatility:  0.21191505818911857
Sharpe Ratio:  2.0678015196575585
['PGAS.JK' 'UNTR.JK' 'TBIG.JK' 'TLKM.JK']
Optimal Weights:  [4.13647417e-01 5.86352583e-01 0.00000000e+00 1.06241346e-14]

# Print Result from Monte Carlos Simulation
print ('\nResult from Monte Carlos Simulation')
print('Stock List: ', selected_stocks)
print('Optimal Weights: ', opt_weights)
print('Optimal Return: ', optimal_return)
print('Optimal Volatility: ', optimal_volatility)
print('Max Sharpe Ratio: ', max_sr)

# Print Result from Scipy Minimized Function
print ('\nResult from Scipy Minimized Function')
print ('Stock List: ', selected_stocks)
print ('Optimal Weights: ', opt_results.x)
print('Return: ', ret_vol_sr(opt_results.x)[0])
print('Volatility: ', ret_vol_sr(opt_results.x)[1])
print('Sharpe Ratio: ', ret_vol_sr(opt_results.x)[2])

Result from Monte Carlos Simulation
Stock List:  ['PGAS.JK' 'UNTR.JK' 'TBIG.JK' 'TLKM.JK']
Optimal Weights:  [0.47205837 0.51425067 0.01118475 0.00250621]
Optimal Return:  0.4819911385679692
Optimal Volatility:  0.20426165776105273
Max Sharpe Ratio:  2.0346018098713943

Result from Scipy Minimized Function
Stock List:  ['PGAS.JK' 'UNTR.JK' 'TBIG.JK' 'TLKM.JK']
Optimal Weights:  [4.13647417e-01 5.86352583e-01 0.00000000e+00 1.06241346e-14]
Return:  0.5045982793617794
Volatility:  0.21191505818911857
Sharpe Ratio:  2.0678015196575585

# Draw Efficient Frontier

plt.figure(figsize = (10,5))
plt.scatter(vol_arr, ret_arr, c = sharp_arr, cmap = 'plasma')
plt.colorbar(label = 'Sharpe Ratio')
plt.xlabel('Volatility', fontweight = 'bold')
plt.ylabel('Return', fontweight = 'bold')

# Optimial point for Monte Carlos Simulation is labelled as red star
plt.scatter(optimal_volatility, optimal_return, c='red', s=200, edgecolors='black', marker='*')
plt.grid(True, ls=':', lw=1)

# Optimial point for Scipy Optimization is labelled as green star
plt.scatter(ret_vol_sr(opt_results.x)[1], ret_vol_sr(opt_results.x)[0], c='green', s=200, edgecolors='black', marker='*')