I like to show you some techniques for statistical evaluations in EMT. I will add a tutorial about this later to the files delivered with EMT.

We compute the value of a European call option. This is simply an option (but not an obligation) to buy an asset at a future time T for a price K. The value of this option is clearly

\(G = \max(0,S(T)-K)\)

if S(T) denotes the price of the asset at time T. This price is not predictable. So we assume a statistical model. For simplicity, we assume the somewhat dubious model of Black and Scholes which predicts that the logarithmic returns are distributed according to a normal distribution plus a trend. We could expand the model to more realistic distributions like some Pareto distribution or a normal distribution with random jumps at random time.

According to the B-S model,

\(S(T) = S_0 \exp((r-\sigma^2/2)T+\sigma\sqrt{T} X), \qquad X \sim N(0,1)\)

With this, we can compute the expected value of a specific option in EMT as follows.

>function C(X,S0,K,r,T,sigma) ...
$  S=S0*exp((r-sigma^2/2)*T+sigma*sqrt(T)*X);
$  return exp(-r*T)*max(0,S-K);
$endfunction
>integrate("C(x,100,110,0.02,10,0.02)*qnormal(x)",-10,10)
 10.0618725582

We have simply integrated the value G (which depends on X) with the normal density. Note that we had to discount the value with the assumed interest rate r.

It is one of the virtues of EMT and similar programs that we can check the result with a numerical simulation.

>mean(C(normal(1000000),100,110,0.02,10,0.02))
 10.0641102799

To get an impression of the distribution, we can plot the simulated values.

>plot2d(C(normal(1000000),100,110,0.02,10,0.02),>distribution):

Let us define a function in EMT for the value.

>function %Cf(x,S0,K,r,T,sigma) := C(x,S0,K,r,T,sigma)*qnormal(x);
>function map Cexact(S0,K,r,T,sigma) ...
$return integrate("%Cf",-10,10;S0,K,r,T,sigma);
$endfunction
>Cexact(100,110,0.02,10,0.02)
 10.0618725582

Using a simulation, we can also easily answer the question, how many values are larger than a given value, or we can compute the standard deviation of the price.

For the simple B-S model, these things can be answered exactly. Here is a formula for the expected value.

>function Cformula(S0,K,r,T,sigma) ...
$  d=log(S0/(K*exp(-r*T)))/(sigma*sqrt(T));
$  dpm=sigma*sqrt(T)/2;
$  return S0*normaldis(d+dpm)-K*exp(-r*T)*normaldis(d-dpm);
$endfunction
>Cformula(100,110,0.02,10,0.02)
 10.0618725582

What changes would have to be done for other statistical models?

For models with a density, we can simply replace qnormal() with the density q(), and define a random generator e.g. using the rejection method. For distributions with discrete jumps at specific times, we have to add a sum. For random jumps, things get a bit more involved. I will go into that in a later blog.