%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Xavier Martin G. Bautista
% Econometrics 2
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This code estimates the probit model using maximum likelihood when there 
% is sample selectivity.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
clc

global x w y z N

% Set the true parameters and placeholders for results.

select = 0;                                                                 % 0 = no control for selectivity
                                                                            % 1 = control for selectivity

N = 300;                                                                    % Number of observations.
beta0 = 0.5;                                                                % Intercept for y0.
beta1 = 3;                                                                  % Coefficient on X.
sigmau = 2;                                                                 % Standard deviation for u. Standard deviation for v is set to 1.
gamma0 = 1;                                                                 % Intercept for z0.
gamma1 = 4;                                                                 % Coefficient on W.
rho = 0.8;                                                                  % Correlation between u and v.

sim = 30;                                                                   % Number of simulations.

if select == 0                                                              % Results matrix.
    results_mat = zeros(sim, 3);                                            
elseif select == 1
    results_mat = zeros(sim, 6);                                            
end

% Begin simulation.

for s = 1:sim 

    % Generate the data. The latent variable model is:
    %           y0(j) = beta0 + beta1*x(j) + u(j)
    %           z0(j) = gamma0 + gamma1*w(j) + v(j)
    % where y(j) = is observed if z0(j) > 0 and unobserved otherwise and
    % z(j) = 1 if z0(j) > 0 and 0 otherwise.
  
    e = mvnrnd([0; 0],[sigmau^2 rho*sigmau; rho*sigmau 1],N);               % Correlated error terms. Standard deviation of v is set to 1 as in Davidson-Mackinnon.
    
    u = e(:,1);                                                             % Error terms for y.                                                     
    v = e(:,2);                                                             % Error terms for q.
                   
    x = ((1:N)'/N).*normrnd(0,1,N,1);
    w = ((1:N)'/N).*normrnd(0,1,N,1);

    y0 = beta0*ones(size(x,1),1) + beta1*x + u;                             % Latent variable y.
    z0 = gamma0*ones(size(w,1),1) + gamma1*w + v;                           % Latent variable z.
    
    y = y0.*double((z0 > 0));                                               % Observed y.
    z = double((z0 > 0));                                                   % Observed z.
    
    options = optimoptions(@fminunc,'Display','off',...                     % Optimization options.
                                               'Algorithm','quasi-newton');
    
    if select == 0

        b0 = [1 1 1];                                                       % Initial values.

        [b_mle, ~, ~, ~, ~, hess] = fminunc('logl_ssbias', b0, options);    % Minimization.
    
    elseif select == 1

        b0 = [1 1 1 1 1 0.05];                                              % Initial values.
    
        [b_mle, ~, ~, ~, ~, hess] = fminunc('logl_ss', b0, options);        % Minimization.
        
    end

    results_mat(s, 1:size(results_mat,2)) = b_mle';                         % Store results.
    vmat = inv(hess);
    
% Display the results of each simulation.
 
    fprintf('--------------------------------------------------------------\n')
    fprintf('\n')
    
    fprintf('Simulation %d \n',s)
    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('sigmau    %0.4f       %0.4f       %0.4f \n', b_mle(3), sqrt(vmat(3,3)), b_mle(3)/sqrt(vmat(3,3)))
   	if select == 1        
    fprintf('g0        %0.4f       %0.4f       %0.4f \n', b_mle(4), sqrt(vmat(4,4)), b_mle(4)/sqrt(vmat(4,4)))
    fprintf('g1        %0.4f       %0.4f       %0.4f \n', b_mle(5), sqrt(vmat(5,5)), b_mle(5)/sqrt(vmat(5,5)))
    fprintf('rho       %0.4f       %0.4f       %0.4f \n', b_mle(6), sqrt(vmat(6,6)), b_mle(6)/sqrt(vmat(6,6)))
    end
    
    fprintf('\n')

end

% Display the results.

fprintf('--------------------------------------------------------------\n')
fprintf('Results of Last Simulation \n')
fprintf('--------------------------------------------------------------\n')
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('sigmau    %0.4f       %0.4f       %0.4f \n', b_mle(3), sqrt(vmat(3,3)), b_mle(3)/sqrt(vmat(3,3)))
if select == 1
fprintf('g0        %0.4f       %0.4f       %0.4f \n', b_mle(4), sqrt(vmat(4,4)), b_mle(4)/sqrt(vmat(4,4)))
fprintf('g1        %0.4f       %0.4f       %0.4f \n', b_mle(5), sqrt(vmat(5,5)), b_mle(5)/sqrt(vmat(5,5)))
fprintf('rho       %0.4f       %0.4f       %0.4f \n', b_mle(6), sqrt(vmat(6,6)), b_mle(6)/sqrt(vmat(6,6)))
end
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)))
fprintf('sigmau  %0.4f       %0.4f     \n', mean(results_mat(:,3)), std(results_mat(:,3)))
if select == 1
fprintf('g0      %0.4f       %0.4f     \n', mean(results_mat(:,4)), std(results_mat(:,4)))
fprintf('g1      %0.4f       %0.4f     \n', mean(results_mat(:,5)), std(results_mat(:,5)))
fprintf('rho     %0.4f       %0.4f     \n', mean(results_mat(:,6)), std(results_mat(:,6)))
end