import numpy as np

def bsm(S0,r,sigma,T,K,R = 100000 , seed=500):
    np.random.seed(seed)
    z = np.random.standard_normal(R) 
    ST = S0 * np.exp(( r - 0.5 * sigma ** 2) * T + sigma * np.sqrt(T) * z)
    hT = np.maximum(ST - K, 0)
    C0 = np.exp(-r * T) * np.sum(hT) / R
    return C0

import time
tm = time.time()
C0 = bsm(S0=105,r=0.06,sigma=0.22,T=1.0,K=109,R = 100000 , seed=500)
pm = time.time() - tm
print("Value of European Call Option: {0:.4g}".format(C0)+" - time[{0:.4g} secs]".format(pm))

"""
Explanation: Warming up
Value of an European Call Option
The Black-Scholes-Merton model is a 
mathematical model applicable in financial markets containing derivative investment instruments for 
the theoretical estimation of European-style options' price. According to such model, 
the index level at maturity of a call option can be modeled as it follows with z normally distributed
$$S_T(z)=S_0e^{(r-0.5\sigma^2)T+\sigma\sqrt{T}z}$$
where 
* $S_0$ is the initial stock index level
* K is the strike price of the European call option
* T is the Time-to-maturity in year 
* r is the riskless short rate 
* $\sigma^2$ is the volatility
Let's assume the following values 
* $S_0 = 105$
* K = 109 
* T = 1 
* r = 6% 
* $\sigma^2$ = 22%
In order to calculate the present value of the call option, it is possible to apply a Monte Carlo approach
according to the following procedure: 
* Let's consider R realizations of the random variable z from a standard normal distribution 
* Compute related index levels at maturity $S_T(i)$ for $z(i)\in{1,2,...,R}$ 
* Compute the call option value $h_T(i)=\max{(S_T(i)–K,0)}$. 
* Compute the present value of the call option, i.e. 
$$C_0=e^{-rT}{R}^{-1}\sum_{i=1}^{R} h_T{(i)}$$
End of explanation
"""


tm = time.time()
C0 = bsm(S0=105,r=0.06,sigma=0.22,T=1.0,K=109,R = 70000000 , seed=500)
pm = time.time() - tm
print("Value of European Call Option: {0:.4g}".format(C0)+" - time[{0:.4g} secs]".format(pm))

"""
Explanation: Let's see how much time is necessary for 70,000,000 iterations intead of 100,000 iterations.
End of explanation
"""


import numexpr as ne
def bsm_ne(S0,r,sigma,T,K,R = 70000000 , seed=500):
    np.random.seed(seed)
    z = np.random.standard_normal(R) 
    ST = ne.evaluate('S0 * exp(( r - 0.5 * sigma ** 2) * T + sigma * sqrt(T) * z)')
    hT = np.maximum(ST - K, 0)
    C0 = np.exp(-r * T) * np.sum(hT) / R
    return C0

tm = time.time()
C0 = bsm_ne(S0=105,r=0.06,sigma=0.22,T=1.0,K=109,R = 70000000 , seed=500)
pm = time.time() - tm
print("Value of the European Call Option: {0:.4g}".format(C0)+" - time[{0:.4g} secs]".format(pm))

"""
Explanation: Let's see how we can speed up the computation with the numexpr package.
End of explanation
"""


import numpy as np
import pandas as pd
import pandas.io.data as web

df_final = web.DataReader(['GOOG','SPY'], data_source='yahoo',
                      start='1/21/2010', end='4/15/2016')

print(df_final)
print(df_final.shape)

df_final.ix[:,:,'SPY'].head()

print(type(df_final.ix[:,:,'SPY']))
print("\n>>> null values:"+str(pd.isnull(df_final.ix[:,:,'GOOG']).sum().sum()))

df_final = web.DataReader(['GOOG','SPY'], data_source='yahoo',
                      start='1/21/1999', end='4/15/2016')
df_final.ix[:,:,'GOOG'].head()

print(type(df_final.ix[:,:,'GOOG']))
print("\n>>> null values:"+str(pd.isnull(df_final.ix[:,:,'GOOG']).sum().sum()))

"""
Explanation: Key Factors for Evaluating the Performance of a Portfolio
The daily return of a stock is easily computable as it follows  $$dr(t)=\frac{P(t)}{P(t-1)}-1$$
Similarly, the cumulative return of a stock  is easily computable as it follows $$cr(t)=\frac{P(t)}{P(0)}-1$$
What is it P(t)? There are basically 2 options for this value, i.e. 
* the adjusted close price of a stock typically indicated in financial feed as Adj Close, or 
* the close price of a stock typically indicated in financial feed as Adj. 
We take the adjusted close price 
(see What Hedge Funds Really Do
to understand why).     
Typically for evaluating the performance of a portfolio key factors to focus on are 
1. Cumulative return 
2. Average daily return 
3. Rsk (Standard deviation of daily return)
4. Sharpe ratio 
We will see how to compute and plot these factors. 
Get financial data
Functions from pandas.io.data and pandas.io.ga extract data from various Internet sources into a DataFrame
The following sources are supported:

Yahoo! Finance
Google Finance
St.Louis FED (FRED)
Kenneth French’s data library
World Bank
Google Analytics

For further info see pandas documentation
End of explanation
"""


import matplotlib.pyplot as plt
def get_data(symbols, 
             add_ref=True,
             data_source='yahoo',
             price='Adj Close',
             start='1/21/2010', 
             end='4/15/2016'):
    """Read stock data (adjusted close) for given symbols from."""
    
    if add_ref and 'SPY' not in symbols:  # add SPY for reference, if absent
        symbols.insert(0, 'SPY')

    df = web.DataReader(symbols, 
                        data_source=data_source,
                        start=start, 
                        end=end)
    
    return df[price,:,:]


get_data(symbols=['GOOG','SPY']).tail()

"""
Explanation: There is a couple of observations to be done:
1. calling pandas.io.data with multiple stocks gets a pandas.core.panel.Panel instead of a pandas.DataFrame
but filtering to specific axis (e.g. Google) we get pandas.core.frame.DataFrame
2. pandas.io.data does not handle missing values 
Hence, we can define the following functions.
End of explanation
"""


df_stock = get_data(symbols=['GOOG','SPY'],start='1/21/1999',end='4/15/2016')
print(">> Trading days from pandas:"+str(df_stock.shape[0]))
dates = pd.date_range('1/21/1999', '4/15/2016')
df = pd.DataFrame(index=dates)
print(">> Calendar days:"+str(df.shape[0]))
df = df.join(df_stock)
print(">> After join:"+str(df.shape[0]))
df = df.dropna(subset=["SPY"])
print(">> After removing non trading days:"+str(df.shape[0]))

"""
Explanation: Also, notice that it is not necessary to perform an initial join with the data range of interest filtering out 
non trading days, as padas does it for us , i.e.
End of explanation
"""


ax = get_data(symbols=['GOOG','SPY','IBM','GLD'],start='1/21/1999', end='4/15/2016').plot(title="Stock Data", fontsize=9)
ax.set_xlabel("Date")
ax.set_ylabel("Price")
plt.show()

"""
Explanation: Plotting stock prices
End of explanation
"""


def fill_missing_values(df_data):
    """Fill missing values in data frame, in place."""
    df_data.fillna(method='ffill',inplace=True)
    df_data.fillna(method='backfill',inplace=True)
    return df_data

ax = fill_missing_values(get_data(symbols=['GOOG','SPY','IBM','GLD'],
                             start='1/21/1999', 
                             end='4/15/2016')).plot(title="Stock Data", fontsize=9)
ax.set_xlabel("Date")
ax.set_ylabel("Price")
plt.show()

"""
Explanation: Imputing missing values
As clear from above plot, we need to handle missing values.
End of explanation
"""


def normalize_data(df):
    return df/df.ix[0,:]

ax = normalize_data(
    fill_missing_values(
        get_data(
            symbols=['GOOG','SPY','IBM','GLD'],
                             start='1/21/1999', 
                             end='4/15/2016'))).plot(title="Stock Data", fontsize=9)
ax.set_xlabel("Date")
ax.set_ylabel("Normalized price")
plt.show()

"""
Explanation: Normalizing prices
End of explanation
"""


df = fill_missing_values(
        get_data(
            symbols=['GOOG','SPY','IBM','GLD'],
                             start='4/21/2015', 
                             end='7/15/2016'))

# 1. Computing rolling mean using a 20-day window 
rm_df = pd.DataFrame.rolling(df, window=20).mean()

ax = rm_df.plot(title="Rolling Mean")
ax.set_xlabel("Date")
ax.set_ylabel("Price")
plt.show()

# 2. Computing rolling standard deviation using a 20-day window 
rstd_df = pd.DataFrame.rolling(df, window=20).std()

ax = rstd_df.plot(title="Rolling Standard Deviation")
ax.set_xlabel("Date")
ax.set_ylabel("Price")
plt.show()


# 3. Compute upper and lower bands
def get_bollinger_bands(rm, rstd):
    """Return upper and lower Bollinger Bands."""
    upper_band, lower_band = rm + 2 * rstd, rm - 2 * rstd
    return upper_band, lower_band

upper_band, lower_band = get_bollinger_bands(rm_df, rstd_df)

# Plot raw SPY values, rolling mean and Bollinger Bands
ax = df['SPY'].plot(title="Bollinger Bands",label='SPY')
rm_df['SPY'].plot(label='Rolling mean', ax=ax)
upper_band['SPY'].plot(label='upper band', ax=ax)
lower_band['SPY'].plot(label='lower band', ax=ax)

# Add axis labels and legend
ax.set_xlabel("Date")
ax.set_ylabel("Price")
ax.legend(loc='lower left')
plt.show()

"""
Explanation: Rolling statistics

Notice that pandas.rolling_mean has been deprecated for DataFrame and will be removed in a future version.
Hence, we will replace it with DataFrame.rolling(center=False,window=20).mean() 
Notice that pd.rolling_std has been deprecated for DataFrame and will be removed in a future version.
Hence, we will replace it with DataFrame.rolling(center=False,window=20).std()
End of explanation
"""


def compute_daily_returns_2(df):
    """Compute and return the daily return values."""
    # Note: Returned DataFrame must have the same number of rows
    daily_returns = df.copy()
    daily_returns[1:] = (df[1:]/df[:-1].values) - 1
    daily_returns.ix[0,:] = 0 
    return daily_returns

def compute_daily_returns(df):
    """Compute and return the daily return values."""
    # Note: Returned DataFrame must have the same number of rows
    daily_returns = (df / df.shift(1)) - 1 
    daily_returns.ix[0,:] = 0 
    return daily_returns

pd.util.testing.assert_frame_equal(
    compute_daily_returns(
        fill_missing_values(
            get_data(
                symbols=['GOOG','SPY','IBM','GLD'],
                start='4/21/2016', 
                end='7/15/2016'))), 
    compute_daily_returns_2(
        fill_missing_values(
            get_data(
                symbols=['GOOG','SPY','IBM','GLD'],
                start='4/21/2016', 
                end='7/15/2016'))))

ax = compute_daily_returns(fill_missing_values(
        get_data(
            symbols=['GOOG','SPY','IBM','GLD'],
                             start='4/21/2016', 
                             end='7/15/2016'))).plot(title="Daily returns")

ax.set_xlabel("Date")
ax.set_ylabel("Daily return")
plt.show()

df = compute_daily_returns(fill_missing_values(get_data(
            symbols=['SPY'], 
            start='4/21/2000', 
            end='7/15/2016')))

plt.hist(df['SPY'],bins=30,color='c',label=['Daily return'])
plt.axvline(df['SPY'].mean(), color='b', linestyle='dashed', linewidth=2 , label='Mean')
plt.axvline(-df['SPY'].std(), color='r', linestyle='dashed', linewidth=2 , label='Std')
plt.axvline(df['SPY'].std(), color='r', linestyle='dashed', linewidth=2 )
plt.title('SPY daily return distribution')
plt.xlabel('Daily return')
plt.grid(True)
plt.legend()

plt.show()

"""
Explanation: Daily returns
There are two ways to compute the daily return of a stock with pandas. We check they produce same 
results and plot them.
End of explanation
"""


def cumulative_returns(df):
    return df/df.ix[0,:] - 1 

ax = cumulative_returns(fill_missing_values(
        get_data(
            symbols=['GOOG','SPY','IBM','GLD'],
                             start='4/21/2016', 
                             end='7/15/2016'))).plot(title="Cumulative returns")

ax.set_xlabel("Date")
ax.set_ylabel("Cumulative return")
plt.show()

"""
Explanation: Cumulative returns
End of explanation
"""


def sharpe_ratio(df,sample_freq='d',risk_free_rate=0.0):
    
    sr = (df - risk_free_rate).mean() / df.std()
    
    if sample_freq == 'd':
        sr = sr * np.sqrt(252)
    elif sample_freq == 'w':
        sr = sr * np.sqrt(52)
    elif sample_freq == 'm':
        sr = sr * np.sqrt(12)
    else:
        raise Exception('unkown sample frequency :'+str(sample_freq))
        
    return sr

# Sharpe ratio
sharpe_ratio(
    compute_daily_returns(
        fill_missing_values(
            get_data(
                symbols=['GOOG','SPY','IBM','GLD'],
                             start='4/21/2015', 
                             end='7/15/2016'))))

"""
Explanation: Sharpe Ratio
Sharpe ratio is a 
way to examine the performance of an investment by adjusting for its risk. The ratio measures the excess 
return (or risk premium) per unit of deviation in an investment asset or a trading strategy, typically 
referred to as risk (and is a deviation risk measure), named after William F. Sharpe.
The ex-ante Sharpe ratio is defined as:
    $$S=\frac{E{R_p-R_f}}{Std{R_p-R_f}}$$
where $R_p$ is the asset return, $R_f$ is the return on a benchmark asset, 
such as the risk free rate of return. $E{R_p-R_f}$ is the expected value of the excess of the asset return 
over the benchmark return, and $Std{R_p-R_f}$ is the standard deviation of the asset return.
The ex-post Sharpe ratio uses the same equation as the one above but with 
realized returns of the asset and benchmark rather than expected returns.
Examples of risk free rates are 
* LIBOR
* Interest rate of 3 month treasury bill 
* 0% (approximation used a lot in last years)
Using daily returns as $R_p$  the ex-post Sharpe ratio is computable as:
    $$S=\frac{mean{R_p^{daily}-R_f^{daily}}}{Std{R_p^{daily}}}$$
as $R_f^{daily}$ is typically costant for several months it is possible to remove from the denominator. 
Also, $R_f^{daily}$ is typically approximable with 0% but it is possible to compute given $R_f^{yearly}$ as 
it follows remembering there are 252 trading days in a year
$$R_f^{daily}=\sqrt[252]{1+R_f^{yearly}}$$
Sharpe ratio is typically an annual measure. This means if we are using different sample frequencies 
need to apply the following formula:
    $$S^{annual}=\sqrt[2]{SPY} \times S$$
where $SPR$ is the number of samples per year considered in computing $S$, e.g. 

$S^{annual}=\sqrt[2]{252} \times S^{daily}$
$S^{annual}=\sqrt[2]{52} \times S^{weekly}$
$S^{annual}=\sqrt[2]{12} \times S^{monthly}$
End of explanation
"""


df = fill_missing_values(get_data(symbols=['GOOG','SPY','IBM','GLD'],
                             start='4/21/2015', 
                             end='7/15/2016'))

# 1. Cumulative return
cumulative_returns(df).ix[-1,:]

# 2. Average daily return
compute_daily_returns(df).mean()

# 3. Rsk (Standard deviation of daily return)
compute_daily_returns(df).std()

# 4. Sharpe ratio
sharpe_ratio(compute_daily_returns(df))

"""
Explanation: Summary
For evaluating the performance of a portfolio key factors to focus on are 
1. Cumulative return 
2. Average daily return 
3. Rsk (Standard deviation of daily return)
4. Sharpe ratio
End of explanation
"""


