# Inference and autoregressive processes

Consider a (stationary) autoregressive process, say of order 2,

for some white noise  with variance . Here is a code to generate such a process,

> phi1=.5
> phi2=-.4
> sigma=1.5
> set.seed(1)
> n=240
> WN=rnorm(n,sd=sigma)
> Z=rep(NA,n)
> Z[1:2]=rnorm(2,0,1)
> for(t in 3:n){Z[t]=phi1*Z[t-1]+phi2*Z[t-2]+WN[t]}

Here, we have to estimate two sets of parameters: the autoregressive coefficients, and the variance of the innovation process . There are (at least) three techniques to estimate those parameters.

• using least square regression

A natural idea is to see here a regression model, and thus, if we consider a matrix formulation,

Here we can run (conditional) ordinary least squares estimation,

> base=data.frame(Y=Z[3:n],X1=Z[2:(n-1)],X2=Z[1:(n-2)])
> regression=lm(Y~0+X1+X2,data=base)
> summary(regression)

Call:
lm(formula = Y ~ 0 + X1 + X2, data = base)

Residuals:
Min      1Q  Median      3Q     Max
-4.3491 -0.8890 -0.0762  0.9601  3.6105

Coefficients:
Estimate Std. Error t value Pr(>|t|)
X1  0.45107    0.05924   7.615 6.34e-13 ***
X2 -0.41454    0.05924  -6.998 2.67e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.449 on 236 degrees of freedom
Multiple R-squared: 0.2561,	Adjusted R-squared: 0.2497
F-statistic: 40.61 on 2 and 236 DF,  p-value: 6.949e-16

> regression\$coefficients
X1         X2
0.4510703 -0.4145365
> summary(regression)\$sigma
[1] 1.449276
• using Yule-Walker equations

As we’ve seen in class, we can easily get the following equations for the autocovariance functions,

which can also be written

So we just have to solve a simple linear system of equations. Note that if we divide by the variance, those equations can be written in terms of the autocorrelation functions

The code is the following

> rho1=cor(Z[1:(n-1)],Z[2:n])
> rho2=cor(Z[1:(n-2)],Z[3:n])
> A=matrix(c(1,rho1,rho1,1),2,2)
> b=matrix(c(rho1,rho2),2,1)
> (PHI=solve(A,b))
[,1]
[1,]  0.4517579
[2,] -0.4155920

Now, we need to extract the estimated innovation process, from this set of parameters (note that it could be possible to include the variance term in Yule-Walker equations, to get a three dimensional linear equation)

> estWN=base\$Y-(PHI[1]*base\$X1+PHI[2]*base\$X2)
> sd(estWN)
[1] 1.445706

This estimator is probably not the best one (we can take into account that we’ve lost two degrees of freedom), but as a starting point, let us consider this one.

• using (conditional) likelihood estimators

Finally, we can assume some distribution for the innovation process. Thestandard model is a Gaussian model, i.e.

In that case, the conditional log likelihood (conditional since we set the first two observations here) is

> CondLogLik=function(A,TS){
+ phi1=A[1];  phi2=A[2]
+ sigma=A[3]	; L=0
+ for(t in 3:length(TS)){
+ L=L+dnorm(TS[t],mean=phi1*TS[t-1]+
+ phi2*TS[t-2],sd=sigma,log=TRUE)}
+ return(-L)}

Now, we can run standard optimization procedures,

> LogL=function(A) CondLogLik(A,TS=Z)
> optim(c(0,0,1),LogL)
\$par
[1]  0.4509685 -0.4144938  1.4430930

\$value
[1] 425.0164

\$counts
88       NA

\$convergence
[1] 0

\$message
NULL

Here, our three estimators are rather close. Actually, if we generate 1,000 time series (of size 240), those are the Box-plots of our three estimators, for the first order autoregressive coefficient

for the second one,

and finally for the standard deviation of the innovation process

All those estimators behave nicely, and are rather close. Note that they all might be biased, but they are consistent (see Davidson and MacKinnon for instance, in their book, for more details).

# In statistics, having too much information might not be a good thing

A common idea in statistics is that if we don’t know something, and we use anestimator of that something (instead of the true value) then there will be some additional uncertainty. For instance, consider a random sample, i.i.d., from a Gaussian distribution. Then, a confidence interval for the mean is

where  is the quantile of probability level  of the standard normal distribution . But usually, standard deviation  (the something is was talking about earlier) is usually unknown. So we substitute an estimation of the standard deviation, e.g.

and the cost we have to pay is that the new confidence interval is

where now  is the quantile of the Student distribution, of probability level , with  degrees of freedom.
We call it a cost since the new confidence interval is now larger (the Student distribution has higher upper-quantiles than the Gaussian distribution).
So usually, if we substitute an estimation to the true value, there is a price to pay.
A few years ago, with Jean David Fermanian and Olivier Scaillet, we were writing a survey on copula density estimation (using kernels,  here). At the end, we wanted to add a small paragraph on the fact that we assumed that we wanted to fit a copula on a sample  i.i.d. with distribution , a copula, but in practice, we start from a sample with joint distribution  (assumed to have continuous margins, and – unique – copula ). But since margins are usually unknown, there should be a price for not observing them.
To be more formal, in a perfect wold, we would consider

but in the real world, we have to consider

where it is standard to consider ranks, i.e.  are empirical cumulative distribution functions.
My point is that when I ran simulations for the survey (the idea was more to give illustrations of several techniques of estimation, rather than proofs of technical theorems) we observed that the price to pay… was negative ! I.e. the variance of the estimator of the density (wherever on the unit square) was smaller on the pseudo sample  than on perfect sample .
By that time, we could not understand why we got that counter-intuitive result: even if we do know the true distribution, it is better not to use it, and to use instead a nonparametric estimator. Our interpretation was based on the discrepancy concept and was related to the latin hypercube construction:

With ranks, the data are more regular, and marginal distributions are exactlyuniform on the unit interval. So there is less variance.
This was our heuristic interpretation.
A couple of weeks ago, Christian Genest and Johan Segers proved that intuition in an article published in JMVA,

Well, we observed something for finite , but Christian and Johan obtained an analytical result. Hence, if we denote

the empirical copula in the perfect world (with known margins) and

the one constructed from the pseudo sample, they obtained that, everywhere

with nice graphs of ,

So I was very happy last week when Christian show me their results, to learn that our intuition was correct. Nevertheless, it is still a very counter-intuitive result…. If anyone has seen similar things, I’d be glad to hear about it !

# Estimateurs des moments

Dans le cours, nous avions vu que si , un estimateur naturel pour le paramètre était . On appelle cette méthode la méthode des moments. Mais il est possible, que, plus généralement , et on pourrait considérer comme estimateur naturel

On avait noté en cours que le premier soucis était un soucis d’inversibilité. Pour utiliser cette méthode, il faut une relation bijective entre le moment et le paramètre. On avait vu avec la loi binomiale – que ce n’est pas le cas pour la variance, par exemple.

Le second soucis est un problème numérique: il faut inverser cette fonction (ce qui peut ne pas être trivial) mais avec une méthode de type sécante (programméeici), ça pourrait marcher. Par exemple pour la loi de Poisson, on sait que

(…etc). On peut numériquement construire des estimateurs de la manière suivante

secant=function(fun, x0, x1, tolerence=1e-07, niter=500){
for ( i in 1:niter ) {
x2 <- x1-fun(x1)*(x1-x0)/(fun(x1)-fun(x0))
if (abs(fun(x2)) < tolerence)
return(x2)
x0 <- x1
x1 <- x2
}}
X=rpois(50,10)
m1=mean(X)
f=function(x){x*(1+x)-mean(X^2)}
secant(f,0,40)
m2=secant(f,0,50)
f=function(x){x*(1+3*x+x^2)-mean(X^3)}
secant(f,0,50)
m3=secant(f,0,50)
f=function(x){x*(1+7*x+6*x^2+x^3)-mean(X^4)}
secant(f,0,100)
m4=secant(f,0,20)

Le dernier soucis est un problème plus statistique: il semble que l’estimateur basé sur le premier moment soit meilleur que les autres. Sur le dessin ci-dessous, on a simulé 10,000 fois des lois de Poisson de paramètre 10 (qui est la vraie valeur du paramètre, ici représentée verticalement en pointillés), et on représente en rouge l’estimateur construit sur le premier moment, avec la croix au centre qui représente la valeur moyenne, et le segment qui est un intervalle de confiance à 80% de notre estimateur (i.e. les quantiles à 10% et 90% obtenus sur les 10,000 simulations),

En bleu, on a l’estimateur construit à l’aide du moment d’ordre deux, en vert du moment d’ordre trois, et en mauve du moment d’ordre quatre. Bref, plus on prend des moments d’ordre élevé, plus le biais semble important, et plus la variance de notre estimateur augmente.

Moralité, l’estimateur de la méthode des moments (et l’unique) sera construit sur les moments les plus faibles possibles…

# Qui a eu l’idée d’élargir les trottoirs en France ?

J’étais tombé l’autre jour sur un joli dessin expliquant pourquoi les policiers et les syndicats n’estimaient jamais le même nombre de manifestants, avec une explication scientifique, et vaguement rationnelle, ici.

C’est  effectivement joli… et on retrouve un rapport de 3 dont pas mal de monde semble parler, sauf que si c’était vrai, le rapport entre les estimations de la police et des syndicats devrait être stable d’une manifestation à l’autre, non ?
Si on compare le 19 mars dernier (ici) et le 16 octobre (), on arrive aux chiffres suivants,

Les traits pointillés sont la droite de régression. La courbe en trait plein et la régression (passant par l’origine car on cherche un rapport de proportionnalité) en prenant en compte un effet taille, i.e. en introduisant des poids différents dans ma régression. J’ai pris la taille des villes, autrement dit, plus la ville est grosse, plus je donne un poids important au rapport… L’idée étant d’intégrer un effet médiatique, même si la pente n’est plus le ratio moyen…  Avec cette dernière méthode, j’obtiens des pentes de 4,06 pour mars, et 6,12 pour octobre… ce qui est à un peu différents des rapports nationaux mentionnés sur le site du monde, i.e. (3 millions selon les syndicats et 1,2 selon la police, en mars, contre  2,5 millions contre 825 milliers l’autre jour). Notons qu’on peut vérifier rapidement que ce n’est pas un problème statistique, les valeurs sont sensiblement différentes, si l’on compare les distributions des estimateurs..

Bref, quelque chose semble avoir changé entre ces deux dates…. la taille des trottoirs ?

# Propriétés des estimateurs dans une régression

Comme les rappels (je devrais plutôt dire “remise à niveau“) ont été particulièrement rapides, je vais prendre un peu de temps pour revenir sur quelques éléments du cours d’économétrie (d’autant plus que j’ai reçu quelques questions par mail).

• sur le calcul des estimateurs

Je vais reprendre les questions que j’ai reçu, ça sera plus simple, “je comprend pas très bien quelle méthode utilise l’opération lm(Y~X…..) sur R pour trouver les coef” ou encore “je regresse un modèle simple, je sors le summary, et j’obtiens un betaREG 3 étoiles. Toujours sur R, je compare ce betaREG avec le betaCalc = (X’X)-1 * X’Y, il est différent ! Comment expliquer cela ? Est-ce possible ?“.
Commençons par la première question: la fonction lm() utilise la méthode dite des moindres carrés, ce qui revient à calculer

en adoptant une écriture matricielle, i.e. la solution est alors

Cet estimateur est celui qui minimise la somme des carrés des erreurs. Et on peut vérifier numériquement que cet estimateur est bien celui calculé par R,

>data(cars)
>X=cbind(rep(1,nrow(cars)),cars\$speed); Y=cars\$dist
> X[1:5,]
[,1] [,2]
[1,]    1    4
[2,]    1    4
[3,]    1    7
[4,]    1    7
[5,]    1    8
> solve(t(X)%*%X)
[,1]         [,2]
[1,]  0.19310949 -0.011240876
[2,] -0.01124088  0.000729927
> t(X)%*%Y
[,1]
[1,]  2149
[2,] 38482
> (solve(t(X)%*%X))%*%(t(X)%*%Y)
[,1]
[1,] -17.579095
[2,]   3.932409

ce qui correspond très précisément à ce que calcule la fonction de R

> lm(dist~speed,cars)
Coefficients:
(Intercept)        speed
-17.579        3.932

Et on peut aussi vérifier qu’il minimise bien la somme des carrés des erreurs,

> b0=seq(-30,10,by=3);
> b1=seq(1,7,by=.5)
> SC=matrix(NA,length(b0),length(b1))
> for(i in 1:length(b0)){
> for(j in 1:length(b1)){
+ SC[i,j]=sum((cars\$dist-(b0[i]+b1[j]*cars\$speed))^2)
}}
> contour(b0,b1,SC)

et je peux certifier que la valeur obtenue est bien le minimum de cette fonction (en fait c’est un résultat théorique d’optimisation).
Ca se généralise bien entendu en dimension plus grande. Alors attention, cette formule permet de calculer un estimateur du vecteur des paramètres, mais la formule est bien sûre fausse si on raisonne composante par composante.

> X=cars\$speed; Y=cars\$dist
> (solve(t(X)%*%X))%*%(t(X)%*%Y)
[,1]
[1,] 2.909132

Mais ça c’est un résultat (trivial) d’algèbre linéaire. En effet, le fait d’avoir

ne signifie pas que l’on ait ce genre de formule composante par composante, i.e.

Dit de manière un peu formelle, le produit matriciel n’est pas un produit terme à terme. La traduction économétrique de cette idée est qu’une régression multiple n’est pas une succession de régressions simples. En effet, dans l’équation précédante, le dernier terme correspond à la régression sans la constante, comme on peut le voir ci-dessous

> lm(dist~0+speed,cars)
Coefficients:
speed
2.909

et le premier est tout simplement la moyenne empirique,

> lm(dist~1,cars)
Coefficients:
(Intercept)
42.98
> mean(cars\$dist)
[1] 42.98
> X=rep(1,nrow(cars)); Y=cars\$dist
> (solve(t(X)%*%X))%*%(t(X)%*%Y)
[,1]
[1,] 42.98
• sur les propriétés des estimateurs

Une question était “doit-on vérifier les propriétés de l’estimateur (sans biais, efficace etc…) du summary de la Reg ? pour chaque régression ?
Bon, la réponse est non car c’est impossible ! Pour répéter tranquillement ce que j’avais dit plusieurs fois oralement: quand on fait un modèle (théorique), on part d’hypothèses (par exemple ici les résidus sont centrées, de variance constante). Sous ces hypothèses on a des propriétés: ce sont des théorèmes, autrement dit de la théorie. Le plus connu étant le théorème de Gauss-Markov sur le modèle linéaire, qui garantie que l’estimateur par moindre carré est BLUE (best linear unbiased estimator).
Autrement dit (si je retraduit ce que dit ce théorème), si les hypothèses sont valides, alors la théorie nous garantie que les estimateurs vérifient des propriétés, dont celle d’être sans biais, par exemple.
On ne
peut pas
vérifier qu’un estimateur est sans biais, on peut juste vérifier ex post que les hypothèses sont valides (ou non), ce qui garantit (ou pas) l’absence de biais.
Mais comme on va le voir par la suite, le biais ce n’est pas forcément gênant, ce qui est important, c’est surtout la convergence. Comme le notait Clive Granger  “if you can’t even get a consistant estimator, you shouldn’t be in this business“.

• et sur des points plus théoriques

sinon la personne qui avait des soucis numériques à calculer les coefficients me demandait si cela avait un “rapport avec les plim X’E/n et plim (X’X)-1/ n“.
Pour revenir sur ce point, une hypothèse forte dans le modèle linéaire (de base) est que le bruit soit soit un vrai bruit, c’est à dire non corrélé avec la variable explicative, i.e.

En effet, sinon

autrement dit l’estimateur par moindre carré est biasé. La solution la plus naturelle pour éviter ce “problème” est d’utiliser un instrument (ou une variable instrumentale).
Mais surtout, sur cet exemple, on voit qu’on ne pourra avoir la convergence de notre estimateur par moindres carrés que si la corrélation entre le bruit et la variable explicative tend vers 0 asymptotiquement.
Il existe en effet une différence fondamentale entre le biais et la consistance (ou la convergence vers la vraie valeur).

• un estimateur sans biais est généralement convergent (et de manière consistante)
• un estimateur biaisé peut ne pas converger

Prenons un petit exemple pour illustrer ces points. On supposera n1 < n2 < n3par la suite. Comme je l’ai expliqué auparavant, la distribution de l’estimateur est une propriété théorique qui ne peut se voir sur un jeu de données (à moins de faire des simulations mais ça sort du cadre de mon billet d’aujourd’hui). La figure ci-dessous correspond au cas sans biais, et convergent

La figure ci-dessous correspond au cas biaisé, mais convergent, et consistant, au sens où asymptotiquement, l’estimateur sera sans biais.

Enfin, un dernier cas correspondant au cas biaisé, convergent, mais non consistant. Asymptotiquement l’estimateur sera toujours biaisé,

Mais revenons un peu sur la formaliation sous-jacente au modèle linéaire. Les propriétés sur les résidus sont conditionnelles aux variables explicatives, en particulier.

(c’est l’hypothèse d’homoscédasticité et d’indépendance des résidus). Bon, il faut aussi l’absence de relation linéaire entre les variables explicatives, ce qui se traduit parfois comme une absence de collinéarité.
On peut s’intéresser aux propriétés asymptotiques de l’estimateur. Pour cela, il faut peut être rappeler ce que signifie la convergence pour une variable aléatoire (en statistique mathématique, les estimateurs sont vus comme des variables aléatoires). On dira que  converge en probabilité vers , parfois noté

si

Par exemple la loi des grands nombres garantie que, pour un échantillon i.i.d.

si l’espérance des variables est finie, alors

Alors sous les hypothèses mentionnées auparavant, on peut montrer que, dans le modèle linéaire,

On parlera alors d’estimateur convergent. On peut aussi montrer que

Mais formellement,

si

Voilà en gros quelques éléments pour répondre aux questions, mais tout ça sera repris en détails dans le cours d’économétrie 1 qui commencera très bientôt… et surtout d’économétrie 2 qui abordera les aspects dynamiques.

# Delta method avec R, ou transformation des paramètres

Un petit billet aujourd’hui sur une question qui m’a été posée (et qui m’a permis de découvrir des nouvelles fonctions sous R). On suppose que l’on a un modèle économétrique de la forme

que l’on peut encore écrire

Supposons que ce soit la formulation en terme de ces derniers coefficients qui nous intéresse, et que l’on cherche un intervalle de confiance.
La Delta Method nous dit que si on dispose d’un estimateur asymptotiquement gaussien (mais c’est le cas pour la plupart des estimateurs que l’on construit en pratique), tout transformation bijective de cet estimateur reste asymptotiquement Gaussien. Bref, on va pouvoir utiliser ces résultats pour en déduire un intervalle de confiance (asymptotique) pour le nouveau paramètre. En effet, l’estimateur par moindre carrés vérifie la normalité asymptotique (sous quelques hypothèses)

ce qui donne avec les notations précédentes

Autrement dit, on peut facilement en déduire un intervalle de confiance (asymptotique). Pour cela, on a trois méthodes (au moins) avec R. Pour la suite, je vais considerer la régression sur la base de données cars (comme toujours), et on s’intéresse à l’estimation de

(ne me demandez pas pourquoi, il me fallait un exemple numérique, j’ai pris un truc nonlinéaire pour qu’au moins ça puisse avoir un vague intérêt).

• En faisant un peu de calcul formel

Pour les calculs symboliques, il y a la fonction D dans le package de base,

> f <- expression(sin(cos(x + y^2)))
> ( Dfx <- D(f, "x") )
-(cos(cos(x + y^2)) * sin(x + y^2))

ce qui pourra alors être utilisé pour faire les calculs de la variance asymptotique.

> x=c(0,1,2); y=1
> eval(D.sc)
[1] -0.7216061 -0.8316919 -0.0774320

On peut également utiliser cette fonction pour faire plein de jolis dessins,

> x <-seq(-5, 5, length=100)
> fx=x^2
> Df <-D(expression(x^2), "x")
> Dfx<-eval(Df)
> plot(x, fx, type = "l", las = 1,ylim=c(-2,2))
> for (i in 1:length(x)) {
+ abline(a = fx[i] - Dfx[i] * x[i], b = Dfx[i],
+ col = rainbow(100)[i])
+ }

Mais bon, c’était juste pour rigoler un peu… On peut utiliser cette fonction pour calculer l’écart type de l’estimateur de la grandeur qui nous intéresse.

> reg1 = lm(dist~speed, data = cars)
> estmean = coef(reg1)
> estvar = summary(reg1)\$cov.unscaled * summary(reg1)\$sigma^2
> f = function(x,y){
+ -y/(x+y)
+ }
> X=f(estmean[1],estmean[2])
> g=expression(-y/(x+y))
> x<-as.numeric(estmean[1]); y<-as.numeric(estmean[2])
> V=eval(c(eval(D(g, "x")),eval(D(g, "y"))))
> sqrt(t(V) %*% estvar %*% V)
[,1]
[1,] 0.1063255

Voilà pour la première méthode. Petite remarque: pour faire le calcul symbolique de la dérivée, j’utilise expression et non pas function.

• En utilisant la fonction deltamethod de la librarie msm

Ici on commence comme avant, sauf que la phase de calcul est immédiat…

> reg2 <- lm(dist~speed, data = cars)
> estmean = coef(reg2)
> estvar = summary(reg2)\$cov.unscaled * summary(reg2)\$sigma^2
> deltamethod (~ -x2 / (x1 + x2), estmean, estvar)
[1] 0.1063255

et d’ailleurs, ça peut s’appliquer pour un paramètre en dimension 2

> deltamethod (list(~ -x2 / (x1 + x2), ~ x1+x2), estmean, estvar)
[1] 0.1063255 6.3664368

En revanche, on obtient des estimateur des écart-types, et pas de la matrice de variance-covariance des estimateurs. Autrement dit on ne peut pas en déduire l’ellipse de confiance, seulement des intervalles de confiance marginaux (en utilisant la normalité asymptotique).

• En utilisant la fonction delta.method de librarie alr3

Pour finir, une dernière fonction. La seule chose qui change c’est la synthaxe. Mais la sortie est plus élégante,

> reg3 = lm(dist~speed, data = cars)
> delta.method(reg3, "-b1/(b1+b0)")
Functions of parameters:  expression(-b1/(b1+b0))
Estimate = 0.2881585 with se = 0.1063255

En fait, l’écart-type indiqué ici est précisément celui qui est utilisé pour écrire l’intervalle de confiance. Dans cet exemple, l’intervalle de confiance approché (car on suppose que la normalité – en théorie seulement asymptotique – est vérifiée) et estimé (car les paramètres sont inconnus et estimés) est de la forme suivant, pour un seuil de confiance à 95%

> DM=delta.method(reg3, "-b1/(b1+b0)")
Functions of parameters:  expression(-b1/(b1+b0))
Estimate = 0.2881585 with se = 0.1063255
> DM\$estimate+qt(.975,df=nrow(cars)-1)*DM\$se
[1] 0.5018277
> DM\$estimate+qt(.025,df=nrow(cars)-1)*DM\$se
[1] 0.07448936

C’est à dire qu’on a un intervalle de confiance de la forme [0.074;0.502].

# Exposé à Toulouse sur l’estimation (nonparamétrique) de quantiles

Exposé à Toulouse 1, sur l’estimation nonparamétrique de quantiles.

In this talk we propose several nonparametric estimators of quantiles based on Beta kernel and applied to transformed data by the generalized Champernowne distribution initially fitted to the data. A Monte-Carlo based study will show that those estimators improve the efficiency, not only for light tailed distributions, but mainly for heavy tailed, when the probability level is close to 1. Another application will be seen, on portfolio optimization in the mean-VaR context.

# Advances in copula density estimation

Exposé aux Gemeinsame Jahrestagung der Deutschen Mathematiker-Vereinigung und der Gesellschaft für Didaktik der Mathematik, dans la session Multivariate Dependence Modelling using Copulas – Applications in Finance, Humboldt-Universität zu Berlin, Mars 2007, à partir du chapitre écrit avec Jean David Fermanian et Olivier Scaillet sur The estimation of copulas : theory and pratice.

Copulas are a way of formalising dependence structures of random vectors. Although they have been known about for a long time (Sklar (1959)), they have been rediscovered relatively recently in applied sciences (biostatistics, reliability, biology, etc). In finance, they have become a standard tool with broad applications: multiasset pricing (especially complex credit derivatives), credit portfolio modelling, risk management, etc. For example, see Li (1999), Patton (2001) and Longin and Solnik (1995). […] In this chapter, we focus on the practical issues practitioners are faced with, in particular concerning estimation and visualisation. In the first section, we give a general setting for the estimation of copulas. Such a framework embraces most of the available techniques. In the second section, we deal with the estimation of the copula density itself, with a particular focus on estimation near the boundaries of the unit square.