%% [0] LE MODELE DE LA NATURE & CONSTRUCTION D'UN ECHANTILLON.

x = rand(50,1);% [ .13 ; .24 ; .33; .45 ; .54 ; .34 ; .12 ]; % On se donne une variable explicative.
X = [x , ones(size(x)) ];                       % On forme la matrice des "variables explicatives" en rajoutant la constante.

a = .1;                                         % La vraie valeur du coefficient associé à x.  
b = .5;                                         % La vraie valeur de la constante.
Coefficients = [a;b];                           % Le vecteur (colonne) des coefficients.

y = X*Coefficients + randn(size(x));            % On joue à Dieu en créant un échantillon (à l'aide de la fonction randn qui
                                                % génère des nombres aléatoirement tités dans une gaussienne).


%% [1] ESTIMATION DU MODELE DE L'ECONOMETRE.

%% [1.1] On code bêtment les formules de l'estimateur des MCO dans le modèle simple (vu en cours) en utilisant des
%%       boucles.

ybar = mean(y);                                 % La moyenne de y.
xbar = mean(x);                                 % La moyenne de x.
n = size(x,1);                                  % Le nombre d'observations dans l'échantillon.

tic;                                            % Pour mesurer le temps de calcul (i).

%% Calcul du dénominateur apparaissant dans l'expression de l'estimateur des MCO du paramètre a.
den = 0;
for i=1:n
    den = den + (x(i)-xbar)^2;
end

%% Calcul du numérateur apparaissant dans l'expression de l'estimateur des MCO du paramètre a.
num = 0;
for i=1:n
    num = num + (y(i)-ybar)*(x(i)-xbar);
end

%% Finalement l'estimateur du paramètre associé à x est:
ahat = num/den;
%% ... Et l'estimateur de la constante:
bhat = ybar-xbar*ahat;

tElapsed = toc(tStart);                          % Pour mesurer le temps de calcul (ii)

%% Affichons les résultats:
disp(['ahat = ' num2str(ahat)])
disp(['bhat = ' num2str(bhat)])
disp(['Elapsed time is ' num2str(tElapsed)])


%% [1.2] On code bêtment les formules de l'estimateur des MCO dans le modèle simple (vu en cours) sans les boucles.

ybar = mean(y);                                 % La moyenne de y.
xbar = mean(x);                                 % La moyenne de x.
n = size(x,1);                                  % Le nombre d'observations dans l'échantillon.

tic;                                            % Pour mesurer le temps de calcul (i).

%% Calcul du dénominateur apparaissant dans l'expression de l'estimateur des MCO du paramètre a.
xc  = (x-xbar); 
den = transpose(xc)*(xc);

%% Calcul du numérateur apparaissant dans l'expression de l'estimateur des MCO du paramètre a.
num = transpose(xc)*(y-ybar);

%% Finalement l'estimateur du paramètre associé à x est:
ahat = num/den;
%% ... Et l'estimateur de la constante:
bhat = ybar-xbar*ahat;

tElapsed = toc(tStart);                          % Pour mesurer le temps de calcul (ii)

%% Affichons les résultats:
disp(['ahat = ' num2str(ahat)])
disp(['bhat = ' num2str(bhat)])
disp(['Elapsed time is ' num2str(tElapsed)])