Uvodné operácie

Balíky

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

Nastavenie pracovného  adresára

(*doma*)

SetDirectory["G:\\acr\\program"]

(*v skole*)

SetDirectory["F:\\acr\\program"]

Nastavenie podmienok pre ListPlot

SetOptions[ListPlot,Joined->True,PlotRange->All,BaseStyle->{10,FontFamily->"Arial"}];

Načítanie dát

Xt=ReadList["SVKINFLACIA.txt",Number]

Základná analýza časového radu

Vykresliť časový rad s popisom osí

n=Length[Xt]     (*pocet merani*)
popisOsiX=Table[{12*i,1999+i},{i,0,n/12}]
graf=ListPlot[Xt,AspectRatio->1/2,AxesLabel->{rok,"%"},AxesOrigin->{0,0},Ticks->{popisOsiX,Automatic}]

stredná hodnota, smerodajná odchýlka a rozptyl

(*stredna hodnota casoveho radu*)
Mean[Xt]

(*smerodajná odchýlka casoveho radu*)
StandardDeviation[Xt]

(*rozptyl casoveho radu*)
Variance[Xt]

Príprava dát na modelovanie a prognózu

Vytvoriť skúšobnú vzorku (5 posledných dát) a testovaciu vzorku (ostávajúce dáta)

(*Xr=Table[(∑_(j=1)^12)flow[[j+12*i]] , {i,0,n1/12-1}] /12;*)
ntv=n-5;
Xsv=Take[Xt,-5]
Xtv=Take[Xt,ntv]

Vykresliť skúšobnú aj testovaciu vzorku s popisom osí

popisOsiX1=Table[{12*i,1999+i},{i,0,ntv/12}]
grafTV=ListPlot[Xtv,AspectRatio->1/2,AxesLabel->{rok,"%"},AxesOrigin->{0,0},Ticks->{popisOsiX1,Automatic}]

popisOsiX2=Table[{i,"2007/".(i+7)},{i,0,5}]
grafSV=ListPlot[Xsv,AspectRatio->1/2,AxesLabel->{rok,"%"},AxesOrigin->{0,0},Ticks->{popisOsiX2,Automatic}]

Uložiť skúšobnú aj testovaciu vzorku spolu s ich grafmi do definičného súboru

Save["skus_test_vzorka.def",Xsv]
Save["skus_test_vzorka.def",Xtv]
Save["skus_test_vzorka.def",grafTV]
Save["skus_test_vzorka.def",grafSV]

Vypočítať strednú hodnotu, smerodajnú odchýlku a rozptyl testovacej vzorky

(*stredna hodnota testovacej vzorky*)
Mean[Xtv]

(*smerodajná odchýlka testovacej vzorky*)
StandardDeviation[Xtv]

(*rozptyl testovacej vzorky*)
Variance[Xtv]

Vypočítať ACF a PACF pre testovaciu vzorku

cortv=CorrelationFunction[Xtv, 25];
partv= PartialCorrelationFunction[Xtv,25];
myplotcorr1[cortv,2/Sqrt[n],AxesLabel->{k,(Subscript[ρ, k])},PlotLabel->"Autocorrelation Function",PlotRange->All]
myplotcorr1[partv,2/Sqrt[n], AxesLabel->{k,(Subscript[ϕ, kk])},PlotLabel->"Partial Autocorrelation Function",PlotRange->All]

Na základe PACF určiť pomocou LevinsonDurbinEstimate a informačných kritérií AIC, BIC vhodný rád p pre model AR(p)

armodels = LevinsonDurbinEstimate[Xtv, 10]

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

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

(* p max =5*)

Nahrať balíky na analýzu časových radov a definičný súbor

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

<<skus_test_vzorka.def;

Zadefinovať funkcie pre:

indikačnú funkciu

IndFunkcia[x_]:=If[x==True,1,0]

premenné pre prvý a druhý režim

Y1[t_,p1_]:=Prepend[Table[Xtv[[t-i]],{i,1,p1}],1]
Y2[t_,p2_]:=Prepend[Table[Xtv[[t-i]],{i,1,p2}],1]

Xtv

skeleton (deterministická časť modelu)

Yt[d_,t_,c_]:=Flatten[Join[Y1[t,p1]*IndFunkcia[Xtv[[(t-d)]]<=c],Y2[t,p2]*IndFunkcia[Xtv[[(t-d)]]>c]]];

odhad autoregresných parametrov

ϕ[c_,d_]:=PseudoInverse[(∑_(t=k+1)^ntv)Outer[Times,Yt[d,t,c], Yt[d,t,c]]].(∑_(t=k+1)^ntv)Yt[d,t,c]*Xtv[[t]];

model SETAR (skeleton +Eps( t))

SETAR[t_,d_]:=ϕ1.Y1[t,p1]*IndFunkcia[Xtv[[t-d]]≤c]+ϕ2.Y2[t,p2]*IndFunkcia[Xtv[[t-d]]>c]

informačné kritériá AIC a BIC

AIC[p1_,p2_]:=N1*Log[σ1^2]+N2*Log[σ2^2]+2*(p1+1)+2*(p2+1)
BIC[p1_,p2_]:=N1*Log[σ1^2]+N2*Log[σ2^2]+(p1+1)/Log[N1]+(p2+1)/Log[N2]

Určiť minimálnu a maximálnu hranicu pre výber parametra c:

usporiadať dáta (označiť napr.du)

datausp=Sort[Xtv];

vypočítať k1 = Round[0.075*n], k2 = Round[0.925*n]

k1=Round[0.075*ntv];
k2=Round[0.925*ntv];

minimálna hranica pre parameter c je hmin = du[[k1]]

cmin=datausp[[k1]]

maximálna hranica pre parameter c je hmax = du[[k2]]

cmax=datausp[[k2]]

h=(cmax-cmin)/100

Do[p1=p2=p;
Do[
k=Max[{p1,p2,d}];
vc={};
vr={};
Do[fii=ϕ[j,d];
pom=(1/(ntv-(k+1)))(∑_(t=k+1)^ntv)(Xtv[[t]]-fii.  Yt[d,t,j])^2;
vr=Append[vr,pom];
vc=Append[vc,j],{j,cmin,cmax,h}];
usporiadane=Sort[Transpose[{vr,vc}]];
sigma1=usporiadane[[1,1]];
c=usporiadane[[1,2]];
cc[d,p]=usporiadane[[1,2]];
prvy={};druhy={};
Do[If[Xtv[[i]]≤c,prvy=Append[prvy,Xtv[[i]]],druhy=Append[druhy,Xtv[[i]]]],{i,ntv}];N1=Length[prvy];N2=Length[druhy];
σ1=StandardDeviationMLE[prvy];
σ2=StandardDeviationMLE[druhy];
Aic[d,p]=AIC[p1,p2];Bic[d,p]=BIC[p1,p2],
{d,1,p}];
optim=Sort[Table[{Bic[d,p],Aic[d,p],cc[d,p],d},{d,1,p}]];
Print["c = ", optim[[1,3]], "   d =",optim[[1,4]], "   p = ",p,   "      BIC = ",optim[[1,1]],"    AIC = ",optim[[1,2]] ],
{p,1,5}]

c=6.206;
d=1;
p1=p2=p=5;
k=Max[{p,d}];

prvy={};druhy={};
Do[
If[Xtv[[i]]≤c,prvy=Append[prvy,Xtv[[i]]],druhy=Append[druhy,Xtv[[i]]]],
{i,ntv}];
N1=Length[prvy];
N2=Length[druhy];
Print[Style["prvý režim: ",Darker[Blue]],N1,Style[" \ndruhý režim: ",Darker[Blue]],N2]

ϕ1=Take[ϕ[c,d],{1,p+1}]
ϕ2=Take[ϕ[c,d],{p+2,2p+2}]
φ=Join[ϕ1,ϕ2];

km=Inverse[(∑_(t=k+1)^ntv)Outer[Times,Yt[d,t,c],Yt[d,t,c]]];
sch=Table[(Sqrt[((km[[i,i]])/ntv)]),{i,Length[km]}];
Column[Transpose[{φ,sch}]]
ϵ=Table[Xtv[[t]]-φ.Yt[d,t,c],{t,k+1,ntv}];n1=Length[ϵ];
s2=Sum[ϵ[[t]]^2,{t,n1}];
sig=s2/n1
σ=Sqrt[sig]

Model=Join[Take[Xtv,k],Table[SETAR[t,d],{t,k+1,ntv}]];
graf2=ListPlot[Model,PlotRange->{Min[Xtv]-1,Max[Xtv]+1},PlotStyle->{Red},Joined->True]
Show[grafTV,graf2]