%% Monte Carlo Simulations - Canen and Song (2021)

clear
clc
rng(12123712) %set seed.

S = 1000; %Simulation number
beta = 0.1449;

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

if Specification1 == 1
id_set = [0,1.04]; %ID Set (Table 1)
LRR_set = [0,0.666];%LRR results (Table 1)
else
id_set = [0,1.637]; %ID Set (Table 1 - Spec 2)
LRR_set = [0,1.157];%LRR results (Table 1 - Spec 2)
end

%% Simulation 

%Random draws from set of parameters
beta_s_id = min(id_set)+(max(id_set)-min(id_set))*rand(S,1); %which parameter to draw in uniform case from id set
beta_s_LRR = min(LRR_set)+(max(LRR_set)-min(LRR_set))*rand(S,1); %which parameter to draw in uniform case from LRR
prob_mixed = rand(S,1); %mixed parameter draw: do we draw from uniform case or extreme case?

%Furthest-case scenario parameters: find which parameters in the LRR or ID
%95% set are furthest away from the true beta. Since the problem is 1
%dimensional, we only compare the extreme values of the ID set.
if max(id_set)-beta>beta-min(id_set)
id_furthest = max(id_set);
else
id_furthest = min(id_set);
end

if max(LRR_set)-beta>beta-min(LRR_set)
LRR_furthest = max(LRR_set);
else
LRR_furthest = min(LRR_set);
end

%Distribution of distances
dist_ID_uniform = abs(beta_s_id-beta);
dist_LRR_uniform = abs(beta_s_LRR-beta);
dist_ID_mixed = abs((prob_mixed>0.5).*abs(beta_s_id-beta)+(prob_mixed<=0.5)*id_furthest);
dist_LRR_mixed = abs((prob_mixed>0.5).*abs(beta_s_LRR-beta)+(prob_mixed<=0.5)*LRR_furthest);

%Store results (i.e. each statistic, across S simulations, and 3 types of
%draws of the parameters)
%uniformly drawn parameters
absdev_id(1) = mean(abs(beta_s_id-beta));
maxdev_id(1) = max(abs(beta_s_id-beta));
absdev_LRR(1) = mean(abs(beta_s_LRR-beta));
maxdev_LRR(1) = max(abs(beta_s_LRR-beta));

%Furthest-away case - largest beta away from the true one (and still within the set)
absdev_id(2) = abs(id_furthest-beta); %Average and Max are the same because beta and furthest draw are constant.
maxdev_id(2) = abs(id_furthest-beta);
absdev_LRR(2) = abs(LRR_furthest-beta);
maxdev_LRR(2) = abs(LRR_furthest-beta);
   
%mixed case scenario drawn parameters (50% uniform case, 50% furthest away)
%uniformly drawn parameters
absdev_id(3) = mean((prob_mixed>0.5).*abs(beta_s_id-beta)+(prob_mixed<=0.5)*abs(id_furthest-beta));
maxdev_id(3) = max((prob_mixed>0.5).*abs(beta_s_id-beta)+(prob_mixed<=0.5)*abs(id_furthest-beta));
absdev_LRR(3) = mean((prob_mixed>0.5).*abs(beta_s_LRR-beta)+(prob_mixed<=0.5)*abs(LRR_furthest-beta));
maxdev_LRR(3) = max((prob_mixed>0.5).*abs(beta_s_LRR-beta)+(prob_mixed<=0.5)*abs(LRR_furthest-beta));
   
%% Figures from Stability Simulations

%Deviation - Uniform draws of parameters from 95% CS
[f,xi] = ksdensity(dist_ID_uniform);
plot(xi,f,'k','LineWidth',2)
hold on
[f,xi] = ksdensity(dist_LRR_uniform);
plot(xi,f,':k','LineWidth',2)
legend('Identified Set', 'LRR')
ylabel('Density')
xlabel('Absolute Deviation - |\psi^s - \psi|')
set(gca, 'LineWidth',1.2, 'fontsize', 11, 'fontweight', 'bold');
xlim([0,2])

%Mean Deviation - (Uniform-case draws of parameters from 95% CS
disp(['Mean Abs. Deviation for (Uniform, Furthest Case, Mixed Case): (i) LRR = ', num2str(absdev_LRR), ', (ii) ID Set = ', num2str(absdev_id)])

%Maximum Deviation - Furthest-case draws of parameters from 95% CS
disp(['Maximum Abs. Deviation for (Uniform, Furthest Case, Mixed Case): (i) LRR = ', num2str(maxdev_LRR), ', (ii) ID Set = ', num2str(maxdev_id)])


%% Second Exercise
%Let us focus on the case in which our parameter of interest is the intercept
S = 1000;
pi = 0.25;

load StabilityExercise2_sets

%In this case, the LRR minimizes the INTERCEPT (i.e. \varphi_0).
id_set = [delta1_in_idset, gamma_ID_set]; %ID Set 
LRR_set = [repmat(delta1_LRR_kappa,size(gamma_ofdelta1_LRR,1),1), gamma_ofdelta1_LRR];

%Distribution of distances
id_draw_list = randi(size(id_set,1),S,1); %which parameter to draw in uniform case from id set
LRR_draw_list = randi(size(LRR_set,1),S,1); %which parameter to draw in uniform case from LRR

%Store vectors
dist_ID_uniform = zeros(S,1);
dist_LRR_uniform = zeros(S,1);

%Predicted data
for s = 1:S
    Y_avg = p_eta*(p_treatment*(delta_1+gamma)+(1-p_treatment)*delta_1)+(1-p_eta)*(p_treatment*log(delta_1+gamma)+(1-p_treatment)*log(delta_1));
    Y_ID_pred = pi*(p_treatment*(id_set(id_draw_list(s),1)+id_set(id_draw_list(s),2))+(1-p_treatment)*id_set(id_draw_list(s),1))+(1-pi)*(p_treatment*log(id_set(id_draw_list(s),1)+id_set(id_draw_list(s),2))+(1-p_treatment)*log(id_set(id_draw_list(s),1)));
    Y_LRR_pred = pi*(p_treatment*(LRR_set(LRR_draw_list(s),1)+LRR_set(LRR_draw_list(s),2))+(1-p_treatment)*LRR_set(LRR_draw_list(s),1))+(1-pi)*(p_treatment*log(LRR_set(LRR_draw_list(s),1)+LRR_set(LRR_draw_list(s),2))+(1-p_treatment)*log(LRR_set(LRR_draw_list(s),1)));

    %Distances
    dist_ID_uniform(s) = abs(Y_avg - Y_ID_pred);
    dist_LRR_uniform(s) = abs(Y_avg - Y_LRR_pred);
end

%Deviation - Uniform draws of parameters from 95% CS
[f,xi] = ksdensity(dist_ID_uniform);
plot(xi,f,'k','LineWidth',2)
hold on
[f,xi] = ksdensity(dist_LRR_uniform);
plot(xi,f,':k','LineWidth',2)
ylabel('Density')
xlabel('|Predicted ASO - True ASO|')
set(gca, 'LineWidth',1.2, 'fontsize', 11, 'fontweight', 'bold');
xlim([0,6])
ylim([0,4])
% legend('Identified Set', 'LRR')
