colordef none crcv5 format compact us_deaths=[crpath,'/brockwel.dat/deaths.dat']; [x,xinf]=getdata(us_deaths); % nacteni a vykresleni n=length(x) n = 72 t=(1:n).'; pause; %--------------------------------------------------------------- % Nalezneme vyznamne periodicke komponenty pomoci periodogramu a % Fisherova-Siegelova testu: [pdg,kF,kS]=tperiod(x,figure); % periodogram pdg a jeho graf kF % indexy vyznamnych harmonickych slozek uzitim Fisherova testu kF = 6 1 12 30 3 18 2 24 kS % indexy vyznamnych harmonickych slozek uzitim Siegelova testu kS = 6 1 12 30 3 18 2 24 % Nejvyraznejsi slozka kF(1)=kS(1)=6 (6-ta harmonicka) indikuje periodu % 12 mesicu sezonni slozky. pause; %--------------------------------------------------------------- % Zlepsime citlivost testovani odectenim parabolickeho trendu: D=[t.^2 ,t,ones(n,1)]; % Cast matice planu pro parabolicky trend [ans,ans,ans,ans,yTr]=linregr(x,D(:,1:3),D(:,1:3)); % Jednoduseji totez uzitim 'polfs' - viz dale xTr=polfs(x,2,[],t); max(abs(xTr(:,1)-yTr(:,1))) % Kontrola shody obou vysledku ans = 1.8190e-011 % Fisherova-Siegelova testu: [pdg,kF,kS]=tperiod(x-xTr(:,1),figure); % periodogram pdg a jeho graf kF % indexy vyznamnych harmonickych slozek uzitim Fisherova testu kF = 6 12 30 18 3 kS % indexy vyznamnych harmonickych slozek uzitim Siegelova testu kS = 6 12 30 18 3 24 % Nejvyraznejsi slozka kF(1)=kS(1)=6 (6-ta harmonicka) indikuje periodu % 12 mesicu sezonni slozky. pause; %-------------------------------------------------------------- % do Fourierovy rady v linearnim modelu zahrneme nejprve jen % nejvyznamnejsi harmonicke komponenty 6,12,18,30: w = 2*pi*[6,12,18,30]./n; % uhlove rychlosti vybranych harm. komponent tt=linspace(1,n).'; N=length(tt); D=[t.^2 ,t,ones(n,1)]; % Cast matice planu pro parabolicky trend C=[tt.^2,tt,ones(N,1)]; % Cast matice param. funkci pro trend D=[D,cos(t*w),sin(t*w)] ; % Pridame do D cast pro Fourierovu radu C=[C,cos(tt*w),sin(tt*w)]; % Pridame do C cast pro Fourierovu radu [bE,rE,sE,PbE,CbE]=linregr(x,D,C); % Provedeme regresi tr=1:3; % indexy koeficientu paraboly trendu sz=4:size(C,2); % indexy Fourier. koeficientu a,a,,..; b,b,... harmon. slozek bE(tr,:) % a jejich odhady spolu s polosirkou intervalu spolehlivosti ans = 1.0e+003 * 0.0008 0.0002 -0.0721 0.0135 9.9619 0.2148 bE(sz,:) ans = -733.7835 98.1169 418.5698 98.0865 156.0846 98.0852 163.2807 98.0850 -758.6635 98.8367 77.3372 98.2018 -198.1549 98.0852 204.4604 98.0310 pause; %--------------------------------------------------------------- figure; plot(tt,CbE(:,1),'y',tt,[CbE(:,1)-CbE(:,2),CbE(:,1)+CbE(:,2)],'r-'); title('trend + sezonni slozka'); pause; %--------------------------------------------------------------- xtrf=C(:,tr)*bE(tr,1); % Jen trend figure; plot(tt,xtrf); title('trend'); pause; %--------------------------------------------------------------- xszf=C(:,sz)*bE(sz,1); % Jen sezonni slozka figure; plot(tt,xszf); title('sezonni slozka'); pause; %--------------------------------------------------------------- % Vyse uvedeny postup je pohodlneji implementovan v procedure 'polfs', % kterou nyni pouzijeme pro dekompozici s indexy kS, kterych je vice nez v kF dt = 1/12; t = (1973:dt:1979-dt).'; % Casovou osu definujeme dle skutecnosti tt = linspace(1973,1979-dt).'; [yE,trE,pE,r,s,bE,cE,sE]=polfs(x,2,n*dt./kS,t,tt,[],figure); bE, % parametry parabolickeho trendu + 95% polosirka intervalu spol. bE = 1.0e+003 * 0.1171 0.0221 -0.1239 0.0357 8.4364 0.0885 cE, % koeficienty u kosinu ve F.r. + 95% polosirka intervalu spol. cE = -904.9425 83.2033 401.2074 83.1757 239.7575 83.1740 -29.7291 83.1743 -67.5980 83.6505 130.6105 83.1741 sE, % koeficienty u sinu ve F.r. + 95% polosirka intervalu spol. sE = -537.3115 83.9643 -139.4185 83.3865 -103.3025 83.2310 -248.4413 83.2803 -146.1104 86.2397 84.2279 83.2449 pause; %--------------------------------------------------------------- % Alternativne muzeme pro dekompozici pouzit proceduru 'polper', ktera % vsak umoznuje predikci na siti tt pouze pro trend. [yE,trE,pE,r,s,bE]=polper(x,2,n/kS(1),t,tt,[],figure); bE, % parametry parabolickeho trendu + 95% polosirka intervalu spol. bE = 1.0e+003 * 0.1190 0.0242 -0.1397 0.0380 7.6078 0.2381 pause; %--------------------------------------------------------------- % Porovname graficky vysledky dekompozice z obou procedur na siti t: [yE1,trE1,pE1,r1,s1,bE1,cE1,sE1]=polfs(x,2,n*dt./kS,t); [yE2,trE2,pE2,r2,s2,bE2]=polper(x,2,n/kS(1),t); bE1, bE2, % Porovnejte koeficienty parabolickeho trendu bE1 = 1.0e+003 * 0.1171 0.0221 -0.1239 0.0357 8.4364 0.0885 bE2 = 1.0e+003 * 0.1190 0.0242 -0.1397 0.0380 7.6078 0.2381 pause; figure; plot(t,[yE1(:,1),yE2(:,1)]); title('Porovnani vyhlazeni'); pause; figure; plot(t,[trE1(:,1),trE2(:,1)]); title('Porovnani trendu'); pause; figure; plot(t,[pE1(:,1),pE2(:,1)]); title('Porovnani periodicke slozky'); pause; figure; plot(t,[r1(:,1),r2(:,1)]); title('Porovnani chybove slozky'); %--------------------------------------------------------------------------- % PROSTUDUJTE ALGORITMUS PROCEDUR 'polfs' a 'polper' pomoci prikazu 'type'. %================================== KONEC ================================== echo off close all help szsmtr ESTIMATING SEASONAL COMPONENT BY THE SMALL TREND METHOD ======================================================================== function [sz,tr]=szsmtr(x,d,model,fig) INPUT: x ... data vector n x 1 or 1 x n to be processed d ... period of the seasonal component (scalar) model ... ==0 ... additive model (default) ~=0 ... multiplicative model fig ... figure handle for the plot of results (optional: default=[]=no plot) OUTPUT: sz ... seasonal component of the same length as x (column vector) tr ... column vector of the same length as x: help szma ESTIMATING SEASONAL COMPONENT BY THE MOVING AVERAGE METHOD ======================================================================== function [sz,tr,msz]=szma(x,d,model,fig) INPUT: x ... data vector n x 1 or 1 x n to be processed d ... period of the seasonal component (scalar), d=2*q or d=2*q+1 model ... ==0 ... additive model (default) ~=0 ... multiplicative model fig ... figure handle for the plot of results (optional: default=[]=no plot) OUTPUT: sz ... seasonal component of the same length as x (column vector) tr ... column vector of the same length as x: tr(q+1:n-q)=preliminary trend estimate obtained from x via moving average. tr(1)=...=tr(q)=tr(q+1), tr(n-q)=tr(n-q+1)=...=tr(n). msz... additive or multiplicative trend correction to obtain mean(sz)=0 for the additive model or mean(sz)=1 for the multiplicative model crcv6 format compact format long us_deaths=[crpath,'/brockwel.dat/deaths.dat']; [x,xinf]=getdata(us_deaths); % nacteni a vykresleni n=length(x) n = 72 t=(1:n).'; pause; %--------------------------------------------------------------- % Jiz znamym postupem identifikujeme periodicke komponenty napr. % pomoci Fisherova testu: xTr=polfs(x,2,[],t); [pdg,kF]=tperiod(x-xTr(:,1),figure); % periodogram pdg a jeho graf kF % indexy vyznamnych harmonickych slozek uzitim Fisherova testu kF = 6 12 30 18 3 % Nejvyraznejsi slozka kF(1)=kS(1)=6 (6-ta harmonicka) indikuje periodu % 12 mesicu sezonni slozky. % pripadnou zakladni periodu kF(#)=1 budeme ignorovat: % kF(#)=[] pause; %-------------------------------------------------------------- % Pro dekompozici pouzijeme proceduru 'polfs'. % Predpokladejme nejprve stupen 0 (konstantni trend): [yE,trE,pE,r,s,bE]=polfs(x,0,n./kF,[],[],[],figure); bE, % hodnota konstatniho trendu + 95% polosirka intervalu spol. bE = 1.0e+003 * 8.78773611111111 0.11494377623347 % Rezidualni chybova slozka nevykazuje nahodne oscilace kolem nuly; % stupen je tedy nutno zvetsit - predpokladejme linearni trend: pause; %--------------------------------------------------------------- close all [yE,trE,pE,r,s,bE]=polfs(x,1,n./kF,[],[],[],figure); bE, % koeficienty linearniho trendu + 95% polosirka intervalu spol. bE = 1.0e+003 * -0.01045148423449 0.00520023639692 8.78773611111111 0.10290280110268 % 95% interval spolehlivosti pro koeficient u prvni mocniny je: [bE(1,1)-bE(1,2),bE(1,1)+bE(1,2)] ans = -15.65172063140605 -5.25124783757277 % V tomto intervalu nelezi hodnota 0, takze s rizikem 5% zamitame % hypotezu o nulovosti tohoto koeficientu a tedy predpoklad konstantniho % trendu; % To potvrzuje i prubeh rezidualni chybove slozky, ktera stale nevykazuje % nahodne oscilace kolem nuly. % Predpokladejme proto parabolicky trend: pause; %--------------------------------------------------------------- close all [yE,trE,pE,r,s,bE]=polfs(x,2,n./kF,[],[],[],figure); bE, % koeficienty linearniho trendu + 95% polosirka intervalu spol. bE = 1.0e+003 * 0.00081364056207 0.00016794206233 -0.01045148423449 0.00325785612783 8.43631119167784 0.09704421262971 % 95% interval spolehlivosti pro koeficient u druhe mocniny je: [bE(1,1)-bE(1,2),bE(1,1)+bE(1,2)] ans = 0.64569849973690 0.98158262439941 % V tomto intervalu nelezi hodnota 0, takze s rizikem 5% zamitame % hypotezu o nulovosti tohoto koeficientu a tedy take predpoklad % linearniho trendu; % Protoze rezidualni chybova slozka jiz nahodne osciluje kolem nuly, % lze se duvodne domnivat, ze trend je parabolicky. % Pro kontrolu zkusime kubicky trend: pause; %--------------------------------------------------------------- close all [yE,trE,pE,r,s,bE,cE,sE]=polfs(x,3,n./kF,[],[],[],figure); bE, % koeficienty linearniho trendu + 95% polosirka intervalu spol. bE = 1.0e+003 * 0.00000225268844 0.00001019969048 0.00081364056207 0.00016915999850 -0.01209771784956 0.00814414867930 8.43631119167784 0.09774798900844 % 95% interval spolehlivosti pro koeficient u treti mocniny je: [bE(1,1)-bE(1,2),bE(1,1)+bE(1,2)] ans = -0.00794700204645 0.01245237891875 % V tomto intervalu lezi hodnota 0, takze nezamitame % hypotezu o nulovosti tohoto koeficientu a definitivne prijimame % predpoklad parabolickeho trendu. pause; close all; %--------------------------------------------------------------- % Jina ponekud mene exaktni metoda odhadu stupne polynomialniho % trendu je zalozena na diferencovani. % Je pouzitelna pouze, pokud data obsahuji nejvyse jednu periodickou slozku % o zname zakladni periode (zpravidla sezonni slozku). % 1.krok: % Jestlize data obsahuji periodickou slozku, odstranime ji diferencovanim % na vzdalenost rovnajici se jeji zakladni periode d; % pro nas soubor je d=12; xd12 = x(d+1:n)-x(1:n-d); figure; plot(xd12); title('Prubeh po odstraneni sezonni slozky diferencovanim'); % Timto diferencovanim jsme vsak snizili o 1 i stupen polynomialniho trendu. % Chybovy prubeh neosciluje kolem konstantni hodnoty. % Trend tedy nemohl byt ani konstantni (oscilace by musely byt kolem nuly), % ani linearni (oscilace by musely byt kolem nenulove konstantni hodnoty). pause; %--------------------------------------------------------------- % 2.krok: % Snizime stupen dalsim diferencovanim, tentokrat jiz standardne na % vzdalenost 1: nd12=length(xd12); xd12d1=x(2:nd12)-x(1:nd12-1); figure; plot(xd12d1); title('Prubeh po dvojim diferencovani'); % Chybovy prubeh jiz osciluje kolem konstantni hodnoty. % Trendovou komponentu proto muzeme pokladat za parabolickou. % Kdyby tomu tak nebylo, zopakovali bychom krok 2 atd. % V pripade, ze data neobsahuji sezonni slozku, pak krok 1 samozrejme % vypustime a zacneme provadet rovnou krok 2. pause; close all; %=========================================================================== % Pro separaci sezonni slozky pouzijeme nyni metodu maleho trendu % implementovanou v procedure 'szsmtr': AM = 0; d=12; % 0=aditivni model; 12=perioda [sz,tr] = szsmtr(x,12,AM,figure); % Metoda maleho trendu pause; %--------------------------------------------------------------- % Alternativne lze pouzit presnejsi metodu klouzaveho prumeru % implementovanou v procedure 'szma': [sz,tr,msz] = szma(x,12,AM,figure); % Metoda klouzaveho prumeru msz, % Aditivni korekce trendu je zanedbatelna msz = -9.53888888888843 pause; close all; %=============================================================== % Predchozi metody lze uzit i primo pro pripad multiplikativniho modelu. % Pro tento ucel uzijeme % SOUBOR: airpass.dat % Cestujici v mezinarodni letecke doprave % x=mesicni souhrny v tisicich cestujicich % t=1,2,...,144: leden 1949 az prosinec 1960 %=============================================================== air_pass=[crpath,'/brockwel.dat/airpass.dat']; [x,xinf]=getdata(air_pass); % nacteni a vykresleni n=length(x); [l,y]=powtr(x,4); >>> LINEAR GROWTH OF VARIANCE (logarithmic transform): y=log(x) <<< pause(2); figure; plot(y); title('Po logaritmicke transformaci'); xlabel(xinf(2,:)); % Zavislost rozptylu na trendu je linearni, takze lze predpokladat, % ze se data ridi multiplikativnim modelem. pause; %--------------------------------------------------------------- % Zjistime periodu sezonni slozky: [pdg,kF]=tperiod(y,figure); % Model musi byt aditivni, tj. uzijeme y % po logaritmicke transformaci. pause(2) % Presneji odectenim parabolickeho trendu: yTr=polfs(y,2,[],(1:n).'); [pdg,kF]=tperiod(y-yTr(:,1),figure); kF kF = 12 24 3 48 36 60 11 13 6 23 % kF(1)=12 je harmonicka slozka odpovidajici rocni sezonni periode, % tj. za periodu zvolime opet 12 mesicu. pause; %--------------------------------------------------------------- MM = 1; d=12; % 1=multiplikativni model; 12=perioda [sz,tr] = szsmtr(x,12,MM,figure); % Metoda maleho trendu pause; %--------------------------------------------------------------- [sz,tr,msz] = szma(x,12,MM,figure); % Metoda klouzaveho prumeru msz, % Multiplikativni korekce trendu je zanedbatelna msz = 0.99823565827641 %--------------------------------------------------------------- % PROSTUDUJTE ALGORITMUS PROCEDUR 'szsmtr' a 'szma' pomoci prikazu 'type'. %================================== KONEC ================================== echo off diary off