%========================================================================
%         BOX-JENKINSOVY ARMA MODELY PRO STACIONARNI CASOVE RADY
% Systematicky postup identifikace, odhadu parametru a overeni modelu MA
%              (vyzaduje SYSTEM IDENTIFICATION TOOLBOX)
%========================================================================

% (C) 1996-98, V. Vesely, Masaryk University of Brno (Czech Republic)

format compact;
echo on
% Bez ujmy na obecnosti se muzeme omezit na stabilni kauzalni a invertibilni
% proces ARMA(p,q) s nulovou stredni hodnotou, ktery vygenerujeme pomoci
% procedury 'armagen'.
% Z teorie totiz vyplyva, ze kazdy stacionarni proces lze s libovolnou
% presnosti (mereno shodou autokorelacni funkce, tj. momentu 2.radu)
% modelovat prave takovym procesem, pricemz nulove stredni hodnoty lze
% vzdy dosahnout vycentrovanim.
%
% Invertibilita znaci, ze proces lze vlastne modelovat jako AR proces,
% ktery ciste teoreticky muze byt az nekonecneho radu.
% Proto analyzu kazdeho neznameho procesu muzeme zahajit pokusem
% identifikovat jej jako AR proces postupem podrobne rozebranym jiz v
% 'crcv12'. Pokud vzniknou problemy s urcenim radu, tj. kriterialni ztratove
% funkce (napr. AIC, MDL apod.) nebudou mit prukazne minimum uvnitr i
% pomerne sirokeho rozsahu testovanych radu (budou se stale jevit spise
% jako klesajici), pak je toto signalem k volbe modelu MA, resp. smiseneho
% modelu ARMA.
%
% Pri analyze nyni budeme davat prednost funkcim ze SYSTEM IDENTIFICATION
% TOOLBOXu, ktere nabizi sirsi moznosti i komfortnejsi pouzivani,
% nebot redukuji mnozstvi vstupnich a vystupnich informaci jejich
% koncentrovanim do specialnich datovych formatu (jako 'theta', 'zepo'
% apod.), o jejichz obsah se uzivatel pri bezne praci nemusi prilis starat.
pause;
%-----------------------------------------------------------------------------
% Zjisteni pritomnosti SYSTEM IDENTIFICATION TOOLBOXu:
if ~exist('ident')==2
   disp('CHYBA: SYSTEM IDENTIFICATION TOOLBOX NENALEZEN');
   return;
end
%---------------------------------------------------------
% DEMONSTRACE ODHADU PARAMETRU SIMULOVANEHO PROCESU MA(2)
%=========================================================
% 1. GENEROVANI PROCESU
%    ------------------
  q=2; n=500;
% Rozptyl, resp. smerodatna odchylka generujiciho gaussovskeho bileho sumu
  sigma2=2.25; sigma=sqrt(sigma2);
% Generovani MA(2) procesu s parametry:
  THETA = -[0.5;0.2];
  [x,phir,thetar]=armagen([],THETA,2*n,sigma);
  thetar, abs(thetar)
% Koreny polynomu THETA(z) lezi vne jednotkoveho kruhu,
% takze proces je invertibilni
  x  = x - mean(x); % Pro jistotu data vycentrujeme
  xv = x(n+1:2*n);  % 2. usek dat pouzijeme pro overeni modelu
  x  = x(1:n);      % 1. usek dat pouzijeme pro odhad modelu
% Proverime nahodnost
  xrnd = trand(x)
% Alespon jeden z testu by mel zamitnout nahodnost (vraci 0).
% Pokud tomu tak neni, je x bily sum a nema smysl pokracovat.
% V takovem pripade spustte davku znovu.
if xrnd
   disp('PREDCASNY KONEC: Vygenerovana data se jevi jako nahodny sum');
%==================================KONEC======================================
echo off
return;
end
pause;
%-----------------------------------------------------------------------------
% 2. IDENTIFIKACE PROCESU
%    --------------------
% a) ODHAD STREDNI HODNOTY A AUTOKOVARIANCNI FUNKCE PROCESU (procedura acf)
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  [xro,xroE,mi,miE,xacf]=acf(x,figure);
  mi                      % Stredni hodnota a jeji konfidencni intervaly
  [mi-miE(1),mi+miE(1)]   % pouze za predpokladu stacionarity
  [mi-miE(2),mi+miE(2)]   % presnejsi za predpokladu stacionarity a normality
% - Nula by mela lezet v intervalu spolehlivosti s polosirkou malou relativne
%   k velikosti dat, kdy nezamitame hypotezu o nulove stredni hodnote.
% - Pro h>2 lezi acf(h) v nulovem konfidencnim pasu a vykazuje
%   zde nahodne chovani.
%   Podle ocekavani tedy usuzujeme na MA proces radu q=2.
pause;
%-----------------------------------------------------------------------------
% b) ODHAD PARCIALNI AUTOKORELACNI FUNKCE A AUTOREGRESNICH KOEFICIENTU phi
%    (procedura pacf)
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  [xpacf,phi]=pacf(xacf,n,figure);
% Obalka parcialni autokorelacni funkce klesa exponencialne k nule, pricemz
% v konfidencnim pasu nulovosti nema nahodny charakter.
% Proto predbezne usuzujeme, ze se asi nejedna o AR proces konecneho radu.
pause;
%-----------------------------------------------------------------------------
close all
% c) POKUS O ODHAD RADU p PRO MODEL AR(p)
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Spocteme odhady rozptylu sigma2 generujiciho bileho sumu pro sirsi rozsah
% radu, napr. pro p=1:10:
  pmax = 10;
  V = arxstruc(x,x,(1:pmax).'); % Ztratova funkce jen z x
  figure; p = selstruc(V)       % graf ztratove funkce, v p nabyva minima
pause;
%-----------------------------------------------------------------------------
% Vykreslime grafy nekolika kriterialnich ztratovych funkci soucasne
  rp = (1:pmax)./n; AIC = 2; MDL = log(n);
  vaic = V(1,1:pmax).*(1+ AIC*rp);
  vmdl = V(1,1:pmax).*(1+ MDL*rp);
  plot(V(1,1:pmax),'w-');
  hold on; plot(vaic,'g--'); plot(vmdl,'g:');
  title('v(p): Ztratova funkce (-), AIC (--) a MDL (:)');
  xlabel('rad p');
  hold off;
  figure(gcf);
% a spocteme minima
  p = selstruc(V,'aic')  % p, kde AIC nabyva minima
  p = selstruc(V,'mdl')  % p, kde MDL nabyva minima
  p = selstruc(V,3)      % p, kde kriterium s vahou c=3 nabyva minima
pause;
%-----------------------------------------------------------------------------
% Predchozi postup zopakujeme pro presnejsi odhady z nezavisleho
% pozorovani xv.
  VV = arxstruc(x,xv,(1:pmax).');  % Ztratova funkce uzitim xv
  figure; p = selstruc(VV) % graf ztratove funkce, v p nabyva minima
pause;
%-----------------------------------------------------------------------------
% Vykreslime grafy nekolika kriterialnich ztratovych funkci soucasne
  vvaic = VV(1,1:pmax).*(1+ AIC*rp);
  vvmdl = VV(1,1:pmax).*(1+ MDL*rp);
  plot(VV(1,1:pmax),'w-');
  hold on; plot(vvaic,'g--'); plot(vvmdl,'g:');
  title('vv(p): Ztratova funkce (-), AIC (--) a MDL (:)');
  xlabel('rad p');
  hold off;
% a spocteme minima
  p = selstruc(VV,'aic') % p, kde AIC nabyva minima
  p = selstruc(VV,'mdl') % p, kde MDL nabyva minima
  p = selstruc(VV,3)     % p, kde kriterium s vahou c=3 nabyva minima
% Body zlomu ztratovych funkci, resp. minima kriterialnich funkci jsou
% nevyrazna, davaji spise vyssi navzajem se lisici hodnoty.
% To potvrzuje nase ocekavani o neadkvatnosti AR modelu, ktere naznacil
% jiz prubeh acf.
% V dalsim se pokusime najit optimalni ARMA model.
pause;
%-----------------------------------------------------------------------------
close all
% d) ODHAD RADU p,q PROCESU ARMA(p,q), q>0
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Pro ARMA model IDENT bohuzel nenabizi (ponekud nelogicky) funkce,
% ktere by byly analogii funkci 'arxstruc' a 'selstruc' pro AR model.
% Pro odhad parametru modelu ARMA vsak muzeme pouzit funkci 'armax',
% kterou budeme opakovane volat pro vsechny kombinace p,q z vhodne
% zvoleneho rozsahu radu, napr. pro p=0:5, q:1:5.
% Pomoci 'th2par', resp. 'pe' pak spocteme odhady rozptylu sigma2
% generujiciho bileho sumu:
% v(p)  ... odhad rozptylu spocteny z tehoz useku x uziteho pro odhad modelu
% vv(p) ... odhad rozptylu spocteny z nezavisleho useku xv
  pmax = 5; qmax = 5;
% Vypocet dvojic [p,q] v uvedenem rozsahu:
  pq  = [kron((0:5),ones(1,5));kron(ones(1,6),(1:5))]
  Npq = sum(pq);  % pocty odhadovanych parametru p+q pro kazdou dvojici [p,q]
  [Npq,I] = sort(Npq); % pq setridime vzestupne podle poctu parametru
  pq = pq(:,I); [pq;Npq]
  pq = pq.';
  kmax = size(pq,1)    % pocet dvojic
% Provedte help 'armax'
keyboard;
%-----------------------------------------------------------------------------
% Odhad v(k) z useku x a vv(k) z useku xv pro kazdou dvojici pq(k,:):
  for k=1:kmax
      th = armax(x,pq(k,:));
      [par,P,v(k)] = th2par(th);
%      z = pe(x ,th); v(k)  = std(z).^2;
      z = pe(xv,th); vv(k) = std(z).^2;
echo off
  end
echo on
% Vykreslime prubeh v(k) a prislusnych kriterialnich funkci AIC a MDL:
  rp   = Npq ./n;
  vaic = v.*(1+ AIC*rp);
  vmdl = v.*(1+ MDL*rp);
  vc3  = v.*(1+   3*rp);
  plot(v,'w-');
  hold on; plot(vaic,'g--'); plot(vmdl,'g:');
  title('v(k): Ztratova funkce (-), AIC (--) a MDL (:)');
  xlabel('k');
  hold off;
pause;
% a spocteme minima
  [vm,I]    = min(v)   ; pq(I,:), vm     % p,q kde v   nabyva minima
  [vaicm,I] = min(vaic); pq(I,:), vaicm  % p,q kde AIC nabyva minima
  [vmdlm,I] = min(vmdl); pq(I,:), vmdlm  % p,q kde MDL nabyva minima
  [vc3m ,I] = min(vc3) ; pq(I,:), vc3m   % p,q kde vc3 nabyva minima
pause;
%-----------------------------------------------------------------------------
% Podobne vykreslime prubeh vv(k) a prislusnych kriterialnich funkci
% AIC a MDL:
  rp   = Npq ./n;
  vvaic = vv.*(1+ AIC*rp);
  vvmdl = vv.*(1+ MDL*rp);
  vvc3  = vv.*(1+   3*rp);
  figure; plot(vv,'w-');
  hold on; plot(vvaic,'g--'); plot(vvmdl,'g:');
  title('vv(k): Ztratova funkce (-), AIC (--) a MDL (:)');
  xlabel('k');
  hold off;
pause;
% a spocteme minima
  [vvm,I]    = min(vv)   ; pq(I,:), vvm     % p,q kde vv   nabyva minima
  [vvaicm,I] = min(vvaic); pq(I,:), vvaicm  % p,q kde AIC  nabyva minima
  [vvmdlm,I] = min(vvmdl); pq(I,:), vvmdlm  % p,q kde MDL  nabyva minima
  [vvc3m ,I] = min(vvc3) ; pq(I,:), vvc3m   % p,q kde vvc3 nabyva minima
% Grafy i minima kriterialnich funkci dobre urcuji bod zlomu a tedy
% rady pro p=0 a q=2.
pause;
%-----------------------------------------------------------------------------
close all
% 3. ODHAD PARAMETRU PROCESU IDENTIFIKOVANEHO JAKO MA(2)
%    ---------------------------------------------------
  q = 2;
  th = armax(x,[0,q]);
  present(th);
  [par,P,LAM]=th2par(th)
% kde
  theta = par.'
% je odhad presnych hodnot
  THETA
% spolu s polosirkami 95%-ho konfidencniho intervalu
  thetaE = 1.96*sqrt(diag(P))
% a
  v = LAM
% je odhad skutecneho rozptylu
  sigma2
pause;
%-----------------------------------------------------------------------------
% Koreny (zeros) a poly ARMA modelu odpovidaji v IDENT prevracenym hodnotam
% korenu polynomu porade theta(z) a phi(z).
% Jejich odhady spolu se smerodatnymi odchylkami vraci
% ve specialnim formatu 'zepo' funkce 'th2zp'.
% 'zepo' = znaci ZEros & POles
% Pokud koreny phi(z) lezi vne jednotkoveho kruhu (tj. v IDENT uvnitr
% jednotkoveho kruhu), pak proces je kauzalni a stabilni, nebot
% 1/phi(z) je pak holomorfni v uzavrenem jednotkovem kruhu a da se zde
% rozvest do absolutne konvergentni Taylorovy rady.
% Podobne koreny theta(z) vne jednotkoveho kruhu (tj. v IDENT uvnitr
% jednotkoveho kruhu) zarucuji invertibilitu procesu.
% Rozmisteni korenu (o) a polu (x) a tedy invertibilitu, resp. kauzalnost
% a stabilitu systemu nejsnaze overime z grafu, ktery kresli funkce 'zpplot'.
% Pro MA proces se zobrazi pouze koreny, nebot phi(z)=1 nema koreny.
% Provedte help th2zp a help zpplot.
keyboard;
%-----------------------------------------------------------------------------
% Tedy koreny jsou ulozeny v
  zepo = th2zp(th)
% pricemz prvy sloupec submatice
  zp = zepo(2:q+1,1:2)
% jsou skutecne (az na poradi) prave prevracene hodnoty korenu
% polynomu theta(z)
  THETAr = roots(flipud([1;THETA])) % koreny pro presne parametry THETA
  thetar = roots(flipud([1;theta])) % koreny pro odhad theta
  zer = 1./thetar                   % koreny dle IDENT
% a navic druhy sloupec zp pro tyto koreny udava smerodatne odchylky
% a tedy polosirky 95% konfidencniho intervalu ziskame jako
  1.96*zp(:,2)
% Rozmisteni korenu vykresli funkce 'zpplot':
  figure; zpplot(zepo)
% Vidime, ze koreny skutecne lezi uvnitr jednotkoveho kruhu.
pause;
%-----------------------------------------------------------------------------
% Na zaver porovname odhad s presnymi hodnotami (1.sloupec):
% Parametry:
  [THETA,theta]
% Polosirky 95%-ho intervalu spolehlivosti:
  [zeros(q,1),thetaE]
% Rozptyl bileho sumu:
  [sigma2,v]
% Vidime dobrou shodu odhadu se skutecnosti.
pause;
%-----------------------------------------------------------------------------
close all
% 4. OVERENI MODELU ODHADNUTEHO JAKO MA(2)
%    -------------------------------------
% a) Nejprve zkontrolujeme, zda hodnoty koncovych parametru jsou
%    vyznamne (jejich konfidencni intervaly neobsahuji nulu):
  [theta-thetaE,theta+thetaE]
% Pokud by napriklad [theta(q)-thetaE(q),theta(q)+thetaEqp)] obsahoval nulu,
% znamenalo by to, ze theta(q) se nelisi vyznamne od nuly a ze tedy
% rad modelu je pravdepodobne nadhodnocen.
% Cely proces odhadu bychom pak zopakovali pro rad q-1 a znovu overili.
pause;
%-----------------------------------------------------------------------------
% b) Proverime, zda odhadnuty model a jeho komponenty splnuji teoreticke
%    predpoklady, za nichz byl nalezen.
%    V nasem pripade staci spocist tzv. rezidualni sum, tj. zbytkovy
%    sum, ktery by mel korespondovat s generujicim bilym sumem.
  [xs,z,vz] = armasim(x,[],theta);  % uzijeme-li usek, z nehoz byl provaden
                                    % odhad modelu
% Nebo jinak uzitim standardni funkce MATLABu 'filter':
  zz = filter(1,[1;theta],x); xxs = x-zz; vzz = std(z).^2;
  max(abs(z-zz)), max(abs(xs-xxs)), vz-vzz
% Provedte help filter.
keyboard;
%-----------------------------------------------------------------------------
  [xvs,zv,vzv] = armasim(xv,[],theta); % uzijeme-li nezavisly usek
  zzv = filter(1,[1;theta],xv); xxvs = xv-zzv; vzzv = std(zzv).^2;
  max(abs(zv-zzv)), max(abs(xvs-xxvs)), vzv-vzzv
% Proverime, vlastnost IID rezidulaniho sumu napr. pomoci 'trand'
% (existuji i specialni testy pouze na nekorelovanost, tj. jen na bily sum)
  trand(z)
  trand(zv)
% V nasem pripade by zadny test nemel zamitnout nahodnost.
% Pokud by v realnem pripade k zamitnuti nekterymi testy doslo,
% znacilo by to, ze v rezidualnim sumu existuji korelace a tedy
% je nutno prejit pravdepodobne k modelu ARMA(p,q) s vyssim radem q.
% Protoze MA(q) = ARMA(0,q), mohli bychom v nasem pripade uvazovat
% o modelu napr. MA(q+1).
pause;
%-----------------------------------------------------------------------------
% Pro vypocet rezidualniho sumu z nabizi IDENT funkci 'pe', resp. 'resid',
% ktera navic zobrazi graf autokorelacni funkce pro z, jejiz prubeh by
% mel indikovat bily sum.
% Provedte help pe a help resid.
keyboard;
%-----------------------------------------------------------------------------
  z1 = pe(x,th);    % vypocet residualniho sumu
  acf(z1,figure);   % graf autokorelacni funkce pomoci 'acf'
pause;
%-----------------------------------------------------------------------------
% podobne, ale trochu jinak pro zv:
  ev1 = pe(xv,th);
  figure; zv1 = resid(xv,th);
  max(abs(ev1-zv1))
% Tedy skutecne 'pe' a 'resid' pocitaji totez.
pause;
  zvro1 = acf(zv1);
  figure; plot(zvro1(1:25));
% Porovnejte posledni dva obrazky - prubehy jsou stejne, i kdyz spocteny
% ruznymi funkcemi.
pause;
%-----------------------------------------------------------------------------
% Vidime, ze rovnez grafy autokorelacnich funkci potvrzuji bily sum pro
% z1 a zv1.
% Az na drobne odchylky (dano jinymi pocatecnimi podminkami pro z)
% se residualni sumy spoctene pomoci 'armasim' a 'resid' shoduji
%
  figure; zoom on; plot([z,z1])      % Pouzijte zoom (viz help zoom)
pause;
  figure; zoom on; plot([zv,zv1])
pause;
% Pro kontrolu jeste porovnejme vsechny odhady sigma2
  [v,vz,vzv]
% s presnou hodnotou
  sigma2
pause;
%-----------------------------------------------------------------------------
close all;
% c) na zaver overime vizualne z grafu soulad pozorovanych hodnot
%    s hodnotami predikovanymi z q hodnot pozorovanych o 1 krok zpet
%    pri vyuziti odhadnuteho modelu.
%    Tuto predikci jiz drive spocetla funkce 'armasim' v 'xs' a v 'xvs'.
%    V IDENT predikce pocita funkce 'predict'.
%    Funkce 'compare' navic zobrazi do grafu ruznymi carami predikci
%    soucasne s pozorovanymi hodnotami.
%    Obe funkce predepsanim tretiho parametru umoznuji i predikci
%    z p hodnot o predepsany pocet kroku k>=1 zpet.
%    Provedte help predict a help compare.
keyboard;
%-----------------------------------------------------------------------------
% Spocteme nejprve predikce a porovname je s xs a xvs:
  xs1  = predict(x, th);
  xvs1 = predict(xv,th);
  max(abs(xs-xs1))
  max(abs(xvs-xvs1))
  figure; zoom on; plot([xs,xs1]);   % Pouzijte zoom (viz help zoom)
  figure; zoom on; plot([xvs,xvs1]);
% Shoda je takrka uplna.
pause;
%-----------------------------------------------------------------------------
% Podobne postupujeme s 'compare', ktery navic zobrazi srovnavaci graf:
  figure; zoom on;
  xs2  = compare(x, th,1);
% Detaily lepe uvidite s vyuzitim zoomovani pomoci mysi (viz help zoom).
pause;
  figure; zoom on;
  xvs2 = compare(xv,th,1);
  max(abs(xs1-xs2))
  max(abs(xvs1-xvs2))
% Vsimnete si, ze graf predikce v tomto pripade nevykazuje tvarove posunuti
% jako u AR modelu v 'crcv12'.
% Je to dano prave tim, ze v modelu zcela chybi autoregresni cast,
% ktera nese zbytky deterministicke casti minule informace, ktera zpravidla
% prevysuje cast stochastickou.
pause;
%-----------------------------------------------------------------------------
% Kvalita modelu se nejvice pozna pri predikci na vetsi vzdalenost:
  figure; zoom on;
  compare(x, th,2);
pause;
  figure; zoom on;
  compare(x, th,3);
% I v techto pripadech je shoda tvaru do jiste miry zachovana, takze
% celkem muzeme nalezeny model MA(2) definitivne prijmout za verohodny.
pause;
%-----------------------------------------------------------------------------
%==================================KONEC======================================
echo off
