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;

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

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"]];
n=Length[test];

Zadefinovať funkcie pre:

indikačnú funkciu

Indication[A_]:=Module[{},If[TrueQ[A],Return[1],Return[0]];]

premenné pre prvý a druhý režim

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

skeleton (deterministická časť modelu)

skeleton[c_,d_,t_]:=Module[{},
Return[Flatten[Join[ypsilon[t,p1]*(1-Indication[test[[t-d]],c]),ypsilon[t,p2]*Indication[test[[t-d]],c]]]];]

parameterk[p1_,p2_,d_]:=Module[{},Return[Max[p1,p2,d]];]

odhad autoregresných parametrov

phic[c_,d_]:=Module[{ϕ,k=parameterk[p1,p2,d]},
Return[PseudoInverse[Sum[Outer[Times,skeleton[c,d,t],skeleton[c,d,t]],{t,k+1,n}]].Sum[skeleton[c,d,t]*test[[t]],{t,k+1,n}]];]

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

rozptyl[c_,d_]:=Module[{k=parameterk[p1,p2,d]},
Return[(1/(n-k))*Sum[(test[[t]]-phic[c,d].skeleton[c,d,t])^2,{t,k+1,n}]];]

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

setar[d_,t_,phi1_,phi2_,c_]:=Module[{xxx},
xxx=phi1.ypsilon[t,p1]*(1-Indication[test[[t-d]],c])+phi2.ypsilon[t,p2]*Indication[test[[t-d]],c];
Return[xxx];]

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]//N;]

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]//N;]

Výpočet

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

p1=6;
p2=1;

cmin=rozsahc[n][[1]];
cmax=rozsahc[n][[2]];
pmax=Max[p1,p2];

samotný výpočet

Vypocet[pmax_,cmin_,cmax_]:=Module[{p,d,k,rozptyly,cecka,c,minrozptyl,minc,vektor1,vektor2,t,pocet1=0,pocet2=0,σ1,σ2},
For[p=1,p≤pmax,p++,
{For[d=1,d≤p,d++,
{k=parameterk[p,p,d];
rozptyly={};cecka={};
For[c=cmin,c≤cmax,c=c+((cmax-cmin)/10),
{AppendTo[rozptyly,rozptyl[c,d]];
AppendTo[cecka,c];}];
minc=cecka[[Flatten[Position[rozptyly,Min[rozptyly]]]]][[1]];
vektor1={};vektor2={};
For[t=1,t≤n,t++,
If[test[[t]]≤c,
vektor1=Append[vektor1,test[[t]]],
vektor2=Append[vektor2,test[[t]]]];];
pocet1=Length[vektor1];
pocet2=Length[vektor2];
σ1=StandardDeviationMLE[vektor1];
σ2=StandardDeviationMLE[vektor2];
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];}];}];]

Vypocet[pmax,cmin,cmax]

Cvičenie 3
Testovanie linearity oproti nelinearite typu SETAR

InformationMatrix[test,ARModel[Phi1,Roz1]]

ARModel[{phi10,phi11,...phi1p},Roz1]

Phi1 su parametre prveho rezimu a Roz1 je rozptyl, co sme pocitali pri AIC a BIC

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

rezim1

p=;
d=;
c=;
phi1=phic

rezim2

p=;
d=;
c=;

rezim3

p=;
d=;
c=;

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

InformationMatrix[test,ARModel[

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