Určenie systematických zložiek časového radu

Úvodné operácie

SetDirectory["C:\\Documents and Settings\\stdn03.B206\\My Documents\\jozo urban\\ACR\\4.3.2009"]
flow1=ReadList["Bebrava.txt",Number];
flow=Drop[flow1,{1,2}];
n=Length[flow]

<<"timeseri\\datasmoo.m"
<<"timeseri\\timeseri.m"
<<"timeseri\\userfunc.m"

<<"LinearRegression`"
<<"NonlinearRegression`"

SetOptions[ListPlot,PlotRange->All,BaseStyle->{10,FontFamily->"Arial"},Joined->True,AxesOrigin->{0,0}];
popis=Table[{60+i*120,1936+i*10},{i,0,6}];
gdata=ListPlot[flow,PlotLabel->"Month Averages, Bebrava",AxesLabel->{"time[month]","flow[(m^3)/s]"},Ticks->{popis,Automatic}];

Trend

Pre modelovanie časového radu priemerných mesačných prietokov sme vybrali  lineárny trend.

reztrend=ReadList["rezlintrend.dat"][[1]];
trend=ReadList["lintrend.def"][[1]]

Sezónna zložka

rezsz=ReadList["kreztrend.dat"][[1]];
sz=ReadList["szmf.def"][[1]];

ZOBRAZENIE  POZOROVANÝCH   DÁT  A  FUNKCIE, KTORÁ JE SÚČTOM TRENDU A SEZÓNNEJ ZLOŽKY

func=Plot[k + lintrend/.t->i,{i,1,852},PlotStyle->RGBColor[1,0,0]];
Show[gdata,func]

Cyklická zložka

1. Odhad výberového spektra a vykreslenie periodogramu:

spec=Spectrum[rezsz];

Length[spec]

gs=plotspectrum[spec,AxesLabel->{"ω","f(ω)"},Joined->True]

2. Vypísanie 20 najväčších hodnôt spektra (spolu s frekvenciami):

dvojice=Transpose[{spec,Table[j,{j,0,n-1}]}];
Take[Sort[dvojice],-20]

3. Určenie vyznamných frekvencií  pomocou Fisherovho testu :

m=Round[n/2];
α=0.05;

Normované výberové spektrum:

mF=m;
Y[j_]:=dvojice[[j,1]]/Sum[dvojice[[i,1]],{i,m}]

Usporiadať a vybrať poslednú dvojicu:

Sort[Table[{Y[j],j-1},{j,m}]];
sY1=Take[Sort[Table[{Y[j],j-1},{j,m}]],{-1}]

Odhad kritickej hodnoty:

1-(α/mF)^(1/(mF-1))

mF2=m-1;
sY12=Take[Sort[Table[{Y[j],j-1},{j,m}]],{-2}]
1-(α/mF2)^(1/(mF2-1))

mF3=m-2;
sY13=Take[Sort[Table[{Y[j],j-1},{j,m}]],{-3}]
1-(α/mF3)^(1/(mF3-1))

Rozhodnúť, či je daná frekvencia na základe Fisherovho testu významná alebo nie.
Opakovať, až kým nezistíme, že nasledujúca frekvencia nie je významná.

4.  Priradiť významné frekvencie do premennej vf a pomocou regresnej analýzy určiť funkciu  modelujúcu cyklickú zložku, výslednú funkciu priradiť  do premennej cz :

Vypočítať  časový ekvivalent významných frekvencií , t.j. dlzku periody:

n/6/12.

n/14.

Časový ekvivalent významnej frekvencie je 6 mesiacov.

Významné funkcie:

vf={2π*6/n,2π*14/n};

Zobrazenie periodogramu s vyznačenými významnými funkciami:

gs=plotspectrum[spec,AxesLabel->{"ω","f(ω)"},Ticks->{vf,Automatic},Joined->True]

Určenie cyklickej zložky regresiou:

funkcie=Flatten[Table[{Sin[vf[[i]]*t],Cos[vf[[i]]*t]},{i,Length[vf]}]];

regcz=Regress[rezsz,funkcie,t,IncludeConstant->False]

regcz=Regress[rezsz,{Sin[((π t)/71)],Cos[((7 π t)/213)]},t,RegressionReport->{  BestFit,FitResiduals}, IncludeConstant->False];

cz=BestFit/.regcz

5. Určiť rezíduá,  priradiť do premennej rez, vykresliť a vypočítať rozptyl rezíduí:

rez=FitResiduals/.regcz;

6. Všetky tri systematické zložky priradiť do premennej syz a zobraziť  spolu s pôvodným časovým radom  pre celé obdobie:

syz=trend+sz+cz;

7. Vykresliť priebeh systematických zložiek počas desiatich rokov:

plsyz=Plot[syz,{t,1,n},PlotStyle->RGBColor[1,0,0]];

Show[gdata,plsyz]

8. Uložiť časový rad reziduí po systematických zložkách do dátového súboru:

(*DeleteFile["rezBebrava.dat"]*)
Save["rezBebrava.dat",rez]

Funkciu pre systematické zložky uložiť do definičného súboru:

(*DeleteFile["syzBebrava.def"]*)
Save["syzBebrava.def",syz]

Testovanie nulovosti autokorelačnej funkcie

corrx=CorrelationFunction[rez,100];

myplotcorr1[corrx,2/Sqrt[n],AxesLabel->{"k","(Subscript[ρ, k])"},PlotLabel->"Autocorrelation Function"]

Z grafu je zrejmé, že reziduá sú ešte korelované, teda má zmysel ich ďalej modelovať napr. lineárnymi ARMA modelmi.