Preliminaries

Off[General::"spell"]; Off[General::"spell1"];

AppendTo[$Path, "/Users/andrzej/Desktop/Finance"];

AppendTo[$Path, "/Users/andrzej/Desktop"];

RateOfInterest = 0.1;

<< "BlackScholes`"

<< ItosLemma`

SetOptions[ItoD, SuppressTime -> False];
{TimeSymbol, TimeIncrement, BrownianIncrement, 
            CorrelationSymbol} = {t, dt, dW, ρ};

dB = ItoMake[B[t], 0, 1];

<< Statistics`

NEC = Security[0.3, 0.1, 550, 0.1]

RateOfInterest = 0.1;

SetAttributes[setVolatility, HoldFirst];
SetAttributes[setDrift, HoldFirst];
SetAttributes[setPrice, HoldFirst];
SetAttributes[setDividend, HoldFirst];
setVolatility[x_, 
          v_] := x[[1]] = 
          v;
setDrift[x_, v_] := x[[2]] = 
          v;
setPrice[x_, v_] := x[[3]] = 
          v;
setDividend[x_, 
          v_] := x[[4]] = 
          v;
getVolatility[x_Security] := 
        x[[1]];
getDrift[x_Security] := 
        x[[2]];
getPrice[x_Security] := 
        x[[3]];
getDividend[x_Security] := 
        x[[4]];

Approximate Method

The stochastic differential equation of the Black-Scholes model is:

d S = 
      S (μ  dt + Z σ (Sqrt[dt]))

Using the Itoslemma package this takes the form:

dS = 
      ItoMake[S[t], S[t] μ, S[t] σ]

We can use this SDE to simulate the random walk performed by the price of the stock:

Simulate[x_Security, T_, n_ : 100] := 
      Module[{μ = getDrift[x], σ = getVolatility[x], dt = (T/n)}, 
        NestList[#1 Random[
                NormalDistribution[1 + μ dt, σ (Sqrt[dt])]] &, 
          getPrice[x], n]]

We use NestList to iterate the function #1 Random[
            NormalDistribution[1 + μ dt, σ (Sqrt[dt])]] &, which multiplies the argument by a random number of the form (1 + 
            1 + μ dt + σ (Sqrt[dt]) 
              Random[Normal[0, 1]])== Random[NormalDistribution[1 + μ dt, σ (Sqrt[dt])]].

Here are

values = Simulate[NEC, 1, 10]

We can plot these values:

ListPlot[values, 
        PlotJoined -> True];

Note that the points on the x-axis range from 0 to 10. If we would like then to range from 0 to T (in our case 1) in steps of (T/n) then we first make a list of points to be plotted on the x-axis:

timesteps = Range[0, 1, 0.1]

Next we use Transpose to make a list of pairs:

pairs = Transpose[{timesteps, values}]

We then plot the pairs:

ListPlot[pairs, 
        PlotJoined -> True];

Note the way Transpose works:

Transpose[{{a, b, c}, {x, y, z}}]

The function below creates the simulation, but without plotting the graph. The option DisplayFunction->Identity prevents plotting. Module is used to isolate the variables form previous code:

PlotSimulate[x_Security, T_, n_ : 100] := 
      Module[{μ = getDrift[x], σ = getVolatility[x], dt = (T/n)}, 
        ListPlot[
          Transpose[{Range[0, T, dt], 
              NestList[#1 Random[
                      NormalDistribution[1 + μ dt, σ (Sqrt[dt])]] &, 
                getPrice[x], n]}], PlotJoined -> True, 
          PlotStyle -> RGBColor[1, 0, 0], PlotRange -> All, 
          Background -> RGBColor[0, 1, 1], 
          DisplayFunction -> Identity]]

Now this will make the simulation but without a picture:

sim1 = PlotSimulate[NEC, 12, 10];

This will make a picture:

Show[sim1, 
        DisplayFunction -> $DisplayFunction];

Here we make 50 simulations using 100 time steps (the default in our program) and then display them:

sims = 
        Table[PlotSimulate[NEC, 12], {50}];

Show[sims, 
        DisplayFunction -> $DisplayFunction];

Exact Method

There is another approach to making a simulation. We first find the solution of the SDE (1). Let's start with the following stochastic process:

Y = Y0 (E^(ν t + σ B[t]));

ItoD[Y] /. Y -> HoldForm[Y]

This means that if we take ν = μ - (1/2) (σ^2), Y will be the solution to (1). In other words, the solution is

(Subscript[S, t])=(Subscript[S, 0]) (E^(t (μ - ((σ^2)/2)) + σ B(t)))

If we take (Subscript[S, t]) as the starting value and consider the time interval Δt we get the following iterative scheme:

(Subscript[S, t + Δt])=(Subscript[S, t]) (E^(Δt (μ - ((σ^2)/2)) + σ B(Δt)))

This gives us a second possible way to simulate the process:

PlotSimulate1[x_Security, T_, n_ : 100] := 
      Module[{μ = getDrift[x], σ = getVolatility[x], dt = (T/n)}, 
        ListPlot[
          Transpose[{Range[0, T, dt], 
              NestList[#1 (E^Random)[NormalDistribution[(μ - ((σ^2)/2)) dt, σ (Sqrt[dt])]] &, getPrice[x], n]}], PlotJoined -> True, 
          PlotStyle -> RGBColor[1, 0, 0], PlotRange -> All, 
          Background -> RGBColor[0, 1, 1], 
          DisplayFunction -> Identity]]

sims1 = 
        Table[PlotSimulate1[NEC, 12], {50}];

Show[sims1, 
        DisplayFunction -> $DisplayFunction];

Unequal intervals

So far all our simulations used equal times intervals. Sometimes it is useful to be able to use intervals of unequal length. Suppose suppose we wish to sample our process only at fixed times, (Subscript[t, 0]), (Subscript[t, 1]), (Subscript[t, 2]), (Subscript[t, 3]), where the intervals (Subscript[t, i ]) - (Subscript[t, i - 1]) are not necessarily equal.

sampletimes = {(Subscript[t, 0]), (Subscript[t, 1]), (Subscript[t, 2]), (Subscript[t, 3])};

We first create a list of time intervals. The easiest way to this with Mathematica is by using the function ListConvolve:

timeintervals = 
      ListConvolve[{1, -1}, sampletimes]

We replace the functions Nest and NestList by more general functions Fold and FoldList. This is how Fold works:

Fold[f, 
      x, {a, b, c, d}]

FoldList[f, x, {a, b, c, d}]

Note that f has to be a function of two variables. Here is anther function that creates a simulation based on a list of sample times, using the second approach based on the solution of the Black-Scholes SDE. As the first argument to FoldList we sue the pure function of two variables:

\(TraditionalForm\`\(\(
\)\(#1\ E\^HoldForm[Random[NormalDistribution[\((μ - σ\^2\/2)\)\ #2, σ\ \@#2]]] &\)\)\)

Note that HoldForm was included to prevent Mathematica attempting to evaluate  the exponent, which it would no t be able to do without numerical values for μ and σ. This is what this function does:

\(TraditionalForm\`\(#1\ E\^HoldForm[Random[NormalDistribution[\((μ - σ\^2\/2)\)\ #2, σ\ \@#2]]] &\)[u, v]\)

Here is what happens when we use FoldList:

\(TraditionalForm\`FoldList[#1\ E\^HoldForm[Random[NormalDistribution[\((μ - σ\^2\/2)\)\ #2, σ\ \@#2]]] &, t\_0, timeintervals]\)

\(PlotSimulate2[x_Security, sampletimes_] := 
      Module[{μ = getDrift[x], σ = getVolatility[x], 
          timeintervals = ListConvolve[{1, \(-1\)}, sampletimes]}, 
        ListPlot[
          Transpose[{sampletimes, 
              FoldList[#1\ E\^Random[NormalDistribution[\((μ - σ\^2\/2)\)\ #2, σ\ \@#2]] &, getPrice[x], timeintervals]}], PlotJoined -> True, 
          PlotStyle -> RGBColor[1, 0, 0], 
          Background -> RGBColor[0, 1, 1], 
          DisplayFunction -> Identity]]\)

sampletimes = Union[Table[Random[], {20}]]

s = 
        Table[PlotSimulate2[NEC, sampletimes], {50}];

Show[s, 
        DisplayFunction -> $DisplayFunction];

Monte Carlo Method

The above method can be used not only for making graphics of stock movements but also for actual computations, using the so called Monte-Carlo method. Let's illustrate it with the example of a European call option. The basic idea is this. Recall the principle of "risk neutral" valuation for options. This says that the value of an option can be computed as the discounted future payoff of the option assuming that the stock performs a random walk given by the Balck-Scholes SDE, with the drift replaced by the risk free interest rate ( -the  continuous dividend rate). Using our first method of simulation this takes the following form:

\(TraditionalForm\`monteCarlo[europeanCall[time_, strike_, security_], 
        timesteps_, iterations_] := 
      Module[{dt = time\/timesteps, r = RateOfInterest, 
          div = getDividend[security], σ = getVolatility[security], 
          S = getPrice[security], }, 
\(E\^\(\(-r\)\ time\)\) 
          Mean[Table[
              Max[0, Nest[#1\ E\^Random[NormalDistribution[\((r - div - σ\^2\/2)\)\ dt, σ\ \@dt]] &, S, timesteps] - 
                  strike], {iterations}]]]\)

Note that the function Mean is a built in function only in Mathematica  version 5.  In previous versions it is defined in one of the statistics packages which we loaded in at the beginning of this session. Note that we use Nest and not NestList, since in this case we only need to know the value of the stock at the strike time.

Here is our security:

NEC

RateOfInterest = 0.1;

We shall compute this using 1000 simulated paths each with 100 steps.

monteCarlo[europeanCall[3, 550, NEC], 100, 
      1000]

Recall from the Introduction that the Black-Scholes formula gives us the

value[europeanCall[3, 550, NEC], Method -> analytic[]]