
function [wagehat, Youtput_test]= SolvingWagesNF_DC(PiFp,Pip,tauXp,tauFp,gammas,alphas,VAn,Surplusnp,J,N,S,wagehat,vfactor)
    %CtyConsp = max(wagehat.*VAn - Surplusnp,1e-18);
    CtyConsp = wagehat.*VAn - Surplusnp;
    F1=1-sum(PiFp./tauFp.*(tauFp-1).*repmat(alphas,N,1));
    F1_test=1-sum(PiFp./tauFp.*(tauFp-1).*repmat(alphas,N,1));
    CtyConsp_F_test=CtyConsp./F1_test';
    FGS = sum(PiFp.*repmat(alphas,N,1).*repmat(CtyConsp',J,1),2);
    FGS_test = sum(PiFp./tauFp.*repmat(alphas,N,1).*repmat(CtyConsp_F_test',J,1),2);
    IGS1 = Pip.*repmat(gammas,N,1);
    IGS21 = sum(Pip./tauXp.*(tauXp-1).*repmat(gammas,N,1));
    IGS22 = reshape(IGS21,56,44);
    IGS23 = PiFp.*repmat(alphas,N,1).*repmat(IGS22,N,1);
    IGS2 = repmat(IGS23',56,1)';
    Youtput_test = inv(eye(J) - IGS1-IGS2)*FGS_test;
    Youtput = inv(eye(J) - IGS1)*FGS;
    
    %Youtput(Youtput <= 0) = 0;

    %Imports
    LHS = zeros(N,1);
    LHS_test=zeros(N,1);
    %Exports
    RHS = zeros(N,1);
    RHS_test=zeros(N,1);
    IGS3=sum(IGS22.*reshape(Youtput,56,44));
    for j = 1:N
        for s = 1:S
            for i = 1:N
                iiss = (i-1)*S+s;
                jjss = (j-1)*S+s;
                for r = 1:S
                    jjrr = (j-1)*S+r;
                    iirr = (i-1)*S+r;
                    LHS(j) = LHS(j) + Pip(iiss,jjrr)*gammas(s,jjrr)*Youtput(jjrr);
                    RHS(j) = RHS(j) + Pip(jjss,iirr)*gammas(s,iirr)*Youtput(iirr);
                    LHS_test(j) = LHS_test(j) + Pip(iiss,jjrr)/tauXp(iiss,jjrr)*gammas(s,jjrr)*Youtput_test(jjrr);
                    RHS_test(j) = RHS_test(j) + Pip(jjss,iirr)/tauXp(jjss,iirr)*gammas(s,iirr)*Youtput_test(iirr);
                end
                LHS(j) = LHS(j) + PiFp(iiss,j)*alphas(s,j)*CtyConsp(j);
                RHS(j) = RHS(j) + PiFp(jjss,i)*alphas(s,i)*CtyConsp(i);
                LHS_test(j) = LHS_test(j) + PiFp(iiss,j)/tauFp(iiss,j)*alphas(s,j)*CtyConsp_F_test(j)+PiFp(iiss,j)*alphas(s,j)*IGS3(j)/F1_test(j);
               RHS_test(j) = RHS_test(j) + PiFp(jjss,i)/tauFp(jjss,i)*alphas(s,i)*CtyConsp_F_test(i)+PiFp(jjss,i)*alphas(s,i)*IGS3(i)/F1_test(i);
            end
        end
    end
    %Exports = Imports + Surplus
    ExcessSurplus = RHS - LHS - Surplusnp;
    ExcessSurplus_test = RHS_test - LHS_test - Surplusnp;
    
    %Normalized Excess Surplus
    %NES = ExcessSurplus./CtyConsp;
    NES = ExcessSurplus./VAn;
    NES_test = ExcessSurplus_test./VAn;
    %ExcessSurplus means wages should go up to kill exports (and drop ES)
    %delta = sign(NES).*min(abs(vfactor*NES),50);
    delta = sign(NES).*abs(vfactor*NES);
    delta_test = sign(NES_test).*abs(vfactor*NES_test);
    %wagehat_cankao =wagehat.*(1+delta./wagehat);
    wagehat =wagehat.*(1+delta_test./wagehat);
    YYp = Youtput;
    YYp_test = Youtput_test;
end 

