%=====================================
%=====================================
% CASOVE RADY SE STRUKTURALNIM ZLOMEM
%=====================================
%=====================================
close all
clear all
clc

load zlom.txt

ep = zlom(:,1);
y1 = zlom(:,2);
y2 = zlom(:,3);

%------------------------------------
figure
plot1=plot(y1);
set(plot1,'LineWidth',2)
title('Stacionarni casova rada se zmenou urovne v case t=51 (y_1)')
xlabel('cas')

figure
plot2=plot(y2,'r');
set(plot2,'LineWidth',2)
title('Casova rada (nestacionarni) se pulsni zmenou v case t=51 (y_2)')
xlabel('cas')

%------------------------------------
figure
subplot(2,2,1)
plotacf(y1,20)
title('acf pro y_1')
subplot(2,2,2)
plotpacf(y1,20)
title('pacf pro y_1')
subplot(2,2,3)
plotacf(diff(y1),20)
title('acf pro \Delta y_1')
subplot(2,2,4)
plotpacf(diff(y1),20)
title('pacf pro \Delta y_1')

figure
subplot(2,2,1)
plotacf(y2,20)
title('acf pro y_2')
subplot(2,2,2)
plotpacf(y2,20)
title('pacf pro y_2')
subplot(2,2,3)
plotacf(diff(y2),20)
title('acf pro \Delta y_2')
subplot(2,2,4)
plotpacf(diff(y2),20)
title('pacf pro \Delta y_2')
%------------------------------------

n=length(y1);
t=(1:n)';
k=ones(size(y1));
alfa = 0.05;

% y_t = a0 + a2 t + e_t
X=[k t];
[b,bi,e,ei,st] =  regress(y1,X);
disp('------------------------------------------')
disp('odhad regresni primky')
disp('------------------------------------------')
disp('parametry a jejich IS')
disp([b bi]) %a0<0, a2>0 
st
y1v = X*b;

pr1=mean(y1(1:50));
pr2=mean(y1(51:end));
pr=[ones(50,1)*pr1;ones(50,1)*pr2;];
figure
plot([y1,y1v,pr])
hold on
title('data a regresni primka')


%============================================================
%============================================================
%============================================================


% AR(1) proces y_t = a0 +a1 y_t-1 + e_t
dat=iddata(y1,ones(size(y1))); % funkce iddata, skladame data (simulovana) a vektor jednicek, ktery se nasobi s parametrem urovnove konstanty delta
th=armax(dat,[1 1 0 0]);     % odhad parametru, zde 2. vstup je [na nb nc nk], nk je zpozdeni na vstupu, trochu jina syntax nez v pripade procesu bez urovnove konstanty
disp('------------------------------------------')
disp('odhad AR(1) procesu pro y1')
disp('------------------------------------------')
present(th);
[p_sig,test,stud] = par_sig(th,alfa);
test
stud
[sig,tests,sk,QQ,chi]=res_sig(th,dat,alfa,0);
disp('QQ a chi')
disp([QQ chi])
% odhad parametru a1 je blizky nule 0.95



% AR(1) proces y_t = a0 +a1 y_t-1 + e_t + D_L
% D_L ... umela promenna (zmena urovne rady od casu 51 o jednotku)
D_L = [zeros(50,1);ones(50,1)];

dat=iddata(y1,[k D_L]); % funkce iddata, skladame data, vektor jednicek a umelou promennou
nb=[1 1]; % mame 2 externi vstupy: k, D_L
nk=[0 0]; % casovy posun externich vstupu (zpozedni)
na = 1; % rad AR casti
nc = 0; % rad MA casti
th=armax(dat,[na nb nc nk]);  % odhad parametru modelu

disp('-------------------------------------------------')
disp('odhad AR(1) procesu pro y1 s umelou promennou D_L')
disp('-------------------------------------------------')
present(th);
[p_sig,test,stud] = par_sig(th,alfa);
test
stud
[sig,tests,sk,QQ,chi]=res_sig(th,dat,alfa,0);
disp('QQ a chi')
disp([QQ chi])
% odhady parametru: a1=0.43, b1=0.035, b2=1.064
% a1 jiz neni blizky nule
% b1 je temer nula, nevyznamny
% b2 temer jednicka - parametr u umele promenne
% takto je to spravne


%============================================================
%============================================================
%============================================================


% y2 je nahodna prochazka s impulsnim skokem

% AR(1) proces y_t = a0 +a1 y_t-1 + e_t
dat=iddata(y2,ones(size(y2))); % funkce iddata, skladame data (simulovana) a vektor jednicek, ktery se nasobi s parametrem urovnove konstanty delta
th=armax(dat,[1 1 0 0]);     % odhad parametru, zde 2. vstup je [na nb nc nk], nk je zpozdeni na vstupu, trochu jina syntax nez v pripade procesu bez urovnove konstanty

disp('------------------------------------------')
disp('odhad AR(1) procesu pro y2')
disp('------------------------------------------')
present(th);
[p_sig,test,stud] = par_sig(th,alfa);
test
stud
[sig,tests,sk,QQ,chi]=res_sig(th,dat,alfa,0);
disp('QQ a chi')
disp([QQ chi])
% odhad parametru a1 je blizky nule 0.97



% AR(1) proces y_t = a0 +a1 y_t-1 + e_t + D_L
D_P  = zeros(size(y2));
D_P(51) = 2;
dat=iddata(y1,[k D_P]); % funkce iddata, skladame data, vektor jednicek a umelou promennou
nb=[1 1]; % mame 2 externi vstupy: k, D_L
nk=[0 0]; % casovy posun externich vstupu (zpozedni)
na = 1; % rad AR casti
nc = 0; % rad MA casti
th=armax(dat,[na nb nc nk]);  % odhad parametru modelu
disp('-------------------------------------------------')
disp('odhad AR(1) procesu pro y2 s umelou promennou D_P')
disp('-------------------------------------------------')
present(th);
[p_sig,test,stud] = par_sig(th,alfa);
test
stud
[sig,tests,sk,QQ,chi]=res_sig(th,dat,alfa,0);
disp('QQ a chi')
disp([QQ chi])


alfa = 0.05;
n = length(y2);
k = ones(size(y2));
t =(1:n)';
dy1 = y1 - lag(y1,1);
dy2 = y2 - lag(y2,1);


%============================================================
%============================================================
%============================================================
disp('')
disp('****************************************')
disp('DF testy bez zahrnuti umelych promennych')
disp('****************************************')
disp('')

%----------------------------------------
% DF testy s konstantou a trendem pro y1
%----------------------------------------
X   = [k t lag(y1,1)];
vys1 = regrese(dy1,X,alfa);

disp('---------------------------------------------------')
disp('1. model pro y1 - konstanta, trend, zadna zpozdeni')
disp('dy1_t = a0 + a2 t + gama y_{t-1} ')
disp('---------------------------------------------------')
disp('t-test (k t y_1)')
disp(vys1.t_test')
%disp('t_krit')
%disp(vys.t_krit)
disp('DF_krit n=100,250, a=0.05. je -3.45')
% vyslo 2.11 => nezam. H0 o unit root

%------------------------------
% TEST, zda tam patri trend
% H0: gamma=a2=0 => phi_3
%------------------------------
Xr    = [k];
vys1r = regrese(dy1,Xr,alfa);

SSE_u = vys1.sse;  %suma ctvercu chyb 1. modelu
SSE_r = vys1r.sse; %suma ctvercu chyb 2. modelu
p = vys1.n - vys1.k; %pocet stupnu volnosti modelu
r = 2; % pocet testovanych omezeni
phi_3 = (SSE_r - SSE_u)/(r*(SSE_u/p));

disp(' ')
disp('H0: gamma = a2 = 0')
disp('phi_3 a phi_3_krit')
disp([phi_3 7.44]) % 3.75 < 7.44

%------------------------------
% TEST, zda tam patri konstanta a0 a trend
%H0: a0=gamma=a2=0 => phi_3 nejde, nic by nezbylo, tj. jinak
%------------------------------
X2   = [lag(y1,1)];
vys2 = regrese(dy1,X2,alfa);

disp('---------------------------------------------------')
disp('2. model - bez konst., trendu, pro y1')
disp('dy1_t = gama y_{t-1} ')
disp('---------------------------------------------------')
disp('t-test ( y_1)')
disp(vys2.t_test')
disp('DF_krit n=100, a=0.05. je -1.95') %-0.98
disp(' ')
disp(' ')


%============================================================
%============================================================
%============================================================


%----------------------------------------
% DF testy s konstantou a trendem pro y2
%----------------------------------------
X   = [k t lag(y2,1)];
vys3 = regrese(dy2,X,alfa);

disp('--------------------------------------------------')
disp('1. model pro y2 - konstanta, trend, zadna zpozdeni')
disp('dy2_t = a0 + a2 t + gama y_{t-1} ')
disp('--------------------------------------------------')
disp('t-test (k t y_1)')
disp(vys3.t_test')
disp('DF_krit n=100,250, a=0.05. je -3.45')
% vyslo  -0.6585 => nezam. H0 o unit root


%------------------------------
% TEST, zda tam patri trend
% H0: gamma=a2=0 => phi_3
%------------------------------
Xr    = [k];
vys3r = regrese(dy2,Xr,alfa);

SSE_u = vys3.sse;  %suma ctvercu chyb 1. modelu
SSE_r = vys3r.sse; %suma ctvercu chyb 2. modelu
p = vys3.n - vys3.k; %pocet stupnu volnosti modelu
r = 2; % pocet testovanych omezeni
phi_3 = (SSE_r - SSE_u)/(r*(SSE_u/p));

disp(' ')
disp('H0: gamma = a2 = 0')
disp('phi_3 a phi_3_krit')
disp([phi_3 7.44]) % -1.8089 < 7.44

%------------------------------
% TEST, zda tam patri konstanta a0 a trend
%H0: a0=gamma=a2=0 => phi_3 nejde, nic by nezbylo, tj. jinak
%------------------------------
X2   = [lag(y2,1)];
vys2 = regrese(dy2,X2,alfa);

disp('---------------------------------------------------')
disp('2. model - bez konst., trendu, pro y2')
disp('dy1_t = gama y_{t-1} ')
disp('---------------------------------------------------')
disp('t-test ( y_1)')
disp(vys2.t_test')
disp('DF_krit n=100, a=0.05. je -1.95') %-0.3722
disp(' ')
disp(' ')


%============================================================
%============================================================
%============================================================
disp('')
disp('****************************************')
disp('DF testy se zahrnutim umelych promennych')
disp('****************************************')
disp('')


%-----------------------------------------
% DF testy se zahrnutim umelych promennych
%-----------------------------------------

% y1 je stac. rada s urovnovym posunem
% D_L = [zeros(50,1);ones(50,1)];
% y2 je nahodna prochazka s impulsnim skokem
% D_P  = zeros(size(y2));
% D_P(51) = 2;

% model :y1_t = a0 + a1 D_L + a2 t + e_t, pak DF test pro rezidua
X1   = [k D_L t];
vys1 = regrese(y1,X1,alfa);
e1 = vys1.e;
de1 = e1 - lag(e1,1);

disp('---------------------------------------------------')
disp('model: y1_t = a0 + a1 D_L + a2 t + e_t')
disp('---------------------------------------------------')
disp('t-test (k D_L t)')
disp(vys1.t_test')
disp('t_krit')
disp(vys1.t_krit)

X1e   = [lag(e1,1)];
vys1e = regrese(de1,X1e,alfa);
disp('DF test rezidui e1 (bez a0, a2), DF krit. je -1.95')
disp(vys1e.t_test) %-5.74, tj. je to jiz bily sum

figure
stem(e1)
title('rezidua modelu y_1_t = a_0 + a_1 D_L + a_2 t + e_t')


% model :y_2t = a0 + a1 D_P + a2 t + e_t, pak DF test pro rezidua
X2   = [k D_P t];
vys2 = regrese(y2,X2,alfa);
e2 = vys2.e;
de2 = e2 - lag(e2,1);

disp('---------------------------------------------------')
disp('model: y2_t = a0 + a1 D_P + a2 t + e_t')
disp('---------------------------------------------------')
disp('t-test (k D_L t)')
disp(vys2.t_test')
disp('t_krit')
disp(vys2.t_krit)

X2e   = [lag(e2,1)];
vys2e = regrese(de2,X2e,alfa);
disp('DF test rezidui e1 (bez a0, a2), DF krit. je -1.95')
disp(vys2e.t_test) %-0.56, tj. je to unit root proces

figure
stem(e2)
title('rezidua modelu y_2_t = a_0 + a_1 D_P + a_2 t + e_t')
