# -*- coding: utf-8 -*-
"""
Created on Fri Apr 22 14:16:24 2022

@author: turinici
"""

import numpy as np
import matplotlib.pyplot as plt

T=1
N=500
h=T/N
M=200
trange=np.linspace(0,T,N)
dW= np.random.randn(N)* np.sqrt(h)
W= np.cumsum(dW)
plt.plot(trange,dW,trange,W)

W = np.zeros((N+1,M))
dW= np.random.randn(N,M)* np.sqrt(h)
W[0,:]=0
W[1:,:]=W[0,:]+ np.cumsum(dW,axis=0)
plt.figure('Brownian')
plt.plot(W)

WT=W[-1,:]
intWtdW = WT**2/2-T/2
sri = np.sum(dW*W[0:-1,:],axis=0)

plt.figure("integral v2")
plt.hist(intWtdW,alpha=0.5,label="int")
plt.hist(sri,alpha=0.5,label=" discrete int")
plt.hist(sri-intWtdW,alpha=0.5,label='diff')

plt.figure("integral v3")
plt.hist([intWtdW,sri,sri-intWtdW],density=True)
plt.legend(['exact',"numeric","diff"])

Y=np.zeros_like(W)
theta=1.
mu=10.
sigma=3
for ii in range(N):
    Y[ii+1,:]=Y[ii,:]+ theta*(mu-Y[ii,:])*h+ sigma*dW[ii,:]
    
plt.figure("OU plot")
plt.hist([Y[0,:], Y[ int(N/4),:], 
          Y[int(N/2),:],Y[-1,:]  ],density=True)
plt.legend(['t=0',"t=T/4","t=T/2","t=T"])

plt.plot(Y)