Úvodné operácie

1. Nastavenie pracovného adresára :

Textovy subor sa ma nachadzat v rovnakom priecinku ako subor .nb

SetDirectory[NotebookDirectory[]]

2. Načítanie dát zo súboru rieka.txt
(hodnoty priemerných mesačných prietokov v(m^3)/sza obdobie od novembra 1930 do decembrom 2001)

flow1=ReadList["Boca.txt",Expression];
(*Odstranenie prebytocnych zatvoriek*)
flow=Flatten[flow1];

Odstránenie prvých dvoch mesiacov (november, december 1930) :

flow=flow[[3;;Length[flow]]];

Vyčíslenie dĺžky časového radu priemerných mesačných prietokov :

n=Length[flow];
(*Vycislenie poctu rokov*)
nyear=n/12;

3. Nahratie balíkov pre analýzu časových radov :

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

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

4. Nastavenie podmienok pre vykresľovanie údajov príkazom ListPlot :

SetOptions[ListPlot,PlotRange->All,Joined->True,AxesOrigin->{0,0}];

5. Vykreslenie časového radu priemerných mesačných prietokov :

Škálovanie časovej osi (s krokom 10 rokov)

popis=Table[{60+i*120,1936+i*10},{i,0,6}];

Zobrazenie

gdata=ListPlot[flow,PlotLabel->Style["Month Averages, Boca",Large,Bold],AxesLabel->{Style["time [month]"],Style["flow [(m^3)/s]"]},Ticks->{popis,Automatic}]

6. Výpočet výberovej strednej hodnoty, smerodajnej odchýlky a rozptylu časového radu priemerných mesačných prietokov :

mi= Mean[flow]
sigma=StandardDeviation[flow]
sigma2=Variance[flow]

Trend

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

trend=ReadList["kvadtrend.dat"];
reztrend=ReadList["rezkvadtrend.dat"];

Sezónna zložka

Pre modelovanie časového radu priemerných mesačných prietokov sme vybrali sezonnu zlozku podla KVALITATIVNYCH PREMENNYCH.

sz=ReadList["sz.dat"];
rezsz=ReadList["rezsz.dat"];

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

sucet=trend+sz;
gsucet=Plot[sucet,{t,0,n},PlotStyle->{Hue[0.05],Thickness[0.005]}];
Show[{gdata,gsucet},PlotRange->{{0,120},{0,10}}]

Cyklická zložka

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

spec=Spectrum[rezsz[[1]]];

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];
alfa=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-(alfa/mF)^(1/(mF-1))

Druha hodnota:

mF=m-1;
sY2=Take[Sort[Table[{Y[j],j-1},{j,m}]],{-2}]
1-(alfa/mF)^(1/(mF-1))

Tretia hodnota:

mF=m-2;
sY3=Take[Sort[Table[{Y[j],j-1},{j,m}]],{-3}]
1-(alfa/mF)^(1/(mF-1))

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/59//N
n/30//N

Významné funkcie:

vf={2*Pi*59/n,2*Pi*30/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]}]];

Regress[rezsz[[1]],funkcie,t]

funkcie=funkcie[[{2,4}]];

regcz=Regress[rezsz[[1]],funkcie,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;

gsyz=Plot[syz,{t,0,n},PlotStyle->{Hue[0.05],Thickness[0.005]}];
Show[{gdata,gsyz},PlotRange->{{120,240},{0,10}}]