Algorithm 4.1 30 step random walk 

import numpy as np
import matplotlib.pyplot as plt
nSteps = 30;
nTrials = 200000;
H = [0]*nTrials;
for j in range(nTrials):
  # do all nSteps coin tosses
  T = 2*(np.random.rand(nSteps) < 0.5 )-1;
  H[j] = sum(T);
plt.hist(H)
plt.show()

Algorithm 4.2 2d random walk

import numpy as np
import random
import matplotlib.pyplot as plt

sig=2; #step size std dev
nSteps = 2000; #steps each walk
nSteps = 2*int(nSteps/2); # assure even}

x = [0]*nSteps; y = [0]*nSteps;
for i in range(1,nSteps):
  #print(x[i-1],y[i-1])
  theta = 2*np.pi*np.random.random()
  N = np.random.normal(0,sig)
  x[i]= x[i-1]+N*np.cos(theta);
  y[i]= y[i-1]+N*np.sin(theta);

#typical walk
plt.subplot(2,2,1); #1st in 2x2 plots, upper left
plt.plot(x,y);
plt.title('(A) Typical walk for '+str(nSteps)+' steps');
plt.xlabel('X position');
plt.ylabel('Y position');

#do the walk from subStart to subEnd}
sub1Start = 500;
sub1End = 1500;
xsub1 = [0]*(sub1End-sub1Start)
ysub1 = [0]*(sub1End-sub1Start)
k = 0;
for i in range(sub1Start,sub1End):
  xsub1[k] = x[i];
  ysub1[k] = y[i];
  k = k+1;
plt.subplot(2,2,2); #2nd in 2x2 plots, upper right
plt.title("Walk (A) from step 500 to 1500");
plt.plot(xsub1,ysub1);

sub2Start = 1000;
sub2End = 1300;
xsub2 = [0]*(sub2End-sub2Start)
ysub2 = [0]*(sub2End-sub2Start)
k = 0;
for i in range(sub2Start,sub2End):
  xsub2[k] = x[i];
  ysub2[k] = y[i];
  k = k+1;
plt.subplot(2,2,3);
plt.title("Walk (A) from step 1000 to 1300");
plt.plot(xsub2,ysub2);

nWalks = 400; # walks for density
Endx = [0]*nWalks; Endy = [0]*nWalks;
for j in range(nWalks):
  x = [0]*nSteps; y = [0]*nSteps;
  for i in range(1,nSteps):
    theta = 2*np.pi*np.random.random()
    N = np.random.normal(0,sig)
    x[i]= x[i-1]+N*np.cos(theta);
    y[i]= y[i-1]+N*np.sin(theta);
  Endx[j] = x[-1];
  Endy[j] = y[-1];

#density plot
plt.subplot(2,2,4);
plt.scatter(Endx,Endy,s=3);
plt.title("Density Plot for "+str(nSteps)+" steps");
plt.xlabel('X end');
plt.ylabel('Y end');

plt.tight_layout();
plt.show()

Algorithm 4.3 Pricing a Call Option

import numpy as np
import matplotlib.pyplot as plt

nDays = 90;
dt = 1/365.0; # use time in years
T = nDays*dt; # expiry in years
S0=20; # starting stock price
X=25; # strike price
r = .031; # risk-free rate
sigma = .6; # volatility
expTerm = r*dt;
stddev = sigma*np.sqrt(dt);
nTrials = 10000;
S90 = [0]*nTrials;
value=0;
for j in range(nTrials):
  # make a vector of N(0,1) random numbers
  N = np.random.randn(nDays);
  S = S0;
  for i in range(nDays):
    dS = S*(expTerm + stddev*N[i]);
    S = S + dS;
  S90[j] = S; # to histogram
  value = value + max(S-X,0);
value = value/float(nTrials);  # expected value
price = np.exp(-r*T)*value # discount back
print("price= "+str(Price))
plt.hist(S90,bins=30)
plt.show()

Algorithm 4.4 Fortune under Kelly Bets

import numpy as np
import matplotlib.pyplot as plt
def kellyEval(nPlays,bf):
  F = 1 # starting fortune
# save the history
  h = [0]*nPlays; h[0] = F;
  for i in range(1,nPlays):
    w = F*bf  # wager}
    r = np.random.random()
    if( r < 0.6):
      F += w
     else:
      F -= w
    h[i] = F
  return h
nPlays = 300
bf = [0.05,0.1,0.2,0.25]; # betting fraction
for i in range(len(bf)):
  h = kellyEval(N,bf[i])
  eF = h[-1]; # ending fortune
  gain = np.exp((1/float(N))*np.log(h[-1]))}
  print("betting fraction: "+str(bf[i])+
" ending fortune= "+str(eF)+" gain= "+str(gain))
  plt.subplot(2,2,i+1)
  plt.plot(range(N),h)
plt.tight\_layout();
plt.show()

Algorithm 4.5 Calculating the Fundamental Matrix

import numpy as np
I = np.identity(5)
p = 0.25
Ph = np.zeros((5,5))
for i in range(5):
  if( i == 0):
    Ph[i][1] = p
  elif( i == 1):
   Ph[i][0] = Ph[i][2] = Ph[i][3] = p
  elif( i == 2):
   Ph[i][1] = Ph[i][4] = p
  elif( i == 3 ):
   Ph[i][1] = Ph[i][4] = p
  else:
   Ph[i][2] = Ph[i][3] = p
N = np.linalg.inv(I-Ph)
R = np.array([[2*p,p],[0,p],[0,2*p],[2*p,0],[2*p,0]])
B = np.matmul(N,R)
print(B)

Algorithm 4.6 Kinetic MC for a Lattice of "Grains"

import numpy as np
def setRates(stateVec):
  rates = [0]*nSites
  for i in range(nSites):
    if( stateVec[i] == 0 ):
      rates[i] = 1 #A to B
    else:
      rates[i] = 5 #B to A
  return rates
def makeBreakPts(rates):
  bPts = [0]*nSites
  bPts[0] = rates[0]
  for i in range(1,nSites):
    bPts[i] = bPts[i-1]+rates[i]
  return bPts
def rouletteSelect(rates):
  bPts = makeBreakPts(rates)
  U = np.random.uniform(0,bPts[-1])
  k=0;
  while( U >= bPts[k]):
    k += 1;
  #k is selected
  return k
#---main---
nSites = 8 # 0 for A, 1 for B
stateVec = [0]*nSites # initial state, all A's
rates = setRates(stateVec) # set the rates
R = sum(rates)
nIts = 64
avgFracB = 0 #average fraction of B sites
# track the accumulated number of B sites
accnBsites = 0
t = 0
for i in range(nIts):
  # time to next event
  deltT = (-1/R)*np.log(np.random.uniform())
  t += deltT
  # choose a site to flip
  k = rouletteSelect(rates)
  stateVec[k] = 1-stateVec[k] # flip site k
  # update rates
  if( stateVec[k] == 0 ):
    rates[k] = 1 # rates[k] used to be 5
    R = R-5+1
  else:
    rates[k] = 5
    R = R-1+5  # rates[k] used to be 1
  nBsites = sum(stateVec)
  accnBsites += nBsites
  avgFracB += nBsites/float(nSites)
print("ending time: ",t)
avgFracB /= nIts
print("avg.Frac. of B sites",avgFracB)
accnAsites = nSites*nIts - accnBsites
print("ratio A to B sites= ",
accnAsites/float(accnBsites))