%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % 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