%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Econometrics 2
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This code must be modified to answer questions 8.10-11 of
% Davidson-Mackinnon.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
clc

% Set the true parameters and placeholders for results.

N = 10;                                                                     % Number of observations.
sim = 1000;                                                                 % Number of simulations.

beta = 0;                                                                   % Coefficient on X.
pi = 1;                                                                     % Coefficient on W.

sigmau = 1;                                                                 % Standard deviation for u.
sigmav = 1;                                                                 % Standard deviation for v.
rho = 0.5;                                                                  % Correlation between u and v.

bhat = zeros(sim,1);                                                        % Placeholder for estimates.
se = zeros(sim,1);                                                          % Placeholder or standard errors.

% Begin simulation.

for s = 1:sim 

    % Generate the data. Modify the code so that the data generating 
    % process is exactly as in (8.40) of the book. Rescaling of the 
    % instrument and plotting of empirical and theoretical distributions 
    % have already been done for you.
  
    e = mvnrnd([0; 0],[sigmau^2 0; 0 sigmav^2],N);                          % Variance-covariance matrix of error terms.
    
    u = e(:,1);                                                             % Error terms for y.                                                     
    v = e(:,2);                                                             % Error terms for x.
                   
    w = normrnd(0,1,N,1);                                                   % Exogenous instrument.
    w = w*sqrt(N/(sum(w.*w)));                                              % Rescale w so that w'*w = N.
    
    x = normrnd(0,1,N,1);                                                   % Regressor.
    y = beta*x + u;            
    
    % OLS estimates.
    
    bhat(s) = inv(x'*x)*x'*y;                                               % Estimate of beta.
    
    uhat = y - bhat(s)*x;
    s2 = (uhat'*uhat)/N-1;
    
    vmat = s2*inv(x'*x);
    se(s) = sqrt(vmat);                                                     % Standard error.
    
end

% Display the results.

fprintf('--------------------------------------------------------------\n')
fprintf('Results of Last Simulation \n')
fprintf('--------------------------------------------------------------\n')
fprintf('\n')

fprintf('N = %d \n', N)
fprintf('Sim = %d \n\n', sim)

fprintf('             b            SE   \n')
fprintf('beta      %0.4f       %0.4f     \n', bhat(sim), se(sim))
fprintf('\n')

fprintf('--------------------------------------------------------------\n')
fprintf('Empirical Results \n')
fprintf('--------------------------------------------------------------\n')
fprintf('\n')

fprintf('           Mean        Std.Dev   \n')
fprintf('bhat      %0.4f       %0.4f     \n', mean(bhat(:)), std(bhat(:)))

clf

[cdf1, b1] = ecdf(bhat);                                                    % Empirical CDF of OLS estimate.
plot(b1,cdf1);

hold on

[cdf2, b2] = ecdf(bhat_iv);                                                 % Empirical CDF of IV estimate.
plot(b2,cdf2);

b3 = -2:0.001:2;
cdf3 = evcdf(b3,0,(sigmau^2/(N*pi^2)));                                     % Theoretical CDF.
plot(b3,cdf3,'m');

title('CDF of $$\hat{b}$$','Interpreter','Latex')
legend('Empirical (OLS)','Empirical (IV)','Theoretical','location','NW');
xlabel('$$\hat{b}$$','Interpreter','Latex');
ylabel('$$F(\hat{b}$$)','Interpreter','Latex');

hold off