%matplotlib inline
import pandas as pd
import statsmodels.formula.api as smf
import statsmodels.api as sm
import matplotlib.pyplot as plt
import numpy as np

"""
ESTIMATING THE COEFFICIENTS
"""
# Dataset: How many people rent capitol bikeshare bikes
bike_dat = pd.read_csv("Data/day.csv")
bike_dat.head(2)

# Plot the data in a scatter plot
plt.scatter(bike_dat.atemp, bike_dat.cnt, alpha=0.3)  # Plot the raw data

# Estimate the model parameters
est_s = smf.ols(formula='cnt ~ temp', data=bike_dat).fit()

# View the model estimates
est_s.summary()

# Plot the data with scatter plot
plt.scatter(bike_dat.temp, bike_dat.cnt, alpha=0.3) 
plt.xlabel("Temperature")
plt.ylabel("Number of Bike Rentals")

# Generate coefficient estimates 
plt.plot([bike_dat.temp.min(), bike_dat.temp.max()],
        [est_s.params.Intercept + est_s.params.temp * bike_dat.temp.min(),
           est_s.params.Intercept + est_s.params.temp * bike_dat.temp.max()],
            linewidth=2)      

# Generate data points
x_prime = pd.DataFrame({'temp' : np.linspace(bike_dat.temp.min(), 
                                           bike_dat.temp.max(), 100)})

# Generate the predictions using the built in method
y_hat = est_s.predict(x_prime)

# Plot the data
plt.plot(x_prime, y_hat, 'r', linewidth=2, alpha=0.9)


# Now let's run a multiple linear regression
# The temp variable is no longer significant. Why? Multicollinearity
est_m = smf.ols(formula='cnt ~ atemp + temp + workingday + windspeed', 
                data=bike_dat).fit()
est_m.summary()



# Scatter plot (observe the (unsurprising) correlation between atemp and temp)
cols = ['cnt','atemp','windspeed','weathersit','temp','workingday','hum']
pd.scatter_matrix(bike_dat[cols])

# Correlation coefficient matrix
corr_matrix = np.corrcoef(bike_dat[cols].T)
sm.graphics.plot_corr(corr_matrix, xnames=cols)


# Let's say we wanted to include an interaction term
# We would do this by including the ':' between interacting variables
est_m = smf.ols(formula='cnt ~ temp + windspeed + temp:windspeed + workingday', 
                data=bike_dat).fit()

est_m.summary()

# An alternate way of specifying interaction terms
# a*b is equivalent to a + b + a:b
est_m = smf.ols(formula='cnt ~ temp*windspeed + workingday',data=bike_dat).fit()

est_m.summary()

"""
Explanation: Bikesharing data
Data Description
* datetime - hourly date + timestamp
* season -  1 = spring, 2 = summer, 3 = fall, 4 = winter 
* holiday - whether the day is considered a holiday
* workingday - whether the day is neither a weekend nor holiday
* weather - 1: Clear, Few clouds, Partly cloudy, Partly cloudy 
-  2: Mist + Cloudy, Mist + Broken clouds, Mist + Few clouds, Mist 
- 3: Light Snow, Light Rain + Thunderstorm + Scattered clouds, Light Rain + Scattered clouds 
- 4: Heavy Rain + Ice Pallets + Thunderstorm + Mist, Snow + Fog 
* temp - temperature in Celsius
* atemp - "feels like" temperature in Celsius
* humidity - relative humidity
* windspeed - wind speed
* casual - number of non-registered user rentals initiated
* registered - number of registered user rentals initiated
* count - number of total rentals
End of explanation
"""


# Dataset: How many people rent capitol bikeshare bikes
bike_hour_dat = pd.read_csv("Data/hour.csv")
bike_hour_dat.head(2)

# Multiple regression 
est_hour_m = smf.ols(formula='cnt ~ temp + hum + workingday + hr + weathersit', data=bike_hour_dat).fit()
est_hour_m.summary()

est_hour_m.rsquared_adj

est_m.rsquared_adj

"""
Explanation: Working with Hours.csv

Read in the 'hour.csv' file
Run the regression with: cnt ~ temp + hum + workingday + hour + C(weathersit)
Evaluate the results, how does this compare with the day 
Create a binary variable for rush hour defined by 6-9a & 4-6p
Run the regression again. Does this new variable improve the results?
End of explanation
"""


bike_hour_dat['rush']=0
bike_hour_dat['rush'][((bike_hour_dat['hr'] >= 6) & (bike_hour_dat['hr'] <= 9 ))|
            ((bike_hour_dat['hr'] >= 16) & (bike_hour_dat['hr'] <= 18 ))] = 1

est_hour_with_rush_m = smf.ols(formula='cnt ~ temp + hum + workingday + weathersit + rush', data=bike_hour_dat).fit()
est_hour_with_rush_m.summary()

# print (est_hour_m.rsquared_adj)
# print (est_m.rsquared_adj)
# print (est_hour_with_rush_m.rsquared_adj)

print ("Day                   :  [%0.6s]"%est_m.rsquared_adj)
print ("Hour                  :  [%0.6s]"%est_hour_m.rsquared_adj)
print ("Hour(With Rush hour)  :  [%0.6s]"%est_hour_with_rush_m.rsquared_adj)

"""
Explanation: Add new variable rush hour
End of explanation
"""


