clear;
clc;
close;

%% load results
cd results
load gas_data_final
cd ..

%% We compute parameters in the real data for June 2010.
% this is necessary for the Monte Carlo simulations;
real_data = consumption(:,(day>=517).*(day<517 + 30)==1);
real_hour = hour((day>=517).*(day<517 + 30)==1);
real_mean = NaN(48,1);
real_sd = NaN(48,1);
real_corr = NaN(48,1);
for a=1:48
    real_mean(a) = mean(mean(real_data(:,real_hour==a)));
    real_sd(a) = mean(std(real_data(:,real_hour==a))); 
end
for a=1:48
    aux_data = [];
    for i=1:size(real_data,1)
        for j=2:size(real_data,2)
            if real_hour(j)==a
                aux_data = [aux_data;[real_data(i,j-1),real_data(i,j)]]; %#ok<AGROW>
            end
        end
    end 
    aux_data2 = corrcoef(aux_data(:,1),aux_data(:,2));
    real_corr(a) = aux_data2(1,2);
end

cd results
save gas_data_OneDay real_mean real_sd real_corr;
cd ..


%% Compute statistics by group
clear
cd results
load gas_data_final
cd ..

% sample size;
n = size(consumption,1);
treat_indicator = treatment>0;
experiment = day >= 517; % 517 day is June 1, 2010;
% monthly pattern: Jan 2009 and onwards;
month_start = [342,366,397,425,456,485,517,547,578,609,639,670,700,731,762,790,821,851]';

% mean and standard deviation
month_mean2 = nan(size(month_start,1)-1,6);
month_sd2 = nan(size(month_start,1)-1,6);
month_zeros2 = nan(size(month_start,1)-1,6); 

for j=1:size(month_start,1)-1
    for a=0:max(treatment)
    aux_matrix = consumption(treatment==a,(day>=month_start(j)).*(day<month_start(j+1))==1);    
    aux_matrix = reshape(aux_matrix,size(aux_matrix,1)*size(aux_matrix,2),1);
    month_mean2(j,a+1) = mean(aux_matrix);
    month_sd2(j,a+1) = std(aux_matrix);
    month_zeros2(j,a+1) = mean(aux_matrix==0);
    end
    
    aux_matrix = consumption(treatment>0,(day>=month_start(j)).*(day<month_start(j+1))==1);
    aux_matrix = reshape(aux_matrix,size(aux_matrix,1)*size(aux_matrix,2),1);
    month_mean2(j,6) = mean(aux_matrix);
    month_sd2(j,6) = std(aux_matrix);
    month_zeros2(j,6) = mean(aux_matrix==0);
end

% correlations
% % Note: the correlations below take a long time.
% corr_vector = NaN(size(consumption,2)-1,6);
% 
% for j=2:size(consumption,2)
%     for a=0:max(treatment)
%     aux_data = consumption(treatment==a,:);
%     aux = corrcoef(aux_data(:,j-1),aux_data(:,j));
%     corr_vector(j-1,a+1)=aux(2,1);
%     end
%     
%     aux_data = consumption(treatment>0,:);
%     aux = corrcoef(aux_data(:,j-1),aux_data(:,j));
%     corr_vector(j-1,6)=aux(2,1);
% end
% cd results
% save gas_data_corr_vector corr_vector;
% cd ..

cd results
    load gas_data_corr_vector corr_vector;
cd ..

month_corr_group = nan(size(month_start,1)-1,size(corr_vector,2));
for j=1:size(month_start,1)-1
    for a=1:size(corr_vector,2)
    aux_matrix = corr_vector((day>=month_start(j)).*(day<month_start(j+1))==1,a);
    aux_matrix = reshape(aux_matrix,size(aux_matrix,1)*size(aux_matrix,2),1);
    month_corr_group(j,a) = mean(aux_matrix);
    end
end

%% Generate graphs in the paper
% Code to generate Figure 1
close
hFig = figure(1);
plot(month_mean2(:,1),'k')
hold on;
plot(month_mean2(:,2),'r')
plot(month_mean2(:,3),'g')
plot(month_mean2(:,4),'m')
plot(month_mean2(:,5),'b')
xlabel('Month of 2010')
ylabel('Consumption (kw)')
title('Average half-hour consumption')
legend('Control','Treatment 1','Treatment 2','Treatment 3','Treatment 4','Location','NorthEast')
axis([1 12 0 max(max(month_mean2(1:12,1:5)))+0.1])
grid on
% focus on months 6-12;
axis([6 12 0 max(max(month_mean2(1:12,1:5)))+0.1])
grid on
set(hFig, 'Position', [0 0 800 600])
cd results
saveas(hFig,'Figure1.png')
cd ..

% Code to generate Figure 2
close
hFig = figure(1);
plot(month_sd2(:,1),'k')
hold on;
plot(month_sd2(:,2),'r')
plot(month_sd2(:,3),'g')
plot(month_sd2(:,4),'m')
plot(month_sd2(:,5),'b')
xlabel('Month of 2010')
ylabel('Consumption (kw)')
title('Sample standard deviation of half-hour consumption')
legend('Control','Treatment 1','Treatment 2','Treatment 3','Treatment 4','Location','NorthEast')
axis([1 12 0 max(max(month_sd2(1:12,1:5)))+0.1])
grid on
% focus on months 6-12;
axis([6 12 0 max(max(month_sd2(1:12,1:5)))+0.1])
grid on
set(hFig, 'Position', [0 0 800 600])
cd results
saveas(hFig,'Figure2.png')
cd ..

% Code to generate Figure 3
close
hFig = figure(1);
plot(month_corr_group(:,1),'k')
hold on;
plot(month_corr_group(:,2),'r')
plot(month_corr_group(:,3),'g')
plot(month_corr_group(:,4),'m')
plot(month_corr_group(:,5),'b')
xlabel('Month of 2010')
ylabel('Consumption (kw)')
title('Sample correlation of consumption between consecutive half-hours')
legend('Control','Treatment 1','Treatment 2','Treatment 3','Treatment 4','Location','NorthEast')
axis([1 12 min(min(month_corr_group(1:12,1:5)))-0.2 max(max(month_corr_group(1:12,1:5)))+0.2])
grid on
set(hFig, 'Position', [0 0 800 600])
% focus on months 6-12;
axis([6 12 min(min(month_corr_group(1:12,1:5)))-0.2 max(max(month_corr_group(1:12,1:5)))+0.2])
set(hFig, 'Position', [0 0 800 600])
cd results
saveas(hFig,'Figure3.png')
cd ..