import numpy as np
import matplotlib.pyplot as plt
from numpy import pi, sqrt, exp, cumsum, sort, cos, diag, sum, linspace, sign, log
from numpy.random import randn, rand
from scipy.stats import norm
from numpy.linalg import cholesky, eig
from european_call import european_call

def wait():
    # python2.7 and python3 compatibility
    if hasattr(__builtins__, 'raw_input'):
        raw_input("Press Enter to continue...")
    else:
        input("Press Enter to continue...")

plt.rc('text', usetex=True)
plt.rc('font', family='serif', size=16)

plt.close("all")
plt.ion()

np.random.seed(123456789)

#
# Q1(a)
#

print('\nQ1(a)')

N = 10**6
U = rand(1,N)
X = randn(1,N)

mean = sum(U)/N
var  = U.dot(U.T)/N - mean**2
print("Uniform: mean = %f, var = %f" % (mean, var))

mean = sum(X)/N
var  = X.dot(X.T)/N - mean**2
print("Normal: mean = %f, var = %f" % (mean, var))

#
# Q1(b)
#

print('\nQ1(b)')

Sig = np.array([[4., 1.], [1., 4.]])

L = cholesky(Sig)
X = randn(2,N)
Y = L.dot(X)

mean = sum(Y,1)/N
var  = Y.dot(Y.T)/N - mean.dot(mean.T)

print("mean = %s \nvar =\n %s" % (mean, var))

# %timeit L.dot(randn(2,N))
# 92.1 ms ± 330 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
# so approximately 600 pairs can be generated in 1 minute


#
# Q1(c)
#

print('\nQ1(d)')

Sig = np.array([[4., 1.], [1., 4.]])

D,V = eig(Sig)
L = V.dot(diag(sqrt(D)))
Y = L.dot(X)

mean = sum(Y,1)/N
var  = Y.dot(Y.T)/N - mean.dot(mean.T)

print("mean = %s \nvar =\n %s" % (mean, var))

# in IPython: %timeit L.dot(randn(2,N))
# 102 ms ± 472 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
# so approximately 600 pairs can be generated in 1 minute

#
# Q2(b)
#

print('\nQ2(b)')

M = 1000
N = 1000

# exact values
fbar = - 2./pi**2 
sig2 = 1./(4.*pi**2) + 1./6. - fbar**2

U = rand(N,M)
f = U*cos(pi*U)

#
# Q2(c)
#

print('\nQ2(c)')

Y = sort(sum(f,0)/N) # the sum here is a numpy.sum
C = (np.arange(1,M+1) - 0.5)/M

plt.figure()
plt.plot(Y,C, fbar+sqrt(sig2/N)*norm.ppf(C),C)
plt.xlabel(r'$Y$'); plt.ylabel(r'$C$')
plt.legend(('Monte Carlo results','Central Limit Theorem'), loc = 2, fontsize = 14)

plt.figure()
plt.plot(C,Y, C,fbar+sqrt(sig2/N)*norm.ppf(C))
plt.ylabel(r'$Y$'); plt.xlabel(r'$C$')
plt.legend(('Monte Carlo results','Central Limit Theorem'), loc = 2, fontsize = 14)


#
# Q2(d)
#

print('\nQ2(d)')

N = 10**6
U = rand(1,N)
f = U*cos(pi*U)

S1 = cumsum(f)
S2 = cumsum(f**2)

NN   = np.arange(10**3,10**6+1,10**3)
mean = S1[NN-1]/NN
sd   = sqrt( (S2[NN-1]/NN - mean**2)/(NN-1) )

plt.figure()
plt.plot(NN,mean, NN,mean-3*sd, NN,mean+3*sd, NN,fbar*np.ones_like(NN))
plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0))
plt.xlabel(r'$N$'); plt.ylabel('value')
plt.legend(('mean','lower','upper','exact'))

plt.figure()
plt.plot(NN,NN*sd**2, NN,sig2*np.ones_like(NN))
plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0))
plt.xlabel('N'); plt.ylabel('variance')
plt.legend(('estimated','exact'))


#
# Q3(b)
#

print('\nQ3(b)')

M = 1000
N = 1000

r     = 0.05
sigma = 0.2
T     = 1.
K     = 100.
S0    = 100.

U = rand(N,1000)
W = sqrt(T)*norm.ppf(U)
f = exp(-r*T)*np.maximum(0,S0*exp((r-0.5*sigma**2)*T+sigma*W)-K)

fbar = european_call(r,sigma,T,S0,K,'value')
sig2 = sum(sum((f-fbar)**2)) / (M*N)

#
# Q3(c)
#

print('\nQ3(c)')

Y = sort(sum(f,0)/N)
C = (np.arange(1,M+1) - 0.5)/M

plt.figure()
plt.plot(Y,C, fbar+sqrt(sig2/N)*norm.ppf(C),C)
plt.xlabel(r'$Y$'); plt.ylabel(r'$C$')
plt.legend(('Monte Carlo results','Central Limit Theorem'),loc=2, fontsize=14)

plt.figure()
plt.plot(C,Y, C,fbar+sqrt(sig2/N)*norm.ppf(C))
plt.ylabel(r'$Y$'); plt.xlabel(r'$C$')
plt.legend(('Monte Carlo results','Central Limit Theorem'),loc=2, fontsize=14)

#
# Q3(d)
#

print('\nQ3(d)\n')

N = 10**6
U = rand(1,N)
W = sqrt(T)*norm.ppf(U)
f = exp(-r*T)*np.maximum(0,S0*exp((r-0.5*sigma**2)*T+sigma*W)-K)

S1 = cumsum(f)
S2 = cumsum(f**2)

NN   = np.arange(10**3,10**6+1,10**3)
mean = S1[NN-1]/NN
sd   = sqrt( (S2[NN-1]/NN - mean**2)/(NN-1) )

plt.figure()
plt.plot(NN,mean, NN,mean-3*sd, NN,mean+3*sd, NN,fbar*np.ones_like(NN))
plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0))
plt.xlabel(r'$N$'); plt.ylabel(r'$value$')
plt.legend(('mean','lower','upper','exact'))

wait()
plt.close("all")

#
# Q4(a)
#

print('\nQ4(a)\n')

N = 10**6

r     = 0.05
sigma = 0.2
T     = 1.0
K     = 100.0
S0    = 100.0

W  = sqrt(T)*randn(N,1)
fp = exp(-r*T)*np.maximum(0.0,S0*exp((r-0.5*sigma**2)*T+sigma*W)-K)
fm = exp(-r*T)*np.maximum(0.0,S0*exp((r-0.5*sigma**2)*T-sigma*W)-K)

fbar = np.sum(fp+fm)/(2.0*N)
sig2 = np.sum(fp**2+fm**2)/(2.0*N) - fbar**2
cov  = np.sum((fp-fbar)*(fm-fbar))/N
corr = cov/sig2

variance_reduction_factor = 0.5*(1.0+corr)

print('correlation:       %f' % corr)
print('variance reduction: %f' % variance_reduction_factor)


#
# Q4(b)
#

print('\nQ4(b)\n')

N = 10**6

r     = 0.05
sigma = 0.2
T     = 1.0
K     = 100.0
S0    = 100.0

W = sqrt(T)*randn(N,1)
S = S0*exp((r-0.5*sigma**2)*T+sigma*W)
f = exp(-r*T)*np.maximum(0.0,S-K)
g = exp(-r*T)*S

fbar  = np.sum(f)/N
gbar  = np.sum(g)/N
sig2f = np.sum(f**2)/N - fbar**2
sig2g = np.sum(g**2)/N - gbar**2
cov   = np.sum(f*g)/N - fbar*gbar
corr  = cov / sqrt(sig2f*sig2g)

variance_reduction_factor = 1.0-corr**2

print('correlation:       %f' % corr)
print('variance reduction: %f' % variance_reduction_factor)

#
# Q5(a)
#

print('\nQ5(a)\n')

N = 10**6

r     = 0.05
sigma = 0.2
T     = 1.0
K     = 50.0
S0    = 100.0

W = sqrt(T)*randn(N,1)
S = S0*exp((r-0.5*sigma**2)*T+sigma*W)
f = exp(-r*T)*0.5*(1.0+sign(K-S))

fbar = np.sum(f)/N
sig2 = np.sum(f**2)/N - fbar**2

fexact = exp(-r*T)*norm.cdf((log(K)-(log(S0)+(r-0.5*sigma**2)*T))/(sigma*sqrt(T)))

N_required = 9.0*sig2 / (0.1*fbar)**2

print('numerical value: %f  +/- %f' % (fbar,3.0*sqrt(sig2/N)))
print('analytic value:  %f' % fexact)
print('number of paths required: %d' % N_required)

#
# Q5(b)
#

print('\nQ5(b)\n')

N = 10**6

r     = 0.05
sigma = 0.2
T     = 1.0
K     = 50.0
S0    = 100.0

W     = sqrt(T)*randn(N,1)
mu1   = (r-0.5*sigma**2)*T
mu2   = log(K) - log(S0)
shift = mu2 - mu1
R     = exp(-(0.5*shift+sigma*W)*shift/(sigma**2.0*T))
S     = K*exp(sigma*W)
f     = exp(-r*T)*0.5*(1+sign(K-S)) * R

fbar = np.sum(f)/N
sig2 = np.sum(f**2)/N - fbar**2

fexact = exp(-r*T)*norm.cdf((log(K)-(log(S0)+(r-0.5*sigma**2)*T))/(sigma*sqrt(T)))

N_required = 9.0*sig2 / (0.1*fbar)**2

print('numerical value: %f  +/- %f' % (fbar,3.0*sqrt(sig2/N)))
print('analytic value:  %f' % fexact)
print('number of paths required: %d' % N_required)

wait()
plt.close("all")

#
# Q6(a)
#

print('\nQ6(a) (bumping) - this might take some time...')

np.random.seed(123456789) # reset random number generator

N  = 10**6
S0 = 100.0
r  = 0.05
sigma = 0.2
K = 100.0
T = 1.0

for option in [1,2]:
    if option == 1:
        gexact = european_call(r,sigma,T,S0,K,'delta')
        bumps  = 10.**linspace(-1,1,20)
    else:
        gexact = european_call(r,sigma,T,S0,K,'vega')
        bumps  = 10.**linspace(-3,-1,20)
  
    for it in [1,2]:
        vals  = []
        vars  = []

        for bump in bumps:
            if it == 2:
                np.random.seed(123456789) # reset random number generator

            W = sqrt(T)*randn(N,1)
            if option == 1:
                S = (S0-bump)*exp((r-0.5*sigma**2)*T+sigma*W)
            else:
                S = S0*exp((r-0.5*(sigma-bump)**2)*T+(sigma-bump)*W)

            fm = exp(-r*T)*np.maximum(0.0,S-K)

            if it == 2:
                np.random.seed(123456789) # reset random number generator

            W = sqrt(T)*randn(N,1)
            if option == 1:
                S = (S0+bump)*exp((r-0.5*sigma**2)*T+sigma*W)
            else:
                S = S0*exp((r-0.5*(sigma+bump)**2)*T+(sigma+bump)*W)

            fp = exp(-r*T)*np.maximum(0.0,S-K)

            Del = (fp-fm)/(2.0*bump)

            val =  np.sum(Del)/N
            var = (np.sum(Del**2)/N - val**2)/(N-1.0)

            vals.append(val)
            vars.append(var)

        fig = plt.figure(2*(option-1)+it)

        plt.subplot(2,1,1)
        plt.plot(bumps,vals-gexact,'-x')
        plt.xlabel('bump')
        if option == 1:
            plt.ylabel('Delta error')
        else:
            plt.ylabel('Vega error')

        if it == 1:
            plt.title('different random numbers')
        else:
            plt.title('same random numbers')

        plt.subplot(2,1,2)
        plt.plot(bumps,sqrt(vars),'-x');
        plt.xlabel('bump'); plt.ylabel('s.d.')

        fig.tight_layout()

        if it == 1:
            plt.savefig('assignment1_1.eps', format='eps', dpi=1200)
        else:
            plt.savefig('assignment1_2.eps', format='eps', dpi=1200)


#
# Q6(b)
#

print('\nQ6(b) (pathwise)\n')

N = 10**6

r     = 0.05
sigma = 0.2
T     = 1.0
K     = 100.0
S0    = 100.0

W = sqrt(T)*randn(N,1)
S = S0*exp((r-0.5*sigma**2)*T+sigma*W)

fp = exp(-r*T)*0.5*(1.0+sign(S-K))

#
# delta
#

g = fp * S/S0

gbar = np.sum(g)/N
sig2 = np.sum(g**2)/N - gbar**2

gexact = european_call(r,sigma,T,S0,K,'delta')

print('delta\n')
print('numerical value: %f  +/- %f' % (gbar,3.0*sqrt(sig2/N)))
print('analytic value:  %f' % gexact)

#
# vega
#

g = fp * S*(W-sigma*T)

gbar = np.sum(g)/N
sig2 = np.sum(g**2)/N - gbar**2

gexact = european_call(r,sigma,T,S0,K,'vega')

print('\nvega\n')
print('numerical value: %f  +/- %f' % (gbar,3.0*sqrt(sig2/N)))
print('analytic value:  %f' % gexact)

#
# Q6 -- extra (LRM)
#

print('\nQ6 -- extra (LRM)\n')

f = exp(-r*T)*np.maximum(0.0,S-K)

#
# delta
#

g = f * W/(S0*sigma*T)

gbar = np.sum(g)/N
sig2 = np.sum(g**2)/N - gbar**2

gexact = european_call(r,sigma,T,S0,K,'delta')

print('delta\n')
print('numerical value: %f  +/- %f' % (gbar,3.0*sqrt(sig2/N)))
print('analytic value:  %f' % gexact)

#
# vega
#

g = f * ( (W**2-T)/(sigma*T) - W)

gbar = np.sum(g)/N
sig2 = np.sum(g**2)/N - gbar**2

gexact = european_call(r,sigma,T,S0,K,'vega')

print('\nvega\n')
print('numerical value: %f  +/- %f' % (gbar,3.0*sqrt(sig2/N)))
print('analytic value:  %f' % gexact)

wait()