Modely s dlhou pamäťou 
pre rieku Poprad- originálne dáta

Zostrojíme modely s dlhou pamäťou pre  rieku Poprad. 
Vypočítame predpovede pre viac typov modelov a porovnáme ich predikčné schopnosti pomocou Diebold - Mariannovho testu

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

Nastavíme adresár, v ktorom budeme pracovať:

SetDirectory["g:\\MPM ACR\data\Poprad"]

Nastavíme podmienky pre vykresľovanie údajov :

SetOptions[ListPlot,Joined->True,PlotRange->All,
AxesLabel->{"time [Month]","flow [(m^3)/s]"},AxesOrigin->{0,0},BaseStyle->{11,FontFamily->"Arial"}];

Originálny časový rad

Načítame dáta :

flow1=ReadList["PopradM.txt",Number];
flow=Drop[flow1,2];
n1=Length[flow]
T1=n1/12

Z originálnych dát vytvoríme testovaciu vzorku (70 rokov)
 a  skúšobnú vzorku (posledných 12 mesiacov). 
 Vypočítame dĺžku časového radu v mesiacoch a rokoch.

sv=Take[flow,-12];
ocr=Drop[flow,-12];
n=Length[ocr]
T2=n/12

Zobrazíme testovaciu vzorku (tento graf priradíme do premennej gd pre porovnanie s modelom)

gd=ListPlot[ocr,Joined->True]

Zobrazíme skúšobnú vzorku (tento graf priradíme do premennej gsv pre porovnanie s predpoveďami)

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

Zadefinujeme funkcie na výpočet predpovedných chýb:

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

Pretože máme originálne dáta, ktoré sú silne sezónne,  vytvorime štandardizovaný časový rad.
 Pre každý kalendárny mesiac vypočítame priemer μ(j) a smerodajnú odchýlku σ(j), ktoré potom použijeme na transformáciu (Subscript[M, i,j])=(((Subscript[Y, ij])-μ(j))/(σ(j))) pre i = 1, ..., (Subscript[T, 1]) , j = 1, ..., 12. 
 Po skončení modelovania musíme urobiť spätnú transformáciu (Subscript[Y, i,j])= (Subscript[M, i,j])  σ(j) + μ(j) pre model aj  predpoveď.

Štandardizovaný časový rad

Y=Transpose[Partition[flow,12]];
μ[i_]:=1/T1*Sum[Y[[i,j]],{j,T1}]
σ[i_]:=Sqrt[Sum[(Y[[i,j]]-μ[i])^2,{j,T1}]/(T1-1)]
M[i_,j_]:=(Y[[i,j]]-μ[i])/σ[i]

Xs1=Flatten[Table[M[i,j],{j,T1},{i,12}]];
Xs=Take[Xs1,n];

gdxs=ListPlot[Xs,AxesLabel->{" ","Stand Flow"},Joined->True]

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

cor[[2]]

spec=Spectrum[Xs];
gs=plotspectrum[spec,PlotRange->All, Joined->True]

Testovanie hypotéz

Budeme testovať štandardizovaný časový rad na hladine významnosti α = 0.05 -> hranice oblasti, 
kde H0 nezamietame sú ⟨ 0.809; 1.862⟩

1. Test pre short-memory

S[a_,m_]:=StandardDeviation[Take[a,m]]
Xp[a_,m_]:=Mean[Take[a,m]]
Q[a_,m_]:=1/S[a,m]*(Max[Table[Sum[a[[i]]-Xp[a,m],{i,k}],{k,m}]]-Min[Table[Sum[a[[i]]-Xp[a,m],{i,k}],{k,m}]])
V[a_,m_]:=Q[a,m]/Sqrt[m]

V1=V[Xs,n/2.]

V2=V[Xs,n/4.]

Jedna hodnota je v intervale, ale druha je  mimo intervalu akceptacie -> aby sme zistili, 
či časový rad je generovaný procesom s  dlhou pamäťou, musíme ešte vypočítať modifikovanú testovaciu štatistiku.

cov[a_,m_]:=CovarianceFunction[a,m];
ω[j_,q_]:=1-j/(q+1)
σ2[a_,m_,q_]:=Variance[Take[a,m]]+2*Sum[ω[j,q]*cov[a,m][[j]],{j,q}]
Q1[a_,m_,q_]:=1/σ2[a,m,q]*(Max[Table[Sum[a[[i]]-Xp[a,m],{i,k}],{k,m}]]-Min[Table[Sum[a[[i]]-Xp[a,m],{i,k}],{k,m}]])
V[a_,m_,q_]:=Q1[a,m,q]/Sqrt[m]

V3=V[Xs,n/4.,n/10.]

V4=V[Xs,n/4.,n/5.]

V5=V[Xs,n/2.,n/10.]

V6=V[Xs,n/2.,n/5.]

Všetky vypočítané testovacie štatistiky sú menšie ako dolný kvantil 0.809 (pre α = 0.05) -> H0, že proces je s krátkou pamäťou zamietame.

Modelovanie časového radu s dlhou pamäťou

V prvom kroku odhadneme Hurstov koeficient H a pomocou neho parameter d pre frakcionálne diferencovanie štandardizovaného časového radu

Hurstov odhad H

Do vektora ind uložíme tie hodnoty k ∈ ⟨3; n⟩, pre ktoré je podiel (n/k) celé číslo (delitele čísla n). 
Tieto potom použijeme pri vytváraní dát pre lineárnu regresiu

ind={};Do[If[IntegerQ[n/k],If[k>2,ind=Append[ind,k]]],{k,n/2}]
Length[ind]

Procedúra pre výpočet RS štatistiky:

data[k_]:=(Y1=Partition[Xs,k];
μ1[m_]:=Mean[Y1[[m]]];
S[m_]:=StandardDeviation[Y1[[m]]];
X1[m_]:=Y1[[m]]-μ1[m];
Z[i_,m_]:=Sum[X1[m][[j]],{j,i}];
R[m_]:=Max[Table[Z[i,m],{i,k}]]-Min[Table[Z[i,m],{i,k}]];
P1[m_]:=R[m]/S[m];
T=n/k;
 λ[k]=1/T*Sum[P1[m],{m,T}])

Zostavíme dvojice {Log k, λ(k)},  ktoré použijeme v lineárnej regresii.

Do[λ[ind[[j]]]=data[ind[[j]]],{j,Length[ind]}]

rd=Table[{Log[ind[[j]]],Log[λ[ind[[j]]]]},{j,Length[ind]}];

reg=Regress[rd,{t},t,RegressionReport->{ParameterTable,AdjustedRSquared,ANOVATable,BestFitParameters}]

Hurstov koeficient H je smernica regresnej priamky, d = H - 0.5

H=(BestFitParameters/.reg)[[2]]
d=H-0.5

Pretože d ∈ (0; 0.5), štandardizovaný časový rad  je generovaný stacionárnym invertibilným procesom s dlhou pamäťou.

Vytvorenie diferencovaného časového radu

V druhom kroku vytvoríme frakcionálne diferencovaný časový rad (((1 - B))^d)(Subscript[X, s]). 
Nekonečnú sumu budeme aproximovať konečným súčtom pre k = 50. 
Vytvoríme diferencovaný časový rad pre celý časový rad. 
Na modelovanie však použijeme len prvých n hodnôt, ktoré označíme dcr.

β[j_]:=Gamma[j-d]/Gamma[-d]/Gamma[j+1]

dcra=Table[Sum[β[j]*Xs1[[t-j]],{j,0,50}],{t,51,n1}];
dcr=Drop[Join[Take[Xs1,50],dcra],-12];
sv1=Take[dcra,-12];
gn=ListPlot[dcr,PlotStyle->RGBColor[0,0,1],Joined->True]

Box - Jenkinsova metodologia

V treťom kroku budeme diferencovaný časový rad modelovať lineárnymi ARMA modelmi. 
Na identifikáciu parametrov modelu použijeme autokorelačnú a parciálnu autokorelačnú funkciu modifikovanú pre časový rad s dlhou pamäťou

ρ[k_]:=Gamma[k+d]*Gamma[1-d]/Gamma[k-d+1]/Gamma[d]
ϕ[k_]:=d/(k-d)

kf=Table[ρ[k],{k,0,30}];
pkf=Table[ϕ[k],{k,0,30}];

myplotcorr1[kf,2/Sqrt[n],AxesLabel->{k,(Subscript[ρ, k])},PlotLabel->"Autocorrelation Function",PlotRange->All]

myplotcorr1[pkf,2/Sqrt[n], AxesLabel->{k,(Subscript[ϕ, kk])},PlotLabel->"Partial Autocorrelation Function",PlotRange->{-0.1,0.1}]

Korelačná funkcia už nie je štatisticky významne rôzna od 0 pre k > 15 a parciálna autokorelačná funkcia pre k > 4.

Použijeme len Hannan Rissanenovu procedúru, v ktorej použijeme hodnoty (Subscript[p, max ])=15 ,   (Subscript[q, max ])= 5 
a necháme vypísať 8 najlepších modelov

hrx=HannanRissanenEstimate[dcr,25,15,5,8]

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,8}]

Nasledujúca časť odteraz závisí od výsledkov Hannan - Rissanenovej procedúry 
a môže sa pre iné časové rady meniť.

Aby sme mali pre všetky modely rovnakú šumovú zložku, vypočítame  ju na začiatku. 
Pre modely ARMA budeme používať rozdelenie N(0, 1)

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

Model AR(1)

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

φ1=ar[[1,1]];
Km1=Inverse[InformationMatrix[dcr,ar]];
sch1=Table[(Sqrt[((Km1[[i,i]])/n)]),{i,Length[Km1]}];
Column[Transpose[{ar[[1]],sch1}]]

reza=Residual[dcr,ar];
mi=Mean[reza]
Variance[reza]
sigma=StandardDeviation[reza]
ListPlot[reza]

Testovanie nahodnosti reziduí po AR(1)

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

Pretože hodnoty autokorelačnej funkcie sú pre k ≥ 1 v intervale [-(2/(Sqrt[n])), (2/(Sqrt[n]))], reziduá 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 = 210, teda kritickú hodnotu počítame pre df = 209 (máme p = 1, q = 0)

n/4

PortmanteauStatistic[reza,210]
Quantile[ChiSquareDistribution[209],0.95]

Pretože Q je menšie ako kritická hodnota, AR( 1) je vhodný model pre diferencovaný časový rad .

y=dcr;
AR[t_]:=φ1*y[[t-1]]

model[t_]:=φ1*y[[t-1]];

Model zobrazíme spolu s diferencovaným časovým radom

mo1=Join[Take[y,1],Table[AR[t]+z[[t]],{t,2,n}]];
gmo=ListPlot[mo1,PlotStyle->RGBColor[1,0,0],Joined->True];
Show[gmo,gn]

Odchýlky medzi diferencovaným časovým radom a jeho modelom:

et=dcr-mo1;
StandardDeviation[et]

Spätná transformácia modelu a zobrazenie s pozorovaným časovým radom:

ψ[k_]:=Gamma[k+d]/Gamma[d]/Gamma[k+1]

cr1=Table[Sum[ψ[k]*mo1[[t-k]],{k,0,50}],{t,51,n}];
cr2=Join[Take[Xs,50],cr1];
xs=Partition[cr2,12];
cr3=Flatten[Table[xs[[i,j]]*σ[j]+μ[j],{i,T2},{j,12}]];
gmo=ListPlot[cr3,PlotStyle->RGBColor[1,0,1],Joined->True];
Show[gd,gmo,AspectRatio->1/3]

et=ocr-cr3;
StandardDeviation[et]

Predpovede

Jednokroková predpoveď

predp1[k1_]:=(y=dcr;
nah1=Join[Take[y,1],Table[model[t],{t,2,n}]];
Do[ For[i=1,i≤k1,i++,z2=Random[nr];pom[i]=model[j]+z2];m1=Sum[pom[i],{i,k1}]/k1;
nah1=Append[nah1,m1];y=nah1,{j,n+1,n+12}];cr=Join[mo1,Take[nah1,-12]];pcr=Table[Sum[ψ[k]*cr[[t-k]],{k,0,50}],{t,n+1,n+12}];pp1=Table[If[pcr[[j]]*σ[j]+μ[j]<0,0,pcr[[j]]*σ[j]+μ[j]],{j,12}])

Timing[L=1;pp1=Table[1,{12}];While[mse[sv,pp1,12]>0.1 && L<100,predp1[20];P[L]=pp1;L=L+1];]

K=L-1

s=Take[Sort[Table[{mse[sv,P[L],12],L},{L,K}]],5]

H=s[[1,2]]

gp=ListPlot[P[H],PlotStyle->Hue[0.05],PlotRange->{0,25},Joined->True];
Show[gp,gsv]

prog=P[H];
mae[sv,prog,12]
mse[sv,prog,12]
rmse[sv,prog,12]

par=P[H];

Model ARMA(1, 1)

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

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

reza=Residual[dcr,arma];
mi=Mean[reza]
Variance[reza]
sigma=StandardDeviation[reza]
ListPlot[reza,Joined->True]

Testovanie nahodnosti reziduí po ARMA(1,1)

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

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

n/4

Portmanteauov test

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 = 208 (máme p = 1, q = 1)

PortmanteauStatistic[reza,210]
Quantile[ChiSquareDistribution[208],0.95]

Pretože Q je menšie ako kritická hodnota, ARMA( 1, 1) je vhodný model pre diferencovaný časovy rad

y=dcr;
ARMA[t_]:=φ1*y[[t-1]]+z[[t]]+θ1*z[[t-1]]

model[t_]:=φ1*y[[t-1]]+z[[t]]+θ1*z[[t-1]];

Model zobrazíme spolu s diferencovaným časovým radom:

mo1=Join[Take[y,1],Table[ARMA[t]+z[[t]],{t,2,n}]];
gmo=ListPlot[mo1,PlotStyle->RGBColor[1,0,0],Joined->True];
Show[gmo,gn]

Odchýlky medzi diferencovaným časovým radom a jeho modelom :

et=dcr-mo1;
StandardDeviation[et]

Spätná transformácia modelu a zobrazenie s pozorovaným časovým radom:

ψ[k_]:=Gamma[k+d]/Gamma[d]/Gamma[k+1]

cr1=Table[Sum[ψ[k]*mo1[[t-k]],{k,0,50}],{t,51,n}];
cr2=Join[Take[Xs,50],cr1];
xs=Partition[cr2,12];
cr3=Flatten[Table[xs[[i,j]]*σ[j]+μ[j],{i,T2},{j,12}]];
gmo=ListPlot[cr3,PlotStyle->RGBColor[1,0,1],Joined->True];
Show[gd,gmo,AspectRatio->1/3]

et=ocr-cr3;
StandardDeviation[et]

Predpovede

Jednokroková predpoveď

predp1[k1_]:=(y=dcr;
nah1=Join[Take[y,1],Table[model[t],{t,2,n}]];
Do[ For[i=1,i≤k1,i++,z2=Random[nr];pom[i]=model[j]+z2];m1=Sum[pom[i],{i,k1}]/k1;
nah1=Append[nah1,m1];y=nah1,{j,n+1,n+12}];cr=Join[mo1,Take[nah1,-12]];pcr=Table[Sum[ψ[k]*cr[[t-k]],{k,0,50}],{t,n+1,n+12}];pp1=Table[If[pcr[[j]]*σ[j]+μ[j]<0,0,pcr[[j]]*σ[j]+μ[j]],{j,12}])

Timing[L=1;pp1=Table[1,{12}];While[mse[sv,pp1,12]>0.1 && L<100,predp1[20];P[L]=pp1;L=L+1];]

K=L-1

s=Take[Sort[Table[{mse[sv,P[L],12],L},{L,K}]],5]

H=s[[1,2]]

gp=ListPlot[P[H],PlotStyle->Hue[0.05],PlotRange->{0,25},Joined->True];
Show[gp,gsv]

prog=P[H];
mae[sv,prog,12]
mse[sv,prog,12]
rmse[sv,prog,12]

parma=P[H];

Predpovede pomocou naivnej metody

Funkcia na vypocet k-krokovej "naivnej predpovede":

lab[k_] := Table[{10 j, 10*j + (k - 1)}, {j, 0, 4}]

naive[k_]:=Module[{sv1,prog},sv1=Drop[sv,k];prog=Drop[sv,-k];
gs=ListPlot[sv1,PlotRange->rozsah,Joined->True];
gp=ListPlot[prog,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]]]]

Vypočítame 1-krokovú pre nasledujúcich 12 mesiacov pomocou  "naivnej" metody. a chyby predpovede mse, rmse a mae :

rozsah={0,25}

naive[1]

Show[gs,gp,Ticks->{lab[1],Automatic}]
naiv = Prepend[Drop[sv,-1], sv[[1]]];

Diebold - Mariano dvojstranny test pre  (α = 0.05)

H0:  (Subscript[e, 1])^2-   (Subscript[e, 2])^2  = 0                  H1:  (Subscript[e, 1])^2-   (Subscript[e, 2])^2  ≠ 0

test1[d1_,d2_,d3_,P_,h_]:=Module[{E1,E2,d,dp,V,DM,nr,qh,qd},g[x_]:=x^2;
E1=d1-d2;E2=d1-d3;d=g[E1]-g[E2];dp=Mean[d];γ[i_]:=1/P*Sum[(d[[t-(h-1)]]-dp)(d[[t-(h-1)-i]]-dp),{t,h+i,P+h-1}];V=γ[0]+2*Sum[γ[i],{i,h-1}];DM=dp/Sqrt[Abs[V/P]];nr=NormalDistribution[0,1];
qh=Quantile[nr,0.975];qd=Quantile[nr,0.025];TS=Which[DM<qd,1,qd≤DM≤qh,0,True,-1]]

Zistenie najlepšieho modelu pre 1-krokovú predpoveď

pocet=3;

k=1;Np=12-(k-1);

Ma[1]=naiv
Ma[2]=par
Ma[3]=parma

Apom={};Do[If[i==j,Apom=Append[Apom,"x"],Apom=Append[Apom,test1[sv,Ma[i],Ma[j],Np,k]]],{i,pocet},{j,pocet}];A=Partition[Apom,pocet];TableForm[A,TableHeadings->Automatic]

Všetky tri modely majú rovnaké predikčné vlastnosti.