%% IV estimation (Michele Piffer, LSE, 11/2/2012)

clc 
clear % get rid of data already in matlab

%% Generate DGP

n = 5000;                             % number of rows
x1 = ones(n,1);

z = randn(n,1);

epsilon = randi([-15 15], n, 1);     % zero in expected value
vareps = sum(epsilon.^2)/n-(sum(epsilon)/n)^2;

rho1 = 100;                            % cov(x,epsilon)
rho2 = 60;                            % cov(x,z)
meanx2 = 5;
x2 = (rho1/vareps)*epsilon + 10*randn(n,1) + meanx2 + rho2*z;    

% Sample covariances
display('Endogeneity of the regressor')
epsilon'*x2/n - (sum(epsilon)/n)*(sum(x2)/n)   % yes, approximately equals rho

display('Validity of the instrument')
epsilon'*z/n - (sum(epsilon)/n)*(sum(z)/n) 

display('Relevance of the instrument')
x2'*z/n - (sum(x2)/n)*(sum(z)/n) 

X = [x1 x2];     
Z = [x1 z];

beta_true = [2 5]';                  % true parameters

y = X*beta_true + epsilon; 

%% Estimate model using X and y

beta_hat = (X'*X)^(-1)*X'*y;         % OLS regressors
beta_hat_iv = (Z'*X)^(-1)*Z'*y;         % IV regressors

display('true')
beta_true'

display('Direction of the bias')
beta_hat'

display('IV estimation')
beta_hat_iv'
