function [rub]=panel_simul_cons(Qylife,copt,sopt,ylook,ylifelog,ye,lvarcons,...
    commonbet,choose,invher,yedtb)

global da dy TL a dye azero



%% Simulation of the process to find BPP coefficients

fast=1;

N=100000;
y_matrix=zeros(N,TL);
ye_matrix=zeros(N,TL);
c_matrix=zeros(N,TL);
a_matrix=zeros(N,TL);
rando=rand(N,TL);
rande=rand(N,TL);

invcum=cumsum(invher);

for i=1:N
    I=find(rando(i,1)<invcum,1,'first');
    y_matrix(i,1)=I;
    for t=1:TL-1
        Iprev=I;
        I=find(rando(i,t+1)<cumsum(Qylife(I,:,t)),1,'first');
        y_matrix(i,t+1)=I;
    end
end

cumue=zeros(dye,1);
cumue(1)=0.5*yedtb(1,1);
for i=2:dye
   cumue(i)=cumue(i-1)+0.5*yedtb(i,1)+0.5*yedtb(i-1,1);
end

if dye>1
for t=1:TL
ye_matrix(:,t)=interp1(cumue,1:dye,rande(:,t),'nearest','extrap');
end
else
   ye_matrix=ones(N,TL); 
end

% First:
ind=sub2ind([da dye dy 1 1 TL],azero*ones(N,1),ye_matrix(:,1),y_matrix(:,1),ones(N,1),ones(N,1),ones(N,1));

c_matrix(:,1)=copt(ind);
a_matrix(:,1)=sopt(ind);

for t=2:TL
    toa=interp1(a,1:da,a_matrix(:,t-1)-1e-6,'next','extrap');
    toa(toa==1)=2;
    toa2=toa-1;
    soa=(a(toa)-a_matrix(:,t-1))./(a(toa)-a(toa2)); % weight on previous gp
    ind=sub2ind([da dye dy 1 1 TL],toa,ye_matrix(:,t),y_matrix(:,t),ones(N,1),ones(N,1),t*ones(N,1));
    ind2=sub2ind([da dye dy 1 1 TL],toa2,ye_matrix(:,t),y_matrix(:,t),ones(N,1),ones(N,1),t*ones(N,1));
    
    c_matrix(:,t)=soa.*copt(ind2)+(1-soa).*copt(ind);
    a_matrix(:,t)=soa.*sopt(ind2)+(1-soa).*sopt(ind);
end


y_exp=zeros(N,TL);
y_log_exp=zeros(N,TL);
ye_exp=zeros(N,TL);
randi=rand(N,TL);

for t=1:TL
    y_exp(:,t)=ylook(y_matrix(:,t),t);
    y_log_exp(:,t)=ylifelog(y_matrix(:,t),t);
    ye_exp(:,t)=ye(ye_matrix(:,t),t);
    
    ymenos1=y_matrix(:,t)-1;
    ymenos1=max(ymenos1,1);
    ymas1=y_matrix(:,t)+1;
    ymas1=min(ymas1,dy);
    
    y_log_exp(randi(:,t)<0.5,t)=(0.5+randi(randi(:,t)<0.5,t)).*ylifelog(y_matrix(randi(:,t)<0.5,t),t)+...
        (0.5-randi(randi(:,t)<0.5,t)).*ylifelog(ymenos1(randi(:,t)<0.5),t);
    
    y_log_exp(randi(:,t)>0.5,t)=(randi(randi(:,t)>0.5,t)).*ylifelog(y_matrix(randi(:,t)>0.5,t),t)+...
        (1-randi(randi(:,t)>0.5,t)).*ylifelog(ymas1(randi(:,t)>0.5),t);
    
end   


% Now we can generate the samples to pass to R

id=kron(ones(TL,1),(1:N)');
age=kron((25:60)',ones(N,1));
y_export=reshape(y_exp,TL*N,1);
ye_export=reshape(ye_exp,TL*N,1);
a_export=reshape(a_matrix,TL*N,1);
c_export=reshape(c_matrix,TL*N,1);
ylog_export=reshape(y_log_exp,TL*N,1);

M=[id age y_export ye_export a_export c_export ylog_export];

if fast==0
    
if commonbet==0
name3 = strcat('Outputs/dfp_simulation_case',num2str(choose),'.csv');
elseif commonbet==1
name3 = strcat('Outputs/dfp_simulation_case',num2str(choose),'.csv');    
end
csvwrite(name3,M)

mivar=var(log(c_matrix),[],1);


%plot(1:36,mivar,1:36,lvarcons(1:36))
%legend('Simu','Agre')


end
rub=0;

end