Určenie systematických zložiek

Ú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]

Trendy

1)  KONŠTANTNÝ  A  LINEÁRNY  TREND

a) Test

Vyčísliť a zobraziť prvé diferencie časového radu
(V prípade, že prvé diferencie sú približne nulové, vhodný je konštantný trend;
  ak sú prvé diferencie približne konštantné, ale nenulové, vhodný je lineárny trend.)

tv1=flow[[2;;n]]-flow[[1;;n-1]];
Variance[tv1]//N
ListPlot[tv1,Joined->False,PlotRange->{-5,5}]

Väčšina hodnôt prvých diferencií je z intervalu (-2; 2), takže ich môžeme považovať za približne konštantné, ale nenulové.Vzhľadom na to by mal byť vhodným trendom lineárny trend.

b) Odhad a zobrazenie funkcie trendu

Odhadneme parametre lineárneho trendu lineárnou regresiou.Najprv použijeme štandardný tvar príkazu Regress bez volieb :

Regress[flow,{t},t]

Hodnoty všetkých PValue z časti Parameter Table sú veľmi malé, preto môžeme tento model trendu považovať za vyhovujúci.

Príkaz Regress s voľbou, ktorá umožní vypísať parametre regresnej funkcie :

reglin=Regress[flow,{t},t,RegressionReport->BestFit]

lintrend=BestFit/.reglin

glintrend=Plot[lintrend,{t,0,n},PlotStyle->{Hue[0.05],Thickness[0.005]}]

c) Zobrazenie pozorovanych dat a trendu

Show[{gdata,glintrend}]

d) Vytvorenie casoveho radu rezidui po trende, vypocet rozptylu, zobrazenie rezidui

Použijeme príkaz Regress s voľbou, ktorá vytvorí časový rad rezíduí po odčítaní hodnôt lineárneho trendu :

reglin2=Regress[flow,{t},t,RegressionReport->FitResiduals];
rezlintrend=FitResiduals/.reglin2;
milin=Mean[rezlintrend]//Chop
sigmalin=StandardDeviation[rezlintrend]
sigma2lin=Variance[rezlintrend]

grtlin=ListPlot[rezlintrend]

Save["rezlintrend.dat",rezlintrend]

2)  KVADRATICKÝ  TREND

a) Test

tv2=flow[[3;;n]]+2*flow[[2;;n-1]]+flow[[1;;n-2]];
Variance[tv2]//N
ListPlot[tv2,Joined->False,PlotRange->{-20,20}]

b) Odhad a zobrazenie funkcie trendu

Regress[flow,{t,t^2},t]

Hodnoty všetkých PValue z časti Parameter Table sú veľmi malé, preto môžeme tento model trendu považovať za vyhovujúci.

Príkaz Regress s voľbou, ktorá umožní vypísať parametre regresnej funkcie :

regkvad=Regress[flow,{t,t^2},t,RegressionReport->BestFit]

kvadtrend=BestFit/.regkvad

gkvadtrend=Plot[kvadtrend,{t,0,n},PlotStyle->{Hue[0.05],Thickness[0.005]}]

c) Zobrazenie pozorovanych dat a trendu

Show[{gdata,gkvadtrend}]

d) Vytvorenie casoveho radu rezidui po trende, vypocet rozptylu, zobrazenie rezidui

regkvad2=Regress[flow,{t,t^2},t,RegressionReport->FitResiduals];
rezkvadtrend=FitResiduals/.regkvad2;
mikvad=Mean[rezkvadtrend]//Chop
sigmakvad=StandardDeviation[rezkvadtrend]
sigma2kvad=Variance[rezkvadtrend]

grtkvad=ListPlot[rezkvadtrend]

Save["rezkvadtrend.dat",rezkvadtrend]
Save["kvadtrend.dat",kvadtrend]

3)  EXPONENCIÁLNY  TREND

a) Test

tv3=flow[[2;;n]]/flow[[1;;n-1]];
Variance[tv3]//N
ListPlot[tv3,Joined->False,PlotRange->{-5,5}]

b) Odhad a zobrazenie funkcie trendu

NonlinearRegress[flow,a*b^t,{a,b},t]

Príkaz Regress s voľbou, ktorá umožní vypísať parametre regresnej funkcie :

regexp=NonlinearRegress[flow,a*b^t,{a,b},t,RegressionReport->BestFit]

exptrend=BestFit/.regexp

gexptrend=Plot[exptrend,{t,0,n},PlotStyle->{Hue[0.05],Thickness[0.005]}]

c) Zobrazenie pozorovanych dat a trendu

Show[{gdata,gexptrend}]

d) Vytvorenie casoveho radu rezidui po trende, vypocet rozptylu, zobrazenie rezidui

regexp2=NonlinearRegress[flow,a*b^t,{a,b},t,RegressionReport->FitResiduals];
rezexptrend=FitResiduals/.regexp2;
miexp=Mean[rezexptrend]//Chop
sigmaexp=StandardDeviation[rezexptrend]
sigma2exp=Variance[rezexptrend]

grtexp=ListPlot[rezexptrend]

Save["rezexptrend.dat",rezexptrend]

4) LOGISTICKÝ TREND

a) Test

tv4=((1./flow[[3;;n]])-(1./flow[[2;;n-1]]))/((1./flow[[2;n-1]])-(1./flow[[1;;n-2]]));
Variance[tv4]//N
ListPlot[tv4,Joined->False,PlotRange->{-5,5}]

b) Odhad a zobrazenie funkcie trendu

NonlinearRegress[flow,c/(1+a*b^t),{a,b,c},t]

Príkaz Regress s voľbou, ktorá umožní vypísať parametre regresnej funkcie :

reglog=NonlinearRegress[flow,c/(1+a*b^t),{a,b,c},t,RegressionReport->BestFit]

logtrend=BestFit/.reglog

glogtrend=Plot[logtrend,{t,0,n},PlotStyle->{Hue[0.05],Thickness[0.005]}]

c) Zobrazenie pozorovanych dat a trendu

Show[{gdata,glogtrend}]

d) Vytvorenie casoveho radu rezidui po trende, vypocet rozptylu, zobrazenie rezidui

reglog2=NonlinearRegress[flow,c/(1+a*b^t),{a,b,c},t,RegressionReport->FitResiduals];
rezlogtrend=FitResiduals/.reglog2;
milog=Mean[rezlogtrend]//Chop
sigmalog=StandardDeviation[rezlogtrend]
sigma2log=Variance[rezlogtrend]

grtlog=ListPlot[rezlogtrend]

Save["rezlogtrend.dat",rezexptrend]

Pre modelovanie časového radu priemerných mesačných prietokov sme vybrali LINEÁRNY trend

Sezonna zlozka

1. Korelogram časového radu

Autokorelačnú funkciu vypočítame a vykreslíme pre 10 rokov (120 údajov).Na popis osi s posunutím k použijeme predpokladané maximá (v násobkoch 12 mesiacov)

cf=CorrelationFunction[rezlintrend,120];

myplotcorr1[cf,2/Sqrt[n],Ticks->{{12,24,36,48,60,72,84,96,108},Automatic},AxesLabel->{k,Subscript[p,k]},PlotLabel->Style["Autocorrelation function",Large,Bold]]

2. Určenie sezónnej zložky pomocou matematických funkcií

Funkciu, ktorá bude modelovať sezónnu zložku určíme regresnou analýzou.
 Modelovať budeme reziduá po lineárnom trende. 
Pretože ich stredná hodnota je nulová, v príkaze Regress použijeme voľbu   IncludeConstant -> False

rezlintrend;
Regress[rezlintrend,{Sin[2*Pi*t/12],Cos[2*Pi*t/12]},t,IncludeConstant->False]

Vypočítame regresnú funkciu rezíduí (ktorú priradíme do premennej sz) a rad rezíduí po sezónnej zložke:

regsez=Regress[rezlintrend,{Sin[2*Pi*t/12],Cos[2*Pi*t/12]},t,IncludeConstant->False,RegressionReport->BestFit];
sz=BestFit/.regsez;
seztrend=Plot[sz,{t,0,12},PlotStyle->{Hue[0.05],Thickness[0.005]}]

Zobrazenie  rezíduí po trende a sezónnej zložke:

regsez2=Regress[rezlintrend,{Sin[2*Pi*t/12],Cos[2*Pi*t/12]},t,IncludeConstant->False,RegressionReport->FitResiduals];
rezseztrend=FitResiduals/.regsez2;
grtsez=ListPlot[rezseztrend]

Odhadneme strednú hodnotu, smerodajnú odchýlku a rozptyl reziduí:

milog=Mean[rezseztrend]//Chop
sigmalog=StandardDeviation[rezseztrend]
sigma2log=Variance[rezseztrend]

3. Určenie sezónnej zložky pomocou kvalitatívnych premenných

Skúsime modelovať sezónnu zložku aj pomocou kvalitatívnych premenných.Tento model budeme ďalej používať len v prípade,ak bude index determinácie(R^2)   väčší (a odhadovaný rozptyl reziduí menší) ako v modeli s goniometrickými funkciami

f[i_,k_,t_]:=If[(i - 1+k*12)<= t<(i+k*12),1,0]
x[i_,t_]:=Sum[f[i,k,t],{k,0,71-1}]
funkcie=Table[x[i,t],{i,2,12}];
Regress[rezlintrend,funkcie,t]
regkp=Regress[rezlintrend,funkcie,t,RegressionReport->BestFit];
regkp2=Regress[rezlintrend,funkcie,t,RegressionReport->FitResiduals];
szp=BestFit/.regkp;
rezszp=FitResiduals/.regkp2;

Zobrazenie  rezíduí po trende a sezónnej zložke:

grtsezp=ListPlot[rezszp]

Save["rezsz.dat",rezszp];
Save["sz.dat",szp];

Odhadneme strednú hodnotu, smerodajnú odchýlku a rozptyl reziduí:

Mean[rezszp]
StandardDeviation[rezszp]
Variance[rezszp]

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

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