Euler–Milstein and Monte Carlo Extensions

Apr 15, 2025 Sharif Haason 24 min read Study note

Computational Methods

Euler Milstein

We will expand on the variance reduction techniques, we will repeat a couple things from last week.

But before we do that, we should talk about Euler Milstein or Euler Maruyama.

𝑃 𝑒 𝑜 𝑝 𝑙 𝑒 𝑘 𝑒 𝑒 𝑝 𝑡 𝑎 𝑙 𝑘 𝑖 𝑛 𝑔 . 𝐼 𝑡 ℎ 𝑢 𝑟 𝑡 𝑠 𝐹 𝑙 𝑜 𝑟 𝑒 𝑠 𝑐 𝑢 𝑏 𝑒 𝑐 𝑎 𝑢 𝑠 𝑒 ℎ 𝑒 ℎ 𝑎 𝑠 𝑎 𝑑 𝑖 𝑠 𝑒 𝑎 𝑠 𝑒 𝑐 𝑎 𝑙 𝑙 𝑒 𝑑 𝐴 𝐷 𝐻 𝐷 .

This is a better approximation. I’m going to show you this because it involves applying Itô’s formula a bunch.

Here is the stochastic process we can approximate. Euler works for any stochastic process, but Milstein has to be homogeneous, like so:

𝑑 𝑋_{𝑡} = 𝛼 (𝑋_{𝑡}) 𝑑 𝑡 + 𝛽 (𝑋_{𝑡}) 𝑑 𝑊_{𝑡}

These coefficients are not functions of time, so they are homogeneous.

The idea is to use Itô for both $𝛼$ and $𝛽$ . They are functions of stochastic process $𝑋_{𝑡}$ . Therefore I can use Itô for both.

I get

𝑑 𝛼 (𝑋_{𝑡}) = 𝛼^{'} (𝑋_{𝑡}) 𝑑 𝑋_{𝑡} + \frac{1}{2} 𝛼^{″} (𝑋_{𝑡}) (𝑑 𝑋_{𝑡})^{2}

This is Itô, so we are lacking the dt, that term would screw up our calculations.

Now we substitute $𝑋_{𝑡}$ here, so we get

The dX^2 gets eliminated

= 𝛼^{'} (𝑋_{𝑡}) 𝛼 (𝑋_{𝑡}) 𝑑 𝑡 + 𝛼^{'} (𝑋_{𝑡}) 𝛽 (𝑋_{𝑡}) 𝑑 𝑊_{𝑡} + \frac{1}{2} 𝛼^{″} (𝑋_{𝑡}) 𝛽^{2} (𝑋_{𝑡}) 𝑑 𝑡

If I just group these terms and drop $𝑋_{𝑡}$ all over the place to make it easier to write:

= (𝛼^{'} 𝛼 + \frac{1}{2} 𝛼^{″} 𝛽^{2}) 𝑑 𝑡 + 𝛼^{'} 𝛽 𝑑 𝑊_{𝑡}

Now, if we do the same calculation/derivation for $𝛽$ , then we basically get something very similar. We’re not going to go through the entire calculation. You can derive this yourself,

𝑑 𝛽 (𝑋_{𝑡}) = (\dots) 𝑑 𝑡 + \dots + 𝑑 𝑊_{𝑡}

Writing in integral from,

Now we substitute $𝛼$ and $𝛽$ from the same formula. So there are two integrals:

X\_{t+\\Delta t} \- X\_t \= \\int\_t^{t+\\Delta t}

When we substitute, we will get four terms, two terms for $𝛼$ and $𝛽$ each.

We will get terms such as

𝑑 𝑠 𝑑 𝑢 \sim 𝑂 (Δ 𝑡^{2})

We will change the letters to make sure they’re right, based on what our dummy variables are

𝑑 𝑠 𝑑 𝑊_{𝑢} \approx 𝑑 𝑢 𝑑 𝑊_{𝑠} \sim 𝑂 (Δ 𝑡^{\frac{3}{2}})

Because $𝑑 𝑠$ is order $Δ 𝑡$ , and $𝑑 𝑊_{𝑡}$ is order $\sqrt{Δ 𝑡}$

Then there is

𝑑 𝑊_{𝑢} 𝑑 𝑊_{𝑠} \sim 𝑂 (Δ 𝑡)

For the reason just discussed.

Once I substitute everything, I will neglect the first two orders, and will be left with terms with just du and ds

𝑑 𝑢

𝑑 𝑠

The existing equation then simplifies a lot

𝑋_{𝑡 + Δ 𝑡} = 𝑋_{𝑡} + 𝛼 (𝑋_{𝑡}) Δ 𝑡 + 𝛽 (𝑋_{𝑡}) Δ 𝑊_{𝑡} + \int_{𝑡}^{𝑡 + Δ 𝑡} \int_{𝑡}^{𝑢} 𝛽_{𝑠}^{'} 𝛽_{𝑠} 𝑑 𝑊_{𝑠} 𝑑 𝑊_{𝑢}

And there’s another integral term you can see here.

You need to take the increments of Brownian motion and such, but then you can show that the integral term is equal to the

= 𝛽_{𝑡}^{'} 𝛽_{𝑡} \frac{1}{2} (Δ 𝑊_{𝑡} - Δ 𝑡)

This gives the Euler Milstein scheme.

The $𝑋_{𝑡 + Δ 𝑡}$ is what you’re approximating. The first part is the regular Euler, and then the integrals are the Milstein part.

If the model

𝑑 𝑋_{𝑡} = 𝛼 (𝑋_{𝑡}) 𝑑 𝑡 + 𝛽 (𝑋_{𝑡}) 𝑑 𝑊_{𝑡}

, then the Euler Milstein scheme is

𝑋_{𝑡 + Δ 𝑡} = 𝑋_{𝑡} + 𝛼 (𝑋_{𝑡}) Δ 𝑡 + 𝛽 (𝑋_{𝑡}) Δ 𝑊_{𝑡} + \frac{1}{2} 𝛽^{'} (𝑋_{𝑡}) 𝛽 (𝑋_{𝑡}) (Δ 𝑊_{𝑡}^{2} - Δ 𝑡)

And by $𝑊$ you introduce a normal variable multiplied by

𝑍 \sim 𝑁 (0, 1)

. The $Δ 𝑊$ is created by

𝑍 \sqrt{Δ 𝑡}

. The other one is

(𝑍^{2} - 1) Δ 𝑡

when the thing factors. For the same increment, you put it in two places, not just one places.

However, there is a possible issue, which might be the derivative. If the $𝛽$ function, the volatility part, if it’s complicated, how do you calculate the derivative? That might be hard. There’s another way to deal with this, a scheme called Runge-Kutta, a generalization where you calculate the Euler part, and plug it into the $𝛽$ , then you calculate the finite difference as an approximation of the derivative.

One more thing to mention. Which I shouldn’t, because it’s from the homework. Let’s have an example.

We have a process

𝑑 𝑌_{𝑡} = 𝜅 (\bar{𝑦} - 𝑌_{𝑡}) 𝑑 𝑡 + 𝜎 \sqrt{𝑦_{𝑡}} 𝑑 𝑊_{𝑡}

This is the CIR process.

In this process, I have my

𝛼 (𝑥) = 𝜅 (\bar{𝑌} - 𝑥)

Then you have $𝛽$

𝛽 (𝑥) = 𝜎 \sqrt{𝑥}

If we now substitute in this formula, we have to calculate the derivative.

𝛽^{'} (𝑥) = \frac{𝜎}{2 \sqrt{𝑥}}

Now if we do Euler Milstein:

𝑌_{𝑡 + Δ 𝑡} = 𝑌_{𝑡} + 𝜅 (\bar{𝑌} - 𝑌_{𝑡}) Δ 𝑡 + 𝜎 \sqrt{𝑌_{𝑡}} Δ 𝑊_{𝑡} + milstein correction: 𝜎 \sqrt{𝑌_{𝑡}} \frac{𝜎}{2 \sqrt{𝑌_{𝑡}}} (Δ 𝑊_{𝑡}^{2} - Δ 𝑡)

Here, the two $\sqrt{𝑌_{𝑡}}$ terms cancel, so it becomes $𝜎^{2} / 2$ .

The other thing to show is that if you look for example, this process.

𝑑 𝑋_{=} 𝑒^{𝑋_{𝑡}} cos 𝑋_{𝑡} 𝑑 𝑡 + 0.7 𝑑 𝑊_{𝑡}

Here, you have this constant

𝛽 (𝑋_{𝑡}) = 0.7

𝛽^{'} (𝑋_{𝑡}) = 0

Therefore there is no Euler-Milstein correction, because it relies on multiplication by the derivative. Like for example in GBM.

Variance Reduction Redux

More of the idea to use on HW 4.

One additional note will be given on antithetic variate.

We discussed the CLT and the basis of the whole thing, and the variance being smaller gives you better estimates.

However, this is more complicated than that.

The better method is not always the one with smaller variability of the sample paths.

We also need to consider the time to generate paths. If you have a Monte Carlo technique with less variability, but it takes one minute to generate a path, then it could be that another method with much higher variability that can generate 1 per second, is better, just by raw brute force. The paper mentioned is really good, the basis of the Monte Carlo book, which expands on the paper. Boyle is a Canadian professor from University of Waterloo, in 1997 they met, he did a summer school. He was drunk all the time in the morning lectures. Glasserman is a friend of the show as well. Third guy Brodie sucks.

If method 1 has variance $𝜎_{1}^{2}$ , and $𝑏_{1}$ , they have a term called “work”, which could be time, but could be other stuff, $𝑏_{1}$ is the work to generate one replication of the final parameter, the one you’re trying to estimate, then we need to look at

𝜎_{1}^{2} 𝑏_{1} < 𝜎_{2}^{2} 𝑏_{2}

He has a better expression, but I’m showing the simple thing. Why multiplication? There’s a reason

If you look at this perspective, and if you rewrite this as

\frac{𝜎_{1}^{2}}{𝜎_{2}^{2}} < \frac{𝑏_{2}}{𝑏_{1}}

That’s when you would prefer one method over another.

Antithetic note:

Reminder: The way it works is you create one normal for one path, and use the negative of that normal for another path.

For antithetic variates, one replicate is

\frac{𝑐_{𝑖} + 𝑐_{𝑖}^{𝑎}}{2} = \bar{𝐶_{𝑖}}

The actual estimate is the average of the 2. This is not necessarily important for the final pricing, because if you take the general average, it doesn’t matter, but it’s important for the variance estimate, for confidence intervals.

This is because these paths are not independent, they’re very related. So you would use something like

Variance = \frac{1}{𝑛 - 1} \sum_{𝑖 - 1}^{𝑛} (\bar{𝐶_{𝑖}} - \bar{𝐶})^{2}

This is important because each sequence of n is independent.

You might think I should iterate over 2n elements instead of n. and do both the regular and antithetic.

But this one will show less than it really is.

Control Variates for the Asian Option

Glasserman’s explanation is better than mine, I will follow him.

The one with $Δ$ -hedging and Asian option. There’s nothing wrong with $Δ$ -hedging or Asian option, it’s just that it needs to be explained where it’s coming from.

The control variate for Monte Carlo, what it does, it should be called “use what you know”. You know, for example, that the option that $Δ$ -hedges, you’re using the concept that if you do this very fast, the two values should be the same.

Similarly, for the Asian option, we’re using what we know.

𝑃_{𝐴}

is the price of an Asian option based on Arithmetic average, which is what is encountered in practice.

𝑃_{𝐺}

is geometric.

This $𝑃_{𝐺}$ has a formula, and for any variation on Asian options it’s very easy to get a geometric formula.

We’ll say we have a formula, if you give me characteristics and $𝜇$ and $𝜎$ , you get an exact number.

With this, let

\hat{𝑃_{𝐴}}

be the value calculated using a single path.

It doesn’t matter what you use, just one single path, where you calculate the value of the Asian option based on this average of all these paths.

And we’ll do the same for geometric

\hat{𝑃_{𝐺}}

We don’t need this for geometric because we have the formula.

We are using the same path for both of these.

I do know that

𝔼 [\hat{𝑃_{𝐴}}] = 𝑃_{𝐴}

On expectation, I get the true vlau eof my option.

I also know

𝔼 [\hat{𝑃_{𝐺}}] = 𝑃_{𝐺}

If we subtract, we get

𝑃_{𝐴} - 𝑃_{𝐺} = 𝔼 [\hat{𝑃_{𝐴}} - \hat{𝑃_{𝐺}}]

This gives you a very natural kind of estimate.

𝑃_{𝐴} = 𝑃_{𝐺} + 𝔼 [\hat{𝑃_{𝐴}} - \hat{𝑃_{𝐺}}]

We can create a Monte Carlo path using this control variate:

{\hat{𝑃_{𝐴}}}^{𝑐 𝑣} = \hat{𝑃_{𝐴}} - \hat{𝑃_{𝐺}} + 𝑃_{𝐺}

If I create a new path and average this, then it should give me on average this difference

You can write it like this

= \hat{𝑃_{𝐴}} + (𝑃_{𝐺} - \hat{𝑃_{𝐺}})

That parentheses statement is the control variate. I have my original path with Euler Milstein, then I control it.

For every path, I obtain the difference between the true value and the value of geometric from that particular path. But is this better than just using

\hat{𝑃_{𝐴}}

If we want, we can calculate the variance. So the variance should be less. And remember that

𝑃_{𝐺}

is just a number, a constant with no variance.

𝕍 [{\hat{𝑃_{𝐴}}}^{𝑐 𝑣}] = 𝕍 [\hat{𝑃_{𝐴}}] + 𝕍 [\hat{𝑃_{𝐺}}] - 2 Cov (\hat{𝑃_{𝐴}}, \hat{𝑃_{𝐺}})

Everything after the first term should be negative to give me a better variance. The covariance should be greater than the variance. It’s only worth it if the covariance is large.

This brings the next idea. This is the original term plus this term. I can control the size of the difference, which puts a $𝛽$ on the coefficient. That $𝛽$ allows me to make the thing smaller. This is all specific to the Asian option, where there is this arithmetic and geometric thing. But there is no assumption about the stochastic model.

$𝛽$ means I’m going to parameterize this:

{\hat{𝑃_{𝐴}}}^{𝛽} = \hat{𝑃_{𝐴}} + 𝛽 (𝑃_{𝐺} - \hat{𝑃_{𝐺}})

I can play around with $𝛽$ in such a way that the resulting variance is the smallest. What is that? We can do the same exact calculation, and minimize the result with respect to $𝛽$ .

𝕍 [{\hat{𝑃_{𝐴}}}^{𝛽}] = 𝕍 [\hat{𝑃_{𝐴}}] + 𝛽^{2} 𝕍 [\hat{𝑃_{𝐺}}] - 2 𝛽 Cov (\hat{𝑃_{𝐴}}, \hat{𝑃_{𝐺}})

This is a quadratic expression. It’s a parabola, so the smallest value is in the vertex. The position of the vertex is obtain for

𝛽_{min} = - \frac{𝑏}{2 𝑎}

Where the varG is a, and the cov stuff is b.

= - \frac{- 2 Cov (\hat{𝑃_{𝐴}}, \hat{𝑃_{𝐺}})}{2 𝕍 [\hat{𝑃_{𝐺}}]}

This cancels to:

= \frac{Cov (\hat{𝑃_{𝐴}}, \hat{𝑃_{𝐺}})}{𝕍 [\hat{𝑃_{𝐺}}]}

If you have two variables, the regression is the covariance divided by the variance, so this is the formula for market $𝛽$ /regression.

So if we regress

\hat{𝑃_{𝐴}} = 𝛼 + 𝛽 \hat{𝑃_{𝐺}} + 𝜖

, the $𝛽$ is the $𝛽$ .

Last week, we learned this was true, but now we know how to get it.

This control variate thing is basically

{\hat{𝑃_{𝐴}}}^{𝛽} = \hat{𝑃_{𝐴}} + \hat{𝛽} (𝑃_{𝐺} - \hat{𝑃_{𝐺}})

Now we have another problem.

This is kinda screwed up. Because you’re using the same paths to estimate $𝛽$ , and the same path to estimate the value of the option. That introduces a bias, and this is calculated in the Glasserman paper.

Typically you have $𝑛$ paths, and you set $𝑛_{1}$ paths out to do regression. Then you use $𝑛 - 𝑛_{1}$ paths for calculation.

The advantage here is to do a proper regression, you don’t need a lot of observations, 100 would be plenty. But for Monte Carlo, you need hundreds of thousands.

Hopefully this is more clear, and what I hope you get is that he used the particular relationship that exists in the Asian option, to reason through the whole thing.

I’m going to make another expansion to this. We can introduce more control!!

It’s not really necessary because this existing control variate already gives a good estimate, but this shows you can introduce as many control as you like!

Under risk-neutral equivalent martingale measure, we have

𝑆_{0} = 𝔼^{𝑄} [𝑆_{𝑇} 𝑒^{- 𝑟 𝑇}]

If you take the stock price as a martingale, and discount it back, you should get $𝑆_{0}$ . Nothing new.

This brings up another way to control.

We use the following, with $𝛽_{1}$ being our original control variate.

{\hat{𝑃_{𝐴}}}^{𝑐 𝑣} = \hat{𝑃_{𝐴}} + 𝛽_{1} (𝑃_{𝐺} - \hat{𝑃_{𝐺}}) + 𝛽_{2} (𝑆_{0} - \hat{𝑆_{𝑇}} 𝑒^{- 𝑟 𝑇})

Now you have the path, you know what $𝑆_{𝑇}$ is, and $𝑆_{0}$ is a constant.

You have the same kind of regression of

\hat{𝑃_{𝐴}} = 𝛼 + 𝛽_{1} \hat{𝑃_{𝐺}} + 𝛽_{2} \hat{𝑆_{𝑇}} + 𝜖

The constants don’t matter because you’re doing a regression, it just changes the y-intercept, doesn’t impact the $𝛽$ .

Technically, it’s $𝜎 \sqrt{𝑡}$ , which is the confidence interval size, the diffusion size. It’s an estimate, work could refer to other things.

Moment Matching Method

This is simple to use, so I’ll mention it, even though it’s useless.

If we have

𝑍_{1}, \dots 𝑍_{𝑛} \sim 𝑁 (0, 1)

They should have theoretical mean 0, but the sample mean is not 0.

The idea is to modify the sample to match the theoretical moments.

The reason I actually have never taught this method is because it’s kinda stupid, as a statistician.

Also, from the raw power of this method, it doesn’t do better than straight Monte Carlo. It does better when you pair it with control variate, where the power comes from the control variate.

In the example here, the random variables don’t have mean 0, so instead use

𝑍_{1} - \bar{𝑍}, 𝑍_{2} - \bar{𝑍}, \dots

These now all have mean 0, but now they’re correlated. So it creates problems when estimating stdev, it becomes bad, very tricky for estimating confidence intervals.

In the example, let’s say I’m going to price a European option based on GBM. When you do GBM, you don’t have to do all the intermediate steps, you can do all in one step. Because the terminal value

\tilde{𝑆_{𝑇}} (𝑖) = 𝑆_{0} 𝑒^{(𝑟 - \frac{𝜎^{2}}{2}) 𝑇 + 𝜎 \sqrt{𝑇} {\tilde{𝑍}}_{𝑖}}

After you modify the normal variable with the thing I said.

This is fine for European options, but it won’t work for path-dependent options.

Confidence intervals are hard to obtain.

This is the first order moment matching. You can also do second order moment matching. Say we want to create

𝑁 (𝜇_{𝑍}, 𝜎_{𝑍}^{2})

I’m trying to create numbers that are normal with this particular target. The usual thing to do here is

𝑍_{𝑖} \sim 𝑁 (0, 1) \to 𝜎_{𝑍} 𝑍_{𝑖} + 𝜇_{𝑍}

Multiply to create the desired distribution.

But the numbers in the sample will have their own sample mean and stdev. So you have to modify this as

\tilde{𝑍_{𝑖}} = \frac{𝜎_{𝑍}}{𝑆_{𝑍}} (𝑍_{𝑖} - \bar{𝑍}) + 𝜇_{𝑍}

Each number is modified by these two numbers $𝑆_{𝑍}$ and $\bar{𝑍}$ , where $𝑆_{𝑍}$ is the sample stdev.

𝑆_{𝑍} = \sqrt{\frac{1}{𝑛 - 1} \sum (𝑍_{𝑖} - \bar{𝑍})^{2}}

For the random variables in the sample, they will have the desired distribution.

Like I said, this is a method from the 90s. I never liked it because it’s slower. All these modifications…

And in order to estimate the sample $𝜇$ and stdev, I have to do all of the simulations first, then calculate the samples, then plug them back, so it’s SLOWER than doing it all at once.

And generally, from my experience, improvement is marginal, it’s not particularly useful.

Is the sample calculated in one path or all paths? It depends. In this example, you need all the paths. But you can do it for one path if you like.

Computationally, when you have numbers that you’re observing, if I have to estimate an average or variance of a certain number of them, theoretically, I can do this.

The idea of the GPU is that you have a lot of matrices and you do these calculations really fast. There is a memory of GPU and a memory of CPU. You want to keep everything in GPU memory as much as possible. The problem is that this method aggregates the paths, and does it again, which is impossible with a GPU.

Additional Stuff

Stratified sampling is not the best explained here. It is a statistics technique, check any statistics book they will explain better there.

Importance sampling, same deal. I have a chapter in my book about importance sampling which is way better.

Conditional MC: very well explained in the paper

Low discrepancy sequences. (Quasi-Monte Carlo techniques). It’s not explained well in the paper, but you can look it up and read better papers.

Briefly:

when you start generating one dimensional random variables, we’ll say for a uniform distribution, we can use testers to see if the numbers are actually uniform.

Then you can generate pairs, two at a time $(𝑋_{1}, 𝑋_{2})$ and plot them. You should definitely do this experiment with a random number generator. It would not look uniform at all. Human mind when you say uniform, thinks that it’s perfectly spread out. The point of the quasi generator is to spread out on purpose, to get a perfect uniform distribution. It’s okay as an exhaustive search method. But if you need more points, you have to quadruple them to have the same spread everywhere. The more detail you want, the finer quasi becomes, and it gets a lot slower.

But there are circumstances in which this is useful, and people in engineering that don’t understand randomness like this thing.

Chapter 5 is about estimating American options using Monte Carlo simulations.

subjective probability

George Calhoun sent an article to me in Nature

Trump’s election is a subjectiv eprobabiolity, it matters what people’s perceptions are.

A stock, TSLA. It’s been going down, so you’r wondering if it should keep going down or should it stabilize? And we don’t know because the stock is not the value of the company, it reflects the perception of people about the company.

It’s the same as poker. Game theory is so close to probability.

When people lose concentration, they react poorly.

Generating Correlated Processes

If you have multiple stock processes that you want to generate that are correlated, how do you do that?

Each asset looks like

𝑑 𝑋_{𝑡}^{𝑖} = 𝑓 (𝑋_{𝑡}^{𝑖}) 𝑑 𝑡 + 𝑔 (𝑋_{𝑡}^{𝑖}) 𝑑 𝑊_{𝑡}^{𝑖}

Each is a separate stochastic process, but they don’t move their own way at random.

This Brownian motion is correlated.

If we take the notation

𝑋_{𝑡}

the vector of components, and same for f(x). This is the simpler case. You can have a more complicated thing. Each function could depend on all the other ones, it doesn’t have to be driven by just one variable.

g(x) is a matrix dxd size, where dW is dx1, the vector of Brownian motion.

When you do the Monte Carlo simulation, the Brownian motion is what you’re interested in simulating.

We will assume that these are correlated. With this notation, we have

𝑑 𝑋_{𝑡} = 𝑓 (𝑋_{𝑡}) 𝑑 𝑡 + 𝑔 (𝑋_{𝑡}) 𝑑 𝑊_{𝑡}

It’s more complicated because it’s a collection of multiple integrals on dW, but we will write like this for conciseness.

We have the covariance matrix

Cov (Δ 𝑊_{𝑡}) = Σ Δ 𝑡

This is Brownian motion. So what exactly is the covariance matrix? The components are the covariances of each of the components.

[\begin{matrix} 1 & 𝜌_{12} & 𝜌_{13} & 𝜌_{1 𝑑} \end{matrix}]

Symmetric and positive definite matrix!

Positive is vector times matrix times vector transpose, which is a 1x1 number. Positive definite means that for a vector u, this is always greater than 0. This is very important because you can take a vector multiply with vector which will give us the variance, which will be always positive.

For example, take Heston

𝑑 𝑆_{𝑡} = 𝑟 𝑆_{𝑡} 𝑑 𝑡 + \sqrt{𝑌_{𝑡}} 𝑆_{𝑡} 𝑑 𝑊_{𝑡}^{1}

𝑑 𝑌_{𝑡} = 𝛼 (\bar{𝑌} - 𝑌_{𝑡}) 𝑑 𝑡 + 𝜎 \sqrt{𝑌_{𝑡}} 𝑑 𝑊_{𝑡}^{2}

If I want to simulate this, which is part of homework, how would you do this?

Generate two random numbers which are correlated with $𝜌$ .

For two, it’s very simple.

We need

𝑋_{1}, 𝑋_{2} \sim 𝑁 (0, 1)

, such that

Corr (𝑋_{1}, 𝑋_{2}) = 𝜌

Then we can multiply by $\sqrt{Δ 𝑡}$ and everything will work.

We will start with uncorrelated $𝑍_{1}, 𝑍_{2} \sim 𝑁 (0, 1)$ . I’ll take $𝑋_{1} = 𝑍_{1}$ .

Then, I’ll take $𝑋_{2} = 𝜌 𝑍_{1}$ . We can do this because the variance of $𝑍_{1}$ is 1, and covariance of $𝑍_{1}$ and $𝑍_{2} = 0$ ;

The combination has to have variance of 1, and if we’re combining linear combinations, then it’s normal.

So total value is

𝑋_{2} = 𝜌 𝑍_{1} + \sqrt{1 - 𝜌^{2}} 𝑍_{2}

If you understand the principle, what do I do when I want to do three? $𝑋_{1}, 𝑋_{2}, 𝑋_{3}$

Take the same idea

𝑋_{1} = 𝑍_{1}

𝑋_{2} = 𝜌_{12} 𝑍_{1} + \sqrt{1 - 𝜌_{12}^{2}} 𝑍_{2} ? ? ? ?

𝑋_{3} = 𝜌_{13} 𝑍_{1} ? ? ? ?

It becomes tricky. So is there a method to do this?

There is! You should know where this is coming from

We have a vector $𝑋$ which is a lot of $𝑋$ s. We are interested in covariance, which is distinct from correlation. Our covariance matrix will be one diagonal, where we multiply by $\sqrt{Δ 𝑡}$ . The values on the diagonal of the matrix are the variance, and the rest are covariance.

So how do I generate a vector with this covariance structure? Here is the idea.

If $𝑋$ is a random vector, (more details you could talk about in 540) with mean $𝜇$ , componentwise for each element in vector, then

Cov (𝑋) = 𝔼 [(𝑋 - 𝜇) (𝑋 - 𝜇)^{𝑇}]

Take $𝑌 = 𝐴 𝑋$ . $𝐴$ is a matrix, but it can be ANY DIMENSION $𝑛 \times 𝑑$ . We can transform 4 components into 15 components with different linear combinations.

The basis of the whole method is this:

Cov (𝑌) = 𝔼 [(𝑌 - 𝜇_{𝑌}) 𝑌 - 𝜇_{𝑦})^{𝑇}]

And we know that the expectation of Y from linear combination si

𝔼 [𝑌] = 𝐴 𝔼 [𝑋]

Because of this, we get

= 𝔼 [(𝐴 𝑋 - 𝐴 𝜇_{𝑋}) (𝐴 𝑋 - 𝐴 𝜇_{𝑋})^{𝑇}]

And you can factor this. Remember that the order is very important.

= 𝔼 [𝐴 (𝑋 - 𝜇_{𝑋}) (𝐴 𝑋 - 𝐴 𝜇_{𝑋})^{𝑇}]

And we can also do this in the transpose

= 𝔼 [𝐴 (𝑋 - 𝜇_{𝑋}) {(𝑋 - 𝜇_{𝑋})}^{𝑇} 𝐴^{𝑇}]

And the inner product is the covariance matrix! So it ends up being this

= 𝐴 Σ 𝐴^{𝑇}

CHOLESKY DECOMPOSITION

Take the $𝑍$ vector of iid normals, which are $𝑁 (0, 𝐼_{𝑑})$ . On the diagonal, you have 1 correlation.

Find a matrix A such that

𝐴 𝐴^{𝑇} = Σ

, the desired covariance.

The reason why you do this is because you have A, and you multiply it, and you get the thing you need.

Cholesky uses eigenvalues and eigenvectors, but it does exactly this.

Covariance vs Correlation

How do you get from the covariance matrix to correlation?

𝜌_{𝑖 𝑗} = \frac{𝜎_{𝑖 𝑗}}{𝜎_{𝑖} 𝜎_{𝑗}}

How do you go from this method to other methods and vice versa?

In terms of matrix operations, it’s not complicated.

If D is the diagonal matrix of the variances (just the diagonal), then the correlation matrix is this. Inverse is regular inverse because it’s just a diagonal.

= {\sqrt{𝐷}}^{- 1} \times Σ ({\sqrt{𝐷}}^{- 1})^{𝑇}

And the transpose is actually the same thing.

And you do non-inverse to get from correlation to covariance.

Not commutative, so be careful.

Markov Chain Monte Carlo

Left unsaid

Back to all writing