%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Xavier Martin G. Bautista
% Econometrics 2
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This code estimates the probit model using maximum likelihood.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
clc

global x z N

% Set the true parameters and placeholders for results.

N = 300;                                                                    % Number of observations.
beta0 = 0.5;                                                                % Intercept.
beta1 = 3;                                                                  % Coefficient on X.
sigma = 1;                                                                  % Standard deviation.

sim = 10;                                                                   % Number of simulations.
results_mat = zeros(sim, 2);                                                % Results matrix.

% Begin simulation.

for s = 1:sim 

    % Generate the data. The latent variable model is:
    %           y(j) = beta0 + beta1*x(j) + u(j)
    % where z(j) = 1 if y(j) is positive and z(j) is zero otherwise.
  
    u = normrnd(0,sigma,N,1);
    x = ((1:N)'/N).*normrnd(0,1,N,1);

    y = beta0*ones(size(x,1),1) + beta1*x + u;
    z = double((y > 0));

    b0 = [1 1];                                                             % Initial values.

    options = optimset('Display','off');                                    % Turn off display.
    [b_mle, ~, ~, ~, ~, hess] = fminunc('logl_prob', b0, options);          % Minimization.

    results_mat(s, 1:size(results_mat,2)) = b_mle';                         % Store results.

end


H = hess;
vmat = inv(H);

% Display the results.

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

fprintf('Hessian Matrix \n')
fprintf('%0.4f       %0.4f     \n', H(1,1), H(1,2))
fprintf('%0.4f       %0.4f     \n', H(2,1), H(2,2))
fprintf('\n')

fprintf('Variance Matrix \n')
fprintf('%0.4f       %0.4f     \n', vmat(1,1), vmat(1,2))
fprintf('%0.4f       %0.4f     \n', vmat(2,1), vmat(2,2))
fprintf('\n')

fprintf('        b            SE          t \n')
fprintf('b0   %0.4f       %0.4f       %0.4f \n', b_mle(1), sqrt(vmat(1,1)), b_mle(1)/sqrt(vmat(1,1)))
fprintf('b1   %0.4f       %0.4f       %0.4f \n', b_mle(2), sqrt(vmat(2,2)), b_mle(2)/sqrt(vmat(2,2)))
fprintf('\n')

fprintf('--------------------------------------------------------------\n')
fprintf('Empirical Results \n')
fprintf('--------------------------------------------------------------\n')
fprintf('\n')

fprintf('         Mean       Std.Dev   \n')
fprintf('b0      %0.4f       %0.4f     \n', mean(results_mat(:,1)), std(results_mat(:,1)))
fprintf('b1      %0.4f       %0.4f     \n', mean(results_mat(:,2)), std(results_mat(:,2)))