(* setting the parameters for the forward binomial tree *)
(* initial stock price *) S0=100;(* volatility *)σ=0.3; (* dividend yield *) δ=0; (* continuously compounded interest rate *)r=0.04; (* length of each period *) h=1; (* number of periods *) n=1;

T=n*h

(* up factor *) u = E^((r-δ)*h+σ*Sqrt[h])

(* down factor *) d = E^((r-δ)*h-σ*Sqrt[h])

(* risk-neutral probability *) pstar=1/(1+E^(σ*Sqrt[h]))

(* strike price *) K=100

(* drawing Bernoulli trials *) LBern=RandomVariate[BernoulliDistribution[pstar],n]

(* cointosses from Bernoulli trials *) LCoinToss=-1+2*LBern

(* prepending the origin *) LWalk=Prepend[Accumulate[LCoinToss],0]

(* plotting the walk *) WalkToPlot[w_]:=Table[{i-1,w[[i]]}, {i,1,Length[w]}]; WalkToPlot[LWalk]; ListLinePlot[WalkToPlot[LWalk]]

(* height of the final position *) ht=LWalk[[n+1]]

(* number of steps up *) upsteps= (ht+n)/2

(* number of steps down *) downsteps= n-upsteps

(* from cointosses to stock prices *) S1=S0* u^upsteps*d^downsteps

(* from stock prices to the payoffs *) V1=Max[S1-K,0]

(* simulated price *) VSim= E^(-r*T)*V1

(* analytic risk-neutral price *) VAn=E^(-r*T)*pstar*Max[S0*u-K,0]

(* this is TERRIBLE -- we need to have MORE than one simulated final stock price !!! *)

(* number of simulations *) nsim=2

(* repeated Bernoulli trials *)LBernSims=Table[RandomVariate[BernoulliDistribution[pstar],n],{nsim}]

(* cointosses from Bernoulli trials *) LCoinTossSims=-1+2*LBernSims

(* prepending the origin *) LWalkSims=Table[ Prepend[Accumulate[LCoinTossSims[[i]]],0],{i,1,nsim}]

(* formatted to plot *) PlotFormat=Table[WalkToPlot[LWalkSims[[i]]],{i,1,nsim}]

(* plotting *) ListLinePlot[PlotFormat]

(* height of the final position *) htSims=Table[LWalkSims[[i,n+1]], {i,1,nsim}]

(* number of steps up *) upstepsSims= (htSims+n)/2

(* number of steps down *) downstepsSims= n-upstepsSims

(* from cointosses to stock prices *) S1Sims=Table[S0*u^upstepsSims[[i]]*d^downstepsSims[[i]], {i,1,nsim}]

(* from stock prices to the payoffs *) V1Sims=Table[Max[S1Sims[[i]]-K,0],{i,1,nsim}]

(* simulated price *) VSim= E^(-r*T)*Total[V1Sims]/nsim

(* we need a LOT of simulations !!!! *)

(* number of simulations *) nsim=200

(* repeated Bernoulli trials *)LBernSims=Table[RandomVariate[BernoulliDistribution[pstar],n],{nsim}];

(* cointosses from Bernoulli trials *) LCoinTossSims=-1+2*LBernSims;

(* prepending the origin *) LWalkSims=Table[ Prepend[Accumulate[LCoinTossSims[[i]]],0],{i,1,nsim}];

(* formatted to plot *) PlotFormat=Table[WalkToPlot[LWalkSims[[i]]],{i,1,nsim}];

(* plotting *) ListLinePlot[PlotFormat]

(* height of the final position *) htSims=Table[LWalkSims[[i,n+1]], {i,1,nsim}];

(* number of steps up *) upstepsSims= (htSims+n)/2;

(* number of steps down *) downstepsSims= n-upstepsSims;

(* from cointosses to stock prices *) S1Sims=Table[S0*u^upstepsSims[[i]]*d^downstepsSims[[i]], {i,1,nsim}] ;

(* from stock prices to the payoffs *) V1Sims=Table[Max[S1Sims[[i]]-K,0],{i,1,nsim}];

(* simulated price *) VSim= E^(-r*T)*Total[V1Sims]/nsim

(* MANY MANY more simulations ..... *)

(* number of simulations *) nsim=10000

(* repeated Bernoulli trials *)LBernSims=Table[RandomVariate[BernoulliDistribution[pstar],n],{nsim}];

(* cointosses from Bernoulli trials *) LCoinTossSims=-1+2*LBernSims;

(* prepending the origin *) LWalkSims=Table[ Prepend[Accumulate[LCoinTossSims[[i]]],0],{i,1,nsim}];

(* formatted to plot *) PlotFormat=Table[WalkToPlot[LWalkSims[[i]]],{i,1,nsim}];

(* plotting *) ListLinePlot[PlotFormat]

(* height of the final position *) htSims=Table[LWalkSims[[i,n+1]], {i,1,nsim}];

(* number of steps up *) upstepsSims= (htSims+n)/2;

(* number of steps down *) downstepsSims= n-upstepsSims;

(* from cointosses to stock prices *) S1Sims=Table[S0*u^upstepsSims[[i]]*d^downstepsSims[[i]], {i,1,nsim}] ;

(* from stock prices to the payoffs *) V1Sims=Table[Max[S1Sims[[i]]-K,0],{i,1,nsim}];

(* simulated price *) VSim= E^(-r*T)*Total[V1Sims]/nsim

(* how about more time periods? *) n=5

(* exercise date *) T=n*h

(* number of simulations *) nsim=1000

(* repeated Bernoulli trials *)LBernSims=Table[RandomVariate[BernoulliDistribution[pstar],n],{nsim}];

(* cointosses from Bernoulli trials *) LCoinTossSims=-1+2*LBernSims;

(* prepending the origin *) LWalkSims=Table[ Prepend[Accumulate[LCoinTossSims[[i]]],0],{i,1,nsim}];

WalkToPlot[w_]:=Table[{i-1,w[[i]]}, {i,1,Length[w]}]; 
(* formatted to plot *) PlotFormat=Table[WalkToPlot[LWalkSims[[i]]],{i,1,nsim}];

(* plotting *) ListLinePlot[PlotFormat]

(* height of the final position *) htSims=Table[LWalkSims[[i,n+1]], {i,1,nsim}];

(* number of steps up *) upstepsSims= (htSims+n)/2;

(* number of steps down *) downstepsSims= n-upstepsSims;

(* from cointosses to stock prices *) S1Sims=Table[S0*u^upstepsSims[[i]]*d^downstepsSims[[i]], {i,1,nsim}] ;

(* from stock prices to the payoffs *) V1Sims=Table[Max[S1Sims[[i]]-K,0],{i,1,nsim}];

(* simulated price *) VSim= E^(-r*T)*Total[V1Sims]/nsim

(* analytic risk-neutral price *)VAn=E^(-r*T)*Expectation[Max[S0*u^x*d^(n-x)-K,0],x \[Distributed]BinomialDistribution[n,pstar]]