Cvičenie 1
Príprava dát na prognózu

Ú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["DoS_USA.txt"];
data=Flatten[data];
n=Length[data];
test=Take[data,n-5];
skus=Take[data,-5];
n=Length[test];

test

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

Vykresliť časový rad s popisom osí
Vypočítať strednú hodnotu, smerodajnú odchýlku a rozptyl
Vytvoriť skúšobnú vzorku (5 posledných dát) a testovaciu vzorku (ostávajúce dáta)
Vykresliť skúšobnú aj testovaciu vzorku s popisom osí
Uložiť skúšobnú aj testovaciu vzorku spolu s ich grafmi do definičného súboru
Vypočítať strednú hodnotu, smerodajnú odchýlku a rozptyl testovacej vzorky
Vypočítať ACF a PACF pre testovaciu vzorku
Na základe PACF určiť pomocou LevinsonDurbinEstimate a informačných kritérií AIC, BIC vhodný rád p pre model AR(p)

Analyza[]:=Module[{gdata,μ,μt,σ1,σ2,σ1t,σ2t,gskus,gtest,autocor,parcor,nn=10,dis,graflev,mod,p1,p2,pmax},
gdata=ListPlot[data,Ticks->{popis,Automatic}, AxesOrigin->{0,0}, AxesLabel->{"time [month]","priemerná inflácia"},PlotLabel->"SVK Inflácia",Joined->True];
μ=Mean[data];
σ1=StandardDeviation[data];
σ2=Variance[data];
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];
Save["skus.def",skus];
Save["test.def",test];
μt= Mean[test];
σ1t= StandardDeviation[test];
σ2t = Variance[test];
autocor=myplotcorr1[CorrelationFunction[test,20],2/Sqrt[n],AxesLabel->{k,(Subscript[ρ, k])},PlotLabel->"Autocorrelation Function"];
parcor=myplotcorr1[PartialCorrelationFunction[test,20],2/Sqrt[n],AxesLabel->{k,(Subscript[ρ, k])},PlotLabel->"Partial Correlation Function"];
mod=LevinsonDurbinEstimate[test,nn];
dis=Map[#[[-1]]&,mod];
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}]]]];
pmax=Max[p1,p2]+1;
Do[Print["p = ", p,"   AIC= ",AIC[mod[[p]],n],"   BIC= ",BIC[mod[[p]],n]],{p,nn}];
Return[{gtest,pmax,gskus}];

]
pmax=Analyza[][[2]];
gtest=Analyza[][[1]];
gskus=Analyza[][[3]];

pmax

Cvičenie 2
Odhad parametrov modelu SETAR

Úvodné operácie

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

<<test.def;
<<skus.def;

Zadefinovať funkcie pre:

indikačnú funkciu

Indication[A_]:=If[TrueQ[A],Return[1],Return[0]]

premenné pre prvý a druhý režim

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

skeleton (deterministická časť modelu)

skeleton[c_,d_,t_,p1_,p2_]:=Flatten[Join[ypsilon[t,p1]*Indication[test[[t-d]]≤c],ypsilon[t,p2]*Indication[test[[t-d]]>c]]]

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

setar[d_,t_,Phi_,c_,p1_,p2_]:=Phi.skeleton[c,d,t,p1,p2]

odhad autoregresných parametrov

phic[c_,d_,p_]:=PseudoInverse[(∑_(t=Max[p,d]+1)^n)Outer[Times,skeleton[c,d,t,p,p], skeleton[c,d,t,p,p]]].(∑_(t=Max[p,d]+1)^n)skeleton[c,d,t,p,p]*test[[t]]
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}];]
rezidua[p_,c_,d_,t_,Phi_]:=test[[t]]-setar[d,t,Phi,c,p,p]
rozptyl[c_,d_,p_,Phi_]:=(1/(n-Max[p,d]))*(∑_(t=Max[p,d]+1)^n)rezidua[p,c,d,t,Phi]^2

informačné kritériá AIC a BIC

kriteriumAIC[p1_,p2_,pocet1_,pocet2_,σ1_,σ2_]:=(pocet1*Log[σ1]+2*(p1+1)+pocet2*Log[σ2]+2*(p2+1))//N
kriteriumBIC[p1_,p2_,pocet1_,pocet2_,σ1_,σ2_]:=(pocet1*Log[σ1]+pocet2*Log[σ2]+(p1+1)*Log[pocet2]+(p2+1)*Log[pocet2])//N

Výpočet

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

cmin=rozsahc[n][[1]];
cmax=rozsahc[n][[2]];
n=Length[test];
vektorBIC={};

samotný výpočet

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

pmax

Vypocet[pmax,cmin,cmax]

vektorBIC

Cvičenie 3
Testovanie linearity oproti nelinearite typu SETAR

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

LRtest=Read["LR_test.def"];

Funkcie

1. vypočítať hodnotu autoregresných parametrov pre jednotlivé režimy
pomocou informačnej matice vypočítať štandardné chyby odhadov  autoregresných parametrov
vypočítať počet dát v prvom a druhom režime
vypočítať rezíduá a ich rozptyl

InformationMatrix[p_,c_,d_]:=(∑_(t=Max[p,d]+1)^n)Outer[Times,skeleton[c,d,t,p,p],skeleton[c,d,t,p,p]]
pocetnost[c_,d_]:=Module[{t,pocet1=0,pocet2=0},
For[t=1,t≤n,t++,If[test[[t-d]]≤c,pocet1++,pocet2++]];
Return[{pocet1,pocet2}];]
rozptylyyy[c_,d_,p_,Phi_]:=Module[{t,k=Max[p,d],pocet1=0,pocet2=0,vektor,vektor1,vektor2,σ1,σ2,σσ},
vektor1={};vektor2={};
For[t=k+1,t≤n,t++,
If[test[[t-d]]≤c,
vektor1=AppendTo[vektor1,rezidua[p,c,d,t,Phi]],
vektor2=AppendTo[vektor2,rezidua[p,c,d,t,Phi]]];];
vektor=Flatten[Table[rezidua[p,c,d,t,Phi],{t,k+1,n}]];
σσ=rozptyl[c,d,p,Phi];
pocet1=Length[vektor1];
pocet2=Length[vektor2];
σ1=StandardDeviationMLE[vektor1];
σ2=StandardDeviationMLE[vektor2];
Return[{vektor1,σ1,vektor2,σ2,vektor,σσ}];]

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

LinearModel[p_]:=Module[{},Return[LevinsonDurbinEstimate[test,10][[p]]];]
LR[s1_,s2_]:=n*((s1-s2)/s2)
pV[x_]:=1-If[x<0,1-f1[-x],f1[x]]
f1[x_]:=1+(Sqrt[(x/(2Pi))])*Exp[-(x/8)]+(3/2)*Exp[x]*CDF[NormalDistribution[0,1],-((3(Sqrt[x]))/2)]-((x+5)/2)*CDF[NormalDistribution[0,1],-((Sqrt[x])/2)]

3.  vytvoriť model a zobraziť ho s originálnymi dátami (testovacou vzorkou) + 7.  vypočítať P-value LR testu

Vypocet2setar[vektorBIC_,NoModel_]:=Module[{t,Model1,pVal,vectormodel,p,d,c,ϕ,minBIC,gpred,s21},
vectormodel={};
minBIC=Sort[vektorBIC,#1[[1]]<#2[[1]]&][[NoModel]];
p=minBIC[[2]];
d=minBIC[[3]];
c=minBIC[[4]];
ϕ=phic[c,d,p];
For[t=1,t≤n,t++,vectormodel=AppendTo[vectormodel,setar[d,t,ϕ,c,p,p]];];
For[t=1,t≤Max[p,d]+1,t++,vectormodel[[t]]=test[[t]];];
gpred=ListPlot[vectormodel,Ticks->{popis,Automatic}, AxesOrigin->{0,0},PlotStyle->Red, AxesLabel->{"time [month]","priemerná inflácia"},Joined->True];
s21=rozptylyyy[c,d,p,ϕ][[6]];
pVal=pV[LR[LevinsonDurbinEstimate[test,10][[p]][[2]],rozptylyyy[c,d,p,ϕ][[6]]]];
Model1={p,d,c,ϕ,s21};
Return[{gpred,Model1,pVal}];]

Výpočet

1. Model

navrat11=Vypocet2setar[vektorBIC,1];
Show[gtest,navrat11[[1]]]
model1dvoj=navrat11[[2]];
Print[Style["1.  Modlel2setar : ",16,Red,Bold]];
TableForm[{{navrat11[[2]][[1]],navrat11[[2]][[2]],navrat11[[2]][[3]],navrat11[[3]],If[navrat11[[3]]>0.05,"H0 nezamietam","H0 zamietam - je potrebne pridat rezim"],pocetnost[navrat11[[2]][[3]],navrat11[[2]][[2]]]}},
TableHeadings->{None,{"p","d","c","pValue","OK>?",pocetnost}}]

2. Model

navrat12=Vypocet2setar[vektorBIC,2];
Show[gtest,navrat12[[1]]]
model2dvoj=navrat12[[2]];
Print[Style["2.  Modlel2setar : ",16,Red,Bold]];
TableForm[{{navrat12[[2]][[1]],navrat12[[2]][[2]],navrat12[[2]][[3]],navrat12[[3]],If[navrat12[[3]]>0.05,"H0 nezamietam","H0 zamietam - je potrebne pridat rezim"],pocetnost[navrat12[[2]][[3]],navrat12[[2]][[2]]]}},
TableHeadings->{None,{"p","d","c","pValue","OK>?",pocetnost}}]

3. Model

navrat13=Vypocet2setar[vektorBIC,3];
Show[gtest,navrat13[[1]]]
model3dvoj=navrat13[[2]]
Print[Style["3.  Modlel2setar : ",16,Red,Bold]];
TableForm[{{navrat13[[2]][[1]],navrat13[[2]][[2]],navrat13[[2]][[3]],navrat13[[3]],If[navrat13[[3]]>0.05,"H0 nezamietam","H0 zamietam - je potrebne pridat rezim"],pocetnost[navrat13[[2]][[3]],navrat13[[2]][[2]]]}},
TableHeadings->{None,{"p","d","c","pValue","OK>?",pocetnost}}]

Cvičenie 4
Diagnostická kontrola modelu SETAR

Testovanie ostávajúcej nelinearity
Spracovať zadanie pre najlepšie modely

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

LRtest=Read["LR_test.def"];
<<test.def;

Funkcie

trojskeleton[c1_,c2_,t_,d_,p_]:=Flatten[Join[
ypsilon[t,p]*Indication[test[[t-d]]≤c1],ypsilon[t,p]*Indication[c1<test[[t-d]]≤c2],ypsilon[t,p]*Indication[test[[t-d]]>c2]]]
trojsetar[Phi_,c1_,c2_,t_,d_,p_]:=Phi.trojskeleton[c1,c2,t,d,p]
reziduatrojsetar[Phi_,c1_,c2_,t_,d_,p_]:=test[[t]]-trojsetar[Phi,c1,c2,t,d,p]
rozptyltrojsetar[Phi_,c1_,c2_,d_,p_]:=(1/(n-Max[p,d]))*(∑_(t=Max[p,d]+1)^n)reziduatrojsetar[Phi,c1,c2,t,d,p]^2
phictrojsetar[c1_,c2_,d_,p_]:=PseudoInverse[(∑_(t=Max[p,d]+1)^n)Outer[Times,trojskeleton[c1,c2,t,d,p], trojskeleton[c1,c2,t,d,p]]].(∑_(t=Max[p,d]+1)^n)trojskeleton[c1,c2,t,d,p]*test[[t]]
kriteriumAICtrojsetar[p_,pocet1_,pocet2_,pocet3_,σ1_,σ2_,σ3_]:=(pocet1*Log[σ1]+2*(p+1)+pocet2*Log[σ2]+2*(p+1)+pocet3*Log[σ3]+2*(p+1))//N
kriteriumBICtrojsetar[p_,pocet1_,pocet2_,pocet3_,σ1_,σ2_,σ3_]:=(pocet1*Log[σ1]+pocet2*Log[σ2]+pocet3*Log[σ3]+(p+1)*Log[pocet2]+(p+1)*Log[pocet2]+(p+1)*Log[pocet3])//N
InformationMatrixtrojsetar[p_,c1_,c2_,d_]:=(∑_(t=Max[p,d]+1)^n)Outer[Times,trojskeleton[c1,c2,d,t,p],trojskeleton[c1,c2,d,t,p]]
pocetnosttrojsetar[c1_,c2_,d_]:=Module[{t,pocet1=0,pocet2=0,pocet3=0},
For[t=1,t≤n,t++,If[test[[t-d]]≤c1,pocet1++,If[c1<test[[t-d]]≤c2,pocet2++,pocet3++]]];
Return[{pocet1,pocet2,pocet3}];]
vektorBICtrojsetar={};
Vypocet3setar[NoMod_]:=Module[{s2,p,d,k,rozptyly,cec1,cec2,c1,c2,phi,minc1,minc2,minphic,vektor1,vektor2,vektor3,t,pocet1,pocet2,pocet3,σ1,σ2,σ3,vector,reziduav,minBIC,cc1,cc2,dd,ϕ,gpred,s3,pVal,Model},
Which[NoMod==1,{s2=model1dvoj[[5]],p=model1dvoj[[1]]},
NoMod==2,{s2=model2dvoj[[5]],p=model2dvoj[[1]]},
NoMod==3,{s2=model3dvoj[[5]],p=model3dvoj[[1]]}];
For[d=1,d≤p,d++,
{k=Max[p,d];
rozptyly={};cec1={};cec2={};
For[c2=cmax,c2>cmin,c2=c2-((cmax-cmin)/10),
{For[c1=cmin,c1<c2,c1=c1+((cmax-cmin)/10),
{phi=phictrojsetar[c1,c2,d,p];

rozptyly=AppendTo[rozptyly,rozptyltrojsetar[phi,c1,c2,d,p]];
cec1=AppendTo[cec1,c1];
cec2=AppendTo[cec2,c2];};];};];
minc1=cec1[[Flatten[Position[rozptyly,Min[rozptyly]]]]][[1]];
minc2=cec2[[Flatten[Position[rozptyly,Min[rozptyly]]]]][[1]];
minphic=phictrojsetar[minc1,minc2,d,p];
vektor1={};vektor2={};vektor3={};
For[t=k+1,t≤n,t++,
{If[test[[t-d]]≤minc1,
vektor1=AppendTo[vektor1,reziduatrojsetar[minphic,minc1,minc2,t,d,p]],
If[minc1<test[[t-d]]≤minc2,
vektor2=AppendTo[vektor2,reziduatrojsetar[minphic,minc1,minc2,t,d,p]],
vektor3=AppendTo[vektor3,reziduatrojsetar[minphic,minc1,minc2,t,d,p]]]];}];
pocet1=Length[vektor1];pocet2=Length[vektor2];pocet3=Length[vektor3];
σ1=StandardDeviationMLE[vektor1];σ2=StandardDeviationMLE[vektor2];σ3=StandardDeviationMLE[vektor3];
AIC=kriteriumAICtrojsetar[p,pocet1,pocet2,pocet3,σ1,σ2,σ3];
BIC=kriteriumBICtrojsetar[p,pocet1,pocet2,pocet3,σ1,σ2,σ3];
vektorBICtrojsetar=AppendTo[vektorBICtrojsetar,{BIC,d,minc1,minc2}];};];
vector={};reziduav={};
Which[p==1,
{dd=vektorBICtrojsetar[[1]][[2]],cc1=vektorBICtrojsetar[[1]][[3]],cc2=vektorBICtrojsetar[[1]][[4]]},
p≠1,
{minBIC=Sort[vektorBICtrojsetar,#1[[1]]<#2[[1]]&][[NoMod]],dd=minBIC[[2]],cc1=minBIC[[3]],cc2=minBIC[[4]]}
];
ϕ=phictrojsetar[cc1,cc2,dd,p];
For[t=1,t≤n,t++,vector=AppendTo[vector,trojsetar[ϕ,cc1,cc2,t,dd,p]];];
For[t=1,t≤Max[p,dd]+1,t++,vector[[t]]=test[[t]];];
gpred=ListPlot[vector,Ticks->{popis,Automatic}, AxesOrigin->{0,0},PlotStyle->Red, AxesLabel->{"time [month]","priemerná inflácia"},Joined->True];
For[t=Max[p,dd]+1,t<n,t++,reziduav=AppendTo[reziduav,reziduatrojsetar[ϕ,cc1,cc2,t,dd,p]];];
s3=rozptyltrojsetar[ϕ,cc1,cc2,dd,p];
pVal=pV[Abs[LR[s2,s3]]];
Model={p,dd,cc1,cc2,ϕ,s3};
Return[{gpred,Model,pVal}];]

Výpočet pre jednotlivé modely

1. MODEL

navrat1=Vypocet3setar[1];
Show[gtest,navrat1[[1]]]
model1troj=navrat1[[2]];
Print[Style["1.  Modlel3setar : ",16,Red,Bold]];
TableForm[{{navrat1[[2]][[1]],navrat1[[2]][[2]],navrat1[[2]][[3]],navrat1[[2]][[4]],navrat1[[3]],If[navrat1[[3]]>0.05,"H0 nezamietam","H0 zamietam - je potrebne pridat rezim"],pocetnosttrojsetar[navrat1[[2]][[3]],navrat1[[2]][[2]],navrat1[[2]][[2]]]}},
TableHeadings->{None,{"p","d","c1","c2","pValue","OK>?",pocetnost}}]

2. MODEL

navrat2=Vypocet3setar[2];
Show[gtest,navrat2[[1]]]
model2troj=navrat2[[2]];
Print[Style["2.  Modlel3setar : ",16,Red,Bold]];
TableForm[{{navrat2[[2]][[1]],navrat2[[2]][[2]],navrat2[[2]][[3]],navrat2[[2]][[4]],navrat2[[3]],If[navrat2[[3]]>0.05,"H0 nezamietam","H0 zamietam - je potrebne pridat rezim"],pocetnosttrojsetar[navrat2[[2]][[3]],navrat2[[2]][[2]],navrat2[[2]][[2]]]}},
TableHeadings->{None,{"p","d","c1","c2","pValue","OK>?",pocetnost}}]

3. MODEL

navrat3=Vypocet3setar[3];
Show[gtest,navrat3[[1]]]
model3troj=navrat3[[2]];
Print[Style["3.  Modlel3setar : ",16,Red,Bold]];
TableForm[{{navrat3[[2]][[1]],navrat3[[2]][[2]],navrat3[[2]][[3]],navrat3[[2]][[4]],navrat3[[3]],If[navrat3[[3]]>0.05,"H0 nezamietam","H0 zamietam - je potrebne pridat rezim"],pocetnosttrojsetar[navrat3[[2]][[3]],navrat3[[2]][[2]],navrat3[[2]][[2]]]}},
TableHeadings->{None,{"p","d","c1","c2","pValue","OK>?",pocetnost}}]

Cvičenie 5
Výpočet predpovedí pomocou modelu SETAR

Pre 3 najlepšie modely (podľa výsledkov testu ostávajúcej nelinearity 2-
alebo 3- režimové)

Vypočítať jednokrokové predpovede na 5 časových jednotiek (naivnou
metódou aj metódou Monte Carlo, resp. bootstrap metódou)

Naivná metóda

Naivna2rez[model_,krok_]:=Module[{p,d,c,ϕ,i,pom,predp={}},
p=model[[1]];d=model[[2]];c=model[[3]];ϕ=model[[4]];
pom=Prepend[Table[{test[[i]]},{i,n-p+1,n}],{1}];
Do[{pom=Take[pom,-p-1];
pom[[1]]=1;
predp=AppendTo[predp,ϕ.Flatten[Join[Flatten[pom]*Indication[pom[[-1]]≤c],Flatten[pom]*Indication[pom[[-1]]>c]]]];
pom=AppendTo[pom,predp[[-1]]];
},{i,5}];
Return[predp];]

Naivna3rez[model_,krok_]:=Module[{c1,p,d,c2,ϕ,i,pom,predp={}},
p=model[[1]];d=model[[2]];c1=model[[3]];c2=model[[4]];ϕ=model[[5]];
pom=Prepend[Table[{test[[i]]},{i,n-p+1,n}],{1}];
Do[
{
pom=Take[Flatten[pom],-p-1];
pom[[1]]=1;
predp=AppendTo[predp,ϕ.Flatten[Join[Flatten[pom]*Indication[pom[[-1]]≤c1],
Flatten[pom]*Indication[c1<pom[[-1]]≤c2],
Flatten[pom]*Indication[pom[[-1]]>c2]]]];
pom=AppendTo[pom,predp[[-1]]];
},{i,5}];
Return[predp];]

Monte Carlo resp. bootstrap metóda

Monte2Carlo[model_,krok_]:=Module[{ϵ,pom,σ,p,d,c,ϕ,i,predp={},dlzka},
p=model[[1]];d=model[[2]];c=model[[3]];ϕ=model[[4]];σ=model[[5]];
ϵ=Table[RandomReal[NormalDistribution[0,(Sqrt[σ])]],{i,1,Length[test]}];
pom=Prepend[Table[{test[[i]]},{i,n-p+1,n}],{1}];
Do[{pom=Take[pom,-p-1];pom[[1]]=1;
predp=AppendTo[predp,(1/krok)*(∑_(i=1)^krok){pom[[-1]]*ϵ[[i]]+ϕ.Flatten[Join[Flatten[pom]*Indication[pom[[-1]]≤c],
Flatten[pom]*Indication[pom[[-1]]>c]]]}];
pom=AppendTo[pom,predp[[-1]]];},{i,5}];
Return[Flatten[predp]];]

Monte3Carlo[model_,krok_]:=Module[{ϵ,σ,c1,p,d,c2,ϕ,i,pom,predp={}},
p=model[[1]];d=model[[2]];c1=model[[3]];c2=model[[4]];ϕ=model[[5]];σ=model[[6]];
ϵ=Table[RandomReal[NormalDistribution[0,(Sqrt[σ])]],{i,n}];
pom=Prepend[Table[{test[[i]]},{i,n-p+1,n}],{1}];
Do[{pom=Take[Flatten[pom],-p-1];
pom[[1]]=1;
predp=AppendTo[predp,(1/krok)(∑_(i=1)^krok){pom[[-1]]*ϵ[[i]]+ϕ.Flatten[Join[Flatten[pom]*Indication[pom[[-1]]≤c1],
Flatten[pom]*Indication[c1<pom[[-1]]≤c2],
Flatten[pom]*Indication[pom[[-1]]>c2]]]}];
pom=AppendTo[pom,predp[[-1]]];},{i,5}];
Return[Flatten[predp]];]

Zobraziť ich spolu so skúšobnou vzorkou a vypočítať chyby predpovedí

Naive2rez[model_,krok_]:=Module[{gnaivna,chyba},
chyba=Variance[skus-Naivna2rez[model,krok]];
gnaivna=ListPlot[Naivna2rez[model,krok], AxesOrigin->{0,0}, AxesLabel->{"time [month]","priemerná inflácia"},PlotLabel->"Naivná jednokroková metóda, Model 1",Joined->True,PlotStyle->Red];
Print["chyba predpovede = ",chyba];
Return[{gnaivna,chyba}];
]

Naive3rez[model_,krok_]:=Module[{gnaivna,chyba},
chyba=Variance[skus-Naivna3rez[model,krok]];
gnaivna=ListPlot[Naivna3rez[model,krok], AxesOrigin->{0,0}, AxesLabel->{"time [month]","priemerná inflácia"},PlotLabel->"Naivná jednokroková metóda, Model 1",Joined->True,PlotStyle->Red];
Print["chyba predpovede = ",chyba];
Return[{gnaivna,chyba}];
]

Naivná jednokroková metóda - 2-režimy

Show[gskus,Naive2rez[model1dvoj,5][[1]]]

Show[gskus,Naive2rez[model2dvoj,5][[1]]]

Show[gskus,Naive2rez[model3dvoj,5][[1]]]

Naivná jednokroková metóda - 3-režimy

Show[gskus,Naive3rez[model1troj,5][[1]]]

Show[gskus,Naive3rez[model2troj,5][[1]]]

Show[gskus,Naive3rez[model3troj,5][[1]]]

Monte Carlo resp. bootstrap metóda - 2-režimy

Monte2C[model_,krok_]:=Module[{a=Length[skus],gmonte,chyba},
chyba=Variance[skus-Monte2Carlo[model,krok]];
gmonte=ListPlot[Monte2Carlo[model,krok], AxesOrigin->{0,0}, AxesLabel->{"time [month]","priemerná inflácia"},PlotLabel->"Monte Carlo, Model 1",Joined->True,PlotStyle->Red];
Print["chyba predpovede = ",chyba];
Return[{gmonte,chyba}];
]

Monte3C[model_,krok_]:=Module[{a=Length[skus],gmonte,chyba},
chyba=Variance[skus-Monte3Carlo[model,krok]];
gmonte=ListPlot[Monte3Carlo[model,krok], AxesOrigin->{0,0}, AxesLabel->{"time [month]","priemerná inflácia"},PlotLabel->"Monte Carlo, Model 1",Joined->True,PlotStyle->Red];
Print["chyba predpovede = ",chyba];
Return[{gmonte,chyba}];
]

Show[gskus,Monte2C[model1dvoj,50][[1]]]

Show[gskus,Monte2C[model2dvoj,50][[1]]]

Show[gskus,Monte2C[model3dvoj,50][[1]]]

Monte Carlo resp. bootstrap metóda - 3-režimy

Show[gskus,Monte3C[model1troj,50][[1]]]

Show[gskus,Monte3C[model2troj,50][[1]]]

Show[gskus,Monte3C[model3troj,50][[1]]]

Cvičenie 6
Testovanie linearity oproti nelinearite typu LSTAR

Vo vonkajšom cykle pre p od 1 po pmax
 Určiť vhodný lineárny model AR(p) pomocou lineárnej regresie
a vypočítať reziduálny súčet štvorcov SSR0
V cykle pre oneskorenie d (od 1 po p) priradiť qt = Xt-d
 Urobiť pomocnú regresiu v tvare a
vypočítať reziduálny súčet štvorcov SSR1
 Vypočítať index determinácie R2 = 1 – SSR1/SSR0
 Vypočítať testovaciu štatistiku LM3 = n R2, ktorá má asymptoticky c2(3p)
 Určiť P-hodnotu

lPval={};sPval={};mins={};minl={};
ypsilon[t_,p_]:=Prepend[Table[{test[[t-i]]},{i,1,p}],{1}]
Funkcia[]:=Module[{x,xx,A,p,b,e,AA,A1,i,j,d,qt,LM,SSSR0,SSSR1,thet,R},
For[p=1,p≤pmax,p++,
{
A=Table[0,{i,n-p},{j,p+1}];
For[i=1,i≤n-p,i++,A[[i]]=Flatten[ypsilon[p+i,p]]];
b=Drop[test,p];
x=Transpose[A].b.PseudoInverse[Transpose[A].A];
e=b-A.x;
SSSR0=e.e;
AA=Table[0,{i,n-p},{j,4p+1}];
A1=Take[A,All,{2,p+1}];
For[d=1,d≤p,d++,
{
qt=Take[test,{p+1-d,n-d}];
For[i=1,i<=n-p,i++,AA[[i]]=Flatten[Join[A[[i]],A1[[i]],A1[[i]]*qt[[i]],A1[[i]]*(qt[[i]]^2),Take[A,All,{2,p+1}][[i]]*(qt[[i]]^3)]];
];
xx=PseudoInverse[Transpose[AA].AA].Transpose[AA].e;
thet=e-AA.xx;
SSSR1=thet.thet;
R=1-(SSSR1/SSSR0);
LM=n*R;
sPval=AppendTo[sPval,{p,d,ChiSquarePValue[LM,2p][[2]]}];
lPval=AppendTo[lPval,{p,d,ChiSquarePValue[LM,3p][[2]]}];
};
];};];]
Funkcia[]

minsPval=Sort[sPval,#1[[3]]<#2[[3]]&];


minlPval=Sort[lPval,#1[[3]]<#2[[3]]&];

Print[Style["Najlepšie parametre pre LSTAR: ",16,Red,Bold]];
TableForm[Take[minlPval,5],TableHeadings->{None,{"p","d","pValue"}}]

Print[Style["Najlepšie parametre pre ESTAR: ",14,Red,Bold]];
TableForm[Take[minsPval,5],TableHeadings->{None,{"p","d","pValue"}}]

Cvičenie 7
Odhad parametrov dvojrežimového modelu STAR

Funkcie

1. Definovať prechodove funkcie

LogF[γ_,q_,c_]:=(1/(1+(E^(-γ(q-c)))))
ExpF[γ_,q_,c_]:=1-(E^(-γ(((q-c))^2)))

ypsilonc[γ_,c_,t_,p_,d_]:=Flatten[Join[Flatten[ypsilon[t,p]*(1-LogF[γ,test[[t-d]],c])],ypsilon[t,p]*LogF[γ,test[[t-d]],c]]]
ϕ[γ_,c_,p_,d_]:=PseudoInverse[(∑_(t=Max[p,d]+1)^n)Outer[Times,ypsilonc[γ,c,t,p,d], ypsilonc[γ,c,t,p,d]]].(∑_(t=Max[p,d]+1)^n)ypsilonc[γ,c,t,p,d]*test[[t]]
skelet[γ_,c_,p_,d_,t_,Phi_]:=Phi.ypsilonc[γ,c,t,p,d]
rezid[γ_,c_,p_,d_,t_,Phi_]:=test[[t]]-skelet[γ,c,p,d,t,Phi]
rozpty[γ_,c_,p_,d_,Phi_]:=(1/(n-Max[p,d]))*(∑_(t=Max[p,d]+1)^n)rezid[γ,c,p,d,t,Phi]^2
AICstar[p_,n1_,n2_,σ1_,σ2_]:=n1*Log[σ1]+n2*Log[σ2]+2*(p+1)+2*(p+1)
BICstar[p_,n1_,n2_,σ1_,σ2_]:=n1*Log[σ1]+n2*Log[σ2]+(p+1)*Log[n1]+(p+1)*Log[n2]

OdhadL[p_,d_]:=Module[{γ,c,phi,rozptyly,t,cec,gamma,minc,minγ,minphi,rezi,n1,n2,σ1,σ2,aic=0,bic=0},
rozptyly={};cec={};gamma={};
For[γ=0.5,γ≤1,γ+=.1,
{For[c=cmin,c<cmax,c++,
{
phi=ϕ[γ,c,p,d];
rozptyly=AppendTo[rozptyly,rozpty[γ,c,p,d,phi]];
cec=AppendTo[cec,c];
gamma=AppendTo[gamma,γ];};];};];
minc=cec[[Flatten[Position[rozptyly,Min[rozptyly]]]]][[1]];
minγ=gamma[[Flatten[Position[rozptyly,Min[rozptyly]]]]][[1]];
minphi=ϕ[minγ,minc,p,d];
n1=(∑_(t=p+1)^n)(1-LogF[minγ,test[[t-d]],minc]);
n2=(∑_(t=p+1)^n)LogF[minγ,test[[t-d]],minc];
σ1=(1/n1)(∑_(t=p+1)^n)rezid[minγ,minc,p,d,t,minphi]^2*(1-LogF[minγ,test[[t-d]],minc]);
σ2=(1/n2)(∑_(t=p+1)^n)rezid[minγ,minc,p,d,t,minphi]^2*LogF[minγ,test[[t-d]],minc];
Return[{AICstar[p,n1,n2,σ1,σ2],BICstar[p,n1,n2,σ1,σ2],minc,minγ}];]

OdhadS[p_,d_]:=Module[{γ,c,phi,rozptyly,t,cec,gamma,minc,minγ,minphi,rezi,n1,n2,σ1,σ2,aic=0,bic=0},
rozptyly={};cec={};gamma={};
For[γ=0.5,γ≤1,γ+=.1,
{For[c=cmin,c<cmax,c++,
{
phi=ϕ[γ,c,p,d];
rozptyly=AppendTo[rozptyly,rozpty[γ,c,p,d,phi]];
cec=AppendTo[cec,c];
gamma=AppendTo[gamma,γ];};];};];
minc=cec[[Flatten[Position[rozptyly,Min[rozptyly]]]]][[1]];
minγ=gamma[[Flatten[Position[rozptyly,Min[rozptyly]]]]][[1]];
minphi=ϕ[minγ,minc,p,d];
n1=(∑_(t=p+1)^n)(1-ExpF[minγ,test[[t-d]],minc]);
n2=(∑_(t=p+1)^n)LogF[minγ,test[[t-d]],minc];
σ1=(1/n1)(∑_(t=p+1)^n)rezid[minγ,minc,p,d,t,minphi]^2*(1-ExpF[minγ,test[[t-d]],minc]);
σ2=(1/n2)(∑_(t=p+1)^n)rezid[minγ,minc,p,d,t,minphi]^2*ExpF[minγ,test[[t-d]],minc];
Return[{AICstar[p,n1,n2,σ1,σ2],BICstar[p,n1,n2,σ1,σ2],minc,minγ}];]

Hodnoty

lstar={};estar={};
For[i=1,i≤Length[minl],i++,
lstar=AppendTo[lstar,OdhadL[minl[[i]][[1]],minl[[i]][[2]]]];];
For[i=1,i≤Length[mins],i++,
estar=AppendTo[estar,OdhadS[mins[[i]][[1]],mins[[i]][[2]]]];];

Print[Style["Hodnoty AIC a BIC pre LSTAR: ",16,Red,Bold]];
TableForm[lstar,TableHeadings->{None,{"AIC","BIC","c","γ"}}]
Print[Style["Hodnoty AIC a BIC pre ESTAR: ",16,Red,Bold]];
TableForm[estar,TableHeadings->{None,{"AIC","BIC","c","γ"}}]

Hessove matice a Kovariancna matica

prvygradient[γ_,c_,p_,d_,t_,phi_]:=Module[{dγ,dc,dphi1,dphi2,gam,cc},
dγ=D[skelet[gam,c,p,d,t,phi],gam]/.gam->γ;
dc=D[skelet[γ,cc,p,d,t,phi],cc]/.cc->c;
dphi1=ypsilon[t,p]*(1-LogF[γ,test[[t-d]],c]);
dphi2=ypsilon[t,p]*LogF[γ,test[[t-d]],c];
Return[Flatten[{"1",dγ,"2",dc,"3",dphi1,"4",dphi2}]];]

prvygradient[0.5,2.3,7,1,8,ϕ[0.5,2.3,7,1]]

Grad2Skel[t_,p_,d_,g_,c_,Phic_]:=Module[{gg,gc,gphi1,gphi2},ggg=D[D[Skel[t,p,d,gx,c,Phic],gx],gx]/.gx->g;
ggc=D[D[Skel[t,p,d,gx,cx,Phic],gx],cx]/.gx->g/.cx->c;
gcc=D[D[Skel[t,p,d,g,cx,Phic],cx],cx]/.cx->c;
gcg=D[D[Skel[t,p,d,gx,cx,Phic],cx],gx]/.cx->c/.gx->g;
Return[SparseArray[{{1,1}->ggg[[1,1]],{1,2}->ggc[[1,1]],{2,1}->gcg[[1,1]],{2,2}->gcc[[1,1]]},{2p+4,2p+4}]];];

druhygradient[γ_,c_,p_,d_,t_,phi_]:=Module[{ddγ,dγdc,ddc,dcdγ,ddphi1dc,ddphi1dγ,ddphi2dc,ddphi2dγ,gam,cc},
ddγ=Flatten[D[D[skelet[gam,c,p,d,t,phi],gam],gam]/.gam->γ];
dγdc=D[D[skelet[gam,cc,p,d,t,phi],gam],cc]/.gam->γ/.cc->c;
ddc=D[D[skelet[γ,cc,p,d,t,phi],cc],cc]/.cc->c;
dcdγ=D[D[skelet[gam,cc,p,d,t,phi],cc],gam]/.cc->c/.gam->γ;
ddphi1dc=D[ypsilon[t,p](1-LogF[γ,test[[t-d]],cc]),cc]/.cc->c;
ddphi1dγ=D[ypsilon[t,p](1-LogF[gam,test[[t-d]],c]),gam]/.gam->γ;
ddphi2dc=D[ypsilon[t,p]LogF[γ,test[[t-d]],cc],cc]/.cc->c;
ddphi2dγ=D[ypsilon[t,p]LogF[gam,test[[t-d]],c],gam]/.gam->γ;
Return[SparseArray[{{1,1}->ddγ,{1,2}->dγdc,{2,1}->dcdγ,{2,2}->ddc},{2p+4,2p+4}]];];

prvy stlpec mam prva derivacia podla gammy a druha podla vsetkych ostatnych (1,1-gamagama,2,1-gamac,3,1 gamaphi1....2p+4,1 gamaphi posledne)

druhygradient[0.5,2.3,7,1,8,ϕ[0.5,2.3,7,1]]

MaticaA[γ_,c_,p_,d_,phi_]:=
(2/n)((∑_(t=p+1)^n)(Outer[Times,prvygradient[γ,c,p,d,t,phi],prvygradient[γ,c,p,d,t,phi]]-rezid[γ,c,p,d,t,phi]*druhygradient[γ,c,p,d,t,phi]))

MatrixForm[MaticaA[0.5,2.3,5,1,ϕ[0.5,2.3,5,1]]]

MaticaB[γ_,c_,p_,d_,phi_]:=(4/n)(∑_(t=1)^n)rezid[γ,c,p,d,t,phi]^2*Outer[Times,prvygradient[γ,c,p,d,t,phi],prvygradient[γ,c,p,d,t,phi]]

MaticaB[γ_,c_,p_,d_,phi_]:=
(1/n)(∑_(t=p+1)^n)(Outer[Times,-2rezid[γ,c,p,d,t,phi]*prvygradient[γ,c,p,d,t,phi],-2rezid[γ,c,p,d,t,phi]*prvygradient[γ,c,p,d,t,phi]])

MatrixForm[MaticaB[0.5,2.3,5,1,ϕ[0.5,2.3,5,1]]]

Kovariancna[γ_,c_,p_,d_,phi_]:=PseudoInverse[MaticaA[γ,c,p,d,phi]].MaticaB[γ,c,p,d,phi].PseudoInverse[MaticaA[γ,c,p,d,phi]]

MatrixForm[Kovariancna[0.5,2.3,5,1,ϕ[0.5,2.3,5,1]]]