# -*- coding: utf-8 -*-
"""
Created on Mon Apr 12 14:51:49 2021

@author: turinici
"""
import numpy as np
import matplotlib.pyplot as plt


#TODO implement a brownian motion

#first idea: use the properties of the Wt:    
T=1.0
N=255
M=10#number of scenarios
dt= T/N
W0=0#standard brownian motion
trange=np.linspace(0,T,N+1,endpoint=True)
# We know that Wt is a normal value of variance t
W1=np.sqrt(trange)*np.random.randn(N,1)
plt.figure(1)
plt.subplot(1,2,1)
plt.plot(trange[1:],W1)
#not working because the covariance is always zero ... and not min(s,t)
#good implementation: with increments

#another idea: use the cummulative increments property of B.M.
dW=np.sqrt(dt)*np.random.randn(N,M)
W=np.zeros((N+1,M)) 
W[0,:]=W0
W[1:,:]=W0+np.cumsum(dW,0)

plt.figure(2)
plt.plot(trange,W)

#compute \int_0^T W_t d W_t : using the Riemann-Ito sums
# in fact we compute sum_n W(t_n) * (increment between $t_n$ and $t_n+h$)
    # also compute integral minus W_T^2/2 and plot for all scenarios

int_WdW=np.zeros_like(W)
int_WdW[0,:]=0.0

for ii in range(N):
    int_WdW[ii+1,:] =int_WdW[ii,:]+ W[ii,:]*dW[ii,:]
    
plt.figure(3,figsize=(15,5))
plt.subplot(1,3,1)
plt.plot(trange,int_WdW)
plt.title('$t \mapsto \int_0^t W_u d W_u$')
plt.xlabel('t')
plt.subplot(1,3,2)
plt.plot(trange,W**2/2-int_WdW)
plt.title('$t \mapsto W_t^2/2-\int_0^t W_u d W_u$')
plt.xlabel('t')
plt.subplot(1,3,3)
plt.plot(trange,W**2/2-int_WdW- trange[:,None]**2/2)
plt.title('$t \mapsto W_t^2/2-t/2-\int_0^t W_u d W_u$')
plt.xlabel('t')
plt.tight_layout() 
plt.show()