Université de Montpellier

Licence 3 de Mathématiques

Modélisation stochastique

Mathieu Ribatet
Table des matières

1 Introduction 1

2 Simulation selon une loi donnée 3

2.1 Le début du voyage : la U (0, 1) . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2 La méthode d’inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.3 La méthode du rejet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.4 Algorithme du ratio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.5 Méthodes spécifiques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.6 Avant de terminer ce chapitre. . . . . . . . . . . . . . . . . . . . . . . . . . . 14

3 Intégration par Monte Carlo 15

3.1 Modèle de Black–Scholes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.2 Approche basique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.3 Variables antithétiques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.4 Variables de contrôle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.5 Échantillonnage préférentiel . . . . . . . . . . . . . . . . . . . . . . . . . . 23

i
Remarques préliminaires

Ce cours est probablement un peu en avance sur un programme classique d’un étudiant
en mathématiques appliquées. Ainsi nous ne verrons pas tous les classiques des méthodes
de Monte–Carlo mais seulement ses éléments incontournables ! Ce cours commence par
la génération de variables aléatoires puis poursuit par l’intégration par Monte–Carlo.
Pour ce deuxième partie, j’ai décidé d’illustrer la thématique avec des applications en
Mathématiques financières. Bien sûr c’est juste « une carotte » afin de rendre ce cours
plus intéressant et il ne vous ait aucunement demandé de connaître les Mathématiques
financières—ni de les apprécier d’ailleurs ;-)
L’évaluation du cours se fera par deux notes :
— une note sur un devoir sur table (avec peut être l’aide de l’ordinateur) ;
— une note sur un projet que vous ferez en binôme et pour lequel un rapport en
LATEXde 5 à 10 pages devra être rédigé.
L’évaluation portera tout autant sur la théorie que vous aurez appris tout au long
de ce cours que sur sa mise en pratique avec R. Votre engagement personnel comptera
également pour une grande partie de la note.
D’ailleurs en passant puisque vous allez utiliser le langage R (et que ce cours n’est pas
un cours sur R), je vous invite très fortement à jeter un oeil à mon cours sur R (ou tout
autre ressource qui vous plaira)

www.math.univ-montp2.fr/~ribatet/docs/LogicielsScientifiques/cours.pdf

Mais aussi d’avoir votre feuille de pompe toujours sur vous (et que je vous mets plus
bas)

http://cran.r-project.org/doc/contrib/Short-refcard.pdf

Sachez que la programmation ne s’apprend pas réellement avec un prof, c’est donc à
vous de vous y mettre. Il faut à tout prix fournir un travail régulier, n’hésitez donc pas si
vous en avez la possibilité de venir en cours avec votre ordinateur portable.

iii
 R Reference Card character or factor columns are surrounded by quotes ("); sep is the Indexing lists
field separator; eol is the end-of-line separator; na is the string for x[n] list with elements n
by Tom Short, EPRI PEAC, [email protected] 2004-11-07 missing values; use col.names=NA to add a blank column header to x[[n]] nth element of the list
Granted to the public domain. See www.Rpad.org for the source and latest get the column headers aligned correctly for spreadsheet input x[["name"]] element of the list named "name"
version. Includes material from R for Beginners by Emmanuel Paradis (with sink(file) output to file, until sink() x$name id.
permission). Most of the I/O functions have a file argument. This can often be a charac- Indexing matrices
ter string naming a file or a connection. file="" means the standard input or x[i,j] element at row i, column j
output. Connections can include files, pipes, zipped files, and R variables. x[i,] row i
On windows, the file connection can also be used with description = x[,j] column j
Getting help "clipboard". To read a table copied from Excel, use x[,c(1,3)] columns 1 and 3
x <- read.delim("clipboard") x["name",] row named "name"
Most R functions have online documentation. Indexing data frames (matrix indexing plus the following)
To write a table to the clipboard for Excel, use
help(topic) documentation on topic x[["name"]] column named "name"
write.table(x,"clipboard",sep="\t",col.names=NA)
?topic id. x$name
For database interaction, see packages RODBC, DBI, RMySQL, RPgSQL, and id.
help.search("topic") search the help system
ROracle. See packages XML, hdf5, netCDF for reading other file formats.
apropos("topic") the names of all objects in the search list matching
the regular expression ”topic” Data creation
help.start() start the HTML version of help c(...) generic function to combine arguments with the default forming a Variable conversion
str(a) display the internal *str*ucture of an R object vector; with recursive=TRUE descends through lists combining all as.array(x), as.data.frame(x), as.numeric(x),
summary(a) gives a “summary” of a, usually a statistical summary but it is elements into one vector as.logical(x), as.complex(x), as.character(x),
generic meaning it has different operations for different classes of a from:to generates a sequence; “:” has operator priority; 1:4 + 1 is “2,3,4,5” ... convert type; for a complete list, use methods(as)
ls() show objects in the search path; specify pat="pat" to search on a seq(from,to) generates a sequence by= specifies increment; length=
pattern specifies desired length Variable information
ls.str() str() for each variable in the search path seq(along=x) generates 1, 2, ..., length(along); useful for for is.na(x), is.null(x), is.array(x), is.data.frame(x),
dir() show files in the current directory loops is.numeric(x), is.complex(x), is.character(x),
methods(a) shows S3 methods of a rep(x,times) replicate x times; use each= to repeat “each” el- ... test for type; for a complete list, use methods(is)
methods(class=class(a)) lists all the methods to handle objects of ement of x each times; rep(c(1,2,3),2) is 1 2 3 1 2 3; length(x) number of elements in x
class a rep(c(1,2,3),each=2) is 1 1 2 2 3 3 dim(x) Retrieve or set the dimension of an object; dim(x) <- c(3,2)
data.frame(...) create a data frame of the named or unnamed dimnames(x) Retrieve or set the dimension names of an object
Input and output arguments; data.frame(v=1:4,ch=c("a","B","c","d"),n=10); nrow(x) number of rows; NROW(x) is the same but treats a vector as a one-
load() load the datasets written with save row matrix
shorter vectors are recycled to the length of the longest
data(x) loads specified data sets ncol(x) and NCOL(x) id. for columns
list(...) create a list of the named or unnamed arguments;
library(x) load add-on packages class(x) get or set the class of x; class(x) <- "myclass"
list(a=c(1,2),b="hi",c=3i);
read.table(file) reads a file in table format and creates a data unclass(x) remove the class attribute of x
array(x,dim=) array with data x; specify dimensions like
frame from it; the default separator sep="" is any whitespace; use attr(x,which) get or set the attribute which of x
dim=c(3,4,2); elements of x recycle if x is not long enough
header=TRUE to read the first line as a header of column names; use attributes(obj) get or set the list of attributes of obj
matrix(x,nrow=,ncol=) matrix; elements of x recycle
as.is=TRUE to prevent character vectors from being converted to fac-
factor(x,levels=) encodes a vector x as a factor
tors; use comment.char="" to prevent "#" from being interpreted as
gl(n,k,length=n*k,labels=1:n) generate levels (factors) by spec-
Data selection and manipulation
a comment; use skip=n to skip n lines before reading data; see the which.max(x) returns the index of the greatest element of x
ifying the pattern of their levels; k is the number of levels, and n is
help for options on row naming, NA treatment, and others which.min(x) returns the index of the smallest element of x
the number of replications
read.csv("filename",header=TRUE) id. but with defaults set for
expand.grid() a data frame from all combinations of the supplied vec- rev(x) reverses the elements of x
reading comma-delimited files
tors or factors sort(x) sorts the elements of x in increasing order; to sort in decreasing
read.delim("filename",header=TRUE) id. but with defaults set
rbind(...) combine arguments by rows for matrices, data frames, and order: rev(sort(x))
for reading tab-delimited files
others cut(x,breaks) divides x into intervals (factors); breaks is the number
read.fwf(file,widths,header=FALSE,sep=""  ,as.is=FALSE) of cut intervals or a vector of cut points
cbind(...) id. by columns
read a table of f ixed width f ormatted data into a ’data.frame’; widths match(x, y) returns a vector of the same length than x with the elements
is an integer vector, giving the widths of the fixed-width fields Slicing and extracting data of x which are in y (NA otherwise)
save(file,...) saves the specified objects (...) in the XDR platform- Indexing vectors which(x == a) returns a vector of the indices of x if the comparison op-
independent binary format x[n] nth element eration is true (TRUE), in this example the values of i for which x[i]
save.image(file) saves all objects x[-n] all but the nth element == a (the argument of this function must be a variable of mode logi-
cat(..., file="", sep=" ") prints the arguments after coercing to x[1:n] first n elements cal)
character; sep is the character separator between arguments x[-(1:n)] elements from n+1 to the end choose(n, k) computes the combinations of k events among n repetitions
print(a, ...) prints its arguments; generic, meaning it can have differ- x[c(1,4,2)] specific elements = n!/[(n − k)!k!]
ent methods for different objects x["name"] element named "name" na.omit(x) suppresses the observations with missing data (NA) (sup-
format(x,...) format an R object for pretty printing x[x > 3] all elements greater than 3 presses the corresponding line if x is a matrix or a data frame)
write.table(x,file="",row.names=TRUE,col.names=TRUE, x[x > 3 & x < 5] all elements between 3 and 5 na.fail(x) returns an error message if x contains at least one NA
sep=" ") prints x after converting to a data frame; if quote is TRUE, x[x %in% c("a","and","the")] elements in the given set
unique(x) if x is a vector or a data frame, returns a similar object but with fft(x) Fast Fourier Transform of an array nchar(x) number of characters
the duplicate elements suppressed mvfft(x) FFT of each column of a matrix
table(x) returns a table with the numbers of the differents values of x filter(x,filter) applies linear filtering to a univariate time series or
Dates and Times
(typically for integers or factors) to each series separately of a multivariate time series The class Date has dates without times. POSIXct has dates and times, includ-
subset(x, ...) returns a selection of x with respect to criteria (..., Many math functions have a logical parameter na.rm=FALSE to specify miss- ing time zones. Comparisons (e.g. >), seq(), and difftime() are useful.
typically comparisons: x$V1 < 10); if x is a data frame, the option ing data (NA) removal. Date also allows + and −. ?DateTimeClasses gives more information. See
select gives the variables to be kept or dropped using a minus sign also package chron.
sample(x, size) resample randomly and without replacement size ele-
Matrices as.Date(s) and as.POSIXct(s) convert to the respective class;
ments in the vector x, the option replace = TRUE allows to resample t(x) transpose format(dt) converts to a string representation. The default string
with replacement diag(x) diagonal format is “2001-02-21”. These accept a second argument to specify a
prop.table(x,margin=) table entries as fraction of marginal table %*% matrix multiplication format for conversion. Some common formats are:
solve(a,b) solves a %*% x = b for x
Math solve(a) matrix inverse of a %a, %A Abbreviated and full weekday name.
sin,cos,tan,asin,acos,atan,atan2,log,log10,exp rowsum(x) sum of rows for a matrix-like object; rowSums(x) is a faster %b, %B Abbreviated and full month name.
max(x) maximum of the elements of x version %d Day of the month (01–31).
min(x) minimum of the elements of x colsum(x), colSums(x) id. for columns %H Hours (00–23).
range(x) id. then c(min(x), max(x)) rowMeans(x) fast version of row means %I Hours (01–12).
sum(x) sum of the elements of x colMeans(x) id. for columns %j Day of year (001–366).
diff(x) lagged and iterated differences of vector x %m Month (01–12).
prod(x) product of the elements of x
Advanced data processing
%M Minute (00–59).
mean(x) mean of the elements of x apply(X,INDEX,FUN=) a vector or array or list of values obtained by
%p AM/PM indicator.
median(x) median of the elements of x applying a function FUN to margins (INDEX) of X
%S Second as decimal number (00–61).
quantile(x,probs=) sample quantiles corresponding to the given prob- lapply(X,FUN) apply FUN to each element of the list X
%U Week (00–53); the first Sunday as day 1 of week 1.
abilities (defaults to 0,.25,.5,.75,1) tapply(X,INDEX,FUN=) apply FUN to each cell of a ragged array given
%w Weekday (0–6, Sunday is 0).
weighted.mean(x, w) mean of x with weights w by X with indexes INDEX
%W Week (00–53); the first Monday as day 1 of week 1.
rank(x) ranks of the elements of x by(data,INDEX,FUN) apply FUN to data frame data subsetted by INDEX
%y Year without century (00–99). Don’t use.
var(x) or cov(x) variance of the elements of x (calculated on n − 1); if x is merge(a,b) merge two data frames by common columns or row names
%Y Year with century.
a matrix or a data frame, the variance-covariance matrix is calculated xtabs(a b,data=x) a contingency table from cross-classifying factors
%z (output only.) Offset from Greenwich; -0800 is 8 hours west of.
sd(x) standard deviation of x aggregate(x,by,FUN) splits the data frame x into subsets, computes
%Z (output only.) Time zone as a character string (empty if not available).
cor(x) correlation matrix of x if it is a matrix or a data frame (1 if x is a summary statistics for each, and returns the result in a convenient
vector) form; by is a list of grouping elements, each as long as the variables
in x Where leading zeros are shown they will be used on output but are optional
var(x, y) or cov(x, y) covariance between x and y, or between the on input. See ?strftime.
columns of x and those of y if they are matrices or data frames stack(x, ...) transform data available as separate columns in a data
cor(x, y) linear correlation between x and y, or correlation matrix if they frame or list into a single column
are matrices or data frames unstack(x, ...) inverse of stack()
round(x, n) rounds the elements of x to n decimals reshape(x, ...) reshapes a data frame between ’wide’ format with
log(x, base) computes the logarithm of x with base base repeated measurements in separate columns of the same record and Plotting
’long’ format with the repeated measurements in separate records; plot(x) plot of the values of x (on the y-axis) ordered on the x-axis
scale(x) if x is a matrix, centers and reduces the data; to center only use
use (direction=”wide”) or (direction=”long”) plot(x, y) bivariate plot of x (on the x-axis) and y (on the y-axis)
the option center=FALSE, to reduce only scale=FALSE (by default
hist(x) histogram of the frequencies of x
center=TRUE, scale=TRUE) Strings barplot(x) histogram of the values of x; use horiz=FALSE for horizontal
pmin(x,y,...) a vector which ith element is the minimum of x[i], paste(...) concatenate vectors after converting to character; sep= is the bars
y[i], . . . string to separate terms (a single space is the default); collapse= is dotchart(x) if x is a data frame, plots a Cleveland dot plot (stacked plots
pmax(x,y,...) id. for the maximum an optional string to separate “collapsed” results line-by-line and column-by-column)
cumsum(x) a vector which ith element is the sum from x[1] to x[i] substr(x,start,stop) substrings in a character vector; can also as- pie(x) circular pie-chart
cumprod(x) id. for the product sign, as substr(x, start, stop) <- value boxplot(x) “box-and-whiskers” plot
cummin(x) id. for the minimum strsplit(x,split) split x according to the substring split sunflowerplot(x, y) id. than plot() but the points with similar coor-
cummax(x) id. for the maximum grep(pattern,x) searches for matches to pattern within x; see ?regex dinates are drawn as flowers which petal number represents the num-
union(x,y), intersect(x,y), setdiff(x,y), setequal(x,y), gsub(pattern,replacement,x) replacement of matches determined ber of points
is.element(el,set) “set” functions by regular expression matching sub() is the same but only replaces stripplot(x) plot of the values of x on a line (an alternative to
Re(x) real part of a complex number the first occurrence. boxplot() for small sample sizes)
Im(x) imaginary part tolower(x) convert to lowercase coplot(x˜y | z) bivariate plot of x and y for each value or interval of
Mod(x) modulus; abs(x) is the same toupper(x) convert to uppercase values of z
Arg(x) angle in radians of the complex number match(x,table) a vector of the positions of first matches for the elements interaction.plot (f1, f2, y) if f1 and f2 are factors, plots the
Conj(x) complex conjugate of x among table means of y (on the y-axis) with respect to the values of f1 (on the
convolve(x,y) compute the several kinds of convolutions of two se- x %in% table id. but returns a logical vector x-axis) and of f2 (different curves); the option fun allows to choose
quences pmatch(x,table) partial matches for the elements of x among table the summary statistic of y (by default fun=mean)
matplot(x,y) bivariate plot of the first column of x vs. the first one of y, mtext(text, side=3, line=0, ...) adds text given by text in lty controls the type of lines, can be an integer or string (1: "solid",
the second one of x vs. the second one of y, etc. the margin specified by side (see axis() below); line specifies the 2: "dashed", 3: "dotted", 4: "dotdash", 5: "longdash", 6:
fourfoldplot(x) visualizes, with quarters of circles, the association be- line from the plotting area "twodash", or a string of up to eight characters (between "0" and
tween two dichotomous variables for different populations (x must segments(x0, y0, x1, y1) draws lines from points (x0,y0) to points "9") which specifies alternatively the length, in points or pixels, of
be an array with dim=c(2, 2, k), or a matrix with dim=c(2, 2) if (x1,y1) the drawn elements and the blanks, for example lty="44" will have
k = 1) arrows(x0, y0, x1, y1, angle= 30, code=2) id. with arrows the same effect than lty=2
assocplot(x) Cohen–Friendly graph showing the deviations from inde- at points (x0,y0) if code=2, at points (x1,y1) if code=1, or both if lwd a numeric which controls the width of lines, default 1
pendence of rows and columns in a two dimensional contingency ta- code=3; angle controls the angle from the shaft of the arrow to the mar a vector of 4 numeric values which control the space between the axes
ble edge of the arrow head and the border of the graph of the form c(bottom, left, top,
mosaicplot(x) ‘mosaic’ graph of the residuals from a log-linear regres- abline(a,b) draws a line of slope b and intercept a right), the default values are c(5.1, 4.1, 4.1, 2.1)
sion of a contingency table abline(h=y) draws a horizontal line at ordinate y mfcol a vector of the form c(nr,nc) which partitions the graphic window
pairs(x) if x is a matrix or a data frame, draws all possible bivariate plots abline(v=x) draws a vertical line at abcissa x as a matrix of nr lines and nc columns, the plots are then drawn in
between the columns of x abline(lm.obj) draws the regression line given by lm.obj columns
plot.ts(x) if x is an object of class "ts", plot of x with respect to time, x rect(x1, y1, x2, y2) draws a rectangle which left, right, bottom, and mfrow id. but the plots are drawn by row
may be multivariate but the series must have the same frequency and top limits are x1, x2, y1, and y2, respectively pch controls the type of symbol, either an integer between 1 and 25, or any
dates polygon(x, y) draws a polygon linking the points with coordinates given single character within ""
ts.plot(x) id. but if x is multivariate the series may have different dates by x and y 1 ● 2 3 4 5 6 7 8 9 10 ● 11 12 13 ● 14 15
and must have the same frequency legend(x, y, legend) adds the legend at the point (x,y) with the sym- 16 ● 17 18 19 ● 20 ● 21 ● 22 23 24 25 * * . X X a a ? ?
qqnorm(x) quantiles of x with respect to the values expected under a nor- bols given by legend ps an integer which controls the size in points of texts and symbols
mal law title() adds a title and optionally a sub-title pty a character which specifies the type of the plotting region, "s": square,
qqplot(x, y) quantiles of y with respect to the quantiles of x axis(side, vect) adds an axis at the bottom (side=1), on the left (2), "m": maximal
contour(x, y, z) contour plot (data are interpolated to draw the at the top (3), or on the right (4); vect (optional) gives the abcissa (or tck a value which specifies the length of tick-marks on the axes as a fraction
curves), x and y must be vectors and z must be a matrix so that ordinates) where tick-marks are drawn of the smallest of the width or height of the plot; if tck=1 a grid is
dim(z)=c(length(x), length(y)) (x and y may be omitted) rug(x) draws the data x on the x-axis as small vertical lines drawn
filled.contour(x, y, z) id. but the areas between the contours are locator(n, type="n", ...) returns the coordinates (x, y) after the tcl a value which specifies the length of tick-marks on the axes as a fraction
coloured, and a legend of the colours is drawn as well user has clicked n times on the plot with the mouse; also draws sym- of the height of a line of text (by default tcl=-0.5)
image(x, y, z) id. but with colours (actual data are plotted) bols (type="p") or lines (type="l") with respect to optional graphic xaxt if xaxt="n" the x-axis is set but not drawn (useful in conjonction with
persp(x, y, z) id. but in perspective (actual data are plotted) parameters (...); by default nothing is drawn (type="n") axis(side=1, ...))
stars(x) if x is a matrix or a data frame, draws a graph with segments or a yaxt if yaxt="n" the y-axis is set but not drawn (useful in conjonction with
star where each row of x is represented by a star and the columns are
Graphical parameters axis(side=2, ...))
the lengths of the segments These can be set globally with par(...); many can be passed as parameters
symbols(x, y, ...) draws, at the coordinates given by x and y, sym- to plotting commands.
bols (circles, squares, rectangles, stars, thermometres or “boxplots”) adj controls text justification (0 left-justified, 0.5 centred, 1 right-justified)
which sizes, colours . . . are specified by supplementary arguments bg specifies the colour of the background (ex. : bg="red", bg="blue", . . . Lattice (Trellis) graphics
termplot(mod.obj) plot of the (partial) effects of a regression model the list of the 657 available colours is displayed with colors())
bty controls the type of box drawn around the plot, allowed values are: "o", xyplot(y˜x) bivariate plots (with many functionalities)
(mod.obj) barchart(y˜x) histogram of the values of y with respect to those of x
The following parameters are common to many plotting functions: "l", "7", "c", "u" ou "]" (the box looks like the corresponding char-
acter); if bty="n" the box is not drawn dotplot(y˜x) Cleveland dot plot (stacked plots line-by-line and column-
add=FALSE if TRUE superposes the plot on the previous one (if it exists) by-column)
axes=TRUE if FALSE does not draw the axes and the box cex a value controlling the size of texts and symbols with respect to the de-
fault; the following parameters have the same control for numbers on densityplot(˜x) density functions plot
type="p" specifies the type of plot, "p": points, "l": lines, "b": points histogram(˜x) histogram of the frequencies of x
connected by lines, "o": id. but the lines are over the points, "h": the axes, cex.axis, the axis labels, cex.lab, the title, cex.main,
and the sub-title, cex.sub bwplot(y˜x) “box-and-whiskers” plot
vertical lines, "s": steps, the data are represented by the top of the qqmath(˜x) quantiles of x with respect to the values expected under a the-
vertical lines, "S": id. but the data are represented by the bottom of col controls the color of symbols and lines; use color names: "red", "blue"
see colors() or as "#RRGGBB"; see rgb(), hsv(), gray(), and oretical distribution
the vertical lines stripplot(y˜x) single dimension plot, x must be numeric, y may be a
xlim=, ylim= specifies the lower and upper limits of the axes, for exam- rainbow(); as for cex there are: col.axis, col.lab, col.main,
col.sub factor
ple with xlim=c(1, 10) or xlim=range(x) qq(y˜x) quantiles to compare two distributions, x must be numeric, y may
xlab=, ylab= annotates the axes, must be variables of mode character font an integer which controls the style of text (1: normal, 2: italics, 3:
bold, 4: bold italics); as for cex there are: font.axis, font.lab, be numeric, character, or factor but must have two ‘levels’
main= main title, must be a variable of mode character splom(˜x) matrix of bivariate plots
sub= sub-title (written in a smaller font) font.main, font.sub
las an integer which controls the orientation of the axis labels (0: parallel to parallel(˜x) parallel coordinates plot
Low-level plotting commands the axes, 1: horizontal, 2: perpendicular to the axes, 3: vertical) levelplot(z˜x*y|g1*g2) coloured plot of the values of z at the coor-
points(x, y) adds points (the option type= can be used) dinates given by x and y (x, y and z are all of the same length)
lines(x, y) id. but with lines wireframe(z˜x*y|g1*g2) 3d surface plot
text(x, y, labels, ...) adds text given by labels at coordi- cloud(z˜x*y|g1*g2) 3d scatter plot
nates (x,y); a typical use is: plot(x, y, type="n"); text(x, y,
names)
In the normal Lattice formula, y x|g1*g2 has combinations of optional con- rpois(n, lambda) Poisson
ditioning variables g1 and g2 plotted on separate panels. Lattice functions rweibull(n, shape, scale=1) Weibull
take many of the same arguments as base graphics plus also data= the data rcauchy(n, location=0, scale=1) Cauchy
frame for the formula variables and subset= for subsetting. Use panel= rbeta(n, shape1, shape2) beta
to define a custom panel function (see apropos("panel") and ?llines). rt(n, df) ‘Student’ (t)
Lattice functions return an object of class trellis and have to be print-ed to rf(n, df1, df2) Fisher–Snedecor (F) (χ2 )
produce the graph. Use print(xyplot(...)) inside functions where auto- rchisq(n, df) Pearson
matic printing doesn’t work. Use lattice.theme and lset to change Lattice rbinom(n, size, prob) binomial
defaults. rgeom(n, prob) geometric
rhyper(nn, m, n, k) hypergeometric
Optimization and model fitting rlogis(n, location=0, scale=1) logistic
optim(par, fn, method = c("Nelder-Mead", "BFGS", rlnorm(n, meanlog=0, sdlog=1) lognormal
"CG", "L-BFGS-B", "SANN") general-purpose optimization; rnbinom(n, size, prob) negative binomial
par is initial values, fn is function to optimize (normally minimize) runif(n, min=0, max=1) uniform
nlm(f,p) minimize function f using a Newton-type algorithm with starting rwilcox(nn, m, n), rsignrank(nn, n) Wilcoxon’s statistics
values p All these functions can be used by replacing the letter r with d, p or q to
lm(formula) fit linear models; formula is typically of the form response get, respectively, the probability density (dfunc(x, ...)), the cumulative
termA + termB + ...; use I(x*y) + I(xˆ2) for terms made of probability density (pfunc(x, ...)), and the value of quantile (qfunc(p,
nonlinear components ...), with 0 < p < 1).
glm(formula,family=) fit generalized linear models, specified by giv-
ing a symbolic description of the linear predictor and a description of
the error distribution; family is a description of the error distribution
and link function to be used in the model; see ?family Programming
nls(formula) nonlinear least-squares estimates of the nonlinear model
function( arglist ) expr function definition
parameters
return(value)
approx(x,y=) linearly interpolate given data points; x can be an xy plot-
if(cond) expr
ting structure
if(cond) cons.expr else alt.expr
spline(x,y=) cubic spline interpolation
for(var in seq) expr
loess(formula) fit a polynomial surface using local fitting
while(cond) expr
Many of the formula-based modeling functions have several common argu-
repeat expr
ments: data= the data frame for the formula variables, subset= a subset of
break
variables used in the fit, na.action= action for missing values: "na.fail",
next
"na.omit", or a function. The following generics often apply to model fitting
Use braces {} around statements
functions:
ifelse(test, yes, no) a value with the same shape as test filled
predict(fit,...) predictions from fit based on input data
with elements from either yes or no
df.residual(fit) returns the number of residual degrees of freedom
do.call(funname, args) executes a function call from the name of
coef(fit) returns the estimated coefficients (sometimes with their
the function and a list of arguments to be passed to it
standard-errors)
residuals(fit) returns the residuals
deviance(fit) returns the deviance
fitted(fit) returns the fitted values
logLik(fit) computes the logarithm of the likelihood and the number of
parameters
AIC(fit) computes the Akaike information criterion or AIC
Statistics
aov(formula) analysis of variance model
anova(fit,...) analysis of variance (or deviance) tables for one or more
fitted model objects
density(x) kernel density estimates of x
binom.test(), pairwise.t.test(), power.t.test(),
prop.test(), t.test(), ... use help.search("test")
Distributions
rnorm(n, mean=0, sd=1) Gaussian (normal)
rexp(n, rate=1) exponential
rgamma(n, shape, scale=1) gamma
Chapitre 1

Introduction

La modélisation stochastique est un terme générique dont le but principal est de

représenter un comportement par des modèles probabilistes. La plupart du temps ce que
nous cherchons à modéliser présente un aléas (telle que la précipitation par exemple) mais
ce n’est pas toujours le cas.
En ce qui nous concerne pour ce cours nous allons nous concentrer sur une petite partie
de la modélisation stochastique : la simulation et les méthodes de Monte Carlo.
Les méthodes de Monte Carlo visent à reproduire le comportement d’un processus (au
sens physique du terme) aléatoire. Cela dit les données générées par ces méthodes n’ont
(la plupart du temps) rien d’aléatoires mais son générées par des algorithmes détermi-
nistes. Les résultats pourraient donc être en principe prévisibles—si nous avions
entièrement connaissance du procédé.
Avant 1920 Utilisation de simulation afin de récupérer ou suggérer la loi théorique
(Student–t, aiguille de Buffon et π).
1940 Approximation d’intégrale pour la 1ère bombe atomique.
1950 Algorithme de Metropolis
1970 Développement d’algorithmes pour la simulation selon des lois standards
1980 Méthodes MCMC et bootstrap
1990 Méthodes MCMC largement répandues (BUGS, . . . )
2000 ABC, SMC, . . .
Petite citation de William Sealy Gosset, alias Student, en 1908 sur l’utilité de la
simulation
« Before I had succeeded in solving my problem analytically, I had endeavoured
to do so empirically. The material used was a correlation table containing the
height and left middle finger measurements of 3000 criminals, from a paper
by W. R. MacDonell (Biometrika, Vol. I., p. 219). The measurements were
written out on 3000 pieces of cardboard, which were then very thoroughly shuf-
fled and drawn at random. As each card was drawn its numbers were written
down in a book, which thus contains the measurements of 3000 criminals in a
random order. Finally, each consecutive set of 4 was taken as a sample—750
in all—and the mean, standard deviation, and correlation of each sample de-
termined. The difference between the mean of each sample and the mean of the
population was then divided by the standard deviation of the sample, giving us
the z of Section III. »
On se doute que William Gosset a dû se faire plaisir en faisant tous ces calculs certes
simples mais vraiment pas intéressants. . . Aujourd’hui nous avons l’ordinateur et il se
démotive bien moins vite que nous sur ce genre de tâches barbantes.

1
1. Introduction

L’objectif de ce cours est donc de vous montrer comment ont peu utiliser l’ordinateur
pour répondre à des problèmes bien précis. Cela dit l’ordinateur est un outil et il faut donc
savoir s’en servir à bon escient. En gros, ce n’est pas parce que l’on utilise l’ordinateur
que l’on n’aura pas besoin de maîtriser la théorie ! ! !
Remarque. Pour être parfaitement honnête, ce sera un premier pas et vous ne connaîtrez
que les méthodes basiques. Cela dit en M1 ça deviendra bien plus sérieux et séduisant ! ! !

2
Chapitre 2

Simulation selon une loi donnée

Avant de commencer ce chapitre, je souhaite reprendre l’introduction faite par Bernard

Bercu et Djalil Chafaï dans leur libre « Modélisation stochastique et simulation ».
Supposons que je joue à pile ou face avec vous mais qu’après 20 lancés, je n’ai fait
que des faces. Vous seriez très sceptique quand à l’honnêteté de ma pièce. Pourtant je
pourrais vous répondre que toutes les suites de vingt lancers sont équiprobables
et de probabilité 2−20 ! ! !
Clairement ce qui vous gêne dans la suite

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1,

c’est qu’elle semble être moins aléatoire que la suite

1 0 0 1 0 1 1 0 1 0 1 0 1 1 1 0 1 0 0 1.

Mais au fait c’est quoi une suite aléatoire ? Nous entrons clairement dans un débat
philosophique ici et nous ne perdrons pas de temps avec. Nous reviendrons juste très
rapidement sur les propriétés souhaitables de telles suites aléatoires lors de la Section 2.1.2.

2.1 Le début du voyage : la U (0, 1)

Nous n’allons pas trop insister sur la génération de variable aléatoires uniformes
sur [0, 1] mais juste parler des notions / approches élémentaires qu’ils faut absolument
connaître au moins pour votre culture générale.

2.1.1 Générateurs linéaires congruentiels

Avant de pouvoir simuler des nombres aléatoires selon une loi spécifique, il faut déjà
pouvoir le faire sur la loi de base : la loi uniforme sur l’intervalle [0, 1], notée U (0, 1).
C’est que nous allons voir dans cette section.
S’il existe différentes méthodes pour arriver à ce but nous allons nous concentrer sur
les générateurs linéaires congruentiels, i.e., dont les nombres pseudo aléatoires
sont donnés par la suite
Xn+1 ≡ aXn + b mod m, (2.1)
où les entiers a, b et m doivent être choisis minutieusement afin que la séquence générée ait
de bonnes propriétés et la valeur initiale X0 est appelée la graine aléatoire. Clairement
les valeurs X1 , X2 , . . . ∈ {0, 1, . . . , m − 1} et l’on obtient donc comme désiré des valeurs
dans [0, 1] en considérant les valeurs Xi /m.
Voici quelques exemples de générateurs :

3
2. Simulation selon une loi donnée

rand() (langage C ANSI) m = 231 , a = 1103515245 et b = 12345 ;

drand48() (langage C ANSI) m = 248 , a = 25214903917 et b = 11 ;
RANDU (qui était implémenté sur les ordinateurs IBM) m = 231 , a = 65539 et
c = 0;
Maple m = 1012 − 11, a = 427419669081 et b = 0.
Notons au passage que les deux derniers générateurs sont multiplicatifs (puisque
b = 0).
Remarque. La logiciel R que nous allons utiliser utilise par défaut le générateur Mersenne
Twister—bien que d’autres générateurs soient disponibles.
Clairement les générateurs de type (2.1) sont périodiques et de période p ≤ m. Es-
sayons de connaître un peu mieux les propriétés de ce type de générateurs mais afin de
nous faciliter la tâche, on supposera pour la suite que b = 0.
Remarque. Dans le cas multiplicatif, on ne peut clairement pas poser X0 ≡ 0 sinon Xi ≡ 0
pour tout i ≥ 0. La période des générateurs congruentiels multiplicatifs est donc au plus
m − 1.
Il est possible de construire des générateurs multiplicatifs de période maximale, i.e.,
égale à m − 1, comme le montre le théorème suivant (donné sans preuve).
Théorème 2.1.1. Une condition nécessaire pour qu’un générateur multiplicatif soit de
période m − 1 est que m soit un nombre premier.
Par le théorème de Lagrange sur les groupes, la période divise alors nécessairement
m − 1. Cela dit elle est égale à m − 1 si et seulement si a est une racine primitive de m,
i.e., a 6= 0 et
a(m−1)/p 6≡ 1 mod m,
pour tout facteur premier p de m − 1.
Exemple 2.1. Considérons le cas où m = 231 − 1, qui est un nombre premier—si si je
vous le dis. Puisque m − 1 = 2 × 32 × 7 × 11 × 31 × 151 × 331, on peut en déduire que 7
est une racine primitive de m et qu’il en est de même pour a = 75 = 16807.
Remarque. Le générateur Mersenne Twister, qui n’est pas un générateur multiplicatif, a
pour période 219937−1 ! ! !

2.1.2 Validation des générateurs

Comme déjà évoqué plus haut, nous aimerions que le générateur possède de bonnes
propriétés. S’il est clair que l’on souhaite avoir une grande période, ce n’est pas le seul
critère à prendre en compte comme le montre le générateur débile suivant
Xn+1 ≡ Xn + 1 mod m, X0 = 0,
et m arbitrairement grand.
Pour faire simple (car c’est une notion compliquée), il faut vérifier deux hypothèses :
— l’uniformité sur [0, 1] ;
— l’indépendance.
Pour tester l’uniformité (du moins marginalement) on peut utiliser des notions de
statistiques que vous verrez en M1 (probablement) appelées tests d’hypothèses (test du
χ2 , Kolmogorov–Smirnov, . . . ) Tester l’indépendance est un peu plus délicat et typique-
ment les tests sont basés sur la manière dont l’hypercube unité [0, 1]d se remplit par
le générateur—où les points remplissant l’hypercube sont de la forme (U` , . . . , U`+d−1 ),
` = 1, . . . , N − d + 1.

4
 2.2 La méthode d’inversion

2.2 La méthode d’inversion

La méthode d’inversion est une technique générique dont l’objectif est de géné-
rer un n–échantillon selon une loi donnée et ceci à partir de réalisations indépendantes
d’une U (0, 1). Cette méthode est de loin la plus simple et se base sur la fonction de
répartition inverse.

Définition 2.2.1 (Fonction de répartition et son inverse généralisée).

Soit X une variable aléatoire à valeur dans R. On appelle fonction de répartition la
fonction
F (x) = Pr(X ≤ x), x ∈ R. (2.2)
La fonction F est une fonction croissante et cadlag—cela dit pour nous elle sera souvent
continue, i.e., sans aucun point de discontinuité.
Ceci nous permet de définir son inverse généralisée, F −1 , définie par

F −1 (u) = inf {x ∈ R : F (x) ≥ u} , u ∈ [0, 1]. (2.3)

La fonction F −1 est croissante puisque F l’est et est parfois appelée fonction quan-
tile.

Lemme 2.2.1. Pour tout x ∈ R et u ∈ [0, 1], on a u ≤ F (x) ⇐⇒ F −1 (u) ≤ x.

L’implication ⇒ est triviale par la définition même

Démonstration.
d’une borne inf. Pour l’autre sens, on note que si F −1(u) ≤ x
alors, par croissance de F , on a, pour tout y ≥ x, F (y) ≥ u—faire
un dessin si nécessaire. Et puisque F est continue à droite, on en
déduit que F (x) ≥ u—en passant par la limite y ↓ x.
Remarque. Évidemment si F est bijective (et ce sera bien souvent le cas pour nous), alors
F −1 correspond à la réciproque de F .

Proposition 2.2.2. Soient F une fonction de répartition quelconque et U ∼ U (0, 1).

Alors F −1 (U ) est une variable aléatoire de fonction de répartition F .

Démonstration. Soit x ∈ R et X = F −1(U ). On a, d’après le lemme

précédent,
Pr(X ≤ x) = Pr{F −1(U ) ≤ x} = Pr{U ≤ F (x)} = F (x),
montrant le résultat annoncé.
Exemple 2.2 (Loi exponentielle).
Rappelons qu’une variable aléatoire exponentielle de paramètre λ > 0, utilisée notamment
pour modéliser les durées de survie, a pour fonction de répartition

F (x) = 1 − exp(−λx), x ≥ 0.

Calculons sa fonction réciproque (F étant clairement bijective).

ln(1 − u)
F (x) = u ⇐⇒ 1 − exp(−λx) = u ⇐⇒ x = − .
λ

5
2. Simulation selon une loi donnée

Ainsi pour générer un n–échantillon selon une Exp(λ) on pose

!
ln U1 ln Un
(X1 , . . . , Xn ) = − ,...,− ,
λ λ

d
où on a utilisé au passage que U = 1 − U .
Voici maintenant une implémentation en R.

> myrexp <- function(n, lambda)

+ return(-log(runif(n)) / lambda)
>
> ## Simulation d'une Exp(5)
> lambda <- 5
> x <- myrexp(500, lambda)
>
> ## Comparaison histogramme / densite
> hist(x, freq = FALSE)##freq = FALSE pour aire = 1
> curve(dexp(x, lambda), xlim = c(0, max(x)), col = "red", add = TRUE)

Notons au passage que pour que cette méthode soit efficace (au sens informatique),
il faut que F −1 soit rapidement calculable.
La méthode d’inversion prend une forme particulière lorsque le domaine de définition
de variable aléatoire X à simuler est un ensemble fini {x1 , . . . , xk }, k ∈ N. Sans perte de
généralité, nous pouvons supposer que x1 < x2 < · · · < xk . Ainsi on a donc



 x1 , si U ≤ p1
...






X = F −1 (U ) = xj ,
Pj−1 Pj
si `=1 p` < U ≤ `=1 p`
.
..





 Pk−1 Pk
p` < U ≤


xk , si `=1 `=1 p` = 1,

avec pj = Pr(X = xj ), j = 1, . . . , k.

Exemple 2.3 (Loi de Bernoulli).

Rappelons que la loi de Bernoulli de paramètre p ∈ (0, 1) et notée Ber(p) est définie par

p, si j = 1
Pr(X = j) =
1 − p, si j = 0.

Ainsi un n-échantillon distribué selon une Ber(p) est générée selon

 
iid
(X1 , . . . , Xn ) = 1{U1 ≤p} , . . . , 1{Un ≤p} , U1 , . . . , Un ∼ U (0, 1).

Cette méthode peut être lente car elle peut nécessiter beaucoup de test de la forme
p` < U ≤
P P
p` et il peut donc être rentable de réordonner les xj par valeurs de
probabilités décroissantes—on augmente ainsi la probabilité que U ≤ p1 par exemple.

6
 2.3 La méthode du rejet

2.3 La méthode du rejet

Cette méthode est très connue et elle doit donc faire partie intégrante de votre réper-
toire de statisticien. Elle est parfois également appelée méthode d’acceptation–rejet.
Pour cette méthode, nous nous plaçons dans le cadre où nous savons simuler des
variables aléatoires de densité de probabilité g ainsi que des U (0, 1) mais que notre but
est de générer des variables aléatoires de densité f dont le support est contenu dans
celui de g.
Remarque. J’appellerai la densité g la « densité outil »—les anglais utilisent quand à eux
l’appelation « instrumental density ».
Supposons qu’il existe une constante M > 0 telle que pour tout x ∈ R on ait
f (x)
≤ M < ∞.
g(x)
La méthode du rejet procède de la manière suivante :
1. Générer Y ∼ g et U ∼ U (0, 1) indépendemment ;
2. Poser X = Y si U ≤ f (Y )/{M g(Y )}, sinon retourner en 1.
3. Renvoyer Y .
Clairement l’algorithme sera performant dès lors que l’événement
( )
f (Y )
U≤
M g(Y )
à une forte probabilité de se réaliser. Nous reviendrons sur ce point plus tard.
Proposition 2.3.1. L’algorithme énoncé plus haut se termine en T itérations où T est
une variable aléatoire de loi géométrique de paramètre 1/M et renvoie une variable aléa-
toire de densité f .

Clairement chaque itération des étapes 1–2 est indé-

Démonstration.
pendante les unes des autres et ont toutes une probabilité de succès
 

 f (Y ) 

p = Pr  U≤ ,
 M g(Y ) 


prouvant ainsi que T suit une loi géométrique de paramètre p—que

nous allons calculer par la suite.
Soit x ∈ R, on a
 
 f (Yt) 
Pr(X ≤ x, T = t) = (1 − p)n−1 Pr 
 
Yt ≤ x, Ut ≤
 M g(Yt) 

ZZ
= (1 − p)n−1 1{y≤x,M g(y)u≤f (y)}g(y)dy du
n−1 x
Z f (y)
= (1 − p) −∞ M g(y)
g(y)dy
−1 x
Z
n−1
= (1 − p) M −∞
f (y)dy.

7
2. Simulation selon une loi donnée

Mg(x)

1.2
f(x)

1.0
0.8
Densite

0.6

●
0.4

●
●
●
●
0.2

●
●
●
●
0.0

0 1 2 3 4

x
Figure 2.1 – Illustration de la méthode du rejet via l’Exemple 2.4. Les points ont pour coor-
données {Yi , M Ui g(Yi )} et les ronds sont acceptés, les croix rejetées.

Du coup, en marginalisant, on obtient

Pr(X ≤ x) = Pr(X ≤ x, T = n)
X

n≥1
−1 x
Z
= M f (y)dy × (1 − p)n−1
X
−∞ n≥1
−1
Z x 1
=M −∞
f (y)dy
1 − (1 − p)
p Z x
= −∞
f (y)dy,
M
montrant le résultat annoncé mais aussi au passage que p = M −1
sinon f ne serait pas une densité de probabilité !

Exemple 2.4 (Loi demi–normale par une exponentielle).

La loi demi–normale de paramètre σ 2 = 1 a pour densité
q  
2
 2
π
exp − x2 , si y > 0
f (x) =
0, sinon.

8
 2.3 La méthode du rejet

Elle porte donc bien son nom !

Pour tout x > 0, on a donc , en prenant pour g une Exp(1),
q  2

f (x)
2
π
exp − x2
=
g(x) exp(−x)
s !
2 x2 − 2x
= exp −
π 2
s ( )
2 (x − 1)2 − 1
= exp −
π 2
s
2
≤ exp(1/2) ≈ 1.316.
π
La Figure 2.1 illustre cet exemple en simulant 10 réalisations indépendantes selon cette
loi demi–normale. Voici le code R associé

> f <- function(x)

+ sqrt(2 / pi) * exp(-0.5 * x^2)
>
> g <- function(x)
+ exp(-x)
>
> M <- sqrt(2 / pi) * exp(0.5)
>
> ratio <- function(x)
+ M * g(x)
>
>
> n.sim <- 10 ## nb de v.a. a generer
> i <- 1 ##nb de v.a. generees
> ans <- rep(NA, n.sim) ## stocke les v.a. generees
> while (i <= n.sim){
+ U <- runif(1)
+ Y <- abs(rnorm(1))

+ if (U <= (f(Y) / (M * g(Y)))){

+ ans[i] <- Y
+ i <- i + 1
+ }
+ }
>
> ans

[1] 1.28825688 1.25588403 1.24018084 1.46666334 0.75421121 0.03724202

[7] 0.39697197 0.46589588 1.11137043 0.18153538

Clairement d’après la Proposition 2.3.1 pour l’algorithme soit le plus efficace, il faut
que M soit le plus petit possible—afin que l’on accepte plus facilement les v.a. géné-
rées. Souvent la borne M sera appelée l’efficacité de l’algorithme car elle correspond au
nombre moyen de v.a. à générer pour en accepter une comme étant générée selon la loi
d’intérêt—nous verrons cela en exercice.

9
2. Simulation selon une loi donnée

2.4 Algorithme du ratio

Nous allons passer vite sur les détails théoriques de cette méthode car j’en ai fait
un gros exercice, cf. Feuille 1 ! Les algorithmes de type ratio sont basés sur le théorème
suivant.

h(x)dx < ∞, posons

R
Théorème 2.4.1. Pour toute fonction positive h telle que R
( s  )
v
Ch = (u, v) : 0 ≤ u ≤ h .
u

Alors Ch est d’aireR finie et si (U, V ) est uniformément distribué sur Ch alors la v.a. V /U
a pour densité h/ h.

L’idée concrète derrière ce théorème est de simuler des v.a. selon une densité f que l’on
pourrait ne connaître qu’à une constante de proportionnalité près, i.e., f ∝ h. Regardons
comment cela fonctionne à travers un exemple.

Exemple 2.5. La loi Beta(a, b), a, b > 0, a pour densité

xa−1 (1 − x)b−1 xa−1 (1 − x)b−1

f (x) = R 1 = , x ∈ [0, 1],
0 ua−1 (1 − u)b−1 du B(a, b)

où B(a, b) est la fonction (spéciale) Bêta.

Soit h(x) = x2 (1 − x)2 , 0 ≤ x ≤ 1, qui est bien sûr positive et d’aire finie, on a donc
 q 
Ch = (u, v) : 0 ≤ u ≤ h(v/u) = {(u, v) : 0 ≤ u ≤ (v/u)(1 − v/u)} .

On générera donc des v.a. Beta(3, 3)

Un peu de calcul (exo à la maison) montre que Ch ⊂ [0, 1/4] × [0, 4/27], on peut donc
simuler des U (Ch ) via la méthode du rejet en simulant d’abord des U ([0, 1/16]×[0, 4/27]).
Voici donc le code associé.

> n.sim <- 10 ##nb de va a simuler

> i <- 1 ##nb de va deja simulees
> ans <- rep(NA, n.sim) ## stocke les v.a. simulees
>
> h <- function(x)
+ x^2 * (1 - x)^2 * (x >= 0) * (x <= 1)
>
> while (i <= n.sim){
+ ## Simulation U(C_h)
+ U <- runif(1, 0, 1/4)
+ V <- runif(1, 0, 4/27)

+ if (U <= sqrt(h(V/U))){
+ ans[i] <- V / U
+ i <- i + 1
+ }
+ }

10
 2.5 Méthodes spécifiques

2.5 Méthodes spécifiques

Les méthodes énoncées plus hauts sont génériques et peuvent donc s’appliquer à un
grand nombre de situations. Ce n’est pas le cas pour les méthodes que nous allons présenter
ici. Elles méritent cependant d’être connues au moins pour votre culture générale.

2.5.1 Méthode de Box–Muller

Cette méthode permet de générer des variables selon une loi normale standard, i.e.,
une N (0, 1).
Proposition 2.5.1. Soient R ∼ Exp(1/2) et Θ ∼ U (0, 2π) mutuellement indépendantes.
Alors √ √
X = R cos Θ, Y = R sin Θ,
sont deux N (0, 1) indépendantes.

Nous allons procéder par identification et utiliser chan-

Démonstration.
gement de variable en coordonnées polaires. Soit donc une fonction
Borélienne test h. On a
2 2
 
1 x + y
h(x, y)f (x, y)dxdy = h(x, y) exp  −  dxdy

2π 2
2
 
1 r 
= h(r cos θ, r sin θ) exp  −  rdrdθ
2π 2
√ √ 1 s!   1 
  

= h( s cos θ, s sin θ)  exp − ds  dθ ,

2 2 2π
soit le produit de la densité d’une Exp(1/2) et d’une U (0, 2π).
L’algorithme de Box–Muller génère donc deux N (0, 1) indépendantes. Bien que sé-
duisant sur le papier cet algorithme a le désavantage de reposer sur les fonc-
tions trigonométriques cos et sin qui sont coûteuses en temps de calcul. On peut
toutefois contourner leur utilisation à l’aide de la méthode du rejet comme le montre
l’algorithme suivant
iid
1. Générer U1 , U2 ∼ U (0, 1) et poser V1 = 2U1 − 1 et V2 = 2U2 − 1 ;
2. Rejeter (V1 , V2 ) si (V1 , V2 ) est en dehors du disque unité, i.e., S = V12 + V22 > 1 ;
3. Poser
√
Z1 = R cos Θ = V1 −2S −1 ln S
√
Z2 = R sin Θ = V2 −2S −1 ln S

qui sont deux N (0, 1) indépendantes.

Remarque. Faisons un petit aparté puisque nous parlons de la simulation de variables
aléatoires gaussiennes. Clairement nous pouvons simuler selon une N (µ, σ 2 ) sans trop de
mal à l’aide d’une N (0, 1) puisque

Z ∼ N (0, 1) =⇒ µ + σZ ∼ N (µ, σ 2 ).

11
2. Simulation selon une loi donnée

Bien plus intéressant est la simulation de vecteurs gaussiens.

Définition 2.5.1 (Rappel).

Un vecteur aléatoire X à valeurs dans Rd , d ≥ 2, est dit gaussien si toute combinaison
linéaire de ses composantes est une variable aléatoire gaussienne.
Si X = (X1 , . . . , Xd )T est un vecteur gaussien, on définit son espérance par

E(X) = {E(X1 ), . . . , E(Xd )}T ,

et sa matrice de variance–(covariance) par

h i
Var(X) = E {X − E(X)}{X − E(X)}T

Ainsi pour simuler un vecteur gaussien de moyenne µ et de matrice de covariance

Σ—rappelons que Σ est alors nécessairement (semi) définie positive, il suffit d’utiliser
l’algorithme suivant :
iid
1. Générer le vecteur gaussien ε = (ε1 , . . . , εd )T où ε1 , . . . , εd ∼ N (0, 1) ;
2. Calculer la matrice C qui est la décomposition de Cholesky de Σ, i.e., C T C = Σ ;
3. Retourner X = µ + C T ε.

Clairement, par linéarité, X est un vecteur gaussien.

Démonstration.
De plus on a bien
 
E(X) = E µ + C ε = µ + C T E(ε) = µ
T
 
T T T
Var(X) = E C ε(C ε) = C T E(εεT )C = C T C = Σ.

2.5.2 Loi de Poisson

Rappelons qu’une variable aléatoire X de loi de Poisson de paramètre λ > 0 est définie
par
λk
Pr(X = k) = exp(−λ), k = 0, 1, . . .
k!
A priori on pourrait être tenté d’utiliser l’algorithme donné lors de la Section 2.2 pour
les variables aléatoires discrètes. C’est une mauvaise idée puisque X possède un nombre
infini d’états et l’algorithme ne peut théoriquement pas s’appliquer—numériquement oui
mais il sera très peu performant.
Une autre approche consiste à utiliser le résultat suivant.

Proposition 2.5.2. Soient E1 , E2 , . . . des v.a.i.i.d. de loi Exp(λ), λ > 0. Alors, pour tout
entier n ≥ 1, on a
n
λn X
Pr(Sn ≤ 1 ≤ Sn+1 ) = exp(−λ), Sn = E` .
n! `=1

12
 2.5 Méthodes spécifiques

Démonstration.
Z
Pr(Sn ≤ 1 ≤ Sn+1) = 1{sn≤1≤sn+1}λn+1 exp{−λ(x1 + · · · + xn+1)}dx1. . . d
Z
n+1
=λ 1{sn≤1≤sn+xn+1} exp{−λ(sn + xn+1)}dsndxn+1
n+1
Z Z ∞
=λ 1{sn≤1} exp(−λsn) 1−sn
exp(−λxn+1)dsndxn+1
Z
= λn 1{sn≤1} exp(−λsn) exp{−λ(1 − sn)}dsn
Z
n
= λ exp(−λ) 1{sn≤1}dsn
λn
= exp(−λ),
n!
où la dernière formule correspond au volume du n–simplexe (for-
mule connue ou qui se retrouve facilement).
On en déduit l’algorithme suivant pour simuler une Poisson(λ), λ > 0,
1. Posez X ← 0 et S ∼ Exp(λ) ;
2. Tant que S < 1 faire
S ← S + Exp(λ) ;
X ← X + 1;
3. Retournez X.
Remarque. Mettre un exo où l’on initialise par x ← 70 par exemple avec une Poisson(100)–
car la proba d’être inférieur à 70 est très faible. . .

2.5.3 Mélanges
En statistique on appelle un mélange une variable aléatoire qui typiquement fait in-
tervenir un nombre fini de loi de probabilité. Plus précisément, la densité d’un mélange
prend deux formes selon que le mélange soit continu ou discret
Z X
f (x) = g(x | y)p(y)dy, f (x) = pi gi (x),
Y i∈Y

où Y est le support de la loi « contrôlant » le mélange p.

L’algorithme permettant de générer des v.a. selon f est plus que naturel :
Mélange continu Mélange discret

1. Générez Y ∼ p 1. Générez Y ∼ p ;
2. Puis générez X ∼ g(· | Y ) ; 2. Puis générez X ∼ gY ;
3. Retournez X. 3. Retournez X.
Exemple 2.6. Loi de Student par mélange continu La loi de Student peut être vue comme
un mélange continu. En effet il est connu (ou du moins à savoir pour vous !) que la variable
aléatoire
ν
r
T =Z , ν > 0, X ∼ χ2ν , Z ∼ N (0, 1),
X
suit une loi de Student à ν degrés de liberté.

13
2. Simulation selon une loi donnée

Avec nos notations, ceci correspond à un mélange continu où p est la densité d’une χ2ν
et g(· | ν, Y ) celle d’une loi Normale de moyenne nulle et de variance ν/Y . L’algorithme
est donc
1. Générez Y ∼ χ2ν ;
2. Générez T ∼ N (0, ν/Y ) ;
3. Retournez T .

Remarque. mettre un exemple sur la loi négative binomiale X ∼ Neg(n, p), i.e., X | Y ∼
Poisson(Y ) avec Y ∼ Gamma(n, p).

2.6 Avant de terminer ce chapitre. . .

. . . je tiens à préciser une chose extrêmement importante. Le langage R possède
de base des fonctions pour générer des variables aléatoires selon de nombreuses loi de
probabilités usuelles. Leurs implémentations ont été soigneusement construites par des
experts reconnus mondialement et optimisées afin qu’elles soient en plus d’être fiable
statistiquement, rapide d’exécution. Il est donc parfaitement inutile d’utiliser les codes
que nous avons produits dans ce chapitre, ces derniers étaient là juste pour illustrer le
cours.
Alors pourquoi avons nous vu tout cela me direz vous ? Tout simplement car R ne
possède pas d’implémentation pour toutes les lois de probabilités et que vous aurez tôt
ou tard besoin de le faire pour des distributions exotiques. . .

14
Chapitre 3

Intégration par Monte Carlo

Tout au long de ce chapitre nous allons nous intéresser à évaluer (je dis bien évaluer
et non pas calculer) des intégrales de la forme
Z
I= h(x)dx.

Ne vous fiez pas à l’apparence simpliste de ce problème car ce problème se rencontre

extrêmement fréquemment en pratique. Alors bien sûr vous pensez tout de suite à l’inté-
gration numérique déterministe (méthodes des trapèzes / de Simpson) et vous auriez bien
raison. Sachez toutefois que ces méthodes bien qu’efficace en petites dimensions, ne le
sont plus du tout pour des dimensions plus grandes. Ainsi une intégration par Monte
Carlo pour le cas de la grande dimension sera bien plus profitable—car la vitesse de
convergence est justement indépendante de cette dernière, contrairement aux méthodes
déterministes ! ! !

3.1 Modèle de Black–Scholes

Afin de motiver ce chapitre nous allons travailler sur le modèle de Black–Scholes et de
sa formule pour établir le prix juste d’une option d’achat européen. Comme ce n’est pas
un cours de maths financières 1 nous n’irons pas trop dans les détails.
Le modèle de Black–Scholes modélise l’évolution temporelle du prix d’un actif, noté
St , à l’aide d’une équation différentielle stochastique, i.e.,

dSt = rSt dt + σSt dWt ,

où σ > 0 est la volatilité (écart type) et Wt un mouvement Brownien 2

Il est facile (mais admis car bien trop tôt pour vous) de montrer à l’aide de la Formule
d’Itô que la solution de cette équation différentielle stochastique est donnée par
( ! )
σ2
St = S0 exp r− t + σWt .
2

Comme le mouvement Brownien vous est inconnu, on supposera juste qu’à temps t > 0
fixé, Wt ∼ N (0, t2 ).
La formule de Black–Scholes établit la valeur théorique d’une option européenne.

1. et que je n’y connais pas grand chose. . .

2. Ce n’est pas grave si vous ne savez pas ce qu’est un mouvement Brownien. . .

15
3. Intégration par Monte Carlo

200
150
St

100
50
0

0 1 2 3 4 5

Annees
Figure 3.1 – Quelques trajectoires du processus issu du modèle de Black–Scholes avec S0 = 100,
r = 0.05, σ = 0.15.

Définition 3.1.1 (Option européenne).

Une option est un contrat entre un acheteur et un vendeur. Cela donne a l’acheteur de
l’option la possibilité, mais pas l’obligation, d’acheter (call) ou de vendre (put) un actif
donné (action, devises. . . ). Les prix fixés à l’avance ainsi que la durée de validité de
l’option sont définis dans le contrat.
Il existe plusieurs type d’options (européenne, asiatique, américaine, . . . ). L’option
européenne signifie que le fait de bénéficier de son droit de vente/d’achat doit se faire à
un instant T fixé par l’option. Ce n’est par exemple pas le cas d’une option américaine
qui autorise cette vente/achat à n’importe quel moment d’une période fixée par l’option,
i.e., sur [0, T ].
La formule dépend de cinq variables :
S0 la valeur actuelle de l’actif sous-jacent ;
T la date de maturité ou échéance, i.e., le temps qui reste à l’option avant son
échéance ;
K le prix d’exercice ou strike, i.e., le prix d’achat/vente de l’actif fixé par l’option ;
r le taux d’intérêt sans risque ;
σ la volatilité du prix de l’actif, i.e., son écart type annuel.
Dans le cas d’un call européen, le bénéfice (payoff in English) que je pourrais faire si
je revendais immédiatement mon actif en t = T serait bien entendu
max{St − K, 0}.

16
 3.2 Approche basique

Ainsi le « prix juste » d’un call européen, en finance on parle alors de prime de l’option,
est alors défini par
C = exp(−rT )E {max(ST − K, 0)} . (3.1)
Remarque. Le terme exp(−rT ) est une correction due aux intérêts composés. Supposons
que je place de M0 euros aujourd’hui à un taux annuel de r et que les intérêts soient versés
toutes les 1/n années, e.g., n = 12, alors au bout de T années j’aurais
nT
r

MT = M0 1 + −→ M0 exp(rT ), n → ∞.
n
Le fait de faire tendre n → ∞ correspond aux intérêts composés continûment et
convient parfaitement aux marchés financiers—hautes fréquences.
Ainsi le fait de multiplier par exp(−rT ) permet donc de « remonter le temps » afin de
déterminer le prix juste de l’option au temps t = 0.
On peut montrer (nous ne le ferons pas) que le prix juste d’un call européen est

1 S0 1 2 √
     
C = S0 Φ(d+ ) − exp(−rT )KΦ(d− ), d± = √ ln + r± σ T ,
σ T K 2

où Φ est la fonction de répartition d’une N (0, 1).

Remarque. Pour un put européen, le prix juste de l’option est alors

P = exp(−rT )E {max(K − ST , 0)} ,

de sorte que l’on a

C − P = exp(−rT )E {max(ST − K, 0) − max(K − ST , 0)}

= exp(−rT )E {max(ST − K, 0) + min(ST − K, 0)}
= exp(−rT )E(ST − K)
" ( ! )#
σ2
= exp(−rT )S0 E exp r − t + σWt − K exp(−rT )
2
( ! )
σ2
= exp(−rT )S0 exp r − T E {exp(σWt )} − K exp(−rT )
2
( ! ) !
Laplace σ2 σ2T
= exp(−rT )S0 exp r − T exp − K exp(−rT )
2 2
= S0 − K exp(−rT ).

Culture générale : Cette égalité est connue sous le nom d’arbitrage call–put.

3.2 Approche basique

L’intégration par Monte Carlo repose sur la loi (faible) des grands nombres.

Théorème 3.2.1 (Loi faible des grands nombres (rappel)).

Soient X1 , X2 , . . . sont des copies indépendantes d’une v.a. X d’espérance et de variance
finies. Alors
n
1X proba
Z
Xi −→ E(X) = xf (x)dx, n → ∞,
n i=1

17
3. Intégration par Monte Carlo

Démonstration.La preuve est directe pour peu que l’on se rappelle de

l’inégalité de Bienaymé–Tchebychev. Posons
1 Xn
X̄n = Xi , n ≥ 1.
n i=1
Par Bienaymé–Tchebychev, on a donc pour tout ε > 0,
 Var(Sn)    Var(X)
Pr |S̄n − E(Sn)| ≥ ε ≤ ⇐⇒ Pr |S̄n − E(X)| ≥ ε ≤ ,
ε2 nε2
montrant bien le résultat désiré en faisant tendre n → ∞.

Remarque. La version forte de la loi des grands nombres est bien plus puissante mais sa
démonstration bien plus compliquée également. De plus elle ne nous sera pas d’une grande
aide dans ce cours. . .
En effet bien souvent l’intégrale I peut s’écrire comme I = E(X) pour une v.a. X
de loi convenable et l’on peut donc, via la loi des grands nombres, évaluer 3 (et donc pas
calculer !) l’intégrale par

n
1X iid
Iˆ = f (Xi ), X1 , . . . , Xn ∼ X.
n i=1

Bien entendu plus n sera grand plus notre évaluation sera précise. . . Le théorème central
limite permet de connaître l’ordre de grandeur des erreurs commises.

Théorème 3.2.2 (Théorème central limite (rappel)).

Soient X1 , X2 , . . . sont des copies indépendantes d’une v.a. X d’espérance et de variance
finies. Alors
√ n1 ni=1 Xi − E(X) loi
P
n q −→ N (0, 1), n → ∞.
Var(X)

Afin de rester sur une démonstration simple, on sup-

Démonstration.
posera acquis le théorème de continuité de Lévy, i.e.,
loi
Xn −→ X ⇐⇒ E{exp(itXn)} −→ E{exp(itX)}, pour tout t ∈ R.

Posons
√ n−1 Pni=1 Xi − E(X) X n Yi Xi − E(X)
Zn = n r = √ , Yi = r ,
Var(X) i=1 n Var(X)
3. En statistiques on préfère alors dire estimer

18
 3.2 Approche basique

alors on a pour tout t ∈ R

 
 n 
E {exp(itZn)} = E exp(it Yi)
X

i=1
n
iid

−1/2
= E exp(itn Y)
 n
∞

 i 1 2 2


∼ 1 + √ tE(Y ) − t E(Y )
 n 2n 


1 2 n
 

= 1 − t 
2n
2
 
t
−→ exp −  ,

n → ∞,
2
qui est la fonction caractéristique d’une N (0, 1), prouvant ainsi le
résultat annoncé.
√
Ce résultat montre que la vitesse de convergence vers I est en n mais dépend égale-
ment de la variance de X. Clairement plus cette variance sera faible plus la précision sera
grande. . .
Remarque. Pour bien enfoncer le clou, remarquez que la vitesse de convergence ne fait au-
cunement apparaître la dimension de l’intégrale. . . (sauf peut–être au travers de Var(X)).
q
Définition 3.2.1. On appellera précision de l’estimateur Iˆ la quantité 1.96 Var(X)/n
que l’on estimera par
s Pn n
(n − 1)2 i=1 (Xi − X̄)2 1X
1.96 , X̄ = Xi .
n n i=1

Note : vous verrez plus tard dans un autre cours que cela correspond à demi longueur
d’un intervalle de confiance à 95%. . .
Revenons donc à nos moutons avec la formule de Black–Scholes (3.1) pour la prime
d’un call européen. Avec tout ce que je viens de dire, on peut l’estimer facilement par
n √
" ( ! ) #
exp(−rT ) X σ2 iid
max S0 exp r− T + σ T Zi − K, 0 , Z1 , . . . , Zn ∼ N (0, 1).
n i=1 2

Exemple 3.1 (Black–Scholes simple).

Application numérique avec S0 = 100, K = 100, σ = 0.15, r = 0.05, T = 1.

> ## Fixons la graine pour pouvoir reproduire

> ## les memes resultats plus tard
> set.seed(123)
> S0 <- 100
> K <- 100
> sigma <- 0.15
> r <- 0.05
> t <- 1

19
3. Intégration par Monte Carlo

>
> ## Calcul du prix theorique
> dplus <- (log(S0 / K) + (r + 0.5 * sigma^2) * t) / (sigma * sqrt(t))
> dminus <- (log(S0 / K) + (r - 0.5 * sigma^2) * t) / (sigma * sqrt(t))
> call.theo <- S0 * pnorm(dplus) - exp(-r * t) * K * pnorm(dminus)
>
> n.sim <- 10^4
> Z <- rnorm(n.sim, 0, 1)
> St <- S0 * exp((r - 0.5 * sigma^2) * t + sigma * sqrt(t) * Z)
>
> call <- exp(-r * t) * pmax(St - K, 0)
>
> precision <- 1.96 * sqrt(var(call) / n.sim)
> c(call.theo, mean(call), precision)

[1] 8.5916583 8.5496363 0.2181641

3.3 Variables antithétiques

Supposons que nous ayons deux estimations 4 Iˆ2 et Iˆ2 de notre fameuse intégrale I.
Ne serait il pas pertinent d’estimer alors I par I˜ = (Iˆ1 + Iˆ2 )/2 ? C’est l’idée de base des
variables antithétiques.
Si de plus, comme se sera souvent le cas, Var(Iˆ1 ) = Var(Iˆ2 ) = Var(I),ˆ alors on a

Var(I)˜ = 1 Var(I)(1
ˆ + Corr(Iˆ1 , Iˆ2 ),
2
˜ ≤ Var(I)
de sorte que la Var(I) ˆ dès lors que Corr(Iˆ1 , Iˆ2 ) ≤ 0.

Lemme 3.3.1. Soient U ∼ U (0, 1) et g une fonction monotone définie sur [0, 1]. Alors
Corr{g(U ), g(1 − U )} < 0.

On va travailler sur la covariance plutôt que sur la

Démonstration.
corrélation car cela ne change strictement rien. . . Sans perte de
généralité on peut supposer que g croissante et on notera µ =
E{g(U )}. Soient U1, U2 deux U (0, 1) indépendantes. Alors puisque
pour tout u1, u2 ∈ [0, 1],
{f (u1) − f (u2)}{f (1 − u1)f (1 − u2)} ≤ 0,
on a donc
E [{f (U1) − f (U2)}{f (1 − U1)f (1 − U2)}] ≤ 0 ⇐⇒ E {f (U1)f (1 − U1)} +
⇐⇒ E{f (U )f (1 − U )} − µ
montrant ainsi le résultat attentu.
4. Je devrais dire estimateurs mais vous ne l’avez pas encore vu alors. . .

20
 3.3 Variables antithétiques

En pratique on se sert du Lemme précédent avec F −1 (U ) et F −1 (1 − U ) qui sont alors

des v.a. négativement corrélées et de fonction de répartition F .
Remarque. Plus g est linéaire plus g(U ) et g(1 − U ) sont négativement corrélées puisque
évidemment Corr(U, 1 − U ) = −1.
Exemple 3.2 (Black–Scholes antithétique).
Du fait de la symétrie en 0 de la loi N (0, 1) on prendra comme variables antithétiques Z
et −Z. Ainsi on a l’estimation suivante
n 
1 X 
Iˆ1 + Iˆ2 ,
2n i=1

où Iˆ1 et Iˆ2 sont comme l’estimateur défini dans la Section 3 avec les Zi et −Zi respecti-
vement. Ceci s’écrit en R de la manière suivante

> ## Fixons la graine pour reproduire les memes resultats qu'avant

> set.seed(123)
> S0 <- 100
> K <- 100
> sigma <- 0.15
> r <- 0.05
> t <- 1
>
> ## Calcul du prix theorique
> dplus <- (log(S0 / K) + (r + 0.5 * sigma^2) * t) / (sigma * sqrt(t))
> dminus <- (log(S0 / K) + (r - 0.5 * sigma^2) * t) / (sigma * sqrt(t))
> call.theo <- S0 * pnorm(dplus) - exp(-r * t) * K * pnorm(dminus)
>
> n.sim <- 10^4
> Z <- rnorm(n.sim, 0, 1)
> St <- S0 * exp((r - 0.5 * sigma^2) * t + sigma * sqrt(t) * Z)
> St.anti <- S0 * exp((r - 0.5 * sigma^2) * t - sigma * sqrt(t) * Z)
>
> call <- exp(-r * t) * pmax(St - K, 0)
> call.anti <- exp(-r * t) * pmax(St.anti - K, 0)
>
> cor(call, call.anti)

[1] -0.5781445

> precision.anti <- 1.96 * sqrt(0.5 * var(call) / n.sim *

+ (1 + cor(call, call.anti)))
>
> c(call.theo, mean(0.5 * (call + call.anti)), precision.anti)

[1] 8.5916583 8.5763419 0.1001959

> precision.anti / precision * 100

[1] 45.92687

On a donc augmenté la précision de notre estimation de près de 50% juste à l’aide

d’un copié collé ! ! ! ! Des fois la fainéantise a du bon non ?

21
3. Intégration par Monte Carlo

3.4 Variables de contrôle

Le principe des variables de contrôle repose sur l’égalité toute bête suivante
I = E{f (X)} = E{f (X) − h(X)} + E{h(X)},
où, pour que la méthode ait un sens, E{h(X)} = µ soit calculable explicitement mais
aussi (et surtout !) que
Var{f (X) − h(X)} ≤ Var{f (X)}.
Définition 3.4.1. La variable h(X) est alors appelée variable de contrôle.
Ainsi on estimera I par une approche Monte–Carlo basique mais faîte sur E{f (X) −
h(X)}, i.e.,
n
1X iid
Iˆcontrole = µ + {f (Xi ) − h(Xi )} , X1 , . . . , Xn ∼ X.
n i=1
Clairement la variance de cet estimateur est alors
Var(Iˆcontrole ) = n−1 Var{f (X) − h(X)}.
Exemple 3.3 (Black–Scholes contrôle).
Nous allons utiliser l’équation d’arbitrage put–call de la Remarque 3.1 pour définir une
variable de contrôle.
En utilisant cette équation le prix juste d’un call européen peut se réécrire sous la
forme
C = S0 − K exp(−rT ) + exp(−rT )E {max(K − ST , 0)} .
On estime alors l’espérance par Monte Carlo classiquement, ce qui en R donne

> ## Fixons la graine pour reproduire les memes resultats qu'avant

> set.seed(123)
> S0 <- 100
> K <- 100
> sigma <- 0.15
> r <- 0.05
> t <- 1
>
> ## Calcul du prix theorique
> dplus <- (log(S0 / K) + (r + 0.5 * sigma^2) * t) / (sigma * sqrt(t))
> dminus <- (log(S0 / K) + (r - 0.5 * sigma^2) * t) / (sigma * sqrt(t))
> call.theo <- S0 * pnorm(dplus) - exp(-r * t) * K * pnorm(dminus)
>
> n.sim <- 10^4
> Z <- rnorm(n.sim, 0, 1)
> St <- S0 * exp((r - 0.5 * sigma^2) * t + sigma * sqrt(t) * Z)
>
>
> put <- exp(-r * t) * pmax(K - St, 0)
> call <- put + S0 - K * exp(-r * t)
>
> precision.control <- 1.96 * sqrt(var(put) / n.sim)
>
> c(call.theo, mean(call), precision.control)

22
 3.5 Échantillonnage préférentiel

[1] 8.5916583 8.5878936 0.1235383

> precision.control / precision * 100

[1] 56.62631

Ce qui nous fait une précision accrue d’environ 40% par rapport à la première approche
mais cependant une performance moindre par rapport à l’approche antithétique.
Cela dit on peut tout à fait combiner les deux ! ! ! Mais le calcul des erreurs sera alors
plus complexe (et non mis)

> St.anti <- S0 * exp((r - 0.5 * sigma^2) * t - sigma * sqrt(t) * Z)

> put.anti <- exp(-r * t) * pmax(K - St.anti, 0)
> call.anti <- put.anti + S0 - K * exp(-r * t)
>
> cor(put, put.anti)

[1] -0.3475422

> precision.final <- 1.96 * sqrt(0.5 * var(put) / n.sim *

+ (1 + cor(put, put.anti)))
>
> c(call.theo, mean(0.5 * (call + call.anti)), precision.final)

[1] 8.59165831 8.64387624 0.07200217

> precision.final / precision * 100

[1] 32.80068

3.5 Échantillonnage préférentiel

L’échantillonnage préférentiel (importance sampling in English) est une technique qu’il
faut à tout prix avoir dans sa trousse à outils. Bien que l’idée soit très simple, elle peut
être d’une efficacité redoutable.
Nous allons motiver son utilisation par un exemple simpliste mais instructif. Supposons
que nous soyons intéressés en l’évaluation de la fonction de répartition d’une N (0, 1) en
un point quelconque x ∈ R, i.e.,
!
Z x
y2
I(x) = (2π) exp − dy.
−∞ 2
ˆ iid
= n−1 ni=1 1{Xi ≤x} avec X1 , . . . , Xn ∼ N (0, 1)
P
Clairement on pourrait l’estimer par I(x)
et sa variance vérifie, comme somme de Bernoulli,

ˆ I(x){1 − I(x)} 1
Var{I(x)} = ≤ .
n 4n
Nous en concluons que pour avoir une précision au moins égale à 10−4 il faut (approxi-
mativement)
1
2 √ ≤ 10−4 ⇐⇒ n ≥ 108 ,
4n

23
3. Intégration par Monte Carlo

impliquant donc que cette approche est numériquement infaisable dès lors que l’on cherche
à estimer de très faibles probabilités 5 .
Essayons donc de faire mieux. . . Ainsi plutôt que de faire du Monte Carlo basique, i.e.,
Z Z
I= h(x)dx = m(x)f (x)dx, f densité,

nous écrivons Z
f (x)
I= m(x) dx,
g(x)
où g est une nouvelle densité dont le support contient celui de m × f .
Remarque. Parfois le rapport w(x) = f (x)/g(x) est appelé poids d’importance ou
poids préférentiels.
Ainsi en utilisant la loi des grands nombre comme précédemment, nous pouvons évaluer
l’intégrale d’une nouvelle manière
n
1X iid
Iˆpref = m(Xi )w(Xi ), X1 , . . . , Xn ∼ g.
n i=1

Proposition 3.5.1. Sous les conditions décrites juste avant, on a

(Z )
2 f (x)
E(Iˆpref ) = I, Var(Iˆpref ) = n −1
m(x) g(x)dx − I 2
.
g(x)

Démonstration. Pour l’espérance

1 Xn
E(Iˆpref ) =
Z Z
E {m(Xi)w(Xi)} = m(x)w(x)g(x)dx = m(x)f (x)dx =
n i=1
Et pour la variance

Var(Iˆpref ) = n−1Var{m(X)w(X)}
   
−1 2 2
=n E {m(X)w(X)} − I
Z 
−1 2 2 2
=n m (x)w (x)g(x)dx − I
 
 f (x) 
= n−1 
Z
2 2

m (x) f (x)dx − I 
 g(x) 

On peut montrer (cf. Feuille 2) que la variance de l’estimateur par échantillonnage

préférentiel est nulle lorsque g(x) = m(x)f (x)/I ce qui est bien entendu d’une inutilité
la plus totale puisque nous cherchons justement à estimer I ! ! ! Cela dit ce résultat nous
suggère qu’un bon choix pour g devrait vérifier g(x) ∝ m(x)f (x).

5. Ce n’est pas juste un exemple théorique, pour l’étape de pré dimensionnement des barrages par
exemple, la loi impose que l’ouvrage résiste à des événements de probabilités inférieures à 10−3 !

24
 3.5 Échantillonnage préférentiel

Exemple 3.4 (Black–Scholes avec une option peu « jouable »).

Rappelons qu’un call européen donne la possibilité droit de vendre un actif au temps T
aux prix ST alors que nous l’avons acheté au prix K. Clairement le call européen est
rentable dès lors que ST > K et le détenteur de l’option se fait un bénéfice de ST − K—je
néglige ici le prix de l’option.
Ainsi l’option a peu de chance d’être utilisée si Pr(ST > K) est petite et nous nous
ramenons donc à la situation expliquée en début de cette Section. C’est notamment le cas
lorsque S0  K.
Plutôt que d’échantillonner selon une N (0, 1), prenons pour loi d’importance une
N (m, 1) pour un m ∈ R bien choisi. On a donc
n √
" ( ! ) # !
exp(−rT ) X σ2 X 2 − (Xi − m)2
Ĉpref = max S0 exp r− T + σ T Xi − K, 0 exp − i
n i=1 2 2
n √
" ( ! ) # !
exp(−rT ) X σ2 2mXi − m2
= max S0 exp r − T + σ T Xi − K, 0 exp − ,
n i=1 2 2
iid
où X1 , . . . , Xn ∼ N (m, 1). Ce qui en R donne

> ## Fixons la graine pour reproduire les memes resultats qu'avant

> set.seed(123)
> S0 <- 50
> K <- 100
> sigma <- 0.15
> r <- 0.05
> t <- 1
>
> ## Calcul du prix theorique
> dplus <- (log(S0 / K) + (r + 0.5 * sigma^2) * t) / (sigma * sqrt(t))
> dminus <- (log(S0 / K) + (r - 0.5 * sigma^2) * t) / (sigma * sqrt(t))
> call.theo <- S0 * pnorm(dplus) - exp(-r * t) * K * pnorm(dminus)
>
> n.sim <- 10^4
> Z <- rnorm(n.sim, 0, 1)
>
> ## Monte-Carlo basique
> St <- S0 * exp((r - 0.5 * sigma^2) * t + sigma * sqrt(t) * Z)
> call <- exp(-r * t) * pmax(St - K, 0)
>
> ## Echantillonnage preferentiel
> m <- log(K / S0) / sigma
> X <- m + Z
>
> St.pref <- S0 * exp((r - 0.5 * sigma^2) * t + sigma * sqrt(t) * X)
> call.pref <- exp(-r * t) * pmax(St.pref - K, 0)
> weights <- exp(- m * X + 0.5 * m^2)
>
> c(call.theo, mean(call), mean(call.pref * weights))

[1] 1.983208e-05 0.000000e+00 2.005198e-05

25
Conclusion

Voilà ce cours est maintenant terminé. Clairement certains aspects manque cruellement
mais les inclure aurait été bien trop ambitieux pour un cours de 3ème année de Licence.
En effet un cours sur les méthodes de Monte–Carlo se donne généralement en 1ère année
de Master. . . Parmi les aspects manquants je pense notamment à
— Théorème de Rao–Blackwell pour réduire la variance
— Les Chaînes de Markov par Monte–Carlo
— Bootstrap, Jackknife
Les plus curieux (et motivés) pourront donc jeter un oeil sur ces thémes !

27