VIEWS: 0 PAGES: 205 POSTED ON: 2/12/2012 Public Domain
SS 2008 LV 1837 Computersimulation for Finance R.Gismondi FIRM-Research Institute for Computational Methods Index • Elements of Probability I 4 • Random Numbers 27 • Elements of Statistics 40 • Elements of Probability II 51 • Generating Random Variables 64 • Brownian Motion I 89 • Geometric Brownian Motion, B&S Model, Pricing of (exotic) options: 1-dim. 113 Folie 2 Riccardo Gismondi Index • Brownian Motion II 151 • Geometric Brownian Motion, pricing of (exotic) options: n-dim. 161 • Square root diffusion process: CIR Model 173 • Heston Model and stochastic Volatility 192 • Pricing of Exotic Options: 205 Structured Products (Equity) Folie 3 Riccardo Gismondi Elements of Probability I Folie 4 Riccardo Gismondi Introduction • Ground-up review of probability necessary to do and understand simulation • Assume familiarity with • Mathematical analysis (especially derivate and integrals) • Some probability ideas (especially probability spaces, random variables, probability distributions) • Outline • Probability – basic ideas, terminology • Random variables, joint distributions Folie 5 Riccardo Gismondi Probability Basics (1) • Experiment – activity with uncertain outcome • Flip coins, throw dice, pick cards, draw balls from urn, … • Drive to work tomorrow – Time? Accident? • Operate a (real) call center – Number of calls? Average customer hold time? Number of customers getting busy signal? • Simulate a call center – same questions as above • Sample space – complete list of all possible individual outcomes of an experiment • Could be easy or hard to characterize • May not be necessary to characterize Folie 6 Riccardo Gismondi Probability Basics (2) • Event – a subset of the sample space • Describe by either listing outcomes, “physical” description, or mathematical description • Usually denote by E, F, E1, E2, etc. • Union, intersection, complementation operations • Probability of an event is the relative likelihood that it will occur when you do the experiment • A real number between 0 and 1 (inclusively) • Denote by P(E), P(E F), etc. • Interpretation – proportion of time the event occurs in many independent repetitions (replications) of the experiment • May or may not be able to derive a probability Folie 7 Riccardo Gismondi Probability Basics (3) • Some properties of probabilities • If S is the sample space, then P(S) = 1 • Can have event E S with P(E) = 1 • If Ø is the empty event (empty set), then P(Ø) = 0 • Can have event E Ø with P(E) = 0 • If EC is the complement of E, then P(EC) = 1 – P(E) • P(E F) = P(E) + P(F) – P(E F) • If E and F are mutually exclusive (i.e., E F = Ø), then P(E F) = P(E) + P(F) • If E is a subset of F (i.e., the occurrence of E implies the occurrence of F), then P(E) P(F) • If o1, o2, … are the individual outcomes in the sample space, then Folie 8 Riccardo Gismondi Probability Basics (4) • Conditional probability • Knowing that an event F occurred might affect the probability that another event E also occurred • Reduce the effective sample space from S to F, then measure “size” of E relative to its overlap (if any) in F, rather than relative to S • Definition (assuming P(F) 0): • E and F are independent if P(E F) = P(E) P(F) • Implies P(E|F) = P(E) and P(F|E) = P(F), i.e., knowing that one event occurs tells you nothing about the other • If E and F are mutually exclusive, are they independent? Folie 9 Riccardo Gismondi Random Variables • One way of quantifying, simplifying events and probabilities • A random variable (RV) “is a number whose value is determined by the outcome of an experiment”. • Technically, a function or mapping from the sample space to the real numbers, but can usually define and work with a RV without going all the way back to the sample space • Think: RV is a number whose value we don’t know for sure but we’ll usually know something about what it can be or is likely to be • Usually denoted as capital letters: X, Y, W1, W2, etc. • Probabilistic behavior described by distribution function Folie 10 Riccardo Gismondi Discrete vs. Continuous RVs • Two basic “flavors” of RVs, used to represent or model different things • Discrete – can take on only certain separated values • Number of possible values could be finite or infinite • Continuous – can take on any real value in some range • Number of possible values is always infinite • Range could be bounded on both sides, just one side, or neither Folie 11 Riccardo Gismondi Discrete Distributions (1) • Let X be a discrete RV with possible values (range) x1, x2, … (finite or infinite list) • Probability mass function (PMF) p(xi) = P(X = xi) for i = 1, 2, ... • The statement “X = xi” is an event that may or may not happen, so it has a probability of happening, as measured by the PMF • Can express PMF as numerical list, table, graph, or formula • Since X must be equal to some xi, and since the xi’s are all distinct, Folie 12 Riccardo Gismondi Discrete Distributions (2) • Cumulative distribution function (CDF) – probability that the RV will be a fixed value x: • Properties of discrete CDFs 0 F(x) 1 for all x As x –, F(x) 0 These four properties are also true of As x +, F(x) 1 continuous CDFs F(x) is nondecreasing in x F(x) is a step function continuous from the right with jumps at the xi’s of height equal to the PMF at that xi Folie 13 Riccardo Gismondi Discrete Distributions (3) • Computing probabilities about a discrete RV – usually use the PMF • Add up p(xi) for those xi’s satisfying the condition for the event • With discrete RVs, must be careful about weak vs. strong inequalities – endpoints matter! Folie 14 Riccardo Gismondi Discrete Expected Values • Data set has a “center” – the average (mean) • RVs have a “center” – expected value • Also called the mean or expectation of the RV X • Weighted average of the possible values xi, with weights being their probability (relative likelihood) of occurring • What expectation is not: The value of X you “expect” to get E(X) might not even be among the possible values x1, x2, … • What expectation is: Repeat “the experiment” many times, observe many X1, X2, …, Xn E(X) is what converges to (in a certain sense) as n Folie 15 Riccardo Gismondi Discrete Variances and Standard Deviations • Data set has measures of “dispersion” – • Sample variance • Sample standard deviation • RVs have corresponding measures • Other common notation: • Weighted average of squared deviations of the possible values xi from the mean • Standard deviation of X is • Interpretation analogous to that for E(X) Folie 16 Riccardo Gismondi Continuous Distributions (1) • Now let X be a continuous RV • Possibly limited to a range bounded on left or right or both • No matter how small the range, the number of possible values for X is always (uncountably) infinite • Not sensible to ask about P(X = x) even if x is in the possible range • Technically, P(X = x) is always 0 • Instead, describe behavior of X in terms of its falling between two values Folie 17 Riccardo Gismondi Continuous Distributions (2) • Probability density function (PDF) is a function f(x) with the following three properties: f(x) 0 for all real values x The total area under f(x) is 1: For any fixed a and b with a b, the probability that X will fall between a and b is the area under f(x) between a and b: • Fun facts about PDFs • Observed X’s are denser in regions where f(x) is high • The height of a density, f(x), is not the probability of anything – it can even be > 1 • With continuous RVs, you can be sloppy with weak vs. strong inequalities and endpoints Folie 18 Riccardo Gismondi Continuous Distributions (3) • Cumulative distribution function (CDF) - probability that the RV will be a fixed value x: F(x) may or may not have a closed-form formula • Properties of continuous CDFs 0 F(x) 1 for all x As x –, F(x) 0 These four properties are also true of As x +, F(x) 1 discrete CDFs F(x) is nondecreasing in x F(x) is a continuous function with slope equal to the PDF: f(x) = F'(x) Folie 19 Riccardo Gismondi Continuous Expected Values, Variances and Standard Deviations • Expectation or mean of X is • Roughly, a weighted “continuous” average of possible values for X • Same interpretation as in discrete case: average of a large number (infinite) of observations on the RV X • Variance of X is • Standard deviation of X is Folie 20 Riccardo Gismondi Joint Distributions (1) • So far: Looked at only one RV at a time • But they can come up in pairs, triples, …, tuples, forming jointly distributed RVs or random vectors • Input: (T, P, S) = (type of part, priority, service time) • Output: {W1, W2, W3, …} = output process of times in system of exiting parts • One central issue is whether the individual RVs are independent of each other or related • Will take the special case of a pair of RVs (X1, X2) • Extends naturally (but messily) to higher dimensions Folie 21 Riccardo Gismondi Joint Distributions (2) • Joint CDF of (X1, X2) is a function of two variables Replace “and” with “,” • Same definition for discrete and continuous • If both RVs are discrete, define the joint PMF • If both RVs are continuous, define the joint PDF f(x1, x2) as a nonnegative function with total volume below it equal to 1, and • Joint CDF (or PMF or PDF) contains a lot of information – usually won’t have in practice Folie 22 Riccardo Gismondi Marginal Distributions • What is the distribution of X1 alone? Of X2 alone? • Jointly discrete • Marginal PMF of X1 is • Marginal CDF of X1 is • Jointly continuous • Marginal PDF of X1 is • Marginal CDF of X1 is • Everything above is symmetric for X2 instead of X1 • Knowledge of joint knowledge of marginals – but not vice versa (unless X1 and X2 are independent) Folie 23 Riccardo Gismondi Covariance Between RVs • Measures linear relation between X1 and X2 • Covariance between X1 and X2 is • If large (resp. small) X1 tends to go with large (resp. small) X2, then covariance > 0 • If large (resp. small) X1 tends to go with small (resp. large) X2, then covariance < 0 • If there is no tendency for X1 and X2 to occur jointly in agreement or disagreement over being big or small, then Cov = 0 • Interpreting value of covariance – difficult since it depends on units of measurement Folie 24 Riccardo Gismondi Correlation Between RVs • Correlation (coefficient) between X1 and X2 is • Has same sign as covariance • Always between –1 and +1 • Numerical value does not depend on units of measurement • Dimensionless – universal interpretation Folie 25 Riccardo Gismondi Independent RVs • X1 and X2 are independent if their joint CDF factors into the product of their marginal CDFs: • Equivalent to use PMF or PDF instead of CDF • Properties of independent RVs: • They have nothing (linearly) to do with each other • Independence uncorrelated • But not vice versa, unless the RVs have a joint normal distribution • Important in probability – factorization simplifies greatly • Tempting just to assume it whether justified or not • Independence in simulation • Input: Usually assume separate inputs are indep. – valid? • Output: Standard statistics assumes indep. – valid?!? Folie 26 Riccardo Gismondi Random Numbers Folie 27 Riccardo Gismondi Introduction • The building block of a simulation study is the ability to generate random numbers. • Random number represent the value of a random variable uniformly distributed on (0,1) Folie 28 Riccardo Gismondi Pseudorandom Number Generation • Multiplicative congruential method x0 ......seed a 75 16807 xn a * xn1 mod m m 231 1 xn un m • Mixed congruential method x0 ......seed a 75 16807 m 231 1 xn a * xn1 c mod m c xn un m Folie 29 Riccardo Gismondi Computer exercises ( MATLAB ) • Let x0 5 a 3 m 150 xn a * xn1 mod m find x1 , x2 ,..., x10 • Let x0 3 a 5 m 200 c=7 xn a * xn1 c mod m find x1 , x2 ,..., x10 Folie 30 Riccardo Gismondi Using random numbers to evaluate Integrals (1) • Let g(x) a function and suppose we want to compute : 1 g ( x)dx 0 • If U is uniform distributed over (0,1), then we can express as E g (U ) • If U1....U k are independent uniform (0,1) random variables, it follows that the random variables g U ....g U 1 k are i.i.d variables with mean Folie 31 Riccardo Gismondi Using random numbers to evaluate Integrals (2) • Therefore, by the strong law of large numbers, it follows that, with probability 1: k 1 g (U i ) lim E g (U ) g ( x)dx k k i 1 0 • Thus, we can approximate by generating a large number of U i and taking as our approximation the average value of g ui . This approach to approximate integrals is called Monte Carlo method Folie 32 Riccardo Gismondi Computer exercises (MATLAB) • Use simulation (Monte Carlo method) to approximate the following integrals : 1 x dx 2 1. 0 1 2t 1 dt 2 2. 0 Folie 33 Riccardo Gismondi Using random numbers to evaluate Integrals (3) • If we want to compute: b g ( x)dx a y x a, dx then, by substitution dy ba ba 1 1 g a [b a ] y b a dy h y dy 0 0 h y b a g a [b a ] y Folie 34 Riccardo Gismondi Computer exercises (MATLAB) • Use simulation (Monte Carlo method) to approximate the following integrals : 1 1. cos( x)dx 3 2 e 1 dt t 2 2. 1 3 Folie 35 Riccardo Gismondi Using random numbers to evaluate Integrals (4) • Similarly, if we want to compute: g ( x)dx then, by substitution , 0 1 dx y dy y 2dx x 1 x 1 2 1 1 g 1 h y dy with h y y 0 y2 Folie 36 Riccardo Gismondi Example: the estimation of • Suppose that the random vector (X,Y) is uniformly distributed in the square of area 4 centered in the origin. Let us consider the probability that this random point is contained within the inscribed circle of radius 1: Area of the circle P X , Y is in the circle P X 2 Y 2 1 Area of the square 4 (-1,1) (1,1) If U is uniform on (0,1) then 2U is uniform on (0,2) , and so 2U-1 (0,0) is uniform on (-1,1) (-1,-1) (1,-1) Folie 37 Riccardo Gismondi Example: the estimation of • Therefore generate random numbers U1 and U2 set X 2U1 1 , Y 2U 2 1 and define : 1, if X 2 Y 2 1 I E I P X Y 1 2 2 0, otherwise 4 Hence we can estimate / 4 and hence by the fraction of pairs for which 2u1 1 2u2 1 1 2 2 Folie 38 Riccardo Gismondi Home assignment (Matlab) • Write a code to generate Uniform Random Numbers • Write a code to generate Uniform Random Vectors • Write a code to estimate and study the convergence Folie 39 Riccardo Gismondi Elements of Statistics Folie 40 Riccardo Gismondi Sampling • Statistical analysis – estimate or infer something about a population or process based on only a sample from it • Think of a RV with a distribution governing the population • Random sample is a set of independent and identically distributed (IID) observations X1, X2, …, Xn on this RV • In simulation, sampling is making some runs of the model and collecting the output data • Don’t know parameters of population (or distribution) and want to estimate them or infer something about them based on the sample Folie 41 Riccardo Gismondi Sampling (cont’d.) • Population parameter • Sample estimate Population mean m = E(X) Sample mean Population variance s2 Sample variance Population proportion Sample proportion • Parameter – need to know • Sample statistic – can be whole population computed from a sample • Fixed (but unknown) • Varies from one sample to another – is a RV itself, and has a distribution, called the sampling distribution Folie 42 Riccardo Gismondi Sampling Distributions • Have a statistic, like sample mean or sample variance • Its value will vary from one sample to the next • Some sampling-distribution results • Sample mean If Regardless of distribution of X, • Sample variance s2 E(s2) = s2 • Sample proportion E( ) = p Folie 43 Riccardo Gismondi Point Estimation • A sample statistic that estimates (in some sense) a population parameter • Properties • Unbiased: E(estimate) = parameter • Efficient: Var(estimate) is lowest among competing point estimators • Consistent: Var(estimate) decreases (usually to 0) as the sample size increases Folie 44 Riccardo Gismondi Confidence Intervals • A point estimator is just a single number, with some uncertainty or variability associated with it • Confidence interval quantifies the likely imprecision in a point estimator • An interval that contains (covers) the unknown population parameter with specified (high) probability 1 – a • Called a 100 (1 – a)% confidence interval for the parameter • Confidence interval for the population mean m : s tn-1,1-a/2 is point below which is area X t n 1,1a / 2 1 – a/2 in Student’s t distribution with n n – 1 degrees of freedom • CIs for some other parameters – in text Folie 45 Riccardo Gismondi Confidence Intervals in Simulation • Run simulations, get results • View each replication of the simulation as a data point • Random input random output • Form a confidence interval • Brackets (with probability 1 – a) the “true” expected output (what you’d get by averaging an infinite number of replications) Folie 46 Riccardo Gismondi Hypothesis Tests • Test some assertion about the population or its parameters • Can never determine truth or falsity for sure – only get evidence that points one way or another • Null hypothesis (H0) – what is to be tested • Alternate hypothesis (H1 or HA) – denial of H0 H 0: m = 6 vs. H1: m 6 H0: s < 10 vs. H1: s 10 H0: m1 = m2 vs. H1: m1 m2 • Develop a decision rule to decide on H0 or H1 based on sample data Folie 47 Riccardo Gismondi Errors in Hypothesis Testing H0 is really true H1 is really true Decide H0 No error Type II error (“Accept” H0) Probability 1 – a Probability a is chosen is not controlled – (controlled) affected by a and n Decide H1 Type I Error No error (Reject H0) Probability a Probability 1 – = power of the test Folie 48 Riccardo Gismondi p-Values for Hypothesis Tests • Traditional method is “Accept” or Reject H0 • Alternate method – compute p-value of the test • p-value = probability of getting a test result more in favor of H1 than what you got from your sample • Small p (like < 0.01) is convincing evidence against H0 • Large p (like > 0.20) indicates lack of evidence against H0 • Connection to traditional method • If p < a, reject H0 • If p a, do not reject H0 • p-value quantifies confidence about the decision Folie 49 Riccardo Gismondi Hypothesis Testing in Simulation • Input side • Specify input distributions to drive the simulation • Collect real-world data on corresponding processes • “Fit” a probability distribution to the observed real-world data • Test H0: the data are well represented by the fitted distribution • Output side • Have two or more “competing” designs modeled • Test H0: all designs perform the same on output, or test H0: one design is better than another • Selection of a “best” model scenario Folie 50 Riccardo Gismondi Elements of Probability II Folie 51 Riccardo Gismondi Discrete Random Variables • Binomial Random Variables Suppose that n independent trials, each of which results in success (1) with probability p, are to be performed. If X represents the number of success in the n trials, then X is said to be a Binomial Random Variable with parameters (p, n). Its PMF is given by n i Pi P X i p (1 p)ni i A binomial (1, p) random variable is called Bernoulli random variable Folie 52 Riccardo Gismondi Discrete Random Variables • Binomial Random Variables For a Binomial Random Variable with parameters (p,n) expectation and variance is given by (Home assignment !) E X np, Var(X)=np 1-p • Poisson Random Variables A random variable that takes on one of the values 0,1,2 3,…. is said to be Poisson Random Variable with parameter if its PMF is given by n 1 i Pi P X i e , e lim 1 i! n n Folie 53 Riccardo Gismondi Discrete Random Variables • Poisson Random Variables Poisson random variables have a wide range of applications. One reason is because such random variables may be used to approximate the distribution of the number of success in a large number of trials, when each trial has a small probability of success. To see why, suppose that X is a binomial random variable (n,p) and let np n i Then Pi P X i p (1 p)ni i n! p i (1 p) ni (n i )!i ! Folie 54 Riccardo Gismondi Discrete Random Variables • Poisson Random Variables ni n! i 1 (n i)!i ! n n n i 1 n(n 1)..(n i 1) n ni i ! i 1 n Folie 55 Riccardo Gismondi Discrete Random Variables • Poisson Random Variable Now for n p n(n 1)...(n i 1)! n i 1 e , 1 , 1 1 n n i n Hence for n p i Pi P X i e i! Since the mean and variance of a binomial random variable Y are given by E Y np, Var(Y)=np 1-p np Then for a Poisson variable with parameter , we get E X Var(X)= Folie 56 Riccardo Gismondi Continuous Random Variables • Uniformly Distributed Random Variables A random variable X is said to be Uniformly Distributed over the interval (a,b) if its pdf is given by 1 if a<x<b f ( x) b a 0 otherwise For the mean we have : b2 a 2 b a b 1 E X xdx 2(b a) 2 ba a Folie 57 Riccardo Gismondi Continuous Random Variables • Uniformly Distributed Random Variables For the variance we have: b3 a3 b2 a 2 ab b 1 ba E X 2 x 2dx 3(b a) 3 a b a b a 2ab 2 2 2 E X 2 2 4 b 2 a 2 ab b 2 a 2 2ab 1 Var ( X ) E X 2 E X b a 2 2 3 4 12 Folie 58 Riccardo Gismondi Continuous Random Variables • Normal Random Variables A random variable X is said to be normally distributed if its pdf is given by x m 2 f x 1 s 2 e 2s 2 - x It is not difficult to show (Home assignment !) that : E X m Var(X)=s 2 Folie 59 Riccardo Gismondi Continuous Random Variables • Normal Random Variables An important fact about normal random variables is that if X is normally distributed then (Home assignment !) X N ( m , s 2 ) aX b N (am b, a 2s 2 ) It follows from this that if X m X N (m ,s ) Z 2 N (0,1) s Such a random variable is said to have standard normal distribution Folie 60 Riccardo Gismondi Continuous Random Variables • Normal Random Variables Let denote the CDF of a X N (0,1) , that is x t2 x 1 2 e 2 dt - x We have the Central Limit Theorem : Let X 1 , X 2 ,... a sequence of i.i.d random variables with finite mean m and finite variances 2: then X1 X 2 .. X n nm lim P x x n s n Folie 61 Riccardo Gismondi Continuous Random Variables • Exponential Random Variables A random variable X is said to be exponentially distributed with parameter if its pdf is given by f x e x - x Its CDF is given by x F x e dt =1-e t xt 0 x 0 It is not difficult to show (Home assignment !) that : 1 1 E X Var(X)= 2 Folie 62 Riccardo Gismondi Continuous Random Variables • Exponential Random Variables The key property of exponential random variable is the “memoryless property”: P X s t X s P X t s,t 0 This equation is equivalent to P X s t P X s P X t which is clearly satisfied whenever X is an exponential random variable, since P X x ex Folie 63 Riccardo Gismondi Generating Random Variables Folie 64 Riccardo Gismondi Discrete Random Variables • The Inverse Transform Method (1) Suppose we want to generate a DRV X having PMF P X x j p j , j=0,1,..., p j 1 j To accomplish this, we generate u U (0,1) and set x , if u p 0 0 x1 , if p0 u p0 p1 . X . j-1 j x j , if pi u pi i=1 i=1 . . Folie 65 Riccardo Gismondi Discrete Random Variables • The Inverse Transform Method (2) Since for 0 a b 1, P a u<b b a we have that j-1 j P X x j = P pi u pi =pi i=1 i=1 And so X has the desired distribution Folie 66 Riccardo Gismondi Computer exercises (MATLAB) 1. Simulate a RV X such that p1 0.2, p2 0.15, p3 =0.25, p4 =0.4 p j P X j 2. Simulate a RV X such that p1 0.05, p1 0.15, p3 0.22, p4 0.08, p5 0.12, p6 0.1, p7 0.1, p8 0.08, p9 0.05, p10 0.03, p11 0.02 p j P X j 3. Simulate a RV X such that 1 p j P X j n Folie 67 Riccardo Gismondi Home assignment (Matlab) • Write a code to generate a DRV with the Inverse Transform Method [U,Hist(U)]=f(data,N) Folie 68 Riccardo Gismondi Continuous Random Variables • The Inverse Transform Algorithm (2) Consider a CRV having PDF F. A general method for generating such a random is based on the following proposition. • Proposition : Let U be a uniform (0,1) RV. For any continuous distribution function F the RV X defined by X F 1 U has a distribution F • Proof : let FX denote the PDF of X F 1 U . Then FX x P X x P F 1 U x Folie 69 Riccardo Gismondi Continuous Random Variables • The Inverse Transform Algorithm (3) Now since F is a PDF it follows that a b F (a) F (b) Hence, we see that FX x P F F 1 U F x = P U F x = F x Folie 70 Riccardo Gismondi Continuous Random Variables • The Inverse Transform Algorithm (4) 1. Supposed we wanted to generate a CRV X with PDF F x xn , 0<x 1 If we let x F 1 u then 1 u F x xn , x u n 2. If X is a exponential RV with rate 1, then is PDF is given by F x 1 e x x F 1 u u F x 1 e x x log 1 u ,or x log u ,since Y U(0,1) 1-Y U(0,1) Folie 71 Riccardo Gismondi Continuous Random Variables • The Box-Muller method for generating Standard Normal Random Variables (1) Let X and Y be independent standard normal RV and let R and denote the polar coordinates of the vector (X,Y) , that is Y R (X,Y) R2 X 2 Y 2 sin Y tan X cos X Folie 72 Riccardo Gismondi Continuous Random Variables • The Box-Muller method for generating Standard Normal Random Variables (2) Since X and Y are independent , their joint PDF is given (why?) by 2 2 x 2 2 y x y 1 1 1 f ( x, y ) e 2 e 2 e 2 2 2 2 To determine the joint density of R 2 and we make the change of variables : d X 2 Y 2 Y Y tan 1 arctan X X Folie 73 Riccardo Gismondi Continuous Random Variables • The Box-Muller method for generating Standard Normal Random Variables (3) Since the Jacobian of the transformation is equal to 2 ( home assignment !) it follows that the joint PDF of R 2 and is d2 1 1 f ( d , ) e 2 2 2 However, this is equal to the product of an exponential density having mean 2 and the uniform density on 0,2 ! It follows that R 2 and are independent with R 2 exp(2), U 0,2 Folie 74 Riccardo Gismondi Continuous Random Variables • The Box-Muller method for generating Standard Normal Random Variables (4) We can now generate a pair of independent standard normal RV with the following algorithm: Step 1 : generate uniform random numbers U 1 and U 2 Step 2: R2 2log U1 , 2U2 X R cos 2log U1 cos 2U 2 Step 3: Y R sin 2log U1 sin 2U 2 Folie 75 Riccardo Gismondi Home assignment (Matlab) • Write a code to generate a standard normal RV with the Box-Muller Method [U,hist(U)]=f(N) Folie 76 Riccardo Gismondi Continuous Random Variables • The Polar method for generating Standard Normal Random Variables (1) Unfortunately, the Box-Mueller Transformation is computationally not very efficient: sine and cosine trigonometric functions Idea: indirect computation of a sine and cosine of a random angle. (-1,1) (1,1) V2 V2 sin 1 R V2 R V1 V2 2 V1 V1 V1 cos R V V 1 (-1,-1) (1,-1) 1 2 2 Folie 77 Riccardo Gismondi Continuous Random Variables • The Polar method for generating Standard Normal Random Variables (2) Thus: V1 X 2 log U cos 2 U 2 log U 1 V1 V2 2 V2 Y 2 log U sin 2 U 2 log U 1 V1 V2 2 If U1 and U 2 are Uniform over (0,1) , then V1 2U1 1 and V2 2U 2 1 are also Uniform over (0,1); setting R 2 V12 V22 S we obtain that Folie 78 Riccardo Gismondi Continuous Random Variables • The Polar method for generating Standard Normal Random Variables (3) 1 V1 2 log S 2 X 2 log S 1 V1 S 2 S 2 log S V2 Y 2 log S 1 V2 S 2 S are independent standard normal when V1,V2 is a randomly chosen point in the circle of radius 1 centered at the origin . Summing up we have the following algorithm to generate a pair of independent standard normal Folie 79 Riccardo Gismondi Continuous Random Variables • The polar method Algorithm STEP 1 : Generate random numbers U and U2 1 STEP 2 : Set V1 2U1 1 V2 2U 2 1 S V 2 V 2 1 2 STEP 3 : if S 1 go to STEP 1 1 2 log S else 2 return the i.s.n X V1 S 1 2 log S 2 Y V2 S Folie 80 Riccardo Gismondi Home assignment (Matlab) • Write a code to generate a standard normal RV with the Polar Method [U,hist(U)]=f(N) Folie 81 Riccardo Gismondi Continuous Random Variables • Generating Multivariate Normal Variables (1) A d-dimensional multivariate normal distribution is characterized by a d-vector m and a d d matrix : we abbreviate it as N , . m must be symmetric and positive semidefinite, meaning that xT x 0 for all x d If is positive semidefinite, then the normal distribution N m, has density 1 xm 1 x m T 1 m , x d 1 e 2 , x d 2 2 2 Folie 82 Riccardo Gismondi Continuous Random Variables • Generating Multivariate Normal Variables (2) The standard d-dimensional multivariate normal distribution N 0, I d with I dthe d d identity matrix, is the special case: 1 T 1 x x 0,I d x d e 2 2 2 If X N m, , then its i-th component X i has distribution N mi , ii with s i ii . The i-th and j-th components have covariance Cov X i X j E X i mi X j m j ij Folie 83 Riccardo Gismondi Continuous Random Variables • Generating Multivariate Normal Variables (3) Any linear transformation of a normal vector is again normal: X N m , AX N Am , AAT , for any d-vector m and d d matrix , and any k d matrix A, for any k. Using any of the methods described before, we can generate standard normal variables Z1 , Z 2 ,...., Z d and assemble them into a vector Z N 0, I . Thus the problem of sampling X from the multivariate normal N m, reduces to finding a matrix A for which AA T Folie 84 Riccardo Gismondi Continuous Random Variables • Generating Multivariate Normal Variables (4) Among all such A, a lower triangular one is particularly convenient because it reduces the calculation of m AZ to the following: X 1 m1 A11Z1 X 2 m2 A21Z1 A22 Z 2 . . X d md Ad 1Z1 Ad 2 Z 2 .. Add Z d A full multiplication of the vector Z by the matrix A would require approximately twice as many multiplications and additions. Folie 85 Riccardo Gismondi Continuous Random Variables • Generating Multivariate Normal Variables (5) A representation of as AAT with A lower triangular is a Cholesky factorization of A . Consider a 2 2 covariance matrix represented as s 2 s 1s 2 1 s 1s 2 s2 2 Assuming s 1 0 and s 2 0 , the Cholesky factor is ( verified ! ) s1 0 A s s 1 2 2 2 Folie 86 Riccardo Gismondi Continuous Random Variables • Generating Multivariate Normals Variables (6) Thus we can generate a bivariate normal distribution N m, by setting X 1 m1 s 1Z1 X 2 m2 s 2 Z1 s 2 1 2 Z 2 For a d-dimensional case we need to solve : A11 A11 A12 . Ad 1 A21 A22 A22 . Ad 2 . . . . . Ad 1 Ad 2 . Add Add Folie 87 Riccardo Gismondi Home assignment (Matlab) • Write a code to generate Multivariate Normal RVs X1, X 2 ,.. X d , hist ( X1 ), hist ( X 2 ),..hist ( X d ) f ( m, ) Folie 88 Riccardo Gismondi Brownian Motion I Folie 89 Riccardo Gismondi Brownian Motion • One Dimension: • a standard one-dimensional Brownian motion on [0, T ] is a stochastic process {W (t ),0 t T } with the following properties: • W (0) 0 • t W (t ) is, with probability 1, a continuous function on [0, T ] • the increments {W (t1 ) W (t0 ), W (t 2 ) W (t1 ),..., W (t k ) W (t k 1 )} are independent for any k and any 0 t0 t1 ...t k T • W (t ) W ( s) ~ N (0, t s) for any 0 s t T Folie 90 Riccardo Gismondi Brownian Motion • One Dimension: • because of its mentioned properties, the following condition W (t ) ~ N (0, t ) is valid for 0 t T • for constants m and s 2 0 , the process X (t ) is a Brownian Motion with drift m and diffusion coefficient s 2 , if X (t ) mt is a Standard Brownian Motion s • X (t ) mt sW (t ) … represents a Brownian Motion with drift m and diffusion coefficient s 2 Folie 91 Riccardo Gismondi Brownian Motion • Standard Brownian Motion W (t ) Folie 92 Riccardo Gismondi Brownian Motion • Brownian Motion X (t ) Folie 93 Riccardo Gismondi Brownian Motion • One Dimension: • Consequently: • X solves the stochastic differential equation (SDE) dX (t ) mdt sdW (t ) • for deterministic, but time-varying m (t ) and s (t ) 0 dX (t ) m (t )dt s (t )dW (t ) • through integration t t X (t ) X (0) m ( s)ds s ( s)dW ( s) 0 0 • we come to the results. Folie 94 Riccardo Gismondi Brownian Motion • One Dimension: • The process X • has continuous sample paths and • independent increments • Each increment X (t ) X ( s) is normally distributed with • mean t E[ X (t ) X ( s )] m (u )du s • and variance t s (u )dW (u ) t s 2 (u )du Var[ X (t ) X ( s)] Var s s Folie 95 Riccardo Gismondi Brownian Motion • Random Walk Construction: • subsequent values for a standard Brownian motion from t 0 0 and W (0) 0 can be generated as follows: W (ti 1 ) W (ti ) ti 1 ti Zi 1 i 0,...,n 1 Z i ~ N (0,1) • for X ~ BM ( m , s 2 ) X (ti 1 ) X (ti ) m (ti 1 ti ) s ti 1 ti Zi 1 i 0,...,n 1 • with time-dependent coefficients ti1 ti1 X (ti 1 ) X (ti ) m ( s)ds s 2 (u )du Z i 1 i 0,...,n 1 ti ti • with Euler approximation X (ti 1 ) X (ti ) m (ti )(ti 1 ti ) s (ti ) ti 1 ti Zi 1 i 0,...,n 1 Folie 96 Riccardo Gismondi Home assignment (Matlab) • Write a code to generate a Standard Brownian Motion W (t ) • Write a code to generate a Brownian Motion X (t ) Folie 97 Riccardo Gismondi Brownian Motion • Random Walk Construction: • for a standard Brownian motion, we know that E W (ti ) 0 and for the covariance matrix we get CovW ( s), W (t ) CovW ( s), W ( s) (W (t ) W ( s)) CovW ( s), W ( s) CovW ( s), W (t ) W ( s) s0 s for the covariance matrix of W (t1 ),..., W (t n ) we denote: Cij min( ti , t j ) Folie 98 Riccardo Gismondi Brownian Motion • Cholesky Factorization • the vector W (t1 ),..., W (t n ) has the distribution N (0, C ) • and we may simulate this as vector AZ, where Z ( Z1 ,..., Z n )T ~ N (0, I ) • and A satisfies AAT C • by applying the Cholesky Factor, we get A as a lower triangular matrix of the form t1 0 ... 0 t1 t 2 t1 ... 0 A ... ... ... ... t1 t 2 t1 ... t n t n 1 Folie 99 Riccardo Gismondi Home assignment (Matlab) • Write a code to generate a Standard Brownian Motion W (t ) with Random Walk construction ( Cholesky factor ) Folie Riccardo Gismondi 100 Brownian Motion • Brownian Bridge Construction • besides generating the vector W (t1 ),..., W (t n ) from left to right, there also exists another method • e.g. • generating the last sample step first: W (t n ) • then generate W (t ) conditional on the value of W (t n ) n / 2 • and proceed with filling in the intermediate values • this method is useful for • variance reduction and • low-discrepancy methods • but it does not give us any computational advantage Folie Riccardo Gismondi 101 Brownian Motion • Brownian Bridge Construction • suppose 0u st • consider the problem, to generate W (s ) conditional on W (u ) x and W (t ) y • the unconditional distribution is given by W (u ) u u u W ( s ) ~ N 0, u s s W (t ) u s t Folie Riccardo Gismondi 102 Brownian Motion • Brownian Bridge Construction • first we permutate the entries W (s) s u s W (u ) ~ N 0, u u u W (t ) s u t • and by applying the conditioning formula to find the distribution of W (s ) conditional on the value of (W (u ),W (t )) we get the mean E[W ( s ) | W (u ) x,W (t ) y ] 1 u u x (t s ) x ( s u ) y u t y 0 (u s ) (t u ) Folie Riccardo Gismondi 103 Brownian Motion • Brownian Bridge Construction • and the variance 1 u u u ( s u )(t s) u t s s (u s) (t u ) • since 1 u u 1 t / u 1 u t t u 1 1 • the conditional mean is obtained by linearly interpolating between (u, x) and (t , y) Folie Riccardo Gismondi 104 Brownian Motion • Brownian Bridge Construction • suppose that W ( s1 ) x1 , W ( s2 ) x2 ,..., S ( sk ) xk • combining the observations of the mean and the variance (W ( s) | W ( s1 ) x1 , W ( s2 ) x2 ,...,W ( sk ) xk ) ( si 1 s) xi ( s si ) xi 1 ( si 1 s )(s si ) N , ( si 1 si ) ( si 1 si ) ... si s si 1 • and finally … ( si 1 s) xi ( s si ) xi 1 ( si 1 s)(s si ) W ( s) Z ( si 1 si ) ( si 1 si ) Folie Riccardo Gismondi 105 Brownian Motion • Brownian Bridge Construction • Brownian Bridge construction of Brownian path. Conditional on W ( s ) x and W ( s ) x , the value is normally distributed i i i 1 i 1 Folie Riccardo Gismondi 106 Home assignment (Matlab) • Write a code to generate a Standard Brownian Motion W (t ) with Brownian Bridge Construction Folie Riccardo Gismondi 107 Brownian Motion • Brownian Bridge Construction • How could the construction be modified for a Brownian motion with drift m • sampling the rightmost point would change • sample W (t m ) from N ( mt m , t m ) instead from N (0, t m ) Folie Riccardo Gismondi 108 Brownian Motion • Principal Components Construction • to visualize the construction of single Brownian path we can write in vector form: W (t1 ) a11 a12 a1n W (t2 ) a21 a22 a2 n ... ... Z1 ... Z 2 ... ... Z n W (t ) a a a n n1 n2 nn • with ai ( a1i ,... ani )T originally derived from AAT C through Cholesky-Factorzation • now, we want to derive a i through principal component construction, such that ai i vi • where 1 2 ... n 0 are the eigenvalues of C and vi are the eigenvectors Folie Riccardo Gismondi 109 Brownian Motion • Principal Components Construction • it can be shown that for an n-step path with equal spacing ti 1 ti t 2 2i 1 vi ( j ) sin j j 1 ,...,n 2n 1 2n 1 t 2i 1 i 1 ,...,n i sin 2 4 2n 1 2 • because of Cv v we can write min(t , t )v( j) v(i) j i j • for the discrete case 1 • and 0 min( s, t ) ( s )ds (t ) for the continuous limit with the eigenfunction on [0,1] Folie Riccardo Gismondi 110 Brownian Motion • Principal Components Construction • the solution to this equation is 2 (2i 1)t 2 i (t ) 2 sin i (2i 1) 2 • Karhounen-Loève expansion of Brownian Motion W (t ) i i (t ) Z i 0 t 1 i 0 Folie Riccardo Gismondi 111 Home assignment (Matlab) • Write a code to generate a Standard Brownian Motion W (t ) with Principal Components Construction (Karhounen-Loève) Folie Riccardo Gismondi 112 Geometric Brownian Motion, B&S Model, Pricing of (exotic) options: 1-dim Folie Riccardo Gismondi 113 Geometric Brownian Motion • a stochastic process S (t ) is a geometric Brownian Motion if log S (t ) is a Brownian motion with initial value log S (0) • in other words: a geometric Brownian motion is simply an exponentiated Brownian motion • all methods for simulating Brownian motion become methods for simulating geometric Brownian motion through exponentiation • a geometric Brownian motion is always positive, because the exponential function takes only positive values S (t 2 ) S (t1 ) S (t3 ) S (t 2 ) S (t n ) S (t n 1 ) , ,..., S (t1 ) S (t 2 ) S (t n 1 ) • are independent for t1 t 2 ...t n rather than the absolute changes S (ti 1 ) S (ti ) Folie Riccardo Gismondi 114 Geometric Brownian Motion • Basic properties • suppose W is a standard Brownian motion and X a Brownian Motion generated from W, such that X ~ BM ( m , s 2 ) • we set S (t ) S (0) exp( X (t )) f ( X (t )) • and with application of Itô’s formula 1 dS (t ) f ' ( X (t ))dX (t ) s 2 f ' ' ( X (t ))dt 2 1 S (0) exp( X (t ))[mdt s dW (t )] s 2 S (0) exp( X (t ))dt 2 1 S (t )(m s 2 )dt S (t )s dW (t ) 2 Folie Riccardo Gismondi 115 Geometric Brownian Motion • Basic properties • in contrast, a geometric Brownian motion is specified through the SDE ( B&S Model ) dS (t ) mdt s dW (t ) S (t ) • there seems to be an ambiguity in the role of m • the drift of a geometric Brownian Motion is mS (t ) and implies 1 2 d log S (t ) ( m s )dt s dW (t ) 2 as can be verified through Itô’s formula • we will use the notation S ~ GBM ( m , s 2 ) Folie Riccardo Gismondi 116 Geometric Brownian Motion • Basic properties • if S has initial value S (0), then 1 S (t ) S (0) exp([m s 2 ]t sW (t )) 2 • and respectively, if u t 1 S (t ) S (u ) exp([m s 2 ](t u ) s (W (t ) W (u ))) 2 • this provides a recursive procedure for simulating values of S at 0 t0 t1 ...t n 1 2 S (ti 1 ) S (ti ) exp([m s ](ti 1 ti ) s ti 1 ti Z i 1 ) i 0,1,...n 1 2 Folie Riccardo Gismondi 117 Home assignment (Matlab) • Write a code to generate a Geometric Brownian Motion S (t ) Folie Riccardo Gismondi 118 Geometric Brownian Motion • B&S Model • if S has initial value S (0), then 1 S (t ) S (0) exp([m s 2 ]t sW (t )) 2 • and respectively, if u t 1 S (t ) S (u ) exp([m s 2 ](t u ) s (W (t ) W (u ))) 2 • this provides a recursive procedure for simulating values of S at 0 t0 t1 ...t n 1 2 S (ti 1 ) S (ti ) exp([m s ](ti 1 ti ) s ti 1 ti Z i 1 ) i 0,1,...n 1 2 Folie Riccardo Gismondi 119 Geometric Brownian Motion • Lognormal Distribution • if S ~ GBM ( m , s ) , then the marginal distribution of 2 S (t ) is that of the exponential of a normal random variable, which is called a lognormal distribution • Y ~ LN ( m , s 2 ) , if the random variable Y has the distribution of exp(m sZ ), Z ~ N (0,1) • the distribution thus is given by log( y) m log( y) m P(Y y) P( Z ) ( ) s s • and the density 1 log( y ) m ( ) ys s Folie Riccardo Gismondi 120 Geometric Brownian Motion • Lognormal Distribution • for Y ~ LN ( m , s 2 ) 1 m s • Expected Value: E[Y ] e 2 2 m s 2 s2 • Variance: Var[Y ] e (e 1) 1 • if S ~ GBM ( m , s 2 ) then ( S (t ) / S (0)) ~ LN ([m s 2 ]t , s 2t ) and 2 E[ S (t )] e m t S (0) 2m t 2 s 2t Var[S (t )] e S (0)(e 1) Folie Riccardo Gismondi 121 Geometric Brownian Motion • Lognormal Distribution • in fact, we have E[ S (t ) | S ( ), 0 u ] E[ S (t ) | S (u )] e m (t u ) S (u ) , u t • the first equality is the Markov property (follows from the fact that S is a one-to-one transformation of a Brownian motion, itself a Markov process) • m acts as an average growth rate for S • for a standard Brownian motion W , we have t 1W (t ) 0 with probability 1 1 1 • for S ~ GBM ( m , s 2 ), we therefore find that log( S ) m s 2 with probability 1 t 2 Folie Riccardo Gismondi 122 Geometric Brownian Motion • Lognormal Distribution 1 2 • if the expression m s (growth rate) is positive 2 S (t ) as t • if it is negative S (t ) 0 as t Folie Riccardo Gismondi 123 Geometric Brownian Motion • Risk-Neutral Dynamics • the difficulty within this model lies in finding the correct probability measure for building the expectation and • for finding the appropriate discount rate • this bears on how the paths of the underlying are to be generated and on how the drift parameter m is chosen • we assume the existence of a continuous compounding interest rate r for riskless borrowing and lending • therefore a dollar investment grows as follows (t ) e r t at time t • we call the numeraire asset Folie Riccardo Gismondi 124 Geometric Brownian Motion • Risk-Neutral Dynamics • suppose the existence of an asset S that does not pay dividends • then under the risk-neutral measure, the discounted price process is a martingal S (u) S (t ) E | {S ( ),0 u} (u) (t ) • if S is a GBM under the risk neutral measure, then it must have m r dS (t ) • and further r dt s dW (t ) S (t ) • in a risk neutral world, all assets would have the same average rate of return Folie Riccardo Gismondi 125 Geometric Brownian Motion • Risk-Neutral Dynamics • therefore the drift parameter for S (t ) will equal the risk free rate r • should an asset pay dividends, then the martingal property continues to hold, but • with S replaced by the sum of S , any dividends paid and any interest earned from investing these dividends at r • let D(t ) be the value of all dividends paid over [0, t] • and assume that the asset pays a continuous dividend yield of such that it pays dividends of S (t ) at time t Folie Riccardo Gismondi 126 Geometric Brownian Motion • Risk-Neutral Dynamics • D grows at dD(t ) S (t ) rD (t ) dt • if S ~ GBM ( m , s 2 ) , then its drift in ( S (t ) D(t )) is mS (t ) S (t ) rD(t ) • now apply the martingal property to the combined process (S (t ) D(t )) , requires that this drift equal r ( S (t ) D(t )) • therefore, we must have m r , i.e. m r Folie Riccardo Gismondi 127 Geometric Brownian Motion • Risk-Neutral Dynamics • Example: Futures Contract • commits the holder to buying an underlying asset at a fixed price at a fixed date in the future • both, the buyer and the payer agree to the futures price, specified in the contract, without either party paying anything immediately to the other • S (T ) ... price of an asset at future time T t F (t , T ) ...futures price at time t for a contract o be settled at T t T Folie Riccardo Gismondi 128 Geometric Brownian Motion • Risk-Neutral Dynamics • for a futures contract with a zero value at the inception time t entails 0 e r (T t ) E[( S (T ) F (t , T )) | Ft ] • where Ft is the history of the market prices up to time t • at t T the spot and futures prices must agree, so S (T ) F (T , T ) • and we may rewrite the condition as F (t , T ) E[ F (T , T ) | Ft ] • the futures price is a martingale under the risk neutral measure • when modeling a futures price with a GBM, then you should set the drift parameter to zero: dF (t , T ) s dW (t ) F (t , T ) ( r )(T t ) • and finally this reveals to F (t , T ) e S (t ) Folie Riccardo Gismondi 129 Geometric Brownian Motion • B&S Model In Black-Scholes model many theoretical assumptions: 1. the short interest rate is known and is constant through time; 2. the underlying asset pays no dividends 3. the option is European 4. no transaction costs 5. is it possible to borrow any fraction of the price of a security to buy it or to hold it, at the short interest rate. 6. trading can be carried on continuously Folie Riccardo Gismondi 130 Geometric Brownian Motion • B&S Model With the assumption that the underlying follows a GBM model, S ~ GBM ( m , s 2 ) , and with no-arbitrage argument, Black and Scholes obtained the amazing formula for the European call option price : C S t , T , r , s , K S t N d1 Ke N d2 T t r S t 1 2 log r s T t K 2 d1 d2 s T t s T t S t 1 log r s 2 T t K 2 d2 Folie s T t Riccardo Gismondi 131 Home assignment (Matlab) • Write a code for pricing a European call/put option simulating S (t ) as S ~ GBM ( m , s 2 ) and study the asymptotic convergence to the theoretical B&S price. • Describe and plot the convergence of the MC estimator. • What is the variance of the estimator? • When will it slow down ? • How many simulations you need to be “sure”? • Write a report about your research. Folie Riccardo Gismondi 132 Geometric Brownian Motion • Path-Dependent Options • the reason behind simulating GBM paths, is that this can be used for pricing options, (e.g. not only Vanillas) • particularly for those whose payoff depend on the path of an underlying asset S and not simply its value S (T ) • the price of an option may be represented as an expected discounted payoff, so • simulate by generating paths of the underlying asset, • evaluate the discounted payoff on each path • and average over all paths Folie Riccardo Gismondi 133 Geometric Brownian Motion • Path-Dependent Payoffs • an option price could in principal depend on the complete path over an interval [0, T ] • e.g. • Asian options – discrete monitoring: is an option on a time average of the underlying asset. Payoff Asian Call: S K Payoff Asian Put: K S with constant strike price K 1 n and S S (ti ) n i 1 Folie Riccardo Gismondi 134 Geometric Brownian Motion • Path-Dependent Payoffs • Asian options – continuous monitoring: just replace the discrete average with the continuous 1 t S t u u S ( )d over an integral [u, t ] Folie Riccardo Gismondi 135 Geometric Brownian Motion • Path-Dependent Payoffs • Geometric average option: replacing the arithmetic average S with 1/ n n S (ti ) i 1 • such options are useful in test cases for computational procedures • they are mathematically convenient to work with because of its geometric average • we find (with m replaced by r ), that 1/ n n 1 2 1 n s n S (ti ) S (0) exp [r s ] ti W (ti ) i 1 2 n i 1 n i 1 Folie Riccardo Gismondi 136 Home assignment (Matlab) • Write a code to price a Asian Option with arithmetic/geometric mean and discrete monitoring : function [Call , Put ] AsianOption S0 , K , T , r , s , ts, MeanType S0 ........................Value of S in t 0 Startvalue K ........................Strike T .........................Espiration, e. g. 1 year r..........................Risk free rate s ..........................Volatility 1 ts..........................timestep, e. g. daily, 252 Aritmetical MeanType............ Geometrical Folie Riccardo Gismondi 137 Geometric Brownian Motion • Path-Dependent Payoffs • Barrier Options: the option gets “knocked out”, if the underlying asset crosses a prespecified level. e.g. a “down-and-out call” with • barrier b • strike K • and expiration T • has the payoff 1{ (b) T }( S (T ) K ) • where (b) inf{ti : S (ti ) b} is the first time that the price of the underlying drops below b • and 1{ } represents the indicator function of the event in the braces Folie Riccardo Gismondi 138 Geometric Brownian Motion • Path-Dependent Payoffs • Barrier Options: respectively, a “down-and-in call” has the payoff 1{ (b) T }( S (T ) K ) and gets “knocked-in” only when the underlying asset crosses the barrier. • Up-and-out and up-and-in calls and puts can be derived analogously! Folie Riccardo Gismondi 139 Home assignment (Matlab) • Write a code to price a Barrier Option (down and in, down and out) with discrete monitoring : function [Call , Put ] BarrierOption S 0 , K , T , r , s , b, ts , Type S0 ........................Value of S in t 0 Startvalue K ........................Strike T .........................Espiration, e.g . 1 year r..........................Risk free rate s ..........................Volatility b...........................Barrier 1 ts..........................timestep , e.g . daily , 252 D _ In Type.................... D _ Out Folie Riccardo Gismondi 140 Geometric Brownian Motion • Path-Dependent Payoffs • Lookback options: like barrier options, lookback options depend on extremal values of the underlying asset price. Such puts and calls expiring at t n have payoffs: ( max S (ti ) S (tn )) and (S (tn ) min S (ti )) i 1,...,n i 1,...,n • e.g. a lookback call can be viewed as • the profit from buying at the lowest price over t1 ,..., t n • and selling at the final price S (t n ) Folie Riccardo Gismondi 141 Geometric Brownian Motion • Incorporating a Term Structure • till now, we assumed that risk free rate r is constant • therefore the price of (riskless) zero-coupon bond maturing at T t can be calculated as follows: B(t , T ) e r (T t ) • now we assume, that on the markets we observe bond prices B(0, T ) that are incompatible with the above mentioned pricing formula • therefore, to use this term structure in option-pricing, we have to introduce a deterministic, but time-varying risk-free rate r (u ) r (u ) log B (0, T ) T T u Folie Riccardo Gismondi 142 Geometric Brownian Motion • Incorporating a Term Structure T r (u )du) 0 • then we get B(0, T ) exp • with this time-varying risk-free rate r (u ) , the dynamics of an asset price S (t ) under the risk-neutral measure can be described by the SDE (stochastic differential equation) dS (t ) r (t )dt s dW (t ) S (t ) t • with solution S (t ) S (0) exp r (u ) du s 2t s W (t ) 1 0 2 • and this process can be simulated over 0 t0 t1 ... t n ti1 1 S (ti 1 ) S (ti ) exp r (u ) du s 2 (ti 1 ti ) s ti 1 ti Z i 1 ti 2 Folie Riccardo Gismondi 143 Home assignment (Matlab) • Write a code to price a Asian Option with discrete monitoring and time-varying interest rate: function [Call , Put ] AsianOption _ IntVar S0 , K , T , r t , s , b, ts, MeanType S0 ........................Value of S in t 0 Startvalue K ........................Strike T .........................Espiration, e. g . 1 year r t .....................Time dipendent Risk free rate e.g in line function s ..........................Volatility b...........................Barrier 1 ts..........................timestep, e.g. daily, 252 Aritmetical MeanType............ Geometrical e.g., you can take : r t 1 100 sin 2 t t 3 , 0 t T Folie Riccardo Gismondi 144 Geometric Brownian Motion • Incorporating a Term Structure • if we observe bond prices B(0, ti ),..., B(0, t n ) • we can calculate the term structure as follows B (0, ti ) ti1 r (u ) du exp B (0, ti 1 ) ti • and we can simulate S (t ) using B (0, ti ) 1 2 S (ti 1 ) S (ti ) exp s (ti 1 ti ) s ti 1 ti Z i 1 B (0, ti 1 ) 2 Folie Riccardo Gismondi 145 Geometric Brownian Motion • Simulating Off a Forward Curve • now, we not just assume the observation of S (0) • but also a series of forward prices F (0, T ) for the asset • under the risk-neutral measure, F (0, T ) E[S (T )] • this implies 1 2 S (t ) F (0, T ) exp - s T s W(t ) 2 • and we can simulate F (0, ti 1 ) 1 S (ti 1 ) S (ti ) exp s 2 (ti 1 ti ) s ti 1 ti Z i 1 F (0, ti ) 2 Folie Riccardo Gismondi 146 Geometric Brownian Motion • Deterministic Volatility Functions • on the markets, it has been widely observed, that option prices are incompatible with a GBM model for the underlying asset • applying the Back-Scholes model for option prices, this would assume to use the same volatility s for all traded options with • different strikes, • different maturities that • are traded on the same underlying asset. • in practice, the observed volatility of these traded assets varies with • strike and • maturity Folie Riccardo Gismondi 147 Geometric Brownian Motion • Deterministic Volatility Functions • Black-Scholes has to be modified to find the real market prices! • consequently, we let our volatility depend on S and t , such that it looks like the following s ( S , t ) • this leads to the model dS (t ) rdt s ( S (t ), t )dW (t ) S (t ) • assuming a deterministic volatility function, an option could be hedged through a position in the underlying asset • this would not be the case in a stochastic volatility model Folie Riccardo Gismondi 148 Geometric Brownian Motion • Deterministic Volatility Functions • to get s ( S , t ) a numerical optimization problem has to be solved • in general, there is no exact simulation procedure for these models and it is necessary to use an Euler scheme of the form S (ti 1 ) S (ti ) 1 r (ti 1 ti ) s (S (ti ), ti ) ti 1 ti Zi 1 • or the Euler scheme for log( S (t )) S (ti 1 ) S (ti ) exp [r 1 s 2 (S (ti ), ti )](ti 1 ti ) s (S (ti ), ti ) ti 1 ti Zi 1 2 Folie Riccardo Gismondi 149 Home assignment (Matlab) • Write a code to price a Asian Option with arithmetic/geometric mean, discrete monitoring and time-varying volatility : function [Call , Put ] AsianOption _ VolVar S 0 , K , T , r , s S t , t , ts, MeanType S0 ........................Value of S in t 0 Startvalue K ........................Strike T .........................Espiration, e.g . 1 year r..........................Risk free rate s S t , t ...........Volatility function In line function 1 ts..........................timestep, e.g. daily, 252 Aritmetical MeanType............ Geometrical Folie Riccardo Gismondi 150 Brownian Motion II Folie Riccardo Gismondi 151 Brownian Motion • Multiple Dimensions • a process W (t ) (W1 (t ),..., Wd (t )) , 0 t T T is called a standard Brownian Motion on d , if it has • W (0) 0 • continuous sample paths • independent increments and • W (t ) W (s) ~ N (0, (t s) I ) for all 0 s t T • I is the d d identity matrix • it follows that each process Wi (t ), i 1,..., d is a standard one- dimensional Brownian motion and that Wi and W j are independent for i j Folie Riccardo Gismondi 152 Brownian Motion • Multiple Dimensions • suppose: m is a vector in d and is a d d matrix, positive definite or semidefinite • X is Brownian Motion with drift m and covariance X ~ BM (m , ) • X (t ) X (s) ~ N ((t s)m , (t s)) • X (t ) mt BW (t ) …. whereby BBT • this process solves the SDE (stochastic differential equation) dX (t ) mdt BdW (t ) • for deterministic, but time-varying m (t ) and (t ): dX (t ) m (t )dt B(t )dW (t ) … whereby B (t ) B T (t ) (t ) Folie Riccardo Gismondi 153 Brownian Motion • Multiple Dimensions • this process has continuous sample paths, independent increments and t m (u)du, t (u)du X (t ) X (s) ~ N s s • Cov X i ( s ), X j (t ) min( s, t ) ij Folie Riccardo Gismondi 154 Brownian Motion • Random Walk Construction • let Z1 , Z 2 ,... be independent N (0, I ) random vectors in d • we can construct a standard d-dimensional Brownian motion at times 0 t0 t1 ... t n by setting W (0) 0 and W (ti 1 ) W (ti ) ti 1 ti Zi 1 i 0,...,n 1 • to simulate X ~ BM (m , ) we first have to find a matrix B for which BBT • set X (0) 0 and X (ti 1 ) X (ti ) m (ti 1 ti ) ti 1 ti BZi i 0,...n 1 Folie Riccardo Gismondi 155 Brownian Motion • Random Walk Construction • thus, simulation of BM ( m , ) is straightforward once has been factored • for the case of time-dependent coefficients, we may set ti1 X (ti 1 ) X (ti ) m (s)ds B(ti , ti 1 )Zi i 0,...n 1 ti with ti1 B(ti , ti 1 ) B(ti , ti 1 ) (u)du T ti Folie Riccardo Gismondi 156 Home assignment (Matlab) • Write a code to generate a d-dimensional Standard Brownian Motion W (t ) • Write a code to generate a d-dimensional Brownian Motion X (t ) Folie Riccardo Gismondi 157 Brownian Motion • Brownian Bridge Construction • we can apply independent one-dimensional constructions also for the multi-dimensional case • to include a drift vector, it suffices to add m i t n to Wi (t n ) at the first step of the construction of the ith coordinate • let W be a standard k-dimensional Brownian motion • and B a d k matrix, with k d , • satisfying BBT • this results into X (t ) m (t ) BW (t ) as X ~ BM (m , ) and X (t1 ),..., X (t n ) are recovered through a linear transformation Folie Riccardo Gismondi 158 Brownian Motion • Principal Components Construction • by application of the principal components construction to the multidimensional case • the covariance matrix can be represented as C11 C12 ... C1n C21 C22 ... C2 n (C ) ... ... ... ... C C ... C n1 n2 nn the eigenvectors of this matrix can be written as (vi w j ) where vi are the eigenvectors from C w j are the eigenvectors from and the eigenvalues i j where i are the eigenvalues from C j are the eigenvalues from Folie Riccardo Gismondi 159 Brownian Motion • Principal Components Construction • ranking the products of the eigenvalues (i j ) (1) (i j ) ( 2 ) ...( i j ) ( nd ) • then for any k 1,...,n r 1 (i j )(r ) i k k i 1 (i j ) ( r ) i nd n r 1 i 1 • the variability of the first k factors is always smaller for a d- dimensional Brownian motion than for a scalar Brownian motion over the same time points • since the d-dimensional process has greater total variability Folie Riccardo Gismondi 160 Geometric Brownian Motion, Pricing of (exotic) options : n-dim Folie Riccardo Gismondi 161 Geometric Brownian Motion • Multiple Dimensions • a multidimensional geometric Brownian motion can be written by the SSDE (system of stochastic differential equations) dSi (t ) mi dt s i dX i (t ) i 1,...,d Si (t ) • where each X i is a standard one-dimensional Brownian motion • and X i and X j have the correlation ij Folie Riccardo Gismondi 162 Geometric Brownian Motion • Multiple Dimensions • by setting ij s is j ij (variance-covariance-matrix) • then (s 1 X 1 ,..., s d X d ) ~ BM (0, ) • and S ( S1 ,..., S d ) ~ GBM ( m , ) with m ( m1 ,..., m d ) • the actual drift vector is given by ( m1S1 (t ),..., m d S d (t )) • and the covariances are given by ( mi m j ) t ijs is j Cov[Si (t ), S j (t )] Si (0)S j (0)e (e 1) • with 1 ( m i s i2 ) t s i X i ( t ) i 1,...,d S i (t ) S i (0)e 2 Folie Riccardo Gismondi 163 Geometric Brownian Motion • Multiple Dimensions • with cholesky factorization from we get a matrix A • such that AAT d dSi (t ) • and we can formulate mi dt Aij dW j (t ) Si (t ) j 1 • this leads to an algorithm for simulating GBM ( m , ) at the times 0 t0 t1 ... t n 1 ( m i s i2 )(t k 1 t k ) t k 1 t k d1 Aij Z k 1, j S i (t k 1 ) S i (t k )e j 2 Folie Riccardo Gismondi 164 Home assignment (Matlab) • Write a code to generate a d-dimensional Geometric Brownian Motion S (t ) Folie Riccardo Gismondi 165 Geometric Brownian Motion • Multiple Dimensions • options depending on multiple assets • e.g. • Spread-Option a call option on the spread between two assets S1 and S 2 has the payoff ([ S1 (T ) S 2 (T )] K ) • Basket Option an option on a portfolio of underlying assets and has the payoff ([ c1S1 (T ) c2 S 2 (T ) ... cd S d (T )] K ) Folie Riccardo Gismondi 166 Home assignment (Matlab) • Write a code to price a Basket Option and discrete monitoring: function [Call , Put ] BasketOption S 1 , S 2 ,.., S d , c1 , c 2 ,.., c d , K , T , r , ts S 1 , S 2 ,.., S d ........................Time Series Matrix of Underlyings c1 , c 2 ,.., c d ..........................Weights Vector K ..........................................Strike T ..........................................Espiration, e.g . 1 year r...........................................Risk free rate vector 1 ts..........................................timestep, e.g. daily, 252 Folie Riccardo Gismondi 167 Geometric Brownian Motion • Multiple Dimensions • Outperformance Option two options on the maximum or minimum of multiple assets; e.g. the payoff can be (max{ c1S1 (T ), c2 S 2 (T ),..., cd S d (T )} K ) • Barrier Options e.g. a two asset barrier – here a down-and-in put – with payoff 1{ min S 2 (ti ) b}(K S1 (T )) i 1,...,n this is a down-and-in put on S1 that knocks in when S 2 drops below the barrier b Folie Riccardo Gismondi 168 Home assignment (Matlab) • Write a code to price a Outperformance Option and discrete monitoring: function [Call , Put ] OutPerfOption S 1 , S 2 ,.., S d , c1 , c 2 ,.., c d , K , T , r , ts S 1 , S 2 ,.., S d ........................Time Series Matrix of Underlyings c1 , c 2 ,.., c d ..........................Weights Vector K ..........................................Strike T ..........................................Espiration, e.g . 1 year r...........................................Risk free rate vector 1 ts..........................................timestep, e.g. daily, 252 Folie Riccardo Gismondi 169 Geometric Brownian Motion • Multiple Dimensions • Quantos options that are sensitive to a stock price and an exchange rate (e.g. different currencies for the option and the stock) with S1 as the stock price and S 2 as the exchange rate (expressed as the quantity of domestic currency required per unit of the foreign currency) the payoff of the option in the domestic currency is given by S 2 (T )( S1 (T ) K ) Folie Riccardo Gismondi 170 Geometric Brownian Motion • Multiple Dimensions • Quantos (cont) respectively, the payoff K S1 (T ) S 2 (T ) corresponds to Quanto, whereby the strike is fixed in the domestic currency and the payoff of the option is made in foreign currency Folie Riccardo Gismondi 171 Home assignment (Matlab) • Write a code to price a Quanto Option and discrete monitoring : function [Call , Put ] QuantosOption S , Currency d / Currency f , K , T , r , ts S .........................................................Time Series of Underlying Currency d / Currency f .......................Time Series of Exchange Rate K ........................................................Strike T .........................................................Espiration, e.g. 1 year r..........................................................Risk free rate vector 1 ts.........................................................timestep, e.g. daily, 252 Folie Riccardo Gismondi 172 Square root diffusion process: CIR Model Folie Riccardo Gismondi 173 Square-Root Diffusions • Basic properties • Consider a class of process that includes the square-root diffusion dr (t ) a (b r (t ))dt s r (t )dW (t ), with W a standard one-dimension Brownian motion. • We consider the case in which: a 0, b 0 • If r (0) 0, then r (t ) 0 for all t • If 2ab s 2 , then r (t ) remains strictly positive for all t • All of the coefficients could in principle be time-dependent • i.e , with b b(t ) dr (t ) a (b(t ) r (t ))dt s r (t )dW (t ) Folie Riccardo Gismondi 174 Square-Root Diffusions • Applications • CIR Model • Heston Model Folie Riccardo Gismondi 175 Square-Root Diffusions • Applications • CIR Model • The square-root diffusion process dr (t ) a (b r (t ))dt s r (t )dW (t ), can be used as a model of the short rate. • This model was done to illustrate the workings of a general equilibrium model and was proposed by Cox-Ingersoll-Ross (CIR) (1985). Folie Riccardo Gismondi 176 Square-Root Diffusions • Applications • Heston Model • Heston proposed a stochastic volatility model in which the price of an asset S (t ) is governed by dS (t ) mdt V (t ) dW1 (t ) S (t ) the squared volatility V (t ) follows a square-root diffusion dV (t ) a (b V (t ))dt s V (t )dW2 (t ) and (W1 ,W2 ) is a two-dimensional Brownian motion with corr (dW1 , dW2 ) dt Folie Riccardo Gismondi 177 Square-Root Diffusions • Simulating r (t ) with Euler Method • Discretization at times t1 , t 2 ,..., t n by setting: r (ti 1 ) r (ti ) a (b r (ti ))[ti 1 ti ] s r (ti ) ti 1 ti Zi 1 with Z1 , Z 2 ,..., Z n independent N (0,1) random variables • Notice: • we have taken the positive part of r (ti ) inside the square root. Some modification of this form is necessary because the values of r (ti ) produces by Euler discretization may become negative. Folie Riccardo Gismondi 178 Square-Root Diffusions • Simulating r (t ) with Transition Density Method • A noncentral chi-square random variable ( ) with '2 degrees of freedom and noncentrality parameter has distribution P( '2 ( ) y) F ' 2 ( ) ( y) e / 2 0.5 j / j! y j 0 2 2( / 2) j ( j ) 0 z ( / 2) j 1e z / 2 dz, for y 0 • The transition law of r (t ) can be expressed as s 2 (1 e a (t u ) ) '2 4ae a ( t u ) r (t ) d ( 2 a ( t u ) r (u )), t u 4a s (1 e ) 4ba where d 2 s Folie Riccardo Gismondi 179 Square-Root Diffusions • Simulating r (t ) with Transition Density Method • This says that, given r (u ), r (t ) is distributed as s 2 (1 e a (t u ) ) /(4a ) times a noncentral chi-square random variable with d degrees of freedom and noncentraly parameter a ( t u ) 4a e 2 a ( t u ) r (u ); s (1 e ) • Equivalently 4a y P(r (t ) y | r (u)) F ' 2 ( ) 2 a ( t u ) , d s (1 e ) Folie Riccardo Gismondi 180 Square-Root Diffusions • Simulating r (t ) with Transition Density Method • Chi-Square • If is a positive integer and Z1 , Z 2 ,..., Z are independent N (0,1) random variables, then the distribution of Z12 Z12 ... Z2 is called the chi-square distribution with degrees of freedom. • The chi-square distribution is given by y 1 2 ( / 2) P( y ) ( / 2) 2 e z / 2 z ( / 2)1dz, 0 where (.)denotes the gamma function and (n) (n 1)! for n Folie Riccardo Gismondi 181 Square-Root Diffusions • Simulating r (t ) with Transition Density Method • Noncentral Chi-Square • For integer and constants a1 , a2 ,..., a , the distribution of (Z i 1 i ai ) 2 is noncentral chi-square with degrees of freedom and noncentrality parameter a12 a2 ... a2 . 2 • if 1, then '2 ( ) 1'2 ( ) 21 • Consequently '2 ( ) (Z )2 21 Folie Riccardo Gismondi 182 Square-Root Diffusions • Simulating r (t ) with Transition Density Method • If N is a Poisson random variable with mean / 2, than ( / 2) j P( N j ) e / 2 , j 0,1,2,... j! • Conditional on N j has the random variable 2 N the ordinary chi- 2 square distribution with 2 j degrees of freedom: y 1 (( / 2) j ) P( 2 2 N y | N j ) ( / 2) j e z / 2 z ( / 2) j 1dz 2 0 The unconditional distribution is thus given by ( / 2) j P( N j)P( j 0 2 2 N y | N j) e j 0 / 2 j! P( 2 2 j y) Folie Riccardo Gismondi 183 Square-Root Diffusions • Simulating r (t ) with Transition Density Method • Simulation of square-root diffusion dr(t ) a (b r (t ))dt s r (t )dW (t ), on time gird 0 t0 t1 ... t n with d 4ba / s 2 1 by sampling from the transition density Folie Riccardo Gismondi 184 Square-Root Diffusions • Simulating r (t ) with Transition Density Method • Simulation of square-root diffusion dr(t ) a (b r (t ))dt s r (t )dW (t ), on time gird 0 t0 t1 ... t n with d 4ba / s 2 1 by sampling from the transition density Folie Riccardo Gismondi 185 Square-Root Diffusions • Simulating r (t ) with Transition Density Method • Comparison of exact distribution (solid) and one-step Euler approximation (dashed) for a square-root diffusion: a 0.2, s 0.1, b 5%, r (0) 4% t 0.25 t 1 Folie Riccardo Gismondi 186 Square-Root Diffusions • Sampling Gamma and Poisson • Gamma Distribution • The gamma distribution with shape parameter a and scale parameter has the density 1 f ( y) f a, ( y) y a 1e y / , y0 (a ) a with mean a and variance a 2 • The chi-square distribution is a special case of gamma distribution with scale parameter 2 Folie Riccardo Gismondi 187 Square-Root Diffusions • Sampling Gamma and Poisson • Gamma Distribution • Algorithms GKM1 for sampling from the gamma distribution with parameters (a,1), a 1 Folie Riccardo Gismondi 188 Square-Root Diffusions • Sampling Gamma and Poisson • Gamma Distribution • Algorithms (Ahrens-Dieter) for sampling from the gamma distribution with parameters (a,1), a 1 Folie Riccardo Gismondi 189 Square-Root Diffusions • Sampling Gamma and Poisson • Poisson Distribution • The poisson distribution with mean 0 is given by k P( N k ) e , k 0,1,2,... k! • We write N ~ Poisson( ) • Poisson is the distribution of the number of events in [0,1] when the times between consecutive events are independent and exponentially distributed with mean 1 / . Folie Riccardo Gismondi 190 Square-Root Diffusions • Sampling Gamma and Poisson • Poisson Distribution • Inverse transformations method for sampling from the Poisson distribution with mean 0 Folie Riccardo Gismondi 191 Heston Model and stochastic Volatility Folie Riccardo Gismondi 192 The Heston Model • Is a stochastic volatility model in which the price of an asset S (t ) is governed by dS (t ) m dt V (t )dW1 (t ) S (t ) the squared volatility V (t ) follows a square-root diffusion dV (t ) a (b V (t ))dt s V (t )dW2 (t ) and (W1 ,W2 ) is a two-dimensional Brownian motion with corr (dW1 , dW2 ) dt Folie Riccardo Gismondi 193 The Heston Model • Equations description • The first equation gives the dynamic of a stock price • The parameters represents • S (t ) the stock price at time t • m the risk neutral drift • V(t) the volatility Folie Riccardo Gismondi 194 The Heston Model • Equations description • The second equation gives the evolution of the variance which follows the square-root process • The parameters represents • b the long vol, or long run average price volatility; as t tends to infinity, the expected value of V (t ) tends to b. • a the mean reversion parameter; rate at which V (t ) reverts to b • s the vol of vol, or volatility of the volatility; this determines the variance of V (t ). Folie Riccardo Gismondi 195 The Heston Model • Simulating S (t ) • 1. Simulate W (t ) (W1 (t ), W2 (t )) ( e.g. with random walk construction ) by setting W (ti 1 ) W (ti ) ti 1 ti BZi 1 with B t B , 1 and 1 Z1 , Z 2 ,... independent standard random normal vectors Folie Riccardo Gismondi 196 The Heston Model Simulating S (t ) • 2. Simulate V (t ) ( e.g. with Euler discretization ) at times t1 , t2 ,..., tn by setting V (ti 1 ) V (ti ) a (b V (ti ))[ti 1 ti ] s V (ti ) (W2 (ti 1 ) W2 (ti )) • 3. Simulate S (t ) at times t1 , t2 ,..., tn by setting S (ti 1 ) S (ti ) m S (ti )[ti 1 ti ] V (ti ) S (ti )(W1 (ti 1 ) W1 (ti )) Folie Riccardo Gismondi 197 The Heston Model • Option pricing • Call Option • The time t price of a European call option with time to maturity (T t ) denoted Call ( S , v, t ), is given by Call ( S , v, t ) S (t ) P KP(t , T ) P2 1 • Put Option • The price of a European put option at time t is obtained through put-call parity and is given by Put (S , v, t ) Call ( S , v, t ) KP(t, T ) S (t ) or Put ( S , v, t ) ( P 1) S (t ) (1 P2 ) KP(t , T ) 1 Folie Riccardo Gismondi 198 The Heston Model • Option pricing • P , P2 are given by 1 1 1 eiu ln a j ( S0 ,V0 , t , T , u ) Pj Re du, j 1, 2. 2 0 iu • P(t , T ) is given by P(t , T ) e ( r q )(T t ) Folie Riccardo Gismondi 199 The Heston Model Option pricing • with j ( S0 ,V0 , ; ) exp{C j (r; ) D j (r; )V0 i S0 }, T t. dC j ( ; ) a bD j ( ; ) (r q) i 0, d dD j ( ; ) s 2 D 2 ( ; ) 2 (b j s i ) D j ( ; ) m j i 0. j d 2 2 C j (0; ) D j (0; ) 0. Folie Riccardo Gismondi 200 The Heston Model • Option pricing • and solutions ab 1 ged C ( ; ) (r q) i 2 (b j s i d ) 2 ln s 1 g b j si d 1 ed D( ; ) s2 1 ged where b j si d g , d ( si b j )2 s 2 (2u ji 2 ) b j si d u1 0.5, u2 0.5, b1 a s , b2 a Folie Riccardo Gismondi 201 The Heston Model • Option pricing • The parameters represents the following • S (t ) the spot price of the asset • K the strike price of the option • r interest rate • q dividend yield • P(t , T ) a discount factor from time t to time T . • P1 and P2 are the probabilities that the call option expires in- the-money. Folie Riccardo Gismondi 202 The Heston Model • Exact Simulation • The stock price at time t given the values of S (u ) and V (u ) for u t can be written as: 1 t t S (t ) S (u ) exp m (t u ) V ( s)ds V ( s)dW2 ( s) 2u u • were the variance at time t is given by the equation: t t V (t ) V (u ) a b(t u ) a V ( s)ds s V ( s)dW1 ( s) u u Folie Riccardo Gismondi 203 The Heston Model • Exact Simulation (Algorithm) • Input: S (u ),V (u ) Output: S (t ) with u t • 1. Generate a sample from the distribution of V (t ) given V (u ) t • 2. Generate a sample from the distribution of V ( s )ds given V (t ) and V (u ). u t t • 3. Recover u V ( s)dW1 ( s) given V (t ),V (u ) and V ( s )ds • 4. Generate a sample from the distribution of S(t) u given t t V ( s)dW1 ( s) and V (s)ds u u Folie Riccardo Gismondi 204 Pricing of Exotic Options: Structured Products (Equity) Folie Riccardo Gismondi 205