import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
import scipy.stats as stats
# If the observations are in a dataframe, you can use statsmodels.formulas.api to do the regression instead
from statsmodels import regression
from statsmodels.stats.stattools import jarque_bera

# Set start and end dates for data
start = '2002-06-14'
end = '2017-06-14'

# Load S&P 500, 30-day Fed Funds Futures and 5-year Treasury Note pricing data 
SPY = get_pricing('SPY', fields='price', start_date=start, end_date=end)
FF = get_pricing('FF', fields='price', start_date=start, end_date=end)
FV = get_pricing('FV', fields='price', start_date=start, end_date=end)

# Import Effective Fed Funds Rate from FRED (URL at bottom of document)
data = local_csv('DFF.csv')

# Because the SPY, FF, FV series only include market days, we will use the ffill to
# fill in missing values for the weekends and holidays which fedfunds already has
calendar_dates = pd.date_range(start=start, end=end, freq='D', tz='UTC')
SPY = SPY.reindex(calendar_dates, method='ffill')
FF = FF.reindex(calendar_dates, method='ffill')
FV = FV.reindex(calendar_dates, method='ffill')

# Let's match the length and dates of the different data sets
fedfunds = data["DFF"][-(len(calendar_dates)):]
fedfundsdates = data["DATE"][-(len(calendar_dates)):]
if (len(SPY) == len(FF)) and (len(FF) == len(FV)) and (len(FV) == len(fedfunds)):
    print "The length of all variables is now " + str(len(SPY))
    print "SPY spans from " + str(SPY.index[0]) + " to " + str(SPY.index[-1])
    print "FedFunds spans from " + str(fedfundsdates.iloc[0]) + " to " + str(fedfundsdates.iloc[-1])

# Get the percentage change values (multiplacative returns)
spy_returns = SPY.pct_change()[1:]
ff_returns = FF.pct_change()[1:]
fv_returns = FV.pct_change()[1:]
fedfunds_returns = fedfunds.pct_change()[1:]
fedfunds = fedfunds[1:]
fedfunds.index = spy_returns.index
fedfunds_returns.index = spy_returns.index

# See what the fed funds rates look like over 15 years
plt.plot(fedfunds)
plt.legend(["Fed Funds Interest Rate"]);

# See what SPY ETF looks like over the past 15 years
plt.plot(SPY)
plt.legend(["S&P 500"]);

# Create a data frame with all the variables inside it indexed to the same timedate
df = pd.DataFrame({'SPYReturns' : spy_returns,
                   'SPY' : SPY,
#                   'FF' : ff_returns, 
#                   'FV' : fv_returns,
                   'FedFundsReturns' : fedfunds_returns,
                   'FedFunds' : fedfunds})
#Let's ignore FF and FV until we do the rest first
# Now we can index all variables along the same datetime or value or boolean index
print df[df['FedFundsReturns']>0].head(3)
print df[df['FedFundsReturns']>0].tail(3)
v1 = df[df['FedFundsReturns']<=0]['SPYReturns']
v2 = df[df['FedFundsReturns']>0]['SPYReturns']

plt.scatter(v1,df[df['FedFundsReturns']<=0]['FedFundsReturns']);
plt.legend(["SPYReturns and FedFundsReturns"])
plt.title(["FedFundsReturns <= 0"]);

plt.scatter(v2,df[df['FedFundsReturns']>0]['FedFundsReturns']);
plt.legend(["SPYReturns and FedFundsReturns"])
plt.title(["FedFundsReturns > 0"]);

def summary(x):
    print "Mean:   " + str(x.mean())
    print "SD:     " + str(x.std())
    print "Skew:   " + str(x.skew())
    print "Kurt:   " + str(x.kurtosis())
    print "Median: " + str(x.median())
    print x.describe()[3:];

plt.hist(v1)
plt.legend(["SPYReturns"])
plt.title(["FedFundsReturns <= 0"]);

plt.hist(v2)
plt.legend(["SPYReturns"])
plt.title(["FedFundsReturns > 0"]);

summary(v1)

summary(v2)

def normal_test(x):
    _, pvalue, _, _ = jarque_bera(x)
    if pvalue > 0.05:
        print 'The variable is likely normally distributed.'
    else:
        print 'The variable is likely not normally distributed.'
normal_test(v1)
normal_test(v2)

alpha = 0.05
n1 = len(v1)
sd1 = v1.std()
n2 = len(v2)
sd2 = v2.std()
print "Group 1, Interest-Rate Change <= 0:"
print "n1 = " + str(n1)
print "SD1 = " + str(sd1)
print "Group 2, Interest-Rate Change > 0:"
print "n2 = " + str(n2)
print "SD2 = " + str(sd2)

var1 = sd1**2
var2 = sd2**2
test_stat = var1/var2
print "Test-statistic, or variance ratio, is: " + str(test_stat)

f_crit = stats.f.ppf(1-alpha, (n1-1), (n2-1))
if test_stat > f_crit:
    print "We reject the null hypothesis because our test_stat of " + str(test_stat) + \
    " is greater than our f_crit of " + str(f_crit) + " at alpha " + str(alpha)
else:
    print "We accept the null hypothesis because our test_stat of " + str(test_stat) + \
    " is less than our f_crit " + str(f_crit) + " at alpha " + str(alpha)

# It is important to note that with the Levene test one does not want to run large
# sample sizes (or the population) because the test_stat gets impacted heavily with 
# the degrees of freedom it accounts for (the test_stat will rise greatly)
# Taking samples of 5% the v1 and v2 lengths, we run the levene test on their variances
test_stat, p_value = stats.levene(v1.sample(len(v1)/20), v2.sample(len(v2)/20), center='mean')
# N is the combined sum length of all samples
N = len(v1)*20
# K is essentially the degrees of freedom allowed for all the 'treatments' (in this case, interest rates 
# going up or staying constant/decreasing)
k = 2

# Find f_crit 
f_crit = stats.f.ppf(alpha,(k-1),(N-k))

if test_stat > f_crit:
    print "We reject the null hypothesis because our test_stat of " + str(test_stat) + \
    " is greater than our f_crit of " + str(f_crit) + " at alpha " + str(alpha) + \
    " with p-value " + str(p_value)
else:
    print "We accept the null hypothesis because our test_stat of " + str(test_stat) + \
    " is less than our f_crit " + str(f_crit) + " at alpha " + str(alpha) + \
    " with p-value " + str(p_value)