import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(0,10,1000)

f = np.exp(-x)/(1+(x-1)**2)

plt.figure(figsize=(8,4), dpi=100)
plt.plot(x,f)
plt.grid()
plt.xlabel(r'$x$')
plt.ylabel(r'f(x)')
plt.show()

y = 6
K = 100
m = np.logspace(1,8,K, dtype='int')
estim = np.zeros(K)

for k in range(K):
  xs = np.random.rand(m[k])*y
  estim[k] = np.mean(np.exp(-xs)/(1+(xs-1)**2))*y

plt.figure(figsize=(8,4), dpi=100)
plt.semilogx(m,estim)
plt.grid()
plt.xlabel('Number of samples')
plt.ylabel('Integral estimation')
plt.show()

plt.figure(figsize=(8,4), dpi=100)
plt.semilogx(m,estim, label='Estimate')
plt.semilogx(m,0.6960918269*np.ones(K), ':', linewidth=2, label='Accurate')
plt.grid()
plt.legend()
plt.xlabel('Number of samples')
plt.ylabel('Integral estimation')
plt.show()

import matplotlib.pyplot as plt

plt.figure(figsize=(6,6), dpi=100)
plt.gcf().gca().add_artist(plt.Circle((0, 0), 1))
plt.xlabel(r'$x_1$')
plt.ylabel(r'$x_2$')
plt.axis('square')
plt.xlim(-1.2,1.2)
plt.ylim(-1.2,1.2)
plt.grid()
plt.show()

def unit_circle_indicator(x):
  return np.floor(np.linalg.norm(x,axis=0)<1)

K = 100
m = np.logspace(1,8,K, dtype='int')
estim = np.zeros(K)

for k in range(K):
  xs = np.random.rand(2,m[k])*2 - 1
  estim[k] = np.mean(unit_circle_indicator(xs))*4

plt.figure(figsize=(8,4), dpi=100)
plt.semilogx(m,estim, label='Estimate')
plt.semilogx(m,np.pi*np.ones(K), ':', linewidth=2, label='Accurate')
plt.grid()
plt.legend()
plt.xlabel('Number of samples')
plt.ylabel(r'$\pi$ estimate')
plt.show()

def unit_ball_indicator(x):
  return np.floor(np.linalg.norm(x,axis=0)<1)

K = 100
m = np.logspace(1,8,K, dtype='int')
estim = np.zeros(K)

for k in range(K):
  xs = np.random.rand(6,m[k])*2 - 1
  estim[k] = np.mean(unit_ball_indicator(xs))*2**6

plt.figure(figsize=(8,4), dpi=100)
plt.semilogx(m,estim)
plt.grid()
plt.xlabel('Number of samples')
plt.ylabel(r'Unit 6-ball volumn estimate')
plt.show()

m = 1000000
N = 17
n_list = np.arange(N, dtype='int')
estim = np.zeros(N)

for i in range(N):
  n = n_list[i]
  xs = np.random.rand(n,m)*2 - 1
  estim[i] = np.mean(unit_ball_indicator(xs))*(2**n)

plt.figure(figsize=(8,4), dpi=100)
plt.plot(n_list,estim)
plt.grid()
plt.xlabel(r'$n$')
plt.ylabel(r'Volume of unit $n$-balls')
plt.show()

from scipy.special import gamma
true_volume = np.pi**(n_list/2)/gamma(n_list/2+1)

plt.figure(figsize=(8,4), dpi=100)
plt.plot(n_list, estim, label='Estimate')
plt.plot(n_list, true_volume, ':', linewidth=2, label='Accurate')
plt.grid()
plt.legend()
plt.xlabel(r'$n$')
plt.ylabel(r'Volume of unit $n$-balls')
plt.show()

!pip install yfinance

import yfinance as yf
alpha = yf.download('GOOGL', start='2004-08-19')
price = alpha['Adj Close'].values

plt.figure(figsize=(10,3), dpi=100)
plt.plot(alpha.index, price)
plt.grid(True)
plt.xlabel('Days')
plt.ylabel('Adjusted close price')
plt.title('Alphabet Inc. (GOOGL)')
plt.show()

alpha['Returns'] = alpha['Adj Close'].pct_change()
v_annual = alpha['Returns'].std()*np.sqrt(252)
print (f'Annual Volatility: {v_annual*100:7.2f}%')

days = (alpha.index[-1]-alpha.index[0]).days
cagr = (price[-1]/price[1])**(365.0/days) - 1
print (f'Compound Annual Growth Rate: {cagr*100:7.2f}%')

T = 252

price_future = np.ones(T)*price[-1]

for t in range(T-1):
  ret = cagr/252 + np.random.randn()*v_annual/np.sqrt(252)
  price_future[t+1] = price_future[t] * (1+ret)

plt.figure(figsize=(10,3), dpi=100)
plt.plot(price, label='Price history')
plt.plot(np.arange(T)+len(price), price_future, label='Predicted price')
plt.grid(True)
plt.xlabel('Days')
plt.ylabel('Adjusted close price')
plt.title('Alphabet Inc. (GOOGL)')
plt.legend()
plt.show()

N = 10000

price_final = np.zeros(N)

plt.figure(figsize=(10,3), dpi=100)
plt.plot(price)

for n in range(N):
  price_future = np.ones(T)*price[-1]

  for t in range(T-1):
    ret = cagr/252 + np.random.randn()*v_annual/np.sqrt(252)
    price_future[t+1] = price_future[t] * (1+ret)
  price_final[n] = price_future[-1]

  plt.plot(np.arange(T)+len(price), price_future)

plt.grid(True)
plt.xlabel('Days')
plt.ylabel('Adjusted close price')
plt.title('Alphabet Inc. (GOOGL)')
plt.show()

Expected_value = np.mean(price_final)
Five_percentile = np.percentile(price_final,5)
Prob_double = np.sum(price_final>price[-1]*2)/N
print (f'Expected value: {Expected_value}')
print (f'5-percentile value: {Five_percentile}')
print (f'Prob_double: {Prob_double}')

plt.figure(figsize=(10,3), dpi=100)
plt.hist(price_final, bins=100)
plt.grid(True)
plt.xlabel('Evaluated price')
plt.ylabel('Frequency')
plt.show()

Expected_value = np.mean(price_final)
Five_percentile = np.percentile(price_final,5)
Prob_double = np.sum(price_final>price[-1]*2)/N
print (f'Expected value: {Expected_value}')
print (f'5-percentile value: {Five_percentile}')
print (f'Prob_double: {Prob_double}')

plt.figure(figsize=(10,3), dpi=100)
plt.hist(price_final, bins=100)
plt.axvline(Five_percentile, color='r', label='5-percentile')
plt.axvline(Expected_value, color='y', label='Mean')
plt.axvline(price[-1]*2, color='k', label='Double price')
plt.grid(True)
plt.legend()
plt.xlabel('Evaluated price')
plt.ylabel('Frequency')
plt.show()