!pip install yfinance
pip install openpyxl
!pip install statsmodels==0.13.0
!pip install lazypredict
import pandas as pd
import numpy as np
import yfinance as yf
import datetime as dt
import matplotlib.pyplot as plt
import seaborn as sns
import statsmodels.api as sm
from statsmodels.formula.api import ols, logit
from lazypredict.Supervised import LazyClassifier, LazyRegressor
from sklearn.model_selection import train_test_split
# set matplotlib.pyplot to default style
plt.style.use('default')
# Random selected companies from 4 industry
Manufacturing = ["TSM", "TXN", "AMAT", "IRBT", "LRCX", "PATK", "ACCO", "AIN", "EL", "ECL", "ZBRA", "RRX", "TSLA", "DOV", "KMB", "JBL", "FLEX", "ENTG", "ACH", "HELE", "MHK", "CAJ", "JJSF", "LULU", "CHD"] #25 companies
Medical = ["EW", "WST", "IDXX", "RMD","VRTX","AMGN", "ALGN", "PFE", "ZTS", "DHR", "TMO", "CERN", "UNH", "COO", "CNC", "QGEN", "CAH", "HCA", "MDT", "MRK", "JNJ", "REGN", "TFX", "LLY", "DVA", "ABC"] #26 companies
Logistic_transportation= ["KSU", "CHRW", "EXPD", "UNP", "FDX", "TIAIY", "ORAN", "EIX", "ASR", "DAL", "KNX", "SKYW", "HUBG", "HTLD", "JBHT", "OKE", "KEX", "NCLH", "ODFL", "PAC", "ZNH", "UAL", "TDW"]#23 companies #"TIGO","OSG",
Energy = ['NEE','ORA','ES','SRE','BWXT', "AEP", "DUK", "ATO", "PAGP", "DVN", "VLO", "PDCE", "HLX", "RRC", "KOS", "FRO", "TRGP", "LIN", "HE", "ENLAY", "CPE", "EC", "WLL", "CRK"]#24 compenies in use, excl. 'PLUG','ENPH'"BTU",
# Stock return
start=dt.datetime(2014,12,31)
end=dt.datetime(2021,9,30)
overall = []
overall.extend(Manufacturing)
overall.extend(Medical)
overall.extend(Logistic_transportation)
overall.extend(Energy)
data_O=yf.download(overall,start=start,end=end)
df_O_all=data_O['Adj Close']
df_O_all=df_O_all.ffill()
# we used 5 year data from 2014 to 2019. The return for 2020 is too volatile due to pandemic
df_O = df_O_all[:dt.datetime(2019,12,31)]
# Benchmark data (market indices & ESG ETFs)
benchmark_df = yf.download(['SPY'], start=start, end=end)
benchmark_SPY_all = benchmark_df['Adj Close']
benchmark_SPY = benchmark_SPY_all[:dt.datetime(2019,12,31)]
benchmark_df = yf.download(['^IXIC'], start=start, end=end)
benchmark_IXIC_all = benchmark_df['Adj Close']
benchmark_IXIC = benchmark_IXIC_all[:dt.datetime(2019,12,31)]
benchmark_df = yf.download(['DJI'], start=start, end=end)
benchmark_DJI_all = benchmark_df['Adj Close']
benchmark_DJI = benchmark_DJI_all[:dt.datetime(2019,12,31)]
benchmark_df = yf.download(['ESGD'], start=start, end=end)
benchmark_ESGD_all = benchmark_df['Adj Close']
benchmark_ESG = benchmark_ESGD_all[:dt.datetime(2019,12,31)]
benchmark_df = yf.download(['ESGV'], start=start, end=end)
benchmark_ESGV_all = benchmark_df['Adj Close']
benchmark_ESG2 = benchmark_ESGV_all[:dt.datetime(2019,12,31)]
benchmark_df = yf.download(['ICLN'], start=start, end=end)
benchmark_ICLN_all = benchmark_df['Adj Close']
benchmark_ESG3 = benchmark_ICLN_all[:dt.datetime(2019,12,31)]
benchmark_df = yf.download(['SUSA'], start=start, end=end)
benchmark_SUSA_all = benchmark_df['Adj Close']
benchmark_ESG4 = benchmark_SUSA_all[:dt.datetime(2019,12,31)]
benchmark_df = yf.download(['ESGHX'], start=start, end=end)
benchmark_ESGHX_all = benchmark_df['Adj Close']
benchmark_ESG5 = benchmark_ESGHX_all[:dt.datetime(2019,12,31)]
benchmark_df = yf.download(['NEXTX'], start=start, end=end)
benchmark_NEXTX_all = benchmark_df['Adj Close']
benchmark_ESG6 = benchmark_NEXTX_all[:dt.datetime(2019,12,31)]
benchmark_df = yf.download(['ATEYX'], start=start, end=end)
benchmark_ATEYX_all = benchmark_df['Adj Close']
benchmark_ESG7 = benchmark_ATEYX_all[:dt.datetime(2019,12,31)]
# ESG Benchmark
_5benchmark =['ICLN','SUSA','ESGHX','NEXTX','ATEYX']
# market indices benchamrk
_3benchmark = ['SPY', 'IXIC','DJI']
pip install xlrd
#Read and clean data in Fama-French 3
path_1 = "FF3.CSV"
FF3=pd.read_csv(path_1, skiprows=4)
#FF3 = FF3.sort_index(ascending = False)
#FF3 = FF3.reset_index()
FF3 = FF3.rename(columns={'Unnamed: 0': 'Year'})
FF3 = FF3[:-1]
FF3["Year"] = pd.to_datetime(FF3['Year'],format='%Y%m%d')
FF3 = FF3[FF3['Year'] >= dt.datetime(2014,12,31)]
df = pd.DataFrame()
for i in overall:
path = "Archive/ESG Table for " + i + ".xlsx"
temp_n=pd.read_excel(path,engine='openpyxl')
temp_n = temp_n[1:14]
temp_n = temp_n.drop(temp_n.index[1:7])
#transpose the column
temp_n = temp_n.T
temp_n.columns = temp_n.iloc[0]
temp_n = temp_n.drop(temp_n.index[0])
# Rename the rows
temp_n.set_index(temp_n.columns[0],inplace=True)
temp_n = temp_n[1:6]
temp_n = temp_n.reset_index()
# Rename column names
temp_n = temp_n.rename(columns={np.nan: 'Year'})
temp_n.columns = ['Year','ESG_C_S', 'ESG_S', 'Env_S', 'Soc_S', 'Gov_S', 'ESG_Controversial']
# Rename the score into integer
temp_n = temp_n.replace(["D-", "D","D+","C-", "C","C+","B-", "B","B+","A-", "A","A+"], [0,1,2,3,4,5,6,7,8,9,10,11])
temp_n['Year'] = pd.to_datetime(temp_n['Year'])
temp_n = temp_n.set_index("Year")
new_date1 = pd.DataFrame(temp_n[-1:].values, index=[pd.to_datetime('2014-12-31')], columns=temp_n.columns)
new_date2 = pd.DataFrame(temp_n[-1:].values, index=[pd.to_datetime('2019-12-31')], columns=temp_n.columns)
temp_n = temp_n.append(new_date1)
temp_n = temp_n.append(new_date2)
temp_n = temp_n.sort_index()
temp_n = temp_n.resample('D').asfreq()
temp_n = temp_n.fillna(method = 'ffill')
# adjust column type
#temp_n["Year"] = pd.to_numeric(temp_n["Year"])
temp_n = temp_n.apply(pd.to_numeric)
#Merge two files
temp_n = temp_n.merge(FF3, how = "inner", left_on = temp_n.index, right_on = FF3['Year'])
temp_n = temp_n.set_index("Year")
temp_n = temp_n.drop(columns = 'key_0')
# Compute Daily Return
temp_return = df_O[i].ffill().pct_change()
temp_return = temp_return.to_frame()
temp_return.columns = [i]
temp_return = temp_return[1:]
temp_return = temp_return.round(3)
temp_return = temp_return.reset_index()
temp_return = temp_return.rename(columns={'Date': 'Year', i: "Ret"})
temp_return['Year'] = pd.DatetimeIndex(temp_return['Year'])
temp_n = temp_n.merge(temp_return, how = "outer", on = "Year")
#Set the company name as index
temp_n["Company"] = i
if i in Manufacturing:
temp_n['Industry'] = 'Manufacturing'
elif i in Medical:
temp_n['Industry'] = 'Medical'
elif i in Logistic_transportation:
temp_n['Industry'] = 'Logistic'
else:
temp_n['Industry'] = 'Energy'
temp_n = temp_n.set_index("Company")
df = pd.concat([df,temp_n], axis = 0)
df_I = df
# get dummy variable for Industry
df = pd.get_dummies(df, prefix=['I'], columns=['Industry'],drop_first=True)
df = df.dropna()
df
df.isna().sum()
# plot regress plot to see if there is any trend in the data
plt.figure(figsize = (16,8))
sns.regplot(x = 'ESG_C_S', y = "Ret", data = df, logistic=True, scatter_kws={"color": "black"}, y_jitter=.05,line_kws={"color": "red"})
plt.ylabel("return")
plt.xlabel("ESG index")
def CI_test(df, i):
bs_replicates = np.empty(10000)
for j in range(10000):
bs_replicates[j] = np.mean(np.random.choice(df[i], len(df[i])))
# Bootstrapped confidence intervals
conf_int = np.percentile(bs_replicates,[2.5, 97.5])
print('95% bootstrapped confidence interval =', conf_int, '%',i)
# Theoretical confidence intervals
conf_int_actual_upper = df[i].mean() + df[i].std()/np.sqrt(df[i].count())*1.96
conf_int_actual_lower = df[i].mean() - df[i].std()/np.sqrt(df[i].count())*1.96
conf_int_actual = [conf_int_actual_lower, conf_int_actual_upper]
print('-'*120)
print('95% theoretical confidence interval =', conf_int_actual, '%',i)
CI_test(df, 'Ret')
# remove Year as it contain multicollinearity problem
df = df.drop(columns = 'Year')
df.columns
# Correlation heatmap --> identify potential correlated factors
corr = df.corr()
plt.figure(figsize = (16,8))
sns.heatmap(corr,annot=True)
plt.show()
# dependent/target/outcome variable
y = df['Ret']
# independent/predictor/explanatory variable
X = df[['ESG_C_S','ESG_S','Env_S', 'Soc_S','Gov_S','ESG_Controversial', 'Mkt-RF','SMB', 'RF' ,'HML','I_Logistic', 'I_Manufacturing', 'I_Medical']]
#
#'ESG_C_S', 'ESG_S', 'Env_S', 'Soc_S','Gov_S','ESG_Controversial'
#'Mkt-RF','SMB', 'RF' ,'HML','I_Logistic', 'I_Manufacturing', 'I_Medical',
# OLS Regression
X = sm.add_constant(X)
ols_model = sm.OLS(y,X.astype(float))
result = ols_model.fit()
print(result.summary2())
#conditional No. should be lower than 20
df_O.index
df_O = df_O.drop(df_O.index[0] ,axis = 0)
df_1 = df_O.fillna(method = 'ffill')
df_1 = df_1.pct_change()
df_1 = df_1.drop(df_1.index[0] ,axis = 0)
df_1
# Stationary test
from statsmodels.tsa.stattools import adfuller
def test_stationarity(timeseries, i):
# Determining rolling statistics
rolmean = timeseries.rolling(14).mean()
rolstd = timeseries.rolling(14).std()
# Plot rolling statistics:
#orig = plt.plot(timeseries, color='blue',label='Original')
#mean = plt.plot(rolmean, color='red', label='Rolling Mean')
#std = plt.plot(rolstd, color='black', label = 'Rolling Std')
#plt.legend(loc='best')
#plt.title('Rolling Mean & Standard Deviation')
#plt.show(block=False)
# Perform Dickey-Fuller test:
#print ('Results of Dickey-Fuller Test:' + i)
timeseries = timeseries.iloc[:].values
dftest = adfuller(timeseries, autolag='AIC')
#dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
#for key,value in dftest[4].items():
# dfoutput['Critical Value (%s)'%key] = value
#print(dfoutput)
if dftest[1] > 0.05:
print(i)
for i in overall:
test_stationarity(df_1[i], i)
# load data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.2, random_state=42)
#LazyPredict on regression models
# fit all models
reg = LazyRegressor(predictions=True)
models, predictions = reg.fit(X_train, X_test, y_train, y_test)
models
predictions
!pip install PyPortfolioOpt==1.2.1
pip install --upgrade scipy
from pypfopt.efficient_frontier import EfficientFrontier
from pypfopt import risk_models
from pypfopt import expected_returns
from pypfopt import cla
from pypfopt.plotting import Plotting
from pypfopt.discrete_allocation import DiscreteAllocation, get_latest_prices
from pypfopt import objective_functions
import scipy.stats as st
# compute efficient frontier & portfolio statistics
class EF:
def __init__(self, asset, asset_list):
self.asset = asset
self.asset_list = asset_list
self.computation(asset, asset_list)
def computation(self, asset, asset_list):
# Check NaN values in the data
nullin_df = pd.DataFrame(self.asset,columns=self.asset_list)
# Calculate portfolio mean return
mu = expected_returns.mean_historical_return(self.asset)
# Calculate portfolio return variance
sigma = risk_models.sample_cov(self.asset)
# weight bounds in negative allows shorting of stocks
ef = EfficientFrontier(mu, sigma, weight_bounds=(0,1))
ef.add_objective(objective_functions.L2_reg, gamma=2)
# optional constraints possible, read pypfopt documentation.
sharpe_portfolio=ef.max_sharpe(risk_free_rate=0.00078) #US Bond
sharpe_portfolio_wt=ef.clean_weights()
latest_prices = get_latest_prices(self.asset)
# Allocate Portfolio Value in $ as required to show number of shares/stocks to buy,
# also bounds for shorting will affect allocation
# Maximum Sharpe Portfolio Allocation $1000000
da = DiscreteAllocation(sharpe_portfolio_wt, latest_prices, total_portfolio_value=1000000)
allocation, leftover = da.greedy_portfolio()
max_sharpe_cla = cla.CLA(mu, sigma)
max_sharpe_cla.max_sharpe()
sharpe_portfolio_wt_list = list(sharpe_portfolio_wt.values())
ret_data1 = self.asset.pct_change()[1:]
weighted_returns1 = (sharpe_portfolio_wt_list * ret_data1)
portfolio_ret1 = pd.DataFrame(weighted_returns1.sum(axis=1))
ret_data1 = ret_data1.merge(portfolio_ret1, on="Date", how="left")
ret_data1 = ret_data1.rename(columns={0: "portfolio_ret"})
ret_data1['cumulative_portfolio_ret'] = (ret_data1['portfolio_ret'] + 1).cumprod()
self.sharpe_portfolio_wt = sharpe_portfolio_wt
self.portfolio_ret1 = portfolio_ret1
self.allocation = allocation
self.max_sharpe_cla = max_sharpe_cla
self.ret_data1 = ret_data1
self.leftover = leftover
def SP_weight(self):
print(self.sharpe_portfolio_wt)
print(self.allocation)
print("Leftover Fund value for the maximum Sharpe portfolio is ${:.2f}".format(self.leftover))
P_weight = Plotting.plot_weights(self.sharpe_portfolio_wt)
return P_weight
def Efficient_frontier(self):
EF_plot = Plotting.plot_efficient_frontier(self.max_sharpe_cla, show_assets="True")
return EF_plot
def cumulative_return(self):
Cum_ret = sns.scatterplot('Date', 'cumulative_portfolio_ret', data=self.ret_data1, color = 'green')
print(self.ret_data1['cumulative_portfolio_ret'][-1])
return Cum_ret
def plot_hist(self):
hist_ret = plt.hist(self.ret_data1['portfolio_ret'], bins=10)
print(st.norm.interval(alpha=0.95, loc=np.mean(self.ret_data1['portfolio_ret']), scale=st.sem(self.ret_data1['portfolio_ret'])))
return hist_ret
def port_stat(self):
geometric_port_return = np.prod(self.portfolio_ret1 + 1) ** (252/self.portfolio_ret1.shape[0]) - 1
annual_std = np.std(self.portfolio_ret1) * np.sqrt(252)
port_sharpe_ratio = geometric_port_return / annual_std
print("Annual Return", geometric_port_return)
print("Stdev", annual_std)
print("Sharpe Ratio", port_sharpe_ratio)
def CI_test(self):
bs_replicates = np.empty(10000)
for i in range(10000):
bs_replicates[i] = np.mean(np.random.choice(self.ret_data1['portfolio_ret'], len(self.ret_data1['portfolio_ret'])))
# Bootstrapped confidence intervals
conf_int = np.percentile(bs_replicates,[2.5, 97.5])
print('95% bootstrapped confidence interval =', conf_int, '%')
# Theoretical confidence intervals
conf_int_actual_upper = self.ret_data1['portfolio_ret'].mean() + self.ret_data1['portfolio_ret'].std()/np.sqrt(self.ret_data1['portfolio_ret'].count())*1.96
conf_int_actual_lower = self.ret_data1['portfolio_ret'].mean() - self.ret_data1['portfolio_ret'].std()/np.sqrt(self.ret_data1['portfolio_ret'].count())*1.96
conf_int_actual = [conf_int_actual_lower, conf_int_actual_upper]
print('-'*120)
print('95% theoretical confidence interval =', conf_int_actual, '%')
# compute efficient frontier & portfolio statistics with benchmarks
class EF_benchmark5:
def __init__(self, asset, asset_list,benchmark1,benchmark2,benchmark3,benchmark4,benchmark5, benchmark_name = None):
self.asset = asset
self.asset_list = asset_list
self.benchmark1 = benchmark1
self.benchmark2 = benchmark2
self.benchmark3 = benchmark3
self.benchmark4 = benchmark4
self.benchmark5 = benchmark5
self.benchmark_name = benchmark_name
self.computation(asset, asset_list,benchmark1,benchmark2,benchmark3,benchmark4,benchmark5)
def computation(self, asset, asset_list,benchmark1,benchmark2,benchmark3,benchmark4,benchmark5):
# Check NaN values in the data
nullin_df = pd.DataFrame(self.asset,columns=self.asset_list)
# Calculate portfolio mean return
mu = expected_returns.mean_historical_return(self.asset)
# Calculate portfolio return variance
sigma = risk_models.sample_cov(self.asset)
# weight bounds in negative allows shorting of stocks
ef = EfficientFrontier(mu, sigma, weight_bounds=(0,1))
ef.add_objective(objective_functions.L2_reg, gamma=2)
# optional constraints possible, read pypfopt documentation.
sharpe_portfolio=ef.max_sharpe(risk_free_rate=0.00078) #US Bond
sharpe_portfolio_wt=ef.clean_weights()
latest_prices = get_latest_prices(self.asset)
# Allocate Portfolio Value in $ as required to show number of shares/stocks to buy,
# also bounds for shorting will affect allocation
# Maximum Sharpe Portfolio Allocation $1000000
da = DiscreteAllocation(sharpe_portfolio_wt, latest_prices, total_portfolio_value=1000000)
allocation, leftover = da.greedy_portfolio()
max_sharpe_cla = cla.CLA(mu, sigma)
max_sharpe_cla.max_sharpe()
sharpe_portfolio_wt_list = list(sharpe_portfolio_wt.values())
ret_data1 = self.asset.pct_change()[1:]
weighted_returns1 = (sharpe_portfolio_wt_list * ret_data1)
portfolio_ret1 = pd.DataFrame(weighted_returns1.sum(axis=1))
ret_data1 = ret_data1.merge(portfolio_ret1, on="Date", how="left")
ret_data1 = ret_data1.rename(columns={0: "portfolio_ret"})
ret_data1['cumulative_portfolio_ret'] = (ret_data1['portfolio_ret'] + 1).cumprod()
ret_benchmark1 = self.benchmark1.pct_change()[1:]
ret_benchmark1['cumulative_benchmark_ret'] = (ret_benchmark1 + 1).cumprod()
ret_benchmark2 = self.benchmark2.pct_change()[1:]
ret_benchmark2['cumulative_benchmark_ret'] = (ret_benchmark2 + 1).cumprod()
ret_benchmark3 = self.benchmark3.pct_change()[1:]
ret_benchmark3['cumulative_benchmark_ret'] = (ret_benchmark3 + 1).cumprod()
ret_benchmark4 = self.benchmark4.pct_change()[1:]
ret_benchmark4['cumulative_benchmark_ret'] = (ret_benchmark4 + 1).cumprod()
ret_benchmark5 = self.benchmark5.pct_change()[1:]
ret_benchmark5['cumulative_benchmark_ret'] = (ret_benchmark5 + 1).cumprod()
self.sharpe_portfolio_wt = sharpe_portfolio_wt
self.portfolio_ret1 = portfolio_ret1
self.allocation = allocation
self.max_sharpe_cla = max_sharpe_cla
self.ret_data1 = ret_data1
self.ret_benchmark1 = ret_benchmark1
self.ret_benchmark2 = ret_benchmark2
self.ret_benchmark3 = ret_benchmark3
self.ret_benchmark4 = ret_benchmark4
self.ret_benchmark5 = ret_benchmark5
self.leftover = leftover
def SP_weight(self):
print(self.sharpe_portfolio_wt)
print(self.allocation)
print("Leftover Fund value for the maximum Sharpe portfolio is ${:.2f}".format(self.leftover))
P_weight = Plotting.plot_weights(self.sharpe_portfolio_wt)
return P_weight
def Efficient_frontier(self):
EF_plot = Plotting.plot_efficient_frontier(self.max_sharpe_cla, show_assets="True")
return EF_plot
def cumulative_return(self):
Cum_ret = sns.scatterplot('Date', 'cumulative_portfolio_ret', data=self.ret_data1, color = 'yellowgreen', size = 1)
Cum_ret = sns.scatterplot('Date', 'cumulative_benchmark_ret', data=self.ret_benchmark1, color = 'olivedrab',size = 1)
Cum_ret = sns.scatterplot('Date', 'cumulative_benchmark_ret', data=self.ret_benchmark2, color = 'darkslategray',size = 1)
Cum_ret = sns.scatterplot('Date', 'cumulative_benchmark_ret', data=self.ret_benchmark3, color = 'darkseagreen',size = 1)
Cum_ret = sns.scatterplot('Date', 'cumulative_benchmark_ret', data=self.ret_benchmark4, color = 'darkgreen',size = 1)
Cum_ret = sns.scatterplot('Date', 'cumulative_benchmark_ret', data=self.ret_benchmark5, color = 'teal',size = 1)
if self.benchmark_name != None:
plt.legend(labels = ["ESG", *self.benchmark_name])
else:
plt.legend(labels = ["ESG", "Benchmark"])
print(self.ret_data1['cumulative_portfolio_ret'][-1])
return Cum_ret
def plot_hist(self):
hist_ret = plt.hist(self.ret_data1['portfolio_ret'], bins=10)
print(st.norm.interval(alpha=0.95, loc=np.mean(self.ret_data1['portfolio_ret']), scale=st.sem(self.ret_data1['portfolio_ret'])))
return hist_ret
def port_stat(self):
geometric_port_return = np.prod(self.portfolio_ret1 + 1) ** (252/self.portfolio_ret1.shape[0]) - 1
annual_std = np.std(self.portfolio_ret1) * np.sqrt(252)
port_sharpe_ratio = geometric_port_return / annual_std
print("Annual Return", geometric_port_return)
print("Stdev", annual_std)
print("Sharpe Ratio", port_sharpe_ratio)
def hypothesis_testing(self):
#import scipy.stats as stats
plt1 = self.ret_data1['portfolio_ret'].hist(bins=100, color='r', alpha=0.5, figsize = (16,8))
plt1 = self.benchmark2.pct_change()[1:].hist(bins=100, color='olivedrab', alpha=0.5)
print("Difference in mean return (" + self.benchmark_name[1] + ", ESG): ")
print((self.ret_data1['portfolio_ret'].mean() - self.benchmark2.pct_change()[1:].mean())*100)
stat, p = st.ttest_ind(self.ret_data1['portfolio_ret'], self.benchmark2.pct_change()[1:], equal_var=False)#, alternative='greater')
alpha = 0.05
print("p value is " + str(p))
if p <= alpha:
print('The difference in mean return is significantly different (reject null hypothesis)')
else:
print('The difference in mean return is not significantly different (failed to reject null hypothesis)')
def CI_test(self):
bs_replicates = np.empty(10000)
for i in range(10000):
bs_replicates[i] = np.mean(np.random.choice(self.ret_data1['portfolio_ret'], len(self.ret_data1['portfolio_ret'])))
# Bootstrapped confidence intervals
conf_int = np.percentile(bs_replicates,[2.5, 97.5])
print('95% bootstrapped confidence interval =', conf_int, '%')
# Theoretical confidence intervals
conf_int_actual_upper = self.ret_data1['portfolio_ret'].mean() + self.ret_data1['portfolio_ret'].std()/np.sqrt(self.ret_data1['portfolio_ret'].count())*1.96
conf_int_actual_lower = self.ret_data1['portfolio_ret'].mean() - self.ret_data1['portfolio_ret'].std()/np.sqrt(self.ret_data1['portfolio_ret'].count())*1.96
conf_int_actual = [conf_int_actual_lower, conf_int_actual_upper]
print('-'*120)
print('95% theoretical confidence interval =', conf_int_actual, '%')
def CI_test_benchmark(self):
bs_replicates = np.empty(10000)
for i in range(10000):
bs_replicates[i] = np.mean(np.random.choice(self.benchmark2.pct_change()[1:], len(self.benchmark2.pct_change()[1:])))
# Bootstrapped confidence intervals
conf_int = np.percentile(bs_replicates,[2.5, 97.5])
print('95% bootstrapped confidence interval =', conf_int, '%')
# Theoretical confidence intervals
conf_int_actual_upper = self.benchmark2.pct_change()[1:].mean() + self.benchmark2.pct_change()[1:].std()/np.sqrt(self.benchmark2.pct_change()[1:].count())*1.96
conf_int_actual_lower = self.benchmark2.pct_change()[1:].mean() - self.benchmark2.pct_change()[1:].std()/np.sqrt(self.benchmark2.pct_change()[1:].count())*1.96
conf_int_actual = [conf_int_actual_lower, conf_int_actual_upper]
print('-'*120)
print('95% theoretical confidence interval =', conf_int_actual, '%')
# compute efficient frontier & portfolio statistics with benchmarks
class EF_benchmark3:
def __init__(self, asset, asset_list,benchmark1,benchmark2,benchmark3, benchmark_name = None):
self.asset = asset
self.asset_list = asset_list
self.benchmark1 = benchmark1
self.benchmark2 = benchmark2
self.benchmark3 = benchmark3
self.benchmark_name = benchmark_name
self.computation(asset, asset_list,benchmark1,benchmark2,benchmark3)
def computation(self, asset, asset_list,benchmark1,benchmark2,benchmark3):
# Check NaN values in the data
nullin_df = pd.DataFrame(self.asset,columns=self.asset_list)
# Calculate portfolio mean return
mu = expected_returns.mean_historical_return(self.asset)
# Calculate portfolio return variance
sigma = risk_models.sample_cov(self.asset)
# weight bounds in negative allows shorting of stocks
ef = EfficientFrontier(mu, sigma, weight_bounds=(0,1))
ef.add_objective(objective_functions.L2_reg, gamma=2)
# optional constraints possible, read pypfopt documentation.
sharpe_portfolio=ef.max_sharpe(risk_free_rate=0.00078) #US Bond
sharpe_portfolio_wt=ef.clean_weights()
latest_prices = get_latest_prices(self.asset)
# Allocate Portfolio Value in $ as required to show number of shares/stocks to buy,
# also bounds for shorting will affect allocation
# Maximum Sharpe Portfolio Allocation $1000000
da = DiscreteAllocation(sharpe_portfolio_wt, latest_prices, total_portfolio_value=1000000)
allocation, leftover = da.greedy_portfolio()
max_sharpe_cla = cla.CLA(mu, sigma)
max_sharpe_cla.max_sharpe()
sharpe_portfolio_wt_list = list(sharpe_portfolio_wt.values())
ret_data1 = self.asset.pct_change()[1:]
weighted_returns1 = (sharpe_portfolio_wt_list * ret_data1)
portfolio_ret1 = pd.DataFrame(weighted_returns1.sum(axis=1))
ret_data1 = ret_data1.merge(portfolio_ret1, on="Date", how="left")
ret_data1 = ret_data1.rename(columns={0: "portfolio_ret"})
ret_data1['cumulative_portfolio_ret'] = (ret_data1['portfolio_ret'] + 1).cumprod()
ret_benchmark1 = self.benchmark1.pct_change()[1:]
ret_benchmark1['cumulative_benchmark_ret'] = (ret_benchmark1 + 1).cumprod()
ret_benchmark2 = self.benchmark2.pct_change()[1:]
ret_benchmark2['cumulative_benchmark_ret'] = (ret_benchmark2 + 1).cumprod()
ret_benchmark3 = self.benchmark3.pct_change()[1:]
ret_benchmark3['cumulative_benchmark_ret'] = (ret_benchmark3 + 1).cumprod()
self.sharpe_portfolio_wt = sharpe_portfolio_wt
self.portfolio_ret1 = portfolio_ret1
self.allocation = allocation
self.max_sharpe_cla = max_sharpe_cla
self.ret_data1 = ret_data1
self.ret_benchmark1 = ret_benchmark1
self.ret_benchmark2 = ret_benchmark2
self.ret_benchmark3 = ret_benchmark3
self.leftover = leftover
def SP_weight(self):
print(self.sharpe_portfolio_wt)
print(self.allocation)
print("Leftover Fund value for the maximum Sharpe portfolio is ${:.2f}".format(self.leftover))
P_weight = Plotting.plot_weights(self.sharpe_portfolio_wt)
return P_weight
def Efficient_frontier(self):
EF_plot = Plotting.plot_efficient_frontier(self.max_sharpe_cla, show_assets="True")
return EF_plot
def cumulative_return(self):
Cum_ret = sns.scatterplot('Date', 'cumulative_portfolio_ret', data=self.ret_data1, color = 'yellowgreen', size = 1)
Cum_ret = sns.scatterplot('Date', 'cumulative_benchmark_ret', data=self.ret_benchmark1, color = 'teal',size =1)
Cum_ret = sns.scatterplot('Date', 'cumulative_benchmark_ret', data=self.ret_benchmark2, color = 'olivedrab',size =1)
Cum_ret = sns.scatterplot('Date', 'cumulative_benchmark_ret', data=self.ret_benchmark3, color = 'darkgreen',size =1)
if self.benchmark_name != None:
plt.legend(labels = ["ESG", *self.benchmark_name])
else:
plt.legend(labels = ["ESG", "Benchmark"])
print(self.ret_data1['cumulative_portfolio_ret'][-1])
return Cum_ret
def plot_hist(self):
hist_ret = plt.hist(self.ret_data1['portfolio_ret'], bins=10)
print(st.norm.interval(alpha=0.95, loc=np.mean(self.ret_data1['portfolio_ret']), scale=st.sem(self.ret_data1['portfolio_ret'])))
return hist_ret
def port_stat(self):
geometric_port_return = np.prod(self.portfolio_ret1 + 1) ** (252/self.portfolio_ret1.shape[0]) - 1
annual_std = np.std(self.portfolio_ret1) * np.sqrt(252)
port_sharpe_ratio = geometric_port_return / annual_std
print("Annual Return", geometric_port_return)
print("Stdev", annual_std)
print("Sharpe Ratio", port_sharpe_ratio)
def hypothesis_testing(self):
#import scipy.stats as stats
plt1 = self.ret_data1['portfolio_ret'].hist(bins=100, color='r', alpha=0.5, figsize = (16,8))
plt1 = self.benchmark1.pct_change()[1:].hist(bins=100, color='olivedrab', alpha=0.5)
print("Difference in mean return (" + self.benchmark_name[0] + ", ESG): ")
print((self.ret_data1['portfolio_ret'].mean() - self.benchmark1.pct_change()[1:].mean())*100)
stat, p = st.ttest_ind(self.ret_data1['portfolio_ret'], self.benchmark1.pct_change()[1:], equal_var=False)#, alternative='greater')
alpha = 0.05
print("p value is " + str(p))
if p <= alpha:
print('The difference in mean return is significantly different (reject null hypothesis)')
else:
print('The difference in mean return is not significantly different (failed to reject null hypothesis)')
def CI_test(self):
bs_replicates = np.empty(10000)
for i in range(10000):
bs_replicates[i] = np.mean(np.random.choice(self.ret_data1['portfolio_ret'], len(self.ret_data1['portfolio_ret'])))
# Bootstrapped confidence intervals
conf_int = np.percentile(bs_replicates,[2.5, 97.5])
print('95% bootstrapped confidence interval =', conf_int, '%')
# Theoretical confidence intervals
conf_int_actual_upper = self.ret_data1['portfolio_ret'].mean() + self.ret_data1['portfolio_ret'].std()/np.sqrt(self.ret_data1['portfolio_ret'].count())*1.96
conf_int_actual_lower = self.ret_data1['portfolio_ret'].mean() - self.ret_data1['portfolio_ret'].std()/np.sqrt(self.ret_data1['portfolio_ret'].count())*1.96
conf_int_actual = [conf_int_actual_lower, conf_int_actual_upper]
print('-'*120)
print('95% theoretical confidence interval =', conf_int_actual, '%')
def CI_test_benchmark(self):
bs_replicates = np.empty(10000)
for i in range(10000):
bs_replicates[i] = np.mean(np.random.choice(self.benchmark1.pct_change()[1:], len(self.benchmark1.pct_change()[1:])))
# Bootstrapped confidence intervals
conf_int = np.percentile(bs_replicates,[2.5, 97.5])
print('95% bootstrapped confidence interval =', conf_int, '%')
# Theoretical confidence intervals
conf_int_actual_upper = self.benchmark1.pct_change()[1:].mean() + self.benchmark1.pct_change()[1:].std()/np.sqrt(self.benchmark1.pct_change()[1:].count())*1.96
conf_int_actual_lower = self.benchmark1.pct_change()[1:].mean() - self.benchmark1.pct_change()[1:].std()/np.sqrt(self.benchmark1.pct_change()[1:].count())*1.96
conf_int_actual = [conf_int_actual_lower, conf_int_actual_upper]
print('-'*120)
print('95% theoretical confidence interval =', conf_int_actual, '%')
#Ranked the data into 3 segment
def gmb(df, subset = "Energy"):
test = df[df["Industry"] == subset]
test = test['ESG_C_S']
test = test.groupby(test.index).mean()
category = pd.qcut(test, 3,labels=["bad", "medium", "good"])
category = pd.DataFrame(category)
bad = list(category[category['ESG_C_S'] == "bad"].index)
medium = list(category[category['ESG_C_S'] == "medium"].index)
good = list(category[category['ESG_C_S'] == "good"].index)
return bad, medium, good
E_bad, E_medium, E_good = gmb(df_I,"Energy")
M1_bad, M1_medium, M1_good = gmb(df_I,"Manufacturing")
M2_bad, M2_medium, M2_good = gmb(df_I,"Medical")
L_bad, L_medium, L_good = gmb(df_I,"Logistic")
bad = []
bad += E_bad
bad += M1_bad
bad += M2_bad
bad += L_bad
medium = []
medium += E_medium
medium += M1_medium
medium += M2_medium
medium += L_medium
good = []
good += E_good
good += M1_good
good += M2_good
good += L_good
sns.violinplot(x="Industry", y="ESG_C_S", data=df_I)# hue = 'Industry',
g_ef = EF(df_O[good], good)
b_ef = EF(df_O[bad], bad)
m_ef = EF(df_O[medium], medium)
g_ef.cumulative_return()
# good ESG stocks in each industry
E_good,M1_good,M2_good,L_good
pip install yahoo_fin
pip install requests_html
# Extract company data
import yahoo_fin.stock_info as si
test = si.get_stats_valuation('AEP')
test['Unnamed: 0']
items = ['Market Cap (intraday) 5', 'Enterprise Value 3', 'Trailing P/E', 'Forward P/E 1', 'PEG Ratio (5 yr expected) 1', 'Price/Sales (ttm)', 'Price/Book (mrq)', 'Enterprise Value/Revenue 3', 'Enterprise Value/EBITDA 7']
def stats1(asset_list, item):
dict_temp = {}
for i in asset_list:
#print(i)
temp = si.get_stats_valuation(i)
temp = temp.iloc[:,:2]
temp.columns = ["Attribute", "Recent"]
dict_temp[i] = temp[(temp["Attribute"] == item)]
combined_stats = pd.concat(dict_temp)
combined_stats = combined_stats.reset_index()
if item == 'Market Cap (intraday) 5' or item == 'Enterprise Value 3':
combined_stats.Recent = combined_stats.Recent.apply(lambda x: float(x[:-1]) if 'B' in x else x)
else:
combined_stats.Recent = combined_stats.Recent.apply(lambda x: float(x))
combined_stats = combined_stats[["level_0", "Recent"]]
combined_stats.columns = ['Firm',item]
combined_stats = combined_stats.sort_values(by=[item],ascending=False)
return combined_stats
def stats2(asset_list, *item):
output = pd.DataFrame(columns=['Firm'])
for i in item:
#print(i)
temp = stats1(asset_list, i)
output = output.merge(temp, how = 'outer', on = "Firm")
return output
#Select stock based on their finacial ratios
output_E = stats2(E_good, items[0],items[3],items[5])
output_M1 = stats2(M1_good, items[0],items[3],items[5])
output_M2 = stats2(M2_good, items[0],items[3],items[5])
output_L = stats2(L_good, items[0],items[3],items[5])
output_E
output_M1
output_M2
output_L
# Select stocks based on the above criteria
Selected_Stock = ["NEE", "LIN", "ENLAY", "ES", "SRE",'UNH', 'TMO', 'WST','RMD','LLY', "TSM", "TXN", "ECL", "AMAT", "KMB", "UNP", "FDX", "ASR", "EIX", "ZNH"]
def alpha_beta_calculation(asset, assetlist, benchmark, start_date = "2014-12-31", end_date = "2019-12-31"):
nullin_df = pd.DataFrame(asset.loc[start_date:end_date],columns=assetlist)
# Calculate portfolio mean return
mu = expected_returns.mean_historical_return(asset.loc[start_date:end_date])
# Calculate portfolio return variance
sigma = risk_models.sample_cov(asset.loc[start_date:end_date])
# weight bounds in negative allows shorting of stocks
ef = EfficientFrontier(mu, sigma, weight_bounds=(0,1))
ef.add_objective(objective_functions.L2_reg, gamma=2)
# optional constraints possible, read pypfopt documentation.
sharpe_portfolio=ef.max_sharpe(risk_free_rate=0.00078) #US Bond
sharpe_portfolio_wt=ef.clean_weights()
latest_prices = get_latest_prices(asset.loc[start_date:end_date])
# Allocate Portfolio Value in $ as required to show number of shares/stocks to buy,
# also bounds for shorting will affect allocation
# Maximum Sharpe Portfolio Allocation $1000000
da = DiscreteAllocation(sharpe_portfolio_wt, latest_prices, total_portfolio_value=1000000)
allocation, leftover = da.greedy_portfolio()
max_sharpe_cla = cla.CLA(mu, sigma)
max_sharpe_cla.max_sharpe()
sharpe_portfolio_wt_list = list(sharpe_portfolio_wt.values())
ret_data1 = asset.loc[start_date:end_date].pct_change()[1:]
weighted_returns1 = (sharpe_portfolio_wt_list * ret_data1)
portfolio_ret1 =weighted_returns1.sum(axis=1)
benchmark1 = benchmark.loc[start_date:end_date].pct_change()[1:]
(Beta, Alpha) = st.linregress(benchmark1.values, portfolio_ret1.values)[0:2]
year = int(end_date[:4]) - int(start_date[:4])
print(str(year) + "-year " + f"Beta: {Beta:.6f}", str(year) + "-year " +f"Alpha: {Alpha:.6f}")
#print(str(year) + "-year " + f"Beta: {res[0]:.6f}", str(year) + "-year " +f"Alpha: {res[1]:.6f}")
#idosycratic risk work in favor to the stock (show as residual)
# --> good ESG & valued stock
alpha_beta_calculation(df_O[Selected_Stock], Selected_Stock,benchmark_SPY)
alpha_beta_calculation(df_O[Selected_Stock], Selected_Stock,benchmark_SPY,"2018-12-31","2019-12-31")
alpha_beta_calculation(df_O[Selected_Stock], Selected_Stock,benchmark_SPY,"2016-12-31","2019-12-31")
alpha_beta_calculation(df_O_all[Selected_Stock], Selected_Stock,benchmark_SPY_all,"2016-09-30","2021-09-30")
alpha_beta_calculation(df_O_all[Selected_Stock], Selected_Stock,benchmark_SPY_all,"2020-09-30","2021-09-30")
alpha_beta_calculation(df_O_all[Selected_Stock], Selected_Stock,benchmark_SPY_all,"2018-09-30","2021-09-30")
selected_EF_3 = EF_benchmark3(df_O[Selected_Stock], Selected_Stock,benchmark_SPY, benchmark_IXIC,benchmark_DJI, _3benchmark)
M_EF = selected_EF_3.cumulative_return()
selected_EF_3.Efficient_frontier()
selected_EF_3.SP_weight()
g_ef.Efficient_frontier()
selected_EF_3.port_stat()
selected_EF_3.CI_test()
selected_EF_3.CI_test_benchmark()
selected_EF_5 = EF_benchmark5(df_O[Selected_Stock], Selected_Stock,benchmark_ESG3, benchmark_ESG4,benchmark_ESG5,benchmark_ESG6,benchmark_ESG7, _5benchmark)
S_EF = selected_EF_5.cumulative_return()
selected_EF_3.hypothesis_testing()
selected_EF_5.hypothesis_testing()
selected_EF_5.port_stat()
a = {'NEE': 0.0652, 'LIN': 0.03926, 'ENLAY': 0.05335, 'ES': 0.04204, 'SRE': 0.03119, 'UNH': 0.07758, 'TMO': 0.06344, 'WST': 0.07308, 'RMD': 0.07704, 'LLY': 0.05329, 'TSM': 0.07588, 'TXN': 0.06643, 'ECL': 0.04255, 'AMAT': 0.06896, 'KMB': 0.02258, 'UNP': 0.0358, 'FDX': 0.0, 'ASR': 0.03637, 'EIX': 0.02486, 'ZNH': 0.05111}
a = pd.DataFrame.from_dict(a,orient='index')
a = a.rename(columns={0: 'Weight'})
a = a.sort_values('Weight', ascending = False)
a.head(10)
b = df_I.loc[Selected_Stock]
test = b[['ESG_C_S','Env_S','Soc_S','Gov_S','ESG_Controversial','Industry']]
test1 = round(test.groupby('Industry').mean(),2)
test1
test1_2 = test1.mean()
test1_2 # portfolio
c = df_I
test = c[['ESG_C_S','Env_S','Soc_S','Gov_S','ESG_Controversial','Industry']]
test2 = round(test.groupby('Industry').mean(),2)
test2
test2_2 = test.mean()
test2_2 # overall (98) companies
features = ['Combined \n Score','Environment','Social','Government','Controversial','Combined \n score']
def plot_rader(value1, value2, features, value1_name = 'Portfolio',value2_name = 'Overall', rader_size = 10):
# 使用ggplot的绘图风格
plt.style.use('ggplot')
# 构造数据
values = value1.values.tolist()
values2 = value2.values.tolist()
feature = features
N = len(values)
# 设置雷达图的角度，用于平分切开一个圆面
angles=np.linspace(0, 2*np.pi, N, endpoint=False)
# 为了使雷达图一圈封闭起来，需要下面的步骤
values=np.concatenate((values,[values[0]]))
values2=np.concatenate((values2,[values2[0]]))
angles=np.concatenate((angles,[angles[0]]))
# 绘图
fig=plt.figure()
ax = fig.add_subplot(111, polar=True)
# 绘制折线图
ax.plot(angles, values, 'o-', linewidth=2, label = value1_name, color = 'limegreen')
# 填充颜色
ax.fill(angles, values, alpha=0.25,color= 'limegreen')
# 绘制第二条折线图
ax.plot(angles, values2, 'o-', linewidth=2, label = value2_name, color = 'darkgreen')
ax.fill(angles, values2, alpha=0.25, color = 'darkgreen')
# 添加每个特征的标签
ax.set_thetagrids(angles * 180/np.pi, feature)
# 设置雷达图的范围
ax.set_ylim(0,rader_size)
# 添加标题
plt.title('ESG_Rader')
# 添加网格线
ax.grid(True)
# 设置图例
plt.legend(loc = 2, prop={'size': 7.4})
# 显示图形
plt.show()
plot_rader(test1_2, test2_2, features, 'Portfolio', 'Overall',10)
plot_rader(test1.loc['Energy'], test2.loc['Energy'], features, 'Portfolio', 'Overall',12)
plot_rader(test1.loc['Medical'], test2.loc['Medical'], features, 'Portfolio', 'Overall',12)
plot_rader(test1.loc['Manufacturing'], test2.loc['Manufacturing'], features, 'Portfolio', 'Overall',12)
plot_rader(test1.loc['Logistic'], test2.loc['Logistic'], features, 'Portfolio', 'Overall',12)
# Proportional to Industrial
I = pd.DataFrame(df_I.loc[Selected_Stock][['Year','Industry']])
I_1 = I[I['Year']=='2014-12-31 00:00:00']
I_2 = a.merge(I_1["Industry"], how = 'inner', left_on = a.index, right_on = I_1.index)
I_2 = I_2.groupby("Industry").sum()
# Pie chart, where the slices will be ordered and plotted counter-clockwise:
labels = I_2.index
sizes = round(I_2['Weight'],2)
fig1, ax1 = plt.subplots()
ax1.pie(sizes, labels=labels, autopct='%1.1f%%',
shadow=False, startangle=90, colors = ['darkseagreen','olivedrab','green', 'mediumseagreen'])
ax1.axis('equal') # Equal aspect ratio ensures that pie is drawn as a circle.
plt.title("Portfolio Weighting by Industry")
plt.show()