function output = rddmmsef(X,Y,D,c)

%--------------------------------------------------------------------------
% rddmmse.m: compute the sharp regression discontinuity estimator
%--------------------------------------------------------------------------
%
% DESCRIPTION:
%   Estimate the fuzzy regression discontinuity estimator with the
%   bandwidths proposed by Arai and Ichimura (2015).
%
% REFERENCE:
%   Arai and Ichimura (2015). "Optimal Bandwidth Selection for the 
%   Fuzzy Regression Discontinuity Estimator," GRIPS Discussion Paper,
%   15-13.
%
% 
% INPUT ARGUMENTS
%   X: Imput vector (n by 1) of the forcing variable 
%   Y: Vector of the outcome variable (n by 1) 
%   D: Treat mentstatus (n by 1)
%   c: Evaluation (discontinuity) point (scalar)
%
% OUTPUT ARGUMENTS:
%   tau:    fuzzy RD estimate
%   se:     standard error
%   h:      bandwidths, the first element for the right of the cut-off
%                   the second element for the left of the cut-off
%   tval:   bias-corrected t-value
%      ci:  bias-corrected confidence interval, See Section 4.4 of
%               Fan and Gijbels (1996)
%       n:  the number of observations with nonzero weight,
%               the first element for the right of the cut-off,
%               the second element for the left of the cut-off
% 
% AUTHOR:
%   This code was written by Yoichi Arai.
%   I will appreciate any feedbacks.
%
% Note:
%   Optimization toolbox is necessary.
% 
%
% 
% DEVELOPMENT: 
%    4 Septermber 2015: Original version of rddmmse.m written.
%--------------------------------------------------------------------------

% Check input arguments
if nargin < 4
    error('more input arguments needed.');
end

if nargin > 4
    error('Too many arguments.');
end

% Step 1. Estimation of density and its first derivative
    % Density Estimation
     n = size(X,1);
     lambda = sqrt(var(X));
     a_pilot = 2.34*lambda/n^(1/5);
     f_k = kdensity_est(X,c,a_pilot); % kernel = 2, Epanechnikov
     h1_pilot = lambda*(105*16*sqrt(pi)/15/n)^(1/7);
     f1_pilot = density1d_est(X,c,h1_pilot);         

% Step 2: Get pilot bandwidths to estimate derivatives for m_y
 % 4th order Polynomial regression
 Xbin_lt_x = (X<c);           % binary indicator for observations less than x
 Xbin_bt_x = (X>=c);          % bigger than x
 num_X_bt_x = sum(Xbin_bt_x); % # of observations of X bigger than x, N+
 num_X_lt_x = sum(Xbin_lt_x); % # of observations of X less than x, N- in p.16
 Xx = X-c;
 Xx2 = Xx.^2;
 Xx3 = Xx.*Xx2;
 Xx4 = Xx2.^2;
 nones = ones(n,1);
 data_mat = [Y nones Xx Xx2 Xx3 Xx4];
 data_mat_p = data_mat(Xbin_bt_x==1,:); % data_mat only for positive X
 data_mat_m = data_mat(Xbin_lt_x==1,:); % data_mat only for negative X

 YT_pc = data_mat_p(:,1);
 TT_pc = data_mat_p(:,2:6);
 lambda_pc_4_4g = TT_pc\YT_pc;
 m41_4 = 24*lambda_pc_4_4g(5,1);
 
 YT_mc = data_mat_m(:,1);
 TT_mc = data_mat_m(:,2:6);
 lambda_mc_4_4g = TT_mc\YT_mc;
 m42_4 = 24*lambda_mc_4_4g(5,1);

 vy_p_plugin_4g = sum((YT_pc-TT_pc*lambda_pc_4_4g).^2)/(num_X_bt_x-5);
 vy_m_plugin_4g = sum((YT_mc-TT_mc*lambda_mc_4_4g).^2)/(num_X_lt_x-5);

 % Pilot bandwidths for 3rd order LPR to estimate 3rd derivative with uniform kernel
 band_comp_p = (vy_p_plugin_4g/f_k/m41_4^2/num_X_bt_x)^(1/9);
 band_comp_m = (vy_m_plugin_4g/f_k/m42_4^2/num_X_lt_x)^(1/9);
 C33K_U = Const_LPR_U_est(3,3);
 h3_p = C33K_U*band_comp_p;
 h3_m = C33K_U*band_comp_m;
 % Pilot bandwidths for 3rd order LPR to estimate 2nd derivative with uniform kernel
 C23K_U = Const_LPR_U_est(2,3);
 h23_p = C23K_U*band_comp_p;
 h23_m = C23K_U*band_comp_m;
   
 % Step 3: Estimate the second and the third derivatives of m_y
 % 3rd order LPR for the third derivatives    
 Xbin_bet_xxh = Xbin_bt_x.*(X<=c+h3_p);     % index for X in [x, x+h3_p]
 data_mat_4_p = data_mat(Xbin_bet_xxh==1,:);
 YT_p = data_mat_4_p(:,1);
 T_p = data_mat_4_p(:,2:5);
 lambda_p_all_33_4 = T_p\YT_p;

 Xbin_bet_xhx = (X>=c-h3_m).*Xbin_lt_x;     % index for X in [x-h3_m, x]
 data_mat_4_m = data_mat(Xbin_bet_xhx==1,:);
 YT_m = data_mat_4_m(:,1);
 T_m = data_mat_4_m(:,2:5);
 lambda_m_all_33_4 = T_m\YT_m;   

 m3c_p = 6*lambda_p_all_33_4(4,1);
 m3c_m = 6*lambda_m_all_33_4(4,1);

 % 3rd order LPR for the second derivatives of m_y
 Xbin_bet_xxh = Xbin_bt_x.*(X<=c+h23_p);     % index for X in [x, x+h3_p]
 data_mat_4_p = data_mat(Xbin_bet_xxh==1,:);
 YT_p = data_mat_4_p(:,1);
 T_p = data_mat_4_p(:,2:5);
 lambda_p_all_23 = T_p\YT_p;
 m2c_p_2 = 2*lambda_p_all_23(3,1);

 YT_p_hat = T_p*lambda_p_all_23;
 vy_p_4_2 = sum((YT_p - YT_p_hat).^2)/(size(YT_p,1)-4);

 Xbin_bet_xhx = (X>=c-h23_m).*Xbin_lt_x;     % index for X in [x-h3_m, x]
 data_mat_4_m = data_mat(Xbin_bet_xhx==1,:);
 YT_m = data_mat_4_m(:,1);
 T_m = data_mat_4_m(:,2:5);
 lambda_m_all_23 = T_m\YT_m;
 m2c_m_2 = 2*lambda_m_all_23(3,1);

 YT_m_hat = T_m*lambda_m_all_23;
 vy_m_4_2 = sum((YT_m - YT_m_hat).^2)/(size(YT_m,1)-4);

 m2_pilot_y = [m2c_p_2 m2c_m_2];
 m3_pilot_y = [m3c_p m3c_m];
 v_pilot_y = [vy_p_4_2 vy_m_4_2];


% Step 4: Get pilot bandwidths to estimate derivatives of m_D
% 4th order Polynomial regression 
data_mat2 = [D nones Xx Xx2 Xx3 Xx4];
data_mat2_p = data_mat2(Xbin_bt_x==1,:); % data_mat only for positive X
data_mat2_m = data_mat2(Xbin_lt_x==1,:); % data_mat only for negative X

DT_pc = data_mat2_p(:,1);
lambda_pc_4_4g_d = TT_pc\DT_pc;
m41_4_d = 24*lambda_pc_4_4g_d(5,1);

DT_mc = data_mat2_m(:,1);
lambda_mc_4_4g_d = TT_mc\DT_mc;
m42_4_d = 24*lambda_mc_4_4g_d(5,1);

vy_p_plugin_4g_d = sum((DT_pc-TT_pc*lambda_pc_4_4g_d).^2)/(num_X_bt_x-5);
vy_m_plugin_4g_d = sum((DT_mc-TT_mc*lambda_mc_4_4g_d).^2)/(num_X_lt_x-5);

% Pilot bandiwidths for 3rd order LPR to estimate 3rd derivative with uniform kernel
band_comp_p_d = (vy_p_plugin_4g_d/f_k/m41_4_d^2/num_X_bt_x)^(1/9);
band_comp_m_d = (vy_m_plugin_4g_d/f_k/m42_4_d^2/num_X_lt_x)^(1/9);
h3_p_d = C33K_U*band_comp_p_d;
h3_m_d = C33K_U*band_comp_m_d;

% Pilot bandwidths3rd order LPR to estimate 2nd derivative with uniform kernel
h23_p_d = C23K_U*band_comp_p_d; % possible to use either v_pilot_24 or v_pilot_23
h23_m_d = C23K_U*band_comp_m_d;


% Step 5: Estimate the second and the third derivatives of m_D
% 3rd order LPR for the third derivatives
Xbin_bet_xxh_d = Xbin_bt_x.*(X<=c+h3_p_d);     % index for X in [x, x+h3_p]
data_mat_4_p_d = data_mat2(Xbin_bet_xxh_d==1,:);
DT_p = data_mat_4_p_d(:,1);
T_p_d = data_mat_4_p_d(:,2:5);
lambda_p_all_33_4_d = T_p_d\DT_p;
m3c_p_d = 6*lambda_p_all_33_4_d(4,1);

Xbin_bet_xhx_d = (X>=c-h3_m_d).*Xbin_lt_x;     % index for X in [x-h3_m, x]
data_mat_4_m_d = data_mat2(Xbin_bet_xhx_d==1,:);
DT_m = data_mat_4_m_d(:,1);
T_m_d = data_mat_4_m_d(:,2:5);
lambda_m_all_33_4_d = T_m_d\DT_m;
m3c_m_d = 6*lambda_m_all_33_4_d(4,1);

Xbin_bet_xxh_d = (X>=c).*(X<=c+h23_p_d);     % index for X in [x, x+h3_p]
data_mat_4_p_d = data_mat2(Xbin_bet_xxh_d==1,:);
DT_p = data_mat_4_p_d(:,1);
T_p_d = data_mat_4_p_d(:,2:5);
lambda_p_all_23_d = T_p_d\DT_p;
m2c_p_2_d = 2*lambda_p_all_23_d(3,1);

DT_p_hat = T_p_d*lambda_p_all_23_d;
vy_p_4_2_d = sum((DT_p - DT_p_hat).^2)/(size(DT_p,1)-4);

Xbin_bet_xhx_d = (X>=c-h23_m_d).*(X<=c);     % index for X in [x-h3_m, x]
data_mat_4_m_d = data_mat2(Xbin_bet_xhx_d==1,:);
DT_m = data_mat_4_m_d(:,1);
T_m_d = data_mat_4_m_d(:,2:5);
lambda_m_all_23_d = T_m_d\DT_m;
m2c_m_2_d = 2*lambda_m_all_23_d(3,1);

DT_m_hat = T_m_d*lambda_m_all_23_d;
vy_m_4_2_d = sum((DT_m - DT_m_hat).^2)/(size(DT_m,1)-4);

v_pilot_t = [vy_p_4_2_d vy_m_4_2_d];

% Compute covariances
Xbin_p_ai = Xbin_bet_xxh.*Xbin_bet_xxh_d;
TT = data_mat(:,2:5);
res_y_p = Y - TT*lambda_p_all_23;
res_d_p = D - TT*lambda_p_all_23_d;
cov_p_ai = (res_y_p.*Xbin_p_ai)'*res_d_p/(num_X_bt_x - 4);

Xbin_m_ai = Xbin_bet_xhx.*Xbin_bet_xhx_d;
res_y_m = Y - TT*lambda_m_all_23;
res_d_m = D - TT*lambda_m_all_23_d;
cov_m_ai = (res_y_m.*Xbin_m_ai)'*res_d_m/(num_X_lt_x - 4);

cov_pilot = [cov_p_ai cov_m_ai];

m2_pilot_t = [m2c_p_2_d m2c_m_2_d];
m3_pilot_t = [m3c_p_d m3c_m_d];


% Step 6: Get the pilot estimate for tau
 options = optimset('Algorithm','interior-point','Display','off');
 Qn_value = @(theta)Qn(theta,f_k,f1_pilot,m2_pilot_y,m3_pilot_y,v_pilot_y,n);
 
 % Lower bound for AI search area
 NN_min_AI = 3;
 X_lt_x = data_mat_m(:,3);
 X_bt_x = data_mat_p(:,3);
 sorted_X_bt_x = sort(X_bt_x);               % ascending order
 sorted_X_lt_x_des_abs = abs(sort(X_lt_x,'descend'));
 NN_p_AI = sorted_X_bt_x(NN_min_AI,1);
 NN_m_AI = sorted_X_lt_x_des_abs(NN_min_AI,1);
 h_lb = [NN_p_AI,NN_m_AI];
 % Upper bound for AI search area
 hu1 = sorted_X_bt_x(num_X_bt_x,1);
 hu2 = sorted_X_lt_x_des_abs(num_X_lt_x,1);
 h_ub = [hu1,hu2];

 h_mat_local = zeros(10,2);
 min_val = zeros(10,1);
 h1_range = sorted_X_bt_x(num_X_bt_x,1) - sorted_X_bt_x(3,1);
 h2_range = sorted_X_lt_x_des_abs(num_X_lt_x,1) - sorted_X_lt_x_des_abs(3,1);
 
 for k=1:10
 for j=1:10
     h_ini_1 = sorted_X_bt_x(3,1) + 0.1*h1_range*j;
     h_ini_2 = sorted_X_lt_x_des_abs(3,1) + 0.1*h2_range*k;
     h_ini = [h_ini_1 h_ini_2];
     [theta_2, fval_2] = fmincon(Qn_value,h_ini,[],[],[],[],h_lb,h_ub,[],options);
     h_mat_local(10*(k-1)+j,:) = theta_2;
     min_val(10*(k-1)+j,1) = fval_2;   
 end
 end
 [~, min_index] = min(min_val);
 h_sharp = h_mat_local(min_index,:);
 [tau_pilot,~,~,~,~] = frdd_est(Y,D,X,c,h_sharp(1),h_sharp(2));
 
 
% Step 7: Get bandwidths that minimize MMSE
 m2_pilot = [m2_pilot_y m2_pilot_t];
 m3_pilot = [m3_pilot_y m3_pilot_t];
 v_pilot = [v_pilot_y v_pilot_t cov_pilot];
 
 Qn_value_f = @(theta)Qn_FRDD(theta,f_k,f1_pilot,m2_pilot,m3_pilot,v_pilot,tau_pilot,n);
 
 h_mat_local_f = zeros(10,2);
 min_val_f = zeros(10,1);
 for k=1:10
 for j=1:10
     h_ini_1 = sorted_X_bt_x(3,1) + 0.1*h1_range*j;
     h_ini_2 = sorted_X_lt_x_des_abs(3,1) + 0.1*h2_range*k;
     h_ini = [h_ini_1 h_ini_2];
     [theta_f, fval_f] = fmincon(Qn_value_f,h_ini,[],[],[],[],h_lb,h_ub,[],options);
     h_mat_local_f(j,:) = theta_f;
     min_val_f(j,1) = fval_f;
 end
 end
   [~, min_index] = min(min_val_f);
   h_opt = h_mat_local_f(min_index,:);

   [tau,se,n_p,n_m,tau_d] = frdd_est(Y,D,X,c,h_opt(1),h_opt(2));    
    
% Compute bias corrected t-stat
 % Compute bias correction terms b11 and b12
 Xbin_bt_xh_opt = (X >= c-h_opt(2));
 Xbin_lt_xh_opt = (X <= c+h_opt(1));
 Xbin_bet_xxh_opt = Xbin_bt_x.*Xbin_lt_xh_opt;    % index for X in [x,x+h]
 Xbin_bet_xhx_opt = Xbin_bt_xh_opt.*Xbin_lt_x;    % index for X in [x-h,x]

 data_trim_mat_p_opt = data_mat(Xbin_bet_xxh_opt==1,:);
 T_p_opt = data_trim_mat_p_opt(:,2:3);
 Kvec_p_opt = (1-abs(data_trim_mat_p_opt(:,3))/h_opt(1)).*(abs(data_trim_mat_p_opt(:,3))<h_opt(1));
 W_p_opt = diag(Kvec_p_opt);

 data_trim_mat_m_opt = data_mat(Xbin_bet_xhx_opt==1,:);
 T_m_opt = data_trim_mat_m_opt(:,2:3);
 Kvec_m_opt = (1-abs(data_trim_mat_m_opt(:,3))/h_opt(2)).*(abs(data_trim_mat_m_opt(:,3))<h_opt(2));
 W_m_opt = diag(Kvec_m_opt);
 
 Sn01_p_tilde_inv_opt = invpd(T_p_opt'*W_p_opt*T_p_opt);
 cn2_p_tilde_opt = data_trim_mat_p_opt(:,4:5)'*Kvec_p_opt;
 cn3_p_tilde_opt = data_trim_mat_p_opt(:,5:6)'*Kvec_p_opt;
 b11_vec_opt = (m2c_p_2 - tau*m2c_p_2_d)*Sn01_p_tilde_inv_opt*cn2_p_tilde_opt/2 ...
    + (m3c_p - tau*m3c_p_d)*Sn01_p_tilde_inv_opt*cn3_p_tilde_opt/6;

 Sn01_m_tilde_inv_opt = invpd(T_m_opt'*W_m_opt*T_m_opt);
 cn2_m_tilde_opt = data_trim_mat_m_opt(:,4:5)'*Kvec_m_opt;
 cn3_m_tilde_opt = data_trim_mat_m_opt(:,5:6)'*Kvec_m_opt;
 b12_vec_opt = (m2c_m_2 - tau*m2c_m_2_d)*Sn01_m_tilde_inv_opt*cn2_m_tilde_opt/2 ...
    + (m3c_m - tau*m3c_m_d)*Sn01_m_tilde_inv_opt*cn3_m_tilde_opt/6;

 bias_opt = (b11_vec_opt(1,1)-b12_vec_opt(1,1))/tau_d;

% Construct bias-corrected t-value
 tstat = (tau-bias_opt)/se;

output.tau = tau;
output.se =se;
output.h = h_opt;
output.tstat = tstat;
output.ci = [tau-bias_opt-1.96*se tau-bias_opt+1.96*se];
output.n = [n_p n_m];
    
%------------------------------------------------------------------
% End function rddmmsef
%------------------------------------------------------------------

function fhat = kdensity_est(X,c,h)

% Purpose: estimate density at a point with Epanechnikov kernel
% --------------------------------------------
% Usage: results = density_est
%       X = a column vector of scalar random variables
%       c = point to estimate
%       h = bandwidth
% --------------------------------------------

nobs = size(X,1);
%fhat = 0;
A_h = (X-c)/h;
ExpA_h = 3*sum((abs(A_h)<1).*(1-A_h.^2))/4;
fhat = ExpA_h/nobs/h;

%------------------------------------------------------------------
% End private function density_est
%------------------------------------------------------------------

function f1hat = density1d_est(X,c,h)

% Purpose: estimate 1st derivative density at a point
% --------------------------------------------
% Usage: results = density_est
%       X = a column vector of scalar random variables
%       c = point to estimate
%       h = bandwidth
% --------------------------------------------

nobs = size(X,1);
A_h = (c-X)/h;
B_h = A_h.^2;
k2 =1/5;
const = 3/4;
ExpA_h = (-1)*const*sum(A_h.*(1-B_h).*(abs(A_h)<1))/k2;
f1hat = ExpA_h/nobs/h^2;

%------------------------------------------------------------------
% End private function density1d_est
%------------------------------------------------------------------

function [C_vp] = Const_LPR_U_est(v,p)

% S = (\mu_{j+k}), 0<=j,k<=p, defined in FG(96,p61), p=6
S5 = [1/2 1/4 1/6 1/8 1/10;
       1/4 1/6 1/8 1/10 1/12;
       1/6 1/8 1/10 1/12 1/14;
       1/8 1/10 1/12 1/14 1/16;
       1/10 1/12 1/14 1/16 1/18];
% Sstar = (\nu_{j+k}), 0<=j,k<=p, defined in FG(96,p61), p=6
K5 = [1/4 1/8 1/12 1/16 1/20;
        1/8 1/12 1/16 1/20 1/24;
        1/12 1/16 1/20 1/24 1/28;
        1/16 1/20 1/24 1/28 1/32;
        1/20 1/24 1/28 1/32 1/36];
        
% Vector of \mu_j,
% [\mu_1, \mu_2, ..., \mu_{13}]
mu_vec = [1/4 1/6 1/8 1/10 1/12 1/14 1/16 1/18 1/20];

Sinv_p = invpd(S5(1:p+1,1:p+1));

% int t^{p+1}K_{v,p}^*(t)dt
int_tp_K = Sinv_p(v+1,:)*mu_vec(1,(p+1):(2*p+1))';

% \int_K_{v,p}^{*2}(t)dt
int_K2 = Sinv_p(v+1,:)*K5(1:p+1,1:p+1)*Sinv_p(v+1,:)';

C_vp = ((factorial(p+1)^2*(2*v+1)*int_K2)/(2*(p+1-v)*int_tp_K^2))^(1/(2*p+3));

%------------------------------------------------------------------
% End private function Const_LPR_U_est
%------------------------------------------------------------------

function [Qn] = Qn(h,f_pilot,f1_pilot,m2_pilot,m3_pilot,v_pilot,n)

% matrix for constants depending on kernel
% [b1 c1 c2 v]
k_const = [-0.1 -0.1 -0.08 24/5];

b2 = zeros(1,2);
m2foverf12 = zeros(1,2);
for j=1:2
    m2foverf12(j) = m2_pilot(j)*f1_pilot/f_pilot/2;
    b2(j) = (-1)^(j+1)*(k_const(2)*(m2foverf12(j) + m3_pilot(j)/6)...   
         - k_const(3)*m2foverf12(j));
end
    
Qn = k_const(1)^2*(m2_pilot(1)*h(1)^2 - m2_pilot(2)*h(2)^2)^2/4 ...
+(b2(1)*h(1)^3 - b2(2)*h(2)^3)^2 ...
+ k_const(4)*(v_pilot(1)/h(1) + v_pilot(2)/h(2))/f_pilot/n;

%------------------------------------------------------------------
% End private function Qn
%------------------------------------------------------------------

function [Qn] = Qn_FRDD(h,f_pilot,f1_pilot,m2_pilot,m3_pilot,v_pilot,tau,n)

m2_pilot_y = m2_pilot(1:2);
m2_pilot_t = m2_pilot(3:4);
m3_pilot_y = m3_pilot(1:2);
m3_pilot_t = m3_pilot(3:4);
v_pilot_y = v_pilot(1:2);
v_pilot_t = v_pilot(3:4);
cov_pilot = v_pilot(5:6);

% matrix for constants depending on kernel
k_const = [-0.1 -0.1 -0.08 24/5]; % [b1 xi_1 xi_2 v], C1 = b_1/2
m2foverf12_y = zeros(1,2);
m2foverf12_t = zeros(1,2);
omega = zeros(1,2);
phi = zeros(1,2);
psi = zeros(1,2);
zeta_y = zeros(1,2);
zeta_t = zeros(1,2);

for j=1:2
     m2foverf12_y(j) = m2_pilot_y(j)*f1_pilot/f_pilot/2;
     m2foverf12_t(j) = m2_pilot_t(j)*f1_pilot/f_pilot/2;
     zeta_y(j) = (-1)^(j+1)*(k_const(2)*(m2foverf12_y(j) + m3_pilot_y(j)/6)...   
         - k_const(3)*m2foverf12_y(j));
     zeta_t(j) = (-1)^(j+1)*(k_const(2)*(m2foverf12_t(j) + m3_pilot_t(j)/6)...   
         - k_const(3)*m2foverf12_t(j));
     phi(j) = k_const(1)/2*(m2_pilot_y(j) - tau*m2_pilot_t(j));
     psi(j) = zeta_y(j) - tau*zeta_t(j);
     omega(j) = v_pilot_y(j) + tau^2*v_pilot_t(j) - 2*tau*cov_pilot(j);
end

Qn = (phi(1)*h(1)^2 - phi(2)*h(2)^2)^2 + (psi(1)*h(1)^3 - psi(2)*h(2)^3)^2 ...
+ k_const(4)*(omega(1)/h(1) + omega(2)/h(2))/f_pilot/n;

%------------------------------------------------------------------
% End private function Qn_FRDD
%------------------------------------------------------------------

function [tau,se,n_p,n_m,tau_d] = frdd_est(Y,D,X,c,h_p,h_m)

N=size(Y,1);
Xx = X - c;
nones = ones(N,1);

data_mat = [Y nones Xx];

Xbin_bt_x = (X>=c);
Xbin_bet_xh = Xbin_bt_x.*(X<=c+h_p);     % index for X in [x, x+h2_p]
n_p = sum(Xbin_bet_xh);
data_trim_mat_p = data_mat(Xbin_bet_xh==1,:);
YT_p = data_trim_mat_p(:,1);
T_p = data_trim_mat_p(:,2:3);
Kvec_p = (1-abs(data_trim_mat_p(:,3))/h_p).*(abs(data_trim_mat_p(:,3))<h_p);
W_p = diag(Kvec_p);
lambda_p_all = (T_p'*W_p*T_p)\(T_p'*W_p*YT_p);
delta_p = lambda_p_all(1,1);
Resid_p = YT_p - T_p*lambda_p_all;
Denom_p = trace(W_p - W_p*T_p*invpd(T_p'*W_p*T_p)*T_p'*W_p);
Sigma_p = sum(Resid_p.^2.*Kvec_p)/Denom_p;
Gam_inv_p = invpd(T_p'*W_p*T_p);
V_comp_p = Gam_inv_p*(T_p'*W_p*W_p*T_p)*Gam_inv_p;
V_mat_p = Sigma_p*V_comp_p(1,1);

Xbin_lt_x = (X<c);
Xbin_bet_hx = Xbin_lt_x.*(X>=c-h_m);     % index for X in [x, x+h2_p]
n_m = sum(Xbin_bet_hx);
data_trim_mat_m = data_mat(Xbin_bet_hx==1,:);
YT_m = data_trim_mat_m(:,1);
T_m = data_trim_mat_m(:,2:3);
Kvec_m = (1-abs(data_trim_mat_m(:,3)/h_m)).*(abs(data_trim_mat_m(:,3))<h_m);
W_m = diag(Kvec_m);
lambda_m_all = (T_m'*W_m*T_m)\(T_m'*W_m*YT_m);
delta_m = lambda_m_all(1,1);
Resid_m = YT_m - T_m*lambda_m_all;
Denom_m = trace(W_m - W_m*T_m*invpd(T_m'*W_m*T_m)*T_m'*W_m);
Sigma_m = sum(Resid_m.^2.*Kvec_m)/Denom_m;
Gam_inv_m = invpd(T_m'*W_m*T_m);
V_comp_m = Gam_inv_m*(T_m'*W_m*W_m*T_m)*Gam_inv_m;
V_mat_m = Sigma_m*V_comp_m(1,1);

data_mat2 = [D nones Xx];

data_trim_mat2_p = data_mat2(Xbin_bet_xh==1,:);
DT_p = data_trim_mat2_p(:,1);
lambda_p_all_d = (T_p'*W_p*T_p)\(T_p'*W_p*DT_p);
delta_p_d = lambda_p_all_d(1,1);
Resid_p_d = DT_p - T_p*lambda_p_all_d;
Sigma_p_d = sum(Resid_p_d.^2.*Kvec_p)/Denom_p;
V_mat_p_d = Sigma_p_d*V_comp_p(1,1);
Cov_p = sum(Resid_p'*W_p*Resid_p_d)/Denom_p*V_comp_p(1,1);

data_trim_mat2_m = data_mat2(Xbin_bet_hx==1,:);
DT_m = data_trim_mat2_m(:,1);
lambda_m_all_d = (T_m'*W_m*T_m)\(T_m'*W_m*DT_m);
delta_m_d = lambda_m_all_d(1,1);
Resid_m_d = YT_m - T_m*lambda_m_all_d;
Sigma_m_d = sum(Resid_m_d.^2.*Kvec_m)/Denom_m;
V_mat_m_d = Sigma_m_d*V_comp_m(1,1);
Cov_m = sum(Resid_m'*W_m*Resid_m_d)/Denom_m*V_comp_m(1,1);

tau_y = delta_p - delta_m;
tau_d = delta_p_d - delta_m_d;
tau = tau_y/tau_d;

omega_p = V_mat_p + tau^2*V_mat_p_d - 2*tau*Cov_p;
omega_m = V_mat_m + tau^2*V_mat_m_d - 2*tau*Cov_m;
variance = (omega_p + omega_m)/tau_d^2;
se = sqrt(variance);

%------------------------------------------------------------------
% End private function frdd_est
%------------------------------------------------------------------