Lineárne ARMA modely  časového radu rezíduí  po systematických zložkách priemerných mesačných prietokov

1. Nastavenie pracovného adresára :

SetDirectory[NotebookDirectory[]]
SetOptions[ListPlot,Joined->True,PlotRange->All,
AxesLabel->{"time [Month]","flow [(m^3)/s]"},AxesOrigin->{0,0},BaseStyle->{11,FontFamily->"Arial"}];
flow1=ReadList["Rajcianka.txt",Number];
flow=Drop[flow1,2];
n1=Length[flow]
T1=n1/12

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

<<"timeseri\\datasmoo.m"
<<"timeseri\\timeseri.m"
<<"timeseri\\userfunc.m"
<<"LinearRegression`"
<<"NonlinearRegression`"

1. Načítanie nameraných a  rezíduí po systematických zložkách,   
vytvorenie skúšobnej a testovacej vzorky

Načítanie   časového radu priemerných mesačných prietokov, 
vytvorenie  testovacej vzorky (70 rokov bez posledných 12 mesiacov) a skúšobnej vzorky ( posledných 12 mesiacov) :

tv=Drop[flow,-12];
sv=Take[flow,-12];

n=Length[tv]
Length[sv]

Vykreslenie a uloženie grafu skúšobnej a testovacej vzorky priemerných mesačných prietokov do premenných:

gd=ListPlot[tv,Joined->True]

gsv=ListPlot[sv, PlotRange->{0,25},Joined->True]

Načítanie funkcie systematických zložiek:

syz=ReadList["syzRajcianka.def"][[1]];

Načítanie  časového  radu rezíduí po systematických zložkách,
 vytvorenie testovacej vzorky (70 rokov bez posledných 12 mesiacov) a skúšobnej vzorky ( posledných 12 mesiacov) :

rez=ReadList["rezRajcianka.dat"];

rezsv=Take[flow,-12];
reztv=Drop[flow,-12];

n=Length[reztv]

Vykresliť a uložiť graf skúšobnej a testovacej vzorky rezíduí do premenných:

greztv=ListPlot[reztv,Joined->True]
Save["greztv.dat",greztv]

grezsv=ListPlot[rezsv, PlotRange->{0,25},Joined->True]
Save["grezsv.dat",grezsv]

Vyčísliť strednú hodnotu, rozptyl a smerodajnú odchýlku skúšobnej vzorky  :

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

Priradenie časového radu rezíduí do premennej y:

y=reztv;

Box - Jenkinsova metodológia

2. Identifikácia modelu na základe autokorelačnej a parciálnej autokorelačnej funkcie

(*myplotcorr1[*)
kf=CorrelationFunction[y,20];
myplotcorr1[kf,2/Sqrt[n],PlotLabel->"AutoCorrelation Function"]

Korelačná funkcia   je štatisticky rovná nule  pre k ≥ 12.

kf=CorrelationFunction[y,11];
myplotcorr1[D[kf],2/Sqrt[n],PlotLabel->"Partial AutoCorrelation Function"]

Parciálna autokorelačná funkcia je štatisticky rovná nule    pre k ≥ 11.

3. Identifikácia modelu pomocou rozptylu jednotlivých AR modelov

LEVINSON-DURBINOVA   PROCEDÚRA

Procedúra, ktorá vypočíta parametre a disperziu  pre daný počet AR modelov časového radu rezíduí

armodels = LevinsonDurbinEstimate[y, 10]

Rozptyl vypočítaných modelov :

dis=Map[#[[-1]]&, armodels]

ListPlot[dis, Joined -> True,AxesLabel -> {"p", "sigma2"}]

Rozptyl sa  stabilizoval pre hodnotu p=5.

4. Odhad parametrov modelu

HANNAN - RISSANENOVA PROCEDURA

Použijeme Hannan Rissanenovu procedúru, s hodnotami parametrov (Subscript[p, max ])=10,  (Subscript[q, max ])= 8 a  vypíšeme päť najlepších modelov :

hrx=HannanRissanenEstimate[y,25,10,8,5]

Pre ľahšiu orientáciu zostavíme tabuľku, ktorá obsahuje nasledujúce položky:
 typ modelu,  AIC a BIC modelu, rozptyl reziduí a či je stacionárny

Do[Print["k=",k,"  AIC=",AIC[hrx[[k]],n],"       BIC=" ,BIC[hrx[[k]],n],"  ",hrx[[k]][[-1]],"      ",StationaryQ[hrx[[k]]]],{k,1,5}]

Najlepší model vzhľadom na hodnotu BIC je AR(1).

5. Vybrať model:  AR(1)

Spresnenie odhadu  parametrov vybraného modelu metódou podmienenej maximálnej vierohodnosti :

ar=ConditionalMLEstimate[y,hrx[[1]]]

Výpočet štandardných chýb odhadov :

Koeficient modelu priradíme do premennej φ1 
a pomocou kovariančnej matice odhadneme jeho štandardnú chybu:

φ1=ar[[1,1]];
M=InformationMatrix[y,ar]//MatrixForm
km=Inverse[M];
Ss=Sqrt[Table[km[[i,i]],{i,Length[km]}]/n];
Column[Transpose[{Flatten[{ar[[1]]}],Ss}]]

Určenie rezíduí modelu, odhad ich strednej hodnoty, rozptylu,  smerodajnej odchýlky:

rez1=Residual[y,ar];
Mean[rez1]
Variance[rez1]
StandardDeviation[rez1]

Zobrazenie grafu rezíduí vybraného modelu:

Diagnostická kontrola modelu rezíduí    -   toto dnes nie - bude sa robit ku kazdemu modelu ...!!! 18.3.09

a) Testovanie  nulovosti autokorelačnej funkcie reziduí modelu:

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

Hodnoty autokorelačnej funkcie nie sú pre k ≥ 1 v intervale [-(2/(Sqrt[n])), (2/(Sqrt[n]))], 
reziduá  modelu AR(1) teda  sú ešte  korelované.

Portmanteauov test

n/4

Portmanteauovým testom otestujeme nulovosť autokorelačnej funkcie reziduí pre prvých h hodnôt posunutia naraz.
 Zvolíme h = n/4 = 210, teda kritickú hodnotu počítame pre df = 209 (máme p =1, q = 0)

Q =PortmanteauStatistic[rez1,210]
Quantile[ChiSquareDistribution[210-1],0.95]

Pretože hodnota testovacej štatistiky Q je menšia ako kritická hodnota, model AR(1)  je vhodný pre časový rad reziduí po systematických zložkách.

b) Testovanie normality  (štandardizovaných) rezíduí Jaques - Berovym testom :

mom[a_,j_]:=1/n*Sum[a[[t]]^j,{t,n}]

srez=(reza-μ)/σ;
sik=Sqrt[n/6.]*Sqrt[mom[srez,3]^2/mom[srez,2]^3]
spi=Sqrt[n/24.]*(mom[srez,4]/mom[srez,2]^2-3)

Quantile[NormalDistribution[0,1],0.99]

TS=sik^2+spi^2

Quantile[ChiSquareDistribution[2],0.99]

Reziduá nemôžeme považovať za náhodnu premennú s normálnym rozdelenim .

c) Testovanie podmienenej heteroskedasticity:

Model prietoku  s AR(1) modelom rezíduí         -  toto nazaver 18.3.

Generovanie náhodnej zložky  modelu  z  s dlžkou  n + 12 (aj  pre   predpovede):

nr=NormalDistribution[0,1];
z=Table[Mean[RandomReal[nr,100]],{n+12}];

Vytvorenie modelu rezíduí po systematických zložkách:

nah=TimeSeries[ar,Take[z,n],Take[y,1]];

Pripočítanie modelu rezíduí  k systematickým zložkám:

syst=Table[sz,{t,n}];
model=syst+nah

Zobrazenie grafu  modelu priemerných mesačných prietokov:

Celý postup opakovať pre dalšie štyri modely

Length[y]