% Estimate STY's model

 clear
 
 path(path,'functions')

 % In order to bootstrap this I need to obtain evka and cevka inside code
 % from given data
 
data=csvread('Outputs/psidsample_boot.csv',1,1);
%data(:,4)=data(:,2).^2;


Nboot=1000;
xboot=zeros(4,Nboot);
xboot_u=zeros(4,Nboot);

for jboot=1:Nboot
    
evka= csvread('Outputs/earnvarprofile_age_kaplan.csv',1,1);
cevka= csvread('Outputs/earncovarprofile_age_kaplan.csv',1,1);

% Bootstrap over data

alea=floor(length(data)*rand(length(data),1))+1;
if jboot==1
alea=1:length(data);
end

data_boot=data(alea,:);

evma=zeros(36,1);
cevma=zeros(36,1);

for i=25:60
    ru=data_boot(:,2);
    ru2=data_boot(:,3);
    evma(i-24)=mean(ru(data_boot(:,1)==i));
    cevma(i-24)=mean(ru2(data_boot(:,1)==i));
end

evi=evka(2:37,2);
cevi=cevka(2:37,2);

myfun=@(x)sty_fun(x,evma,cevma);

options=optimoptions('fmincon','StepTolerance',1e-10);
x0=[0.01 0.12 0.01 0.99];
lb=[1e-4 0 0 0];
ub=[1 1 1 1];
[x,val] = fmincon(myfun,x0,[],[],[],[],lb,ub);
while norm(x-x0)>1e-4
 x0=x;
[x,val] = fmincon(myfun,x0,[],[],[],[],lb,ub);
end

[val2,fparts]=sty_fun(x,evma,cevma);

xboot(:,jboot)=x;

jboot

if jboot==1
    csvwrite('Outputs/variance_earn_canonical.csv',fparts)
end

end



for jboot=1:Nboot
    
% Bootstrap over data

alea=floor(length(data)*rand(length(data),1))+1;
if jboot==1
alea=1:length(data);
end

data_boot=data(alea,:);

evma=zeros(36,1);
cevma=zeros(36,1);

for i=25:60
    ru=data_boot(:,2);
    ru2=data_boot(:,3);
    evma(i-24)=mean(ru(data_boot(:,1)==i));
    cevma(i-24)=mean(ru2(data_boot(:,1)==i));
end

evi=evka(2:37,2);
cevi=cevka(2:37,2);

myfun=@(x)sty_fun(x,evma,cevma);

options=optimoptions('fmincon','StepTolerance',1e-10);
x0=[0.01 0.12 0.01 0.99];
lb=[1e-4 0 0 0];
ub=[1 1 1 2];
[x,val] = fmincon(myfun,x0,[],[],[],[],lb,ub);
while norm(x-x0)>1e-4
 x0=x;
[x,val] = fmincon(myfun,x0,[],[],[],[],lb,ub);
end

[val2,fparts]=sty_fun(x,evma,cevma);

xboot_u(:,jboot)=x;

jboot

end

% First result are the parameter estimates (bootstrap forced to take each
% observation only once)

disp(xboot(:,1))


% Of the constrained problem:

pc25=quantile(xboot,0.025,2);
pc75=quantile(xboot,0.975,2);

pc_u25=quantile(xboot_u,0.025,2);
pc_u75=quantile(xboot_u,0.975,2);