If you want to choose randomly your next holidays destination, you are likely to process in a way which is certainly biased. Especially if you choose randomly the latitude and the longitude. A bit like they do in this lovely advertising (For those of you who do not speak French, this is about a couple who have won the national gamble prize and have to decide their next travel. The husband randomly picks Australia and the wife is complaining : "Not again!"). So let me help you to choose your next travel!


If we were able to generate uniformly distributed variables on [0 ; 1], we could easily generate variables on other spaces such as [0 ; 1] x [0 ; 1], which is a square of side 1. It can be simply done by using two independant variables X1 and X2, uniformly distributed on [0 ; 1], and by considering X = (X1 , X2). Then X is a random variable evenly distributed on the square [0 ; 1] x [0 ; 1]

However, this generation may be a bit more complicated if we work with more complicated shapes, such as a sphere for example. Indeed, for a sphere of a given radius, say R = 1, it is possible to project the variables of a square on the sphere.

The first method I thought of, is wrong. I was tempted to generate two variables Y1 and Y2 where Y1 follows a uniform on  [0 ; Pi] and Y2 follows a uniform on [-Pi ; Pi]. In other words :
- Y1 = Pi * X1
- Y2 = 2Pi * (X2 - 0.5).
 Then I wanted to consider (Y1, Y2) as the spherical coordinated of the sphere.

However, this method does NOT generate a uniform on the sphere. Indeed, it has a tendency to over generate north and south poles. The reason is simple, this method generates, in average the same amount of variable in each latitude. North and South poles are of smaller area than Equator. Therefore, the closer we are from the Equator latitude the less variables there are.

On the following graph, we observe a high density of point in the "middle" of the sphere, this is the North Pole. It shows that this method does not offer a uniform distribution. By the way, the graph is computed with the library rgl, which provides a display device for 3 dimensions objects. And then a little code provided in the rgl help allows to move automatically the device and take a snapshot for each view. You can eventually generate a GIF on gifmake (like for the map in Spatial segregation in cities - An explanation by a neural network model).
A wrong method to generate evenly distributed variables on the sphere. The North Pole and the South Pole are over represented.

Many methods have been proposed to generate evenly distributed random variables on a sphere. We propose one of them here. We consider the couple z = (u,v) defined as :

- u = 2 * Pi * X1
- v = arc-cos(2 * (X2 - 0.5))

In this case, a theorem shows that z = (u , v) generates evenly distributed variables. It can be observed on the next graph. There is no irregularity in the distribution of the random variables.
A correct simulation of a uniform distribution on a sphere. There is no over represented area.

Other methods exist. I like this one since it is really simple and uses a uniform distribution at the beginning. The idea of this post is to show that generating a uniform distribution can be adapted to many shapes and cases. However, to do so, a previous analytical study has to be done to find the correct transformation. 

So this method would help you to avoid being too many times in the chilly places such as North Pole and South Pole since it does not overrepresent the extreme latitudes.

The program (R) :

# import package to plot in 3D
install.packages("rgl", dependencies = TRUE)
library(rgl)

################################################################
# Uniform distribution in a square
################################################################

size = 10000
x1 = runif(size)
x2 = runif(size)

# the option pch = '.' change the symbol for the graph into dot.
# cex = 2 doubles the size of the dots
plot(x1,x2, col = 'blue', pch = '.', cex = 2)

################################################################
# Wrong solution for the sphere
################################################################

y1 = pi * x1
y2 = 2* pi * (x2-0.5)
y = matrix(0, nrow = 2, ncol = size)
y[1,] = y1
y[2,] = y2
plot(y1, y2)

# This function transform the spherical coordinates into cartesian coordinates
sphereToCartesian = function(matrice){
  x= matrix(0,nrow = 3, ncol = length(matrice[1,]))
  x[1,] = sin(matrice[2,]) * cos(matrice[1,])
  x[2,] = sin(matrice[2,]) * sin(matrice[1,])
  x[3,] = cos(matrice[2,])
  return(x)
}

a = sphereToCartesian(y)
plot3d(x = a[1,], y = a[2,], z = a[3,])

#you should enlarge the device window, before running this, if you want to have a meaningful graph
rgl.bringtotop()
rgl.viewpoint(0,20)

for (i in 1:45) {
  rgl.viewpoint(i,20)
  filename <- paste("pic",i,".png",sep="")
  rgl.snapshot(filename, fmt="png")
}


################################################################
# Correct solution for the sphere
################################################################

uniformSphere = function(length){
  x1 = runif(length, 0,1)
  x2 = runif(length, 0,1)
  u = 2*pi*x1
  v = acos(2*x2- 1)
  z =matrix(0, ncol = length, nrow = 2)
  z[1,] = u
  z[2,] = v
  return(z)
}

z = uniformSphere(size)
b = sphereToCartesian(z)
plot3d(x = b[1,], y = b[2,], z = b[3,])


4

View comments

The financial market is not only made of stock options. Other financial products enable market actors to target specific aims. For example, an oil buyer like a flight company may want to cover the risk of increase in the price of oil. In this case it is possible to buy on the financial market what is known as a "Call" or a "Call Option".

A Call Option is a contract between two counterparties (the flight company and a financial actor). The buyer of the Call has the opportunity but not the obligation to buy a certain  quantity of a certain product (called the underlying) at a certain date (the maturity) for a certain price (the strike).

If the flight company doesn't want to suffer from the risk of an increase of oil on the market within the next year, a Call Option of maturity 1 year, with a certain strike K may be bought. In this case, if the oil price on the market is over the strike, the flight company will use the Call option to buy oil for the price of the strike. If the price of the oil is lower than the strike, then the flight will not use the option and will buy oil on the market.

As you can see, this contract is not symmetric  The flight company has an option (to buy or not to buy) while the second counterparty just follows the decision of the flight company. Therefore, the flight company will have to pay something to the second counterparty. One of the great question is how much should it be charged.

We will therefore use a pricing algorithm to estimate the price of this Call Option. We first build an algorithm to simulate the value of an asset in the model of Black-Scholes. Then we simulate the historical value of this asset in order to simulate the final value of the Call Option.

Indeed, if we know the value of the asset at the maturity (A(T)), we get directly the value of the Call Option of strike K at maturity : max(A(T) - K , 0).

Our asset is randomly distributed as:
A(t+dt) = A(t) (1 + r dt + v(B(t+dt) -B(t))

Where r is the interest rate, v the volatility of the asset and B a Brownian process (for more detail you can look at the post on Brownian motion).

After repeating the process many times, we estimate the mean of the Call Option for different strike in order to estimate the price of the Call Option.


As we can see the price of the Call decreases with the strike. Actually it converges towards 0. Indeed, the higher the strike is, the lower the chance are such that the asset goes over the strike.

What is done here can be done for many, many, many other financial products in order to price them by Monte Carlo.



The code (R):


sample.size <- 365
mu <- 0.1
sigma <- 0.2

a0 <- 1
Asset <- function(sample.size = 365, mu = 0.1, sigma = 0.2, a0 = 1){
  dt <- 1/sample.size
  sdt =sigma*sqrt(dt)
  gauss <- rnorm(sample.size)
  asset <- NULL
  asset[1] <- a0 + mu + sigma * gauss[1]
  test.default <- FALSE
  for(i in 2:365){
    if (!test.default){
      asset[i] <- asset[i-1] * (1 + dt*mu + sdt * gauss[i])
    }
    else{
      asset[i] <- 0
    }
    if (asset[i] <= 0){
      asset[i] <- 0
      test.default <- TRUE
    }
  }
  return(asset)
}


PriceEstimation <- function(t, tf = 365, r = 0.1, strike = 1, n = 1000, mu = 0.1, sigma = 0.2, a0 = 1){
  mean <- 0
  for(i1 in 1:n){
    mean <- mean + exp(-r*(tf-t)/tf) * max(0, (Asset(tf, mu, sigma, a0)-strike))
  }
  mean <- mean/n
  return(mean)
}

res <- NULL
for(i in 1:length(seq(0, 3, 0.05))){
  res[i] <- PriceEstimation(t = 0, tf = 365, r = 0.1, strike = seq(0, 3, 0.05)[i], n = 1000, mu = 0.1, sigma = 0.2, a0 = 1)
}

plot(seq(0,3, 0.05), res,type = 'l', xlab = "Strike Value",  ylab = "Price of the Call")



0

Add a comment

Blog Archive
Translate
Translate
Loading