% Shrnuti k ARMAX modelum
y = [.....]; %data
alfa = 0.05; % zvolene riziko

th=armax(y,[2 3]); % odhad arma(2,3) modelu bez ur. konst.
present(th); % vysledky odhadu

% nebo ARMAX(2,3) s urovnovou konstantou
dat=iddata(y,ones(size(y))); % funkce iddata, skladame data (simulovana) a vektor jednicek, ktery se nasobi s parametrem urovnove konstanty delta
th=armax(dat,[2 1 3 0]);     % odhad parametru, zde 2. vstup je [na nb nc nk], nk je zpozdeni na vstupu, trochu jina syntax nez v pripade procesu bez urovnove konstanty
present(th);

% vlastnosti odhadu
[p_sig,test,stud] = par_sig(th,alfa) 
[sig,tests,sk,QQ,chi]=res_sig(th,dat,alfa)

% stabilita systemu
figure
pzmap(th)

%-----------------------------------
% predikce napr. pro model ARMA(1,4)
%-----------------------------------
% vytazeni vektoru parametru z promenne "th", kde mame vysledky odhadu z funkce armax
param = th2par(th); 
a = param(1); %AR %cast ma 1 par.
c = param(2:end); %MA cast ma 4 par.

hp = 20; % horizont predpovedi
%n=n-1; %pokud mame zdiferencovana data, nezapomeneme snizit "n" o jednicku

%PE Computes prediction errors.
%E = PE(Model,DATA)
e=pe(y,th4);

% pro vypocet je treba 4 predchozi hodnoty, vypocteme je bokem, zbytek v cyklu, 
% ocekavana stredni hodnota sumu epsilon je nulova, proto z vypoctu predikci postupne mizi
% a zustava jen clen AR casti
yp(1,:) = a(1)*y(n)  + c(1)*e(n) + c(2)*e(n-1) + c(3)*e(n-2) + c(4)*e(n-3);
yp(2,:) = a(1)*yp(1) + c(2)*e(n) + c(3)*e(n-1) + c(4)*e(n-2);
yp(3,:) = a(1)*yp(2) + c(3)*e(n) + c(4)*e(n-1);
yp(4,:) = a(1)*yp(3) + c(4)*e(n);

for g=5:hp
    yp(g,:) = a(1)*yp(g-1);
end
