montecarlo simulation
import pandas as pd

import numpy as np

import matplotlib.pyplot as plt

from scipy.optimize import minimize

Load the dataset

file_path = '/Users//Downloads/F-F_Research_Data_Factors-2.xlsx'

df = pd.read_excel(file_path)

Convert 'Date' to datetime type and filter for the last 30 years

df['Date'] = pd.to_datetime(df['Date'], format='%m/%Y')

df_filtered = df[df['Date'] >= (df['Date'].max() - pd.DateOffset(years=30))]

Extracting the factors for analysis

factors = ['Mkt-RF', 'SMB', 'HML', 'Mom', 'LT_Rev']

returns = df_filtered[factors]

Normalizing the returns (assuming they are in percentage format)

returns_normalized = returns / 100

Monte Carlo Simulation for Portfolio Optimization

num_portfolios = 30000

all_weights = np.zeros((num_portfolios, len(factors)))

ret_arr = np.zeros(num_portfolios)

vol_arr = np.zeros(num_portfolios)

np.random.seed(42) # For reproducibility

for i in range(num_portfolios):

weights = np.random.random(len(factors))

weights /= np.sum(weights)

all_weights[i, :] = weights

ret_arr[i] = np.sum((returns_normalized.mean() * weights * 12))  # Annualizing returns

vol_arr[i] = np.sqrt(np.dot(weights.T, np.dot(returns_normalized.cov() * 12, weights)))  # Annualizing standard deviation

Identifying the portfolio with minimum variance

min_vol_idx = np.argmin(vol_arr)

min_vol_weights = all_weights[min_vol_idx]

min_vol_portfolio_return = ret_arr[min_vol_idx]

min_vol_portfolio_volatility = vol_arr[min_vol_idx]

Plotting the portfolios

plt.figure(figsize=(12, 8))

plt.scatter(vol_arr, ret_arr, c=ret_arr/vol_arr, marker='o')

plt.colorbar(label='Sharpe Ratio')

plt.xlabel('Expected Volatility')

plt.ylabel('Expected Return')

plt.scatter(min_vol_portfolio_volatility, min_vol_portfolio_return, c='red', s=50, marker='*', label='Minimum Variance Portfolio')

plt.legend()

plt.title('Monte Carlo Simulation of Portfolios')

plt.show()

Outputting the minimum variance portfolio details

min_vol_portfolio_details = {

'Return': min_vol_portfolio_return,

'Volatility': min_vol_portfolio_volatility,

'Weights': dict(zip(factors, min_vol_weights))

}

print(min_vol_portfolio_details)

Continuing from the previous analysis, we will now find portfolios that satisfy the additional criteria:

2. Maximum Sharpe Ratio

3. Maximum Return

4. Minimum Value-at-Risk (VaR) at 1%

Calculating the Sharpe Ratio for each portfolio

risk_free_rate = 0 # Assuming a risk-free rate of 0 for simplicity

sharpe_ratios = ret_arr / vol_arr # Sharpe Ratio = Return / Volatility

Finding the portfolio with the maximum Sharpe Ratio

max_sharpe_idx = np.argmax(sharpe_ratios)

max_sharpe_portfolio = {

'Return': ret_arr[max_sharpe_idx],

'Volatility': vol_arr[max_sharpe_idx],

'Sharpe Ratio': sharpe_ratios[max_sharpe_idx],

'Weights': all_weights[max_sharpe_idx]

}

Finding the portfolio with the maximum Return

max_ret_idx = np.argmax(ret_arr)

max_ret_portfolio = {

'Return': ret_arr[max_ret_idx],

'Volatility': vol_arr[max_ret_idx],

'Weights': all_weights[max_ret_idx]

}

Value-at-Risk (VaR) Calculation

Assuming normal distribution of returns, we can calculate VaR using the Z-score for the 1% confidence interval

from scipy.stats import norm

Z-score for 1% confidence interval

z_score_1pc = norm.ppf(1 - 0.01)

Calculating VaR for each portfolio

var_1pc = -z_score_1pc * vol_arr # Negative sign because VaR is typically expressed as a positive number

Finding the portfolio with the minimum VaR

min_var_idx = np.argmin(var_1pc)

min_var_portfolio = {

'Return': ret_arr[min_var_idx],

'Volatility': vol_arr[min_var_idx],

'VaR 1%': var_1pc[min_var_idx],

'Weights': all_weights[min_var_idx]

}