编程代写 Lectures 6 and 7: Simulation

Lectures 6 and 7: Simulation
Computational Finance

Building Models

A model is an approximate mathematical description of real-world
E.g. model describing gravity
122 • 𝐹𝐹=𝐺𝐺𝑚𝑚 𝑚𝑚 and𝐹𝐹=𝑚𝑚𝑚𝑚imply𝑚𝑚=𝐺𝐺𝑚𝑚
• The acceleration of an object towards another object (e.g. Earth) is proportional to the mass of the other object and inversely proportional to the square of the distance between them
• Allows you to predict how fast a ball will fall to the ground on Earth vs. on the moon E.g. model describing stock prices (CAPM)
•𝐸𝐸𝑅𝑅−𝑅𝑅=𝛽𝛽𝑅𝑅−𝑅𝑅 𝑖𝑖𝑓𝑓𝑖𝑖𝑚𝑚𝑓𝑓
• The excess return on a stock is proportional to the return on the market
• Allows you to predict a particular stock’s return given the market’s
• Much less accurate than the gravity model! Human behavior is hard to model…

Solving Models
Up to now: solve the model by solving an equation, or a system of equations, describing what happens
• Sometimes with pen and paper, sometimes on a computer (e.g. tangency portfolio in HW 1)
Many models cannot be solved this way
Simulation: create many different representative scenarios of what could happen
• In finance, Monte Carlo simulation: sample from a known distribution of possible scenarios, calculate the outcome for that scenario, to get distribution of outcomes
• Useful when the distribution is complex and high-dimensional • Also used in physics and engineering

Examples of Simulations Outside of Finance

Problem: find E[g(X)]
Example: if X is a standard normally distributed random variable, what is the expected value of
𝑒𝑒𝑋𝑋 (1+𝑋𝑋2)
• Complicated function
• What is an expected value? Integral with
E[𝑌𝑌] = �∞ 𝑒𝑒𝑥𝑥 −∞ (1 + 𝑥𝑥2)
∗ 𝑒𝑒−𝑥𝑥2 𝑑𝑑𝑥𝑥
respect to the probability density:
• Didn’t get any easier…
• Need a computer to solve

Solving for the expectation numerically
Limit your range to likely values of X (where the integrand is different enough from zero)
Divide the x axis into n equal segments 𝑚𝑚𝑚𝑚𝑚𝑚𝑛𝑛 𝑚𝑚𝑚𝑚𝑚𝑚
• Length: 𝑑𝑑 = 𝑥𝑥 −𝑥𝑥
• Midpoints 𝑥𝑥𝑖𝑖 = 𝑥𝑥𝑚𝑚𝑖𝑖𝑛𝑛 + 𝑑𝑑 𝑖𝑖 Calculate ∑𝑖𝑖 𝑔𝑔 𝑥𝑥𝑖𝑖 𝜙𝜙 𝑥𝑥𝑖𝑖 2to get
your approximation of E[g(X)]

𝑔𝑔𝑋𝑋,𝑋𝑋 12
What if is a function of two variables
Example: if 𝑋𝑋1and 𝑋𝑋2 are jointly standard normally distributed random variable, what is the expected value of
𝑌𝑌=𝑔𝑔𝑋𝑋,𝑋𝑋 =
1 2 ( 1 + 𝑋𝑋 12 + | 𝑋𝑋 2 | ) 2 2
E[𝑌𝑌] = ∫∫∞ 𝑒𝑒𝑋𝑋1+𝑋𝑋2 ∗ 𝑒𝑒− 𝑚𝑚1+𝑚𝑚2 𝑑𝑑𝑥𝑥 𝑑𝑑𝑥𝑥 − ∞ ( 1 + 𝑋𝑋 12 + | 𝑋𝑋 2 | ) 2 𝜋𝜋 1 2
where the second term is the PDF of the multivariate standard normal distribution

Solving for the expectation numerically
Partition both axes into equal intervals forming squares
Compute the value of the function for the center point of each square
Add up the volume of each bar to approximate the volume under the curve i.e. the value of the double integral
Already getting hard
What if you had a function of 200 random variables e.g. value of a portfolio of 200 stocks?
Want to solve a 200-dimensional integral?

Monte Carlo Simulation is Easier and Faster
• Sample random variables from their joint distribution using known methods of generating random numbers
• Evaluate the function at each sample and take the mean • MC_integration.m
• By sampling from the underlying distribution of X, the Monte Carlo method evaluates the function at values that are most likely to arise
• Avoids computing the function for values that will rarely occur
Valuable when
• Computing the function is costly
• The function is of many variables
• Both situations arise in finance all the time!

Generating
Sampling from a univariate distribution

Why generate random numbers?
Because we want unpredictability
• Example: random number inside a slot machine in a casino • Not very relevant for this class
Simulation: because we want a distribution
• Each random number represents a scenario
• Example: surface temperatures in an area of the Pacific Ocean where typhoons may form
• Example: prices of stocks in our portfolio on a future day

Generating Random Numbers
Computers are deterministic machines, so can’t generate truly random numbers
Programming languages generate pseudorandom numbers
• These are numbers that appear and behave random enough that we can treat them as if they were random
• With a pseudorandom number generator (PRNG), once a seed value is passed into the algorithm, the sequence of numbers is deterministic
• seeAppendix:PseudorandomNumberGeneratorsforhowthisworks • We can set the seed value for MATLAB’s PRNG using the rng()
• rng(n) sets the seed to the integer n
• rng(‘shuffle’) uses the current date and time to set the seed
• rng(‘default’) sets the seed to the value used when MATLAB starts

Pseudo- vs. quasi-random numbers
Do we want unpredictability, or do we want a distribution?
• Differentrandomnumbergeneratorsarebetterforoneversustheother
Pseudo-random number algorithms are designed to be as unpredictable as possible
• Goodifyou’rerunningacasino
• Butmayproduce”clustering”(bunching)whenmultiplenumbersaresampled • Requireverylargesamplestoaccuratelymatchadesireddistribution
• Matlab:rand(),randn()andmostotherfunctionswewilltalkabouttoday
Quasi-random numbers are designed to “fill in” a particular distribution
• Knowingthefirstn-1generatednumbershelpspredictthen’thnumbernotgoodfor casinos at all!
• Butcanmatchadesireddistributionquicker(willbeusefulforusstartingnextweek)
• LatinHypercubeSamplingisthemostpopularmethodofgeneratingquasi-random
• Matlab:lhsdesign()andlhsnorm()

Inverse transform sampling
The uniform distribution is the distribution directly generated by most PRNGs
• rand() function in MATLAB
The cumulative probability of sampling
a number between 0 and u1 is u1
Most simulations entail sampling random variables from something other than the uniform distribution
To get a sample from a different distribution, start with standard uniform and transform into the distribution you need.
This works because CDFs are themselves uniformly distributed

Inverse Transform sampling
Suppose we want to sample x from a distribution with the CDF given
The inverse transform method sets: byF𝑥𝑥=𝐹𝐹−1 𝑢𝑢, 𝑢𝑢~𝑈𝑈0,1
• where 𝐹𝐹−1 𝑢𝑢 is the inverse of F, and 𝑈𝑈(0,1) denotes the uniform distribution from [0,1]
• So, if we generate a u from U(0,1), then we can calculate the corresponding x that has cumulative probability u
•𝑥𝑥 =𝐹𝐹−1 𝑢𝑢 and𝑥𝑥 =𝐹𝐹−1 𝑢𝑢 1122
Example: distribution with [−∞. +∞] support 18

For continuous distributions
MATLAB often has functions for directly sampling from common continuous distributions
• E.g. exponential: exprnd() • E.g. normal: randn()
Let’s see how to do inverse transforms ourselves anyway
• Can use built-in Matlab functions to check
• Good practice with known and relatively easy distributions: prepares us for when we have to sample from a distribution that does not have a built-in Matlab random number generator

Inverse Transform Sampling in Practice
Given a uniformly distributed random sample u, solve F(x)=u to get a random sample distributed with CDF F(x)
Example: exponential with mean 1/𝜆𝜆 Let1−𝑒𝑒−𝜆𝜆𝑥𝑥 =𝑢𝑢 u0 Solve for x
𝒙𝒙=−𝟏𝟏 𝐥𝐥𝐥𝐥𝐥𝐥(𝟏𝟏−𝒖𝒖) 𝝀𝝀
We found the Inverse CDF of the Eexponential distribution This is what MATLAB does! Let’s look inside exprnd()

The Normal Distribution
Normal random variables are the building blocks of many financial simulation models
Basic properties • PDFforN(0,1):
𝜑𝜑 𝑧𝑧 = 1 𝑒𝑒−𝑧𝑧2 2𝜋𝜋 2
• PDFforN(𝜇𝜇,𝜎𝜎)
• CDFdoesnothaveananalyticalexpression.
𝑓𝑓 𝑥𝑥 =𝜑𝜑 𝑋𝑋−𝜇𝜇 i.e.𝑧𝑧=𝑧𝑧−𝜇𝜇or𝑥𝑥=𝜇𝜇+𝜎𝜎𝑧𝑧 𝜎𝜎𝜎𝜎
MATLAB functions for the normal PDF, CDF, and inverse CDF • PDF:normpdf(x) ornormpdf(x,mu,sigma)
• CDF:normcdf(x) ornormcdf(x,mu,sigma)
• Inverse CDF: norminv(x) or norminv(x,mu,sigma)

In-Class Exercise: E[g(X)]
Use simulation with 10,000 samples to find the expected value of a random variable 𝑌𝑌 = 1 , where random variable X is exponentially
distributed with rate parameter 𝜆𝜆 = 2 (in other words, 𝐸𝐸 𝑋𝑋 = 𝜆𝜆1 = 0.5)
• Part a: generate a random sample using MATLAB’s built-in exprnd()
• Part b: generate a random uniform sample using rand() and transform it to
exponential Display your results.
We will reconvene in 15 minutes
If you’re done early, think about how you would try to solve this problem without MATLAB. How would you set it up? What difficulties would you run into?

For discrete distributions
Example: the daily price fluctuation X of an asset
Prob𝑋𝑋=Δ1 =𝑝𝑝1
𝑋𝑋∈ −Δ2,0,Δ1 with�Prob𝑋𝑋=0 =1−𝑝𝑝1−𝑝𝑝2=𝑝𝑝3
So, we generate a random variable U(0,1) and X is given by
0 p1 p1+p2 1 Transform
Prob 𝑋𝑋=−Δ2 =𝑝𝑝2

For discrete distributions: CDF representation
Cumulative probability of X (HLS 2.1, p. 49) delta1 = 5,delta2 = 2,p1 = 0.3,p2 = 0.4
The “line segment” visualization makes it easy to figure out how to code this up, but we can also represent the problem the same way as we do for continuous variables 2 2
1 0.8 0.6 0.4 0.2 0
-4 -3 -2 -1 0 1 2 3 4 5 6 7 X
For𝑥𝑥<−Δ ,𝐹𝐹 𝑥𝑥 =Prob𝑋𝑋<−Δ =0: impossible For𝑥𝑥<0,𝐹𝐹 𝑥𝑥 =Prob𝑋𝑋<0 =Prob[𝑋𝑋= −Δ2] = 𝑝𝑝2 For𝑥𝑥<Δ ,𝐹𝐹 𝑥𝑥 =Prob𝑋𝑋<Δ =Prob[𝑋𝑋= 11 −Δ ] + Prob 𝑋𝑋 = 0 = 𝑝𝑝 + 𝑝𝑝 = 1 − 𝑝𝑝 2231 For𝑥𝑥≥Δ,𝐹𝐹𝑥𝑥 =1−Prob𝑋𝑋>Δ =1: 11
guaranteed
Cumluative Probability

In-Class Exercise: Benford’s Law
Problem: Benford’s Law says that naturally- occurring numbers tend to have a leading digit (i.e., first digit) that is not uniformly distributed since smaller leading digits tend to occur more often than large leading digits. Instead, the probability that the first digit of a number X is equal to d, where d is 1,𝑑𝑑2, … 9 is given as:
Benford’s Law can be used as a basic fraud detector
Say you’re analyzing a large dataset e.g.
• Sales records as part of a firm audit
• Election results: vote totals by region
If the first digits of the numbers deviate significantly from Benford’s Law, it may be that the numbers are fake
𝑃𝑃𝑋𝑋=𝑑𝑑 =log (1+1). 10
Write a script file that generates 1000 random leading digits from this distribution using the inverse transform method. Plot a histogram of the values you generated.
Use MATLAB functions log10() and histogram()

Generating correlated random numbers

Generating Multivariate Normals
Case of two jointly normally distributed random variables 𝑋𝑋~𝑁𝑁(𝜇𝜇,𝜎𝜎), 𝑋𝑋~𝑁𝑁(𝜇𝜇,𝜎𝜎),𝜌𝜌 =𝜌𝜌
1 1 12 2 2 2 𝑥𝑥1,𝑥𝑥2
𝜎𝜎2 𝜌𝜌𝜎𝜎 𝜎𝜎 𝑋𝑋𝜇𝜇112
• We can write down the joint distribution of these variables in vector form
𝑋𝑋1 𝜇𝜇 1 𝜌𝜌 𝜎𝜎1 𝜎𝜎2 𝜎𝜎2 2
Goal: generate a 1 x 2 vector
• Firstelementisasampleofrandom
variable 𝑋𝑋1
• Secondelementisacorrespondingsample
of random variable 𝑋𝑋2
Built-in MATLAB function:
mvnrnd(mean_vector,covariance_matr ix)
• See optional Appendix to learn how mvnrnd() works “under the hood”
mu = [mu1, mu2];
Sigma = [sig1^2, rho*sig1*sig2;
rho*sig1*sig2, sig2^2];
X = mvnrnd(mu,Sigma);

Generating Multivariate Normals: N x 2
But we rarely want to sample a distribution just once.
New goal: generate a N x 2 matrix
• First column: samples of X1
• Second column: corresponding samples of X2
Let’s generate N pairs:
mvnrnd(mean_vector,covariance_matrix,N)
Produces a N x 2 matrix
• First column: samples of X1
• Second column: corresponding samples of X2
• Again, transpose to get what we want
X = mvnrnd(mu,Sigma,N);

Generating Multivariate Normals: n x m
Specify 1 x m vector means (zeros if standard normal) Specify m x m covariance matrix
• correlation matrix if standard normal Specify length of sample n
Get n x m sample

Demo: Simulating Stock Returns

Simulating Stock Returns
Example: United Airlines
• Today’s stock price: 36.85 | expected daily return: 0.02% | daily vol: 1.4%
What is the expected stock price a year from now? What is the standard deviation? If we buy the stock today, how likely are we to lose money?
• 𝑆𝑆 = 36.85;𝜇𝜇 = 0.02;𝜎𝜎 = 0.014; 0
Define some variables
• A year is 365 calendar days, but only 252 trading days, so set T = 252
•𝑆𝑆 =𝑆𝑆𝑒𝑒𝑟𝑟 𝑇𝑇0
Stock price at time T can be represented as

Modeling Log Returns as Normally Distributed
Define a one-period log return 𝑟𝑟𝑖𝑖 = log 𝑅𝑅𝑖𝑖 ~𝑁𝑁 𝜇𝜇 − 𝜎𝜎2 , 𝜎𝜎2
• Identically normally distributed with mean 𝜇𝜇 − 𝜎𝜎2 and variance 𝜎𝜎2 2
• Why the − 𝜎𝜎2 term? So that 𝔼𝔼 𝑅𝑅𝑖𝑖 = 𝔼𝔼 exp(𝑟𝑟𝑖𝑖) can be written as 𝑒𝑒𝜇𝜇 Thenlog𝑆𝑆 =log𝑆𝑆 +𝑟𝑟 +𝑟𝑟 +⋯+𝑟𝑟 +⋯𝑟𝑟 =log𝑆𝑆 +𝑟𝑟
𝑇𝑇2 0 1 2 𝑡𝑡 𝑇𝑇 0 0→𝑇𝑇 Further assuming that they are independent (so iid),
• 𝔼𝔼 log𝑆𝑆 =log𝑆𝑆 +∑ 𝔼𝔼[𝑟𝑟]=log𝑆𝑆 +(μ− )𝑇𝑇
𝑇𝑇 0 𝑡𝑡=1 𝑡𝑡 0
•𝕍𝕍log𝑆𝑆 =∑𝑇𝑇 𝕍𝕍[𝑟𝑟]=𝜎𝜎2𝑇𝑇and𝜎𝜎log𝑆𝑆 =𝜎𝜎𝑇𝑇 𝑇𝑇𝑡𝑡=1𝑡𝑡 𝑇𝑇
Sum of iid normally distributed random variables is normal, so

Working with log returns, cont’d
The cumulative T-day log return 𝑟𝑟0→𝑇𝑇 is a normally distributed random variable with a mean of (μ − 𝜎𝜎2) 𝑇𝑇 and a standard deviation of 𝜎𝜎 𝑇𝑇
• Wecanwriteitas𝑟𝑟0→𝑇𝑇 =(μ−𝜎𝜎2)𝑇𝑇+𝜎𝜎 𝑇𝑇 ∗ε
• Where ε is a standard normal random variable (we know how to generate these in
Simulate annual log returns: 2
• generate N normally distributed random numbers with mean (μ − 𝜎𝜎 ) 𝑇𝑇 and standard deviation 𝜎𝜎 𝑇𝑇 and store them in a vector X 2
Compute year-ahead stock prices 𝑆𝑆𝑇𝑇 = 𝑆𝑆0𝑒𝑒𝑋𝑋
𝑆𝑆𝑇𝑇 now contains a distribution of end-of-year stock prices

What about a portfolio of stocks?
Profit or loss on a portfolio (e.g. portfolio of stocks) depends on the future prices of all assets (e.g. stocks) in the portfolio
Example: portfolio with 1 share each of Google (G), Apple (A), and
Facebook (F)
• Valuetoday:𝑆𝑆 +𝑆𝑆
• To get value at time T, we need 𝑆𝑆𝐺𝐺,𝑇𝑇, 𝑆𝑆𝐴𝐴,𝑇𝑇 , and 𝑆𝑆𝐹𝐹,𝑇𝑇
• 𝑆𝑆𝑖𝑖,𝑇𝑇 =𝑆𝑆𝑖𝑖,0exp 𝜇𝜇𝑖𝑖 −𝜎𝜎𝑚𝑚2 𝑇𝑇+𝜎𝜎𝑖𝑖 𝑇𝑇∗𝜀𝜀𝑖𝑖 foreachstock𝑖𝑖=𝐺𝐺,𝐴𝐴,𝐹𝐹
• What is the joint distribution of noise [𝜀𝜀𝐴𝐴,𝜀𝜀𝐴𝐴, 𝜀𝜀𝐹𝐹]′? In other words, how are the stocks
Once we simulate the correlated noise, we can a distribution of stock prices and hence portfolio value and profit/loss
correlated?
• Sample from the multivariate standard normal distribution

[𝜀𝜀 𝜀𝜀 ,𝜀𝜀 ]′ 𝐴𝐴,𝐴𝐴 𝐹𝐹
Sampling noise
1 𝜌𝜌𝐺𝐺,𝐴𝐴 𝜌𝜌𝐺𝐺,𝐹𝐹 𝜀𝜀𝐴𝐴~Normal0,𝜌𝜌𝐺𝐺,𝐴𝐴 1𝜌𝜌𝐴𝐴,𝐹𝐹
𝜀𝜀𝐹𝐹 0 𝜌𝜌𝐺𝐺,𝐹𝐹 𝜌𝜌𝐴𝐴,𝐹𝐹 1
• Need to sample a 3×1 vector of multivariate normally distributed random variables with mean 0 and covariance Σ
Sample correlated normally distributed random numbers
• noise = mvnrnd(zeros(1,3), correlation_matrix)

Capturing Rare Big Movements
Remember: normal distribution doesn’t do a great job capturing rare large (either positive or negative) returns
Our solution: add “jumps”
• On most days (probability 1 − 𝑝𝑝 − 𝑝𝑝 ), log returns are sampled from the
normal distribution 𝑟𝑟𝑡𝑡~𝑁𝑁(𝜇𝜇 − 𝜎𝜎2 , 𝜎𝜎2) 2
distribution (e.g. Pareto distribution)
o Good jumps: Positive jumps with probability 𝑝𝑝 : 𝑟𝑟 ~𝑃𝑃𝑚𝑚𝑟𝑟𝑒𝑒𝑃𝑃𝑃𝑃 𝑥𝑥 , 𝛼𝛼
𝐺𝐺𝑡𝑡 𝑃𝑃𝑃𝑃 o Bad jumps: Negative jumps with probability 𝑝𝑝𝐵𝐵 : −𝑟𝑟𝑡𝑡 ~𝑃𝑃𝑚𝑚𝑟𝑟𝑒𝑒𝑃𝑃𝑃𝑃 𝑥𝑥𝑁𝑁 , 𝛼𝛼𝑁𝑁
• On other days, log returns are sampled from an alternative fat-tailed
Let’s simulate the combined distribution (“mixture”)

How does mvnrnd() work?

Appendix: How Does mvnrnd() work?
Case of a two-dimensional random
1𝜌𝜌 𝑋𝑋20 𝜌𝜌1
, Σ where Σ =
1) Generate 2 x 1 vector Z of standard
normal random numbers (uncorrelated)
2) Compute the lower Cholesky decomposition L of Σ
3) Compute𝑋𝑋=𝐿𝐿 ∗𝑍𝑍
* indicates matrix multiplication
Sigma = [1, rho; rho, 1];
Z = randn(2,1);
L = chol(Sigma,’lower’);
X = L * Z;

Lower Cholesky decomposition
The lower Cholesky decomposition of a square positive-definite matrix Σ yields a matrix L of the same size such that
• Σ = 𝐿𝐿 ∗ 𝐿𝐿′ and L is lower triangular (all elements above the diagonal are zero)
• For2x2matrix,𝐿𝐿=
𝐿𝐿11 0 𝐿𝐿21 𝐿𝐿22
• Think of Cholesky decomposition as a matrix analog of the square root for scalars
o If I asked you to generate a sample of univariate random numbers and gave you the variance, you’d take the square root to get standard deviation.
o Square roots are unique: For any scalar, only one number exists such that its square is equal to that scalar
o Same isn’t true for matrices. Many matrices B exist such that Σ = 𝐵𝐵 ∗ 𝐵𝐵′ o But only one is lower triangular!

Constructing Cholesky decomposition
𝐿𝐿0𝐿𝐿0𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿 11 111121111121
For 2 x 2 matrix, 𝐿𝐿 = so Σ = ∗ = 2
22 𝐿𝐿21𝐿𝐿22 0𝐿𝐿22 𝐿𝐿𝐿𝐿𝐿𝐿+𝐿𝐿
o Σ is symmetric as any covariance matrix is, so off-diagonal elements the same so extra equation is redundant
• But we also know that Σ =
• System of 3 equations in 3 unknowns is easy to solve
o1=𝐿𝐿2 𝐿𝐿 =1 11 11
o 𝜌𝜌 = 𝐿𝐿11𝐿𝐿21 𝐿𝐿21 = 𝜌𝜌
o 1 = 𝐿𝐿2 + 𝐿𝐿2 𝐿𝐿 = 1 − 𝜌𝜌2
•𝑋𝑋=𝐿𝐿∗𝑍𝑍soX=𝐿𝐿𝑍𝑍=𝑍𝑍andX=𝐿𝐿𝑍𝑍+𝐿𝐿𝑍𝑍=𝑍𝑍𝜌𝜌+𝑍𝑍 1−𝜌𝜌 1 111 1 2 211 222 1 2
o Check means: E[X1] = E[𝑍𝑍1] = 0 and E[X2] = 𝐿𝐿21E[𝑍𝑍1] + 𝐿𝐿22E[𝑍𝑍2] = 0
o Checkvariances:Var[X ]=Var[𝑍𝑍 ]=1andVar[X ]=𝐿𝐿2 Var[𝑍𝑍 ] +𝐿𝐿2 Var[𝑍𝑍 ] =𝜌𝜌2 + 1−𝜌𝜌2 =1
11 2211222
o Check correlation: Cov[𝑋𝑋1,𝑋𝑋2] = Var[𝑍𝑍1]𝐿𝐿11𝐿𝐿21 = 𝜌𝜌 Var 𝑋𝑋1 Var[𝑋𝑋2]
• With more than 2 variables, this gets much harder to do by hand
o With lots of variables, this can get hard even for a computer e.g. 10,000 variables creates a 50M-dimensional system of
equations1 11 1 2
𝑋𝑋 = 𝐿𝐿 ∗ 𝑍𝑍 so X = 𝐿𝐿 𝑍𝑍 and X = 𝐿𝐿 𝑍𝑍 + 𝐿𝐿 𝑍𝑍
o Intuition: the correlation between X1 and X2 is created by the presence of Z1 in both

Generating Multivariate Normals
N random samples from the joint distribution of 𝑋𝑋 and 𝑋𝑋
12 • First row are samples of 𝑋𝑋1
• Second row are corresponding samples of 𝑋𝑋
• Each column vector is sample pair Still𝑋𝑋=𝐿𝐿 ∗ 𝑍𝑍
• 𝑋𝑋and𝑍𝑍arenow2xnmatrices • 2x2times2xnis2xn
Sigma = [1, rho; rho, 1];
Z = randn(2,n);
L = chol(Sigma,’lower’);
X = L * Z;

Generating m Multivariate Normals
Foronesample 𝑋𝑋1 ∗ 0 … 0 𝑍𝑍1 𝑋𝑋=𝐿𝐿∗𝑍𝑍→ 𝑋𝑋2 = ∗ ∗ … 0 𝑍𝑍2
𝑋𝑋⋮ ∗ ⋮ .⋮. ⋮ 𝑍𝑍⋮ 𝑚𝑚∗∗…∗𝑚𝑚
ForNsamples(XismxN),it’sstill𝑋𝑋=𝐿𝐿 ∗𝑍𝑍

Appendix: Generating Pseudorandom Numbers
𝑥𝑥𝑖𝑖+1 = 𝑚𝑚𝑥𝑥𝑖𝑖 mod 𝑚𝑚
A linear congruential generator (generate u)
• mod is modulus (remainder) after division o Ex: 7 mod 3 = 1, 11 mod 3 = 2, 13 mod 5 = 3
• 0<𝑥𝑥𝑖𝑖+1 <𝑚𝑚;0<𝑢𝑢𝑖𝑖+1 <1 • Example: a=6; m=11; 𝑥𝑥0 = 1 • Period length is (m-1) and the gap is (1/m) • Reproducibility and unpredictability based on 𝑢𝑢𝑖𝑖+1 = 𝑥𝑥𝑖𝑖+1/𝑚𝑚 • Usesaninitialvalue(seed)𝑥𝑥 ,0<𝑥𝑥 <𝑚𝑚 00 • RandStream.list lists all the generator algorithms in Matlab • See lecture code: Simulation_randomNumberGenerator.m • Note: The period of MATLAB’s default random number generator (the Mersenne Twister) is 2^19937-1! o See https://www.mathworks.com/help/matlab/ref/randstream.list.html 程序代写 CS代考加微信: powcoder QQ: 1823890830 Email: powcoder@163.com

Related Posts