.chapter22<-function(i=0){ " i Chapter 22: Monte Carlo Simulation - ------------------------------------- 1 Why Monte Carlo Simulation is important 2 Concepts of distributions and cumulative distributions 3 Uniform distribution 4 R functions for a uniform distribution 5 Using rand() function to estimate pi 6 Why generate the same random numbers, concept of seed 7 How to generate the SAME random numbers 8 Standard normal distribution 9 Generating random numbers from a standard normal distribution 10 Normal distributions 11 Why a normal distributions is so important for finance? 12 Generating ransom numbers follow a normal distribution 13 Simulating terminal stock prices and stock price paths 14 Using Monte Carlo Simulation to price a European call (put) 15 Normality test in R 16 Using Monte Carlo Simulation to price Asian options 17 Generating 2 set of correlated random numbers 18 Generating correlated n-stock random numbers 19 Poisson, Binomial distributions 20 Links Example #1:>.c22 # see the above list Example #2:>.c22() # the same as the above Example #3:>.c22(1) # see the 1st explanation ";.zchapter22(i)} .n22chapter<-20 .zchapter22<-function(i){ if(i==0){ print(.c22) }else{ .printEachQ(22,i,.n22chapter) } } .c22<-.chapter22 .C22EXPLAIN1<-"Why Monte Carlo Simulation is important ///////////////////////////////// Monte Carlo simulation produces distributions of possible outcome values. By using probability distributions, variables can have different probabilities of different outcomes occurring. Probability distributions are a much more realistic way of describing uncertainty in variables of a risk analysis. What is Monte Carlo simulation? Monte Carlo simulation is a computerized mathematical technique that allows people to account for risk in quantitative analysis and decision making. The technique is used by professionals in such widely disparate fields as finance, project management, energy, manufacturing, engineering, research and development, insurance, oil & gas, transportation, and the environment. Monte Carlo simulation furnishes the decision-maker with a range of possible outcomes and the probabilities they will occur for any choice of action. It shows the extreme possibilities the outcomes of going for broke and for the most conservative decision along with all possible consequences for middle-of-the-road decisions. The technique was first used by scientists working on the atom bomb; it was named for Monte Carlo, the Monaco resort town renowned for its casinos. Since its introduction in World War II, Monte Carlo simulation has been used to model a variety of physical and conceptual systems. http://www.palisade.com/risk/monte_carlo_simulation.asp ///////////////////////////////// " .C22EXPLAIN2<-"Concepts of distributions and cumulative distributions ///////////////////////////////// Properties of a probability: 1) prob(i) >=0 2) summation of all is 1 Case #1: Toss a coin, what is the probability of a head? head is 50%, i.e., 50% vs. 50% Case #2: Toss a dice with values from 1 to 6 Each number has 1/6 probability What is the prob. to get a number less than 4? 1/6 + 1/6 + 1/6 = 3/6 =1/2 1/6 is probability of each number, while 1/2 is the cumulative probability (distribution) Case #3: uniform distribution: assume density function is 1/a and choose a=5 from 1 to 5, the density is 1/5. total sum is 5*1/5 =1 what does it mean that the density function is 1/5 (0.2)? Take a value between 0 and 5, such as x=1.5 and choose a small number about called deltaX. The probability that a random number falls into this small strip = f(x) * deltaX= 0.2*detlaX Case #4: Standard normal distribution The density function is 1 f(x) = --------- * exp(-0.5 *x^2) sqrt(2*pi) When x=0, f(0) = 1/sqrt(2*3.1425926)=0.3988788 Using Excel =normdist(0,0,1,FALSE) = 0.398942280401433 Case #5: Cumulative standard normal distribution. Let's choose a value of 1 meaning: summation of all small strips from negative infinity up to x=1 Using Excel =normdist(1,0,1,TRUE) -> 0.841344746068543 ///////////////////////////////// " .C22EXPLAIN3<-"Uniform distribution ///////////////////////////////// Uniform distribution is the distribution evenly within a range. Uniform probability distributions arise when every outcome in the sample space has the same probability. Examples: toss a dice ///////////////////////////////// " .C22EXPLAIN4<-"R functions for a uniform distribution ///////////////////////////////// dunif(x, min = 0, max = 1, log = FALSE) punif(q, min = 0, max = 1, lower.tail = TRUE, log.p = FALSE) qunif(p, min = 0, max = 1, lower.tail = TRUE, log.p = FALSE) runif(n, min = 0, max = 1) Get a set of random numbers ------------------------- runif(n) > set.seed(123) > runif(10) [1] 0.2875775 0.7883051 0.4089769 0.8830174 0.9404673 0.0455565 0.5281055 [8] 0.8924190 0.5514350 0.4566147 runif(n,min,max) > set.seed(123) > runif(10,1,50) [1] 15.091298 39.626952 21.039869 44.267853 47.082897 3.232268 26.877169 [8] 44.728533 28.020316 23.374122 > Density function -------------------- dunif(x) Cumulative density distribution ----------------------------- punif(x) ///////////////////////////////// " .C22EXPLAIN5<-"Using rand() function to estimate pi ///////////////////////////////// Methodology (logic): 1) Let's draw a square first, assume that the size is d. we know that the area will be S(square)= d * d 2) Let's draw a circle inside the square We know that the area of the circle will be S(circle) = pi * r^2 where d=2r S(square) = d * d =(2r)* (2r) = 4 *r^2 3) S(circle)/ S(square) = pi *r^2 / (4*r^2) = pi/4 pi = 4 * s(circle) /s(square) Procedure: Step 1: we generate a pair of random number between 0 and 1. let call them x1 and y1 Step 2: if x1^2 + y1^2 <= 1 then n1=n1+1 Step 3: repeat this 500 times Step 4: the value of pi pi=n1/n ///////////////////////////////// " .C22EXPLAIN6<-"why generate the same random numbers, concept of seed ///////////////////////////////// It seems a condition that people want to generate the same random numbers. However, many times we need to have the same ransom numbers. 1) An instructor wants to show how to calculate the mean and standard deviation of 10 random numbers drawn from a uniform distribution. If all students have the same 10 values, it is less problematic. 2) If we design a complex experiment. It is a good idea to fix our input variables to have the same intermediate values. Otherwise, we would not know the difference is caused by our next step or by our input values. Concept of seed --------------- When choose the same seed, different users would get the same set of random numbers. >set.seed(12345) ///////////////////////////////// " .C22EXPLAIN7<-"How to generate random numbers with a seed ///////////////////////////////// For example, we want to generate 10 random numbers from a uniform distribution (from 0 to 1) with a seed of 1234. >set.seed(1234) >x<-runif(10) > x [1] 0.113703411 0.622299405 0.609274733 0.623379442 [5] 0.860915384 0.640310605 0.009495756 0.232550506 [9] 0.666083758 0.51425114 ///////////////////////////////// " .C22EXPLAIN8<-"standard normal distributions: ///////////////////////////////// The standard normal distribution is a probability distribution that associates the normal random variable X with a cumulative probability . The normal distribution is defined by the following equation: 1 f(x) = --------- * exp(-0.5 *x^2) sqrt(2*pi) where f(x) is the density function, x is the input pi is 3.1415926 Manually, we could calculate a few values. x=0 -> f(x) = 1/sqrt(2*3.1415925) -> 0.3989423 x=1 -> f(1) = 1/sqrt(2*3.1415926)*exp(-0.5) -> 0.2419707 x=-1 -> f(-1) = 1/sqrt(2*3.1415926)*exp(-0.5*(-1)^2) -> 0.2419707 ///////////////////////////////// " .C22EXPLAIN9<-"generate random numbers from a normal distributions: ///////////////////////////////// >set.seed(1111) >x<-rnorm(100) ///////////////////////////////// " .C22EXPLAIN10<-"normal distributions: ///////////////////////////////// The normal distribution is a probability distribution that associates the normal random variable X with a cumulative probability . The normal distribution is defined by the following equation: 1 f(x) = ---------------- * exp(-0.5 *[(x-mean)/(sigma))]^2) sqrt(2*pi*sigma^2) ///////////////////////////////// " .C22EXPLAIN11<-"why the normal distributions is so important for finance? ///////////////////////////////// There are most important reasons: 1) if our model is good, then the error terms should follow a normal distribution with mean error is zero. 2) In finance, we assume that stock price follow a log-normal distribution while stock returns follow a normal distribution. ///////////////////////////////// " .C22EXPLAIN12<-"How to generate a set of random numbers with a seed? ///////////////////////////////// Generate a standard normal distribution (mean=0,std=1) ---------------------------------- >set.seed(1111) >x<-rnorm(100) > head(x,10) [1] -0.08658011 1.32252443 0.63970204 1.17478657 [5] 0.11629031 -2.93084636 0.67750806 1.11777194 [9] 1.38404752 1.28394086 > Generate a normal distribution (mean=mu,std=std) ---------------------------------- >mu<-0.1 >td<-0.3 >set.seed(22) >x<-rnorm(100,mean=mu,sd=std) > head(x,10) [1] -0.05364173 0.84555510 0.40234785 0.18784437 [5] 0.03731219 0.65742772 0.08019208 0.05117051 [9] 0.04004180 0.19016852 ///////////////////////////////// " .C22EXPLAIN13<-"Simulate stock prices (path) ///////////////////////////////// Given parameters. st : stock price at time t T : maturity (in years) n : number of steps r : mean return or (risk-free rate) sigma : volatility of the stock deltaT= T/n a small interval 1) simulate s(t+1), s(t+2) etc. --------------------------- s(t+1)=st* exp((r-0.5*sigma^2 )*deltaT + sigma*random*sqrt(deltaT)) where: s(t+1) is the stock price at time t+1 random is a random number drawn from a standard normal distribution A better formula: s(t+1) = s(t) * exp( a + b) a = (r-0.5*sigma^2 )*deltaT b = sigma*random*sqrt(deltaT) 2) R command lines to simulate terminal stock price (sT) at time T ------------------------------------------------------------- sT=st* exp((r-0.5*sigma^2 )*T + sigma*rnorm(1)*sqrt(T)) or a = (r-0.5*sigma^2 )*T b = sigma*rnorm(1)*sqrt(T) sT = st * exp( a + b) # R program sT_f<-function(S,r,T,sigma){ S*exp((r-0.5*sigma*sigma)*T+sigma*rnorm(1)*sqrt(T)) } 3) R program to simultae stock prices ------------------------------------- sT_path_f<-function(S,r,T,n,sigma){ delta_T<-T/n ST<-seq(1:n)*0 ST[1]<-S for (i in 2:n) ST[i]<-ST[i-1]*exp((r-0.5*sigma*sigma)*delta_T+sigma*rnorm(1)*sqrt(delta_T)) return(ST) } ///////////////////////////////// " .C22EXPLAIN14<-"Using Monte Carlo Simulation to price a European call (put) ///////////////////////////////// -------------------------------------------------------------- To price a call, we need the following 5 parameters: s: today's stock price X: exercise price T: maturity date in years r: risk-free rate (compounded continuously) sigma : volatility payoff (call) = Max (sT-X,0) where ST is the terminal stock price at time T bsCallSimulation<-function(s,x,r,T,sigma,n=100){ sT<-s*exp((r-sigma*sigma/2.0)*T + sigma*rnorm(n)*sqrt(T)) payoff<-((sT-x) + abs(sT-x))/2 mean(payoff)*exp(-r*T) } bsCallSimulation(42,40,0.1,0.5,0.2,5000) #[1] 4.728094 bsCall<-function(s,x,T,r,sigma){ d1 = (log(s/x)+(r+sigma*sigma/2.)*T)/(sigma*sqrt(T)) d2 = d1-sigma*sqrt(T) s*pnorm(d1)-x*exp(-r*T)*pnorm(d2) } > bsCall(42,40,0.1,0.5,0.2) [1] 4.015032 ///////////////////////////////// " .C22EXPLAIN15<-"Normality test in R ///////////////////////////////// Test 1: we generate a set of data from a normal ------------------- library(dplyr) set.seed(12345) x<-rnorm(500) shapiro.test(x) Shapiro-Wilk normality test data: x W = 0.99656, p-value = 0.3628 Test 2: Market return minus Risk-free from Fama-French data set --------------------------- library(dplyr) x<-.showff3Monthly(0) > head(x) DATE MKT_RF SMB HML RF 1 1926-07-01 0.0296 -0.0230 -0.0287 0.0022 2 1926-08-01 0.0264 -0.0140 0.0419 0.0025 3 1926-09-01 0.0036 -0.0132 0.0001 0.0023 4 1926-10-01 -0.0324 0.0004 0.0051 0.0032 5 1926-11-01 0.0253 -0.0020 -0.0035 0.0031 6 1926-12-01 0.0262 -0.0004 -0.0002 0.0028 shapiro.test(x$MKT_RF) Shapiro-Wilk normality test data: x$MKT_RF W = 0.91694, p-value < 2.2e-16 ///////////////////////////////// " .C22EXPLAIN16<-"Using Monte Carlo Simulation to price Asian options ///////////////////////////////// -------------------------------------------------------------- Asian options are type of path dependent options while European and American options are path independent. Here is an example of European call. Assume that terminal stock price is $50 while the exercise price is $45. Since the payoff of a call is Max(sT-X,0), we have Max(50-45,0) = 5 In other words, our payoff is $5. This payoff has nothing to do with the path. Payoff(Asian options, average price) = Max(s_average - X, 0) where s_average is the average price before maturity (T) For many cases, Asian options are more meaningful for hedging. Let's use a simple example. A refinery uses huge amount of crude oil. If the price increase, the refinery would suffer a loss if no hedge. The company could buy a call option. Let's look at the two payoffs: Payoff( call) = Max(sT-X,0) payoff( Asian call) = Max(s_average -X,0) Since, the company use crude oil every day, they should be more concerned about the average daily price they paid instead of the end of year price. R program ----------- asian_call_ave_price<-function(S,r,T,n_steps,sigma,x,n_trial){ payoff<-payoff<-rep(0,time=n_trial) # initialization for(i in 1:n_trial){ Saverage<- mean(sT_path_f(S,r,T,n_steps,sigma)) payoff[i]<-max(Saverage-x,0) } return(mean(payoff)*exp(-r*T)) } > asian_call_ave_price(42,0.1,0.5,1000,0.2,40,100) [1] 3.247508 ///////////////////////////////// " .C22EXPLAIN17<-"Generate 2 sets of correlated random numbers ///////////////////////////////// Assume that the correlation (coefficient) between two random series is rho. Step 1: Generate two uncorrelated random time series (x1, x2) Step 2: Two correlated random time series are y1 and y2,where: y1=x1 y2= x`*rho + x2*sqrt ( 1- rho^2 ) Below is one example. Assume that rho=0.3 > set.seed(12345) > x1<-rnorm(1000) > x2<-rnorm(1000) > cor(x1,x2) [1] 0.03877533 > rho<-0.3 > y1<-x1 > y2<-x1*rho + x2*sqrt(1-rho^2) > cor(y1,y2) [1] 0.3307948 # method #2: # --------------- require(MASS) out <- mvrnorm(50, mu = c(0,0), Sigma = matrix(c(1,0.3,0.3,1), ncol = 2),empirical = TRUE) cor(out[,1],out[,2]) [1] 0.3 # Method #3 library(mvtnorm) rho<-0.3 n<-1000 out <- rmvnorm(n,c(0,0),matrix(c(1,rho,rho,1),2,2)) dim(out) [1] 10000 2 > cor(out[,1],out[,2]) [1] 0.2857123 ///////////////////////////////// " .C22EXPLAIN18<-"Generating correlated n sets of random numbers ///////////////////////////////// Assume that we have three stocks and their correlation matrix is shown below. 1 rho12 rho13 rho = rho21 1 rho23 rho31 rho32 1 library(mvtnorm) rho12<-0.2 rho21<-rho12 rho13<-0.3 rho31<-rho13 rho32<-0.1 rho23<-rho32 rho<-c(1,rho12, rho13, rho21, 1, rho23, rho31, rho32, 1) rho2<-matrix(rho,3,3) n<-1000 out <- rmvnorm(n,c(0,0,0),rho2) > cor(out) [,1] [,2] [,3] [1,] 1.0000000 0.2336049 0.2643637 [2,] 0.2336049 1.0000000 0.1313078 [3,] 0.2643637 0.1313078 1.0000000 ///////////////////////////////// " .C22EXPLAIN19<-"Poisson, Binomial distributions ///////////////////////////////// An event can occur 0, 1, 2, times in an interval. The average number of events in an interval is called \"lambda\". Lambda is the event rate, also called the rate parameter. The probability of observing k events in an interval is given by the equation lambda ^k * exp(-lambda) P(k events in interval) = --------------------------- k! where \"lambda\" is the average number of events per interval. exp() is the exponential function. exp(x) could be written 2.71828^x k takes values 0, 1, 2, ... k! = k * (k-1)* (k-2) * ... 3*2*1 ///////////////////////////////// " .C22EXPLAIN20<-"Links ///////////////////////////////// Boyle, P., Phelim, 1977, Options: A Monte Carlo Approach, JFE 4,323-338, http://canisius.edu/~yany/boyle.shtml. Kurtverstegen,D.,2013,Simulation: simulating uncorrelated & correlated random variables, https://kurtverstegen.wordpress.com/2013/12/07/simulation/ Hertz D. B.,1964,Risk Analysis in Capital Investment, Harvard Business Re view 46(1),95-106. Joe, Stephen and Frances Kuo,2005, Sobol sequence generator http://web.maths.unsw.edu.au/~fkuo/sobol/index.html Miller, Craig S., 2003, Normal Quantile Plots in Excel http://facweb.cs.depaul.edu/cmiller/it223/normQuant.html Xiong, Changwei, 2015 http://www.cs.utah.edu/~cxiong/ Wolfram, Mathworld, Random numbers, http://mathworld.wolfram.com/topics/RandomNumbers.html Even Simpler Multivariate Correlated Simulations, 2010 https://cerebralmastication.com/2010/08/even-simpler-multivariate-correlated-simulations/ ///////////////////////////////// " .C22EXPLAIN21<-"After midTerm survey ///////////////////////////////// http://canisius.edu/~yany/doc/afterMidTermSurvey.txt x<-read.table('http://canisius.edu/~yany/doc/afterMidTermSurvey.txt',sep='*',nrows=9) print(x,right=F) y<-read.csv('http://canisius.edu/~yany/doc/afterMidTermSurvey.txt',skip=11) ///////////////////////////////// "