We can create a file called asian.jl by using the following code:

function run_asian(N = 100000, PutCall = 'C';)
# European Asian option.
# Uses geometric or arithmetic average.
# Euler and Milstein discretization for Black-Scholes.
# Option features.
  println("Setting option parameters");
  S0  = 100;      # Spot price
  K   = 100;      # Strike price
  r   = 0.05;     # Risk free rate
  q   = 0.0;      # Dividend yield
  v   = 0.2;      # Volatility
  tma = 0.25;     # Time to maturity

  Averaging = 'A';  # 'A'rithmetic or 'G'eometric
  OptType = (PutCall == 'C' ? "CALL" : "PUT");
  println("Option type is $OptType");
# Simulation settings.
  println("Setting simulation parameters");
  T = 100;         # Number of time steps
  dt = tma/T;      # Time increment

# Initialize the terminal stock price matrices
# for the Euler and Milstein discretization schemes.
S = zeros(Float64,N,T);
  for n=1:N
    S[n,1] = S0;
  end

# Simulate the stock price under the Euler and Milstein schemes.
# Take average of terminal stock price.
  println("Looping $N times.");
  A = zeros(Float64,N);
  for n=1:N
    for t=2:T
      dW = (randn(1)[1])*sqrt(dt);
      z0 = (r - q - 0.5*v*v)*S[n,t-1]*dt;
      z1 = v*S[n,t-1]*dW;
      z2 = 0.5*v*v*S[n,t-1]*dW*dW;
      S[n,t] = S[n,t-1] + z0 + z1 + z2;
     end
    if cmp(Averaging,'A') == 0
      A[n] = mean(S[n,:]);
    elseif cmp(Averaging,'G') == 0
      A[n] = exp(mean(log(S[n,:])));
    end
  end

# Define the payoff
  P = zeros(Float64,N);
  if cmp(PutCall,'C') == 0
    for n = 1:N
      P[n] = max(A[n] - K, 0);
    end
  elseif cmp(PutCall,'P') == 0
    for n = 1:N
      P[n] = max(K - A[n], 0);
    end
  end
# Calculate the price of the Asian option
AsianPrice = exp(-r*tma)*mean(P);
@printf "Price: %10.4f\n" AsianPrice;
end

We have wrapped the main body of the code in a function run_asian(N, PutCall) .... end statement. The reason for this is to be able to time the execution of the task in Julia, thereby eliminating the startup times associated with the Julia runtime when using the console.

The stochastic behavior of the stock is modeled by the randn function; randn(N) provides an array of N elements, normally distributed with zero mean and unit variance.

All the work is done in the inner loop; the z-variables are just written to decompose the calculation.

To store the averages for each track, use the zeros function to allocate and initialise the array.

The option would only be exercised if the average value of the stock is above the "agreed" prices. This is called the payoff and is stored for each run in the array P.

It is possible to use arithmetic or geometric averaging. The code sets this as arithmetic, but it could be parameterized.

The final price is set by applying the mean function to the P array. This is an example of vectorized coding.

So, to run this simulation, start the Julia console and load the script as follows:

julia> include("asian.jl")
julia> run_asian()

To get an estimate of the time taken to execute this command, we can use the tic()/toc() function or the @elapsed macro:

include("asian.jl")
tic(); run_asian(1000000, 'C'); toc();
Option Price: 1.6788 elapsed: time 1.277005471 seconds

If we are not interested in the execution times, there are a couple of ways in which we can proceed.

The first is just to append to the code a single line calling the function as follows:

run_asian(1000000, 'C')

Then, we can run the Asian option from the command prompt by simply typing the following: julia asian.jl.

This is pretty inflexible since we would like to pass different values of the number of trials N and to determine the price for either a call option or a put option.

Julia provides an ARG array when a script is started from the command line to hold the passed arguments. So, we add the following code at the end of asian.jl:

nArgs = length(ARGS)
if nArgs >= 2
    run_asian(ARGS[1], ARGS[2])
elseif nArgs == 1
    run_asian(ARGS[1])
else
    run_asian()
end

Julia variables are case-sensitive, so we must use ARGS (uppercase) to pass the arguments.

Because we have specified the default values in the function definition, this will run from the command line or if loaded into the console.

Arguments to Julia are passed as strings but will be converted automatically although we are not doing any checking on what is passed for the number of trials (N) or the PutCall option.

We will discuss the package system in the next section, but to get going, let's look at a simple example that produces a plot by using the ASCIIPlots package. This is not dependent on any other package and works at the console prompt.

You can find it at http://docs.julialang.org and choose the "Available Packages" link. Then, click on ASCIIPlots, which will take you to the GitHub repository.

It is always a good idea when looking at a package to read the markdown file: README.md for examples and guidance.

Let's use the same parameters as before. We will be doing a single walk so there is no need for an outer loop or for accumulating the price estimates and averaging them to produce an option price.

By compacting the inner loop, we can write this as follows:

using ASCIIPlots;

S0  = 100;      # Spot price
K   = 102;      # Strike price
r   = 0.05;     # Risk free rate
q   = 0.0;      # Dividend yield
v   = 0.2;      # Volatility
tma = 0.25;     # Time to maturity
T = 100;        # Number of time steps
dt = tma/T;     # Time increment
S = zeros(Float64,T);
S[1] = S0;
dW = randn(T)*sqrt(dt)
[ S[t] = S[t-1] * (1 + (r - q - 0.5*v*v)*dt + v*dW[t] + 0.5*v*v*dW[t]*dW[t]) for t=2:T ]
x = linspace(1,T);
scatterplot(x,S,sym='*');

Note that when adding a package the statement using its name as an ASCIIPlots string, whereas when using the package, it does not.

using Winston;

S0  = 100;      # Spot price
K   = 102;      # Strike price
r   = 0.05;     # Risk free rate
q   = 0.0;      # Dividend yield
v   = 0.2;      # Volatility
tma = 0.25;     # Time to maturity
T = 100;        # Number of time steps
dt = tma/T;     # Time increment
S = zeros(Float64,T)
S[1] = S0;
dW = randn(T)*sqrt(dt);
[ S[t] = S[t-1] * (1 + (r - q - 0.5*v*v)*dt + v*dW[t] + 0.5*v*v*dW[t]*dW[t]) for t=2:T ]

x = linspace(1, T, length(T));
p = FramedPlot(title = "Random Walk, drift 5%, volatility 2%")
add(p, Curve(x,S,color="red"))
display(p)

We compared the Asian option code above with similar implementations in the "usual" data science languages discussed earlier.

The point of these benchmarks is to compare the performance of specific algorithms for each language implementation. The code used is available to download.

The runs were executed on a Samsung RV711 laptop with an i5 processor and 4GB RAM running CentOS 6.5 (Final).