% main program to run procedure with or without misreporting
clear all;
format long;
tic;

%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% specifications to change

% Outcome
% 1: Extramarital affairs
% 2: more than 1 partner
outcome = 2


% Smoothing and inverse transformation
% 1: smoothing, no cox, misreporting
% 2: no smoothing, cox, misreporting
% 3: no misreporting, cox
trans = 3;

% version of variables
% 0: age, belief
% 1: age, belief, age_sq
% 2: age, bealief, belief_sq
% 3: age, belief, age_sq, belief_sq
% 4: no squares, imposing no constant
v = 0;  

% using lagged behavior as a regressor
% 0: no
% 1: yes
lag = 0;

% Using also perceptions about prevalence in the 
% 0: no
% 1: yes
prev = 0;

% aggregating beans;
% 2: 0,1, 2 and above
% 3: 0,1,2-4, 5-10
group_bean = 2;



% Grid of alphas
A = [0:0.05:0.6];
if trans == 3
    A = 0;
end
A_num = length(A);

% Grid of tuning parameters
if trans ==1
    T = [0.1 1 10];
else
    T = 0;
end
T_num = length(T);



% whether to calculate marginal effects
% 0: no
% 1: yes
marginal = 1;

%joint significance test
test = 0


%% Other specifications

% number of minimal observations per cell
c = 3;

% program parameters
a = 1; % age intervals
b = 1;  % beans intervals for reported beliefs

% grouping ages
% 1: fixed intervals
% 2: 0.33 and 0.66 percentiles
% 3: .25, 0.5,0.75 percentiles
% 4: 5 groups
% 5: deciles
group_age = 4;

%% reading the row data


% order of variables is: (1:age, 2:belief06, 3:belief08, 4:extra04,
% 5:extra06, 6:extra08, 7:verbal06, 8:out_of_ten06, 9:out_of_ten08)
if prev ==1 
	if outcome == 1
	    X = dlmread('extra_perc.out');
	elseif outcome == 2
	    X = dlmread('number_perc.out');
	end
elseif prev == 0
	if outcome == 1
	    X = dlmread('extra.out');
	elseif outcome == 2
	    X = dlmread('number.out');
	end
	X(:,8:9) = 0;
end

    
% number of observations
N = size(X,1);

if group_age == 1
    X(:,1) = floor(X(:,1)/a);
    age = X(:,1);
elseif group_age == 2
    age_1 = (X(:,1) > 39);
    age_2 = (X(:,1) > 51);
    age = double(age_1+age_2);
    clear age_1 age_2;
elseif group_age == 3
    age_1 = (X(:,1) > 37);
    age_2 = (X(:,1) > 45);
    age_3 = (X(:,1) > 54);
    age = double( age_1 + age_2 + age_3 );
    clear age_1 age_2 age_3;
elseif group_age == 4
    age_1 = (X(:,1) > 36);
    age_2 = (X(:,1) > 42);
    age_3 = (X(:,1) > 49);
    age_4 = (X(:,1) > 56);
    age = double( age_1 + age_2 + age_3 + age_4);
    clear age_1 age_2 age_3 age_4;    
elseif group_age == 5
    age_1 = (X(:,1) > 32);
    age_2 = (X(:,1) > 36);
    age_3 = (X(:,1) > 39);
    age_4 = (X(:,1) > 42);
    age_5 = (X(:,1) > 45);
    age_6 = (X(:,1) > 49);
    age_7 = (X(:,1) > 52);
    age_8 = (X(:,1) > 56);
    age_9 = (X(:,1) > 61);
    age = double( age_1 + age_2 + age_3 + age_4 + age_5 + age_6 + age_7 + age_8 + age_9 );
    clear age_1 age_2 age_3 age_4 age_5 age_6 age_7 age_8 age_9;            
end

% grouping beliefs 
if group_bean == 2
    bean_1 = (X(:,2)>=1);
    bean_2 = (X(:,2)>=2);
    bean06 = bean_1 + bean_2;
    clear bean_1 bean_2;
    bean_1_08 = (X(:,3)>=1);
    bean_2_08 = (X(:,3)>=2);
    bean08 = bean_1_08 + bean_2_08;
    clear bean_1_08 bean_2_08;    
elseif group_bean == 3
    bean_1 = (X(:,2)>=1);
    bean_2 = (X(:,2)>=2);
    bean_3 = (X(:,2)>=5);
    bean06 = bean_1 + bean_2 + bean_3;
    clear bean_1 bean_2 bean_3; 
    bean_1_08 = (X(:,3)>=1);
    bean_2_08 = (X(:,3)>=2);
    bean_3_08 = (X(:,3)>=5);
    bean08 = bean_1_08 + bean_2_08 + bean_3_08;
    clear bean_1_08 bean_2_08 bean_3_08;  
end

% grouping percieved prevalence
prev_1 = (X(:,8) > 1);
prev_2 = (X(:,8) > 2);
prev_3 = (X(:,8) > 3);
prev06 = prev_1 + prev_2 + prev_3;
clear prev_1 prev_2 prev_3;
prev_1_08 = (X(:,9) > 1);
prev_2_08 = (X(:,9) > 2);
prev_3_08 = (X(:,9) > 3);
prev08 = prev_1_08 + prev_2_08 + prev_3_08;
clear prev_1 prev_2 prev_3;




%% Information for "lag" cells (conditional on age, behavior i 2004 and
% beliefs in 2006)

% creating a vector of unique cell IDs
if prev == 0
	id = bean08 + X(:,5)*100 + age*10000 + bean06*1000000 + X(:,4)*100000000;
elseif prev == 1
	id = prev08 + bean08*100 + X(:,5)*10000 + age*1000000 + prev06*100000000 + bean06*10000000000 + X(:,4)*1000000000000;
end


clear age bean06 bean08 prev06 prev08;

% order of variables is: (1:id, 2:age, 3:belief06, 4:belief08, 5:extra04, 6:extra06, 7:extra08, 8:verbal06, 9:out_of_ten06, 10:out_of_ten08)
X = [ id X];
clear id;
X = sortrows(X);

% creating a unique id for lag cells only;
if prev == 0 
    id_lag = floor(X(:,1)./1000);
elseif prev == 1 
    id_lag = floor(X(:,1)./100000);
end
    
    
% number of observations in each cell, sorted by id, length of vector is
% number of cells
obs_lag_cells = diff([ 0; find(id_lag(1:end-1) ~= id_lag(2:end)); length(id_lag) ]);
% number of cells (with any number of observations
J_06 = length(obs_lag_cells);

% number of lag cells used in the estimation (with at least 'c'
% respondents)
J_06_c = sum(obs_lag_cells>=c);

% Creating a block diagonal matriX from square matrices of ones, in the
% size of the different cells;
diago_06 = blkdiag(ones(obs_lag_cells(1)),ones(obs_lag_cells(2)));
for i = 3:J_06
    diago_06 = blkdiag(diago_06, ones(obs_lag_cells(i)));
end
clear i;

% creating the vector of number of people in the cell to which the
% respondent belongs
temp_1 = ones(1,N)*diago_06;

% creating the vector of proportion of people cheating in 2006 in each cell
temp_2 = (X(:,6)'*diago_06)./temp_1;

X = [X temp_1' temp_2'];  % Columns : (1:id, 2:age, 3:belief 06, 4:belief08, 5:extra04, 6:extra06, 7:extra08, 8:verbal06, 9:out_of_ten06, 10:out_of_ten08, 11:# obs in lag cell, 12:propensity in lag cell)
clear temp_1 temp_2 ;

% Vector marking respondents belonging to cells of at least 'c'
% observations
mark = (X(:,11) >= c);

% marking the cells with at least 'c' observations
index = double(obs_lag_cells>=c);

% sample size
% the whole sample: 488
N_all = N;
% number of observations in cells with at least 'c' observations
N = sum(mark);



%% Creating the Z's of rlag cells

% vector of lag cell ids, each ID appearing once
val = id_lag([ (find(id_lag(1:end-1) ~= id_lag(2:end)))' length(id_lag(:)) ]);

% according to the order in 'val', determining to which cell the respondent
% belongs
% A 'J_06xN_all' matrix
Z_06 = (val*ones(1,N_all)-(id_lag(:)*ones(1,J_06))') == 0;
Z_06 = double(Z_06);

clear id_lag val;

%% Information for current cells (conditional on age, behavior in 04 & 06,
% beliefs in 06 & 08)

% number of obsevations in each current cell
obs_cells = diff([ 0; find(X(1:end-1,1) ~= X(2:end,1)); length(X(:,1)) ]);
% number of current cells
J_08 = length(obs_cells);

diago = blkdiag(ones(obs_cells(1)),ones(obs_cells(2)));
for i = 3:J_08
    diago = blkdiag(diago, ones(obs_cells(i)));
end

% creating the vector of number of people in the cell to which the
% respondent belongs
temp_1 = ones(1,N_all)*diago;

% creating the vector of proportion of people cheating in 2006 in each cell
temp_2 = (X(:,7)'*diago)./temp_1;



% vector of current cell ids, each ID appearing once
val = X([ (find(X(1:end-1,1) ~= X(2:end,1)))' length(X(:,1)) ], 1);
% Creating the current Z's
Z_08 = (val*ones(1,N_all)-(X(:,1)*ones(1,J_08))') == 0;
Z_08 = double(Z_08);

Z = [ Z_06' Z_08' ]';


% moving the ID column to the end
X = [X(:,2:end-2) temp_1' temp_2' X(:,end-1:end) X(:,1)]; 
% Columns : (1:age, 2:belief 06, 3:belief08, 4:extra04, 5:extra06, 6:extra08, 7:verbal06, 8:out_of_ten06, 9:out_of_ten08, 10:# obs in current cell, 11:propensity in current cell, 12:# obs in lag cell, 13:propensity in lag cell, 14:id)
clear temp_1 temp_2 diago val;

%% verbal score dummies;
no = (X(:,7) == 0);
low = (X(:,7) == 1);
med = (X(:,7) == 2);
high = (X(:,7) == 3);
V = [ no' ; low' ; med' ; high'];


clear no low med high;

%% OBtaining the estimation to center the moments

% Creating the weighting matriX for the GMM
W = zeros(J_06);
for i=1:N_all
    W = W + Z_06(:,i)*(Z_06(:,i)');
end

W = (1/N)*blkdiag(W, eye(6));


%% given a tuning parameter, calculating the objective function at different
% alphas 


center = zeros(J_06+6,1);

    
for i = 1:T_num %loop over tuning parameters
    for j = 1:A_num % loop over alphas

        t = T(i);
        alpha = A(j);
        objective(((i-1)*A_num + j),:)= objective_mis(X,center,index,V, v, trans, alpha, t, group_bean, J_06, N_all, N, Z_06, W, mark, J_06_c, lag, prev);
                                        
    end
end

csvwrite('out_objective', objective); 


%% CALCULATING Marginal effect, only without misreporting
% Columns : (1:age, 2:belief 06, 3:belief08, 4:extra04, 5:extra06,% 6:extra08, 7:verbal06, 8:out_of_ten06, 9:out_of_ten08, 10:# obs in current cell, 11:propensity in current cell, 12:# obs in lag cell, 13:propensity in lag cell, 14:id)
% mean age in 2008 is 46
% mean age in 2006 is 44
age06 = 44;
age08 = 46;

if marginal ==1 && trans == 3


   theta = objective(6:end-1);

   % Marginal effects
if prev == 0 
   change = [ 0 10; 1 10; 2 10; 3 10; 4 10; 5 10; 6 10; 7 10; 8 10; 9 10; 1 0; 2 0; 3 0; 4 0; 5 0; 6 0; 7 0; 8 0; 9 0;];
   for jjj = 1:length(change(:,1))
       if prev == 0 
           first_06 = theta(3)*(change(jjj,1)*ones(N_all,1) - X(:,2)) + log( (X(:,13)+( (100*X(:,12)).^(-1) ))./ (ones(N_all,1)-X(:,13)+( (X(:,12).*100).^(-1) )));
           second_06 = theta(3)*(change(jjj,2)*ones(N_all,1) - X(:,2)) + log( (X(:,13)+( (100*X(:,12)).^(-1) ))./ (ones(N_all,1)-X(:,13)+( (X(:,12).*100).^(-1) )));
           % 2008
           first_08 = theta(3)*(change(jjj,1)*ones(N_all,1) - X(:,3)) + log( (X(:,11)+( (100*X(:,10)).^(-1) ))./ (ones(N_all,1)-X(:,11)+( (X(:,10).*100).^(-1) )));
           second_08 = theta(3)*(change(jjj,2)*ones(N_all,1) - X(:,3)) + log( (X(:,11)+( (100*X(:,10)).^(-1) ))./ (ones(N_all,1)-X(:,11)+( (X(:,10).*100).^(-1) )));
           if v == 2
               first_06 = first_06 + theta(4)*((change(jjj,1)^2)*ones(N_all,1) - (X(:,2).^2));
               second_06 = second_06 + theta(4)*((change(jjj,2)^2)*ones(N_all,1) - (X(:,2).^2));
               % 2008
               first_08 = first_08 + theta(4)*((change(jjj,1)^2)*ones(N_all,1) - (X(:,3).^2));
               second_08 = second_08 + theta(4)*((change(jjj,2)^2)*ones(N_all,1) - (X(:,3).^2));
           end
           if v == 3
               first_06 = first_06 + theta(5)*((change(jjj,1)^2)*ones(N_all,1) - (X(:,2).^2));
               second_06 = second_06 + theta(5)*((change(jjj,2)^2)*ones(N_all,1) - (X(:,2).^2));
               % 2008
               first_08 = first_08 + theta(5)*((change(jjj,1)^2)*ones(N_all,1) - (X(:,3).^2));
               second_08 = second_08 + theta(5)*((change(jjj,2)^2)*ones(N_all,1) - (X(:,3).^2));
           end
       
        end
       logi_first_06 = ones(N_all,1)./(ones(N_all,1) + exp(-first_06));
       logi_first_08 = ones(N_all,1)./(ones(N_all,1) + exp(-first_08));
       logi_second_06 = ones(N_all,1)./(ones(N_all,1) + exp(-second_06));
       logi_second_08 = ones(N_all,1)./(ones(N_all,1) + exp(-second_08));

       me_06 = (1/sum(X(:,12)>=c))*sum((logi_second_06 - logi_first_06).*(X(:,12)>=c));
       me_08 = (1/sum(X(:,10)>=c))*sum((logi_second_08 - logi_first_08).*(X(:,10)>=c));

       prob_06 = (1/sum(X(:,12)>=c))*sum(ones(N_all,1)./(ones(N_all,1) + exp(-first_06)).*(X(:,12)>=c));
       prob_06_2 = (1/sum(X(:,12)>=c))*sum(ones(N_all,1)./(ones(N_all,1) + exp(-second_06)).*(X(:,12)>=c));
       prob_08 = (1/sum(X(:,10)>=c))*sum(ones(N_all,1)./(ones(N_all,1) + exp(-first_08)).*(X(:,10)>=c));
       prob_08_2 = (1/sum(X(:,10)>=c))*sum(ones(N_all,1)./(ones(N_all,1) + exp(-second_08)).*(X(:,10)>=c));


       me_mat(jjj,:) = [ group_bean age06 age08 change(jjj,1) change(jjj,2) me_06 me_08 prob_06 prob_06_2 prob_08 prob_08_2];
    
       
   end
   csvwrite('me_mat', me_mat);
   clear change first_06 first_08 second_06 second_08 logi_first_06 logi_first_08 logi_second_06 logi_second_08;
end


end


%% Calculating the SE of the estimator
% Columns : (1:age, 2:belief 06, 3:belief08, 4:extra04, 5:extra06, 6:extra08, 7:verbal06, 8:out_of_ten06, 9:out_of_ten08, 10:# obs in current cell, 11:propensity in current cell, 12:# obs in lag cell, 13:propensity in lag cell, 14:id)





if trans ==3
    theta = objective(6:end-1);
    J = J_06 + J_08;

    % Creating vector of propensities by cells
    if prev == 0
        X(:,15) = floor(X(:,14)./1000);
    elseif prev == 1 
       X(:,15) = floor(X(:,14)./100000);
    end
    
    p_06 = X( [ (find(X(1:end-1,15) ~= X(2:end,15)))' length(X(:,15)) ] , 13);
    p_08 = X( [ (find(X(1:end-1,14) ~= X(2:end,14)))' length(X(:,14)) ] , 11);
    p = [p_06' p_08']';

    f_08 = log( (X(:,11)+( (100*X(:,10)).^(-1) ))./ (ones(N_all,1)-X(:,11)+( (X(:,10).*100).^(-1) )));
    f_06 = log( (X(:,13)+( (100*X(:,12)).^(-1) ))./ (ones(N_all,1)-X(:,13)+( (X(:,12).*100).^(-1) ))); 

    % vector of number of observations in each cell;
    nj_06 = obs_lag_cells;
    nj_08 = obs_cells;
    nj = [nj_06' nj_08']';

    % derivatives of epsilon w.r.t p_08 & p_06
    deps_p08 = ( (ones(N_all,1) + 2*( (2*X(:,10)).^(-1)))./((X(:,11) + ((2*X(:,10)).^(-1))).*( ones(N_all,1) - X(:,11) + ( (2*X(:,10)).^(-1))))).*mark;
    deps_p06 = -( (ones(N_all,1) + 2*( (2*X(:,12)).^(-1)))./((X(:,13) + ((2*X(:,12)).^(-1))).*( ones(N_all,1) - X(:,13) + ( (2*X(:,12)).^(-1))))).*mark;   

    if prev == 0 
        if v == 0 || v == 5
            % calculating epsilon_hat for each respondent
            epsilon_hat = f_08 - f_06 - 2*theta(2)*ones(N_all,1) - theta(3)*(X(:,3)-X(:,2)) ;
            % calculating etta for 06
            eta_06 = f_06 - theta(1)*ones(N_all,1) - theta(2)*(X(:,1)-2*ones(N_all,1)) - theta(3)*X(:,2); 
            % calculating eta for 08
            eta_08 = f_08 - theta(1)*ones(N_all,1) - theta(2)*X(:,1) - theta(3)*X(:,3);
        elseif v == 1
            % calculating epsilon_hat for each respondent
            epsilon_hat = f_08 - f_06 - 2*theta(2) - theta(3)*(X(:,3)-X(:,2)) - theta(4)*( (X(:,1).^2)-((X(:,1)-2*ones(N_all,1)).^2));
            % calculating etta for 06
            eta_06 = f_06 - theta(1)*ones(N_all,1) - theta(2)*(X(:,1)-2*ones(N_all,1)) - theta(3)*X(:,2) - theta(4)*((X(:,1)-2*ones(N_all,1)).^2); 
            % calculating eta for 08
            eta_08 = f_08 - theta(1)*ones(N_all,1) -  theta(2)*X(:,1) - theta(3)*X(:,3) - theta(4)*(X(:,1).^2);
        elseif v == 2
            % calculating epsilon_hat for each respondent
            epsilon_hat = f_08 - f_06 - 2*theta(2)*ones(N_all,1) - theta(3)*(X(:,3)-X(:,2)) - theta(4)*((X(:,3).^2)-(X(:,2).^2));
            % calculating etta for 06
            eta_06 = f_06 - theta(1)*ones(N_all,1) - theta(2)*(X(:,1)-2*ones(N_all,1)) - theta(3)*X(:,2) - theta(4)*(X(:,2).^2); 
            % calculating eta for 08
            eta_08 = f_08 - theta(1)*ones(N_all,1) - theta(2)*X(:,1) - theta(3)*X(:,3) - theta(4)*(X(:,3).^2);
        elseif v == 3
            % calculating epsilon_hat for each respondent
            epsilon_hat = f_08 - f_06 - 2*theta(2)*ones(N_all,1) - theta(3)*(X(:,3)-X(:,2)) - theta(4)*( (X(:,1).^2)-((X(:,1)-2*ones(N_all,1)).^2)) - theta(5)*((X(:,3).^2)-(X(:,2).^2));
            % calculating etta for 06
            eta_06 = f_06 - theta(1)*ones(N_all,1) - theta(2)*(X(:,1)-2*ones(N_all,1)) - theta(3)*X(:,2) - theta(4)*((X(:,1)-2*ones(N_all,1)).^2) - theta(5)*(X(:,2).^2); 
            % calculating eta for 08
            eta_08 = f_08 - theta(1)*ones(N_all,1) - theta(2)*X(:,1) - theta(3)*X(:,3) - theta(4)*(X(:,1).^2) - theta(5)*(X(:,3).^2);
        elseif v == 4
            % calculating epsilon_hat for each respondent
            epsilon_hat = f_08 - f_06 - 2*theta(1)*ones(N_all,1) - theta(2)*(X(:,3)-X(:,2)) ;
            % calculating etta for 06
            eta_06 = f_06 - theta(1)*(X(:,1)-2*ones(N_all,1)) - theta(2)*X(:,2); 
            % calculating eta for 08
            eta_08 = f_08 - theta(1)*X(:,1) - theta(2)*X(:,3);   
        end
    elseif prev == 1 
         if v == 0 
            % calculating epsilon_hat for each respondent
            epsilon_hat = f_08 - f_06 - 2*theta(2)*ones(N_all,1) - theta(3)*(X(:,3)-X(:,2)) -theta(4)*(X(:,9)-X(:,8));
            % calculating etta for 06
            eta_06 = f_06 - theta(1)*ones(N_all,1) - theta(2)*(X(:,1)-2*ones(N_all,1)) - theta(3)*X(:,2) - theta(4)*X(:,8); 
            % calculating eta for 08
            eta_08 = f_08 - theta(1)*ones(N_all,1) - theta(2)*X(:,1) - theta(3)*X(:,3) - theta(4)*X(:,9);
        elseif v == 1
            % calculating epsilon_hat for each respondent
            epsilon_hat = f_08 - f_06 - 2*theta(2) - theta(3)*(X(:,3)-X(:,2)) -theta(4)*(X(:,9)-X(:,8)) - theta(5)*( (X(:,1).^2)-((X(:,1)-2*ones(N_all,1)).^2));
            % calculating etta for 06
            eta_06 = f_06 - theta(1)*ones(N_all,1) - theta(2)*(X(:,1)-2*ones(N_all,1)) - theta(3)*X(:,2) - theta(4)*X(:,8) - theta(5)*((X(:,1)-2*ones(N_all,1)).^2); 
            % calculating eta for 08
            eta_08 = f_08 - theta(1)*ones(N_all,1) -  theta(2)*X(:,1) - theta(3)*X(:,3) - theta(4)*X(:,9) - theta(5)*(X(:,1).^2);
        elseif v == 2
            % calculating epsilon_hat for each respondent
            epsilon_hat = f_08 - f_06 - 2*theta(2)*ones(N_all,1) - theta(3)*(X(:,3)-X(:,2)) -theta(4)*(X(:,9)-X(:,8)) - theta(5)*((X(:,3).^2)-(X(:,2).^2)) - theta(6)*((X(:,9).^2)-(X(:,8).^2));
            % calculating etta for 06
            eta_06 = f_06 - theta(1)*ones(N_all,1) - theta(2)*(X(:,1)-2*ones(N_all,1)) - theta(3)*X(:,2) - theta(4)*X(:,8) - theta(5)*(X(:,2).^2) - theta(6)*(X(:,8).^2); 
            % calculating eta for 08
            eta_08 = f_08 - theta(1)*ones(N_all,1) - theta(2)*X(:,1) - theta(3)*X(:,3) - theta(4)*X(:,9) - theta(5)*(X(:,3).^2) - theta(6)*(X(:,9).^2);
        elseif v == 3
            % calculating epsilon_hat for each respondent
            epsilon_hat = f_08 - f_06 - 2*theta(2)*ones(N_all,1) - theta(3)*(X(:,3)-X(:,2)) -theta(4)*(X(:,9)-X(:,8)) - theta(5)*( (X(:,1).^2)-((X(:,1)-2*ones(N_all,1)).^2)) - theta(6)*((X(:,3).^2)-(X(:,2).^2)) - theta(7)*((X(:,9).^2)-(X(:,8).^2));
            % calculating etta for 06
            eta_06 = f_06 - theta(1)*ones(N_all,1) - theta(2)*(X(:,1)-2*ones(N_all,1)) - theta(3)*X(:,2) - theta(4)*X(:,8) - theta(5)*((X(:,1)-2*ones(N_all,1)).^2) - theta(6)*(X(:,2).^2) - theta(7)*(X(:,8).^2); 
            % calculating eta for 08
            eta_08 = f_08 - theta(1)*ones(N_all,1) - theta(2)*X(:,1) - theta(3)*X(:,3) - theta(4)*X(:,9) - theta(5)*(X(:,1).^2) - theta(6)*(X(:,3).^2) - theta(7)*(X(:,9).^2);           
        end
    end
        
        
        
    if lag == 1
        epsilon_hat = epsilon_hat - theta(end)*(X(:,5) - X(:,4));
        eta_06 = eta_06 - theta(end)*X(:,4);
        eta_08 = eta_08 - theta(end)*X(:,5);
    end

    sum_deps_p = zeros(J_06,J);
    sum_deta_06_p = zeros(1,J);
    sum_deta_08_p = zeros(1,J);
    sum_d_v_eps_p = zeros(4,J);
    if v ~= 4
        sum_zeta_zeta_prime = zeros(J_06+6+J);
    elseif v ==4
       sum_zeta_zeta_prime = zeros(J_06+4+J);
    end

    for i=1:N_all
        % sum of the derivatives of Z*epsilon w.r.t p
        sum_deps_p = sum_deps_p + Z_06(:,i)*[ Z_06(:,i)'*deps_p06(i) Z_08(:,i)'*deps_p08(i) ];

        if v ~= 4 && v~=5
            % sum of the derivatives of eta_06 w.r.t p
            sum_deta_06_p = sum_deta_06_p + [ Z_06(:,i)'*(-1)*deps_p06(i) zeros(1,J_08)];

            % sum of the derivatives of eta_08 w.r.t p
            sum_deta_08_p = sum_deta_08_p + [ zeros(1,J_06) Z_08(:,i)'*deps_p08(i) ];
        end

        % sum of the derivatives of V*eps w.r.t p
        sum_d_v_eps_p = sum_d_v_eps_p + V(:,i)*[ Z_06(:,i)'*deps_p06(i) Z_08(:,i)'*deps_p08(i) ]; 

        % sum of zeta*zeta'
        if v ~= 4 && v~= 5    
            zeta(:,i) = [(Z_06(:,i)*epsilon_hat(i)*mark(i))' eta_06(i)*mark(i) eta_08(i)*mark(i) (V(:,i)*epsilon_hat(i)*mark(i))' (diag(Z(:,i))*([ones(1,J_06)*X(i,5)*mark(i) ones(1,J_08)*X(i,6)*mark(i)]' - p))']';    
        elseif v== 4
            zeta(:,i) = [(Z_06(:,i)*epsilon_hat(i)*mark(i))' (V(:,i)*epsilon_hat(i)*mark(i))' (diag(Z(:,i))*([ones(1,J_06)*X(i,5)*mark(i) ones(1,J_08)*X(i,6)*mark(i)]' - p))']'; 
        end
        sum_zeta_zeta_prime = sum_zeta_zeta_prime + zeta(:,i)*zeta(:,i)';     
    end;
    clear zeta;

    % sum of derivatives of psi w.r.t p
    if v ~= 4
        sum_dpsi_p = [ sum_deps_p' sum_deta_06_p' sum_deta_08_p' sum_d_v_eps_p']';
    elseif v == 4
        sum_dpsi_p = [ sum_deps_p' sum_d_v_eps_p']';
    end


    % sum of the derivatives of the m function;
    sum_dm = - diag(nj);  

    Q_N = sum_dpsi_p*inv(sum_dm);

    if v ~= 4
        M_hat = [eye(J_06 + 6) -Q_N]*(1/N)*sum_zeta_zeta_prime*[eye(J_06 + 6) -Q_N]';
    elseif v == 4
        M_hat = [eye(J_06 + 4) -Q_N]*(1/N)*sum_zeta_zeta_prime*[eye(J_06 + 4) -Q_N]';
    end

    % sum of derivatives of psi w.r.t theta
    % W.R.T CONSTANT
    sum_dpsi_theta = [ zeros(1,J_06) -N -N zeros(1,4)]';
    % w.r.t age coefficient
    sum_dpsi_theta(:,2) = [ (-Z_06*(2*ones(N_all,1).*mark))' sum((-(X(:,1) - 2*ones(N_all,1))).*mark) sum(-X(:,1).*mark) (-V*(2*ones(N_all,1).*mark))']';
    % w.r.t belief coefficient
    sum_dpsi_theta(:,3) = [ (-Z_06*((X(:,3)-X(:,2)).*mark))' sum(-X(:,2).*mark) sum(-X(:,3).*mark) (-V*((X(:,3)-X(:,2)).*mark))']';

    if prev == 0 
        if v ==1 || v==3
            % w.r.t to age_sq
            sum_dpsi_theta(:,4) = [ (-Z_06*(((X(:,1).^2)-((X(:,1)-2*ones(N_all,1)).^2)).*mark))' sum(-((X(:,1)-2*ones(N_all,1)).^2).*mark) sum(-(X(:,1).^2).*mark) (-V*(((X(:,1).^2)-((X(:,1)-2*ones(N_all,1)).^2)).*mark))' ]';
        elseif v == 2 
            sum_dpsi_theta(:,4) = [ (-Z_06*((((X(:,3).^2)-(X(:,2).^2))).*mark))' sum(-(X(:,2).^2).*mark) sum(-(X(:,3).^2).*mark) (-V*((((X(:,3).^2)-(X(:,2).^2))).*mark))' ]';
        end 

        if v == 3
            sum_dpsi_theta(:,5) = [ (-Z_06*((((X(:,3).^2)-(X(:,2).^2))).*mark))' sum(-(X(:,2).^2).*mark) sum(-(X(:,3).^2).*mark) (-V*((((X(:,3).^2)-(X(:,2).^2))).*mark))' ]';
        end

        if v ==4
            sum_dpsi_theta = [ sum_dpsi_theta(1:J_06,2); sum_dpsi_theta(J_06+3:J_06+6,2) ];

            sum_dpsi_theta(:,2) = [ (-Z_06*((X(:,3)-X(:,2)).*mark.*(X(:,3) > X(:,2))))'  (-V*((X(:,3)-X(:,2)).*mark.*(X(:,3) > X(:,2))))']';
            sum_dpsi_theta(:,3) = [ (-Z_06*((X(:,3)-X(:,2)).*mark.*(X(:,3) < X(:,2))))'  (-V*((X(:,3)-X(:,2)).*mark.*(X(:,3) < X(:,2))))']';
        end
    elseif prev ==1 
        % w.r.t prev
        sum_dpsi_theta(:,4) = [ (-Z_06*((X(:,9)-X(:,8)).*mark))' sum(-X(:,8).*mark) sum(-X(:,9).*mark) (-V*((X(:,9)-X(:,8)).*mark))']';

        
        if v ==1 || v==3 
            % w.r.t to age_sq
            sum_dpsi_theta(:,5) = [ (-Z_06*(((X(:,1).^2)-((X(:,1)-2*ones(N_all,1)).^2)).*mark))' sum(-((X(:,1)-2*ones(N_all,1)).^2).*mark) sum(-(X(:,1).^2).*mark) (-V*(((X(:,1).^2)-((X(:,1)-2*ones(N_all,1)).^2)).*mark))' ]';
        elseif v == 2 
            % w.r.t belief sq
            sum_dpsi_theta(:,5) = [ (-Z_06*((((X(:,3).^2)-(X(:,2).^2))).*mark))' sum(-(X(:,2).^2).*mark) sum(-(X(:,3).^2).*mark) (-V*((((X(:,3).^2)-(X(:,2).^2))).*mark))' ]';
            % w.r.t prev sq
            sum_dpsi_theta(:,6) = [ (-Z_06*((((X(:,9).^2)-(X(:,8).^2))).*mark))' sum(-(X(:,8).^2).*mark) sum(-(X(:,9).^2).*mark) (-V*((((X(:,9).^2)-(X(:,8).^2))).*mark))' ]';
        end 

        if v == 3 
            % w.r.t belief sq
            sum_dpsi_theta(:,6) = [ (-Z_06*((((X(:,3).^2)-(X(:,2).^2))).*mark))' sum(-(X(:,2).^2).*mark) sum(-(X(:,3).^2).*mark) (-V*((((X(:,3).^2)-(X(:,2).^2))).*mark))' ]';
        end
        
        if v == 3 
            % w.r.t prev sq
            sum_dpsi_theta(:,7) = [ (-Z_06*((((X(:,9).^2)-(X(:,8).^2))).*mark))' sum(-(X(:,8).^2).*mark) sum(-(X(:,9).^2).*mark) (-V*((((X(:,9).^2)-(X(:,8).^2))).*mark))' ]';
        end
        
    end
        
 
    if lag == 1
        sum_dpsi_theta(:,end+1) = [ (-Z_06*((X(:,5)-X(:,4)).*mark))' sum(-X(:,4).*mark) sum(-X(:,5).*mark) (-V*((X(:,5)-X(:,4)).*mark))']';

    end



    D_hat = sum_dpsi_theta/N;

    A_N = W; 

    V_hat = (1/N)*inv(D_hat'*A_N*D_hat)*D_hat'*A_N*M_hat*A_N*D_hat*inv(D_hat'*A_N*D_hat);
    
    out_with_se = [objective(3:end-1) diag(V_hat)'.^0.5];
    
    
    % joint significance test
    if test == 1
        m_test =  [0 0 1 0 0; 0 0 0 0 1; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0]*out_with_se(4:8)'
        wald = m_test'*inv(V_hat)*m_test
    end
    
end
 


































