clear all
close all
clc

n = 1000;
%
% exogenni promenne: mozno generovat libovolně(veřejné výdaje,úroková míra)
%
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

% skutečné hodnoty regresních koeficientů
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;

% maticový zápis strukturního tvaru modelu
% Ay = By_1 + C_X + eps; 
%
% Matice koeficientů mezi běžnými endogenními proměnnými
%
AA = [1 -a2 0
     -b2 1 -b3
     0   0   1];
%
% Matice koeficientů mezi zpožděnými endogenními proměnnými
%
BB = [a3 0  0
     0  0 b4
     0  1  0];
%
% Matice koeficientů mezi exogenními proměnnými
%
CC = [a1 a4 0
      b1 0 b5
      0  0 0];
% 
% nulování proměnných modelu
%
Y = zeros(n,1);
Y_1 = zeros(n,1);
Y_2 = zeros(n,1);
C = zeros(n,1);
C_1 = zeros(n,1);
%
% generování zpožděných proměnných
%
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
%
%   Redukovaná forma modelu
%   
    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 parametrů
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);
%
% odhad 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 modelu s parametry OLS
% dosazení parametrů a1,a2,a3,a4 resp. b1,b2,b3,b4,b5 do obou rovnic 
%
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);
%
% Přepočet modelu pomocí aktualizovaných hodnot parametrů
% 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);
%
% generování zpožděných proměnných
%
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
%
% uložení trajektorií C a Y vypočtených dle OLS
%
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 přes redukovanou formu modelu s parametry 2SLS
%
   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
%
% uložení trajektorií C a Y vypočtených dle 2SLS
%
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')