%% Problem Set 6, question 3 (Michele Piffer, LSE, 30/11/2011)
% Dummy Variables

clc 
clear                 % get rid of data already in matlab

%% Generate Regressors

n = 100;               % number of rows
k = 2;                 % number of regressors
X = 10*rand(n,k);      % random values assigned

epsilon = randn(n,1);  % true errors

n1 = 20;               % number of people in each group
n2 = 30;
n3 = 10;
n4 = n - n1 - n2 - n3;

d1 = ones(n1,1);
d2 = ones(n2,1);
d3 = ones(n3,1);
d4 = ones(n4,1);

D = [d1 zeros(n1,3); zeros(n2,1) d2 zeros(n2,2); zeros(n3,2) d3 zeros(n3,1);
    zeros(n4,3) d4];

W = [D X]; 
k = size(W,2);         %  adjust for the increase in the number of regressors

% Note, I assume so far that the DGP does not have a constant

%% Generating the DGP for y

beta_true = [5 8 4 -2 16 3]'; % true parameters, which you want to estimate

y = W*beta_true + epsilon; 

%% Estimate the model regressing without the intercept

theta_hat = (W'*W)^(-1)*W'*y;

%to compute variance covariance
MW = eye(n) - W*(W'*W)^(-1)*W'; 
RSS = y'*MW*y;
Ssquare = RSS/(n-k);
Vtheta = Ssquare*(W'*W)^(-1); 

%% Estimate the model regressing with the intercept and dropping the first
% dummy

WW = W(:, 2:k);              % drop the first dummy
WW = [ones(n,1) WW];         % add a constant    

beta_hat = (WW'*WW)^(-1)*WW'*y;

%to compute variance covariance
MWW = eye(n) - WW*(WW'*WW)^(-1)*WW'; 
RSS = y'*MWW*y;
Ssquare = RSS/(n-k);
Vbeta = Ssquare*(WW'*WW)^(-1); 

%% Compare what you get

beta_true'
display('thetas')
theta_hat'
display('betas')
beta_hat'

Vtheta - Vbeta              % makes sense

%% Test if group 1 or 2 have the same effect
% compute statistics and compare to critical value 1.66

t_beta = beta_hat(2)/Vbeta(2,2)       

t_theta = (theta_hat(2)-theta_hat(1))/(Vtheta(1,1)+Vtheta(2,2)-2*Vtheta(1,2)) % the same

%% Play with a constant

% What we cannot do is use both all dummies and add a constant. But this
% does not mean that the DGP did not generate the data with a constant!
% Assume for instance that DGP has a constant, and that coefficient on the
% constant is 1. Can we estimate the coefficient on the constant? No, the
% effect is taken by the dummies, which now include both the effect of
% belonging to the group and the effect of the constant

ybis = ones(n,1) + W*beta_true + epsilon;       % new y variable

% Estimate the model without the constant and with all dummies

theta_hatbis = (W'*W)^(-1)*W'*ybis;              

% Estimate the model dropping one dummy and adding the constant

beta_hatbis = (WW'*WW)^(-1)*WW'*ybis;

display('thetas without a constant in the DGP')
theta_hat' 

display('thetas with a constant in the DGP')
theta_hatbis'     % the effect of the constant is take in all dummies

display('betas without a constant in the DGP')
beta_hat'

display('betas with a constant in the DGP')
beta_hatbis'      % the effect of the constant is taken in beta0

