function [wel_region,wel_sector] = wel3(NIM,W1_2017_C1)
%NIM 表示流入矩阵（省际调入和进口）961*17；Wel_C 表示反事实福利效应 30*17
Outflow = NIM;
DTrade = zeros(900,17);
FTrade = zeros(30,17);
for i=1:1:30
  DTrade(1+(i-1)*30:i*30,:)= Outflow(1+(i-1)*31:i*31-1,:); %省际调入矩阵 900*17
end
for i=1:1:30
    FTrade(i,:)=Outflow(i*31,:);   %进口矩阵 30*17
end
MDTrade = zeros(30,17);
for i=1:1:30
    SBTrade = DTrade(1+(i-1)*30:i*30,:);
    MDTrade(i,:)=sum(SBTrade);
end
% 各省份各部门总消费
Total_CS = MDTrade+FTrade;
% 调入
Trade_pro = sum(Total_CS,2);           % 各省消费加总 **** 30*1
Trade_na = sum(Trade_pro);             % 全国消费加总 1*1
Share_pro_sec = Total_CS./Trade_pro;   % 各省各部门消费份额 30*17
Share_pro = Trade_pro./Trade_na;       % 各省消费份额 30*1
Trade_sec = sum(Total_CS);             % 各部门消费总额 ****  1*17
Share_sec_pro = Total_CS./Trade_sec;   % 各部门各省消费份额 30*17
%% 各地区流入和消费份额
[Coa_tra,Inl_tra]=Reg_Di(Trade_pro);
Share_coa = Coa_tra./sum(Coa_tra);  % 沿海地区消费份额
Share_inl = Inl_tra./sum(Inl_tra);  % 内陆地区消费份额
%% 各部门流入和消费份额
[Lab_tra,Cap_tra]=Sec_Di(Trade_sec);
Share_lab = Lab_tra./sum(Lab_tra);  % 劳动部门消费份额
Share_cap = Cap_tra./sum(Cap_tra);  % 资本部门消费份额
%% 福利加权
% 地区层面
W1_2017_C1(W1_2017_C1==0)=1;
wel_pro_0 = W1_2017_C1.^Share_pro_sec;
wel_pro = prod(wel_pro_0,2);        % 各省福利加权 ****
wel_na_0 = wel_pro.^Share_pro;     
wel_na = prod(wel_na_0);          % 全国福利加权 ****
[Coa_wel,Inl_wel]=Reg_Di(wel_pro);
wel_coa_0 = Coa_wel.^Share_coa;
wel_coa = prod(wel_coa_0);        % 沿海地区福利加权 ***
wel_inl_0 =Inl_wel.^Share_inl;
wel_inl = prod(wel_inl_0);        % 内陆地区福利加权 ***
wel_region = [wel_na;wel_coa;wel_inl];
% 部门层面
wel_sec_0 =W1_2017_C1.^Share_sec_pro;
wel_sec = prod(wel_sec_0);        % 各部门福利加权 ****
[Lab_wel,Cap_wel]=Sec_Di(wel_sec);
wel_lab_0 = Lab_wel.^Share_lab;
wel_lab = prod(wel_lab_0);        % 劳动部门福利加权 ***
wel_cap_0 = Cap_wel.^Share_cap;
wel_cap = prod(wel_cap_0);        % 资本部门福利加权 ***
Indus_trade =[Lab_tra;Cap_tra];
Share_indus = Indus_trade./sum(Indus_trade);
Indus_wel = [Lab_wel;Cap_wel];
wel_indus_0 = Indus_wel.^Share_indus;
wel_indus = prod(wel_indus_0);   % 工业部门福利加权
wel_sector = [wel_indus;wel_lab;wel_cap];
end