> restart:

TP 07

Série temporelle

Rappel de consignes

(1) titre en "Times New Roman", corps 36,

NOM Prénom - groupe

TP 07 - date

(2) sauvegarder "souvent"

(3) imprimer : paginer (Format/Page_Number)

puis File/Preview.

enfin imprimer en deux colonnes

>

On charge quelques bibliothèques

Vérifier les numéros de version des bibliothèmes disponibles

> with(pldx); # il faut la version 24
ever_loaded;

> with(simul); # il faut la version 16
ever_loaded;

Warning, the name ever_loaded has been redefined

> xima:=proc(pl); op(1,pl) : convert(%,list): map(op,%):
map2(op,2,%): max(op(%)): end:

>

Lecture du fichier

Télécharger le fichier "www.douillet.info/~douillet/mao_tp/mao_dat_07.txt"
ABS, beer.dat, from jan 1956

Déclarer le nom de la copie locale du fichier

> fich:= "F:/docs/Ensait/maotp/maotp07/dat_maotp07.txt";

On lit tous les nombres

> lire:= table(): tmp:= 'tmp':
hdl:= fopen(fich, READ) :
for i to 5000 while tmp <> [ ] do
tmp:= fscanf(hdl,"%a");
lire[i]:= op (tmp);
od:
fclose(fich); n:= i-2 ; i:='i':

On en fait une liste

> li:= [seq(lire[k],k=1..n)]:

>

Paramètres de dispersion

> macro(moy= stats[describe, mean] ); m:= moy(li);

> macro(var= stats[describe, variance] ); s, s2:= sqrt(var(li)), var(li);

> macro(range= stats[describe, range] ); vals:= range(li);
macro(count= stats[describe, count] ); count(li);

>

Visualisation

> an0:= 1956; macro(ptplo=plots[pointplot]);

> ptplo([seq([an0+k/12, li[k]], k=1..n)], style= line, view=[DEFAULT, 0..rhs(vals)], color=magenta); pl_data:= %:

> pl2:= plot([m-s, m, m+s], t=an0..an0+n/12, color=red, linestyle=[8,24,8]) :

> display(pl_data, pl2);

>

>

Auto-corrélation (E2)

> nn:=20*12:

> for j from 0 to nn-2 do
sx,sy:= add(li[k], k=1..nn-j)/(nn-j),add(li[k+j], k=1..nn-j)/(nn-j);
sxx,syy:= add(li[k]^2, k=1..nn-j)/(nn-j),add(li[k+j]^2, k=1..nn-j)/(nn-j);
sxy:= add(li[k]*li[k+j], k=1..nn-j)/(nn-j);
rr[j]:= (sxy-sx*sy)^2/(sxx-sx^2)/(syy-sy^2);
od:

> ptplo([seq([k,rr[k]], k=0..nn-2)], style=line, labels=[decalage, rr]);

>

Transformée de Fourier

> nn:=12*2: ft:= table([]);

> ft:= [seq( add(li[k]*exp(2*I*evalf(Pi)*k*j/nn), k=1..nn)/nn, j=0..nn)];

> ptplo( [seq([k, abs(ft[k+1])], k=0..nn)], style=line, view=[DEFAULT, 0..20]);

>

Attention : cela dure "un certain temps"

> nn:=n : ft:= table([]);

> stamp:= time() :
ft:= evalf([seq( add(li[k]*exp(2*I*Pi*k*j/nn), k=1..nn)/nn, j=0..nn)]) :
elapsed_time = time() - stamp ;

> ptplo( [seq([k, abs(ft[k+1])], k=0..nn)], style=line, view=[DEFAULT, 0..20]);

> seq([k, abs(ft[k+1])], k=37..42); pic:= add(k * abs(ft[k+1]), k=37..42)/add(abs(ft[k+1]), k=37..42);

> # pic := 39.63866159;

> n/pic; per:= round(%); lon:= floor(n/per);

> perg := floor(per/2); perd := per-%;

>

Moyenne mobile (trend)

> moymob:= proc(ma::list, per) global perg, perd; local k, kx, mb, tmp;
kx := nops(ma):
for k from 1 to perg-1 do mb[k] := missing od:;
for k from kx-perd+1 to kx do mb[k] := missing od:
tmp := add(ma[k], k=1..per); mb[perg] := tmp/per;
for k from perg+1 to kx-perd do
tmp := tmp-ma[k-perg]+ma[k+perd]: mb[k] := %/per
od: [seq(mb[k],k=1..kx)];
end:

> trend:= moymob(li, per):

> pl_trend:= plot([seq([an0 + k/12, trend[k]],k=1..n)], color=blue):

> display(pl_data, pl_trend);

>

>

Correction des Variations saisonnières

Méthodes possibles (E2) :

correction simple

correction logarithmique

correction trend simple

correction trend logarithmique

> for u from 1 to per do meu[u]:= moy([seq(li[u+per*j],j=0..lon-1)]); od:
meu_:= moy([seq(meu[u],u=1..per)]);
for u from 1 to per do meu[u]:= meu[u]-meu_: od:

> licor:= [seq(li[k] -meu[1+irem(k-1,per)], k=1..n)] :

> [count, moy, var](licor); [count, moy, var](li);

>

> pl_sai_sim:= plot([seq( [k/12+1956, licor[k] ], k=1..n)], linestyle=16):

> display( pl_data, pl_sai_sim);

>

Résiduels

> resid := [missing$(perg-1), seq( licor[k] - trend[k], k=perg .. n-perd), missing$perd] : nops(%), n;

> plot([seq( [k/12+1956, resid[k]], k=1..n)], linestyle=8) ;

Transformée de Fourier des résiduels

Attention : cela dure "un certain temps"

> size:= 12*(lon-1) ; nn:=perg-1+size ; ftr := table([]) ;

> stamp:= time() :
ftr := evalf([seq( add(resid[k]*exp(2*I*Pi*k*j / size), k=perg..nn)/ size, j=0..size)]) :
elapsed_time = time() - stamp ;

> ptplo( [seq([k, abs(ftr[k+1])], k=0..size)], style=line);

> ran:= 157..161; seq([k, abs(ftr[k+1])], k=ran);
pic2:= add(k * abs(ftr[k+1]), k=ran)/add(abs(ftr[k+1]), k=ran); LL:= size/%;

>

> ran:= 37..40; seq([k, abs(ftr[k+1])], k=ran);
pic3:= add(k * abs(ftr[k+1]), k=ran)/add(abs(ftr[k+1]), k=ran); size/%;

>

Modèle y = polynome + trigo (frv=8.57)

> phi:= 2*Pi/12 ;
add(a[j]*(t/12)^j, j=0..4)+ add(cos(phi*t*j)*cc[j]+sin(phi*t*j)*ss[j], j=1..4) - y :
fu:= unapply(%, t,y);
vars:= op( indets(fu(t,y), indexed) );

> Digits:=20:
chi2_f := collect(add( (evalf@fu)(k, licor[k])^2, k=1..n), [vars], distributed):
soo := solve(map2(diff,chi2_f, {vars})) : fuu := unapply(subs(%, fu(t,y)+y),t):
Digits :=10 : evalf(fuu(t));

> resid1:= [seq([k/12+1956, fuu(k)-licor[k]], k=1..n)] :

> subs(missing=NULL, moymob(resid1,12)): display(plot(%), plot(resid1, color=magenta));

> Digits:=10:
evalf(map2(op,2,resid1)): [count, moy, var](%); frv_cor:= var(licor)/%[3];

>

Modèle ln(y) = polynome + trigo (frv=8.57)

> phi:= 2*Pi/12 ;
add(a[j]*(t/12)^j, j=0..4)+ add(cos(phi*t*j)*cc[j]+sin(phi*t*j)*ss[j], j=1..4) - ln(y) :
gu:= unapply(%, t, y);
vars:= op( indets(gu(t,y), indexed) );

> Digits:=20:
chi2_g := collect(add( (evalf@gu)(k, licor[k])^2, k=1..n), [vars], distributed) :
sog := solve(map2(diff,chi2_g, {vars})) : guu := unapply(subs(%, gu(t,y)+ln(y)),t):
Digits :=10 : evalf(guu(t));

> resid_g:= [seq([k/12+1956, exp(guu(k)) - licor[k]], k=1..n)] :

> subs(missing=NULL, moymob(resid_g,12)): display(plot(%), plot(resid_g, color=magenta));

> Digits:=10:
evalf(map2(op,2,resid_g)): [count, moy, var](%); frv_cor:= var(licor)/%[3];

> resid_g_ln := [seq( evalf(guu(k)) - ln(licor[k]), k=1..n)] :

> frv_cor_ln := var(map(ln, licor))/var(resid_g_ln);

>

>

>