function [y,z1,z2,z3,z4,z5,z6,z7]=fn_main3(w,q)
% A function used to solve the system of F1-7 and the labor market.
% Bing Ye, 20161224

global rhoG bet J del alp epsbar b h0 h1 %!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

etaG=(1-bet)*bet^( bet/(1-bet) )*q^( 1/(1-bet) )*w^( -bet/(1-bet) );

%% Compute V, alp and k1
% initial guess
alpmat=[zeros(J-1,1);1];
LHS=1-rhoG^J*(1-del)^J*alpmat(J);
RHS=0;
for j=0:1:(J-1)
	RHS=RHS+rhoG^(j+1)*(1-del)^( alp/(1-bet)*j );
end
RHS=RHS*etaG*alp/(1-bet);
k1=( LHS/RHS )^( 1/( alp/(1-bet)-1 ) ); 
%fn=@(k1) fn_init_k1(k1,etaG);    % old version
%k1=fsolve(fn,1);    % old version

% iterate until the convergence is reached
err=1000;
tol=1e-8;
iter=0;
iter_max=1000;
alpmatnew=zeros(J,1);
A=zeros(J,J);
B=zeros(J,1);
while err>tol & iter<iter_max
    %
    if iter>0
        alpmat=alpmatnew;
        k1=k1new;
    end
        
    % update V
    A(1,1)=1-rhoG*alpmat(1);
	for i=2:J
	    A(i,i)=1;
		A(i,1)=-rhoG*alpmat(i);
	end
	for i=1:(J-1)
	    A(i,i+1)=-rhoG*(1-alpmat(i));
	end
    for j=1:J
	    B(j,1)=rhoG*( etaG*k1^( alp/(1-bet) )*(1-del)^( alp/(1-bet)*(j-1) )...
               -alpmat(j)*( 1-(1-del)^j )*k1-w*fn_dist_cost(alpmat(j)) );  
	end
	Vaug=A\B;
    
    % update k1 and alpmat
	LHS=1;
	RHS=0;
	prod=1;
	for j=0:1:(J-1)
	    LHS=LHS-rhoG^(j+1)*prod*(1-del)^(j+1)*alpmat(j+1);
		RHS=RHS+rhoG^(j+1)*prod*(1-del)^( alp/(1-bet)*j );
        prod=prod*( 1-alpmat(j+1) );
	end
	RHS=RHS*etaG*alp/(1-bet);
	k1new=( LHS/RHS )^( 1/( alp/(1-bet)-1 ) ); 
		
    epsmat=zeros(J-1,1);
    for j=1:(J-1)
        epsmat(j)=( Vaug(1)-Vaug(j+1)-k1*( 1-(1-del)^j ) )/w;
        epsmat(j)=min(epsbar,epsmat(j));
        epsmat(j)=max(0,epsmat(j));
        alpmatnew(j)=fn_dist_prob(epsmat(j));
    end
    alpmatnew(J)=1;
	    
    %
    err=max([max(abs(alpmatnew-alpmat)),abs(k1new-k1)]);
    iter=iter+1;   
end
if err>tol
    disp('fail to converge')
    % return
end

%% Compute theta
themat=zeros(J,1);
themat(1)=1;
for j=2:J
    themat(j)=themat(j-1)*( 1-alpmat(j-1) );
end
themat=themat/sum(themat);

%% Compute the excess demand
% labor demand
zet=( bet*q/w )^( 1/(1-bet) );
labD=0;
for j=1:J
    labD=labD+themat(j)*zet*( k1*(1-del)^(j-1) )^( alp/(1-bet) )+themat(j)*fn_dist_cost(alpmat(j));
end

% labor supply
labS=( (1-b*rhoG)*w/h0 )^(1/h1);

% excess labor demand
y=labD-labS;

if nargout>1
    z1=themat;
    z2=Vaug;
    z3=k1;
    z4=alpmat;
    z5=epsmat;
    z6=labD;
    z7=iter;
end
