%=================================================================
%      BOX-JENKINSOVY ARMA MODELY PRO STACIONARNI CASOVE RADY
%       Predbezna identifikace modelu dle prubehu acf a pacf
%=================================================================

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

format compact;
echo on
%=================================================================
%      BOX-JENKINSOVY ARMA MODELY PRO STACIONARNI CASOVE RADY
%       Predbezna identifikace modelu dle prubehu acf a pacf
%=================================================================
% Libovolny gaussovsky proces ARMA(p,q) vygenerujeme pomoci
% procedury 'armagen'.
% Jeho identifikaci pak provadime pomoci odhadu autokovariancni funkce
% (procedura 'acf') a odhadu parcialni autokorelacni funkce (procedura 'pacf').
% Provedte help armagen, help acf, help pacf (navod k pouziti), resp.
%          type armagen, type acf, type pacf (algoritmus)
keyboard;
%-----------------------------------------------------------
% 1. DEMONSTRACE ODHADU PARAMETRU SIMULOVANEHO PROCESU AR(2)
%===========================================================
% 1.1 GENEROVANI PROCESU
  p=2; n=1000;
% Rozptyl, resp. smerodatna odchylka generujiciho gaussovskeho bileho sumu
  sigma2=2.25; sigma=sqrt(sigma2);
% Generovani AR(2) procesu s parametry:
  PHI = [0.5,0.2];
  [x,phir,thetar]=armagen(PHI,[],n+10,sigma);
  phir              % Koreny polynomu PHI(z) lezi vne jednotkoveho kruhu,
                    % takze proces je kauzalni
  xx=x(n-p+1:n+10); % Posledni dve + 10 budoucich pozorovanych hodnot
  x=x(1:n)  ;       % Dosud pozorovane hodnoty
pause;
%-----------------------------------------------------------------------------
% 1.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 lezi v intervalu spolehlivosti s polosirkou malou relativne
%   velikosti dat. Proto nezamitame hypotezu o nulove stredni hodnote.
% - Obalka autokorelacni funkce klesa exponencialne k nule, pricemz
%   v konfidencnim pasu nulovosti nema nahodny charakter.
%   Proto usuzujeme, ze se nejedna o MA proces.
pause;
%-----------------------------------------------------------------------------
% b) ODHAD PARCIALNI AUTOKORELACNI FUNKCE A AUTOREGRESNICH KOEFICIENTU phi
%    (procedura pacf)
  [xpacf,phi]=pacf(xacf,n,figure);
% Vidime, ze pro h>2 lezi pacf(h) v nulovem konfidencnim pasu a vykazuje
% zde nahodne chovani.
% Podle ocekavani tedy usuzujeme na AR proces radu p=2.
pause;
%-----------------------------------------------------------------------------
% POZNAMKA: Pokud nepozadujeme vystupni parametr v, lze do pacf dosadit
%           misto autokovariancni funkce xacf autokorelacni funkci xro.
  [xpacf0,phi0]=pacf(xro,n,figure);
  max(abs(xpacf-xpacf0))
  max(abs(phi-phi0))
% Obdrzene vysledky jsou stejne.
pause;
%-----------------------------------------------------------------------------
% c) ODHAD PARAMETRU PROCESU IDENTIFIKOVANEHO JAKO AR(2)
  [xpacf,phi,phiE,v]=pacf(xacf(1:p+1),n,gcf);
  [phi,phiE]  % ZAVERECNY ODHAD PARAMETRU a jejich konfidencni polosirky
  v           % ZAVERECNY ODHAD ROZPTYLU sigma2
% Vidime, ze ziskany odhad phi se blizi idealnim hodnotam PHI, ktere jsou
% zahrnuty v konfidencnim intervalu.
% Podobne v je dobrym odhadem sigma2 = 2.5.
pause;
%-----------------------------------------------------------------------------
% 1.3. POUZITI MODELU PRO PREDIKCI VYHLAZENYCH I POZOROVANYCH HODNOT
% Vyhlazeni modelujeme zadanim nulove smerodatne odchylky bileho sumu:
  xs = [x(n-p+1:n);armagen(PHI     ,[],10,0,x(n-p+1:n))]; % IDEALNI   VYHLAZENI
  xse= [x(n-p+1:n);armagen(phi(1:p),[],10,0,x(n-p+1:n))]; % ODHADNUTE VYHLAZENI
% Do xs, resp. xse zarazujeme i p predchozich hodnot, pomoci nichz se
% predikce spocetla
  t=(n-p+1:n+10).';
  figure;plot(t,xs,'w',t,xse,'y--')
  title('Vyhlazena predikce: idelani (-), odhadnuta (--)');
pause;
%-----------------------------------------------------------------------------
% ODHAD PREDIKCE BUDOUCICH POZOROVANYCH HODNOT, TJ. VCETNE SUMU
% (xx ... ocekavane hodnoty spoctene drive)
  xxe=[x(n-p+1:n);armagen(phi(1:2),[],10,sqrt(v),x)];
  figure;plot(t,xx,'w',t,xxe,'y--')
  title('Odhad predikce vcetne sumu c.1');
pause;
%-----------------------------------------------------------------------------
% JESTE JEDNOU SE ZNOVU GENEROVANYMI BUDOUCIMI HODNOTAMI
  xx = [x(n-p+1:n);armagen(PHI     ,[],10,sigma  ,x)];
  xxe= [x(n-p+1:n);armagen(phi(1:2),[],10,sqrt(v),x)];
  figure;plot(t,xx,'w',t,xxe,'y--')
  title('Odhad predikce vcetne sumu c.2');
% OPAKUJTE PREDIKCE NEKOLIKRAT (diky nahode se nekdy podari lepe, nekdy hure).
% POTOM ZOPAKUJTE CELY POSTUP ad 1 PRO AR-PROCESY JINYCH RADU (1,3, atd.).
% PHI VOLTE TAK, ABY PROCES BYL KAUZALNI.
% (prikaz armagen vypise varovani, pokud tomu tak neni)
keyboard;
%===========================================================
% 2. DEMONSTRACE ODHADU PARAMETRU SIMULOVANEHO PROCESU MA(2)
%===========================================================
% 2.1 GENEROVANI PROCESU
  q=2; n=1000;
% 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,n,sigma);
  thetar,             % Koreny polynomu THETA(z) lezi vne jednotkoveho kruhu,
                      % takze proces je invertibilni
pause;
%-----------------------------------------------------------------------------
% 2.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 lezi v intervalu spolehlivosti s polosirkou malou relativne
%   velikosti dat. Proto 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 usuzujeme, ze se nejedna o AR proces.
pause;
%-----------------------------------------------------------------------------
% c) ODHAD PARAMETRU PROCESU IDENTIFIKOVANEHO JAKO MA(2)
%    Odhad parametru theta a rozptylu v je pro MA-proces komplikovanejsi a
%    vede na system nelinearnich rovnic.
%    Je vhodne pouzit specialni software, napr.:
%    SYSTEM IDENTIFICATION TOOLBOX (viz demo 'crcv12') nebo
%    NAG FOUNDATION TOOLBOX (viz demo 'crnag12')
keyboard; % ZOPAKUJTE VYPOCET ad 2. PRO INVERTIBILNI MA-PROCESY JINYCH RADU
%===========================================
% 3. DEMONSTRACE SMISENEHO PROCESU ARMA(2,2)
%===========================================
% 3.1 GENEROVANI PROCESU
  p=2; q=2; n=1000;
% Rozptyl, resp. smerodatna odchylka generujiciho gaussovskeho bileho sumu
  sigma2=2.25; sigma=sqrt(sigma2);
% Generovani ARMA(2,2) procesu s parametry:
  PHI = [0.5,0.2]; THETA = [-0.4,0.3];
  [x,phir,thetar]=armagen(PHI,THETA,n,sigma);
  thetar, phir
  abs(thetar), abs(phir)
% Koreny polynomu PHI(z) i THETA(z) lezi vne jednotkoveho kruhu,
% takze proces je kauzalni a invertibilni.
pause;
%-----------------------------------------------------------------------------
% 3.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 lezi v intervalu spolehlivosti s polosirkou malou relativne
%   velikosti dat. Proto nezamitame hypotezu o nulove stredni hodnote.
% - Obalka autokorelacni funkce klesa exponencialne k nule, pricemz
%   v konfidencnim pasu nulovosti nema nahodny charakter.
%   Proto usuzujeme, ze se nejedna o MA proces.
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 usuzujeme, ze se nejedna o AR proces.
% Bude se tedy asi jednat o smiseny proces ARMA.
% Identifikace radu p a q je v tomto pripade obtiznejsi a vyzaduje specialni
% software.
pause;
%-----------------------------------------------------------------------------
% c) ODHAD PARAMETRU PROCESU IDENTIFIKOVANEHO JAKO ARMA(?,?)
%    Je opet nutne pouzit specialni software.
%    Pro identifikaci a odhad parametru ARMA modelu a jejich slozitejsich
%    variant (ARMAX, aj.) lze pouzit napr.
%    SYSTEM IDENTIFICATION TOOLBOX (viz demo 'crcv12' az 'crcv14') nebo
%    NAG FOUNDATION TOOLBOX (viz demo 'crnag12' az crnag14')
% ZOPAKUJTE VYPOCET ad 3. PRO KAUZALNI A INVERTIBILNI ARMA-PROCESY JINYCH RADU
keyboard;
%================================================
% 4. DEMONSTRACE CHOVANI NEKAUZALNIHO AR PROCESU
%================================================
% GENEROVANI PROCESU
  p=2; n=1000;
% Rozptyl, resp. smerodatna odchylka generujiciho gaussovskeho bileho sumu
  sigma2=2.25; sigma=sqrt(sigma2);
% Generovani AR(2) procesu s parametry:
  p = poly([1,1]), % Polynom s koreny na jednotkove kruznici necht
                   % reprezentuje PHI(z) = 1 - (p(2)/p(3))z - (p(1)/p(3))z^2
  PHI = -[p(2),p(1)]./p(3);
  [x,phir,thetar]=armagen(PHI,[],n,sigma);
        % Koreny polynomu PHI(z) i THETA(z) nelezi vne jednotkoveho
  phir  % kruhu, takze proces neni kauzalni
  figure; plot(x);
  title('Nekauzalni proces');
% Proces vykazuje jasne nestacionarni chovani jiz pri tesnem poruseni
% kauzality.
%==================================KONEC======================================
echo off
