In [1]:
from histogram1d import histogram1d

In [3]:
import matplotlib.pyplot as plt
%config InlineBackend.figure_format = 'retina'
%matplotlib inline
from numpy import *

In [9]:
def p(v): return 3/8*v**2
def g(u): return (8*u)**(1/3)

In [10]:
n = 1000000
u = random.rand(n)  # uniform random numbers
v = g(u)
# make histogram of v
leftends,heights = histogram1d(v,0,2,100)
w = leftends[1] - leftends[0] # bin width
plt.bar(leftends,heights,w,align='edge')

# add a plot of the density we're supposed to be getting:
v = linspace(0,2,100)
plt.plot(v,p(v),'r',alpha=0.5,lw=3);

In [14]:
def g(u,N,lam): return -log(1-u)/(N*lam)
def p(tau,N,lam): return N*lam*exp(-N*lam*tau)

n = 1000000
u = random.rand(n)  # uniform random numbers

N = 10000
taumax = 0.001
lam = 1
tau = g(u,N,lam)
# make histogram of tau
leftends,heights = histogram1d(tau,0,taumax,100)
w = leftends[1] - leftends[0] # bin width
plt.bar(leftends,heights,w,align='edge')

# add a plot of the density we're supposed to be getting:
tau = linspace(0,taumax,100)
plt.plot(tau,p(tau,N,lam),'r',alpha=0.5,lw=3);

In [20]:
## Standard normal distribution
from scipy.special import erfinv
def g(u): return sqrt(2)*erfinv(2*u-1)
def p(v): return 1/sqrt(2*pi)*exp(-v**2/2)

n = 1000000
u = random.rand(n)  # uniform random numbers
v = g(u)
# make histogram of v
leftends,heights = histogram1d(v,-3,3,100)
w = leftends[1] - leftends[0] # bin width
plt.bar(leftends,heights,w,align='edge')

# add a plot of the density we're supposed to be getting:
vv = linspace(-3,3,100)
plt.plot(vv,p(vv),'r',alpha=0.5,lw=3);

In [21]:
v.std()

Out[21]:
0.99998894703807051
In [23]:
from histogram1d import histogram1d
plt.figure(figsize=(16,8))
n = 100000
colors = 'rgbcymk'
kvals = [1,2,3,4,10,100]
isub = 0
for k,color in zip(kvals,colors):  # k is the number of random numbers we're adding together

u = random.rand(k,n)  #replace this with something non-uniform
def g(u): return (8*u)**(1/3)
u = g(u)
u -= u.mean()

v = u.sum(axis=0)/sqrt(k)
lefts, p = histogram1d(v,-2,2,99)
w = lefts[1]-lefts[0]
isub += 1
plt.subplot(2,6,isub)
plt.bar(lefts,p,w,align='edge',alpha=0.5,color=color)
sigma = v.std()
#print('sigma',sigma)
vv = linspace(-2,2,300)
pp = 1/sqrt(2*pi*sigma**2)*exp(-vv**2/2/sigma**2)
plt.plot(vv,pp,'k',alpha=0.75,lw=2)
plt.title(str(k))
plt.axis('off')
plt.xlim(-1,1);    #break
plt.suptitle('Shrink v by factor of k');


# Stock price simulation¶

In [27]:
n = 365
r = 0 #.03
L = empty(n+1)
L[0] = 0
for t in range(n):
L[t+1] = L[t] + r + sigma*random.randn()

plt.plot(L,alpha=0.2);

In [34]:
ntrials = 40
for i in range(ntrials):
n = 365
r = .02
L = empty(n+1)
L[0] = 0
for t in range(n):
L[t+1] = L[t] + r + sigma*random.randn()
plt.plot(L,'b',alpha=0.4);

In [41]:
sigma

Out[41]:
0.38659983160790179
In [ ]:
# Look at histogram of outcomes after 1 year

In [39]:
ntrials = 4000
finalLs = []
for i in range(ntrials):
n = 365
r = .02
L = empty(n+1)
L[0] = 0
for t in range(n):
L[t+1] = L[t] + r + sigma*random.randn()
#plt.plot(L,'b',alpha=0.4);
finalLs.append(L[-1])
#finalLs
finalLs = array(finalLs)
leftends,heights = histogram1d(finalLs,min(finalLs),max(finalLs),100)
w = leftends[1] - leftends[0] # bin width
plt.bar(leftends,heights,w,align='edge')
plt.xlabel('Log price')

Out[39]:
Text(0.5,0,'Log price')
In [45]:
sigma = .01
r = .02
ntrials = 4000
finalLs = []
for i in range(ntrials):
n = 365
L = empty(n+1)
L[0] = 0
for t in range(n):
L[t+1] = L[t] + r + sigma*random.randn()
#plt.plot(L,'b',alpha=0.4);
finalLs.append(L[-1])
#finalLs
finalLs = exp(array(finalLs))  # actual amounts leftends,heights = histogram1d(finalLs,min(finalLs),max(finalLs),100) w = leftends[1] - leftends[0] # bin width plt.bar(leftends,heights,w,align='edge') plt.xlabel('')

Out[45]:
Text(0.5,0,'\$')