function W_hat = wel1(X,tau_hat,epsilon,ls)
% N 表示双边贸易矩阵;tau_hat 表示贸易成本变动; epsilon表示弹性;ls劳动投入份额
[N,~]=size(X);
Y = sum(X,2); 
w_hat_0 = ones(N,1); 
    diff = 1; 
    tol = 0.0001; 
    update = 0.1;
    while diff>tol 
        P_hat = (sum((X'./(Y*ones(1,N))).*(tau_hat'.^-epsilon).*(ones(N,1)*(w_hat_0.^-epsilon)'),2)).^(-1/epsilon); 
        w_hat_1 = sum((X./(Y*ones(1,N))).*(tau_hat.^-epsilon).*((w_hat_0.^-epsilon)*ones(1,N)).*(ones(N,1)*((P_hat.^epsilon).*w_hat_0)'),2); 
        w_hat_1 = w_hat_1./w_hat_1(1); 
        diff = sum(((w_hat_0 - w_hat_1).^2).^0.5); 
        w_hat_0 = (1-update)*w_hat_0 + update*w_hat_1;
    end
%% 福利变化
w_hat = w_hat_0;
w_hat(isnan(w_hat))=1;
w_hat(isinf(w_hat))=1;
Pi_hat = ones(N,1).*(w_hat.^-epsilon)'.*(tau_hat'.^-epsilon).*(ones(N,1)*(P_hat.^epsilon)');
Pi_jj_hat = zeros(N,1);
for i=1:1:N
    for j=1:1:N
    if i==j
        Pi_jj_hat(i,1)=Pi_hat(i,j);
    end
    end
end
nels = (epsilon*ls').^(-1);
W_hat = Pi_jj_hat.^nels; 
end