Applied Statistics for finance - Appunti completi (ENG)

Applied statistics for finance – lecture notes

Random numbers are not really random because computers execute instructions in a deterministic way,

hence computers cannot produce any randomness. We can simulate numbers from only one distribution

which is Uniform [0, 1], a continuous distribution with minimum 0 and maximum 1 and with equal

probabilities. Each software allows us to reshape a pre-specified sequence in a certain distribution (like the

normal one) and we get different values every time because the starting point/seed of the sequence is

!

different case by case (an algorithm tells R to do it). Pseudo-random numbers are such that, given the

generated value it is possible to predict the next pseudo-random number via the iterative formula:

" "#$

= ( )

"#$ "

Provided that we know the algorithm or the Random Number Generator (RNG). Remark: if we start

(. ),

from the same initial value , called the seed of the RNG, we get the same deterministic sequence.

!

Knowing the starting value and the algorithm which we are using we can obtain the entire sequence

(forward looking algorithm) but this is not backward looking (given we can’t obtain ).

"#$ "

When is enough? The inverse-transform method

(, ) %&$

This is the first method. I generate random numbers from a uniform and I use the inverse to obtain

()

random numbers from a r. v. with cumulative distribution function For discrete random variables,

().

&$

the inverse is always easy to obtain. For continuous r.v., the inverse function is not always available.

The only simple and very special case is the exponential random variable Indeed:

~().

&'%

() = 1 −

Set () = : 1

&'%

(1 − ) = → − ln(1 − ) =

If also Therefore:

~(0,1) 1 − ~(0,1). 1

= − ln()

With is such that Basically we start from the desired distribution, we put it equal to

~(0,1) ~().

the uniform and then we invert it to find the .

When is not enough? The acceptance/rejection method

(, )

If we can’t isolate the because the (CD) function is not invertible (think at the normal) and we can’t use

numerical methods (trying different value for the till we make the equality with true) we use the so-

called acceptance/rejection method. Let be a random variable with density for which no RNG exists.

Assume you know how to simulate from with density and there exists a constant such that:

> 0

()

max ≤ 1, ∀

()

We want it to be close as much as possible to 1. Simulate and If:

~ ~(0,1).

()

≤ ()

Then is taken as simulated value of Otherwise reject and iterate. Basically we create an envelope

.

(container) called and we multiply it for a constant (to make it bigger) and we extract numbers from

()

this. If these numbers fall inside the area of we accept them, otherwise we reject.

()

Finding the right size of the envelope

We can’t take an envelope too big otherwise we have a computational loss, it must be the right size that

contains my function. It’s therefore a maximization problem: I am looking for a and for a that make

()

my ratio as close as possible to 1:

• I am therefore looking for a constant which is big enough so that can contain but I

() (),

want the bare minimum for that to happen: the bigger the higher the risk of falling into the

(),

rejection area and therefore of simulating values that I cannot use;

• I want a which is similar as possible in shape to my so that I do not lose simulations

() ():

due to big areas in which the two distributions and do not overlap.

() ()

Remember that the envelope is not forcedly a uniform distribution (it can also be an exponential one).

The function set.seed()

The function set.seed() allows you the reproducibility of your results, but only if the seed is always the

same (the same number in the brackets). If you change number inside the brackets (e.g. set.seed(111) vs

set.seed(123)) the starting value from the pseudo-random sequence from the U(0, 1) random variable will

be different, and therefore the generated pseudo-random values will be different.

Monte Carlo method

The Monte Carlo method is a numerical method based on statistical arguments which can be used to

generate draws from a probability distribution and to evaluate integrals. In particular, the expected value of

a random variable with density is an integral of this form:

(. ) () = I ∙ ()

In finance, we use the Monte Carlo Method to evaluate the expected payoff of some derivative in order to

price it. Suppose we are interested in calculating the expected value with is any function and

[()],

is a given random variable. Assume we know how to simulate pseudo-random numbers

, … ,

$ "

according to the distribution of Then, we can think to approximate with the arithmetic mean of

. [()]

the numbers ( ):

( "

1 )

[()] ≅ Q ( =

( "

()$

Is this a good approximation? Yes! This is guaranteed by the Law of Large Numbers. In summary, the

number we estimate by simulations, will have a precision, around the true value of order

[()]

"

equal to 1/√. Accuracy and speed of convergence

In some cases, Monte Carlo intervals are not very informative if the variance of is too large, as

= ()

variability affects stability of the Monte Carlo method. There are many techniques to constrain instability:

such techniques are called variance reduction techniques. Examples of existing techniques are preferential

sampling, control variables, antithetic sampling, etc.

Antithetic sampling

The idea of antithetic sampling can be applied when it is possible to find transformations of that leave its

distribution unchanged (for example, if is Gaussian, then is Gaussian as well). Basically in option

−

pricing the idea is to extract half of the values out of our and the others out of

() (−).

Relationship with option pricing

In the Black & Scholes theory of option pricing, we model an underlying asset using the gBm with this SDE:

= +

* * * *

The price of a derivative is given by the general formula:

&+(-&*) )}

= {(

* -

Where is a payoff function. We need to estimate via Monte Carlo the expected value and, to this end,

(. )

we need to be able to simulate the values of . So we need to know how to simulate stochastic processes.

-

Simulation of stochastic processes

A stochastic process is the rule that tells us how a random variable develops through time. The difference

with a distribution (density function) is that this one doesn’t take into account time, it can only describe

what happens to a random variable at a certain point in time. The stochastic process has inside whatever

distribution but it also includes the factor of time. A first classification of stochastic processes is:

We will focus on continuous stochastic

processes. The fact that the process is

continuous through time means that for every

point in time we should be able to get a value

for that stochastic process. Unfortunately, we

need to simulate the process so technically it

can’t be continuous because we will discretize

time (setting different discretization points).

The distance between two discretization points doesn’t necessarily need to be always the same (if we

choose to be equal, given for example 100 discretization points the length of each one is 1/100).

Brownian motion

Brownian motion (or Wiener process but different from the geometric one) is a stochastic process (with

starting from zero at time 0, i.e. with the following properties:

, ≥ 0) = 0,

* !

1. It has independent increments: i.e. and are independent for

− − (, ) ∩ (, ) = ∅.

* / 0 1

What happens between two grid points is independent of what happens between other two points;

2. It has stationary increments: i.e. the distribution of for depends only on

− > ≥ 0 −

* /

and not by and/or separately. The distributional assumption that governs two grid points is the

same one that governs other two grid points;

3. It has Gaussian increments: i.e. − ∼ (0, − ).

* /

In particular We also use the notation . The grid point at is obtained as:

∼ (0, ). =

* * * $

1

f

)

( = 0 + ∙ (0,1)

$ 100

Where 0 is the starting point. At we simply substitute 0 with the value obtained. Problems: starting value

2

in zero and it can also go negative. Given a fixed time increment one can easily simulate a

∆ > 0,

trajectory of the Wiener process in the time interval [0,]. Indeed, for it holds true that:

∆*

(∆) = (∆) − (0)~(0, ∆)~√∆ ∙ (0,1)

And the same is also true for any other increment ( + ∆) − ()~(0, ∆)~√∆ ∙ (0,1).

How to simulate the path for a Brownian motion

The idea is to simulate a lot of paths in order to apply the LLN (because Monte Carlo method relies on it).

Divide the interval [0,] into a grid such as with Then

0 = < < ⋯ < < = − = ∆.

$ 2 4&$ 4 (#$ (

)

we set and and iterate as:

= 1 (0) = ( = 0

$

1. Generate a (new) random number from the standard Gaussian distribution;

2. = + ;

) )

3. Set ( = ( + ∙ √∆;

&

4. If iterate from step 1: in R the output will be a numerical sequence of length (1 path).

≤ , + 1

Be very clear at the exam about the output of a code. We can’t use the Wiener process for simulating a

path for an underlying (it’s a component of a more complicated process we will use) because its problem is

that it considers extractions from a standard normal (which are positive and negative) so it behaves around

the mean of zero (the average simulated path because some paths will be positive but others will be

negative). We use the Wiener process to represent the volatility parts of the underlying.

Some of the codes shown are inefficient as they implement

an iteration using a for statement when this is not needed

at all. But let’s see this one:

W < −c(0, cumsum(sqrt(Delta) × rnorm(N))

In this simple case, the whole trajectory can be simulated

in one line of R code as follows due to the fact that the

algorithm is just simulating a random walk with Gaussian

increments of size We do all together.

∙ (0,1).

√∆

Brownian motion paths are nowhere differentiable

If we compute the derivative of a Brownian motion, the limit explodes to infinite. This happens because of

the independence of the increments of the Brownian motion and also, more importantly, because the

increments behave like instead of Thus, in the limit one can expect:

( + ∆) − () ∆.

√∆

|( + ∆) − ()| •√∆•

lim = lim = +∞

∆ ∆

∆*→! ∆*→!

Stochastic differential equations

These are the mathematical formulas that describe a stochastic process. They are called like this because

they have derivatives inside. Every stochastic process will have its behaviour explained through an SDE.

Given a smooth, non-stochastic dynamical system its evolution with respect to time can be

= (),

*

represented as: −

*#8* * * ) )

= = ( = (

* * *

A stochastic differential equation models the noise (or the stochastic part) of this system by adding the

variation of some stochastic process to the above dynamics, e.g. the Wiener process:

) )

= ( + ( = + ℎ

* * * *

Every SDE/process we will is see is the sum of the deterministic trend plus some noise around that trend.

So it’s the sensitivity to time (drift, derivative with respect to time plus the sensitivity to the variation

)

(derivative with respect to the Wiener process itself ).

*

Geometric Brownian motion

The problem of the Wiener process was the absence of drift so it’s good to represent only the stochastic

noise. In the Black & Scholes theory of option pricing, we model an underlying asset using a stochastic

process called geometric Brownian motion which satisfies the stochastic differential equation:

* = +

* * * *

With volatility interest rate and drift ( is the Wiener process). The solution of the SDE is explicit:

, * 2

= exp ˆ‰ − ‹ + Œ , > 0

* *

2

With is the initial value (and not 0) like the value of the stock today. There is no but it

= , ∈

!

appears the (we will explain why in two lectures). The only random part is the Wiener process, the others

are deterministic information. Also in R it’s similar to code. This nice solution is possible only with the gBm.

Cox-Ingersoll-Ross (CIR) is different because the drift component has two parameters and it’s able to

represent a reversion towards the mean (mean reversion, the tendency to return to a long-term average

over time like for interest rates). Also the Vasicek and the OU have this property.

Markov property

A discrete time stochastic process is said to be Markovian (follow a Markov process) if the

{ , ≥ 1}

"

value of an underlying in any point in time conditional on the filtration (which is a way to say information

the time series of the past) is equal to the value that it would have if

= ( , , … ), [ | ]

"&$ "&$ "&2 " "&$

instead of the time series I only have the last point in time . The Markovian process doesn’t

[ | ] 1

" "&$

care about the past, we don’t need the entire time series but just the last value available. The Wiener

process and the gBm (although here the time series is used to estimate parameters so it’s not true that it

does not use the past at all) are Markovian. All examples in slide 18 are Markovian.

Failure of the Markov property

The problem with Markov processes is that with the fact that they only consider the last point in time and

therefore have no memory of the past, they are not able to capture patterns, i.e. the repetition over time

of the same shapes in different size (the so-called fractals by Benoit Mandelbrot).

Martingale

In probability theory, a martingale is a sequence of random variables (i.e., a stochastic process) for which,

at a particular time, the conditional expectation of the next value in the sequence, given all prior values, is

equal to the present value. A basic definition of a discrete-time martingale is a discrete-time stochastic

process (i.e., a sequence of random variables) that satisfies for any time

, , , … :

$ 2 9

| )

( , … , =

"#$ $ " "

This property means that there is no drift. A Wiener process is martingale (we add to 0 a standard normal

extraction which has mean 0). The gBm is a martingale only if (because otherwise it has drift):

= 0

• | )

If : it is a submartingale (drift is positive, remember). Increasing trend;

( , … , ≥

#

• | )

If : it is a supermartingale. Decreasing trend.

( , … , ≤

#

All examples in slide 18 are Martingale only if the parameters that govern their trend are zero.

1 Or also if the conditional distribution of given the past filtration equals the conditional distribution of given

" "

solely.

"#$ Simulation strategies for SDEs of two components

Simulating out of a SDE means simulating from a stochastic process. The only process we don’t need to

simulate is the gBm. There are many (although not unrelated) types of approaches in simulation for SDEs:

1. Via exact sampling (acceptance/rejection method): knowing the explicit solution for the SDE (valid

only for the gBm) using the acceptance/rejection method;

2. Via conditional distribution: I take the distribution that describes the behaviour of that SDE and I

condition it on the value of my underlying today (Markovian) getting a closed form distribution.

However we don’t arrive to a closed form unless in the case of gBm, CIR or OU;

3. Discretization of the SDE: for any SDE (gBm can use all the 3 methods!). We can’t compute the

derivative so we calculate the value of the deterministic trend and stochastic noise component

trying to evaluate the area (differentiating). This means that the next grid point is equal to the

value of today + deterministic trend + stochastic noise component. Any SDE can be rewritten as:

*#8* *#8*

= + I ( , ) + I ( , )

*#8* * 0 0 0

* *

Integrate and derivativ