Cvičenie 1

Úvodné operácie

Nahratie balíkov na analýzu časových radov

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

Nastavenie pracovného  adresára

SetDirectory[NotebookDirectory[]];

Nastavenie podmienok pre ListPlot

SetOptions[ListPlot,PlotRange->All,BaseStyle->{11,FontFamily->"Arial"}];
popis=Table[{6+i*12,1999+i*1},{i,0,10}];

Načítanie dát (v prípade potreby odstránenie prvých dvoch údajov)

data=ReadList["SVKINFLACIA.txt"];
data=Flatten[data];
n=Length[data];

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

Vykresliť časový rad s popisom osí

gdata=ListPlot[data,Ticks->{popis,Automatic}, AxesOrigin->{0,0}, AxesLabel->{"time [month]","priemerná inflácia"},PlotLabel->"SVK Inflácia",Joined->True]

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

μ=Mean[data]
σ=StandardDeviation[data]
σ2=Variance[data]

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)

skus=Take[data,-5];

test=Take[data,n-5];

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

gskus=ListPlot[skus,Ticks->{popis,Automatic}, AxesOrigin->{0,0}, AxesLabel->{"time [month]","priemerná inflácia"},PlotLabel->"skusobna vzorka",Joined->True]

gtest=ListPlot[test,Ticks->{popis,Automatic}, AxesOrigin->{0,0}, AxesLabel->{"time [month]","priemerná inflácia"},PlotLabel->"testovacia vzorka",Joined->True]

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

Save["skus.def",skus]
Save["test.def",test]

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

μ = Mean[test]
σ = StandardDeviation[test]
σ2 = Variance[test]

Vypočítať ACF a PACF pre testovaciu vzorku

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

parcorrtest=PartialCorrelationFunction[test,40];
myplotcorr1[parcorrtest,2/Sqrt[n],AxesLabel->{k,(Subscript[ρ, k])},PlotLabel->"Partial Correlation Function"]

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

mod=LevinsonDurbinEstimate[test,20];

dis=Map[#[[-1]]&,mod];

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

nn=20;

mod=LevinsonDurbinEstimate[test,nn];
Do[Print["p = ", p,"   AIC= ",AIC[mod[[p]],n],"   BIC= ",BIC[mod[[p]],n]],{p,nn}]

p1=Flatten[Position[Table[AIC[mod[[i]],n],{i,1,nn}],Min[Table[AIC[mod[[i]],n],{i,1,nn}]]]];
p2=Flatten[Position[Table[BIC[mod[[i]],n],{i,1,nn}],Min[Table[BIC[mod[[i]],n],{i,1,nn}]]]];
If[p1==p2,Print["Vysledny model je AR"p1]Print["Vysledny model je AR",Max[p1,p2]]];
pmax=Max[p1,p2];

Cvičenie 2
Odhad parametrov modelu SETAR

Úvodné operácie

Načítanie dát a definičných súborov

skus=Flatten[ReadList["skus.def"]];
test=Flatten[ReadList["test.def"]];

Zadefinovať funkcie pre:

indikačnú funkciu

Indication[A_,B_]:=Module[{cc},
If[A>B,cc=1,cc=0];
Return[cc];]

premenné pre prvý a druhý režim

ypsilon[t_,p_]:=Module[{y,i},
y=Prepend[Table[{test[[t-i]]},{i,1,p}],{1}];
Return[y];]

parameterk[p_,p_,d_]:=Module[{k},
k=Max[p,p,d];
Return[k];]

ypsilonc[p_,t_,c_,d_]:=Module[{y,yy},
y=ypsilon[t,p];
If[Indication[test[[t-d]],c],
y=Flatten[{y*0,y*1},1],
y=Flatten[{y*1,y*0},1]];
Return[y];]

phic[p_,c_,d_,n_]:=Module[{ϕ,k},
k=parameterk[p,p,d];
ϕ=Inverse[Sum[ypsilonc[t,p,c,d].Transpose[ypsilonc[t,p,c,d]],{t,k+1,n}]].Sum[ypsilonc[t,p,c,d]*test[[t]],
{t,k+1,n}];
Return[ϕ];]

skeleton (deterministická časť modelu)

skeleton[p_,d_,c_,n_]:=Module[{t,x,pocet1=0,pocet2=0},
x=Table[0,{i,n}];
For[t=0,t<n,t++,
If[Indication[test[[t-d]],c],
{x[[t]]=ypsilonc[t,p,c,d]*phic[p,c,d,n],pocet1++},
{x[[t]]=ypsilonc[t,p,c,d]*phic[p,c,d,n],pocet2++}];
Return[{x,pocet1,pocet2}];
]]

odhad autoregresných parametrov

rozsahc[n_]:=Module[{k1,k2,du,hmin,hmax},
du=Sort[test];
k1=Round[0.075*n];
k2=Round[0.925*n];
hmin=du[[k1]];
hmax=du[[k2]];
Return[{hmin,hmax}];]

rezidua[t_,p_,c_,d_,n_,Phic_]:=Module[{ec},
ec=test[[t]]-Transpose[Phic].ypsilonc[t,p,c,d];
Return[ec];]

rozptyl[p_,c_,d_,n_,Phic_]:=Module[{sigma,k},
k=parameterk[p,p,d];
sigma=(1/(n-k))*Sum[rezidua[t,p,c,d,n,Phic]^2,{t,k+1,n}];
Return[sigma];]

informačné kritériá AIC a BIC

kriteriumAIC[p1_,p2_,pocet1_,pocet2_,σ1_,σ2_]:=Module[{},
AIC=pocet1*Log[σ1]+2*(p1+1)+pocet2*Log[σ2]+2*(p2+1);
Return[{AIC}];]

kriteriumBIC[p1_,p2_,pocet1_,pocet2_,σ1_,σ2_]:=Module[{},
BIC=pocet1*Log[σ1]+(p1+1)*Log[pocet1]+pocet2*Log[σ2]+(p2+1)*Log[pocet2];
Return[{BIC}];]

model SETAR (skeleton + (Subscript[ϵ, t]))

SETAR[pmax_,hmin_,hmax_,n_]:=Module[{p,c,d,k,minc,minphi,σ1,σ2,cecka,pocet1=0,pocet2=0,t,pphhii,rozptyly,minrozptyl,rezidua1,rezidua2},
For[p=1,p≤pmax,p++,
{
For[d=1,d≤p,d++,
{
k=parameterk[p,p,d];
rozptyly={};
cecka={};
For[c=hmin,c≤hmax,c=c+((hmax-hmin)/10),
{
Print["1"];
pphhii=phic[p,c,d,n];
Print[pphhii];
AppendTo[rozptyly,rozptyl[p,c,d,n,pphhii]];
AppendTo[cecka,c];
}
];
Print[rozptyly];
minrozptyl=Min[rozptyly];
minc=cecka[[Flatten[Position[rozptyly,minrozptyl]]]];
Print[minc];
minc=minc[[1]];
Print[minc];
rezidua1={};
rezidua2={};
Print["6"];
minphi=phic[p,minc,d,n];
For[t=k+1,t≤n,t++,
If[Indication[test[[t-d]],minc]==0,
{rezidua1=AppendTo[rezidua1,rezidua[t,p,minc,d,n,minphi]];pocet1++;rezidua1=Flatten[rezidua1]},
{rezidua2=AppendTo[rezidua2,rezidua[t,p,minc,d,n,minphi]];pocet2++;rezidua2=Flatten[rezidua2]}];];
Print["7"];
σ1=Variance[rezidua1];
σ2=Variance[rezidua2];
AIC=kriteriumAIC[p,p,pocet1,pocet2,σ1,σ2];
BIC=kriteriumBIC[p,p,pocet1,pocet2,σ1,σ2];
Print["p: ",p," d: ",d," c: ",minc," AIC: ",AIC," BIC: ",BIC];
}
];
}
];
]

Výpočet

zadanie premenných, pre ktoré budeme rátať

n=Length[test];
hmin=rozsahc[n][[1]];
hmax=rozsahc[n][[2]];

poznámka: Ja budem rátať pre pmax = 3, lebo nám stačia len tri hodnoty a mne sa AIC a BIS s rastúcimi parametrami výrazne menia

pmax=3;

samotmý výpočet

SETAR[pmax,hmin,hmax,n]

Cvičenie 3
Testovanie linearity oproti nelinearite typu SETAR

Nahrať súbor LR_test.def a vlastný definičný súbor

LRtest=ReadList["LR_test.def"];

Samotný postup

1.   vypočítať hodnotu autoregresných parametrov pre jednotlivé režimy

vezmem si p=1, d=1, c=2.3, AIC=-123.761, BIC=-116.247

p=1;d=1;c=2.3;AIC=-123.761;BIC=-116.247;

phic[p,c,d,n]

2. pomocou informačnej matice vypočítať štandardné chyby odhadov  autoregresných parametrov

rozptyl[5,1,1,1]

n=Length[test]

rozptyl[p,c,d,n]

InformationMatrix[

3.   vypočítať počet dát v prvom a druhom režime

4.   vypočítať rezíduá a ich rozptyl

5.   vytvoriť model a zobraziť ho s originálnymi dátami (testovacou vzorkou)

6. odhadnúť parametre lineárneho AR(p) modelu s absolútnym členom (pre príslušné p) a vypočítať rozptyl jeho rezíduí

7.  vypočítať P-value LR testu