clear all
close all
clc

load smoke.mat;

%Obs:    807

data = smoke;

%   1. educ                     years of schooling
%   2. cigpric                  state cigarette price, cents per pack
%   3. white                    =1 if white
%   4. age                      in years
%   5. income                   annual income, $
%   6. cigs                     cigs. smoked per day
%   7. restaurn                 =1 if state restaurant smoking restrictions
%   8. lincome                  log(income)
%   9. agesq                    age^2
%  10. lcigpric                 log(cigprice)

educ    = data(:,1);
cigpric = data(:,2);
white   = data(:,3);
age     = data(:,4);
income  = data(:,5);
cigs    = data(:,6);
restaurn= data(:,7);
lincome = data(:,8);
agesq   = data(:,9);
lcigpric= data(:,10);



yy = cigs;
XX = [ones(size(yy)) log(income) log(cigpric) educ age age.^2 restaurn];

res_cigs = ols(yy,XX);
vnames = strvcat('cigs', 'konst', 'log_income', 'log_cigpric', 'educ', 'age', 'age^2', 'restaurn');
prt_reg(res_cigs,vnames)


%test heteroskedasticity
White = white_test(res_cigs,XX)
GQ = GQ_test(yy,XX,2,300,300)
    
Z = XX;
BP = BP_test(res_cigs,Z) 
BPKB = BPKB_test(res_cigs,Z) 
Glejser = glejser_test(res_cigs,Z)

%FGLS
log_u2 = log(res_cigs.resid.^2);
res_log_u2 = ols(log_u2,XX);
u2_hat = exp(res_log_u2.yhat); %1/vahy pro WLS resp. diagonalni prvky matice Omega v GLS

%wls - "rucne"
w = 1./u2_hat;
yyw = sqrt(w).*yy;
XXw = [sqrt(w).*ones(size(yyw)) sqrt(w).*log(income) sqrt(w).*log(cigpric) sqrt(w).*educ sqrt(w).*age sqrt(w).*age.^2 sqrt(w).*restaurn];
res_cigsw = ols(yyw,XXw);
prt_reg(res_cigsw,vnames)


%wls - funkce wls
res_cigsw2 = wls(yy,XX,diag(u2_hat));
my_prt(res_cigsw2,vnames)


%gls - funkce gls
res_cigsw3 = gls(yy,XX,diag(u2_hat));
my_prt(res_cigsw3,vnames)