# -*- coding: utf-8 -*-
"""
simulation of 

- a brownian motion and
- Riemann-Ito sums
- log-normal dynamics 
- risk-neutral dynamics
- computing the price of an option
@author: turinici
"""
import numpy as np
import matplotlib.pyplot as plt
T,N=1,500
h=T/N
M=1000
trange=np.linspace(0,T,N+1,endpoint=True)

W=np.zeros((N+1,M))
W[1:,:]= np.sqrt(h)*np.cumsum(np.random.randn(N,M),
                              axis=0)

exact_result=(W[-1,:]**2-T)/2
riemann_ito_sums=np.sum(W[:-1,:]*(W[1:,:]-W[:-1,:]),
                        axis=0)

_=plt.figure('W',figsize=(10,5))
plt.subplot(1,2,1)
_=plt.plot(trange,W[:,:np.minimum(49,M-1)])
plt.subplot(1,2,2)
_=plt.hist(exact_result-riemann_ito_sums,
           bins=int(np.sqrt(M)))

## compute the historical evolution for S_t
S0=100
mu,sigma=0.1,0.25
S=S0*np.exp( (mu-sigma**2/2)*trange[:,None] + sigma*W)
plt.figure('S')
_=plt.plot(trange,S[:,:np.minimum(49,M-1)])


#compute the risk-neutral evolution to obtain the option price
r=0.05
K=110
M=10000
W_risk_neutral=np.zeros((N+1,M))
W_risk_neutral[1:,:]=np.sqrt(h)*np.cumsum(np.random.randn(N,M),
                              axis=0)

S_T_risk_neutral=S0*np.exp( (r-sigma**2/2)*T 
                               + sigma*W_risk_neutral[-1,:])
#compute price as average
results=np.exp(-r*T)*np.maximum(S_T_risk_neutral-K,0)
price=np.mean(results)
print('price=',price)


def bootstrap_mean_confidence_interval(data, num_iterations=10000, alpha=0.05):
    """
    A function to compute the average and a confidence interval around it.

    Use example 
    data = np.array([0.2, 0.5, 0.7, 0.8, 1.1, 1.3, 1.5, 1.8, 2.0, 2.2])
    print(bootstrap_mean_confidence_interval(data))
    

    Parameters
    ----------
    data : 1D array of data
    num_iterations : number of bootstrap iterations, default is 10000.
    alpha : confidence level, the default is 0.05.

    Returns : the average and a confidence interval around it as a tuple
    -------
    """

    n = len(data)
    means = np.zeros(num_iterations)

    for i in range(num_iterations):
        means[i] = np.mean(np.random.choice(data, size=n, replace=True))

    means.sort()
    lower_bound = means[int(num_iterations * (alpha / 2))]
    upper_bound = means[int(num_iterations * (1 - alpha / 2))]
    mean_estimate = np.mean(means)

    return mean_estimate, (lower_bound, upper_bound)

print('price, CI=',price,bootstrap_mean_confidence_interval(results))