import numpy as np
import math
import scipy.stats as scs
import matplotlib.pyplot as plt
import pandas as pd

plt.style.use('seaborn')

""" I am going to looking at the statistical characterisitics of financial instruments in this notebook so I am going to
write a function for describing there key charachterisitics """

def print_statistics(array): 
    """ Returns min max mean and stdev of array
    args:
        
        array(np.array): the financial returns of the instrument
        
    Returns
    
        min(float), max(float), mean(float), std(float)
    """
    import scipy.stats as scs   
    stats=scs.describe(array)
    
    print('%14s %15s' % ('statistic','value'))
    print(30*'-')
    print('%14s %15s' % ('min',stats[1][0]))
    print('%14s %15s' % ('max',stats[1][1]))
    print('%14s %15s' % ('mean',stats[2]))
    print('%14s %15s' % ('standard dev',stats[3]))
    print('%14s %15s' % ('skew',stats[4]))
    print('%14s %15s' % ('kurtosis',stats[5]))
    
def import_returns(ticker):
    """ returns daily data of a stock I have storedd locally"""
    
    import pandas as pd
    data = pd.read_csv(r"C:\Users\Adam\Desktop\Personal\Python Syntax\Python Input Files\finance\{}.csv".format(ticker), parse_dates=True ,index_col=0)
    data.sort_index(ascending=True, inplace=True)
    data['Log_Return'] = np.log((data['Price'])/data['Price'].shift(1))
    data.dropna(inplace=True)
    return data 

"""    
First of all we are going we are going to test the theory that financial instruments can 
be modelled using normal distirbutions by actually looking at the returns of some instruments that
we will later put into a portfolio. The stocks are Apple stock, BP stock, Macdonalds, Stock and 
"""

SBUX = import_returns('SBUX')['Log_Return']
APPL = import_returns('APPL')['Log_Return']
MCD = import_returns('MCD')['Log_Return']
GLD = import_returns('GLD')['Log_Return']

All_data = pd.concat([SBUX,APPL,MCD,GLD],axis=1,join='inner')
All_data.columns=['SBUX','AAPL','MCD','GLD']

""" upon first inspection the graphs suggest that each if the instruments do seem to have a normal esque."""

All_data.hist(bins=50,figsize=(16,8))

""" we can plot a hypothetical normaldistibution using the normal and stdeviation of the returns to get a 
better look. The Gold returns distribuions is narrower than a normal distribution  with a longer tail it 
exhibits positive kurtosis also known as leptokurtic""" 

plt.figure(figsize=(16,8))
GLD.hist(bins=50,  label ='observed', density = True)
plt.xlabel('Distribution of GLD ETF daily returns')
plt.xlabel('GLD daily returns')
plt.ylabel('frequency')

x=np.linspace(plt.axis()[0],plt.axis()[1])
plt.plot(x, scs.norm.pdf(x, GLD.mean(),GLD.std()), color='r',lw=2.0, label = 'pdf')
plt.legend()

""" we can confirm this by printing the statistics - and as shown from the graph
the kurtosis is not close to zero showing the distribution does not follow that of a normal distbution.
It also exhibits some negative skew - the left tail is longer than the right"""

print_statistics(GLD)

""" 
We can test simply for normally using the scs functions I have built some standard nomrlaity tests into the custom functino below
"""

def normality_test(array):
    """ test the normality of a series or array
    args:
        
        array(np.array or pd.series) : data set to be tested for normality
    """
    print('Skew of data set: %8.3f' % scs.skew(array))
    print('Skew test p-value: %7.3f' % scs.skewtest(array)[1])
    print('Kurt of data set: %8.3f' % scs.kurtosis(array))
    print('Kurt test p value: %7.3f' % scs.kurtosistest(array)[1])
    print('normal test p value: %4.3f' % scs.normaltest(array)[1])


""" what you can see us that gold fails all 3 test showing that its returns are not accurately represented by a normal distribution """

print("Gold returns distribution characterisitcs and p values")
normality_test(GLD)
print('\n' + "Starbucks returns distribution characterisitcs and p values")
normality_test(SBUX)
print('\n' + "APPL returns distribution characterisitcs and p values")
normality_test(APPL)
print('\n' + "Mcdonalds returns distribution characterisitcs and p values")
normality_test(MCD)

def import_annual_returns(ticker):
    
    import pandas as pd
    data = pd.read_csv(r"C:\Users\Adam\Desktop\Personal\Python Syntax\Python Input Files\finance\{}.csv".format(ticker), parse_dates=True ,index_col=0)
    data.sort_index(ascending=True, inplace=True)
    data = data.resample('1Y').last()
    data['Annual_Log_Return'] = np.log((data['Price'])/data['Price'].shift(1))
    data.dropna(inplace=True)
    return data 


SBUX = import_annual_returns('SBUX')['Annual_Log_Return']
APPL = import_annual_returns('APPL')['Annual_Log_Return']
MCD = import_annual_returns('MCD')['Annual_Log_Return']
GLD = import_annual_returns('GLD')['Annual_Log_Return']

annual_returns = pd.concat([SBUX,APPL,MCD,GLD],axis=1,join='inner')
annual_returns.columns=['SBUX','APPL','MCD','GLD']

""" In portfolio thory all position must be long so our first basic portolfio will just be a 25% split across all 4 instruments """


portfolio_weights = np.array([0.25,0.25,0.25,0.25])

mean_annual_returns = np.array(annual_returns.mean())

Expected_Portfolio_Return = np.dot(portfolio_weights.T,mean_annual_returns)

#the .cov methdo returns the covariance between each instrument in a 2d array 
Expected_Portfolio_Volatility = math.sqrt(np.dot(portfolio_weights.T,np.dot(annual_returns.cov(),portfolio_weights.T)))

import math

print( "The expected portfolio return from a portfolio split 25% each \
is: {0: 2.2f}% and has volatility:{1: 2.2f}%".format(Expected_Portfolio_Return*100,Expected_Portfolio_Volatility*100))
      
      
      

def portfolio_return(portfolio_weights):
     return np.dot(portfolio_weights.T,np.array(annual_returns.mean()))
     
def portfolio_volatility(portfolio_weights):
     return math.sqrt(np.dot(portfolio_weights.T,np.dot(annual_returns.cov(),portfolio_weights.T)))

portfolio_returns = []
portfolio_volatilities= [] 

import numpy.random as random 

for portfolio_mix in range(10000):
    
    weights = random.rand(len(annual_returns.columns))
    weights = weights / weights.sum()
    
    portfolio_returns.append(portfolio_return(weights))
    portfolio_volatilities.append(portfolio_volatility(weights))

portfolio_returns = np.array(portfolio_returns)
portfolio_volatilities = np.array(portfolio_volatilities)   
    
"""
The sharpe ratio gives us an indication of risk adjusted returns of a portfolio, the higher the sharpe ratio
the greater the risk adjusted returns 

sharpe ratio = (portfolio return - riskfree rate)/portfolio volatility

the riskfree at the moment in the UK = 0.25% """

Rf = 0.0025

Sharpe_Ratios = (portfolio_returns-Rf)/portfolio_volatilities

plt.figure(figsize = (18,6))
plt.scatter(portfolio_volatilities,portfolio_returns, c = Sharpe_Ratios, marker ='o' , cmap ='coolwarm')
plt.title('Expected portfolio returns vs volatilities for different combination weights of SBUX/BP/MCD/GLD')
plt.ylabel('Expected portfolio return')
plt.xlabel('Expected portfolio volatility')
plt.colorbar(label = 'Sharpe Ratio')

import scipy.optimize as sco 

#number 1) the output to minimize 

def minimum_sharpe(weights):
    return - (portfolio_return(weights) - Rf)/portfolio_volatility(portfolio_weights)

#number 2) the parameters to iterate throught - inputted as an intial values array

weights = np.array([0.25,0.25,0.25,0.25])

#number 3) the minimizsing algo 

minimizing_algo = 'SLSQP'

#4) establishing paramter boundaries, the weight for each instrument must reside between 0 and 1 entered as a tuple of tuples

boundaries = tuple((0.0,1) for x in range(len(annual_returns.columns)))

#step 5 ) optimizing constraint - also entered as a tuple - input paramters must sum to 1

constraints = ({'type' : 'eq', 'fun': lambda x: np.sum(x)-1})


optimized_params_sharpe = sco.minimize(minimum_sharpe, weights, method = minimizing_algo , bounds = boundaries , constraints = constraints) 

print(optimized_params_sharpe)

#access the weights attribute of the optimization and round to 3 dp
sharpe_optimized_weights = optimized_params_sharpe['x'].round(3)

""" enter this back into our portfolio return and volatility function to establish the return and vol"""

print(' The optimal portfolio has: {}% SBUX {}% APPL {}% MCD {}% GLD'.format(sharpe_optimized_weights[0]*100,sharpe_optimized_weights[1]*100,sharpe_optimized_weights[2]*100,sharpe_optimized_weights[3]*100))

max_sharpe_return = portfolio_return(sharpe_optimized_weights)
max_sharpe_volatility = portfolio_volatility(sharpe_optimized_weights)
max_sharpe_ratio = (max_sharpe_return - Rf)/max_sharpe_volatility

print(' Its Sharpe Ratio is:{}'.format(max_sharpe_ratio))

EF_constraints = ({'type' : 'eq', 'fun': lambda x: portfolio_return(x)-target_return},
                  {'type' : 'eq', 'fun': lambda x: np.sum(x)-1})

target_returns = np.linspace(0.02, 0.175, 50)

min_vols = [] 

# this time we are minimizing the portfolio volatility output for a constrained portfolio return
for target_return in target_returns:
    optimized_params = sco.minimize(portfolio_volatility, weights, method = minimizing_algo , bounds = boundaries , constraints = EF_constraints)
    # the fun attribute of the optimized params variable is the minimized output value in this case the volatility
    min_vols.append(optimized_params['fun'])
    
    
min_vols = np.array(min_vols)

plt.figure(figsize =(18,6))
plt.scatter(portfolio_volatilities,portfolio_returns, c = Sharpe_Ratios, marker ='o' , cmap ='coolwarm')
plt.plot(min_vols,target_returns, 'b', lw=4.0)
plt.plot(max_sharpe_volatility,max_sharpe_return, 'y*', markersize=15.0)
plt.ylabel('expected portfolio return')
plt.xlabel('expected portfolio volatility')
plt.colorbar(label = 'Sharpe Ratio')

""" This is another 'solver' situation so we can use the scipy package! """

def equations(p, rf=Rf): 
    eq1 = p[1] - rf
    eq2 = rf + p[0]*p[2] - f(p[2])
    eq3 = p[0] - df(p[2])
    return eq1, eq2, eq3
    
"""to find the efficient frontier f(x) and f'(x) we can use cubic spline interplotation which is part of the scipy package it returns the function of 
the line of best fit to data, we will feed it the data point of the efficient frontier to get back a fucntion that
decribes it"""

import scipy.interpolate as sci

""" we use the np.argmin function to remove repeat vol indices by ideniying the index of the first repeat value so that we 
can only keep volatility points from there on after (the efficient frontier) we are interested in the tangiential crossing with the top of the curve. maximum returns """


repeated_index = np.argmin(min_vols)

# nreturn only the coords for the top half of the curve the efficient frontier also need to sort the data for interpolation
# ive added in -5 as the last few data points were duplicate values  and the interpolate algorithmn cant handle them
EFrontier_vols = min_vols[repeated_index:-5]
EFrontier_returns = target_returns[repeated_index:-5]

#fitting the cubic spline to the frontier data points
CubicSplineInterpolation = sci.splrep(EFrontier_vols,EFrontier_returns)

""" define a function for the efficent frontier function and its 1st derivation"""

def f(x):
    ''' Efficient frontier function (the apporximation from the splines interpolation)'''
    return sci.splev(x,CubicSplineInterpolation,der=0)

def df(x):
    ''' Efficient frontier function  1st derivative'''
    return sci.splev(x,CubicSplineInterpolation,der=1)


opt = sco.fsolve(equations, [0.5,0.01,0.15])

tangent_x = opt[2]
tangent_y = opt[0]*opt[2]+opt[1]

print("The equation of the capital market line is y = {0: 2.4f}x + {1: 2.4f} and the tangent it at ({2: 2.4f},{3: 2.4f})".format(opt[0],opt[1],tangent_x,tangent_y)+ '\n')


""" plotting the capitla market line  """ 

plt.figure(figsize=(18,6))
plt.scatter(portfolio_volatilities, portfolio_returns, c = Sharpe_Ratios , marker ='.', cmap ='coolwarm' )
plt.colorbar(label = 'Sharpe Ratio')
plt.scatter(EFrontier_vols, EFrontier_returns, color='b',lw=4.0)

cx = np.linspace(0,0.3)
plt.plot(cx,opt[0]*cx + opt[1], markersize = 15.0)
plt.grid(True)
plt.ylabel('expected portfolio return')
plt.xlabel('expected portfolio volatility')
plt.xlim(0,0.2)
plt.ylim(0,0.25)

print('The optimal portfolio has: {}% SBUX {}% APPL {}% MCD {}% GLD'.format(sharpe_optimized_weights[0]*100,sharpe_optimized_weights[1]*100,sharpe_optimized_weights[2]*100,sharpe_optimized_weights[3]*100))

print('\n' +'The portfolio is predominantly APPL as apple for the last ten years has very high risk adjusted returns compared to the other 3 instruments show below: ' + '\n')

print_statistics(APPL)

# Utilising scipy.solve another theoretical example 

""" just to drill the above in we are going to find the crossing point of two lines using the scipy.solve() function """

""" 
this time we want to find the the crossing point of these 2 lines, to do this we need to solve the python function p(a,x) this time we have 2 equations that are just described by
2 unknowns so we need to write a function we 2 simulatanoues equations:

1) y = 2x + a

2) y = 20 - 3x
    
"""

import scipy.optimize as sco
    
def equations(p): 
    eq1 = 2*p[1]+ p[0]
    eq2 = 20 - 3*p[1]
        
    return eq1, eq2

opt = sco.fsolve(equations, [0,0])

print ("paramter a = : {0: 2.2f}  and paramter b :{1: 2.2f}".format(opt[0],opt[1]))