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

Úvodné operácie

1. Nastavenie pracovného adresára :

SetDirectory[NotebookDirectory[]]

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

flow1=ReadList["Rajcianka.txt",Number];
flow1[[1]]

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

flow=Drop[flow1,{1,2}];

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

n=Length[flow]

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,BaseStyle->{10,FontFamily->"Arial"},Joined->True,AxesOrigin->{0,0}];

5.Vykreslenie časového radu priemerných mesačných prietokov a uloženie grafu do premennej:

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

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

Zobrazenie a uloženie  grafu do premennej gd

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

μ=Mean[flow]
σ=StandardDeviation[flow]
σ2=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=Drop[flow,1]-Drop[flow,-1];
Variance[tv1]//N
ListPlot[tv1, PlotRange->{-20,20},Joined->False]

Väčšina hodnôt prvých diferencií je z intervalu (-5; 5), 
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 :

reg=Regress[flow,{t},t,RegressionReport->BestFit]
lintrend=BestFit/.reg
glintrend=Plot[lintrend,{t,0,n},PlotStyle->{Hue[0.05],Thickness[0.005]}]

ltr=Table[lintrend,{t,1,825}];
Save["lintrend.def",lintrend]

c) ZOBRAZENIE  POZOROVANÝCH   DÁT  A  TRENDU

Show[ gdata,glintrend]

d) VYTVORENIE  ČASOVÉHO RADU  REZÍDUÍ  PO  TRENDE,  VÝPOČET   ROZPTYLU,  ZOBRAZENIE REZÍDUÍ

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

reg=Regress[flow,{t},t,RegressionReport->FitResiduals];
rezlintrend=FitResiduals/.reg;
μ=Mean[rezlintrend]//Chop
σ=StandardDeviation[rezlintrend]
σ2=Variance[rezlintrend]

grt=ListPlot[rezlintrend,Joined->True]

Save["rezlintrend.dat",rezlintrend]

2)  KVADRATICKÝ  TREND

Needs["NonlinearRegression`"]

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

NonlinearRegress[flow,{a t^2 + b t + c},{a,b,c},t]

kvreg=Regress[flow,{ t^2,t},t,RegressionReport->BestFit]
kvtrend=BestFit/.kvreg
gkvtrend=Plot[kvtrend,{t,0,n},PlotStyle->{Hue[0.05],Thickness[0.005]}]
Show[ gdata,gkvtrend]

kvreg=Regress[flow,{ t^2,t},t,RegressionReport->FitResiduals];
rezkvtrend=FitResiduals/.kvreg;
μ=Mean[rezkvtrend]//Chop
σ=StandardDeviation[rezkvtrend]
σ2=Variance[rezkvtrend]
grt=ListPlot[rezkvtrend,Joined->True]
Save["rezkvtrend.dat",rezkvtrend]

3)  EXPONENCIÁLNY  TREND

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

expreg=NonlinearRegress[flow,{a b^t},{a,b},t,RegressionReport->BestFit]
exptrend=BestFit/.expreg
gexptrend=Plot[exptrend,{t,0,n},PlotStyle->{Hue[0.05],Thickness[0.005]}]
Show[ gdata,gexptrend]

expreg=NonlinearRegress[flow,{a b^t},{a,b},t,RegressionReport->FitResiduals];
rezexptrend=FitResiduals/.expreg;
μ=Mean[rezexptrend]//Chop
σ=StandardDeviation[rezexptrend]
σ2=Variance[rezexptrend]
grt=ListPlot[rezexptrend,Joined->True]
Save["rezexptrend.dat",rezexptrend]

4) LOGISTICKÝ  TREND

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

logreg=NonlinearRegress[flow,{c /(1 + a b ^t)},{a,b,c},t,RegressionReport->BestFit]
logtrend=BestFit/.logreg
glogtrend=Plot[logtrend,{t,0,n},PlotStyle->{Hue[0.05],Thickness[0.005]}]
Show[ gdata,glogtrend]

logreg=NonlinearRegress[flow,{c /(1 + a b ^t)},{a,b,c},t,RegressionReport->FitResiduals];
rezlogtrend=FitResiduals/.logreg;
μ=Mean[rezlogtrend]//Chop
σ=StandardDeviation[rezlogtrend]
σ2=Variance[rezlogtrend]
grt=ListPlot[rezlogtrend,Joined->True]
Save["rezlogtrend.dat",rezlogtrend]

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

reztrend=   rezlintrend ;

Sezónna zložka

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)

kf=CorrelationFunction[reztrend,120];
myplotcorr1[kf,2/Sqrt[n],Ticks->{{12,24,36,48,60,72,84,96,108},Automatic},AxesLabel->{k,(Subscript[ρ, k])},PlotLabel->"AutoCorrelation Function"]

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

Regress[reztrend,{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:

REG=Regress[reztrend,{Sin[2 Pi t/12],Cos[2 Pi t/12]},t,IncludeConstant->False,RegressionReport->BestFit]
sz=BestFit/.REG

regdata=ListPlot[reztrend,PlotLabel->"",AxesLabel->{"time[month]","flow[(m^3)/s]"},Ticks->{popis,Automatic}]

szplot=Plot[sz,{t,0,n},PlotStyle->{Hue[0.05],Thickness[0.005]}]
Show[regdata,szplot]

szreg=Regress[reztrend,{Sin[2 Pi t/12],Cos[2 Pi t/12]},t,IncludeConstant->False,RegressionReport->FitResiduals];
szrez=FitResiduals/.szreg;

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

szreztrend=ListPlot[szrez,Joined->True]
Save["szreztrend.dat",szreztrend]

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

μ=Mean[szrez]//Chop
σ=StandardDeviation[szrez]
σ2=Variance[szrez]

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

L=12;
M=71;
f[i_,k_,t_]:= If [(i-1+k*L)≤t<(i+k*L),1,0]

x[i_,t_]:= Sum[f[i,k,t],{k,0, M-1}]
funkcie = Table[ x[i,t], {i,2, L}];

Regress[reztrend,funkcie,t,Weights->(#^2&)]

REGK=Regress[reztrend,funkcie,t,RegressionReport->BestFit(*,IncludeConstant->False*)];
k=BestFit/.REGK;
Kplot=Plot[k,{t,0,n},PlotStyle->{Hue[0.05],Thickness[0.005]}]

Show[regdata,Kplot]

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

kreg=Regress[reztrend,funkcie,t,(*IncludeConstant->False,*)RegressionReport->FitResiduals];
krez=FitResiduals/.kreg;
kreztrend=ListPlot[krez,Joined->True]
Save["szmf.def",k];
Save["kreztrend.dat",krez]

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

μ=Mean[krez]//Chop
σ=StandardDeviation[krez]
σ2=Variance[krez]

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]