Page 166

Matlab: Coin Toss Random Walk

nSteps = 30;
nTrials = 200000;
H=zeros(1,nTrials);
for j=1:nTrials
  % do all nSteps coin tosses
  T = 2*(rand(1,nSteps) < 0.5 ) -1;
  H(j) = sum(T);
end
hist(H)

Page 172

Matlab: 2-D Brownian Motion

sig=2; % step size std dev
nSteps = 2000; % # steps in each walk
nSteps = 2*fix(nSteps/2); % assure even
% box-muller for N(0,1) random variables
u1=rand(nSteps/2,1); %U[0,1) vector nSteps/2 long
u2=rand(nSteps/2,1);
x1 = cos(2*pi*u1).*sqrt(-2*log(u2)); %N(0,1) vector
x2 = sin(2*pi*u1).*sqrt(-2*log(u2));
% concatenate x1 and x2 into h
h(1:nSteps/2) = x1;
h(1+nSteps/2:nSteps) = x2;
x(1)=0; y(1)=0; %start at the origin
for i=2:nSteps
  theta = 2*pi*rand; % random direction
  x(i)=x(i-1)+sig*h(i)*cos(theta);
  y(i)=y(i-1)+sig*h(i)*sin(theta);
end
%typical walk
figure(1)
subplot(2,2,1);
plot(x,y);
title('Typical walk for nSteps steps');
xlabel('X position');
ylabel('Y position');
% do the walk from subStart to subEnd
sub1Start = 500;
sub1End = 1500;
k = 1;
for i=sub1Start:sub1End
  xsub1(k) = x(i);
  ysub1(k) = y(i);
  k = k+1;
end
figure(1)
subplot(2,2,2);
plot(xsub1,ysub1);
sub2Start = 1000;
sub2End = 1300;
k = 1;
for i=sub2Start:sub2End
  xsub2(k) = x(i);
  ysub2(k) = y(i);
  k = k+1;
end
figure(1)
subplot(2,2,3);
plot(xsub2,ysub2);
nWalks = 400; % # walks for density
for j=1:nWalks
  u1=rand(nSteps/2,1);
  u2=rand(nSteps/2,1);
  x1 = cos(2*pi*u1).*sqrt(-2*log(u2));
  x2 = sin(2*pi*u1).*sqrt(-2*log(u2));
  h(1:nSteps/2) = x1;
  h(1+nSteps/2:nSteps) = x2;
  x(1)=0; y(1)=0;
  for i=2:nSteps
    theta = 2*pi*rand;
    x(i)=x(i-1)+sig*h(i)*cos(theta);
    y(i)=y(i-1)+sig*h(i)*sin(theta);
  end
  Endx(j) = x(nSteps);
  Endy(j) = y(nSteps);
end
%density plot
figure(1);
subplot(2,2,4);
scatter(Endx,Endy,3);
title('Density Plot for nSteps steps');
xlabel('X end');
ylabel('Y end');

Page 178

Matlab: Option Price via Simulation

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*sqrt(dt);
nTrials = 10000;
value=0;
for j=1:nTrials
  n = randn(1,nDays); % vector of N(0,1) random numbers
  S = S0;
  for i=1:nDays
    dS = S*(expTerm + stddev*n(i));
    S = S + dS;
  end
  S90(j) = S; % to histogram
  value = value + max(S-X,0);
end
value = value/nTrials;  % expected value
Price = exp(-r*T)*value % discount back
hist(S90,0:.5:65)

Page 188

Matlab: Kelly Problem Evaluation via Simulation

function fitnesses = kellyEval(popSize,chromPop);
popSize = 1; % remove for use in GA
chromPop=0.05; % remove for use in GA
F=linspace(1,1,popSize); % initial fortune
N=1000; % iterations of the game
h=linspace(1,1,N); % histogram, remove for GA 
for i=2:N
  w = F.*chromPop; %term-by-term wager
    % play 60/40 game
  r = (rand(1,popSize) < 0.6); % 1 if <0.6 else 0
  r = 2*r-1; % +1 if <0.6 else -1
  w = r.*w; % +bet if <0.6, else -bet
  h(i)=F; % remove for use in GA
  F = F+w; % new fortune
end
fitnesses = exp((1/N)*log(F));

Page 195

Matlab: Kelly Problem Solution via Genetic Algorithms

% a "chromosome" here is a fraction between 0 and 1
% a population is a vector of chromosomes
popSize=32;
chromPop= rand(1,popSize); %initalize randomly
 % space for their fitnesses = exp(growth rate)
fitnesses=linspace(0,0,popSize);
newPop=linspace(0,0,popSize);
upper = linspace(1,1,popSize); % keep fractions < 1
lower = linspace(0,0,popSize); % keep fractions > 0
bp = linspace(0,0,popSize); % roulette wheel breakpts
pm = 0.1; % mutation probability
orelax = 1.4; % over relaxation factor for mating
epsilon=0.07 % mutation range
bestFitness=0;
%chromPop % for debugging
for g=1:100  % 100 generations
   % mutate each w/prob pm
   m = (rand(1,popSize)<pm); % mutation indices
    % mutate each chromoPop(i) within epsilon
   a = max(lower,chromPop-epsilon);
   b = min(upper,chromPop+epsilon);
   r = m.*(a+(b-a).*rand(1,popSize)); % new value
   chromPop = (1-m).*chromPop; %null out old value
   chromPop = chromPop+r; %replace by new value
  %chromPop % for debugging
  % mate each chromosome via a permutation of chromPop
   permPop = chromPop;
   for i=1:popSize-1
     k = i+floor((popSize+1-i)*rand); %random in remaining
     swap = permPop(i);
     permPop(i) = permPop(k);
     permPop(k)=swap;
  end
%permPop % for debugging
  newPop = chromPop +...
   orelax*(permPop-chromPop).*rand(1,popSize);
    % stay within 0 and 1
  newPop = max(lower,newPop);
  newPop = min(upper,newPop);
%newPop % for debugging
    % evaluate newPop
  fitnesses = kellyEval(popSize,newPop);
  oldBest = bestFitness;
  bestFitness = max(oldBest,max(fitnesses));
  if bestFitness > oldBest
   m = fitnesses>=bestFitness; % where is it
   solution = m.*newPop;
  end
% roulette wheel select, make up the wheel
  sum = 0;
  for i=1:popSize
   sum = sum + fitnesses(i);
   bp(i) = sum;
  end
%bp % for debugging
    % now pick popSize times
  for i=1:popSize
   U=rand; k=1;
   while( U>= bp(k))
     k=k+1;
   end; %k is selected
   chromPop(i) = newPop(k);
  end
%chromPop % for debugging
end
bestFitness
solution