[ Retour aux Travaux Dirigés de SAS ]

/* Séries temporelles avec SAS, Corrigé du TD 2 */
/* Jean-Sebastien Roy (js@jeannot.org), 2000-2001 */

/* Premier jeu de données */
data serie;
  do t=1 to 10*12;
    x=sin(t*(2*3.14)/24)+sin(t*(2*3.14)/(24*5))*3;
    xb=x+sin(t*(2*3.14)/3.3)*5;
    output;
  end;
run;

symbol1 interpol=join color=black value=none;
symbol2 interpol=join color=red value=none;
symbol3 interpol=join color=green value=none;
symbol4 interpol=join color=blue value=none;
symbol5 interpol=join color=orange value=none;
symbol6 interpol=join color=pink value=none;

proc gplot data=serie;
  plot x*t xb*t / overlay;
run;
quit;

/* Moyenne mobile */
/* Avec des lags */
data mlag;
  set serie;
  a3=(xb+lag(xb)+lag2(xb))/3;
  xb=lag(xb);
  x=lag(x);
  t=lag(t);
run;

proc gplot data=mlag;
  plot (x xb a3)*t/ overlay legend;
run;
quit;

/* Avec expand */
proc expand data=serie out=xlag;
  id t;
  convert xb=a3 / transform=(cmovave 3);
  convert xb=a5 / transform=(cmovave 5);
  convert xb=a13 / transform=(cmovave 13);
  convert xb=a23 / transform=(cmovave 23);
run;

proc gplot data=xlag;
  plot (x xb a3 a5 a13 a23)*t / overlay legend;
run;
quit;

/* Moyenne mobile de henderson */
data xlagh;
  set serie;
  h13=(
    -325*xb
    -468*lag1(xb)
    +1100*lag3(xb)
    +2475*lag4(xb)
    +3600*lag5(xb)
    +4032*lag6(xb)
    +3600*lag7(xb)
    +2475*lag8(xb)
    +1100*lag9(xb)
    -468*lag11(xb)
    -325*lag12(xb))/16796;
  t=lag6(t);
  x=lag6(x);
  xb=lag6(xb);
  if h13 ne .;
run;

proc gplot data=xlagh;
  plot (x xb h13)*t / overlay legend;
run;
quit;

/* Moyenne mobile de Henderson avec expand */
proc expand data=serie out=xlag2;
  id t;
  convert xb=a5_5_5 / transform=(cmovave 5 cmovave 5 cmovave 5);
  convert xb=h13_p / transform=(cmovave (0 0 0 1100 2475 3600 4032 3600 2475 1100 0 0 0));
  convert xb=h13_n / transform=(cmovave (325 468 0 0 0 0 0 0 0 0 0 468 325));
run;

data xlag2;
  set xlag2;
  h13=(h13_p*(+1100 +2475 +3600 +4032 +3600 +2475 +1100)
       -h13_n*(+325 +468 +468 +325))/16796;
run;

proc gplot data=xlag2;
  plot (x xb a5_5_5 h13)*t      / overlay legend;
run;
quit;

/* Annulation d'une composante */
data serie2;
  do t=1 to 10*12;
    x=sin(t*(2*3.14)/(17));
    xb=x+sin(t*(2*3.14)/7)*2;
    output;
  end;
run;

proc expand data=serie2 out=xlag3;
  id t;
  convert xb=a3 / transform=(cmovave 3);
  convert xb=a5 / transform=(cmovave 5);
  convert xb=a7 / transform=(cmovave 7);
  convert xb=a13 / transform=(cmovave 13);
run;

proc gplot data=xlag3;
  plot (x xb a3 a5 a7 a13)*t / overlay legend;
run;
quit;

/* Analyse d'une moyenne mobile */

%let width=30;
data one;
  x=0;
  do t=-&width to -1;
    output;
  end;
  t=0;x=1;output;
  x=0;
  do t=1 to &width;
    output;
  end;
run;

proc expand data=one out=weights;
  id t;
  convert x=a5_7_13 / transform=(cmovave 5 cmovave 7 cmovave 13);
  convert x=s / transform=(cmovave (1 2 2 2 1));
  convert x=a17 / transform=(cmovave 17);
  convert x=h13_p / transform=(cmovave (0 0 0 1100 2475 3600 4032 3600 2475 1100 0 0 0));
  convert x=h13_n / transform=(cmovave (325 468 0 0 0 0 0 0 0 0 0 468 325));
run;

data weights;
  set weights;
  h13=(h13_p*(+1100 +2475 +3600 +4032 +3600 +2475 +1100)
       -h13_n*(+325 +468 +468 +325))/
      (-325 -468 +1100 +2475 +3600 +4032 +3600 +2475 +1100 -468 -325);
  drop h13_p h13_n;
run;

proc gplot data=weights;
  plot (a5_7_13 s a17 h13)*t / overlay legend;
run;
quit;

/* Reduction d'un bruit blanc */

data bb;
  do t=1 to 100;
    x=rannor(0);output;
  end;
run;

proc expand data=bb out=xlag4;
  id t;
  convert x=a3 / transform=(cmovave 3);
  convert x=a44 / transform=(cmovave 4 cmovave 4);
  convert x=a5 / transform=(cmovave 5);
  convert x=a7 / transform=(cmovave 7);
  convert x=a17 / transform=(cmovave 17);
  convert x=m5 / transform=(cmovmed 5);
  convert x=s / transform=(cmovave (1 2 2 2 1));
run;

proc gplot data=xlag4;
  plot (x a3 a44 a5 a7 m5)*t / overlay vref=0 legend;
  plot (x s a17)*t / overlay vref=0 legend;
run;
quit;

proc summary data=xlag4 print var;
  var x a3 a44 a5 a7 m5 s;
run;

/* Tau (Slutsky-Yule) */

proc transpose data=weights out=tweights
  (rename=(_name_=nom col1=w));
  by t;
run;

proc sort data=tweights;
  by nom t;
run;

data tau;
  set tweights;
  ttm1 = w*lag(w);
  t2 = w*w;
run;

proc sql;
  select nom,2*3.1415/arcos(sum(ttm1)/sum(t2))
  from tau
  group by nom;
quit;