close all
clear all
clc
format short g

%rhdp USA, billions of 1996 dollars
load gdp_usa.txt;
gdp = gdp_usa;
roky = (1960.25:0.25:2002.25)';

figure
plot(roky,gdp)
title('rgdp USA')
xlabel('cas')
ylabel('billions of 1996 dollars')

n = length(roky);
t = (1:n)';
X=[ones(size(t)) t t.^2 t.^3];
[b,bint,e,eint,stat] = regress(gdp,X);
gdp_vyr = X*b;
af = acf(e,5);

disp('************************************')
disp('Odhad polynomialniho trendu rGDP USA')
disp('************************************')
disp('parametry a jejich IS')
disp([b bint])
disp('R^2, F_R, p-value, sig^2_e')
disp(stat)
disp('autokorelacni funkce rezidui')
disp(af')

figure
subplot(2,1,1)
plotacf(e,20)
title('autokorelacni funkce detrendovaneho rGDP')
subplot(2,1,2)
plotpacf(e,20)
title('parc. autokorelacni funkce detrendovaneho rGDP')

figure
p1 = plot(roky,gdp,roky,gdp_vyr,'--');
set(p1,'LineWidth',2)
title('rgdp USA - data a vyrovnane hodnoty')
xlabel('cas')
ylabel('billions of 1996 dollars')
legend('data','modelova predikce')

% predpoved do budoucna
% h = 8;
% tp = (n+1:n+h)';
% Xp = [ones(size(tp)) tp tp.^2 tp.^3];
% gdp_pred = Xp*b;
% rokypr = (roky(end)+0.25:0.25:roky(end)+h*0.25)';
% 
% figure
% p2 = plot(roky,gdp,roky,gdp_vyr,'--',rokypr,gdp_pred);
% set(p2(1),'LineWidth',2)
% set(p2(3),'LineWidth',2)
% title('rgdp USA - data, vyrovnane hodnoty, predikce')
% xlabel('cas')
% ylabel('billions of 1996 dollars')
% legend('data','modelova predikce','budoucnost','Location','best')



%-------------------------------------------------------
%-------------------------------------------------------
% 1. MODEL S KONSTANTOU a S TRENDEM
%-------------------------------------------------------
%-------------------------------------------------------
% overime, zdali jsme meli zlogaritmovana data jeste diferencovat, tj.
% provedeme DF test, nejprve pro rovnici s konstantou i trendem
% y = log(gdp)
% dy_t = a0 + a2 t + gama y_t-1 + b1 dy_t-1 + b2 dy_t-2
y   = log(gdp); % 
dy = y - lag(y,1);
n = length(y);
t = (1:n)';
k = ones(size(y)); 
alfa = 0.05;

X=[k t lag(y,1) lag(dy,1) lag(dy,2)];
vys = regrese(dy,X,alfa);
% funkci "regrese" je treba mit ve svem pracovnim adresari nebo pripojenou
% na ceste, jeji pouziti - viz "help regrese"

disp('*****************************************************')
disp('1. MODEL S KONSTANTOU A TRENDEM (y = log(gdp))')
disp('dy_t = a0 + a2 t + gama y_t-1 + b1 dy_t-1 + b2 dy_t-2')
disp('*****************************************************' )
disp('parametry'); disp(vys.b')
disp('t-test (k t y_1 dy_1 dy_2))'); disp(vys.t_test')
disp('t_krit'); disp(vys.t_krit)
disp('DF_krit n=100, a=0.1, 0.05 je -3.15, -3.45')
% test gamma = 0
% t-test je -3.37, krit hodnota pro 100 dat a 10% a 5% jr -3.15 1-3.45.
% pro a=10% unit root neni
% pro a=5% nezamitame H0, ze je unit root



% nyni otestujeme opravnenost pritomnosti clenu s trendem
% t-test je pro a2 roven 3.27, ale neporovnavame se Stud. kvantilem!!! ale 
% vypocteme statistiku phi_3 pro H0: gamma=a2=0

%-----------------------------------
% 1. F-TEST: H0: gamma=a2=0 => PHI_3
%-----------------------------------
% unrestricted model uz mame vyse
% restricted model
Xr=[k lag(dy,1) lag(dy,2)];

vys2 = regrese(dy,Xr,alfa);

SSE_u = vys.sse;  %suma ctvercu chyb modelu bez omezeni
SSE_r = vys2.sse; %suma ctvercu chyb modelu s omezenimi
p = vys.n - vys.k; %pocet stupnu volnosti modelu
r = 2; % pocet testovanych omezeni
phi_3 = (SSE_r - SSE_u)/(r*(SSE_u/p));

disp('  ')
disp('phi_3 a phi_3_krit (H0: gamma=a2=0)')
disp([phi_3 7.44])
% vyslo 6.14, krit. je 7.44 => H0 nezamitame, tj. ze trend do modelu nepatri a model
% odhadneme znova, tentokrate bez trendu



%-------------------------------------------------------
%-------------------------------------------------------
% 2. MODEL S KONSTANTOU A BEZ TRENDU
%-------------------------------------------------------
%-------------------------------------------------------
% y = log(gdp)
% dy_t = a0 + gama y_t-1 + b1 dy_t-1 + b2 dy_t-2
X=[k lag(y,1) lag(dy,1) lag(dy,2)];
vys3 = regrese(dy,X,alfa);

disp('*************************************************')
disp('2. MODEL S KONSTANTOU A BEZ TRENDU (y = log(gdp))')
disp('dy_t = a0 + gama y_t-1 + b1 dy_t-1 + b2 dy_t-2')
disp('*************************************************' )
disp('parametry'); disp(vys3.b')
disp('t-test (k y_1 dy_1 dy_2)'); disp(vys3.t_test')
disp('t_krit'); disp(vys3.t_krit)
disp('DF_krit n=100, a=0.1, 0.05 je -2.89, -2.58')
% H0: gamma=0: vyslo -1.22, krit. h. tau_mu je -2.89, tj. nezam. H. unit root


% nyni otestujeme opravnenost pritomnosti urovnove konstanty
%-----------------------------------
% 2. F-TEST: H0: gamma=a0=0 => PHI_1
%-----------------------------------
% H0: a0=gamma=0:

% restricted model
Xr=[lag(dy,1) lag(dy,2)];
vys4 = regrese(dy, Xr, alfa);

SSE_u = vys3.sse;  %suma ctvercu chyb modelu bez omezeni
SSE_r = vys4.sse; %suma ctvercu chyb modelu s omezenimi
p = vys.n - vys.k; %pocet stupnu volnosti modelu
r = 2; % pocet testovanych omezeni
phi_1 = (SSE_r - SSE_u)/(r*(SSE_u/p));

disp('  ')
disp('phi_1 a phi_1_krit (H0: gamma=a0=0)')
disp([phi_1 5.57])

% 13.35>5.57 (krit. h. phi_1), tj. konstanta a0 do modelu patri
% Jelikoz nezamitame a0=0, muzeme testovat gamma=0 s norm. rozlozenim, 
%(vysledek jiz nebude zkreslen neadekvatnosti tvaru modelove rovnice) 
% t-test=-1.22 => neni st. v. ruzny od nuly. 
%
% Tj. casova rada ma unit root a drift term.



%--------------------------------------------------------
%--------------------------------------------------------
% 3. MODEL S KONST. A BEZ TRENDU PRO ZDIFERENCOVANOU RADU
% (overeni, ze stacilo diferencovat 1x)
%--------------------------------------------------------
%--------------------------------------------------------
% y = log(gdp)
% ddy_t = a0 + gama dy_t-1 + b1 ddy_t-1 + b2 ddy_t-2
ddy = dy - lag(dy,1);
X=[k lag(dy,1) lag(ddy,1) lag(ddy,2)];
vys5 = regrese(ddy,X,alfa);

disp('*************************************')
disp('3. ZDIFERENCOVANA DATA MODELU 2')
disp('(overeni, ze se melo diferencovat 1x)')
disp('*************************************' )
disp('parametry'); disp(vys5.b')
disp('t-test (k dy_1 ddy_1 ddy_2)'); disp(vys5.t_test')
disp('t_krit'); disp(vys5.t_krit)
disp('DF_krit n=100, a=0.1, 0.05 je -3.15, -3.45')
disp('  ')
%t-test pro gamu je -6.25, tj. zamitame hypotezu jednotkoveho korene,
% rada "dy" je jiz zrejme stacionarni.



%-------------------------------------------------------
%-------------------------------------------------------
%-------------------------------------------------------
% Tedy pri modelovani B-J pristupem budeme pracovat s tempem rustu GDP,
% jelikoz log(GDP) obsahuje jednotkovy koren a rada diff(log(GDP)) uz ne.

% Napr. odhad modelu AR(2) s urovnovou konstantou by vypadal takto:
% pri oznaceni: y = log(gdp), dy = diff(log(GDP))
% POZOR, v "dy" je na zacatku 1 NaN, tj. musime ho oriznout.

disp('******************************************' )
disp('Odhad AR(2) procesu s urovnovou konstantou')
disp('******************************************' )
dy = dy(2:end);
alfa=0.05;
dat=iddata(dy,ones(size(dy))); % funkce iddata, skladame data (simulovana) a vektor jednicek, ktery se nasobi s parametrem urovnove konstanty delta
th=armax(dat,[2 1 0 0]);
present(th);

[p_sig,test,stud] = par_sig(th,alfa); 
[sig,tests,sk,Q,chi]=res_sig(th,dat,alfa,0); 

disp('p_sig     test'); disp([p_sig   test])
disp('stud'); disp([stud])
disp('Q           chi'); disp([Q chi])
%-------------------------------------------------------