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)

test=Take[data,n-5];
skus=Take[data,-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]

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,20];
myplotcorr1[corrtest,2/Sqrt[n],AxesLabel->{k,(Subscript[ρ, k])},PlotLabel->"Autocorrelation Function"]

parcorrtest=PartialCorrelationFunction[test,20];
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)

nn=10;
mod=LevinsonDurbinEstimate[test,nn];
dis=Map[#[[-1]]&,mod];
ListPlot[dis,Joined->True,AxesLabel->{"p","(σ^2)"}]

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}]]]];
pmax=Max[p1,p2]+1;

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=Read["test.def"];*)
<<test.def
n=Length[test];

n

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_,p1_,p2_]:=Module[{},
Return[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_]:=Module[{ϵ,i},
ϵ=Table[RandomReal[NormalDistribution[0,1]],{i,n}];
Return[Phi.skeleton[c,d,t,p1,p2]];
Print[Phi.skeleton[c,d,t,p1,p2]];]

odhad autoregresných parametrov

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

informačné kritériá AIC a BIC

kriteriumAIC[p1_,p2_,pocet1_,pocet2_,σ1_,σ2_]:=Module[{},
Return[(pocet1*Log[σ1]+2*(p1+1)+pocet2*Log[σ2]+2*(p2+1))//N];]
kriteriumBIC[p1_,p2_,pocet1_,pocet2_,σ1_,σ2_]:=Module[{},
Return[(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=parameterk[p,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}];
Print["p: ",p," d: ",d," c: ",minc," AIC: ",AIC," BIC: ",BIC];
}];}];]

Vypocet[pmax,cmin,cmax]

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"];

Samotný postup

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

autoregresneparametre[p_,c_,d_]:=Module[{ϕ,ϕ1,ϕ2},
ϕ=phic[c,d,p];
ϕ1={ϕ[[1]],ϕ[[2]]};
ϕ2={ϕ[[3]],ϕ[[4]]};
Return[{ϕ,ϕ1,ϕ2}];]

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

InformationMatrix[p_,c_,d_]:=Module[{k=parameterk[p,p,d]},Return[(∑_(t=k+1)^n)Outer[Times,skeleton[c,d,t,p,p],skeleton[c,d,t,p,p]]];]

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

pocetnost[c_,d_,p_]:=Module[{t,k=parameterk[p,p,d],pocet1=0,pocet2=0},
For[t=k+1,t≤n,t++,If[test[[t-d]]≤c,pocet1++,pocet2++]];
Return[{pocet1,pocet2}];]

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

rozptylyyy[c_,d_,p_,Phi_]:=Module[{t,k=parameterk[p,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,σσ}];]

6.  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)]

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

1. Model

vectormodel={};
minBIC=Sort[vektorBIC,#1[[1]]<#2[[1]]&][[1]];
p=minBIC[[2]];
d=minBIC[[3]];
c=minBIC[[4]];
Print["p = ",p];
Print["d = ",d];
Print["c = ",c];
ϕ=autoregresneparametre[p,c,d][[1]];
ϕ1=autoregresneparametre[p,c,d][[2]];
ϕ2=autoregresneparametre[p,c,d][[3]];
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]]];
gtest=ListPlot[test,Ticks->{popis,Automatic}, AxesOrigin->{0,0}, AxesLabel->{"time [month]","priemerná inflácia"},PlotLabel->"vykreslenie modelu a jeho zobrazenie s originálnymi dátami",Joined->True];gpred=ListPlot[vectormodel,Ticks->{popis,Automatic}, AxesOrigin->{0,0},PlotStyle->Red, AxesLabel->{"time [month]","priemerná inflácia"},Joined->True];
Show[{gtest,gpred}]
Print["Informačná matica: ",MatrixForm[InformationMatrix[p,c,d]]];
Print["početnosť prvého režimu je ",pocetnost[c,d,p][[1]]];
Print["početnosť druhého režimu je ",pocetnost[c,d,p][[2]]];
Print["celkový rozptyl = ",s21=rozptylyyy[c,d,p,ϕ][[6]]];
Print["pvalue = ",pV[LR[LevinsonDurbinEstimate[test,10][[p]][[2]],rozptylyyy[c,d,p,ϕ][[6]]]]];

2. Model

vectormodel={};
minBIC=Sort[vektorBIC,#1[[1]]<#2[[1]]&][[2]];
p=minBIC[[2]];
d=minBIC[[3]];
c=minBIC[[4]];
Print["p = ",p];
Print["d = ",d];
Print["c = ",c];
ϕ=autoregresneparametre[p,c,d][[1]];
ϕ1=autoregresneparametre[p,c,d][[2]];
ϕ2=autoregresneparametre[p,c,d][[3]];
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]]];
gtest=ListPlot[test,Ticks->{popis,Automatic}, AxesOrigin->{0,0}, AxesLabel->{"time [month]","priemerná inflácia"},PlotLabel->"vykreslenie modelu a jeho zobrazenie s originálnymi dátami",Joined->True];gpred=ListPlot[vectormodel,Ticks->{popis,Automatic}, AxesOrigin->{0,0},PlotStyle->Red, AxesLabel->{"time [month]","priemerná inflácia"},Joined->True];
Show[{gtest,gpred}]
Print["Informačná matica: ",MatrixForm[InformationMatrix[p,c,d]]];
Print["početnosť prvého režimu je ",pocetnost[c,d,p][[1]]];
Print["početnosť druhého režimu je ",pocetnost[c,d,p][[2]]];
Print["celkový rozptyl = ",s22=rozptylyyy[c,d,p,ϕ][[6]]];
Print["pvalue = ",pV[LR[LevinsonDurbinEstimate[test,10][[p]][[2]],rozptylyyy[c,d,p,ϕ][[6]]]]];

3. Model

vectormodel={};
minBIC=Sort[vektorBIC,#1[[1]]<#2[[1]]&][[3]];
p=minBIC[[2]];
d=minBIC[[3]];
c=minBIC[[4]];
Print["p = ",p];
Print["d = ",d];
Print["c = ",c];
ϕ=autoregresneparametre[p,c,d][[1]];
ϕ1=autoregresneparametre[p,c,d][[2]];
ϕ2=autoregresneparametre[p,c,d][[3]];
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]]];
gtest=ListPlot[test,Ticks->{popis,Automatic}, AxesOrigin->{0,0}, AxesLabel->{"time [month]","priemerná inflácia"},PlotLabel->"vykreslenie modelu a jeho zobrazenie s originálnymi dátami",Joined->True];gpred=ListPlot[vectormodel,Ticks->{popis,Automatic}, AxesOrigin->{0,0},PlotStyle->Red, AxesLabel->{"time [month]","priemerná inflácia"},Joined->True];
Show[{gtest,gpred}]
Print["Informačná matica: ",MatrixForm[InformationMatrix[p,c,d]]];
Print["početnosť prvého režimu je ",pocetnost[c,d,p][[1]]];
Print["početnosť druhého režimu je ",pocetnost[c,d,p][[2]]];
Print["celkový rozptyl = ",rozptylyyy[c,d,p,ϕ][[6]]];
Print["pvalue = ",pV[LR[LevinsonDurbinEstimate[test,10][[p]][[2]],rozptylyyy[c,d,p,ϕ][[6]]]]];

Cvičenie 4
Diagnostická kontrola modelu SETAR

Testovanie ostávajúcej nelinearity

Nahrať súbor LR_test.def a vlastný definičný súbor
Spracovať zadanie pre najlepšie modely

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

Funkcie

trojskeleton[c1_,c2_,t_,d_]:=Module[{},
Return[Flatten[Join[
ypsilon[t,1]*Indication[test[[t-d]]≤c1],ypsilon[t,2]*Indication[c1<test[[t-d]]≤c2],ypsilon[t,1]*Indication[test[[t-d]]>c2]]]];]
trojsetar[Phi_,c1_,c2_,t_,d_]:=Module[{},
Return[Phi.trojskeleton[c1,c2,t,d]];]
reziduatrojsetar[Phi_,c1_,c2_,t_,d_]:=Module[{},
Return[test[[t]]-trojsetar[Phi,c1,c2,t,d]];]
rozptyltrojsetar[Phi_,c1_,c2_,d_]:=Module[{k=parameterk[p,p,d]},
Return[(1/(n-k))*(∑_(t=k+1)^n)reziduatrojsetar[Phi,c1,c2,t,d]^2];]
phictrojsetar[c1_,c2_,d_]:=Module[{k=parameterk[p,p,d]},
Return[PseudoInverse[(∑_(t=k+1)^n)Outer[Times,trojskeleton[c1,c2,t,d], trojskeleton[c1,c2,t,d]]].(∑_(t=k+1)^n)trojskeleton[c1,c2,t,d]*test[[t]]];]
kriteriumAICtrojsetar[p_,pocet1_,pocet2_,pocet3_,σ1_,σ2_,σ3_]:=Module[{},
Return[(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_]:=Module[{},
Return[(pocet1*Log[σ1]+pocet2*Log[σ2]+pocet3*Log[σ3]+(p+1)*Log[pocet2]+(p+1)*Log[pocet2]+(p+1)*Log[pocet3])//N];]
vektorBICtrojsetar={};
Vypocettrojsetar[p_,cmin_,cmax_]:=Module[{rozptyly,k,t,cecka1,cecka2,vektor1,vektor2,vektor3,phi,c1,c2,minc1,minc2,minphic,pocet1,pocet2,pocet3,σ1,σ2,σ3},
For[d=1,d≤p,d++,
{k=parameterk[p,p,d];
rozptyly={};cecka1={};cecka2={};
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];
rozptyly=AppendTo[rozptyly,rozptyltrojsetar[phi,c1,c2,d]];
cecka1=AppendTo[cecka1,c1];
cecka2=AppendTo[cecka2,c2];};];};];
minc1=cecka1[[Flatten[Position[rozptyly,Min[rozptyly]]]]][[1]];
minc2=cecka2[[Flatten[Position[rozptyly,Min[rozptyly]]]]][[1]];
minphic=phictrojsetar[minc1,minc2,d];
vektor1={};vektor2={};vektor3={};
For[t=k+1,t≤n,t++,
{If[test[[t-d]]≤minc1,
vektor1=AppendTo[vektor1,reziduatrojsetar[minphic,minc1,minc2,t,d]],
If[minc1<test[[t-d]]≤minc2,
vektor2=AppendTo[vektor2,reziduatrojsetar[minphic,minc1,minc2,t,d]],
vektor3=AppendTo[vektor3,reziduatrojsetar[minphic,minc1,minc2,t,d]]]];}];
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}];
Print["p: ",p," d: ",d," c1: ",minc1," c2: ",minc2," AIC: ",AIC," BIC: ",BIC];}];]
autoregresneparametretrojsetar[c1_,c2_,d_]:=Module[{ϕ,ϕ1,ϕ2},
ϕ=phictrojsetar[c1,c2,d];
ϕ1={ϕ[[1]],ϕ[[2]]};
ϕ2={ϕ[[3]],ϕ[[4]]};
Return[{ϕ,ϕ1,ϕ2}];]
InformationMatrixtrojsetar[c1_,c2_,d_]:=Module[{k=parameterk[p,p,d]},Return[(∑_(t=k+1)^n)Outer[Times,trojskeleton[c1,c2,d,t],trojskeleton[c1,c2,d,t]]];]
pocetnosttrojsetar[c1_,c2_,d_]:=Module[{t,k=parameterk[p,p,d],pocet1=0,pocet2=0,pocet3=0},
For[t=k+1,t≤n,t++,If[test[[t-d]]≤c1,pocet1++,If[c1<test[[t-d]]≤c2,pocet2++,pocet3++]]];
Return[{pocet1,pocet2,pocet3}];]

Výpočet pre jednotlivé modely

1. MODEL

1. odhadnúť parametre trojrežimového modelu (pre dané p)

Vypocettrojsetar[7,cmin,cmax]
vectormodeltrojsetar={};
minBICtrojsetar=Sort[vektorBICtrojsetar,#1[[1]]<#2[[1]]&][[1]];
d=minBICtrojsetar[[2]];
c1=minBICtrojsetar[[3]];
c2=minBICtrojsetar[[4]];
Print["p = ",p];
Print["d = ",d];
Print["c1 = ",c1];
Print["c2 = ",c2];

2. vypočítať hodnotu autoregresných parametrov pre jednotlivé reži-my a pomocou informačnej matice vypočítať ich štandardné chyby

ϕ=autoregresneparametretrojsetar[c1,c2,d][[1]];
ϕ1=autoregresneparametretrojsetar[c1,c2,d][[2]];
ϕ2=autoregresneparametretrojsetar[c1,c2,d][[3]];
For[t=1,t≤n,t++,vectormodeltrojsetar=AppendTo[vectormodeltrojsetar,trojsetar[ϕ,c1,c2,t,d]];];
For[t=1,t≤Max[p,d]+1,t++,vectormodeltrojsetar[[t]]=test[[t]]];
gtest=ListPlot[test,Ticks->{popis,Automatic}, AxesOrigin->{0,0}, AxesLabel->{"time [month]","priemerná inflácia"},PlotLabel->"vykreslenie modelu a jeho zobrazenie s originálnymi dátami",Joined->True];gpredtrojsetar=ListPlot[vectormodeltrojsetar,Ticks->{popis,Automatic}, AxesOrigin->{0,0},PlotStyle->Red, AxesLabel->{"time [month]","priemerná inflácia"},Joined->True];
Show[{gtest,gpredtrojsetar}]
Print["Informačná matica: ",MatrixForm[InformationMatrixtrojsetar[c1,c2,d]]];
Print["početnosť prvého režimu je ",pocetnosttrojsetar[c1,c2,d][[1]]];
Print["početnosť prvého režimu je ",pocetnosttrojsetar[c1,c2,d][[2]]];
Print["početnosť prvého režimu je ",pocetnosttrojsetar[c1,c2,d][[3]]];

3. vypočítať reziduá trojrežimového modelu a ich rozptyl S3

reziduavector={};
For[t=Max[p,d]+1,t<n,t++,reziduavector=AppendTo[reziduavector,reziduatrojsetar[ϕ,c1,c2,t,d]];];
Print["rezíduá = ",reziduavector];
Print["rozptyl rezíduí = ",s3=rozptyltrojsetar[ϕ,c1,c2,d]];

4. vypočítať P - value LR testu, kde prvý parameter je rozptyl reziduí dvojrežimového modelu S2 a druhý rozptyl S3

Print["pvalue = ",pV[LR[s21,s3]]];

2. MODEL

1. odhadnúť parametre trojrežimového modelu (pre dané p)

Vypocettrojsetar[6,cmin,cmax]

vectormodeltrojsetar={};
minBICtrojsetar=Sort[vektorBICtrojsetar,#1[[1]]<#2[[1]]&][[1]];
d=minBICtrojsetar[[2]];
c1=minBICtrojsetar[[3]];
c2=minBICtrojsetar[[4]];
Print["p = ",p];
Print["d = ",d];
Print["c1 = ",c1];
Print["c2 = ",c2];

2. vypočítať hodnotu autoregresných parametrov pre jednotlivé reži-my a pomocou informačnej matice vypočítať ich štandardné chyby

ϕ=autoregresneparametretrojsetar[c1,c2,d][[1]];
ϕ1=autoregresneparametretrojsetar[c1,c2,d][[2]];
ϕ2=autoregresneparametretrojsetar[c1,c2,d][[3]];
For[t=1,t≤n,t++,vectormodeltrojsetar=AppendTo[vectormodeltrojsetar,trojsetar[ϕ,c1,c2,t,d]];];
For[t=1,t≤Max[p,d]+1,t++,vectormodeltrojsetar[[t]]=test[[t]]];
gtest=ListPlot[test,Ticks->{popis,Automatic}, AxesOrigin->{0,0}, AxesLabel->{"time [month]","priemerná inflácia"},PlotLabel->"vykreslenie modelu a jeho zobrazenie s originálnymi dátami",Joined->True];gpredtrojsetar=ListPlot[vectormodeltrojsetar,Ticks->{popis,Automatic}, AxesOrigin->{0,0},PlotStyle->Red, AxesLabel->{"time [month]","priemerná inflácia"},Joined->True];
Show[{gtest,gpredtrojsetar}]
Print["Informačná matica: ",MatrixForm[InformationMatrixtrojsetar[c1,c2,d]]];
Print["početnosť prvého režimu je ",pocetnosttrojsetar[c1,c2,d][[1]]];
Print["početnosť prvého režimu je ",pocetnosttrojsetar[c1,c2,d][[2]]];
Print["početnosť prvého režimu je ",pocetnosttrojsetar[c1,c2,d][[3]]];

3. vypočítať reziduá trojrežimového modelu a ich rozptyl S3

reziduavector={};
For[t=Max[p,d]+1,t<n,t++,reziduavector=AppendTo[reziduavector,reziduatrojsetar[ϕ,c1,c2,t,d]];];
Print["rezíduá = ",reziduavector];
Print["rozptyl rezíduí = ",s3=rozptyltrojsetar[ϕ,c1,c2,d]];

4. vypočítať P - value LR testu, kde prvý parameter je rozptyl reziduí dvojrežimového modelu S2 a druhý rozptyl S3

Print["pvalue = ",pV[LR[s22,s3]]];