Lineárne ARMA modely 
pre  časový rad počtu
 novozaregistrovaných áut v USA

Zostrojíme lineárne ARMA modely pre reziduá po systematických zložkách.

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

SetDirectory["F:\\anlyzaCR"]

SetOptions[ListPlot,PlotRange->All,AxesLabel->{"time [Month]","počet * 1000"},
BaseStyle->{11,FontFamily->"Arial"},Joined->True];

NAČÍTANIE REZÍDUÍ,    VYTVORENIE  SKUŠOBNEJ VZORKY

Načítame dáta , reziduá a funkciu popisujúcu systémové zložky modelu:

data=ReadList["AUTA_MES.dat",Number];

Length[data]

rez=<<rezauta.dat ;

sz=<<syzauta.def

Z originálnych dát a z rezíduí po systematických zložkách vytvoríme  testovaciu vzorku (bez posledných 12 mesiacov) a skúšobnú vzorku ( posledných 12 mesiacov).

ocr=Drop[data,-12];
sv=Take[data,-12];

Z reziduí vytvoríme rad bez posledných 12 mesiacov .
Vypočítame dĺžku takto upraveného časového radu reziduí

rez1=Drop[rez,-12];
n=Length[rez1]

Skúšobnú vzorku sv uložíme do  súboru:

Save["auta_sv.dat",sv]

gd=ListPlot[ocr,Joined->True]

rozsah={500,1200}

gsv=ListPlot[sv, PlotRange->{500,1100},AxesOrigin->{0,500},
AxesLabel->{"time [Month]","počet * 1000"},Joined->True]

gdr=ListPlot[rez1,Joined->True,AxesLabel->{"Month","počet * 1000"}]

PRIRADENIE  RADU MODELOVANÝCH  REZÍDUÍ  DO PREMENNEJ  y

y=rez1;

DEFÍNÍCIE  FUNKCIÍ

Definícia funkcie na na vypocet predpov. chyb

mae[a_,F_,k_]:=N[(((∑_(i=1)^k)Abs[a[[i]]-F[[i]]])/k)]
mse[a_,F_,k_]:=N[(((∑_(i=1)^k)(a[[i]]-F[[i]])^2)/k)]
rmse[a_,F_,k_]:=N[(Sqrt[(((∑_(i=1)^k)(a[[i]]-F[[i]])^2)/k)])]
lab[k_]:=Table[{10 j, 10*j+(k-1)},{j,0,4}]

Definicia funkcie na vypocet k - krokovej "naivnej predpovede"

naive[k_,rozsah_]:=Module[{sv1,prog},
sv1=Drop[sv,k];
prog=Drop[sv,-k];
gs=ListPlot[sv1,PlotRange->rozsah,PlotStyle->{Dashing[{0.02,0.02,0.01}]},Joined->True];
Print["mse   ",mse[sv1,prog,Length[sv1]],"rmse  ",rmse[sv1,prog,Length[sv1]],"  mae  ",mae[sv1,prog,Length[sv1]]]]

Vygenerovanie nahodnej zlozky modellu

vygenerovanie nahodnej zlozky s dlzkou radu n+12

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

Odhad parametrov modelu

HANNAN - RISSANENOVA PROCEDURA

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

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

hlav=(If[Head[#1]===ARMAModel,
		Head[#1][Length[#1[[1]]],Length[#1[[2]]]],
		Head[#1][Length[#1[[1]]]]]&)/@hrx;

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

Model MA(1)

Najlepší model je MA(1). 
Spresníme odhad jeho parametrov podmienenou metódou maximálnej vierohodnosti :

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

Vypočítame reziduá modelu, odhad ich strednej hodnoty, disperzie,  smerodajnej odchýlky a zobrazíme graf rezíduí :

reza=Residual[y,mam];
μ=Mean[reza]
Variance[reza]
σ=StandardDeviation[reza]
ListPlot[reza,Joined->True]

Jaques - Berovym testom otestujeme normalitu rezíduí (štandardizovaných):

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]

Rezíduá môžeme považovať za náhodnú premennú s normálnym rozdelením 
(na hladine vyznamnosti 0.01).

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

θ1=mam[[1,1]];
Km1=Inverse[InformationMatrix[y,mam]];φ
sch1=Table[(Sqrt[((Km1[[i,i]])/n)]),{i,Length[Km1]}];
Column[Transpose[{Flatten[{mam[[1]]}],sch1}]]

3. Tretia časť je testovanie vhodnosti modelu.

Testovanie náhodnosti reziduí po MA(1)

Ako prvé otestujeme nulovosť autokorelačnej funkcie reziduí MA modelu:

corrx=CorrelationFunction[reza, 25];
myplotcorr1[corrx,2/Sqrt[n],AxesLabel->{k,(Subscript[ρ, k])},PlotLabel->"Autocorrelation Function",PlotRange->All]

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

Portmanteauov test

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

n/4

Q =PortmanteauStatistic[reza,18]
Quantile[ChiSquareDistribution[17],0.95]

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

Vytvoríme model MA(1) pre reziduá po systematických zložkách,
 pripočítame k nemu systematické zložky a
zobrazíme graf kompletného modelu spolu s grafom  originálnych dát.

nah1=Prepend[TimeSeries[mam,Take[z,n],Take[y,10]],y[[1]]];
syst=Table[sz,{t,n}];
model=syst+nah1;

gmo=ListPlot[model,PlotStyle->RGBColor[1,0,0],Joined->True];

Show[gd,gmo]

Rezíduá po modeli:

et=Take[data,n]-model;
StandardDeviation[et]
ListPlot[et,Joined->True]

Predpovede

Vypocitame 1 kirokovu predpoved na nsledujucich 12 mesiacov pomocou modelu a chyby predpovede nse,rmse,mae

predp=TimeSeries[mam,Take[z,-13]]+Table[sz,{t,n+1,n+12}]

Chyby predpocede :

Print["mse  ",mse[sv,predp,12],"     rmse  ",rmse[sv,predp,12],"mae  ",mae[sv,predp,12]]

gs=ListPlot[sv,PlotRange->rozsah,Joined->True];
gp=ListPlot[predp,PlotRange->rozsah,PlotStyle->{Dashing[{0.02,0.02,0.01}]},Joined->True];

Show[gs,gp]

Jednokrokova predpoved pomocou naivnej predpovede

rozsah={500,1100};
naive[1,rozsah]

Jednokrokova predpoved pomocou prediktora

BestLinearPredictor[y,mam,12,Exact->False][[1]]

xhat=BestLinearPredictor[y,mam,12,Exact->False][[1]];
predpoved=Table[sz,{t,n+1,n+12}]+xhat;
gp=ListPlot[predpoved,PlotStyle->Hue[0.05],PlotRange->{600,1300},Joined->True];
Show[gp,gsv]

chyby predpovede mse rmse, mae

prog=predpoved;
Print["mse  ",mse[sv,predp,12],"     rmse  ",rmse[sv,predp,12],"   mae  ",mae[sv,predp,12]]

3. Tretia časť je testovanie vhodnosti modelu.

Testovanie náhodnosti reziduí po AR(1)

Ako prvé otestujeme nulovosť autokorelačnej funkcie reziduí MA modelu:

corrx=CorrelationFunction[reza, 25];
myplotcorr1[corrx,2/Sqrt[n],AxesLabel->{k,(Subscript[ρ, k])},PlotLabel->"Autocorrelation Function",PlotRange->All]

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

Portmanteauov test

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

n/4

Q =PortmanteauStatistic[reza,18]
Quantile[ChiSquareDistribution[17],0.95]

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

Vytvoríme model MA(1) pre reziduá po systematických zložkách,
 pripočítame k nemu systematické zložky a
zobrazíme graf kompletného modelu spolu s grafom  originálnych dát.

nah1=TimeSeries[arm,Take[z,n],Take[y,1]];

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

gmo=ListPlot[model,PlotStyle->RGBColor[1,0,0],Joined->True];

Show[gd,gmo]

Rezíduá po modeli:

et=Take[data,n]-model;
StandardDeviation[et]
ListPlot[et,Joined->True]

Model ARMA(1,1)

2 Najlepší model je AR(1). 
Spresníme odhad jeho parametrov podmienenou metódou maximálnej vierohodnosti :

arma=ConditionalMLEstimate[y,hrx[[3]]]

Vypočítame reziduá modelu, odhad ich strednej hodnoty, disperzie,  smerodajnej odchýlky a zobrazíme graf rezíduí :

reza=Residual[y,arma];
μ=Mean[reza]
Variance[reza]
σ=StandardDeviation[reza]
ListPlot[reza,Joined->True]

Jaques - Berovym testom otestujeme normalitu rezíduí (štandardizovaných):

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]

Rezíduá môžeme považovať za náhodnú premennú s normálnym rozdelením 
(na hladine vyznamnosti 0.01).

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

φ1=arma[[1,1]];
θ1=arma[[2,1]];
Km1=Inverse[InformationMatrix[y,arma]];
sch1=Table[(Sqrt[((Km1[[i,i]])/n)]),{i,Length[Km1]}];
Column[Transpose[{Flatten[{arma[[1,1]],arma[[2,1]]}],sch1}]]

3. Tretia časť je testovanie vhodnosti modelu.

Testovanie náhodnosti reziduí po AR(1)

Ako prvé otestujeme nulovosť autokorelačnej funkcie reziduí MA modelu:

corrx=CorrelationFunction[reza, 25];
myplotcorr1[corrx,2/Sqrt[n],AxesLabel->{k,(Subscript[ρ, k])},PlotLabel->"Autocorrelation Function",PlotRange->All]

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

Portmanteauov test

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

n/4

Q =PortmanteauStatistic[reza,18]
Quantile[ChiSquareDistribution[17],0.95]

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

Vytvoríme model MA(1) pre reziduá po systematických zložkách,
 pripočítame k nemu systematické zložky a
zobrazíme graf kompletného modelu spolu s grafom  originálnych dát.

nah1=TimeSeries[arma,Take[z,n+1],Take[y,1]];

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

gmo=ListPlot[model,PlotStyle->RGBColor[1,0,0],Joined->True];

Show[gd,gmo]

Rezíduá po modeli:

et=Take[data,n]-model;
StandardDeviation[et]
ListPlot[et,Joined->True]