%% Problem Set 6, question 2a (Michele Piffer, LSE, 22/11/2011)
% Test if beta2 - 3beta3 = 4 AND beta1 = 2beta4

clc 
clear % get rid of data already in matlab

%% Generate DGP

n = 100;               % number of rows
k = 3;                 % number of regressors
X = 2*rand(n,k);       % random values assigned
X = [ones(n,1), X];    % add a constant term

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

beta_true = [8 13 3 4]'; % true parameters, which you want to estimate

y = X*beta_true + epsilon; 

if beta_true(2) - 3*beta_true(3) == 4 && beta_true(1) == 2*beta_true(4)
    display('H0 true')
else
    display('H0 false')
end

beta_true(2) - 3*beta_true(3)  
beta_true(1) - 2*beta_true(4)

%% Estimate model using X and y

display('True parameters and estimates')

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

display(beta_true')
display(beta_hat')    

MX = eye(n) - X*(X'*X)^(-1)*X'; 
% ones(1,n)*MX*y        %check that the sum of residuals is zero
RSS = y'*MX*y;
Ssquare = RSS/(n-(k+1));
V = Ssquare*(X'*X)^(-1);  % variance covariance of estimators                
                                      
%% Test your hypothesis using Wald Principle

R = [0, 1, -3, 0; 1, 0, 0, -2];
q = [4 0]';

WaldStatistic = 0.5*(R*beta_hat - q)'*(R*V*R')^(-1)*(R*beta_hat - q)
critical_value = 3.08;  %F(2,n-k) at 5 percent 

if WaldStatistic > critical_value
    display('reject H0')
else
    display('fail to reject H0')
end
                  
%% Test your hypothesis using LR Principle

yrest = y - 4*X(:,2);
x1rest = ones(n,1) + 0.5*X(:,4);
x2rest = 3*X(:,2) + X(:,3);
Xrest = [x1rest x2rest];

MXrest = eye(n) - Xrest*(Xrest'*Xrest)^(-1)*Xrest'; 

RSSrest = yrest'*MXrest*yrest;

Fstatistic = (RSSrest - RSS)*(n-(k+1))/(2*RSS); % exactly the same statistic
    

%% Simulate the model p times

display('Simulate your hypothesis testing')

if beta_true(2) - 3*beta_true(3) == 4 && beta_true(1) == 2*beta_true(4)
    display('H0 true')
else
    display('H0 false')
end

p = 10000;         % number of simulations
Reject = 0;        % number of rejections of H0

for i=1:p;
    
X = 2*rand(n,k);       % random values assigned
X = [ones(n,1), X];    % add a constant term
epsilon = randn(n,1);  % true errors
y = X*beta_true + epsilon;

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

MX = eye(n) - X*(X'*X)^(-1)*X'; 
RSs = y'*MX*y;
Ssquare = RSS/(n-(k+1));
V = Ssquare*(X'*X)^(-1);  % variance covariance of estimators                
                                      
WaldStatistic = 0.5*(R*beta_hat - q)'*(R*V*R')^(-1)*(R*beta_hat - q);

if WaldStatistic > critical_value
    Reject = Reject + 1;
end

end

rejection_ratio = Reject/p;

rejection_ratio 

% Change beta_true so that H0 not true any more and see what happens

