%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Homework 6 - Solution
% Econometrics 2
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


clear
clc

% Set the true parameters and placeholders for results.

T = 150;                                                                   % Number of observations.

beta0 = 0;                                                                  % Intercept in equation 1.
beta1 = 1.5;                                                                % Coefficient on Y2 in equation 1.
beta2 = 0.3;                                                                % Coefficient on X1 in equation 1.
beta3 = 0.0;                                                                % Intercept in equation 2.
beta4 = 0.3;                                                                % Coefficient on x1 in equation 2.
beta5 = 2;                                                                % Coefficient on x2 in equation 2.
beta6 = 0.015;                                                               % Coefficient on x3 in equation 2.
beta7 = 0.0;                                                                % Coefficient on x4 in equation 2
beta8 = 0.0;
beta9 = 0.0;
beta10 = 0.0;



sigma1 = 1;                                                                 % Standard deviation of u1.
sigma2 = 1;                                                                 % Standard deviation of u2.

sim = 200;                                                                 % Number of simulations.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The simultaneous equations model is
%           Y1 = beta0 + beta1*Y2 + beta2*x1 + u1
%           Y2 = beta3 + beta4*x1 + beta5*x2 + beta6*x3 + beta7*x4 + beta8*x5 + beta9*x6 + beta10*x7+ u2
%
% This code estimates the coefficients in the first equation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Generate the data.
x1 = 2 + normrnd(0,1,T,1);                                    % x1.
x2 = 3 + sigma1.*normrnd(0,1,T,1);                                 % x2.
x3 = -2 + 0.4*x2 + sigma1.*normrnd(0,1,T,1);
x4 = normrnd(0,1,T,1);
x5 =  0.9*sigma2.*normrnd(0,1,T,1);
x6 = 1 + sigma1.*normrnd(0,1,T,1);
x7 =  0.3*sigma1.*normrnd(0,1,T,1);
x8 =  0.9*sigma2.*normrnd(0,1,T,1);
x9 = 1 + sigma1.*normrnd(0,1,T,1);
x10 =  0.3*sigma1.*normrnd(0,1,T,1);
x11 =  0.9*sigma2.*normrnd(0,1,T,1);
x12 = 1 + sigma1.*normrnd(0,1,T,1);
x13 =  0.3*sigma1.*normrnd(0,1,T,1);

X = [ones(T,1) x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13];                                        % X.
 

B_2sls = zeros(sim, 3);                                             % 2SLS results.
B_liml = zeros(sim, 3);                                             %liml results.
B_ols = zeros(sim, 3);                                              % OLS results.
B_Fuller = zeros(sim, 3);                                           % Fuller results

for s = 1:sim
    
    u1 = normrnd(0,sigma1,T,1);                                             % Residuls for equation 1.
    u2 = normrnd(0,sigma2,T,1) + 5*u1;                                      % Residuals for equation 2.    
    
    %Generate y2
    y2 = beta3 + beta4*x1 + beta5*x2 + beta6*x3 + beta7*x4 + beta8*x5 + beta9*x6 + beta10*x7 + u2;            
    
    y1 = beta0 + beta1*y2 + beta2*x1 + u1;
    
    Y = [y1 y2];  
    Y1 = Y(:,1);   %same as y1                      
    Y2 = Y(:,2);   %endogenous regressors, same as generated y2
    
     
    X_exo1 = X(:,1:2); %exogenous regressors in the reduced form
    X_OLS=[ones(T,1) Y2 x1];
    
    % Two-stage Least Squares: Stage 1
    B2_hat = inv(X'*X)*X'*Y2;                                               % First stage estimates.
    Y2_hat = X*B2_hat;                                                      % Fitted values of Y2.

    % Two-stage Least Squares: Stage 2
    X1_hat = [ones(T,1) Y2_hat x1];       
    B_2sls(s,:) = inv(X1_hat'*X1_hat)*X1_hat'*Y1;                           % 2SLS estimates.
    
    %LIML estimation(taken form http://web.mit.edu/nagaraj/files/liml.m)
    N = length(Y2);
    Mexo = eye(N) - X*inv(X'*X)*X';                                    %projection matrix, X is the instrument for X1   
    Mexo1 = eye(N) - X_exo1*inv(X_exo1'*X_exo1)*X_exo1';
    W = [Y1 Y2]'*(Mexo)*[Y1 Y2];                                       %2x2 matrix
    W1 = [Y1 Y2]'*(Mexo1)*[Y1 Y2];   %2x2 matrix     
    lambda = min( eig(inv(W)*W1 )) ;  
    lambda_fuller = min(eig(inv(W)*W1))-(1/(T-(size(X,2)-size(X_exo1,2))));                                           %finds the min kappa from the book
    %finds the min kappa from the book
    B_liml(s,:) = inv(X_OLS'*(eye(N)-(lambda*Mexo))*X_OLS)*(X_OLS'*(eye(N)-(lambda*Mexo))*Y1); 
    B_Fuller(s,:) = inv(X_OLS'*(eye(N)-(lambda_fuller*Mexo))*X_OLS)*(X_OLS'*(eye(N)-(lambda_fuller*Mexo))*Y1); 
    
    
    % OLS esimation
    
    X1_ols = [ones(T,1) Y2 x1];    
    B_ols(s,:) = inv(X1_ols'*X1_ols)*X1_ols'*Y1;                            % OLS estimates.

end

% Display the results of the last estimation.

fprintf('--------------------------------------------------------------\n')
fprintf('Empirical Results: 2SLS estimates for Equation 1\n')
fprintf('--------------------------------------------------------------\n\n')

fprintf('         Mean   Median    Std.Dev   \n')
fprintf('b0      %0.4f   %0.4f      %0.4f     \n', mean(B_2sls(:,1)), median(B_2sls(:,1)), std(B_2sls(:,1)))
fprintf('b1      %0.4f   %0.4f      %0.4f     \n', mean(B_2sls(:,2)), median(B_2sls(:,2)), std(B_2sls(:,2)))
fprintf('b2      %0.4f   %0.4f      %0.4f     \n',mean(B_2sls(:,3)), median(B_2sls(:,3)), std(B_2sls(:,3)))

fprintf('--------------------------------------------------------------\n')
fprintf('Empirical Results: LIML estimates for Equation 1\n')
fprintf('--------------------------------------------------------------\n\n')

fprintf('         Mean   Median    Std.Dev   \n')
fprintf('b0      %0.4f   %0.4f      %0.4f     \n', mean(B_liml(:,1)), median(B_liml(:,1)), std(B_liml(:,1)))
fprintf('b1      %0.4f   %0.4f      %0.4f     \n', mean(B_liml(:,2)), median(B_liml(:,2)), std(B_liml(:,2)))
fprintf('b2      %0.4f   %0.4f      %0.4f     \n',mean(B_liml(:,3)), median(B_liml(:,3)), std(B_liml(:,3)))

fprintf('--------------------------------------------------------------\n')
fprintf('Empirical Results: Fuller estimates for Equation 1\n')
fprintf('--------------------------------------------------------------\n\n')

fprintf('         Mean   Median    Std.Dev   \n')
fprintf('b0      %0.4f   %0.4f      %0.4f     \n', mean(B_Fuller(:,1)), median(B_Fuller(:,1)), std(B_Fuller(:,1)))
fprintf('b1      %0.4f   %0.4f      %0.4f     \n', mean(B_Fuller(:,2)), median(B_Fuller(:,2)), std(B_Fuller(:,2)))
fprintf('b2      %0.4f   %0.4f      %0.4f     \n',mean(B_Fuller(:,3)), median(B_Fuller(:,3)), std(B_Fuller(:,3)))

fprintf('--------------------------------------------------------------\n')
fprintf('Empirical Results: OLS estimates for Equation 1\n')
fprintf('--------------------------------------------------------------\n\n')

fprintf('         Mean   Median    Std.Dev   \n')
fprintf('b0      %0.4f   %0.4f      %0.4f     \n', mean(B_ols(:,1)), median(B_ols(:,1)), std(B_ols(:,1)))
fprintf('b1      %0.4f   %0.4f      %0.4f     \n', mean(B_ols(:,2)), median(B_ols(:,2)), std(B_ols(:,2)))
fprintf('b2      %0.4f   %0.4f      %0.4f     \n', mean(B_ols(:,3)), median(B_ols(:,3)), std(B_ols(:,3)))

