<jupyter_start><jupyter_text>Porfolio Optimization  Section 1Choose 10 stocks you would like to analyze, at least 30 monthly observations, and download historical prices<jupyter_code># Import required libraries.
import yfinance as yf # for querying stock data
import pandas as pd # process CSV and manipulations
import numpy as np # for numerical computing
import quandl # for finance yield data
from pprint import pprint # pretty print
import seaborn as sn # for visualizations
import matplotlib.pyplot as plt # for plotting
import statsmodels.api as sm # for statistical models

pd.set_option('display.max_rows', None) # set dataframe display options to display full data content
quandl.ApiConfig.api_key = '<KEY>' # here we are using Quadl Python API for some data query tasks<jupyter_output><empty_output><jupyter_text>Selected stocks for observationsHere we are choosing the following 10 stocks' price data for our research task:| Stock Name                      	| NASDAQ/NYSE Code 	||---------------------------------	|------------------	|| 3M Co                           	|              MMM 	|| Netflix Inc                     	|             NFLX 	|| Amazon.com, Inc.                	|             AMZN 	|| Tesla Inc                       	|             TSLA 	|| Apple Inc                       	|             AAPL 	|| Facebook, Inc. Common Stock     	|               FB 	|| Electronic Arts Inc.            	|               EA 	|| Alibaba Group Holding Ltd - ADR 	|             BABA 	|| GrubHub Inc                     	|             GRUB 	|| Yelp Inc                        	|             YELP 	|<jupyter_code>stocks = ['MMM', 'NFLX', 'AMZN', 'TSLA', 'AAPL', 'FB', 'EA', 'BABA', 'GRUB', 'YELP'] # 10 selected stocks
# stocks_df = pd.DataFrame(columns=stocks) # dataframe to hold stock historical data
# Download historical stock price data with over 30 monthly observations
# for stock in stocks:
#     data = yf.download(stock, '2018-01-01', '2020-07-14', interval='1mo')['Adj Close']
#     stocks_df[stock] = data.dropna()

stocks_df = pd.read_csv('stock_data.csv')
stocks_df.set_index('Date', inplace=True)<jupyter_output><empty_output><jupyter_text>Stock Data<jupyter_code>display(stocks_df)
# stocks_df.to_csv('stocks.csv') # save the stock data into a csv file.<jupyter_output><empty_output><jupyter_text>Section 2Calculate the historical average returns, standard deviations and correlations.<jupyter_code># Define a function to compute monthly returns.
def compute_monthly_returns(previous: float, current: float)-> float:
    return current / previous - 1<jupyter_output><empty_output><jupyter_text>Stock Historical Average Returns and Standard Deviations<jupyter_code># Compute historical average returns, standard deviations and correlations for each stock.
monthly_returns = [] # store the result in a nested list.
for stock in stocks:
    row = []
    stock_arr = stocks_df[stock].values # get the stock price np array
    previous = stock_arr[0]
    for i in range(1, stock_arr.size):
        current = stock_arr[i]
        row.append(compute_monthly_returns(previous, current))
        previous = current
    monthly_returns.append(row)

monthly_returns_df = pd.DataFrame(list(map(list, zip(*monthly_returns))), columns=stocks)
display(monthly_returns_df)
print("Historical Average Returns:")
display(monthly_returns_df.mean())
print("Standard Deviation of the Returns:")
display(monthly_returns_df.std())<jupyter_output>Historical Average Returns:<jupyter_text>Stock Correlations<jupyter_code>monthly_returns_df.corr()<jupyter_output><empty_output><jupyter_text>Visualization of the Correlation Matrix<jupyter_code>correlation_matrix = stocks_df.corr()
plt.figure(figsize = (16,8))
sn.heatmap(correlation_matrix, annot=True, fmt='g', linewidths=.5)
plt.show()<jupyter_output><empty_output><jupyter_text>Section 3Build a covariance matrix.<jupyter_code>covariance_matrix = monthly_returns_df.cov()
display(covariance_matrix)<jupyter_output><empty_output><jupyter_text>Visualization of the Covariance Matrix<jupyter_code>plt.figure(figsize = (16,8))
sn.heatmap(covariance_matrix, annot=True, fmt='g', linewidths=.5)
plt.show()<jupyter_output><empty_output><jupyter_text>Section 4Look up the one month interest rate from Treasury bill quotes (make sure you convert to monthly)<jupyter_code>monthly_interest_rate = 0.17 / 12<jupyter_output><empty_output><jupyter_text>Question 5<jupyter_code># sp = yf.download('spy', '2018-01-01', '2020-07-14', interval='1mo').dropna()<jupyter_output><empty_output><jupyter_text>S&P 500 Index Data Using `spy` Agent<jupyter_code># sp_df = pd.DataFrame(sp['Adj Close'])
# display(sp_df)
sp_df = pd.read_csv('sp.csv')
sp_df.set_index('Date', inplace=True)
display(sp_df)<jupyter_output><empty_output><jupyter_text>Monthly Returns of Index Data<jupyter_code>sp_monthly_returns = []
previous = sp_df['Adj Close'].values[0]
for i in range(1, len(sp_df['Adj Close'].values)):
    current = sp_df['Adj Close'].values[i]
    sp_monthly_returns.append(compute_monthly_returns(previous, current))
    previous = current
pprint(sp_monthly_returns)<jupyter_output>[-0.03636024193828624,
 -0.03129011897430245,
 0.009193838163753076,
 0.024309121008824208,
 0.0012548519155881621,
 0.04170329851153287,
 0.03191970569435809,
 0.0014124522348848423,
 -0.06489033285677448,
 0.018549222701129464,
 -0.09334303153693013,
 0.0863730209435849,
 0.032415579179747134,
 0.013636087444746092,
 0.04543689264241135,
 -0.06377122103707777,
 0.0644093853844665,
 0.020056785373939068,
 -0.016743349021436593,
 0.014771730900268176,
 0.026824706297821743,
 0.03619811100853454,
 0.02402092196012262,
 0.00451066716903914,
 -0.07916569197267742,
 -0.12998723151340352,
 0.13361053327894612,
 0.04764525713784029,
 0.013275427295073028,
 0.03446486952874084,
 -0.008753420483494034]<jupyter_text>Porfolio Optimization Using Monte Carlo Simulation<jupyter_code>num_portfolios = 10000
log_monthly_returns = np.log(stocks_df / stocks_df.shift(1)) # use logarithm for better numerical performance
weights_tangency_porfolio = np.zeros((num_portfolios, len(stocks)))
return_arr = np.zeros(num_portfolios)
std_arr = np.zeros(num_portfolios)
sharpe_arr = np.zeros(num_portfolios)

for portfolio in range(num_portfolios):
    # Get current random weights.
    weights = np.random.randn(len(stocks))
    weights /= np.sum(weights)
    weights_tangency_porfolio[portfolio, :] = weights
    
    # Compute current expected return.
    return_arr[portfolio] = np.sum((np.dot(log_monthly_returns.mean(), weights)) * 12)
    
    # Compute current standard deviation.
    std_arr[portfolio] = np.sqrt(np.dot(weights.T, np.dot(log_monthly_returns.cov() * 12, weights)))
    
    # Sharpe Ratio.
    sharpe_arr[portfolio] = return_arr[portfolio] / std_arr[portfolio]
optimal_weights = weights_tangency_porfolio[sharpe_arr.argmax(),:]
expected_return = return_arr[sharpe_arr.argmax()]
standard_deviation = std_arr[sharpe_arr.argmax()]
optimal_sharpe_ratio = sharpe_arr.max()

print('Weights of the tangency portfolio:')
print(optimal_weights)
print(f'Expected return: {expected_return}')
print(f'Standard Deviation: {standard_deviation}')
print(f'Sharpe Ratio: {optimal_sharpe_ratio}')<jupyter_output>Weights of the tangency portfolio:
[-9.03872035  2.49603837  5.41679138  0.6817262   7.67141131 -4.410834
 -2.47610113  3.43687883 -2.87640417  0.09921357]
Expected return: 6.394356159892053
Standard Deviation: 3.076721089707106
Sharpe Ratio: 2.0783021838683386<jupyter_text>Generate Visualizations of the Markowitz Bullet.<jupyter_code>plt.figure(figsize=(12, 6))
plt.scatter(std_arr, return_arr, c=sharpe_arr, cmap='plasma')
plt.scatter(standard_deviation, expected_return, c='red', s=60, edgecolors='black')
plt.xlim(-5, 60)
plt.ylim(-20, 20)
plt.colorbar(label='Sharpe Ratio')
plt.xlabel('Standard Deviation')
plt.ylabel('Return')<jupyter_output><empty_output><jupyter_text>This is the result I get using Monte Carlo Simulation, which appears to be a better result from our Excel version, however, we cannot achive the similarly high Sharpe Ratio using Excel, hence I'm keeping this result here.  Question 6Find the implied betas, using the tangency portfolio as the “market”.<jupyter_code># Compute the implied betas.
avg_returns = monthly_returns_df.mean() # compute average monthly return for each stock

betas = dict()
for stock in stocks:
    beta = (avg_returns[stock] - monthly_interest_rate) / (expected_return - monthly_interest_rate)
    betas[stock] = beta

pprint(betas)<jupyter_output>{'AAPL': 0.0028210990961154096,
 'AMZN': 0.002279686663194948,
 'BABA': -0.00045573967583096355,
 'EA': -0.0013266196007789843,
 'FB': -0.00030882622998449445,
 'GRUB': -0.00015453196885890265,
 'MMM': -0.0038930902107982145,
 'NFLX': 0.0019149017081609547,
 'TSLA': 0.008108014105713867,
 'YELP': -0.00405126183866965}<jupyter_text>Question 7Also estimate historical betas using the S&P 500 as the “market”1. You can use SPX or SPY as proxies for the S&P 5002. Using the market risk premium of 6%, calculate expected returns with the CAPM  Estimate the Historical Beta ValuesThe historical betas can be estimated as the closed form OLS solution:$$\mathbf{\beta} = (\mathbf{X'}\mathbf{X})^{-1}\mathbf{X'}\mathbf{Y},$$where $\mathbf{X}$ is the Index Return and the $\mathbf{Y}$ is the Stock Return.<jupyter_code>stock_betas = dict()
Y = np.array(monthly_returns_df.values)
X = np.array(sp_monthly_returns)
XXX = np.dot(1/np.dot(X.T, X), X.T)
beta = np.dot(XXX, Y)
historical_betas = np.column_stack((stocks, beta))
historical_betas_df = pd.DataFrame(historical_betas, columns=['NASDAQ/NYSE CODE', 'Estimated Beta'])
historical_betas_df.set_index('NASDAQ/NYSE CODE', inplace=True)
display(historical_betas_df)
# NOTE: the beta can also be computed using OLS regression model provided by statsmodel library, take 3M for example:

X = sp_monthly_returns
Y = monthly_returns_df['MMM']
# X1 = sm.add_constant(X)

# make regression model 
model = sm.OLS(Y, X)

# fit model and print results
results = model.fit()
print(results.summary())<jupyter_output>OLS Regression Results                                
=======================================================================================
Dep. Variable:                    MMM   R-squared (uncentered):                   0.451
Model:                            OLS   Adj. R-squared (uncentered):              0.433
Method:                 Least Squares   F-statistic:                              24.63
Date:                Sat, 18 Jul 2020   Prob (F-statistic):                    2.59e-05
Time:                        23:32:20   Log-Likelihood:                          49.036
No. Observations:                  31   AIC:                                     -96.07
Df Residuals:                      30   BIC:                                     -94.64
Df Model:                           1                                                  
Covariance Type:            nonrobust                                                  
================================[...]<jupyter_text>Compute the Expected Returns<jupyter_code># Compute expected returns in percentage form.
expected_return_dict = dict()
risk_premium = 0.06
for stock in stocks:
    beta = float(historical_betas_df['Estimated Beta'][stock])
    expected_return_dict[stock] = (monthly_interest_rate + beta * risk_premium) * 100

expected_return_df = pd.DataFrame(expected_return_dict.items(), columns=['NASDAQ/NYSE CODE', 'Expected Return'])
expected_return_df.set_index('NASDAQ/NYSE CODE', inplace=True)
pd.options.display.float_format = '{:,.2f}%'.format
display(expected_return_df)<jupyter_output><empty_output>