About

In this worksheet we look at different distribution functions, sampling methods, and probability calculations. Next we consider a calculation of a Europen option using Monte Carlo simulation, and compare results to calculation using Black-Scholes.

Setup

Remember to always set your working directory to the source file location. Go to ‘Session’, scroll down to ‘Set Working Directory’, and click ‘To Source File Location’. Read carefully the below and follow the instructions to complete the tasks and answer any questions. Submit your work to RPubs as detailed in previous notes.

Note

Always read carefully the instructions on Sakai. For clarity, tasks/questions to be completed/answered are highlighted in red color (visible in preview) and numbered according to their particular placement in the task section. Quite often you will need to add your own code chunk.

Execute all code chunks, preview, publish, and submit link on Sakai follwoing the naming convention. Make sure to add comments to your code where appropriate. Use own language!

Task 1: Distribution, Sampling & Probability

#Install package quantmod 
if(!require("quantmod",quietly = TRUE))
  install.packages("quantmod",dependencies = TRUE, repos = "https://cloud.r-project.org")
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Version 0.4-0 included new data defaults. See ?getSymbols.

Consider the free cash flow example with an estimated range set by a minimum of -$25M and maximum $275M. We start by generating a sample random set. Following is an example to generate a random set of 10 data points from a uniform distribution with probability \(\frac{1}{b-a}\), a minimum of a, and a maximum of b.

ud = runif(10, min=5, max=15) # runif is in reference to uniform distribution

# Other R functions needed for this task are hist(), mean(), sd(), punif(), runif(),  pnorm(), and rnorm(). Check the Help command in R for more details about these functions.

##### 1A) Assuming a uniform probability distribution, generate a sample of 100 random observations (or deviates) representing the cash flows 

ud = runif(100, min=-25, max=275)

Next we plot a histogram to describe the generated sample data, and calculate the mean and standard deviation of the sample using the proper R-functions

##### 1B) Plot a histogram and calculate the mean and standard deviation of your sample. Compare the mean and standard deviation of the sample to the theoretical values obtained using the proper formulas. 

hist(ud)

mean(ud)
## [1] 129.8297
sd(ud)
## [1] 91.17722

Mean = (1/600)*(275^2 - 25^2) = 125 SD = [(1/900)*(6,750,000)]^(1/2) = 86.60

The proper R formula generates a mean of 128.1435, which is slightly higher than the manually calculated value of 125. The R formula for standard deviation generated a value of 86.28194, which is slightly lower than the calculated value of 86.60.

##### 1C) Repeat 1A & 1B above by increasing the sample size to 1000. Share your insights. 

ud2 = runif(1000, min=-25, max=275)
hist(ud2)

mean(ud2)
## [1] 123.1715
sd(ud2)
## [1] 85.2861

Once again the R-generated mean value is slightly higher (126.4429 > 125), and the R-generated SD value is slightly lower (86.4416 < 86.60). However, increasing the number of observations (deviates) has brought both R-generated values closer the manually calculated values.

Given the characteristics of a probability distribution we should be able to compute various probability scenarios using the proper functions in R.

##### 1D) Calculate the probability that the free cash flow is negative 

punif(ud,min = -25, max = 275)
##   [1] 0.922445024 0.383631490 0.988438698 0.595904435 0.849784181
##   [6] 0.690483740 0.540643191 0.020614671 0.817021230 0.068361017
##  [11] 0.028765324 0.810722947 0.572781981 0.987591160 0.979849273
##  [16] 0.380034859 0.445267439 0.392599422 0.175083074 0.492855433
##  [21] 0.077083161 0.366275597 0.507001918 0.664608405 0.583232434
##  [26] 0.050940605 0.050424765 0.035294121 0.305347376 0.486080161
##  [31] 0.954857417 0.340759546 0.950416612 0.842973220 0.003301476
##  [36] 0.181231878 0.283491711 0.171390339 0.012167083 0.651725037
##  [41] 0.021278373 0.319628625 0.620492682 0.323914157 0.239290486
##  [46] 0.461931637 0.752578802 0.481364705 0.189941621 0.458547495
##  [51] 0.935393898 0.974501021 0.959687050 0.620175711 0.582293891
##  [56] 0.793912342 0.671684026 0.542862226 0.828572135 0.553717180
##  [61] 0.185581565 0.062915435 0.568671996 0.591927521 0.967568667
##  [66] 0.349253671 0.857178868 0.835794065 0.185461662 0.773594770
##  [71] 0.038112570 0.293770561 0.787585434 0.587048904 0.729331765
##  [76] 0.266431581 0.327008594 0.043967220 0.462268322 0.477183342
##  [81] 0.737035132 0.283915763 0.813906264 0.786377291 0.579840575
##  [86] 0.978816877 0.430627710 0.872224408 0.986544249 0.332720335
##  [91] 0.808247064 0.695660428 0.038104434 0.100950168 0.733596477
##  [96] 0.946021136 0.216655124 0.763908899 0.483081591 0.637755092

Prob Neg = [(0 - (-25))/(275-(-25))] = 0.0833

We will now repeat the above exercises 1A, 1B, 1C using instead a normal distribution.

##### 1E) Repeat steps A-C for the case of a portfolio daily returns with normal probability distribution,a mean=1.2% and a standard deviation= 3.7% 

nd = rnorm(100, mean = 1.2, sd = 3.7)
hist(nd, breaks = 20)

mean(nd)
## [1] 0.8707093
sd(nd)
## [1] 3.890054
nd2 = rnorm(1000, mean = 1.2, sd = 3.7)
hist(nd2, breaks = 20)

mean(nd2)
## [1] 1.104776
sd(nd2)
## [1] 3.581415

Similarly we should be able to compute various probability senarios with our obtained normal distribution.

##### 1F) Calculate the probability that returns will be negative using the values for mean and standard deviation as in 1E 

negnd = nd[nd<0]
p = length(negnd)/length(nd)

Probability of negative returns is 0.39

negnd2 = nd2[nd2<0]
p2 = length(negnd2)/length(nd2)

Probability of negative returns is 0.375

pnorm(0, mean = 1.2, sd = 3.7, lower.tail = "true")
## Warning in pnorm(0, mean = 1.2, sd = 3.7, lower.tail = "true"): NAs
## introduced by coercion
## [1] 0.3728463

0.3728463 is the most accurate probability of negative returns, given 1000 observations.

##### 1G) Repeat the calculation in 1F using instead the standard Z-value. Share insights. P(r < 0) = P(z < –0.324) = N(–0.324) = 0.3745

Task 2: MC Simulation & European Option Pricing

Follow the Algorithm 5.2 example on p 167 (*) to calculate the price of a European option using a MC simulation. Note that the code in the book example is missing one detail and a correction. Those are left for your investigation.

##### 2A) Identify and explain the nature of the missing detail and the correction needed.

Missing: Sold = S This will reset the initial value for the simulated path to the value of S, 100 in this case. Otherwise the initial value for subsequent simulations will be the final value of the previous simulation.

Correction: E = rnorm(0,1) –> E = rnorm(1, mean =0, sd = 1) This correction clarifies for the program the size of the input vector is 1 and mean is 0, whereas the original form would have established the size of the input vector as 0, yielding a value of 0

Given the MC simulation we should be able to price a European option. Note the introduction of user defined function in the book. Below is an example of a user defined function and usage. The function returns the squared value.

myfunction <- function(x=2){
  y=x^2
return(y)}

myfunction(x=5)
## [1] 25

##### 2B) Use MC simulation to price a European Call option on a stock with initial price of $155, strike price $140, a time to expiration equal to six months, a risk-free interest rate of 2.5% and a volatility of 23%. Consider a number of periods n=100 and a number of simulations m=1000 

#User defined function for Euler MC
EulerMC = function(Type = "c", S=100, K=100, T=1/12, r=0.1, sigma=0.25, n=300, m=100, dt=T/n){
  sum=0
  for(i in 1:m){# of paths
    Sold=S
    for(j in 1:n){# of periods in path
      E=rnorm(1,mean = 0,sd = 1);
      Snew=Sold+r*Sold*dt+sigma*Sold*sqrt(dt)*E;}
    if(Type == "c"){payoff = max(Snew-K,0)}
    else if(Type == "p"){payoff = max(K-Snew,0)}
    else{payoff=max(Snew-K,0)} #default
    sum = sum + payoff}
  OptionValue = (sum*exp(-r*T))/m;
  return(OptionValue)}

f <- EulerMC("c", 155, 140, 1/2, 0.025, 0.23, 100, 1000)

f = 14.8445

##### 2C) Write the mathematical representation of the discrete pricing equation modeled in the MC simulation. Explain what each variable in the equation represents, and provide the associated numerical value.

##### 2D) Compare the price obtained from the MC simulation to the option price using the Black-Scholes function pricing GBSOption(). Share insights. 

install.packages("fOptions")
## Installing package into '/home/rstudio-user/R/x86_64-pc-linux-gnu-library/3.5'
## (as 'lib' is unspecified)
library("fOptions")
## Loading required package: timeDate
## Loading required package: timeSeries
## 
## Attaching package: 'timeSeries'
## The following object is masked from 'package:zoo':
## 
##     time<-
## Loading required package: fBasics
## 
## Attaching package: 'fBasics'
## The following object is masked from 'package:TTR':
## 
##     volatility
GBSOption(TypeFlag = "c", S=155, X=140, Time = 1/2, r = 0.025, b = 0.00, sigma = 0.23)
## 
## Title:
##  Black Scholes Option Valuation 
## 
## Call:
##  GBSOption(TypeFlag = "c", S = 155, X = 140, Time = 1/2, r = 0.025, 
##      b = 0, sigma = 0.23)
## 
## Parameters:
##           Value:
##  TypeFlag c     
##  S        155   
##  X        140   
##  Time     0.5   
##  r        0.025 
##  b        0     
##  sigma    0.23  
## 
## Option Price:
##  18.63254 
## 
## Description:
##  Wed Jan 30 22:52:23 2019

The BS option price is found to be 18.63254, whereas the MC simulation calculated a value of 14.8445 using the same parameters.

*http://computationalfinance.lsi.upc.edu