clear all;
diary qeffectsfinal.out;
disp('                   QEFFECTSFINAL.OUT');
disp(' ');

% LOAD DATA AND CONSTRUCT VARIABLES

load jtpa.raw; 

sex = jtpa(:,5);
male = find(sex==1);
female = find(sex==0);

ym = jtpa(male,2);
zm = jtpa(male,3);
dm = jtpa(male,4);
om = ones(size(ym));
xm = jtpa(male,[6 7 8 9 10 17 18 12 13 14 15 16 19]);
nm =length(male);

% Covariates for men: hsorged   black    hispanic   married   wkless13 
%                     class_tr  ojt_jsa  age2225    age2629    age3035 
%                     age3644   age4554  f2sms.  

yf = jtpa(female,2);
zf = jtpa(female,3);
df = jtpa(female,4);
of = ones(size(yf));
xf = jtpa(female,[6 7 8 9 10 11 17 18 12 13 14 15 16 19]);
nf =length(female);

% Covariates for women: hsorged   black    hispanic   married   wkless13 
%                       afdc      class_tr ojt_jsa    age2225   age2629    
%                       age3035   age3644  age4554    f2sms.  


clear jtpa sex;
pack;

% CREATE LABELS

ls   =['training';
       'constant'];

lcovm =['training';
       'hsorged ';
       'black   ';
       'hispanic';
       'married ';
       'wkless13';
       'constant'];

lcovf =['training';
       'hsorged ';
       'black   ';
       'hispanic';
       'married ';
       'wkless13';
       'afdc    ';
       'constant'];



Q = [0.15 0.25 0.50 0.75 0.85];


disp('WITHOUT COVARIATES');
disp(' ');
disp('      QUANTILE REGRESSION');
disp(' ');
disp('FOR MEN');
disp(' ');
for i=1:length(Q);
 q=Q(i);
 disp(sprintf('Quantile %1.2f :\n',q));
 [b,se]=qregkb(ym,[om dm],q);
 disp('variable        coeff.          s.e.');
 disp('------------------------------------')
 disp([ls num2str([b(2);b(1)],'\f %12.4f') num2str([se(2);se(1)],'\f %12.4f')])
 disp('------------------------------------')
% disp(sprintf('Effect of training in percentage: %4.2f%%', 100*b(2)/b(1)));
 disp(' ');
end;

disp('FOR WOMEN');
disp(' ');
for i=1:length(Q);
 q=Q(i);
 disp(sprintf('Quantile %1.2f :\n',q)); 
 [b,se]=qregkb(yf,[of df],q);
 disp('variable        coeff.          s.e.');
 disp('------------------------------------')
 disp([ls num2str([b(2);b(1)],'\f %12.4f') num2str([se(2);se(1)],'\f %12.4f')])
 disp('------------------------------------')
% disp(sprintf('Effect of training in percentage: %4.2f%%', 100*b(2)/b(1)));
 disp(' ');
end;

disp(' ')
disp('     QUANTILE TREATMENT EFFECTS');
disp(' ')
disp('FOR MEN');
disp(' ')

tau = mean(zm);

% Cross validation without covariates
CVm=[];
X = [dm (1-dm)];
K=6;
sc = 10000;
for k=1:K
 X = [X ((ym/sc).^k).*dm ((ym/sc).^k).*(1-dm)];
 S=inv(X'*X);
 b=S*(X'*zm);
 u=(zm-X*b)./(1-sum((X*S)'.*X'))';
 CVm=[CVm sum(u.^2)];
end




K=5;
sc = 10000;
X = [dm (1-dm)];
for k=1:K
 X = [X ((ym/sc).^k).*dm ((ym/sc).^k).*(1-dm)];
end
S=inv(X'*X);
b=S*(X'*zm);
v = X*b;
clear X;


kappa = 1 - (dm.*(1-v)./(1-tau)) - ((1-dm).*v./tau);
kappa = kappa.*(kappa>0);
kappa_z = 1 - (dm.*(1-zm)./(1-tau)) - ((1-dm).*zm./tau);
for i=1:length(Q);
 q=Q(i);
 disp(sprintf('Quantile %1.2f :\n',q))
 b=quantlsf(spdiags(kappa,0,nm,nm)*[om dm],kappa.*ym,q);
 e = ym - [om dm]*b;
 c = (nm^(-1/5))*std(e);
 t = (15/(16*c))*((1-((e/c).^2)).^2).*(abs((e/c))<1);
 J = [om dm]'*spdiags(kappa.*t,0,nm,nm)*[om dm];
 der = ((1-dm).*zm/(tau^2)) - (dm.*(1-zm)/((1-tau)^2));
 M_pi = mean(spdiags((q-(e<0)).*der,0,nm,nm)*[om dm] )';
 S = [om dm]'*spdiags((kappa_z.*(q-(e<0))).^2,0,nm,nm)*[om dm] + ...
     sum(spdiags(kappa_z.*(q-(e<0)).*(zm-tau),0,nm,nm)*[om dm])'*M_pi' + ...
     M_pi* sum(spdiags(kappa_z.*(q-(e<0)).*(zm-tau),0,nm,nm)*[om dm]) + ...
     M_pi*sum((zm-tau).^2)*M_pi';
 V = inv(J)*S*inv(J);
 se = sqrt(diag(V));
 disp('variable        coeff.          s.e.');
 disp('------------------------------------')
 disp([ls num2str([b(2);b(1)],'\f %12.4f') num2str([se(2);se(1)],'\f %12.4f')])
 disp('------------------------------------')
% disp(sprintf('Effect of training in percentage: %4.2f%%', 100*b(2)/b(1)));
 disp(' ');  
end;

disp('FOR WOMEN');
disp(' ')
tau = mean(zf);

% Cross validation without covariates
CVf=[];
X = [df (1-df)];
K=7;
sc = 10000;
for k=1:K
 X = [X ((yf/sc).^k).*df ((yf/sc).^k).*(1-df)];
 S=inv(X'*X);
 b=S*(X'*zf);
 u=(zf-X*b)./(1-sum((X*S)'.*X'))';
 CVf=[CVf sum(u.^2)] 
end



X = [df (1-df)];
K=6;
sc = 10000;
for k=1:K
 X = [X ((yf/sc).^k).*df ((yf/sc).^k).*(1-df)];
end


S=inv(X'*X);
b=S*(X'*zf);
v = X*b;
clear X;

kappa = 1 - (df.*(1-v)/(1-tau)) - ((1-df).*v/tau);
kappa = kappa.*(kappa>0);
kappa_z = 1 - (df.*(1-zf)/(1-tau)) - ((1-df).*zf/tau);
for i=1:length(Q);
 q=Q(i);
 disp(sprintf('Quantile %1.2f :\n',q)); 
 b=quantlsf(spdiags(kappa,0,nf,nf)*[of df],kappa.*yf,q);
 e = yf - [of df]*b;
 c = (nf^(-1/5))*std(e);
 t = (15/(16*c))*((1-((e/c).^2)).^2).*(abs((e/c))<1);
 J = [of df]'*spdiags(kappa.*t,0,nf,nf)*[of df];
 der = ((1-df).*zf/(tau^2)) - (df.*(1-zf)/((1-tau)^2));
 M_pi = mean(spdiags((q-(e<0)).*der,0,nf,nf)*[of df] )';
 S = [of df]'*spdiags((kappa_z.*(q-(e<0))).^2,0,nf,nf)*[of df] + ...
     sum(spdiags(kappa_z.*(q-(e<0)).*(zf-tau),0,nf,nf)*[of df])'*M_pi' + ...
     M_pi* sum(spdiags(kappa_z.*(q-(e<0)).*(zf-tau),0,nf,nf)*[of df]) + ...
     M_pi*sum((zf-tau).^2)*M_pi';
 V = inv(J)*S*inv(J);
 se = sqrt(diag(V));
 disp('variable        coeff.          s.e.');
 disp('------------------------------------')
 disp([ls num2str([b(2);b(1)],'\f %12.4f') num2str([se(2);se(1)],'\f %12.4f')])
 disp('------------------------------------')
% disp(sprintf('Effect of training in percentage: %4.2f%%', 100*b(2)/b(1)));
 disp(' ');
end;




disp('WITH COVARIATES');
disp(' ');
disp('      QUANTILE REGRESSION');
disp(' ');
disp('FOR MEN');
disp(' ');
for i=1:length(Q);
 q=Q(i);
 disp(sprintf('Quantile %1.2f :\n',q));
 [b,se]=qregkb(ym,[om dm xm],q);
 disp('variable        coeff.          s.e.');
 disp('------------------------------------')
 disp([lcovm num2str(b([2:size(lcovm,1) 1]),'\f %12.4f')...
       num2str(se([2:size(lcovm,1) 1]),'\f %12.4f')])
 disp('------------------------------------')
 disp(sprintf('Y0: %4.2f', ...
              ([1 0 mean(xm(find(dm==1),:))]*b)));
% disp(sprintf('Effect of training in percentage: %4.2f%%', ...
%              100*b(2)/([1 0 mean(xm(find(dm==1),:))]*b)));
 disp(' '); 
end;

disp('FOR WOMEN');
disp(' ');
for i=1:length(Q);
 q=Q(i);
 disp(sprintf('Quantile %1.2f :\n',q)); 
 [b,se]=qregkb(yf,[of df xf],q);
 disp('variable        coeff.          s.e.');
 disp('------------------------------------')
 disp([lcovf num2str(b([2:size(lcovf,1) 1]),'\f %12.4f')...
       num2str(se([2:size(lcovf,1) 1]),'\f %12.4f')])
 disp('------------------------------------')
 disp(sprintf('Y0: %4.2f', ...
              ([1 0 mean(xf(find(df==1),:))]*b)));
% disp(sprintf('Effect of training in percentage: %4.2f%%', ...
%              100*b(2)/([1 0 mean(xf(find(df==1),:))]*b)));
 disp(' ');  
end;






disp(' ')
disp('     QUANTILE TREATMENT EFFECTS');
disp(' ')
disp('FOR MEN');
disp(' ')

tau = mean(zm);

K=5;
sc = 10000;
X = [dm (1-dm)];
for k=1:K
 X = [X ((ym/sc).^k).*dm ((ym/sc).^k).*(1-dm)];
end
S=inv(X'*X);
b=S*(X'*zm);
v = X*b;
clear X;


kappa = 1 - (dm.*(1-v)./(1-tau)) - ((1-dm).*v./tau);
kappa = kappa.*(kappa>0);
kappa_z = 1 - (dm.*(1-zm)./(1-tau)) - ((1-dm).*zm./tau);
for i=1:length(Q);
 q=Q(i);
 disp(sprintf('Quantile %1.2f :\n',q))
 b=quantlsf(spdiags(kappa,0,nm,nm)*[om dm xm],kappa.*ym,q);
 e = ym - [om dm xm]*b;
 c = (nm^(-1/5))*std(e);
 t = (15/(16*c))*((1-((e/c).^2)).^2).*(abs((e/c))<1);
 J = [om dm xm]'*spdiags(kappa.*t,0,nm,nm)*[om dm xm];
 der = ((1-dm).*zm/(tau^2)) - (dm.*(1-zm)/((1-tau)^2));
 M_pi = mean(spdiags((q-(e<0)).*der,0,nm,nm)*[om dm xm] )';
 S = [om dm xm]'*spdiags((kappa_z.*(q-(e<0))).^2,0,nm,nm)*[om dm xm] + ...
     sum(spdiags(kappa_z.*(q-(e<0)).*(zm-tau),0,nm,nm)*[om dm xm])'*M_pi' + ...
     M_pi* sum(spdiags(kappa_z.*(q-(e<0)).*(zm-tau),0,nm,nm)*[om dm xm]) + ...
     M_pi*sum((zm-tau).^2)*M_pi';
 V = inv(J)*S*inv(J);
 se = sqrt(diag(V));
 
 disp('variable        coeff.          s.e.');
 disp('------------------------------------')
 disp([lcovm num2str(b([2:size(lcovm,1) 1]),'\f %12.4f')...
       num2str(se([2:size(lcovm,1) 1]),'\f %12.4f')])
 disp('------------------------------------')
% disp(sprintf('Effect of training in percentage: %4.2f%%', ...
%              100*b(2)/([1 0 mean(xm(find(dm==1),:))]*b)));
 disp(sprintf('Y0: %4.2f', ...
              ([1 0 mean(xm(find(dm==1),:))]*b)));
 disp(' '); 

end;

disp('FOR WOMEN');
disp(' ')
tau = mean(zf);


X = [df (1-df) df.*xf(:,7) (1-df).*xf(:,7)];
K=3;
sc = 10000;
for k=1:K
 X = [X ((yf/sc).^k).*df ((yf/sc).^k).*(1-df) ((yf/sc).^k).*df.*xf(:,7) ((yf/sc).^k).*(1-df).*xf(:,7) ];
end


S=inv(X'*X);
b=S*(X'*zf);
v = X*b;
clear X;

kappa = 1 - (df.*(1-v)/(1-tau)) - ((1-df).*v/tau);
kappa = kappa.*(kappa>0);
kappa_z = 1 - (df.*(1-zf)/(1-tau)) - ((1-df).*zf/tau);
for i=1:length(Q);
 q=Q(i);
 disp(sprintf('Quantile %1.2f :\n',q)); 
 b=quantlsf(spdiags(kappa,0,nf,nf)*[of df xf],kappa.*yf,q);
 e = yf - [of df xf]*b;
 c = (nf^(-1/5))*std(e);
 t = (15/(16*c))*((1-((e/c).^2)).^2).*(abs((e/c))<1);
 J = [of df xf]'*spdiags(kappa.*t,0,nf,nf)*[of df xf];
 der = ((1-df).*zf/(tau^2)) - (df.*(1-zf)/((1-tau)^2));
 M_pi = mean(spdiags((q-(e<0)).*der,0,nf,nf)*[of df xf] )';
 S = [of df xf]'*spdiags((kappa_z.*(q-(e<0))).^2,0,nf,nf)*[of df xf] + ...
     sum(spdiags(kappa_z.*(q-(e<0)).*(zf-tau),0,nf,nf)*[of df xf])'*M_pi' + ...
     M_pi* sum(spdiags(kappa_z.*(q-(e<0)).*(zf-tau),0,nf,nf)*[of df xf]) + ...
     M_pi*sum((zf-tau).^2)*M_pi';
 V = inv(J)*S*inv(J);
 se = sqrt(diag(V));
 
 disp('variable        coeff.          s.e.');
 disp('------------------------------------')
 disp([lcovf num2str(b([2:size(lcovf,1) 1]),'\f %12.4f')...
       num2str(se([2:size(lcovf,1) 1]),'\f %12.4f')])
 disp('------------------------------------')
 disp(sprintf('Y0: %4.2f', ...
              ([1 0 mean(xf(find(df==1),:))]*b)));
% disp(sprintf('Effect of training in percentage: %4.2f%%', ...
%              100*b(2)/([1 0 mean(xf(find(df==1),:))]*b)));
 disp(' '); 
 
end;

diary off;