close all;
clear all;
clc;

% USGas_Greene.dat dokumentace

% data ke spotrebe benzinu v USA
% Greene (2000)

load USGas_Greene.dat;

[n nvar] = size(USGas_Greene);

Year =USGas_Greene(1:n,1);  % Rok
G   = USGas_Greene(1:n,2);  % Celkova spotreba benzinu v USA za rok      [celkove vydaje/cenovy index]
Pg   = USGas_Greene(1:n,3); % Cenovy index benzinu                       [1967 = 1]
Y   = USGas_Greene(1:n,4);  % Disponibilni prijem na obyvatele
Pnc  = USGas_Greene(1:n,5); % Cenovy index pro nova auta                 [1967 = 1]
Puc  = USGas_Greene(1:n,6); % Cenovy index pro stara auta                [1967 = 1]
Ppt = USGas_Greene(1:n,7);  % Cenovy index verejne dopravy               [1967 = 1]
Pd  = USGas_Greene(1:n,8);  % Agregovany cenovy index pro durables       [1982=1]
Pn  = USGas_Greene(1:n,9);  % Agregovany cenovy index pro nondurables    [1982=1]
Ps  = USGas_Greene(1:n,10); % Agregovany cenovy index pro sluzby         [1982=1]
Pop  = USGas_Greene(1:n,11);% Americka populace                          [mil.]

const = ones(n,1);

%(a)
YY = G./(Pop);
XX = [const Pg Y Pnc Puc Ppt Pd Pn Ps Year];
result_a = ols(YY,XX);
vnames_a = strvcat('G_Pop','Konstanta','Pg','Y','Pnc','Puc','Ppt','Pd','Pn','Ps','Trend');
prt_reg(result_a,vnames_a)



%(b)
cov_b = result_a.sige*inv(XX'*XX);
t_stat = (result_a.beta(4)-result_a.beta(5))/sqrt(cov_b(4,4)+cov_b(5,5)-2*cov_b(4,5))
t_krit = tinv(0.025,result_a.nobs-result_a.nvar)

%(c)
eta_Pg = result_a.beta(2)*(XX(:,2)./YY);
eta_Y = result_a.beta(3)*(XX(:,3)./YY);
eta_Ppt = result_a.beta(6)*(XX(:,6)./YY);

%prumerne elasticity
eta_Pg_mean = result_a.beta(2)*(mean(XX(:,2))./mean(YY))
eta_Y_mean = result_a.beta(3)*(mean(XX(:,3))./mean(YY))
eta_Ppt_mean = result_a.beta(6)*(mean(XX(:,6))./mean(YY))



figure
plot(Year,[eta_Pg,eta_Y,eta_Ppt]);
legend('eta Pg','eta Y','eta Ppt');

%(d)
yy = log(YY);
xx = [const log(XX(:,2:end-1)) Year];
result_d = ols(yy,xx);
vnames_d = strvcat('g_pop','Konstanta','p_g','y','p_nc','p_uc','p_pt','p_d','p_n','p_s','Trend');
prt_reg(result_d,vnames_d)

%(e)
nPg = Pg./Pg(1982-1967+1);
nPnc = Pnc./Pnc(1982-1967+1);
nPuc = Puc./Puc(1982-1967+1);
nPpt = Ppt./Ppt(1982-1967+1);
YY = G./(Pop);
XX = [const nPg Y nPnc nPuc nPpt Pd Pn Ps Year];
result_e = ols(YY,XX);
vnames_e = strvcat('G_Pop','Konstanta','nPg','Y','nPnc','nPuc','nPpt','Pd','Pn','Ps','Trend');
prt_reg(result_e,vnames_e)