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,'$')