%========================================================================
%         BOX-JENKINSOVY ARMA MODELY PRO STACIONARNI CASOVE RADY
% Systematicky postup identifikace, odhadu parametru a overeni smiseneho
%                      modelu SARIMA(p,d,q,P,D,Q,s)
%                   (vyzaduje NAG FOUNDATION TOOLBOX)
%========================================================================
% ANALYZA SOUBORU deaths.dat:
% Mesicni pocty umrti pri automobilovych nehodach v USA
% x=absolutni pocet umrti
% t=1,2,...,72: 1973 az 1978
%=================================================================
% (C) 1999, V. Vesely, Masaryk University of Brno (Czech Republic)
%=================================================================
echo on
format compact
% Zjisteni pritomnosti NAG FOUNDATION TOOLBOXu:
if ~exist('nagg')==3
   disp('CHYBA: NAG FOUNDATION TOOLBOX NENALEZEN');
   return;
end
%========================================================================
% ANALYZA SOUBORU deaths.dat:
% Mesicni pocty umrti pri automobilovych nehodach v USA
% x=absolutni pocet umrti
% t=1,2,...,72: 1973 az 1978
%========================================================================
% 1. NACTENI DAT
%    -----------
  us_deaths=[crpath,'/brockwel.dat/deaths.dat'];
  [x,xinf]=getdata(us_deaths);   % nacteni a vykresleni
  n=length(x)
% Proverime homoskedasticitu:
  [lambdaE,y] = powtr(x,4,figure);
% Transformace dat neni nutna.
pause;
%========================================================================
% 2. IDENTIFIKACE PROCESU
%    --------------------
% a) ODHAD PERIODY s DOMINANTNI SEZONNI SLOZKY
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Data naznacuji existenci vyznmane periodicke (sezonni) slozky.
% Dominantni sezonni periodu s stanovime pomoci autokorelacni
% funkce (acf) a periodogramu (tperiod):
  acf(x,figure);
  [I,kF,kS] = tperiod(x,figure);
  kF, kS  % Vyznamne harmonicke komponenty dle Fisherova a Siegelova testu:
% Vidime, ze nejsilneji je zastoupena 6-ta harmonicka a jeji celociselne
% nasobky 2*6=12, 3*6=18 a 4*6=24, 5*6=30, jimz odpovidaji mesicni periody:
  n./kF(rem(kF,6)==0)
% Dominantni periodou je tedy
  s = n./kF(1)
pause;
close all
%---------------------------------------------------------------
% b) ODHAD RADU DIFERENCOVANI, STREDNI HODNOTY A AUTOKOVARIANCNI
%    FUNKCE DIFERENCOVANEHO PROCESU (NAGdiff,acf)
%---------------------------------------------------------------
% Rady d,D urcime vyberem diferencovane rady s nejmensim rozptylem,
% pricemz z prubehu diferencovane rady a jeji 'acf' budeme vizualne
% overovat stacionaritu jejiho prubehu.
  std(x).^2
  d=1, D=0
  xd=NAGdiff([d,D,s],x); acf(xd,figure); std(xd)^2
% Podle ocekavani neni stacionarni, nebot v 'acf' se vyskytuji
% silne korelace s priodou s=12.
pause;
  d=0, D=1
  xd=NAGdiff([d,D,s],x); acf(xd,figure); std(xd)^2
% Zjevne rada neni stacionarni (vyskyt trendu).
% Linearni pokles 'acf' pro h=1,...,s-1 indikuje nutnost zvyseni radu d:
pause;
  d=1, D=1
  xd=NAGdiff([d,D,s],x); acf(xd,figure); std(xd)^2
% Rada jiz vypada jako stacionarni (rychly pokles 'acf'),
% pricemz rozptyl je mensi nez v predchozich pripadech.
% Vyzkousime jeste dale zvysit rad diferencovani.
pause;
  d=1, D=2
  xd=NAGdiff([d,D,s],x); acf(xd,figure); std(xd)^2
  d=2, D=1
  xd=NAGdiff([d,D,s],x); acf(xd,figure); std(xd)^2
% V obou pripadech jiz dalsim diferencovanim vzrusta rozptyl.
pause;
% Zvolime proto rady:
  d=1; D=1;
  close all
%---------------------------------------------------------------
% c) PREDBEZNY ODHAD NESEZONNICH RADU p,q a SEZONNICH RADU P,Q
%    analyzou grafu autokorelacni (acf) a parcialni autokorelacni
%    funkce (pacf).
%---------------------------------------------------------------
  xd = NAGdiff([d,D,s],x);
  nd = length(xd)
  xdro = acf(xd,figure);
% Podrobnejsi rozbor korelaci z grafu 'acf' ukazuje:
% Vyznamne sezonni korelace:
%   xdro(h), priblizne pro h=k*s je vyznamne nenulova jen pro k=1,
%   pro k>1 vykazuje prudky pokles => asi P=0, Q=1.
% Vyznamne nesezonni (sousedni) korelace:
%   xdro(h), h=1,...,s-1, je vyznamne nenulova jen pro h=1,
%   pro h>1 vykazuje prudky pokles => asi p=0, q=1.
pause;
  [xdpacf,phi]=pacf(xdro,nd,figure);
% Podrobnejsi rozbor parcialnich korelaci z grafu 'pacf' ukazuje:
% Vyznamne sezonni parcialni korelace:
%   xdpacf(h), priblizne pro h=ks klesa pro k=1,2,... velmi pozvolna,
%   tj. sezonne se asi nejedna o AR proces => potvrzuje P=0, Q=1.
% Vyznamne nesezonni (sousedni) parcialni korelace:
%   xdpacf(h), h=1,...,s-1, je vyznamne nenulova jen pro h=1,2, pricemz
%   xdpacf(2) je na hranici vyznamnosti.
%   To svedci spise pro pomaly pokles => opet potvrzeni p=0, q=1.
pause;
%---------------------------------------------------------------
% d) DEFINITIVNI ODHAD NESEZONNICH RADU p,q a SEZONNICH RADU P,Q
%    uzitim 'NAGpqest'.
%---------------------------------------------------------------
% Vykreslime nejprve prubeh sezonnich komponent:
  na = ceil(nd/s)*s - nd  % pocet hodnot NaN k doplneni na nasobek period
  np = (nd+na)/s          % pocet period po doplneni
  Xd = reshape([ones(na,1)*NaN;xd],s,np);
% prumerna diferencovana sezonni slozka:
  xdsz = [mean(Xd(na+1:s,1:(na>0))),mean(Xd(:,(na>0)+1:np))].';
% Vykreslime prubehy vsech sezonnich komponent spolu s jejich prumerem (tucne)
  figure; plot(Xd.'); hold on; plot(xdsz,'w','LineWidth',2); hold off
pause;
% Sezonni slozky maji jen 5 hodnot, takze je mozno uzit nejvyse 3 parametry,
% tj. Pmax+Qmax=3:
  Pmax=1; Qmax=2;
  [PQo,v,vv,ke,pq,rp] = NAGpqest(ARMA(Pmax,Qmax),xdsz,[],figure);
  PQo      % Odhad optimalnich radu P,Q z 'xdsz' dle FPE,AIC a MDL
pause;
% Spocteme jeste optima z jednotlivych sezonnich slozek:
  for j=2:s
      PQj = NAGpqest(ARMA(Pmax,Qmax),Xd(j,:));
      PQo = [PQo,PQj];
echo off
  end
echo on
% Procenta uspesnosti pro kombimace P,Q:
  100*sum(PQo(1,:)==0 & PQo(2,:)==1)/(3*s) % pro P=0, Q=1
  100*sum(PQo(1,:)==0 & PQo(2,:)==2)/(3*s) % pro P=0, Q=2
  100*sum(PQo(1,:)==1 & PQo(2,:)==1)/(3*s) % pro P=1, Q=1
  100*sum(PQo(1,:)==1 & PQo(2,:)==2)/(3*s) % pro P=1, Q=2
% Nejvetsi uspesnost ma P=0 a Q=1 v souladu s predbeznym odhadem.
  P=0; Q=1;
pause;
  close all
% Zbyva urcit nesezonni rady p,q:
% Nejprve zkusime AR model:
  pmax=8;
  [pqo,v,vv,ke,pq,rp] = NAGpqest(SARIMA(pmax,d,0,P,D,Q,s),x,[],figure);
  pqo        % optimalni rady p,q spoctene porade z FPE,AIC a MDL
  pq(:,1:ke) % rady p,q, z nichz mohl byt uspesne proveden optimalni vyber
             % v(ke+1) = NaN, pokud ke<pmax.
pause;
% Nyni proverime MA modely a smisene ARMA modely:
  pmax=4; qmax=4;
  [pqo,v,vv,ke,pq,rp] = NAGpqest(SARIMA(pmax,d,qmax,P,D,Q,s),x,[],figure);
  pqo        % optimalni rady p,q spoctene porade z FPE,AIC a MDL
  pq(:,1:ke) % rady p,q, z nichz mohl byt uspesne proveden optimalni vyber
             % v(ke+1) = NaN, pokud ke<(pmax+1)*qmax.
% Nabizi se volba p=2,q=0 nebo p=1,q=2.
% Vyzkousime obe varianty.
pause;
  close all
%=========================================================================
% 3. ODHAD PARAMETRU PROCESU MODELU SARIMA(p,d,q,P,D,Q,s) VCETNE VERIFIKACE
%    ----------------------------------------------------------------------
% Pro odhad parametru a jeho verifikaci nabizi NAG funkce:
% NAGpest ... predbezny odhad parametru
% NAGest  ... finalni odhad parametru minimalizaci souctu ctvercu odchylek,
%             ktery je obecne nelinearni funkci neznamych parametru.
% NAGcheck... verifikace, zda rezidualni sum je bily
%-----------------------------------------------------------------------------
% a) Predbezny odhad pomoci 'NAGpest':
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  p1  = 2; q1 = 0;     % varianta 1
  mr1 = SARIMA(p1,d,q1,P,D,Q,s) % prislusny vektor radu pro SARIMA model
  p2  = 1; q2 = 2;     % varianta 2
  mr2 = SARIMA(p2,d,q2,P,D,Q,s) % prislusny vektor radu pro SARIMA model
  [NAGxdro,NAGxdv] = NAGacf(xd,nd-1);
  [par1_0,v,isf] = NAGpest(mr1,NAGxdro,NAGxdv);
% Predbezne odhady varianty 1
  par1_0{:}   % porade parametru phi0, theta0, PHI0, THETA0
  v           % rozptylu sigma2
  isf         % hodnota 1 by mela indikovat uspesnost odhadu
pause;
  [par2_0,v,isf] = NAGpest(mr2,NAGxdro,NAGxdv);
% Predbezne odhady varianty 2
  par2_0{:}   % porade parametru phi0, theta0, PHI0, THETA0
  v           % rozptyl sigma2
  isf         % hodnota 1 by mela indikovat uspesnost odhadu
pause;
%-----------------------------------------------------------------------------
% b1) Finalni odhad pomoci 'NAGest' pro variantu 1:
%     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  [par1,par1S,v1,est1INFO] = NAGest(mr1,x,par1_0);
% Takto jsme pro variantu 1 spocetli finalni odhady pomoci metody
% nejmensich ctvercu z NAG procedury 'G13AFF':
  par1{:}
  par1S{:}
  v1                 % rozptyl sigma2
  est1INFO
pause;
% Definitivni odhad pro variantu 1 spocteme metodou maximalni verohodnosti
% z NAG procedury 'G13BEF':
  [par1,par1S,v1,est1INFO] = NAGest(mr1,x,par1_0,'MLE');
  par1{:}
  abs(par1{1})./(1.96*par1S{1})  % => phi(2) je nevyznamny
  abs(par1{4})./(1.96*par1S{4})
  v1                 % rozptyl sigma2
  est1INFO
pause;
% Snizime rad p:
  mr1(1)=mr1(1)-1
  [par1_0,v,isf] = NAGpest(mr1,NAGxdro,NAGxdv);
  isf
  [par1,par1S,v1,est1INFO] = NAGest(mr1,x,par1_0,'MLE');
  par1{:}
  abs(par1{1})./(1.96*par1S{1})
  abs(par1{4})./(1.96*par1S{4})
  v1                 % rozptyl sigma2
  est1INFO
% Vsechny parametry jsou vyznamne
pause;
% Proverime, zda je rezidualni proces bilym sumem
  z1 = est1INFO.z;
  trand(z1)
  acf(z1,figure);
  [zro1,zro1cm,stat1,idf1,siglev1] = NAGcheck(mr1,par1,z1);
  idf1, siglev1
% Jen tesne nezamitame hypotezu o bilem sumu.
  0.05>siglev1
% Tedy varianta 1 vede prozatim na
  p1 = mr1(1), q1 = mr1(3)
pause;
%-----------------------------------------------------------------------------
% b2) Finalni odhad pomoci 'NAGest' pro variantu 2:
%     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  [par2,par2S,v2,est2INFO] = NAGest(mr2,x,par2_0);
% Takto jsme pro variantu 2 spocetli finalni odhady pomoci metody
% nejmensich ctvercu z NAG procedury 'G13AFF':
  par2{:}
  par2S{:}
  v2                 % rozptyl sigma2
  est2INFO
pause;
% Definitivni odhad pro variantu 2 spocteme metodou maximalni verohodnosti
% z NAG procedury 'G13BEF':
  [par2,par2S,v2,est2INFO] = NAGest(mr2,x,par2_0,'MLE');
  par2{:}
  abs(par2{1})./(1.96*par2S{1})
  abs(par2{2})./(1.96*par2S{2})  % => parametr theta(2) je nevyznamny
  abs(par2{4})./(1.96*par2S{4})
  v2                 % rozptyl sigma2
  est2INFO
pause;
% Snizime rad q:
  mr2(3)=mr2(3)-1
  [par2_0,v,isf] = NAGpest(mr2,NAGxdro,NAGxdv);
  isf
  [par2,par2S,v2,est2INFO] = NAGest(mr2,x,par2_0,'MLE');
  par2{:}
  abs(par2{1})./(1.96*par2S{1})
  abs(par2{2})./(1.96*par2S{2})
  abs(par2{4})./(1.96*par2S{4})
% Vsechny parametry jsou vyznamne
  v2                 % rozptyl sigma2
  est2INFO
pause;
% Proverime, zda je rezidualni proces bilym sumem
  z2 = est2INFO.z;
  trand(z2)
  acf(z2,figure);
  [zro2,zro2cm,stat2,idf2,siglev2] = NAGcheck(mr2,par2,z2);
  idf2, siglev1, siglev2
% Tesne zamitame hypotezu o bilem sumu.
  0.05>siglev2
% Pritom hladina vyznamnosti 'siglev1' je jen nepatrne vetsi, nez 'siglev2',
% coz ponekud svedci o uprednostneni varianty 1 oproti 2.
% Oba modely se ale nezdaji byt prukazne potvrzeny.
pause;
% Vyzkousime, jeste variantu 3, jez odpovida predbezne volbe p=0, q=1 ad 2c),
% tj. snizime rad p jeste jednou az na p=0:
  mr3 = mr2;
  mr3(1)=mr3(1)-1
  [par3_0,v,isf] = NAGpest(mr3,NAGxdro,NAGxdv);
  isf
  [par3,par3S,v3,est3INFO] = NAGest(mr3,x,par3_0,'MLE');
  par3{:}
  abs(par3{2})./(1.96*par3S{2})
  abs(par3{4})./(1.96*par3S{4})
% Vsechny parametry jsou vyznamne
  v3                 % rozptyl sigma2
  est3INFO
pause;
% Proverime, zda je rezidualni proces bilym sumem
  z3 = est3INFO.z;
  trand(z3)
  acf(z3,figure);
  [zro3,zro3cm,stat3,idf3,siglev3] = NAGcheck(mr3,par3,z3);
  idf3, siglev1, siglev3
% Nasledujici relace spolehlive neplati, takze nezamitame hypotezu o
% bilem sumu.
  0.05>siglev3
% Hladina vyznamnosti 'siglev3' je totiz vyrazne vetsi, nez 'siglev1'
% u dosud nejlepsi varianty 1. Definitivne proto prijimame variantu 3:
  p3 = mr3(1), q3 = mr3(3)
% Prohledneme jeste korelacni strukturu odhadnutych parametru:
  est3INFO.corrmx
% vsechny odhady jsou mezi sebou temer nezavisle, coz svedci
% o optimalni parametrizaci modelu.
pause;
  close all
%-----------------------------------------------------------------------------
% 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.
%-----------------------------------------------------------------------------
% SARIMA model pretransformujeme nejdrive na ekvivalentni ARMA model
% pomoci 'ARMApar':
  xdmi3 = est3INFO.xdmi;
  [phi3x,theta3x,xc3,phi3d] = ARMApar(mr3,par3,xdmi3);
% Porovname 1-krokovou predikci 'xdp' s diferencovanymi daty 'xd':
  [xdp,zd,vzd] = armasim(xd-xdmi3,phi3d,theta3x);
  xdp=xdp+xdmi3;
  plot(xd,'y'); hold on; plot(xdp,'m'); hold off;
  title('Diferencovana data xd (zluta), predikce (fialova)');
  zoom on
% Detaily lepe uvidite s vyuzitim zoomovani pomoci mysi (viz help zoom).
pause;
%-----------------------------------------------------------------------------
% Jeste jednou totez s tim, ze 1-krokovou predikci provedeme v puvodni casove
% rade. Nulove pocatecni podminky nelze u SARIMA procesu predpokladat,
% proto pocatecni usek rady o delce px pouzijeme pro tento ucel.
  px = length(phi3x);
  [xp,z,vz] = armasim(x(px+1:n),phi3x,theta3x,x(1:px));
  xp=xp+xc3;
% Porovname 1-krokovou predikci 'xp' s pozorovanymi daty 'x':
  figure; plot(x(px+1:n),'y'); hold on; plot(xp,'m'); hold off;
  title('Pozorovana data x (zluta), predikce (fialova)');
  zoom on
% Detaily lepe uvidite s vyuzitim zoomovani pomoci mysi (viz help zoom).
pause;
%=============================================================
% 4. PREDIKCE UZITIM ODHADNUTEHO MODELU SARIMA(0,1,1,0,1,1,12)
%    ---------------------------------------------------------
% Uzijeme funkci 'NAGpred'.
%-----------------------------------------------------------------------------
% Spocteme predikci o celou sezonu s=12:
  npv = s;   % delka predikce
  [xp3,xp3S] = NAGpred(mr3,par3,npv,est3INFO);
% polosirka intervalu spolehlivosti pri normalni aproximaci chyb:
  xp3E = 1.96*xp3S;
  [xp3,xp3S,xp3E].' % Tabelace vysledku predikce:
  plot(xp3,'y'); hold on; plot([xp3-xp3E,xp3+xp3E],'r--'); hold off
  title('Predikce');
pause;
% Znazornime jeste predikci spolu s puvodnimi daty:
  t=1:n; tp=n+1:n+npv;
  figure; plot(t,x,'w',[n,tp],[x(n);xp3],'y',tp,[xp3-xp3E,xp3+xp3E],'r--');
  title('Pozorovana data (bila), predikce (zluta)');
%================================== KONEC ==================================
echo off
