Odhad parametrov modelu SETAR:

Ú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[]];

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

data=ReadList["SVKINFLACIA.txt"];
data=Flatten[data];
n=Length[data];
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[p1_,p2_,d_]:=Module[{k},
k=Max[p1,p2,d];
Return[k];]

ypsilonc[t_,p_,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],pocet1++},
{x[[t]]=ypsilonc[t,p,c,d]*phic[p,c,d],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_]:=Module[{sigma,k},
k=parameterk[p,p,d];
sigma=(1/(n-k))*Sum[rezidua[t,p,c,d,n]^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,pocet=0,t,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),
{
AppendTo[rozptyly,rozptyl[p,c,d,n]];
AppendTo[cecka,c];
}
];
minrozptyl=Min[rozptyly];
minc=cecka[[Flatten[Position[rozptyly,minrozptyl]]]];
minc=minc[[1]];
rezidua1={};
rezidua2={};
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]}];];
σ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]];
pmax=6;

samotmý výpočet

setar[pmax,hmin,hmax,n]

Testovanie linearity oproti nelinearite typu SETAR

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

LRtest=ReadList["LR_test.def"];

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

InformationMatrix[