%% Monte Carlo Simulations - Canen and Song (2020)
%Based on LRR draft version B55, December 20 - 2020

clear
clc
rng(12123712) %set seed.

%Least Favorable Specification? (i.e. moment inequalities close to equalities)
Specification = 1; %1 or 2? (1 is close to moment equalities)

%Grid
grid_jump = 0.05; 
range_delta1 = 2.75;
range_gamma = 2.5;

%Parameters (Data)
p_eta = 0.5; %prob. of high type.
p_treatment = 0.3; %prob of treatment (~ what we observe in the data)
sigma_e = 1; %variance of error

%parameters (Linear Model)
%parameters fixed across specifications
delta_1 = 2.5;
gamma = 0.15; %e.g. ~5% return to college

theta = [delta_1;gamma;p_eta];
d = size(theta,1);

%True Parameter of interest:
beta = p_eta*gamma+(1-p_eta)*(log(delta_1 + gamma)-log(delta_1));

%bounds of Intervals
%Depending on spec (how much top coding)
if Specification==1
Z1 = 2.4; %top-coding in ~ 5% of obs.
Z2 = 7; %Under the parameters above, Z2 = 9 guarantees Y<=Z2 at all 1000000 draws.

%Smaller grid
range_delta1 = 6;
range_gamma = 5;

else
Z1 = 1.9;  %top-coding in ~20% of obs.
Z2 = 7;
end

%Grid for theta (has to be positive, make it large enough to search) 
delta1_grid = 0.01:grid_jump:delta_1+range_delta1; %all values must have lower grid above 0 due to log() parametrization.
gamma_grid = max(0,gamma-range_gamma):grid_jump:gamma+range_gamma;
p_eta_grid = 0.01:grid_jump*2:0.99;


%% Construct the Identified set:

%Set-up
n = 500000;
R = 1;

%Drawing Data
X = ones(1,n); %observables
D = (rand(n,1)<p_treatment);  %treatment occurs at around 10% of obs. (e.g. college), independent of X.
e = sigma_e*randn(n,R); %unobservables
eta = (rand(n,1)>p_eta);%unobserved ability: high in 50% of pop.

%Outcomes
Y_star = repmat((delta_1+gamma*D).*eta+log(delta_1+gamma*D).*(1-eta),1,R) + e; %if eta=1, linear reduced form; if eta = 0, log reduced form
Ztilde1 = Y_star.*(Y_star<=Z1)+Z1.*(Y_star>Z1); %top-coded Ztilde_1
Ztilde2 = Y_star.*(Y_star<=Z1)+Z2.*(Y_star>Z1); %topc coded Ztilde_2

%Indicator functions for D.
D_eq_1 = (D==1);
D_eq_0 = (D==0);
beta_ID_set = [];
delta1_in_idset = [];
gamma_ID_set = [];
%% Find ID Set

for m = 1:size(p_eta_grid,2)
for j = 1:size(delta1_grid,2)
for l = 1:size(gamma_grid,2)

m1 = (Ztilde1.*D_eq_1 - (p_eta_grid(m)*(delta1_grid(j)+gamma_grid(l)) + (1-p_eta_grid(m))*log(delta1_grid(j)+gamma_grid(l))).*D_eq_1);
m2 = (Ztilde1.*D_eq_0 - (p_eta_grid(m)*(delta1_grid(j)) + (1-p_eta_grid(m))*log(delta1_grid(j))).*D_eq_0);
m3 = (p_eta_grid(m)*(delta1_grid(j)+gamma_grid(l)) + (1-p_eta_grid(m))*log(delta1_grid(j)+gamma_grid(l))).*D_eq_1-Ztilde2.*D_eq_1;
m4 = (p_eta_grid(m)*(delta1_grid(j)) + (1-p_eta_grid(m))*log(delta1_grid(j))).*D_eq_0-Ztilde2.*D_eq_0;

%Average Moment.
mbar = [mean(m1),mean(m2),mean(m3),mean(m4)]; %mean over n.

if sum(mbar<=0)==4
beta_temp = p_eta_grid(m)*gamma_grid(l)+(1-p_eta_grid(m))*mean(log(delta1_grid(j)+gamma_grid(l))-log(delta1_grid(j)));
beta_ID_set = [beta_ID_set; beta_temp];
gamma_ID_set = [gamma_ID_set; gamma_grid(l)];
delta1_in_idset = [delta1_in_idset;delta1_grid(j)];
end
end
end
end

delta1_in_idset_unique = sort(unique(delta1_in_idset));
gamma_ID_set_unique = sort(unique(gamma_ID_set));

%% Find LRR set

kappa = 0.005;
delta1_LRR_kappa = [];
gamma_LRR_kappa = [];

%We are going to fix \tilde \beta and p_eta at its true value of 0.14..., and look at a
%certain value of \gamma...
for k = 1:size(gamma_ID_set_unique,1)
Q_LRR = zeros(size(delta1_in_idset_unique,2),1);
for j = 1:size(delta1_in_idset_unique,1)
Q_LRR(j) = mean(D)*(delta1_in_idset_unique(j)+gamma_ID_set_unique(k)-log(delta1_in_idset_unique(j)+gamma_ID_set_unique(k))).^2+(1-mean(D))*(delta1_in_idset_unique(j)-log(delta1_in_idset_unique(j))).^2;
end
delta_temp = delta1_in_idset_unique(Q_LRR <= min(Q_LRR)+kappa);
delta1_LRR_kappa = [delta1_LRR_kappa; delta_temp];
end

% LRR for delta_1. 
delta1_LRR_kappa = unique(delta1_LRR_kappa);
gamma_ofdelta1_LRR = [];

% Now find the values of gamma in the ID set consistent with this delta.
for j = 1:size(delta1_LRR_kappa,1) 
gamma_temp = unique(gamma_ID_set(delta1_LRR_kappa(j)==delta1_in_idset));
gamma_ofdelta1_LRR = [gamma_ofdelta1_LRR; gamma_temp];
end

save StabilityExercise2_sets p_eta p_treatment delta_1 delta1_LRR_kappa gamma gamma_ofdelta1_LRR delta1_in_idset gamma_ID_set
