Modely s dlhou pamäťou ARFIMA

Nastavenie pracovného adresára a volby pre tlač:

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

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

Model ARFIMA pre originálny časový rad priemerných mesačných prietokov

Načítanie   časového radu priemerných mesačných prietokov:

Vykreslenie  priemerných mesačných prietokov a uloženie grafu do premennej:

1.Vytvorenie štandardizovaného  časového  radu

Originálne dáta priemerných mesačných prietokov sa vyznačujú silne sezónnosťou,  
urobíme preto z neho štandardizovaný časový rad nasledovne:
 - z daného radu vytvoríme 12 radov obsahujúcich hodnoty v jednotlivých mesiacoch roka (Subscript[Y, i,j]), i = 1, ..., T= n/12 , j = 1,..., 12, 
 - pre jednotlivé kalendárne mesiace vypočítame priemer (Subscript[μ, j]) a smerodajnú odchýlku (Subscript[σ, j]), j = 1, ..., 12, 
 - urobíme transformáciu hodnôt (Subscript[M, i,j])=(((Subscript[Y, i,j])-(Subscript[μ, j]))/(Subscript[σ, j])) pre i = 1, ..., T , j = 1, ..., 12,
 - zoradíme  štandardizované hodnoty časového radu v pôvodnom poradí.
 
 Po skončení modelovania musíme urobiť spätnú transformáciu (Subscript[Y, i,j])= (Subscript[M, i,j])  .(Subscript[σ, j]) + (Subscript[μ, j]) .

Y=Transpose[Partition[flow,12]];
T=n/12

μ[i_]:=1/T*Sum[Y[[i,j]],{j,T}]
σ[i_]:=Sqrt[Sum[(Y[[i,j]]-μ[i])^2,{j,T}]/(T-1)]
M[i_,j_]:=(Y[[i,j]]-μ[i])/σ[i]

Xstand=Flatten[Table[M[i,j],{j,T},{i,12}]];

X=Xstand;

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

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

2. Testovanie hypotézy
	H0: rad sa vyznačuje krátkou pamäťou (short-memory)

1. testovacia štatistika (Subscript[V, m]):
, kde
, pre m = n/4, resp. m = n/2

Na hladine významnosti α = 0.05 hypotézu H0 nezamietame, 
ak je hodnota štatistiky (Subscript[V, m]) v intervale ⟨ 0.809; 1.862⟩.

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[X,n/2.]

V2=V[X,n/4.]

Jedna hodnota je v intervale, ale druhá je  mimo intervalu.
Aby sme zistili, či časový rad je generovaný procesom s  dlhou pamäťou, 
treba vypočítať modifikovanú testovaciu štatistiku (Subscript[V, m])(q).

2. testovacia štatistika (Subscript[V, m])(q):
				,      kde		
   
 
 pre m = n/4, resp. m = n/2,  q = n/10, resp. q= n/5

Na hladine významnosti α = 0.05 hypotézu H0 nezamietame, 
ak je hodnota štatistiky (Subscript[V, m])(q) v intervale ⟨ 0.809; 1.862⟩.

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[X,n/4.,n/10.]

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

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

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

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

3. Odhad Hurstovho koeficientu H a parametra diferencovania d

<<"LinearRegression`"

Do vektora ind uložiť  celočíselné delitele k čísla n,  k ∈ ⟨3; n⟩ (podiel (n/k) je celé číslo) . 
Použijeme ich pri vytváraní dát pre lineárnu regresiu hodnôt závisle premennej  ln(N) a závisle premennej (Subscript[ln(R/S), N]). 
Koeficient pri lineárnom člene vypočítaný v tejto regresii je Hurstov koeficient: 
	(Subscript[ln(R/S)), N])=ln(c) + H.ln(N)

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

ind

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

data[k_]:=
(Y1=Partition[X,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}])

Zostaviť 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,BestFitParameters}]

Hurstov koeficient H je smernica regresnej priamky, parameter diferencovania je 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.

4. Vytvorenie frakcionálne diferencovaného časového radu (Subscript[Y, t]) s parametrom diferencovania d

Vytvoríme frakcionálne diferencovaný časový rad (Subscript[Y, t])=(((1 - B))^d)(Subscript[X, stand]). 
Nekonečnú sumu budeme aproximovať konečným súčtom pre k = 50. 
Koeficienty autoregresnej reprezentácie sú označené v programe  β.

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

dcra=Table[Sum[β[j]*X[[t-j]],{j,0,50}],{t,51,n}];
dcr=Drop[Join[Take[X,50],dcra]];

gdcr=ListPlot[dcr,PlotStyle->RGBColor[0,0,1],Joined->True]

5. Model ARMA frakcionálne diferencovaného časového radu s dlhou pamäťou

Na identifikáciu parametrov p,q 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"]

Korelačná funkcia je štatisticky významne rôzna od 0 pre k < 11 
a parciálna autokorelačná funkcia pre k < 4.

Hannan Rissanenova procedúra s hodnotami   (Subscript[k, max ])= 25,  (Subscript[p, max ])=12 ,   (Subscript[q, max ])= 5 ,
 vypísať 10 najlepších modelov  :

hrx=HannanRissanenEstimate[dcr,25,12,5,10]

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

Model AR(1)

Generovanie náhodnej zložky  modelu v počte n  :

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

hrx[[1]]

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

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

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

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

Výpočet rezíduí modelu, odhad ich strednej hodnoty, rozptylu,  smerodajnej odchýlky:

rezar1=Residual[dcr,ar1];
mi=Mean[rezar1]
Variance[rezar1]
sigma=StandardDeviation[rezar1]

Testovanie nulovosti autokorelačnej funkcie   :

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

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

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

n/4

p=1;
q=0;

PortmanteauStatistic[rezar1,n/4]
Quantile[ChiSquareDistribution[n/4-p-q],0.95]

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

Vytvorenie modelu štandardizovaného diferencovaného časového radu:

modcr=TimeSeries[ar1,Take[z,n],Take[dcr,1]];

Zobrazenie modelu spolu s diferencovaným časovým radom:

gmodcr=ListPlot[modcr,PlotStyle->RGBColor[1,0,0],Joined->True];
Show[gdcr,gmodcr]

6. Spätná transformácia modelu : oddiferencovanie a odštandartizovanie

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

T=n/12

cr1=Table[Sum[ψ[k]*modcr[[t-k]],{k,0,50}],{t,51,n}];
cr2=Join[Take[X,50],cr1];
xs=Partition[cr2,12];
cr3=Flatten[Table[xs[[i,j]]*σ[j]+μ[j],{i,T},{j,12}]];
gmo=ListPlot[cr3,PlotStyle->RGBColor[1,0,1],Joined->True];

Length[cr3]

Zobrazenie  modelu s originálnym časovým radom:

Show[gdata,gmo,AspectRatio->1/3]

Diagnostická kontrola výsledného modelu

Vyčísliť základné štatistické charakteristiky rezíduí :

Zobraziť autokorelačnú funkciu  výsledného modelu :

Celý postup opakovať pre dalšie ARMA modely