function [PiHat, PiFHat, PrsjHat, PrFjHat, cost]= PriceHatNF_DC0424(wagehat,tauXhat,tauFhat,VAcoeff,gammas,Pi,PiF,Theta1,Theta2,J,N,S)

    wagehat0 = wagehat;
    PrsjHat0 = ones(S,J);
    PrFjHat0 = ones(S,N);
    tol = 1E-3;
    maxit = 100;
    pfmax = 1;
    iterph  = 1;    

    while (iterph <= maxit) && (pfmax > tol) %This gets me the prices
        cost = ones(S,N);
        for s = 1:1:S,
            for j = 1:1:N,
                cost(s,j) = wagehat0(j)^VAcoeff(s,j);
                for r = 1:1:S,
                    cost(s,j) = cost(s,j)*PrsjHat0(r,(j-1)*S+s)^gammas(r,(j-1)*S+s);
                end
            end
        end
        %This is equation 12 and 11
        PrsijHat = zeros(J,J);
        for i = 1:1:N
            for r = 1:1:S
                for j = 1:1:N
                    for s = 1:1:S
                        iirr = (i-1)*S+r;
                        jjss = (j-1)*S+s;
                        PrsijHat(iirr,jjss) = (tauXhat(iirr,jjss)*cost(r,i))^(-Theta1(r));
                    end
                    % PrFijHat(iirr,j) = PiF(iirr,j)*(tauFhat(iirr,j)*cost(r,i))^(-Theta(r));
                    % PrFijHat(iirr,j) = (tauFhat(iirr,j)*cost(r,i))^(-Theta(r));
                end
            end
        end
        PrsjHat = zeros(S,J);
        for r=1:1:S
            for j = 1:1:N
                for s = 1:1:S
                    jjss = (j-1)*S+s;
                    PrsjHat(r,jjss) = (sum(Pi(r:S:(N-1)*S+r,jjss).*PrsijHat(r:S:(N-1)*S+r,jjss)))^(-1/Theta1(r));
                end
                % PrFjHat(r,j) = (sum(PiF(r:S:J,j).*PrFijHat(r:S:J,j)))^(-1/Theta(r));
            end
        end
        
        pfdev = [abs(PrsjHat - PrsjHat0)]; % Checking tolerance of convergence
        PrsjHat0 = PrsjHat;
        pfmax = max(max(pfdev));
        iterph = iterph + 1;
    end
    iterph
    
    PrFijHat = zeros(J,N);
    for i = 1:1:N
        for r = 1:1:S
            for j = 1:1:N
                iirr = (i-1)*S+r;
                PrFijHat(iirr,j) = (tauFhat(iirr,j)*cost(r,i))^(-Theta2(r));
            end
        end
    end
    PrFjHat = zeros(S,N);      
    for r=1:1:S
        for j = 1:1:N
            PrFjHat(r,j) = (sum(PiF(r:S:(N-1)*S+r,j).*PrFijHat(r:S:(N-1)*S+r,j)))^(-1/Theta2(r));
        end
    end
    
    Tall = repmat(Theta1,N,J);
    TF = repmat(Theta2,N,N);
    %Verified that sum(Pi*Pip) == repmat(S,1,J)
    PiHat = PrsijHat./(repmat(PrsjHat,N,1).^(-Tall));
    PiFHat = PrFijHat./(repmat(PrFjHat,N,1).^(-TF));

end
