clear all
close all
clc

n = 1000;

%exogenni promenne - mozno generovat libovolne
G = 10 + randn(n,1);
R = 4 + 2*randn(n,1);

%nahodne slozky - dobre pohrat si s rozptylem pro odliseni odhadu OLS a
%2SLS!!!
eps_C = 12*randn(n,1);
eps_Y = 4*randn(n,1);

%endogenni promenne C a Y
%soustava:  1. rovnice C = a1 + a2*Y + a3*C_1 + a4*G + eps_C
%           2. rovnice Y = b1 + b2*C + b3*Y_1 + b4*Y_2 + b5*R + eps_Y
%           3. identita Y_1 = Y_1

a1 = 1; a2 = 0.2; a3 = 0.8; a4 = 0.5; 
b1 = 1; b2 = 0.4; b3 = 0.3; b4 = 0.2; b5 = -2.3;

%Ay = By_1 + C_X + eps;
AA = [1 -a2 0
     -b2 1 -b3
     0   0    1];

BB = [a3 0  0
     0  0 b4
     0  1  0];

CC = [a1 a4 0
      b1 0 b5
      0  0 0];
 
Y = zeros(n,1);
Y_1 = zeros(n,1);
Y_2 = zeros(n,1);
C = zeros(n,1);
C_1 = zeros(n,1);


for i=1:n
    if i > 1
        Y_1(i)=Y(i-1);
        C_1(i)=C(i-1);
    end
    if i>2
        Y_2(i)=Y(i-2);
    end
    

    YY = inv(AA)*BB*[C_1(i); Y_1(i); Y_2(i)] + inv(AA)*CC*[1; G(i); R(i)] + inv(AA)*[eps_C(i); eps_Y(i); 0];
    C(i) =YY(1);
    Y(i) =YY(2);
    Y_1(i) =YY(3);
end

C_model = C;
Y_model = Y;

delta = inv(AA)*BB;
lambda = eig(delta);
abs_lambda = abs(lambda);
disp('Lambda')
disp(lambda)
disp('Abs. lambda')
disp(abs_lambda)




%Odhady
iota = ones(n,1);
%odhad
%1. rovnice
 resultC_OLS = ols(C,[Y iota C_1 G]);
 resultC_2SLS = tsls(C,Y,[iota C_1 G],[C_1 Y_1 Y_2 iota G R]);
 vnamesC = strvcat('C','a2','a1','a3','a4');
 prt_reg(resultC_OLS,vnamesC);
 prt_reg(resultC_2SLS,vnamesC);
%2. rovnice
 resultY_OLS = ols(Y,[C iota Y_1 Y_2 R]);
 resultY_2SLS = tsls(Y,C,[iota Y_1 Y_2 R],[C_1 Y_1 Y_2 iota G R]);
 vnamesY = strvcat('Y','b2','b1','b3','b4','b5');
 prt_reg(resultY_OLS,vnamesY);
 prt_reg(resultY_2SLS,vnamesY);

 
%SIMULACE OLS
a1 = resultC_OLS.beta(2); a2 = resultC_OLS.beta(1); a3 = resultC_OLS.beta(3); a4 = resultC_OLS.beta(4); 
b1 = resultY_OLS.beta(2); b2 = resultY_OLS.beta(1); b3 = resultY_OLS.beta(3); b4 = resultY_OLS.beta(4); b5 = resultY_OLS.beta(5);

%Ay = By_1 + C_X + eps;
AA = [1 -a2 0
     -b2 1 -b3
     0   0    1];

BB = [a3 0  0
     0  0 b4
     0  1  0];

CC = [a1 a4 0
      b1 0 b5
      0  0 0];
 
Y = zeros(n,1);
Y_1 = zeros(n,1);
Y_2 = zeros(n,1);
C = zeros(n,1);
C_1 = zeros(n,1);


for i=1:n
    if i > 1
        Y_1(i)=Y(i-1);
        C_1(i)=C(i-1);
    end
    if i>2
        Y_2(i)=Y(i-2);
    end
    

    YY = inv(AA)*BB*[C_1(i); Y_1(i); Y_2(i)] + inv(AA)*CC*[1; G(i); R(i)] + inv(AA)*[eps_C(i); eps_Y(i); 0];
    C(i) =YY(1);
    Y(i) =YY(2);
    Y_1(i) =YY(3);
end

C_OLS = C;
Y_OLS = Y;


delta = inv(AA)*BB;
lambda = eig(delta);
disp('Lambda - OLS')
disp(lambda)
abs_lambda = abs(lambda);
disp('Abs. lambda  - OLS')
disp(abs_lambda)


%SIMULACE 2SLS

a1 = resultC_2SLS.beta(2); a2 = resultC_2SLS.beta(1); a3 = resultC_2SLS.beta(3); a4 = resultC_2SLS.beta(4); 
b1 = resultY_2SLS.beta(2); b2 = resultY_2SLS.beta(1); b3 = resultY_2SLS.beta(3); b4 = resultY_2SLS.beta(4); b5 = resultY_2SLS.beta(5);

%Ay = By_1 + C_X + eps;
AA = [1 -a2 0
     -b2 1 -b3
     0   0    1];

BB = [a3 0  0
     0  0 b4
     0  1  0];

CC = [a1 a4 0
      b1 0 b5
      0  0 0];
 
Y = zeros(n,1);
Y_1 = zeros(n,1);
Y_2 = zeros(n,1);
C = zeros(n,1);
C_1 = zeros(n,1);


for i=1:n
    %Zajisteni pocatecnich podminek
    if i > 1
        Y_1(i)=Y(i-1);
        C_1(i)=C(i-1);
    end
    if i>2
        Y_2(i)=Y(i-2);
    end
    
    %simulace systemu
    YY = inv(AA)*BB*[C_1(i); Y_1(i); Y_2(i)] + inv(AA)*CC*[1; G(i); R(i)] + inv(AA)*[eps_C(i); eps_Y(i); 0];
    C(i) =YY(1);
    Y(i) =YY(2);
    Y_1(i) =YY(3); %neni nutno pocitat
end

C_2SLS = C;
Y_2SLS = Y;

%Overeni stability dle vlastnich cisel prechodove matice Delta ("dynamicke
%jadro?"
delta = inv(AA)*BB;
lambda = eig(delta);
disp('Lambda - 2SLS')
disp(lambda)
abs_lambda = abs(lambda);
disp('Abs. lambda - 2SLS')
disp(abs_lambda)



figure
subplot(2,1,1)
plot([C_model C_OLS C_2SLS])
title('Simulace C')
legend('C','C\_OLS','C\_2SLS')
subplot(2,1,2)
plot([Y_model Y_OLS Y_2SLS])
title('Simulace Y')
legend('Y','Y\_OLS','Y\_2SLS')