Sequence of return risk analysis for the S&P 500 Index from 1871 to 2020¶

# Using sp500 data from R. Schiller. Source: http://www.econ.yale.edu/~shiller/data.htm
df = pd.read_excel('ie_data.xls', sheet_name='Data', header=7)
df = df[['Date', 'P', 'D', 'CPI', 'Price', 'Dividend', 'Price.1']]\
.rename(columns={'Date':'Time', 'Price':'RP', 'Dividend':'RD', 'Price.1':'RTRP' })
'''
Time: year.month
P: Price, 
D: Dividend
CPI: Consumer price index, 
RP: real price (cpi adj. rel. to last date)
RD: real dividend (cpi adj. rel. to last date)
RTRP: Real total return price ((RP[n]+RD[n])/RP[n-1] )
'''
df = df.apply(lambda x: pd.to_numeric(x, errors = 'coerce')).dropna()
print(type(df.loc[df.index[-2], 'Time']))
df['Time'] = pd.to_datetime(df['Time'].astype(str), format='%Y.%m')
display("Head and tail of modified dataset")
display(df.head(), df.tail())

# plt.hist(df['RTR%'], bins=25, density=True, alpha=0.6, color='#1f77b4')
p = plt.hist(df['RTR'].dropna(), bins=100, density=True)

mu, std = stats.norm.fit(df['RTR'].dropna())
print(f'mu:{mu}',f'std:{std}')

# Plot the PDF.
xmin, xmax = plt.xlim()
x = np.linspace(xmin, xmax, 100)
pdf = stats.norm.pdf(x, mu, std)
plt.plot(x, pdf, 'k', linewidth=2)
title = f"Fit results: mu = {mu:.2f},  std = {std:.2f}"
display('Dataset Histogram - Monthly real total return')
plt.title(title)
plt.show()
display(df['RTR'].describe())
display(df['RTR%'].describe())

df = df_m.copy()
# Gouping by year
df.head()
df = df.groupby(df.Time.dt.year).mean().reset_index()
df['Time'] = pd.to_datetime(df['Time'].astype(str), format='%Y')
df['RTR'] = df['RTRP'].divide(df['RTRP'].shift(1))
df['RTR%'] = (df['RTR']-1)*100
df.dropna(inplace=True)
display('Full dataset grouped by year')
display(df)

## Plotting ten year returns
display('Plotting 10 year trailing returns for each year')
s = pd.Series(df['RTRP'].shift(-10).divide(df['RTRP']).values, index=df.Time)\
.subtract(1).multiply(100).dropna()
ax = s.plot(figsize=(12,8), grid=True)
# xtick = pd.date_range( start=ts.index.min( ), end=ts.index.max( ), freq='W' )
# ax.set_xticks( xtick, minor=True )
ytick = np.linspace(s.min()*0.8,s.max()*1.2,40)
_ = ax.set_yticks(ytick, minor=True)
# ax.grid('on', which='minor', axis='x' )
# ax.grid('off', which='major', axis='x' )

##RTR-percent norm distr
# plt.hist(df['RTR%'], bins=25, density=True, alpha=0.6, color='#1f77b4')
plt.hist(df['RTR%'].dropna(), bins=60, density=True)

mu, std = stats.norm.fit(df['RTR%'].dropna())
print(f'mu:{mu}',f'std:{std}')

# Plot the PDF.
xmin, xmax = plt.xlim()
x = np.linspace(xmin, xmax, 100)
pdf = stats.norm.pdf(x, mu, std)
plt.plot(x, pdf, 'k', linewidth=2)
title = f"Fit results: mu = {mu:.2f},  std = {std:.2f}"
display('Plotting histogram of real total returns percent for each year with fitted normal distribution')
plt.title(title)
plt.show()
df['RTR%'].describe()

##RTR norm distr
# plt.hist(df['RTR%'], bins=25, density=True, alpha=0.6, color='#1f77b4')
p_h = plt.hist(df['RTR'].dropna(), bins=60, density=True)
mu, std = stats.norm.fit(df['RTR'].dropna())
print(f'mu:{mu}',f'std:{std}')

# Plot the PDF.
xmin, xmax = plt.xlim()
x = np.linspace(xmin, xmax, 1000)
pdf = stats.norm.pdf(x, mu, std)
plt.plot(x, pdf, 'k', linewidth=2)
display('Plotting histogram of real total returns with fitted normal distribution')
title = f"Fit results: mu = {mu:.2f},  std = {std:.2f}"
plt.title(title)
plt.show()
display(df['RTR'].describe())

# Norm distr fit to log(x) values
# print(df['RTR'])
rtr_ln = np.log(df['RTR'])
# print(rtr_ln)
p = plt.hist(rtr_ln, bins=60, density=True)
mu, std = stats.norm.fit(rtr_ln)
print(f'mu:{mu}',f'std:{std}')
# # print(np.trapz(p_h[0], p_h[1]))
# # print(p_h)
# print(len(p_h[0]))
# print(len(p_h[1]))
# delta = p_h[1][1]-p_h[1][0]
# print(delta)
# print(f'sum:{np.trapz(p[0],dx=delta)}')

# Plot the PDF.
display('Histogram of ln(x) and fitted normal probability distribution (for simplicity)')
xmin, xmax = plt.xlim()
x = np.linspace(xmin, xmax, 400)
delta = x[1]-x[0]
# x = np.arange(xmin, xmax, delta)
pdf = stats.norm.pdf(x, mu, std)
plt.plot(x, pdf, 'k', linewidth=2)
title = f"Fit results: mu = {mu:.2f},  std = {std:.2f}"
plt.title(title)
plt.show()
# print(pdf)
len(pdf)
print(f'sum:{np.trapz(pdf,dx=delta)}')

## Setting up convolution calc params
# mu_check = []
# for res in [1000, 20000, 50000, 100000, 150000, 200000]:
# print(res)
# x = np.linspace(-40, 40, res)
x = np.linspace(-20, 20, 10000) # for the convolution to work, need to have to be symmetric, large x vector and cover practial range as n->
# need high resolution to avoid border effects during convolution
delta = x[1]-x[0]
pdf = stats.norm.pdf(x, mu, std)
pmf = pdf*delta
print(sum(pmf))
print(sum(pdf)*delta)
print(len(pmf))
# plt.plot(x, pmf, 'k', linewidth=2)
# title = f"Fit results for pmf"
# plt.title(title)
print()

n=5
pmf_conv = convolute_list_of_pmfs(n*[pmf])
print(sum(pmf_conv))
plt.plot(x, pmf_conv, 'k', linewidth=2)
plt.xlim([-2, 15])
title = f"Convolution of {n} pdf(ln(x))'s with mu = {mu:.4f},  std = {std:.4f}"
plt.title(title)
# mu_, std_ = mean_and_std_from_pmf(x, pmf_conv)
# print(mu_, std_)
# mu_check.append(mu_)
# print(mu_check)

mu_conv, std_conv = mean_and_std_from_pmf(x, pmf_conv)
# mu_test, std_test = create_pdf_x_params_from_norm_pdf_lnx_params(mu_conv, std_conv)
cdf_log_test = stats.lognorm.cdf(x, std_conv, scale=np.exp(mu_conv))
print(mu_conv, std_conv )
# print(mu_test, std_test )
pdf_log_test = stats.lognorm.pdf(x, std_conv, scale=np.exp(mu_conv))
pmf_log_test = pdf_log_test*delta
plt.plot(x, pdf_log_test, 'k', linewidth=2)
title = f"Fit results: mu = {mu:.2f},  std = {std:.2f}"
plt.title(title)

# stats.rv_continuous.cdf

cd = interp1d(cdf_log_test, x)(0.5)
print(cd)
# fig4, ax4 =  plt.subplots()
plt.plot(x, cdf_log_test, linewidth=2, label=f'n: {n}')
plt.xlim([-2, 10])

# Checking the sums
n = 15 # periods
d = create_all_pmfs_pdfs_from_periods_and_pmf_lnx(x, pmf, n)

## Checking sums etc
# print(d['pdfs_x_params'])
print([sum(x) for x in d['pmfs_x']])
print([sum(x) for x in d['pmf_lognorms']])
print(len(d['pdfs_x'][0]))
print(len(d['pmf_lognorms']))
print(d['pdfs_x_params'])
# cd = interp1d(d['cdf_lognorms'][4], x)(0.5)
# print(cd)

print(d['pdfs_x_params'])
print('Checking last value for cumulative distribution')
print([x[-1] for x in d['cdf_lognorms']])

## Plotting probability distributions
fig,ax = plt.subplots()
fig2, ax2 =  plt.subplots()
fig5, ax5 =  plt.subplots()
plot_nums = 14
for i in range(0,plot_nums-1,1):
    ax.plot(x, d['pmfs_x'][i], linewidth=2, label=f'n: {i}')
    ax.set_xlim([-2, 3])
    ax2.plot(x, d['pmf_lognorms'][i], linewidth=2, label=f'n: {i}')
    ax2.set_xlim([0, 8])
    ax5.plot(x, d['cdf_lognorms'][i], linewidth=2, label=f'n: {i}')
    ax5.set_xlim([0, 10])
#     ax2.set_yscale('log')

fig.suptitle("pmf from convolution of n number of pdf(ln(x)) distributions")
fig2.suptitle("Resulting probability mass function\n for the SP500 real total return over n periods")
fig.legend(loc='right')
fig2.legend(loc='right')
fig5.legend(loc='right')
fig5.suptitle("Cumulative probability distribution \n for the SP500 real total return over n periods")

# Single plots example
plot_num = 0
fig3, ax3 =  plt.subplots()
ax3.plot(x, d['pmf_lognorms'][plot_num], linewidth=2, label=f'n: {n}')
ax3.set_xlim([0, 20])
fig3.suptitle(f"Single plot of probability mass function for {plot_num} periods")

fig4, ax4 =  plt.subplots()
ax4.plot(x, d['cdf_lognorms'][plot_num], linewidth=2, label=f'n: {n}')
ax4.set_xlim([0, 2])
fig4.suptitle(f"Single plot of cumulative probability distribution for {plot_num} periods")

# Check cumulative distribution
cd = interp1d(d['cdf_lognorms'][9], x)
print(cd(0.5))
print(cd(0.05))

## DATAFRAME for norm(pdf(ln(x))) dataset
# params_log = [mean_and_std_from_pmf(x, d['pmf_lognorms'][n]) for i in range(1,n)]
df_res = pd.DataFrame(index = range(1,n+1))
df_res.index.name= 'year'
# display(df_res)
# print(n)
all_cd = [interp1d(d['cdf_lognorms'][i], x) for i in range(0,n)]
# df_res['ev_rtr'] = [np.exp(x[0]) for x in d['pdfs_x_params']]
df_res['ev_rtr'] = [x[0] for x in d['pdfs_y_params']]
df_res['ev_return%'] = df_res['ev_rtr'].subtract(1).multiply(100)
# s_ev = pd.Series([x[0] for x in d['pdfs_y_params']])
df_res['gmean_ev_rtr_pa'] = [x**(1/(i+1)) for i,x in enumerate(df_res['ev_rtr'].values)]
df_res['gmean_ev_rtr%_pa'] = df_res['gmean_ev_rtr_pa'].subtract(1).multiply(100)
df_res['ar_mean_rtr_pa'] = [((x-1)/(i+1))+1 for i,x in enumerate(df_res['ev_rtr'].values)]
# df_res['median_rtr'] = [interp1d(d['cdf_lognorms'][i], x)(0.5) for i in range(0,n)]
df_res['median_rtr'] = [all_cd[i](0.5) for i in range(0,n)]
df_res['median_rtr%'] = df_res['median_rtr'].subtract(1).multiply(100)
df_res['gmean_median_rtr'] = [x**(1/(i+1)) for i,x in enumerate(df_res['median_rtr'].values)]
df_res['5_conf_rtr'] = [all_cd[i](0.05) for i in range(0,n)]
df_res['95_conf_rtr'] = [all_cd[i](0.95) for i in range(0,n)]
df_res['5_conf_gain%'] = df_res['5_conf_rtr'].subtract(1).multiply(100)
df_res['95_conf_gain%'] = df_res['95_conf_rtr'].subtract(1).multiply(100)
df_res['5_conf_gmean_rtr_pa'] = [x**(1/(i+1)) for i,x in enumerate(df_res['5_conf_rtr'].values)]
df_res['95_conf_gmean_rtr_pa'] = [x**(1/(i+1)) for i,x in enumerate(df_res['95_conf_rtr'].values)]
df_res['5_conf_gmean%'] = df_res['5_conf_gmean_rtr_pa'].subtract(1).multiply(100)
df_res['95_conf_gmean%'] = df_res['95_conf_gmean_rtr_pa'].subtract(1).multiply(100)

#Finne x hvor cd(x)=1.04 (hvor sannsynlig det er å gå under 4% swr)?
display('Summary table of S&P 500 sequence of return risk analysis')
display('Including confidence intervals for real total returns and geometric means (annualized returns)')
display(df_res)
# display(pd.DataFrame(d['cdf_lognorms']))