% input reduced original EVS data
% data : size 49,720 times 37

n=length(data(:,1));

% calender time

%J_q4 - J_q4 : quarter 1 - 4

J_q1=zeros(n,1);
J_q2=J_q1;
J_q3=J_q2;
J_q4=J_q3;

I_q1=find(data(:,1)==1);
I_q2=find(data(:,1)==2);
I_q3=find(data(:,1)==3);
I_q4=find(data(:,1)==4);

J_q1(I_q1)=1;
J_q2(I_q2)=1;
J_q3(I_q3)=1;
J_q4(I_q4)=1;

clear I_q1 I_q2 I_q3 I_q4

%reduce data to summer quarters 

I=find( J_q2==1 | J_q3==1 ); 
data=data(I,:);
clear I J_q1 J_q2 J_q3 J_q4

n=length(data(:,1));

% employment status

I_emp=find(data(:,2)==1);
J_emp=zeros(n,1);

J_emp(I_emp)=1;
J_ne=ones(n,1);
J_ne=J_ne-J_emp;

clear I_emp

% separating households into demographic subgroubs
% J_s0 single without child
% J_s1 single with 1 child
% J_c0 couple without...
% J_c1 couple with 1 child
% J_c2 couple with two children
% J_c3 couple with three children

I_s0=find(data(:,15)==1 | data(:,15)==2);
I_s1=find(data(:,15)==3);
I_c0=find(data(:,15)==7 | data(:,15)==8 | data(:,15)==20 | data(:,15)==21);
I_c1=find(data(:,15)==9 | data(:,15)==10 | data(:,15)==22 | data(:,15)==23);
I_c2=find(data(:,15)==11 | data(:,15)==12 | data(:,15)==24 | data(:,15)==24);
I_c3=find(data(:,15)==13 | data(:,15)==14);

J_s0=zeros(n,1);
J_s1=J_s0;
J_c0=J_s0;
J_c1=J_s0;
J_c2=J_s0;
J_c3=J_s0;

J_s0(I_s0)=1;
J_s1(I_s1)=1;
J_c0(I_c0)=1;
J_c1(I_c1)=1;
J_c2(I_c2)=1;
J_c3(I_c3)=1;

clear I_s0 I_s1 I_c0 I_c1 I_c2 I_c3



% columns 3 - 14 : original commodity groups
% columns 16 - 36: durables to be substracted

% create aggregated non-durables

s(:,1)=data(:,3);
s(:,2)=data(:,4);
s(:,3)=data(:,5);
s(:,4)=data(:,6);
s(:,5)=data(:,7)-sum(data(:,16:25),2);
s(:,6)=data(:,8)-data(:,26);
s(:,7)=data(:,9)-sum(data(:,27:30),2);
s(:,8)=data(:,10);
s(:,9)=data(:,11)-sum(data(:,31:34),2)-data(:,35);
s(:,10)=data(:,12);
s(:,11)=data(:,13);
s(:,12)=data(:,14)-data(:,36);

% compute total expenditure
x=sum(data(:,3:14),2);

for i=1:12
s(:,i)=s(:,i)./x;
end

lnx=log(x);

clear i

% indices for demographic groups
I_s0=find( J_s0==1 ); 
I_s1=find( J_s1==1 );
I_c0=find( J_c0==1 );
I_c1=find( J_c1==1 ); 
I_c2=find( J_c2==1 );
I_c3=find( J_c3==1 );

I_s0_emp=find( J_s0==1 & J_emp==1 ); 
I_s1_emp=find( J_s1==1 & J_emp==1 );
I_c0_emp=find( J_c0==1 & J_emp==1 );
I_c1_emp=find( J_c1==1 & J_emp==1 ); 
I_c2_emp=find( J_c2==1 & J_emp==1 );
I_c3_emp=find( J_c3==1 & J_emp==1 );

I_s0_ne=find( J_s0==1 & J_ne==1 ); 
I_s1_ne=find( J_s1==1 & J_ne==1 );
I_c0_ne=find( J_c0==1 & J_ne==1 );
I_c1_ne=find( J_c1==1 & J_ne==1 ); 
I_c2_ne=find( J_c2==1 & J_ne==1 );
I_c3_ne=find( J_c3==1 & J_ne==1 );

q=0.25;
p=quantile(data(I_s0,37),q);
J=find(data(I_s0,37)>p);
I_s0_q25=I_s0;
I_s0_q25(J)=[];

p=quantile(data(I_s1,37),q);
J=find(data(I_s1,37)>p);
I_s1_q25=I_s1;
I_s1_q25(J)=[];

p=quantile(data(I_c0,37),q);
J=find(data(I_c0,37)>p);
I_c0_q25=I_c0;
I_c0_q25(J)=[];

p=quantile(data(I_c1,37),q);
J=find(data(I_c1,37)>p);
I_c1_q25=I_c1;
I_c1_q25(J)=[];

p=quantile(data(I_c2,37),q);
J=find(data(I_c2,37)>p);
I_c2_q25=I_c2;
I_c2_q25(J)=[];

p=quantile(data(I_c3,37),q);
J=find(data(I_c3,37)>p);
I_c3_q25=I_c3;
I_c3_q25(J)=[];

q=0.75;
p=quantile(data(I_s0,37),q);
J=find(data(I_s0,37)<=p);
I_s0_q75=I_s0;
I_s0_q75(J)=[];

p=quantile(data(I_s1,37),q);
J=find(data(I_s1,37)<=p);
I_s1_q75=I_s1;
I_s1_q75(J)=[];

p=quantile(data(I_c0,37),q);
J=find(data(I_c0,37)<=p);
I_c0_q75=I_c0;
I_c0_q75(J)=[];

p=quantile(data(I_c1,37),q);
J=find(data(I_c1,37)<=p);
I_c1_q75=I_c1;
I_c1_q75(J)=[];

p=quantile(data(I_c2,37),q);
J=find(data(I_c2,37)<=p);
I_c2_q75=I_c2;
I_c2_q75(J)=[];

p=quantile(data(I_c3,37),q);
J=find(data(I_c3,37)<=p);
I_c3_q75=I_c3;
I_c3_q75(J)=[];

q=0.25;
p=quantile(data(I_s0_emp,37),q);
J=find(data(I_s0_emp,37)>p);
I_s0_emp_q25=I_s0_emp;
I_s0_emp_q25(J)=[];

p=quantile(data(I_s1_emp,37),q);
J=find(data(I_s1_emp,37)>p);
I_s1_emp_q25=I_s1_emp;
I_s1_emp_q25(J)=[];

p=quantile(data(I_c0_emp,37),q);
J=find(data(I_c0_emp,37)>p);
I_c0_emp_q25=I_c0_emp;
I_c0_emp_q25(J)=[];

p=quantile(data(I_c1_emp,37),q);
J=find(data(I_c1_emp,37)>p);
I_c1_emp_q25=I_c1_emp;
I_c1_emp_q25(J)=[];

p=quantile(data(I_c2_emp,37),q);
J=find(data(I_c2_emp,37)>p);
I_c2_emp_q25=I_c2_emp;
I_c2_emp_q25(J)=[];

p=quantile(data(I_c3_emp,37),q);
J=find(data(I_c3_emp,37)>p);
I_c3_emp_q25=I_c3_emp;
I_c3_emp_q25(J)=[];

q=0.75;
p=quantile(data(I_s0_emp,37),q);
J=find(data(I_s0_emp,37)<=p);
I_s0_emp_q75=I_s0_emp;
I_s0_emp_q75(J)=[];

p=quantile(data(I_s1_emp,37),q);
J=find(data(I_s1_emp,37)<=p);
I_s1_emp_q75=I_s1_emp;
I_s1_emp_q75(J)=[];

p=quantile(data(I_c0_emp,37),q);
J=find(data(I_c0_emp,37)<=p);
I_c0_emp_q75=I_c0_emp;
I_c0_emp_q75(J)=[];

p=quantile(data(I_c1_emp,37),q);
J=find(data(I_c1_emp,37)<=p);
I_c1_emp_q75=I_c1_emp;
I_c1_emp_q75(J)=[];

p=quantile(data(I_c2_emp,37),q);
J=find(data(I_c2_emp,37)<=p);
I_c2_emp_q75=I_c2_emp;
I_c2_emp_q75(J)=[];

p=quantile(data(I_c3_emp,37),q);
J=find(data(I_c3_emp,37)<=p);
I_c3_emp_q75=I_c3_emp;
I_c3_emp_q75(J)=[];

clear J

clear data J_s0 J_s1 J_c0 J_c1 J_c2 J_c3 J_emp J_ne p