<<"timeseri\\datasmoo.m";
<<"timeseri\\timeseri.m";
<<"timeseri\\userfunc.m";

SetDirectory[NotebookDirectory[]];
SetOptions[ListPlot,PlotRange->All,BaseStyle->{11,FontFamily->"Arial"}];
popis=Table[{60+i*120,1936+i*10},{i,0,6}];
data=ReadList["PLN.txt"];
data=Flatten[data];
n=Length[data];
gdata=ListPlot[data,Ticks->{popis,Automatic}, AxesOrigin->{0,0}, AxesLabel->{"time [month]","data [(m^3)/s]"},PlotLabel->"PLN",Joined->True]

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

skus=Take[data,-5]
gskus=ListPlot[skus,PlotLabel->"priemerny kurz euro/skk",AxesLabel->{"time[month]","priemerny kurz"},Joined->True,AxesOrigin->{0,0}]
test=Take[data,n-5]
gtest=ListPlot[test,PlotLabel->"priemerny kurz euro/skk",AxesLabel->{"time[month]","priemerny kurz"},Joined->True]

Save["s_t_vzorka.def",skus]
Save["s_t_vzorka.def",test]
Save["s_t_vzorka.def",gskus]
Save["s_t_vzorka.def",gtest]

Mean[test]
StandardDeviation[test]
Variance[test]
cortv=CorrelationFunction[test, 25];
partv= PartialCorrelationFunction[test,25];
myplotcorr1[cortv,2/Sqrt[n],AxesLabel->{k,(Subscript[ρ, k])},PlotLabel->"Autocorrelation Function",PlotRange->All]
myplotcorr1[partv,2/Sqrt[n], AxesLabel->{k,(Subscript[ϕ, kk])},PlotLabel->"Partial Autocorrelation Function",PlotRange->All]

armodels = LevinsonDurbinEstimate[test, 20];
dis=Map[#[[-1]]&,armodels];
ListPlot[dis,Joined->True,AxesLabel->{"p","sigma2"}]

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

a1=Flatten[Position[Table[AIC[mod[[i]],n],{i,1,20}],Min[Table[AIC[mod[[i]],n],{i,1,20}]]]]
a2=Flatten[Position[Table[BIC[mod[[i]],n],{i,1,20}],Min[Table[BIC[mod[[i]],n],{i,1,20}]]]];
If[a1==a2,Print[mod[[a1]]]]

Print["Vysledny model je AR2"]

<<"LinearRegression`"

indikacna[x_]:=If[x==True,1,0]

(*1 a 2 rezim*)
Y1[t_,p_]:=Prepend[Table[test[[t-i]],{i,1,p}],1]
Y2[t_,p_]:=Prepend[Table[test[[t-i]],{i,1,p}],1]

(*SKELETON - DETERMINISTICKA cast*)
Yt[d_,t_,c_,p_]:=
Flatten[Join[Y1[t,p]*indikacna[test[[(t-d)]]<=c],Y2[t,p]*indikacna[test[[(t-d)]]>c]]];

phiF[c_,d_,n_,k_,p_]:=
Module[{},Return[PseudoInverse[(∑_(t=k+1)^n)Outer[Times,Yt[d,t,c,p], Yt[d,t,c,p]]].(∑_(t=k+1)^n)Yt[d,t,c,p]*test[[t]]];]

AIC[p1_,p2_,N1_,N2_,σ1_,σ2_]:=N1*Log[σ1^2]+N2*Log[σ2^2]+2*(p1+1)+2*(p2+1)
BIC[p1_,p2_,N1_,N2_,σ1_,σ2_]:=N1*Log[σ1^2]+N2*Log[σ2^2]+(p1+1)/Log[N1]+(p2+1)/Log[N2]

rozsahC[]:=Module[{tvSort,k1,k2,hmin,hmax},tvSort=Sort[test];
k1=Round[0.075*Length[test]];
k2=Round[0.925*Length[test]];
hmin=tvSort[[k1]];
hmax=tvSort[[k2]];
Print["Minimum parametra c, cmin = ",hmin];
Print["Maximum parametra c, cmax = ",hmax];];

SETAR[pmax_,cmin_,cmax_]:=
Module[
{n,p,d,k,rozptyly,cecka,c,pp,phi,sorted,σ1,C,CC,prvy,druhy,i,N1,N2,sigma1,sigma2,Aic,Bic,optim},
n=Length[test];
For[p=1,p≤pmax,p++,(*.vonkajší cyklus je pre rád p (od 1 po pmax)*)
{

For[d=1,d≤p,d++,(*cyklus pre parameter oneskorenia d (od 1 po p)*)
{
k=Max[p,d];


rozptyly={};
cecka={};


For[c=cmin,c<=cmax,c=c+(cmax-cmin)/100,(*cyklus je pre parameter c od hmin po hmax*)
{
phi=phiF[c,d,n,k,p];
pp=(1/(n-(k+1)))(∑_(t=k+1)^n)(test[[t]]-phi.  Yt[d,t,c])^2;
rozptyly=Append[rozptyly,pp];
cecka=Append[cecka,c];
}];
sorted=Sort[Transpose[{rozptyly,cecka}]];
σ1=sorted[[1,1]];
C=sorted[[1,2]];
CC[d,p]=sorted[[1,2]];
prvy={};
druhy={};
For[i=1,i≤n,i++,
{
If[test[[i]]≤c,prvy=Append[prvy,test[[i]]],druhy=Append[druhy,test[[i]]]]
}];
N1=Length[prvy];N2=Length[druhy];
sigma1=StandardDeviationMLE[prvy];
sigma2=StandardDeviationMLE[druhy];
Aic[d,p]=AIC[p,p,N1,N2,sigma1,sigma2];
Bic[d,p]=BIC[p,p,N1,N2,sigma1,sigma2];
}];
optim=Sort[Table[{Bic[d,p],Aic[d,p],CC[d,p],d},{d,1,p}]];
Print["c = ", optim[[1,3]], "   d =",optim[[1,4]], "   p = ",p,   "      BIC = ",optim[[1,1]],"    AIC = ",optim[[1,2]] ]
}]
]

rozsahC[]

SETAR[5,3.5374,4.4686]

(*Testovanie linearity*)

c=3.5374`;
d=2;
p1=p2=p=5;
k=Max[{p,d}];

prvy={};druhy={};
Do[
If[test[[i]]≤c,prvy=Append[prvy,test[[i]]],druhy=Append[druhy,test[[i]]]],
{i,Length[test]}];
N1=Length[prvy];
N2=Length[druhy];
Print[Style["prvý režim: ",Darker[Blue]],N1,Style[" \ndruhý režim: ",Darker[Blue]],N2]

ϕ1=Take[phiF[c,d,Length[test],k,p],{1,p+1}]
ϕ2=Take[phiF[c,d,Length[test],k,p],{p+2,2p+2}]
φ=Join[ϕ1,ϕ2]

km=Inverse[(∑_(t=k+1)^(Length[test]))Outer[Times,Yt[d,t,c,p],Yt[d,t,c,p]]]

sch=Table[(Sqrt[((km[[i,i]])/(Length[test]))]),{i,Length[km]}];

Column[Transpose[{φ,sch}]]
ϵ=Table[test[[t]]-φ.Yt[d,t,c,p],{t,k+1,Length[test]}];n1=Length[ϵ];
s2=Sum[ϵ[[t]]^2,{t,n1}];
sig=s2/n1
σ=Sqrt[sig]

SETARmodel[t_,d_,p1_,p2_]:=ϕ1.Y1[t,p1]*indikacna[test[[t-d]]≤c]+ϕ2.Y2[t,p2]*indikacna[test[[t-d]]>c]

Model=Join[Take[test,k],Table[SETAR[t,d,p,p],{t,k+1,Length[test]}]];
graf2=ListPlot[Model,PlotRange->{Min[test]-1,Max[test]+1},PlotStyle->{Red},Joined->True];
Show[gtest,graf2]