Analyse statistique de données climatiques

Marie-Cécile Selosse
(last modiﬁcation date: October 10, 2017)

Version pdf de ce document
Version sans bandeaux

1 La Terre se réchauﬀe-t-elle ?
2 À quelle date le changement climatique a-t-il commencé ?
3 Test de nullité d’un coeﬃcient dans un modèle de régression linéaire simple

1 La Terre se réchauﬀe-t-elle ?
2 À quelle date le changement climatique a-t-il commencé ?
2.1 Modèle statistique
2.2 Position du problème
3 Test de nullité d’un coeﬃcient dans un modèle de régression linéaire simple

1 La Terre se réchauﬀe-t-elle ?

On trouve sur le site du Climatic Research Unit les moyennes globales de température de 1856 à 2004 (mois par mois et annuelles). Une copie locale où on n’a conservé que les données de température est disponible localement temperature_globe.txt

Question 1 Saisir le ﬁchier de données de température.

Extraire les moyennes annuelles.

Tracer la courbe tT(t) donnant l’évolution de ces dernières.

  //saisie de données de température
  Tg=fscanfMat('temperatures_globe.txt')
  //on extrait la dernière colonne qui contient les moyennes annuelles
  y=Tg(:,$)
  // t représente le temps (années)
  [m,n]=size(Tg);t=[1:m]';
  // graphique des données
  xbasc();plot2d(t',y');

Question 2 Calculer les coeﬃcients de la droite de régression y = 𝜃₁ + 𝜃₂t qui est la plus proche de la courbe des températures au sens des moindres carrés.

  //on programme une régression linéaire
  //formule de régression a la main
  x=[ones(t),t];
  M=(x'*x)^(-1);
  theta=M*x'*y;
  // On peut remarquer que Scilab effectue la résolution
  // au sens des moindres carrés si x est << tall >>
  // i.e. est une matrice qui a plus de lignes que de colonnes
  // Le problème est bien posé si son rang est égal à son nombre de colonnes,
  // sinon Scilab donne une solution parmi les solutions possibles
  theta=x\y
  // on peut aussi utiliser reglin
  // [a,b,sig]=reglin(t',y')

Question 3 Superposer la courbe des données et la droite de régression.

  //superposition graphique des données et de la droite de régression
  xbasc();
  plot2d(t,y);
  plot2d(t,(x*theta),2,"000");

Question 4 La température T(t) est supposée être une réalisation du modèle T(t) = 𝜃₁ + 𝜃₂t + 𝜀_t, où 𝜀₁, 𝜀₂... est une suite de variables aléatoires indépendantes de même loi normale 𝒩(0,σ²).

Interpréter en terme de paramètres l’hypothèse (H₀) ¡¡ la température est stationnaire ¿¿. Tester cette hypothèse. (On pourra consulter la section 3 a la ﬁn de ce document).

On calcule l’estimateur de σ² appellé s²

  p=2
  Res=y-x*theta;
  s2=Res'*Res/(n-p);

On calcule l’écart type de ^
𝜃 ₂ à partir de la covariance de ^
𝜃 :

sigtheta2=sqrt(s2*M(2,2))

Alors theta(2)/sigtheta2 suit une loi de Student à (n − 2) degrés de liberté :

  T=theta(2)/sigtheta2
  // Utiliser la fonction cdft pour calculer la
  // probabilite qu'une v.a suivant une loi de Student a (n-2)
  // degré de liberté dépasse T

Question 5 Examiner graphiquement la courbe des résidus _t = T(t) − ₁ − ₂t.

Subsiste-t-il une structure dans les résidus ?

Tracer la courbe t( _t,T(t)).

À l’examen de l’histogramme des résidus, l’hypothèse de normalité vous semble-t-elle raisonnable ?

  //examen des résidus
  r=y-x*theta
  plot(r)
  histplot(10,r)

2 À quelle date le changement climatique a-t-il commencé ?

2.1 Modèle statistique

( { Yt = 𝜃1 + 𝜀t 1 ≤ t < τ ( Yt = 𝜃1 + 𝜃2(t − τ ) + 𝜀t τ ≤ t ≤ n

(1)

On suppose les 𝜀_t indépendants et identiquement distribués (i.i.d.) suivant une loi normale 𝒩(0,σ²).

Le vecteur des paramètres est noté

ϕ = (𝜃1,𝜃2,σ2,τ ).

(2)

2.2 Position du problème

On estime les quatre paramètres par la méthode du maximum de vraisemblance.

Déﬁnition de la fonction de vraisemblance

La fonction de vraisemblance de la loi normale est déﬁnie comme suit :

( ) [ { τ−1 n }] ---1-- n -1-- ∑ 2 ∑ 2 L (ϕ,Y ) = √ --- exp − 2σ2 (yt − 𝜃1) + (yt − 𝜃1 − 𝜃2(t − τ)) 2π σ t=1 t=τ

(3)

Pour calculer le maximum de vraisemblance, c’est-à-dire les valeurs de 𝜃₁, 𝜃₂, σ² et τ qui maximisent log L(ϕ,Y ), on choisit de :

ﬁxer d’abord τ ∈{1,...,n} ;
calculer les valeurs ₁^∗(τ), ₂^∗(τ) et ^∗(τ) qui maximisent log L(ϕ,Y ) ;
calculer enﬁn la valeur τ^∗ qui maximise log L^∗(τ) = log L( ₁^∗(τ), ₂^∗(τ),^∗(τ),τ) .

On eﬀectue donc les opérations dans l’ordre suivant :

( ) max log L(ϕ, Y) = max max log L(ϕ, Y) = max log L(^𝜃∗(τ),^𝜃∗(τ),σ^2∗(τ),τ ). τ,𝜃1,𝜃2,σ2 τ=1,...,n 𝜃1,𝜃2,σ2 τ=1,...,n 1 2

Calcul de ₁^∗(τ), ₂^∗(τ) et ^∗(τ)

Question 6 Montrer, par dérivation de (3), que ₁^∗(τ) et ₂^∗(τ) sont solution de

( | ∑ n ∑n ||| ^𝜃∗1(τ)n + ^𝜃∗2(τ) (t − τ) = yt { t= τ t=1 | ∑n ∑ n ∑n ||| ^𝜃∗1(τ) (t − τ ) + ^𝜃∗2(τ) (t − τ)2 = (t − τ)yt ( t=τ t= τ t=τ

(4)

et que

{ } ∗ 1 τ∑−1 ∑n ^σ2 (τ) = -- (yt − ^𝜃∗1(τ ))2 + (yt − ^𝜃∗1(τ) − (t − τ)^𝜃∗2(τ))2 . n t=1 t=τ

(5)

Calcul de τ^∗

Question 7 Programmer en Scilab le calcul de τ^∗ solution de

^∗ ^∗ ^2∗ τm=a1x,...,n log L(𝜃1(τ),𝜃2(τ),σ (τ),τ )

  function [ttheta1,ttheta2,ssigma2]=rupture(tau,y)
    n=size(y,"*");// taille de y
    A=[n,sum(1:n-tau);sum(1:n-tau),sum((1:n-tau)^2)]
    B=[sum(y);(1:n-tau)*y(tau+1:$)'];
    x=A\B;
    ttheta1=x(1,:);
    ttheta2=x(2,:);
    ssigma2=1/n* ...
            (sum((y(1:tau-1)-ttheta1)^2)+sum((y(tau+1:$)-ttheta1-ttheta2*(1:n-tau))^2));
  endfunction