import numpy as np
import pandas as pd
from quantopian.pipeline.data import morningstar
import statsmodels.api as sm
from statsmodels import regression
import matplotlib.pyplot as plt
import scipy.stats

"""
Explanation: Fundamental factor models
By Evgenia "Jenny" Nitishinskaya, Delaney Granizo-Mackenzie, and Maxwell Margenot.
Part of the Quantopian Lecture Series:

www.quantopian.com/lectures
github.com/quantopian/research_public

Notebook released under the Creative Commons Attribution 4.0 License.

Fundamentals are data having to do with the asset issuer, like the sector, size, and expenses of the company. We can use this data to build a linear factor model, expressing returns on any asset as
$$R_t = a_t + b_{t1} F_1 + b_{t2} F_2 + \ldots + b_{tK} F_K + \epsilon_t$$
There are two different approaches to computing the factors $F_j$, which represent the returns associated with some fundamental characteristics, and the factor sensitivities $b_{tj}$.
Approach 1: Portfolio Construction
In the first, we start by representing each characteristic of interest by a portfolio: we sort all assets by that characteristic, then build the portfolio by going long the top quantile of assets and short the bottom quantile. The factor corresponding to this characteristic is the return on this portfolio. Then, the $b_{ij}$ are estimated for each asset $i$ by regressing over the historical values of $R_i$ and of the factors.
We'll use the canonical Fama-French factors for this example, which are the returns of portfolios constructred based on fundamental factors.
We start by getting the fundamentals data for all assets and constructing the portfolios for each characteristic:
Import some libraries.
End of explanation
"""


start_date = '2011-1-1'
end_date = '2012-1-1'

"""
Explanation: Set the date range for which we want data.
End of explanation
"""


import numpy as np
from quantopian.pipeline import Pipeline
from quantopian.pipeline.data import morningstar
from quantopian.pipeline.data.builtin import USEquityPricing
from quantopian.pipeline.factors import CustomFactor, Returns

# Here's the raw data we need, everything else is derivative.

class MarketCap(CustomFactor):
    # Here's the data we need for this factor
    inputs = [morningstar.valuation.shares_outstanding, USEquityPricing.close]
    # Only need the most recent values for both series
    window_length = 1
    
    def compute(self, today, assets, out, shares, close_price):
        # Shares * price/share = total price = market cap
        out[:] = shares * close_price
        
        
class BookToPrice(CustomFactor):
    # pb = price to book, we'll need to take the reciprocal later
    inputs = [morningstar.valuation_ratios.pb_ratio]
    window_length = 1
    
    def compute(self, today, assets, out, pb):
        out[:] = 1 / pb
        
def make_pipeline():
    """
    Create and return our pipeline.
    
    We break this piece of logic out into its own function to make it easier to
    test and modify in isolation.
    
    In particular, this function can be copy/pasted into research and run by itself.
    """
    pipe = Pipeline()

    # Add our factors to the pipeline
    market_cap = MarketCap()
    # Raw market cap and book to price data gets fed in here
    pipe.add(market_cap, "market_cap")
    book_to_price = BookToPrice()
    pipe.add(book_to_price, "book_to_price")
    
    # We also get daily returns
    returns = Returns(inputs=[USEquityPricing.close], window_length=2)
    pipe.add(returns, "returns")
    
    # We compute a daily rank of both factors, this is used in the next step,
    # which is computing portfolio membership.
    market_cap_rank = market_cap.rank()
    pipe.add(market_cap_rank, 'market_cap_rank')
    
    book_to_price_rank = book_to_price.rank()
    pipe.add(book_to_price_rank, 'book_to_price_rank')

    # Build Filters representing the top and bottom 1000 stocks by our combined ranking system.
    biggest = market_cap_rank.top(1000)
    smallest = market_cap_rank.bottom(1000)
    
    highpb = book_to_price_rank.top(1000)
    lowpb = book_to_price_rank.bottom(1000)
    
    # Don't return anything not in this set, as we don't need it.
    pipe.set_screen(biggest | smallest | highpb | lowpb)
    
    # Add the boolean flags we computed to the output data
    pipe.add(biggest, 'biggest')
    pipe.add(smallest, 'smallest')
    
    pipe.add(highpb, 'highpb')
    pipe.add(lowpb, 'lowpb')
    
    return pipe

"""
Explanation: Using the Pipeline API to Fetch Data
The pipeline API is a very useful tool for factor analysis. We use it here to get data for our analysis. Specifically, we want the daily values of book to price ratio and market cap for every security. But we also do several other useful filtering steps which are detailed in code comments.
End of explanation
"""


pipe = make_pipeline()

"""
Explanation: Now we initialize the pipeline.
End of explanation
"""


pipe.show_graph('png')

"""
Explanation: We can visualize the dependency graph of our data computations here.
End of explanation
"""


from quantopian.research import run_pipeline

"""
Explanation: This function will allow us to run the pipeline.
End of explanation
"""


# This takes a few minutes.
results = run_pipeline(pipe, start_date, end_date)
results

"""
Explanation: Now let's actually run it and check out our results.
End of explanation
"""


results[results.biggest]

"""
Explanation: Great, we have all the data. Now we need to compute the returns of our portfolios over time. We have the daily returns for each equity, plus whether or not that equity was included in any given portfolio on any given day. We can combine that information in the following way to yield daily portfolio returns.
Step 1: Subset our results into only data belonging to our 'biggest' portfolio.
End of explanation
"""


results[results.biggest]['returns']

"""
Explanation: Step 2: Get returns.
End of explanation
"""


results[results.biggest]['returns'].groupby(level=0).mean()

"""
Explanation: Step 3: Group by day and take the mean. This is pretty deep into pandas logic, so if you don't understand this on first pass it is recommended to check out pandas' documentation on all the functions used. Especially groupby, which is very useful. Keep in mind that the index in our results is a MultiIndex rather than a regular Index, that can complicate things.
End of explanation
"""


R_biggest = results[results.biggest]['returns'].groupby(level=0).mean()
R_smallest = results[results.smallest]['returns'].groupby(level=0).mean()

R_highpb = results[results.highpb]['returns'].groupby(level=0).mean()
R_lowpb = results[results.lowpb]['returns'].groupby(level=0).mean()

SMB = R_smallest - R_biggest
HML = R_highpb - R_lowpb

"""
Explanation: Now run through this computation for each portfolio and get our final results.
End of explanation
"""


plt.plot(SMB.index, SMB.values)
plt.ylabel('Daily Percent Return')
plt.legend(['SMB Portfolio Returns']);

plt.plot(HML.index, HML.values)
plt.ylabel('Daily Percent Return')
plt.legend(['HML Portfolio Returns']);

"""
Explanation: What were the daily returns?
End of explanation
"""


plt.plot(SMB.index, np.cumprod(SMB.values+1))
plt.ylabel('Cumulative Return')
plt.legend(['SMB Portfolio Returns']);

"""
Explanation: And what would it look like to hold these portfolios over time?
End of explanation
"""


M = get_pricing('SPY', start_date='2011-1-1', end_date='2012-1-1', fields='price').pct_change()[1:]

plt.plot(M.index, M.values)
plt.ylabel('Daily Percent Return')
plt.legend(['Market Portfolio Returns']);

"""
Explanation: The last data we need are the daily returns on the broad market.
End of explanation
"""


# Get returns data for our portfolio
portfolio = get_pricing(['MSFT', 'AAPL', 'YHOO', 'FB', 'TSLA'], 
                        fields='price', start_date=start_date, end_date=end_date).pct_change()[1:]
R = np.mean(portfolio, axis=1)

"""
Explanation: Actually Running the Regression
Now that we have returns series representing our factors, we can compute the factor model for any return stream using a linear regression. Below, we compute the factor sensitivities for returns on a tech portfolio.
End of explanation
"""


# Define a constant to compute intercept
constant = pd.TimeSeries(np.ones(len(R.index)), index=R.index)

df = pd.DataFrame({'R': R,
              'M': M,
              'SMB': SMB,
              'HML': HML,
              'Constant': constant})
df = df.dropna()

"""
Explanation: Put all the data into one dataframe for convenience.
End of explanation
"""


# Perform linear regression to get the coefficients in the model
b1, b2, b3 = regression.linear_model.OLS(df['R'], df[['M', 'SMB', 'HML']]).fit().params

# Print the coefficients from the linear regression
print 'Historical Sensitivities of portfolio returns to factors:\nMarket: %f\nMarket cap: %f\nB/P: %f' %  (b1, b2, b3)

"""
Explanation: Perform the regression. You'll notice that these are the sensitivities over an entire year. It can be valuable to look at the rolling sensitivities as well to determine how stable they are.
End of explanation
"""


model = pd.stats.ols.MovingOLS(y = df['R'], x=df[['M', 'SMB', 'HML']], 
                             window_type='rolling', 
                             window=100)
rolling_parameter_estimates = model.beta
rolling_parameter_estimates.plot();
plt.title('Computed Betas');
plt.legend(['Market Beta', 'SMB Beta', 'HML Beta', 'Intercept']);

"""
Explanation: Let's perform a rolling regression to look at how the estimated sensitivities change over time.
End of explanation
"""


BTP = results['book_to_price']['2011-1-3']
zscore = (BTP - np.mean(BTP)) / np.std(BTP)
zscore.dropna(inplace=True)

plt.hist(zscore)
plt.xlabel('Z-Score')
plt.ylabel('Frequency');

"""
Explanation: Approach 2: Factor Value Normalization
This is also known as cross-sectional factor analysis.
Another approach is to normalize factor values each bar and see how predictive of that bar's returns they were. We do this by computing a normalized factor value $b_{aj}$ for each asset $a$ in the following way. 
$$b_{aj} = \frac{F_{aj} - \mu_{F_j}}{\sigma_{F_j}}$$
$F_{aj}$ is the value of factor $j$ for asset $a$ during this bar, $\mu_{F_j}$ is the mean factor value across all assets, and $\sigma_{F_j}$ is the standard deviation of factor values over all assets. Notice that we are just computing a z-score to make asset specific factor values comparable across different factors.
The exceptions to this formula are indicator variables, which are set to 1 for true and 0 for false. One example is industry membership: the coefficient tells us whether the asset belongs to the industry or not. 
After we calculate all of the normalized scores during bar $t$, we can estimate factor $j$'s returns $F_{jt}$, using a cross-sectional regression (i.e. at each time step, we perform a regression using the equations for all of the assets). Specifically, once we have returns for each asset $R_{at}$, and normalized factor coefficients $b_{aj}$, we construct the following model and estimate the $F_j$s and $a_t$
$$R_{at} = a_t + b_{a1}F_1 + b_{a2}F_2 + \dots + b_{aK}F_K$$
You can think of this as slicing through the other direction from the first analysis, as now the factor returns are unknowns to be solved for, whereas originally the coefficients were the unknowns. Another way to think about it is that you're determining how predictive of returns the factor was on that day, and therefore how much return you could have squeezed out of that factor.
Following this procedure, we'll get the cross-sectional returns on 2011-01-03, and compute the coefficients for all assets:
Getting the Data
We already have the results of the previous pipeline call, so we can grab book to price information for 2011-1-3 pretty easily.
End of explanation
"""


# Get some random data
X = np.random.normal(0, 1, 100)

# Put in some outliers
X[0] = 1000
X[1] = -1000

# Perform winsorization
print 'Before winsorization', np.min(X), np.max(X)
scipy.stats.mstats.winsorize(X, inplace=True, limits=0.01)
print 'After winsorization', np.min(X), np.max(X)

"""
Explanation: Problem: The Data is Weirdly Distributed
Notice how there are big outliers in the dataset that cause the z-scores to lose a lot of information. Basically the presence of some huge book to price datapoints causes the rest of the data to seem to occupy a relatively small area. We need to get around this issue using some data cleaning technique, here we're use winsorization.
Winsorization
Winzorization takes the top $n\%$ of a dataset and sets it all equal to the least extreme value in the top $n\%$. For example, if your dataset ranged from 0-10, plus a few crazy outliers, those outliers would be set to 0 or 10 depending on their direction. Here is an example.
End of explanation
"""


BTP = results['book_to_price']['2011-1-3']
scipy.stats.mstats.winsorize(BTP, inplace=True, limits=0.01)
BTP_z = (BTP - np.mean(BTP)) / np.std(BTP)
BTP_z.dropna(inplace=True)

plt.hist(BTP_z)
plt.xlabel('Z-Score')
plt.ylabel('Frequency');

"""
Explanation: This looks good, let's see how our book to price data looks when winsorized.
End of explanation
"""


R_day = results['returns']['2011-1-3']

"""
Explanation: We need the returns for that day as well.
End of explanation
"""


constant = pd.TimeSeries(np.ones(len(R_day.index)), index=R_day.index)

df_day = pd.DataFrame({'R': R_day,
              'BTP_z': BTP_z,
              'Constant': constant})
df_day = df_day.dropna()

# Perform linear regression to get the coefficients in the model
F1 = regression.linear_model.OLS(df_day['R'], df_day['BTP_z']).fit().params
print F1

"""
Explanation: Now set up our data and estimate $F_j$ using linear regression.
End of explanation
"""


MKT = results['market_cap']['2011-1-3']
scipy.stats.mstats.winsorize(MKT, inplace=True, limits=0.01)
MKT_z = (MKT - np.mean(MKT)) / np.std(MKT)

constant = pd.TimeSeries(np.ones(len(R_day.index)), index=R_day.index)

df_day = pd.DataFrame({'R': R_day,
              'BTP_z': BTP_z,
              'MKT_z': MKT_z,
              'Constant': constant})
df_day = df_day.dropna()

# Perform linear regression to get the coefficients in the model
F1, F2 = regression.linear_model.OLS(df_day['R'], df_day[['BTP_z', 'MKT_z']]).fit().params
print F1, F2

"""
Explanation: Finally, let's add another factor so you can see how the code changes.
End of explanation
"""


