# always yielding a real result, even dividing two integers
from __future__ import division

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

# plot inline
%matplotlib inline

"""
Explanation: Pricing and Hedging - Exercise Two
start coding
importing modules, and some basic stuff
End of explanation
"""


def mcpath(nsteps, S0, r, vol, T):
    
    # generate nsteps random numbers
    sample=pd.Series(np.random.standard_normal(nsteps))
    
    # define time increment
    dt=T/nsteps
    
    # initialize vectors
    initial = np.zeros(1)
    ts=np.repeat(dt,nsteps)
    
    # calculate s1 and s2
    s1=(r-(vol**2)/2)*dt
    s2=vol*np.sqrt(dt)
    steps=s1*np.repeat(1,nsteps)+s2*sample
    path=steps.cumsum()
    Sj=S0*np.exp(path)
    
    # add initial time t=0
    S = np.concatenate(([S0],Sj))
    t = np.concatenate(([0],ts.cumsum()))
    
    return pd.Series(S,index=t)

"""
Explanation: function for a monte carlo path sample, based on Generalized geometric Brownian motion:
$$S(t) = S(0)\exp\left {\sigma W(t) + (r - \frac{1}{2}\sigma^2)t\right }$$
End of explanation
"""


def mcpaths(npaths, nsteps, S0, r, vol, t):
    # generate npaths using mcpath
    paths=[mcpath(nsteps, S0, r, vol, t) for j in range(npaths)]
    return paths

"""
Explanation: An auxiliary function to generate npaths
End of explanation
"""


npaths=10
nsteps=100
S0=100
r=0.1
vol=0.2
T=1
test_mcpaths=mcpaths(npaths,nsteps,S0,r,vol,T)

[test_mcpaths[i].plot() for i in range(npaths)]

"""
Explanation: Let´s test
End of explanation
"""


def bsv(phi,S,K,r,q,vol,t):
    if(t>0):
        # calculate F
        F=S*np.exp((r-q)*t)

        # calculate d1 and d2
        sigma_rt = vol*np.sqrt(t)
        d1=(np.log(S/K)+(r-q+vol**2/2)*t)/(sigma_rt)
        d2=d1-sigma_rt

        # calculate N(d1) and N(d2)
        Nd1=st.norm.cdf(phi*d1)
        Nd2=st.norm.cdf(phi*d2)

        # calculate delta and premium
        delta=phi*np.exp(-q*t)*Nd1
        premium=S*delta-phi*K*np.exp(-r*t)*Nd2
    else:
        delta=0
        premium=max(phi*(S-K),0)
    
    return [premium,delta]

"""
Explanation: Defining the Black formula, according to:
$$C(S_0,t) = e^{-rT} \left[FN(d_1) - KN(d_2)\right]$$
Where:
$$F = S_0e^{\left(r-q\right)t}$$
$$d_1 = \frac{1}{\sigma\sqrt{t}} \left[ \ln{\left( \frac{S_0}{K} \right) } + \left( r -q + \frac{\sigma^2}{2} \right) \left(t \right) \right]$$
$$d_2 = d_1 - \sigma\sqrt{t}$$
End of explanation
"""


call = 1
put = -1
S0=100
r=0.1
q=.02
vol=0.2
T=1
K=100
[bsv(call,S0,K,r,q,vol,T),bsv(put,S0,K,r,q,vol,T)]

"""
Explanation: Testing the function we can verify that results match with this Black-Scholes Calculator from Drexel University
End of explanation
"""


def calcportfolio(path,phi,K,r,q,vol,T):
    nstp=len(path)-1

    # calculate t,S(t),premium(t),delta(t) using Black Scholes function
    bsvpath=[[path.index[i],path.values[i]]+bsv(call,path.values[i],K,r,q,vol,1-path.index[i])\
            for i in range(len(path))]
    df = pd.DataFrame(bsvpath,columns=['time','spot','premium','delta'])
       
    # generate cfwprem column to include cash flows regarding premium
    df['cfwprem']=0
    
    # the first cash flow is the price payed to enter the option
    df.loc[0,'cfwprem']=-df['premium'][0]
    
    # at the maturity:
    # - all premium, if any, becomes cash flow regarding option settlement
    # - there is no delta as option expires and we should withdraw the position on asset
    df.loc[nstp,'cfwprem']=df['premium'][nstp]
    df.loc[nstp,'premium']=0
    df.loc[nstp,'delta']=0
       
    # calculate time intervals dt
    df['timechg']=df['time'].diff()
    df.loc[0,'timechg']=0
    
    # calculate changes in delta
    df['dltchg']=df['delta'].diff()
    df.loc[0,'dltchg']=0
    
    # calculate changes in spot price
    df['spotchg']=df['spot'].diff()
    df.loc[0,'spotchg']=0

    # cashflows for heding the portfolio buying/selling delta quantities of the asset
    df['cfwspot']=0
    df.loc[0,'cfwspot']=df['delta'][0]*df['spot'][0]
    df.loc[1:,'cfwspot']=df['dltchg'][1:]*df['spot'][1:]

    # dividend cashflows
    df['cfwdivid']=0
    df.loc[1:,'cfwdivid']=-((df['delta'][0:nstp]*df['spot'][0:nstp]).values)*(np.exp(q*df['timechg'][1:].values)-1)

    # cashflows before interest
    df['cfwprer']=df['cfwprem']+df['cfwspot']+df['cfwdivid']

    # interest and consolidation of cashflows
    df['balance']=0
    df.loc[0,'balance']=df['cfwprer'][0]
    for j in range(1,nstp+1):
        df.loc[j,'balance']=df['balance'][j-1]*(np.exp(r*df['timechg'][j]))+df['cfwprer'][j]

    # portfolio
    df['portf']=df['premium']-df['delta']*df['spot']+df['balance']

    # consolidated discount factor
    return df

"""
Explanation: Function for calculating portfolio P&L and cash flow value given a path
End of explanation
"""


call=1
put=-1
nsteps=10
S0=100
r=0.1
q=0
vol=0.2
T=1
K=100
mypath = mcpath(nsteps,S0,r,vol,T)
calcportfolio(mypath,call,K,r,q,vol,T)

"""
Explanation: some function test, just to see the data frame´s face
End of explanation
"""


npaths=500
pls=pd.Series([calcportfolio(mcpath(nsteps, S0, r, vol, T),call,K,r,q,vol,T)['portf'][nsteps] for i in range(npaths)])
pls.plot(kind='hist',bins=30)

"""
Explanation: generating a histogram of portfolio P&L
End of explanation
"""


def bsv2(phi,S,K,r,q,vol_pricing,vol_delta,t):
    if(t>0):
        # calculate F
        F=S*np.exp((r-q)*t)

        # calculate d1 and d2
        sigma_rt = vol_pricing*np.sqrt(t)
        d1=(np.log(S/K)+(r-q+vol_pricing**2/2)*t)/(sigma_rt)
        d2=d1-sigma_rt

        # calculate N(d1) and N(d2)
        Nd1=st.norm.cdf(phi*d1)
        Nd2=st.norm.cdf(phi*d2)

        # calculate premium
        delta=phi*np.exp(-q*t)*Nd1
        premium=S*delta-phi*K*np.exp(-r*t)*Nd2

        # calculate delta
        sigma_rtdelta=vol_delta*np.sqrt(t)
        d1delta=(np.log(S/K)+(r-q+vol_delta**2/2)*t)/(sigma_rtdelta)
        Nd1delta=st.norm.cdf(phi*d1delta)
        delta=phi*np.exp(-q*t)*Nd1delta
    else:
        delta=0
        premium=max(phi*(S-K),0)
        
    return [premium,delta]

def calcportfolio2(path,phi,K,r,q,vol_pricing,vol_delta,T):
    nstp=len(path)-1

    # calculate t,S(t),premium(t),delta(t) using Black Scholes function
    bsvpath=[[path.index[i],path.values[i]]+bsv2(call,path.values[i],K,r,q,vol_pricing,vol_delta,1-path.index[i])\
            for i in range(len(path))]
    df = pd.DataFrame(bsvpath,columns=['time','spot','premium','delta'])
       
    # generate cfwprem column to include cash flows regarding premium
    df['cfwprem']=0
    
    # the first cash flow is the price payed to enter the option
    df.loc[0,'cfwprem']=-df['premium'][0]
    
    # at the maturity:
    # - all premium, if any, becomes cash flow regarding option settlement
    # - there is no delta as option expires and we should withdraw the position on asset
    df.loc[nstp,'cfwprem']=df['premium'][nstp]
    df.loc[nstp,'premium']=0
    df.loc[nstp,'delta']=0
       
    # calculate time intervals dt
    df['timechg']=df['time'].diff()
    df.loc[0,'timechg']=0
    
    # calculate changes in delta
    df['dltchg']=df['delta'].diff()
    df.loc[0,'dltchg']=0
    
    # calculate changes in spot price
    df['spotchg']=df['spot'].diff()
    df.loc[0,'spotchg']=0

    # cashflows for heding the portfolio buying/selling delta quantities of the asset
    df['cfwspot']=0
    df.loc[0,'cfwspot']=df['delta'][0]*df['spot'][0]
    df.loc[1:,'cfwspot']=df['dltchg'][1:]*df['spot'][1:]

    # dividend cashflows
    df['cfwdivid']=0
    df.loc[1:,'cfwdivid']=-((df['delta'][0:nstp]*df['spot'][0:nstp]).values)*(np.exp(q*df['timechg'][1:].values)-1)

    # cashflows before interest
    df['cfwprer']=df['cfwprem']+df['cfwspot']+df['cfwdivid']

    # interest and consolidation of cashflows
    df['balance']=0
    df.loc[0,'balance']=df['cfwprer'][0]
    for j in range(1,nstp+1):
        df.loc[j,'balance']=df['balance'][j-1]*(np.exp(r*df['timechg'][j]))+df['cfwprer'][j]

    # portfolio
    df['portf']=df['premium']-df['delta']*df['spot']+df['balance']

    # consolidated discount factor
    return df

"""
Explanation: Exercise 1
Change the the Black&Scholes and portfolio function to receive implied and realized volatity
End of explanation
"""


call=1
put=-1
S0=100
r=0
q=0
vol_pricing=.2
vol_delta=.3
T=1
K=100
nsteps=500
npaths=500

"""
Explanation: Define parameters
End of explanation
"""


paths = [mcpath(nsteps, S0, r, vol_delta, T) for i in range(npaths)]

"""
Explanation: Generate paths
End of explanation
"""


portfsVolReal=[calcportfolio2(paths[i],call,K,r,q,vol_pricing,vol_delta,T) for i in range(npaths)]
df=pd.DataFrame([pd.Series(portfsVolReal[i]['portf'].values,index=portfsVolReal[i]['time']) for i in range(npaths)])
df.transpose().plot(legend=None)

"""
Explanation: Hedging using Realized Vol
P&L portfolio envolution
End of explanation
"""


pls=pd.Series([portfsVolReal[i]['portf'][nsteps] for i in range(npaths)])
pls.plot(kind='hist',bins=10)

"""
Explanation: Final P&L distribution
End of explanation
"""


portfsVolImpl=[calcportfolio2(paths[i],call,K,r,q,vol_pricing,vol_pricing,T) for i in range(npaths)]
df=pd.DataFrame([pd.Series(portfsVolImpl[i]['portf'].values,index=portfsVolImpl[i]['time']) for i in range(npaths)])
df.transpose().plot(legend=None)

"""
Explanation: Hedging using Implied Vol
End of explanation
"""


pls=pd.Series([portfsVolImpl[i]['portf'][nsteps] for i in range(npaths)])
pls.plot(kind='hist',bins=10)

"""
Explanation: Final P&L distribution
End of explanation
"""


import scipy.optimize as opt

def bev(vol,path,phi,K,r,q,T):
    plfinal = calcportfolio(path,phi,K,r,q,vol,T)['portf'][len(path)-1]
    return plfinal

bevs = pd.Series([opt.newton(bev, 0.2, args=(paths[i],call,K,r,q,T)) for i in range(npaths)])

"""
Explanation: Exercise 2
End of explanation
"""


bevs.plot(kind='hist',bins=20)

def calcvolreal(path,T):
    return np.std(path.values[1:] / path.values[0:-1] - 1)*np.sqrt(len(path)/T)


vols_real = [calcvolreal(paths[i],T) for i in range(npaths)]

"""
Explanation: The Break Event Vol distribution
End of explanation
"""


plt.plot(vols_real,bevs.values,'o')
plt.xlabel('rvol')
plt.ylabel('bev')

"""
Explanation: Scatter plot between Realized Vol and Break Event Vol
End of explanation
"""


def strike_from_delta(phi,S,delta,r,q,vol,t):
    Nd1=delta/(phi*np.exp(-q*t))
    d1=st.norm.ppf(phi*Nd1)
    sigma_rt=vol*np.sqrt(t)
    K=S/(np.exp((d1*sigma_rt)-(r-q+vol**2/2)*t))
    return K

"""
Explanation: Exercise 3
Function that returns Strike value given a Delta value
End of explanation
"""


strikes=[strike_from_delta(call,S0,0.1,r,q,vol,T),\
         strike_from_delta(call,S0,0.25,r,q,vol,T),\
         S0,\
         strike_from_delta(put,S0,0.75,r,q,vol,T),\
         strike_from_delta(put,S0,0.9,r,q,vol,T)]
strikes

"""
Explanation: Calculate strike values for delta 25, delta 10, delta -10 and delta -25
End of explanation
"""


bevsmile = pd.DataFrame([pd.Series([opt.newton(bev, 0.3, args=(path,call,strike,r,q,T)) for strike in strikes], index=strikes) \
                         for path in paths[:20]])

"""
Explanation: Calculate the break even smile for the first 10 paths
End of explanation
"""


bevsmile.transpose().plot(legend=None)

"""
Explanation: Skew of break even vol
End of explanation
"""


plbystrike = pd.DataFrame([pd.Series([calcportfolio(path,call,strike,r,q,vol_pricing,T)['portf'][len(path)-1] \
                                      for strike in strikes], index=strikes) \
                           for path in paths[:5]])


plbystrike.plot()

plbystrike

"""
Explanation: Exercise 4
End of explanation
"""


import pandas.io as web
import pandas.io.data as webdata
from pandas.io.data import Options
import datetime

ticker="PETR4.SA"
start = datetime.datetime(2007, 1, 1)
end = datetime.datetime(2007,6,30)
bullQuotes = webdata.DataReader(ticker, 'yahoo', start, end)
start = datetime.datetime(2015, 7, 1)
end = datetime.datetime(2015,12,31)
bearQuotes = webdata.DataReader(ticker, 'yahoo', start, end)

plt.plot(bullQuotes['Close'])

plt.plot(bearQuotes['Close'])

"""
Explanation: Exercise 5
End of explanation
"""


