%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Xavier Martin G. Bautista
% Econometrics 2
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This code estimates the parameters of a duration model which has and 
% exponential distribution using maximum likelihood.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
clc

global x1 x2 T t

% Set the true parameters and placeholders for results.

T = 100;                                                                    % Number of observations.
beta1 = 3;                                                                  % Coefficient on X1.
beta2 = 1;                                                                  % C0efficient on X2.

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

% Begin simulation.

for s = 1:sim 

    % Generate the data. 
    %              t(i) = exp(x1(i)*beta1 + x2(i)*beta2)*u(i)
    % where the distribution of u(i) is the exponential distribution.
    
    x1 = normrnd(0,1,T,1);                                                
    x2 = normrnd(0,1,T,1);
   
    u = exp(exprnd(ones(T,1)));                                             % Exponential draws.
    t = exp(-(x1.*beta1 + x2.*beta2)).*u;                                   % Duration data. 

    b0 = [0 0];                                                             % Initial values.

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

    results_mat(s,1:size(results_mat,2)) = b_mle';                          % Store results.
  
% Display the results of each simulation.
 
    fprintf('--------------------------------------------------------------\n')
    fprintf('\n')
    
    fprintf('Simulation %d \n',s)
    fprintf('        b          SE             t \n')
    fprintf('b1   %0.4f       %0.4f       %0.4f \n', b_mle(1), sqrt(inv(hess(1,1))), b_mle(1)/sqrt(inv(hess(1,1))))
    fprintf('b2   %0.4f       %0.4f       %0.4f \n', b_mle(2), sqrt(inv(hess(2,2))), b_mle(2)/sqrt(inv(hess(2,2))))
    fprintf('\n')

end

vmat = inv(hess);

% Display the results.

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

fprintf('Hessian Matrix \n')
fprintf('%0.4f       %0.4f     \n', hess(1,1), hess(1,2))
fprintf('%0.4f       %0.4f     \n', hess(2,1), hess(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('b1   %0.4f       %0.4f       %0.4f \n', b_mle(1), sqrt(vmat(1,1)), b_mle(1)/sqrt(vmat(1,1)))
fprintf('b2   %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('b1      %0.4f       %0.4f     \n', mean(results_mat(:,1)), std(results_mat(:,1)))
fprintf('b2      %0.4f       %0.4f     \n', mean(results_mat(:,2)), std(results_mat(:,2)))