%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Xavier Martin G. Bautista
% Econometrics 2
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The following performs likelihood ratio, Wald, Lagrange multiplier and 
% tests. 
%                       The null is H0: beta2 = 0.
%                 Th alternative is HA: beta2 != 0.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
clc

global y x1 x2 T

% Set the true parameters and placeholders for results.

T = 1000;                                                                   % Sample size.

beta0 = 0.5;                                                                % Intercept.
beta1 = 3;                                                                  % Coefficient on X1.
beta2 = 0;                                                                  % Coefficient on X2.

sigma = 1;                                                                  % Standard deviation.

sim = 1000;                                                                 % Number of simulations.

b_ols = zeros(sim,3);                                                       % OLS estimates.

t = zeros(sim,1);                                                           % t-stats.
pval_t = zeros(sim,1);                                                      % p-values.
reject_t = zeros(sim,1);                                                    % Rejections.

LR = zeros(sim,1);                                                          % Likelihood ratio test statistic.
pval_LR = zeros(sim,1);                                                     % p-values.
reject_LR = zeros(sim,1);                                                   % Rejections.

Wald = zeros(sim,1);                                                        % Wald statistic.
pval_Wald = zeros(sim,1);                                                   % p-values.
reject_Wald = zeros(sim,1);                                                 % Rejections.

LM = zeros(sim,1);                                                          % Lagrange multiplier test statistic.
pval_LM = zeros(sim,1);                                                     % p-values.
reject_LM = zeros(sim,1);                                                   % Rejections.

for s = 1:sim 
    
    % Generate the data.

    x1 = ((1:1:T)'/T).*normrnd(0,1,T,1);
    x2 = ones(T,1) + 2.*x1 + 5.*normrnd(0,1,T,1);
    x = [ones(T,1) x1 x2];

    u = sigma.*randn(T,1);

    y = beta0 + beta1.*x1 + beta2.*x2 + u;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % OLS estimates.
    
    bhat = inv(x'*x)*x'*y;                                                  % OLS estimates.
    b_ols(s,:) = bhat;                                                     
    
    % t-stats and p-vlaues of beta2.
    
    uhat = y - x*bhat;                                                      % Residuals.
    s2 = (uhat'*uhat)/(T-3);                                                % s squared.

    vmat_ols = inv(x'*x).*s2;                                               % Variance matrix,
    t(s) = bhat(3)/sqrt(vmat_ols(3,3));                                     % t-stat of bhat2.
    pval_t(s) = 2*tcdf(abs(t(s)),T - 3,'upper');                            % p-value.

    if pval_t(s) >= 0.05                                                    % Rejection decision; 5% significance. 0 = fail to reject, 1 = reject.
        reject_t(s) = 0;
    else
        reject_t(s) = 1;
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % Likelihood Ratio Test.

    options = optimset('Display','off');
    
    %( MISSING CODE %}

    if pval_LR(s) >= 0.05                                                   % Rejection decision; 5% significance. 0 = fail to reject, 1 = reject.
        reject_LR(s) = 0;
    else
        reject_LR(s) = 1;
    end

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % Wald test.

    vmat_w = inv(hess);                                                     % Empirical hessian.

    Wald(s) = ((b1_mle(3))^2)/vmat_w(3,3);                                  % Wald statistic.
    pval_Wald(s) = chi2cdf(Wald(s),1,'upper');                              % p-value.

    if pval_Wald(s) >= 0.05                                                 % Rejection decision; 5% significance. 0 = fail to reject, 1 = reject.
        reject_Wald(s) = 0;
    else
        reject_Wald(s) = 1;
    end

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % Lagrange Multiplier Test.
    
    grad_LM = zeros(4,1);                                                   % Gradient of unrestricted loglikelihood evaluated at restricted ML estimates.
    grad_LM(1,1) = sum((y - ones(size(y,1),1).*b2_mle(1) - ...
                                                b2_mle(2)*x1)/b2_mle(3)^2);
    grad_LM(2,1) = sum(x1.*(y - ones(size(y,1),1).*b2_mle(1) - ...
                                               b2_mle(2)*x1)./b2_mle(3)^2);
    grad_LM(3,1) = sum(x2.*(y - ones(size(y,1),1).*b2_mle(1) - ...
                                                b2_mle(2)*x1)/b2_mle(3)^2);
    grad_LM(4,1) = sum(ones(size(y,1),1).*(-1/b2_mle(3)) + ...
                            (y - ones(size(y,1),1).*b2_mle(1) - ...
                                             b2_mle(2)*x1).^2/b2_mle(3)^3);

    vmat_IM = [b2_mle(3)*inv(x'*x) zeros(3,1);                              % IM estimator evaluated at restricted ML estimates.
               zeros(1,3)    b2_mle(3)/(2*T)];

    LM(s) = grad_LM'*vmat_IM*grad_LM;                                       % LM statistic.
    pval_LM(s) = chi2cdf(LM(s),1,'upper');                                  % p-values.

    if pval_LM(s) >= 0.05                                                   % Rejection decision; 5% significance. 0 = fail to reject, 1 = reject.
        reject_LM(s) = 0;
    else
        reject_LM(s) = 1;
    end

end

% Display the results of the last estimation.

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

fprintf('Unrestricted ML Estimates \n')
fprintf('\n')
fprintf('            b         \n')
fprintf('b0        %0.4f       \n', b1_mle(1))
fprintf('b1        %0.4f       \n', b1_mle(2))
fprintf('b2        %0.4f       \n', b1_mle(3))
fprintf('sigma     %0.4f       \n', b1_mle(4))
fprintf('\n')

fprintf('Restricted ML Estimates \n')
fprintf('Restriction: beta2 = 0\n')
fprintf('\n')
fprintf('            b         \n')
fprintf('b0        %0.4f       \n', b2_mle(1))
fprintf('b1        %0.4f       \n', b2_mle(2))
fprintf('sigma     %0.4f       \n', b2_mle(3))
fprintf('\n')

fprintf('--------------------------------------------------------------\n')
fprintf('t-test. \n')
fprintf('--------------------------------------------------------------\n')
fprintf('\n')

fprintf('t-statistic \n')
fprintf('\n')
fprintf('The most recent t-statistic is:   %0.4f \n', t(sim))
fprintf('\n')
fprintf('p-value \n')
fprintf('\n')
fprintf('The most recent p-value is:   %0.4f \n', pval_t(sim))

if pval_t(sim)>=0.05
    fprintf('The p-value is at least 0.05 so the t-test does not reject the null that beta2 is zero. \n')
    fprintf('\n')
else
    fprintf('The p-value is less than 0.05 so the t-test rejects the null that beta2 is zero.')
    fprintf('\n')
end
fprintf('\n')

fprintf('--------------------------------------------------------------\n')
fprintf('Empirical size of t-test. \n')
fprintf('--------------------------------------------------------------\n')
fprintf('\n')

%{ MISSING CODE %}

fprintf('--------------------------------------------------------------\n')
fprintf('Likelihood Ratio Test. \n')
fprintf('--------------------------------------------------------------\n')
fprintf('\n')

fprintf('LR Statistic \n')
fprintf('\n')
fprintf('The most recent LR statistic is:   %0.4f \n', LR(sim))
fprintf('\n')
fprintf('p-value \n')
fprintf('\n')
fprintf('The most recent p-value is:   %0.4f \n', pval_LR(sim))

if pval_LR(sim)>=0.05
    fprintf('The p-value is at least 0.05 so the LR test does not reject the null that beta2 is zero. \n')
    fprintf('\n')
else
    fprintf('The p-value is less than 0.05 so the LR test rejects the null that beta2 is zero.')
    fprintf('\n')
end
fprintf('\n')

fprintf('--------------------------------------------------------------\n')
fprintf('Empirical size of Likelihood Ratio Test. \n')
fprintf('--------------------------------------------------------------\n')
fprintf('\n')

%{ MISSING CODE %}

fprintf('--------------------------------------------------------------\n')
fprintf('Wald Test. \n')
fprintf('--------------------------------------------------------------\n')
fprintf('\n')

fprintf('Wald Statistic \n')
fprintf('\n')
fprintf('The most recent Wald statistic is:   %0.4f \n', Wald(sim))
fprintf('\n')
fprintf('p-value \n')
fprintf('\n')
fprintf('The most recent p-value is:   %0.4f \n', pval_Wald(sim))

if pval_Wald(sim)>=0.05
    fprintf('The p-value is at least 0.05 so the Wald test does not reject the null that beta2 is zero. \n')
    fprintf('\n')
else
    fprintf('The p-value is less than 0.05 so the Wald test rejects the null that beta2 is zero.')
    fprintf('\n')
end
fprintf('\n')

fprintf('--------------------------------------------------------------\n')
fprintf('Empirical size of Wald Test. \n')
fprintf('--------------------------------------------------------------\n')
fprintf('\n')

%{ MISSING CODE %}

fprintf('--------------------------------------------------------------\n')
fprintf('Lagrange Multiplier Test. \n')
fprintf('--------------------------------------------------------------\n')
fprintf('\n')

fprintf('LM Statistic \n')
fprintf('\n')
fprintf('The most recent LM statistic is:   %0.4f \n', LM(sim))
fprintf('\n')
fprintf('p-value \n')
fprintf('\n')
fprintf('The most recent p-value is:   %0.4f \n', pval_LM(sim))

if pval_LM(sim)>=0.05
    fprintf('The p-value is at least 0.05 so the LM test does not reject the null that beta2 is zero. \n')
    fprintf('\n')
else
    fprintf('The p-value is less than 0.05 so the LM test rejects the null that beta2 is zero.')
    fprintf('\n')
end
fprintf('\n')

fprintf('--------------------------------------------------------------\n')
fprintf('Empirical size of Lagrange Multiplier Test. \n')
fprintf('--------------------------------------------------------------\n')
fprintf('\n')

%{ MISSING CODE %}