clc
clear all
close all

%{
vlastnosti odhadu OLS pro dvourovnicovy model model
    y1 = a1 + a2*x1 + e
    y2 = b1*y1 + b2 + b3*x2 + e
-konzistence
-nestrannost
-asymptoticka konzistence
%}

a1 = 1.0; a2 = 1.0;
b1 = 1.0; b2 = 1.0; b3 = 1.0;


disp('KONZISTENCE');
maxnobs = 1000;
aOLS = zeros(2,maxnobs-10);
bOLS = zeros(3,maxnobs-10);
disp('probihaji regrese pro konzistenci!');
tic;
for i=10:maxnobs
    x1 = randn(i,1); x2 = randn(i,1);
    iota = ones(i,1);
    y1 = zeros(i,1); y2 = zeros(i,1); evec = randn(i,1);
        for j = 1:i
            y1(j,1) = iota(j,1)*a1 + x1(j,1)*a2 + evec(j,1);
            y2(j,1) = iota(j,1)*b1 + y1(j,1)*b2 + x2(j,1)*b3 +evec(j,1);
        end;
    resulta = ols(y1,[iota x1]);
    resultb = ols(y2,[iota y1 x2]);
    aOLS(:,i-9) = resulta.beta;
    bOLS(:,i-9) = resultb.beta;
end;
fprintf('\n Celkovy cas vypoctu pro konzistenci %6.2f sekund \n \n',toc);

figure
    subplot(2,1,1)
        plot(aOLS(1,:));
        title('a1');
    subplot(2,1,2)
        plot(aOLS(2,:));
        title('a2');
figure
    subplot(3,1,1)
        plot(bOLS(1,:));
        title('b1');
    subplot(3,1,2)
        plot(bOLS(2,:));
        title('b2');
    subplot(3,1,3)
        plot(bOLS(3,:));
        title('b3');

disp('NESTRANNOST');
nobs = 50;
niter = 1000;
disp('probihaji regrese pro nestrannost!');
tic;
aOLS = zeros(2,niter);
bOLS = zeros(3,niter);
x1 = randn(nobs,1); x2 = randn(nobs,1);
iota = ones(nobs,1);
for iter = 1:niter;
        y1 = zeros(nobs,1); y2 = zeros(nobs,1); evec = randn(nobs,1);
        for j = 1:nobs
            y1(j,1) = iota(j,1)*a1 + x1(j,1)*a2 + evec(j,1);
            y2(j,1) = iota(j,1)*b1 + y1(j,1)*b2 + x2(j,1)*b3 +evec(j,1);
        end;
        resulta = ols(y1,[iota x1]);
        resultb = ols(y2,[iota y1 x2]);
        aOLS(:,iter) = resulta.beta;
        bOLS(:,iter) = resultb.beta;
end;
aOLSm = mean(aOLS,2);
bOLSm = mean(bOLS,2);

fprintf('\n Celkovy cas vypoctu pro nestrannost %6.2f sekund \n \n',toc);

fprintf(['OLS vysledky prvni rovnice pres ', num2str(niter),' iteraci \n']);
in.rnames = strvcat('Koeficient','a1','a2');
in.cnames = strvcat('mean');
mprint([aOLSm],in);
fprintf(['OLS vysledky druhe rovnice pres ', num2str(niter),' iteraci \n']);
in.rnames = strvcat('Koeficient','b1','b2','b3');
mprint([bOLSm],in);

disp('ASYMPTOTICKA NESTRANNOST')
niter = 50;
maxnobs = 500;
aOLS = zeros(2,niter);
bOLS = zeros(3,niter);
aOLSm=zeros(2,maxnobs-10);
bOLSm=zeros(3,maxnobs-10);
disp('probihaji regrese pro asymptotickou nestrannost!');
tic;
for i=10:maxnobs
    x1 = randn(i,1); x2 = randn(i,1);
    iota = ones(i,1);
    for iter = 1:niter;
        y1 = zeros(i,1); y2 = zeros(i,1); evec = randn(i,1);
        for j = 1:i
            y1(j,1) = iota(j,1)*a1 + x1(j,1)*a2 + evec(j,1);
            y2(j,1) = iota(j,1)*b1 + y1(j,1)*b2 + x2(j,1)*b3 +evec(j,1);
        end;
        resulta = ols(y1,[iota x1]);
        resultb = ols(y2,[iota y1 x2]);
        aOLS(:,iter) = resulta.beta;
        bOLS(:,iter) = resultb.beta;
    end;
    aOLSm(:,i-9) = mean(aOLS,2);
    bOLSm(:,i-9) = mean(bOLS,2);
end;
fprintf('\n Celkovy cas vypoctu pro asymptotickou nestrannost %6.2f sekund \n \n',toc);

figure
    subplot(2,1,1)
        plot(aOLSm(1,:));
        title('mean a1');
    subplot(2,1,2)
        plot(aOLSm(2,:));
        title('mean a2');
figure
    subplot(3,1,1)
        plot(bOLSm(1,:));
        title('mean b1');
    subplot(3,1,2)
        plot(bOLSm(2,:));
        title('mean b2');
    subplot(3,1,3)
        plot(bOLSm(3,:));
        title('mean b3');