%simulace heteroscedasticita a jeji testovani

clear all
close all
clc
 
niter=100; %pocet iteraci
n = 100; k=3; %pocet pozorovani, pocet vysvetlujicich promennych

White_hom = 0; %pocet nezamitnuti hypotezy o homoskedasticite (White-test)
GQ_hom = 0;    %pocet nezamitnuti hypotezy o homoskedasticite (GQ-test)
BP_hom = 0;    %pocet nezamitnuti hypotezy o homoskedasticite (BP-test)
BPKB_hom = 0;  %pocet nezamitnuti hypotezy o homoskedasticite (BPKP-test)
Glejser_hom = 0;  %pocet nezamitnuti hypotezy o homoskedasticite (Glesjer-test)
        Glejser_hom_e2 = 0;
        Glejser_hom_abse = 0;
        Glejser_hom_logabse = 0;  

tic %spusteni pocitadla
for iter=1:niter
    rand('seed',10); %vyber generatoru
    
    x1 = (randn(n,1)*10+6);   %1. vysvetlujici promenna ~N(6,100)
    x2 = randn(n,1)*4+3;    %2. vysvetlujici promenna ~N(3,16)
    
    %volba typu nahodne slozky
    e = randn(n,1).*x1;  % heteroscedastic error term (podle x1)
    %e = randn(n,1);        % homoscedastic error term
    
    b = ones(k,1);
    iota = ones(n,1);
    x = [iota x1 x2];

    % generate y-data
    y = x*b + e;

    % do ols regression
    reso = ols(y,x);

    White = white_test(reso,x);%viz. help white_test
    GQ = GQ_test(y,x,2,0.3*n,0.3*n); %viz. help GQ_test
    
    Z = [iota x1 x2];
    BP = BP_test(reso,Z); %viz. help BP_test
    BPKB = BPKB_test(reso,Z); %viz. help BPKB_test
    Glejser = glejser_test(reso,Z); %viz. help Glejser_test
    
    White_hom = White_hom + White.homosc;
    GQ_hom = GQ_hom + GQ.homosc;
    BP_hom = BP_hom + BP.homosc;
    BPKB_hom = BPKB_hom + BPKB.homosc;
    Glejser_hom = Glejser_hom + Glejser.homosc;
        Glejser_hom_e2 = Glejser_hom_e2 + Glejser.homosc_e2;
        Glejser_hom_abse = Glejser_hom_abse + Glejser.homosc_abse;
        Glejser_hom_logabse = Glejser_hom_logabse + Glejser.homosc_abse;
    
    clc  
    fprintf('Probihaji iterace ... \n');
    fprintf(['Iterace cislo ', num2str(iter),'/',num2str(niter)]);
end

%zobrazeni vysledku
fprintf('\n Celkovy cas iterace %6.2f sekund \n',toc);
fprintf('\n Testovani hypotezy o homoskedasticite z %6.0f iteraci \n',niter);
fprintf('\n Delka vzorku %6.0f \n',n);
fprintf('\n Whiteuv test:          prijata %6.0f        zamitnuta %6.0f \n',[White_hom, niter-White_hom]);
fprintf('\n Golden-Quandtuv test:  prijata %6.0f        zamitnuta %6.0f \n',[GQ_hom, niter-GQ_hom]);
fprintf('\n Breusch-Paganuv test:  prijata %6.0f        zamitnuta %6.0f \n',[BP_hom, niter-BP_hom]);
fprintf('\n Breusch-Paganuv test:  prijata %6.0f        zamitnuta %6.0f \n',[BPKB_hom, niter-BPKB_hom]);
fprintf(' (Koenker-Bassettova varianta) \n');
fprintf('\n Glejseruv test:            prijata %6.0f        zamitnuta %6.0f \n',[Glejser_hom, niter-Glejser_hom]);
fprintf('\n Glejseruv test e^2:        prijata %6.0f        zamitnuta %6.0f \n',[Glejser_hom_e2, niter-Glejser_hom_e2]);
fprintf('\n Glejseruv test abs(e):     prijata %6.0f        zamitnuta %6.0f \n',[Glejser_hom_abse, niter-Glejser_hom_abse]);
fprintf('\n Glejseruv test:log(abs(e)) prijata %6.0f        zamitnuta %6.0f \n',[Glejser_hom_logabse, niter-Glejser_hom_logabse]);