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

% (C) 1999, V. Vesely, Masaryk University of Brno (Czech Republic)

format compact;
echo on
%========================================================================
%         BOX-JENKINSOVY ARMA MODELY PRO STACIONARNI CASOVE RADY
% Systematicky postup identifikace, odhadu parametru a overeni modelu MA
%                   (vyzaduje NAG FOUNDATION TOOLBOX)
%========================================================================
% Bez ujmy na obecnosti se muzeme omezit na invertibilni
% proces MA(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 kauzalnim AR procesem.
%
% Invertibilita znamena, 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
% 'crnag12'. Pokud vzniknou problemy s urcenim radu, tj. kriterialni ztratove
% funkce (napr. FPE, 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(q) nebo
% smiseneho modelu ARMA(p,q), q>0.
pause;
%-----------------------------------------------------------------------------
%
% V dalsim jiz budeme take vyuzivat funkce z NAG FOUNDATION TOOLBOXu (NAG).
%
% Zjisteni pritomnosti NAG FOUNDATION TOOLBOXu:
if ~exist('nagg')==3
   disp('CHYBA: NAG FOUNDATION TOOLBOX NENALEZEN');
   return;
end
% Libovolny gaussovsky proces MA(q) vygenerujeme pomoci procedury 'armagen'.
% Knihovna NAG nabizi funkce 'NAGacf' a 'NAGpacf' jako alternativu k
% funkcim 'acf' a 'pacf' z knihovny 'cas_rady', coz nam umozni porovnat
% vysledky.
% Z knihovny 'cas_rady' dale vyuzijeme take procedury 'armagen' a 'armasim'.
% Provedte help armagen, help armasim, help acf, help NAGacf,
%          help pacf, help NAGpacf (navod k pouziti),
% resp.
%          type armagen, type armasim, type acf, type NAGacf,
%          type pacf, type NAGpacf (algoritmus).
keyboard;
%---------------------------------------------------------
% DEMONSTRACE ODHADU PARAMETRU SIMULOVANEHO PROCESU MA(2)
%=========================================================
% 1. GENEROVANI PROCESU
%    ------------------
  q=2; n=300;
% Rozptyl, resp. smerodatna odchylka generujiciho gaussovskeho bileho sumu
  sigma2=5e-3; 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
  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-mean(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 by acf(h) mela lezet v nulovem konfidencnim pasu a vykazovat
%   zde nahodne chovani.
%   Podle ocekavani tedy usuzujeme na MA proces radu q=2.
pause;
%-----------------------------------------------------------------------------
% A nyni alternativne s vyuzitim procedury 'NAGacf':
  nq = length(xro)-1;
  [NAGxro,NAGxv,NAGmi,stat] = NAGacf(x,nq);
% vidime, ze ziskane vysledky jsou shodne v ramci strojove presnosti:
  max(abs(xro(2:nq+1)-NAGxro))
  mi-NAGmi
pause;
%-----------------------------------------------------------------------------
% Navic statistika 'stat' ma priblizne rozdeleni 'chi2(nq)' a umoznuje testovat
% hypotezu o nulovosti NAGxro (tj. x je bily sum) s rizikem zamitnuti alpha:
  alpha = 0.05;
  kvantil=chi2inv(1-alpha,nq), stat
% hypotezu o nulovosti NAGxro bychom meli zamitnout s rizikem alpha:
  stat>kvantil
pause;
%-----------------------------------------------------------------------------
% b) ODHAD PARCIALNI AUTOKORELACNI FUNKCE A AUTOREGRESNICH KOEFICIENTU phi
%    (procedura pacf)
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  [xpacf,phi]=pacf(xacf,n,figure);
% Obalka parcialni autokorelacni funkce by mela klesat exponencialne k nule,
% pricemz v konfidencnim pasu nulovosti by nemela mit nahodny charakter.
% Proto predbezne usuzujeme, ze se asi nejedna o AR proces konecneho radu.
pause;
%-----------------------------------------------------------------------------
% Porovname vysledky jeste s vystupy procedury 'NAGpacf':
  [NAGxpacf,NAGphi,vr]=NAGpacf(NAGxro,nq);
  max(abs(NAGxpacf-xpacf(2:nq+1)))
  max(abs(NAGphi-phi))
pause;
%-----------------------------------------------------------------------------
close all
% c) POKUS O ODHAD RADU p PRO MODEL AR(p)
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Odhad v(p) z useku x a vv(p) z useku xv pro kazdy rad p=1,...,pmax a
% nalezeni optimalnich radu pro nektera kriteria spocteme pomoci 'NAGpqest'.
% v(p)  ... odhad rozptylu spocteny z tehoz useku x uziteho pro odhad modelu
% vv(p) ... odhad rozptylu spocteny z nezavisleho useku xv
% v(p)=NaN a vv(p)=NaN, pokud nemohly byt ziskany spolehlive odhady.
% To zpravidla indikuje preparametrizovani (prilis velky rad p).
%
% Provedte help NAGpqest (navod k pouziti), resp.
%          type NAGpqest (algoritmus).
%
keyboard;
%-----------------------------------------------------------------------------
  pmax = 10;
  [pqo,v,vv,ke,pq,rp] = NAGpqest(AR(pmax),x,xv,figure);
  pqo(1,:)   % optimalni rady p spoctene porade z FPE,AIC,MDL a vv
  pq(:,1:ke) % rady p,q, z nichz mohl byt uspesne proveden optimalni vyber
             % v(ke+1) = vv(ke+1) = NaN, pokud ke<pmax.
pause;
  rp  % pocty parametru pripadajici na jedno pozorovani:
      % umoznuje spocist jine kriterium, napr. v.*(1+c*rp) pro vhodne c.
  v10 = v.*(1+10*rp)  % pro c=10.
  figure; plot(v10);
  title('Silne modifikujici kriterium (c=10)');
  xlabel('rad p');
  [v10min,kmin] = min(v10)   % kmin=po ... optimalni rad p dle kriteria v10

pause;
%-----------------------------------------------------------------------------
% Body zlomu vetsiny ztratovych funkci, resp. minima kriterialnich funkci
% by mely byt nevyrazna a davat spise vyssi navzajem se lisici hodnoty.
% Proto se na zaklade prubehu 'acf a 'pacf' ve 2a)b) priklonime
% k neadkvatnosti AR modelu a pokusime se najit spise optimalni MA model.
pause;
%-----------------------------------------------------------------------------
close all
% d) ODHAD RADU q PROCESU MA(q), q>0
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Odhad v(q) z useku x a vv(q) z useku xv pro kazdy rad q=1,...,qmax a
% nalezeni optimalnich radu pro nektera kriteria
% spocteme opet pomoci 'NAGpqest'.
% v(q)  ... odhad rozptylu spocteny z tehoz useku x uziteho pro odhad modelu
% vv(q) ... odhad rozptylu spocteny z nezavisleho useku xv
% v(q)=NaN a vv(q)=NaN, pokud nemohly byt ziskany spolehlive odhady (errc=7).
% To zpravidla indikuje preparametrizovani (prilis velky rad q).
  qmax = 8;
  [pqo,v,vv,ke,pq,rp] = NAGpqest(MA(qmax),x,xv,figure);
  pqo(2,:)   % optimalni rady q spoctene porade z FPE,AIC,MDL a vv
  pq(:,1:ke) % rady p,q, z nichz mohl byt uspesne proveden optimalni vyber
             % v(ke+1) = vv(ke+1) = NaN, pokud ke<pmax.
pause;
  rp  % pocty parametru pripadajici na jedno pozorovani:
      % umoznuje spocist jine kriterium, napr. v.*(1+c*rp) pro vhodne c.
  v10 = v.*(1+10*rp)  % pro c=10.
  figure; plot(v10);
  title('Silne modifikujici kriterium (c=10)');
  xlabel('rad q');
  [v10min,kmin] = min(v10)   % kmin=qo ... optimalni rad dle kriteria v10
%
% Grafy i minima kriterialnich funkci by mely dobre urcovat bod zlomu a tedy
% rad q=2.
pause;
%-----------------------------------------------------------------------------
close all
% 3. ODHAD PARAMETRU PROCESU IDENTIFIKOVANEHO JAKO MA(2)
%    ---------------------------------------------------
% Pro odhad parametru nabizi NAG funkce:
% NAGpest ... predbezny odhad parametru
% NAGest  ... finalni odhad parametru minimalizaci souctu ctvercu odchylek,
%             ktery je obecne nelinearni funkci neznamych parametru.
% Provedte help NAGpest a help NAGest.
keyboard;
%-----------------------------------------------------------------------------
% a) Predbezny odhad pomoci 'NAGpest':
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  q  = 2;
  mr = MA(q)      % vektor radu pro SARIMA model redukovany na MA(q).
  [par0,v0,isf] = NAGpest(mr,NAGxro,NAGxv);
% Predbezne odhady
  theta0 = par0{2}  % parametru theta
  v0                % rozptylu sigma2
  isf(2)            % hodnota 1 by mela indikovat uspesnost odhadu
pause;
%-----------------------------------------------------------------------------
% Autokorelacni funkci nemusime pocitat predem, pokud 'NAGpest' zavolame
% takto:
  [par0,v0,isf] = NAGpest(mr,x);
% Pri tomto volani zaridi vse procedura NAGpest sama a obdrzene vysledky
% jsou opet stejne
  theta0 = par0{2}
  v0
  isf(2)
pause;
%-----------------------------------------------------------------------------
% b) Finalni odhad pomoci 'NAGest':
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  [par1,par1S,v1,est1INFO] = NAGest(mr,x,par0);
% Takto jsme spocetli finalni odhady pomoci metody nejmensich ctvercu
% z NAG procedury 'G13AFF':
  theta1  = par1{2}  % parametru theta
  theta1S = par1S{2} % jejich smerodatnych odchylek
  v1                 % rozptylu sigma2
% Z theta1S spocteme 95% konfidencni polosirky s vyuzitim
% normalni aproximace chyb odhadu parametru.
% Obecne: theta1E=norminv(1-alpha/2)*theta1S, kde 'norminv' je kvantilova funkce
%         normalniho rozlozeni a alpha riziko.
% Odtud:
  theta1E = 1.96*theta1S % 1.96 = norminv(1-0.05/2) odpovida riziku 5%
pause;
%-----------------------------------------------------------------------------
% Prohledneme jeste obsah struktury estINFO, ktera obsahuje dodatecne
% informace:
  est1INFO
% kde je:
         % xdmi        ... odhad stredni hodnoty
         % xdmiS       ... jeji smerodatna odchylka
  est1INFO.corrmx   %  ... korelacni matice odhadu vsech parametru
                    %      (celkem tri parametry vcetne 'xdmi')
  est1INFO.isf(2)   %  ... hodnota 1 by mela indikovat uspesnost odhadu
         % stset       ... stavova mnozina pro dalsi zpracovani
         % z           ... rezidua (mela by reprezentovat pozorovani
         %                          generujiciho bileho sumu Z(t))
         % rss         ... rezidualni soucet ctvercu
         % ndf         ... pocet stupnu volnosti
         % itc         ... skutecny pocet iteraci provedenych pri minimalizaci
         %                 souctu ctvercu odchylek.
pause;
%-----------------------------------------------------------------------------
% Do NAGest lze pridat dalsi vstupni parametry, napriklad vlastni
% pocatecni odhad 'xmi', resp. potlacit odhadovani 'xmi' zadanim xdmiEST = 0:
  xmi0 = mean(x);
  xmiEST = 0;
  [par2,par2S,v2,est2INFO] = NAGest(mr,x,par0,[],xmi0,xmiEST);
% Finalni odhady
  theta2  = par2{2}       % parametru theta
  theta2E = 1.96*par2S{2} % jejich konfidencnich polosirek
  v2                      % rozptylu sigma2
  est2INFO
pause;
%-----------------------------------------------------------------------------
% Definitivni odhad spocteme metodou maximalni verohodnosti z NAG procedury
% 'G13BEF':
  [par3,par3S,v3,est3INFO] = NAGest(mr,x,par0,'MLE');
% Finalni odhady
  theta3  = par3{2}       % parametru theta
  theta3E = 1.96*par3S{2} % jejich konfidencnich polosirek
  v3                      % rozptylu sigma2
  est3INFO
  xmi3 = est3INFO.xdmi;   % odhad stredni hodnoty
pause;
%-----------------------------------------------------------------------------
% 4. OVERENI MODELU ODHADNUTEHO JAKO MA(2)
%    -------------------------------------
% a) Nejprve zkontrolujeme, zda hodnoty koncovych parametru jsou
%    vyznamne (jejich konfidencni intervaly neobsahuji nulu):
  theta3 = theta3'; theta3E = theta3E';
  [theta3-theta3E,theta3+theta3E]
% Pokud by napriklad [theta3(q)-theta3E(q),theta3(q)+theta3E(q)]
% obsahoval nulu, znamenalo by to, ze theta3(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.
%    Uzijeme pro verifikaci nejprve usek dat 'x', z nehoz byl nalezen odhad
%    parametru modelu:
  [xs,z,vz] = armasim(x-xmi3,[],theta3);
  xs = xs + xmi3;                    % korekce o stredni hodnotu
% Porovname s estINFO.z:
  max(abs(z-est3INFO.z))
  plot([z,est3INFO.z]); zoom on;
% odhady se lisi nejvice na pocatku v dusledku rozdilnych pocatecnich
% podminek pro Z(t). Totiz 'armasim' dosazuje nulove, zatimco 'NAGest'
% spocte Z(t) z prislusne invertibilni reprezentace.
%
% Spocteme 'z' jinak uzitim standardni funkce MATLABu 'filter':
  zz = filter(1,[1;theta3],x-xmi3); xxs = x-zz; vzz = (norm(z).^2)./n;
  max(abs(z-zz)), max(abs(xs-xxs)), vz-vzz
% Provedte help filter.
keyboard;
  zoom off;
%-----------------------------------------------------------------------------
% Nyni pro verifikaci uzijeme nezavisly usek 'xv':
  [xvs,zv,vzv] = armasim(xv-xmi3,[],theta3,x-xmi3);
  xvs = xvs + xmi3;
  zzv = filter(1,[1;theta3],xv-xmi3);
  xxvs = xv-zzv; vzzv = (norm(zzv).^2)./n;
  max(abs(zv-zzv)), max(abs(xvs-xxvs)), vzv-vzzv
% Jeste jednou presneji s doplnenim pocatecnich podminek pro 'z':
  [xvs,zv,vzv] = armasim(xv-xmi3,[],theta3,x-xmi3,z);
  xvs = xvs + xmi3;
% Vysledky se nepatrne zmenily:
  max(abs(zv-zzv)), vzv-vzzv
  plot([xvs,xxvs]); zoom on;
pause;
% Proverime vlastnost IID rezidulaniho sumu napr. pomoci 'trand'
  z = est3INFO.z;
  xs = x - z;
  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 p nebo q.
% Protoze MA(q) = ARMA(0,q), museli bychom v nasem pripade uvazovat
% o modelu napr. MA(q+1) nebo o smisenem modelu pro p,q>0.
pause;
%-----------------------------------------------------------------------------
% Spocteme jeste autokorelacni funkci rezidualniho sumu z, resp. zv pomoci
% 'acf' nebo 'NAGacf' a proverime, zda se jedna o bily sum:
  acf(z ,figure);
  acf(zv,figure);
% V obou pripadech by vsechny hodnoty pro h>0 mely byt v 95% konfidencnim
% pasu nulovosti, takze hypotezu o bilem sumu nezamitneme.
pause;
%-----------------------------------------------------------------------------
% Podobne uzitim 'NAGacf':
  nqz = round(sqrt(n));
  [zro,ans,ans,stat] = NAGacf(z,nqz);
  stat, kvantil=chi2inv(0.95,nqz-q)
% tato relace by nemela platit, pokud nezamitame hypotezu o bilem sumu.
  stat>kvantil
pause;
close all;
%-----------------------------------------------------------------------------
% Pro vypocet rezidualniho sumu z nabizi NAG funkci 'NAGcheck'
% ktera navic spocte autokorelacni funkci pro 'z' vcetne presnejsi
% portmanteau statistiky pro testovani, zda z je bily sum.
% Provedte help 'NAGcheck'.
keyboard;
%-----------------------------------------------------------------------------
  [zro3,zro3cm,stat3,idf3,siglev3] = NAGcheck(mr,par3,z);
  max(abs(zro3-zro))
  idf3, siglev3
  stat3, kvantil=chi2inv(0.95,idf3)
% tato relace by nemela platit, pokud nezamitame hypotezu o bilem sumu.
  stat3>kvantil
% nebo alternativne:
% tato relace by opet nemela platit, pokud nezamitame hypotezu o bilem sumu.
  0.05>siglev3
% spocteme jeste polosirky intervalu spolehlivosti pro 'zro3' pri normalni
% aproximaci rozlozeni odhadu:
  zro3E = 1.96*diag(zro3cm);
  errorbar(1:length(zro3),zro3,zro3E,zro3E,'o');
  title('Rezidualni autokorelacni funkce');
pause;
%-----------------------------------------------------------------------------
% c) na zaver overime vizualne z grafu soulad pozorovanych hodnot
%    s hodnotami predikovanymi z p hodnot pozorovanych o 1 krok zpet
%    pri vyuziti odhadnuteho modelu.
%    Tuto predikci jiz drive spocetla funkce 'armasim' v 'xs' a v 'xvs'.
%    V NAG predikce pocita funkce 'NAGpred'.
%    Provedte help 'NAGpred'.
keyboard;
%-----------------------------------------------------------------------------
% Porovname 1-krokovou predikci 'xs' se skutecnym pozorovnim puvodnich dat x:
   plot(x,'y'); hold on; plot(xs,'m'); hold off;
   title('Pozorovana data x (zluta), predikce (fialova)');
   zoom on
% Detaily lepe uvidite s vyuzitim zoomovani pomoci mysi (viz help zoom).
% Vsimnete si, ze graf predikce v tomto pripade nevykazuje tvarove posunuti
% jako u AR modelu v 'crnag12'.
% 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;
%-----------------------------------------------------------------------------
% Porovname 1-krokovou predikci 'xvs' se skutecnym pozorovnim nasledujiciho
% nezavisleho useku xv puvodnich dat:
   plot(xv,'y'); hold on; plot(xvs,'m'); hold off;
   title('Nezavisly usek pozorovanych dat (zluta), predikce (fialova)');
   zoom on
% Detaily opet lepe uvidite s vyuzitim zoomovani pomoci mysi (viz help zoom).
pause;
%-----------------------------------------------------------------------------
  zoom off
% Spocteme nyni predikci o 10 kroku pomoci 'NAGpred' a porovname ji s xv(1:10):
  npv = 10;   % delka predikce
  [xp3,xp3S] = NAGpred(mr,par3,npv,est3INFO);
% polosirka intervalu spolehlivosti pri normalni aproximaci chyb:
  xp3E = 1.96*xp3S;
  plot(xv(1:npv),'y'); hold on;
  plot(xp3,'m'); plot([xp3-xp3E,xp3+xp3E],'m--'); hold off;
  title('Skutecny usek pozorovanych dat (zluta), jeho predikce (fialova)');
pause;
%-----------------------------------------------------------------------------
% Kvalita odhadu modelu se nejvice pozna pri predikci na vetsi vzdalenost.
% Porovname idealni predikci s odhadnutou:
  estINFO = est3INFO; estINFO.xdmi = 0;
  xp = NAGpred(mr,{[],THETA,[],[]},npv,estINFO);  % IDEALNI PREDIKCE
  t=(n+1:n+npv).';
  figure;plot(t,xp,'w',t,xp3,'y--')
  title('Predikce: idealni (-), odhadnuta (--)');
%==================================KONEC======================================
echo off
