%% 20200901修改
%利用Estimation_baseline_0824.m的均衡数据（即校准过的数据）作为模拟的初步输入变量（与cp处理方式一致）
%20220302第一次修改，修改所有的贸易弹性为4.其中修改的地方为line28和line29（新增两行），其他不变，得到的结果保存为20220302.mat
%20220302第二次修改，修改所有的贸易弹性为Tombe（2015）估计θ=4.1（农业）和4.6（非农业）的情况。.其中修改的地方为line31和line32、33（新增三行），其他不变，得到的结果保存为20220302.mat
%20220302第三次修改，修改所有的贸易弹性为自己估的2014年的弹性的情况。.其中修改的地方为line34（新增两行），其他不变，得到的结果保存为20220302_2014.mat
close all
clear all
clc
userpath('');
savepath
vfactor  = -0.1;%cp设定为-0.2
tol      = 1E-7;%cp设定为1E-7
maxit    = 500;%cp设定为1E10
wfmax = 1E7;

N = 44;
S = 56;
J = N*S;

%%%%%%%%%%%%%%%%%
% Modify the following to specify the set of theta's 
% Either constant theta or Caliendo-Parro theta's
% cons = 5;
% Constheta = cons*ones(S,1);
% % CPtheta = importdata('CPtheta.txt');
% Theta1 = Constheta;
%%修改弹性20200824，再修改贸易弹性为所有弹性均为4，依据是Simonovska和Waugh（2014）利用跨国价格数据估计θ≈4。2022.3.2日修改
T = xlsread('T.csv','B2:B24');
%T_xiugai = ones(23,1)*4;%新增，20220302，结果保存为20220302.mat
%T = T_xiugai;%新增，20220302
%T_xiugai = ones(23,1)*4.6;%新增，20220302，结果保存为20220302_4.mat
%T_xiugai(1:3,1) = 4.1;%修改弹性20200824，再修改农业贸易弹性为所有弹性均为4，依据是Tombe（2015）估计θ=4.1（农业）和4.6（非农业）的情况。。2022.3.2日修改,20220302，结果保存为20220302_4.6.mat
%T = T_xiugai;%新增，20220302
%T = xlsread('T.csv','A2:A24');
%T=[1./T; ones(33,1)*(1/8.22)];
% Theta = [T; ones(33,1)*8.22];
%%修改弹性20200424
% T1 = csvread('T1.csv');%T=[1./T; ones(33,1)*(1/8.22)];
% T2 =  csvread('T2.csv');
Theta1 = [T; ones(33,1)*8.22];
Theta2 = [T; ones(33,1)*8.22];
%%%%%%%%%%%%%%%%%%

%load the equilibrim data to counterfacual;
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%两种方式，一种是利用原始数据模拟变化的均衡，一种是利用均衡后的数据模拟变化的均衡，先实验第一种。
k=1;%2000年
head = 'WIO2016T';
WIOT = [head,num2str(k),'.xlsb'];
XXread = xlsread(WIOT,'E7:CPX2470');
%M = mat2cell(Ma,[56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56],[56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56]);
Ya =  xlsread(WIOT,'CYK7:CYK2470');
%Y = mat2cell(Ya,[56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56],1);
FFread = xlsread(WIOT,'CPY7:CYJ2470');
%fb = mat2cell(fa,[56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56 56],[5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5]);
Mc = sum(XXread); %每一个国家的每一个行业的中间投入品加总，即Mj（s），得到一个1*1435的向量
VVread = Ya'-Mc;
%第二种方式
% load WIOT_14_nfbaseline0626.mat
% clear FFdip XXip%为了不和新得到的变量重复，同时保留位置指标GVC_cf变量，与变化的结果作对比。
%%%%%%%%%%%%%%%%%
% Modify the following to get the tauhat's and tauFhat's of interest 
% (some examples given)
tariff_int2000 = xlsread('inter_tariff_recount20200814.csv','C2:AF1681');%注意C2和s953不能为NaN,不然就会丢失数据，注意列为关税征收国.确定最后一个国家是欧盟还是世界其他国家啊，核验20200824
tariff_fin2000 = xlsread('fin_tariff_recount20200814.csv','C2:AF1681');%注意C2和s953不能为NaN,不然就会丢失数据，注意列为关税征收国
%读入2007年关税数据
tariff_int2007 = xlsread('inter_tariff_recount_2007_20200827.csv','C2:AF1681');%注意C2和s953不能为NaN,不然就会丢失数据，注意列为关税征收国.确定最后一个国家是欧盟还是世界其他国家啊，核验20200824
tariff_fin2007 = xlsread('fin_tariff_recount_2007_20200827.csv','C2:AF1681');%注意C2和s953不能为NaN,不然就会丢失数据，注意列为关税征收国

%最终品第4个行业无关税，用中间品的代替。
%中间品第3个行业无关税，用最终品的代替。
for j =1:30    
       tariff_fin2000(4+(j-1)*S,:)=tariff_int2000(4+(j-1)*S,:);
       tariff_int2000(3+(j-1)*S,:)=tariff_fin2000(3+(j-1)*S,:); 
       tariff_fin2007(4+(j-1)*S,:)=tariff_int2007(4+(j-1)*S,:);
       tariff_int2007(3+(j-1)*S,:)=tariff_fin2007(3+(j-1)*S,:); 
end
%修改为2007年中国变化关税，即taup

%扩展为44国56行业
%C1=[3,4,9,10,11,12,13,14,15,16,17,18,19,20,23,24,27,28,29,31,32,34,35,36,38,39,40];
%C2=[1,2,5,6,7,8,21,22,25,26,30,33,37,41,42,43,44];
C1=[2,3,11,12,13,15,16,17,18,23,24,28,32,35,40];
C2=[1,4,5,6,7,8,9,10,14,19,20,21,22,25,26,27,29,30,31,33,34,36,37,38,39,41,42,43,44];%注意，欧盟为30.44正好排名第29
for c1 =1:29
      c11= C2(1,c1);
      tariff_fin2000_1(:,c11)=tariff_fin2000(:,c1);
      tariff_int2000_1(:,1+(c11-1)*S:c11*S)=repmat(tariff_int2000(:,c1),1,56);
      tariff_fin2007_1(:,c11)=tariff_fin2007(:,c1);
      tariff_int2007_1(:,1+(c11-1)*S:c11*S)=repmat(tariff_int2007(:,c1),1,56);
end
for c1 =1:15
      c11= C1(1,c1);
      tariff_fin2000_1(:,c11)=tariff_fin2000(:,30);
      tariff_int2000_1(:,1+(c11-1)*S:c11*S)=repmat(tariff_int2000(:,30),1,56);
      tariff_fin2007_1(:,c11)=tariff_fin2007(:,30);
      tariff_int2007_1(:,1+(c11-1)*S:c11*S)=repmat(tariff_int2007(:,30),1,56);
end

for c1 =1:29
      c11= C2(1,c1);
      tariff_fin2000_2(1+(c11-1)*S:c11*S,:)=tariff_fin2000_1(1+(c1-1)*S:c1*S,:);
      tariff_int2000_2(1+(c11-1)*S:c11*S,:)=tariff_int2000_1(1+(c1-1)*S:c1*S,:);
      tariff_fin2007_2(1+(c11-1)*S:c11*S,:)=tariff_fin2007_1(1+(c1-1)*S:c1*S,:);
      tariff_int2007_2(1+(c11-1)*S:c11*S,:)=tariff_int2007_1(1+(c1-1)*S:c1*S,:);
end
for c1 =1:15
      c11= C1(1,c1);
       tariff_fin2000_2(1+(c11-1)*S:c11*S,:)=tariff_fin2000_1(29*S+1:30*S,:);
       tariff_int2000_2(1+(c11-1)*S:c11*S,:)=tariff_int2000_1(29*S+1:30*S,:);
       tariff_fin2007_2(1+(c11-1)*S:c11*S,:)=tariff_fin2007_1(29*S+1:30*S,:);
       tariff_int2007_2(1+(c11-1)*S:c11*S,:)=tariff_int2007_1(29*S+1:30*S,:);
      
end
%还有一个问题需要解决，就是欧盟国家对欧盟国家的关税应为0，检验是否有，已检验，没问题。
inttariff=tariff_int2000_2;
fintariff_L = tariff_fin2000_2;
inttariff(isnan(inttariff))=0;
fintariff_L(isnan(fintariff_L))=0;

% 构建    tauXp和 tauFp 
% inttariffp=tariff_int2007_2;
% fintariffp_L = tariff_fin2007_2;
inttariffp=inttariff;
inttariffp(:,1+7*56:8*56) = tariff_int2007_2(:,1+7*56:8*56);
fintariff_Lp=fintariff_L;
%fintariff_Lp(:,8)=tariff_fin2007_2(:,8);
inttariffp(isnan(inttariffp))=0;
fintariff_Lp(isnan(fintariff_Lp))=0;
%构建    tauhat和 tauFhat     

tauX=1+inttariff/100;                           % actual tariff vector 中间品和最终品关税  
tauF=1+fintariff_L/100;                          % counterfactual tariff vector
tauXp=1+inttariffp/100;
tauFp=1+fintariff_Lp/100;
tauXhat=tauXp./tauX;
tauFhat=tauFp./tauF;
%把数字变小，wiot中的单位是million，借鉴cp的，要缩小两个100000，因此要再缩小10000.20200515新加程序

 %将关税纳入到数据校准中，修改WIOTcleaning_f_DC程序
[XXi, FFdi,YY, CtyCons, gammas, VAcoeff, alphas, Surplusn, VAn, Pi, PiF,XXi_gs,FFdi_gs] = WIOTcleaning_f0515_DC(FFread, XXread, VVread, tauX, tauF, J, S, N);



%%%%%%%%%%%%%%%%%
% Final surplus value - can be adjusted.
% Should take this to be either 0 or the final surplus value in 2011.
Surplusnp = Surplusn;
%%%%%%%%%%%%%%%%%
%计算Ij,解决思路见“记录”的20200515的第3条。20200515新加程序
Ij = (VAn+XXi_gs+FFdi_gs-Surplusn); %VAn和Surplusn已知。第一种方法,Ijp的计算在Welfarelinal.m中。
%Initialize wages
wagehat=ones(N,1); 
iter = 1;
Price_Index0 = ones(N,1);
for s = 1:S
    Price_Index0 = Price_Index0.*(((ones(1,N)./alphas(s,:)).^alphas(s,:))');
end


while (iter <= maxit) && (wfmax > tol)
     [PiHat, PiFHat, PrsjHat, PrFjHat, cost] =PriceHatNF_DC0424(wagehat,tauXhat,tauFhat,VAcoeff,gammas,Pi,PiF,Theta1,Theta2,J,N,S);

    PiHat(isinf(PiHat))=0;
    PiFHat(isinf(PiFHat))=0;
    Pip=PiHat.*Pi;
    PiFp = PiFHat.*PiF;

    %Labor Market Clearing and solving for wages
    
    wagehatprev = wagehat;
    [wagehat, Youtput_test] = SolvingWagesNF_DC(PiFp,Pip,tauXp,tauFp,gammas,alphas,VAn,Surplusnp,J,N,S,wagehat,vfactor);      
    wfmax = max(abs(wagehat-wagehatprev))
    iter = iter+1
     
    if min(wagehat) < 0,
        wagehat = wagehatprev;
    end
end

Price_Index_hat = ones(N,1);
for s = 1:S
    Price_Index_hat = Price_Index_hat.*(((PrFjHat(s,:)./alphas(s,:)).^alphas(s,:))');
end

realwagehatp = wagehat.*Price_Index0./Price_Index_hat;


%真实工资的变化 0515新增
P=prod(PrFjHat.^(alphas)); %price index
W_all=wagehat./P';      %real wages

%福利的变化  0515新增
 [Welfare ToTF VoTF ToTEjsF ToTMjsF VoTMjsF  ToTX VoTX ToTEjsX ToTMjsX VoTMjsX votEijsF votEijsX EjisF_dc EjisX_dc Ijp] = Welfarelineal0515(Ij,PiF,Pi,alphas,gammas,tauF,tauX,cost,Surplusnp,VAn,wagehat,tauFp,PiFp,tauXp,Pip,YY,Youtput_test,J,S,N);
%计算中国的福利 20201029新增
aus = 1; aut = 2; bel = 3; bgr = 4; bra = 5; can = 6; che = 7;chn = 8; cyp  = 9; cze = 10; deu = 11;
dnk =12; esp = 13; est =14;fin = 15;fra = 16;gbr = 17; grc = 18;hrv = 19; hun = 20; idn =21; ind = 22;
irl = 23;ita = 24; jpn =25; kor = 26;ltu =27; lux = 28; lva =29; mex = 30; mlt =31;nld = 32;nor = 33;
pol  =34; prt =35; rou = 36;rus =37;svk = 38; svn = 39;swe =40; tur =41;twn = 42;usa = 43; 
row = 44;

Welfare_chn=[Welfare(8,1);ToTX(8,1);ToTF(8,1);VoTX(8,1);VoTF(8,1);W_all(8,1)-1];
Welfare_ame=[Welfare(43,1);ToTX(43,1);ToTF(43,1);VoTX(43,1);VoTF(43,1);W_all(43,1)-1];
share_Ijp=Ijp/sum(Ijp);
Welfare_world=[sum(Welfare.*share_Ijp);sum(ToTX.*share_Ijp);sum(ToTF.*share_Ijp);sum(VoTX.*share_Ijp);sum(VoTF.*share_Ijp);sum((W_all-1).*share_Ijp);];%通过收入占比加权得到
Welfare_world_c=[Welfare(8,1)*share_Ijp(8,1);ToTX(8,1)*share_Ijp(8,1);ToTF(8,1)*share_Ijp(8,1);VoTX(8,1)*share_Ijp(8,1);VoTF(8,1)*share_Ijp(8,1);(W_all(8,1)-1)*share_Ijp(8,1)];%通过收入占比加权得到

Welfare_world_chn=[sum(Welfare.*share_Ijp);sum(ToTX.*share_Ijp);sum(ToTF.*share_Ijp);sum(VoTX.*share_Ijp);sum(VoTF.*share_Ijp);sum((W_all-1).*share_Ijp);]-Welfare_world_c;%去掉中国所占比例
Welfare_jpn=[Welfare(25,1);ToTX(25,1);ToTF(25,1);VoTX(25,1);VoTF(25,1);W_all(25,1)-1];
Welfare_kor=[Welfare(26,1);ToTX(26,1);ToTF(26,1);VoTX(26,1);VoTF(26,1);W_all(26,1)-1];
Welfare_mex=[Welfare(30,1);ToTX(30,1);ToTF(30,1);VoTX(30,1);VoTF(30,1);W_all(30,1)-1];
Welfare_ind=[Welfarnde(22,1);ToTX(22,1);ToTF(22,1);VoTX(22,1);VoTF(22,1);W_all(22,1)-1];
%计算欧盟
eu= [aut bel bgr cyp cze deu dnk esp est fin fra gbr grc hrv hun irl ita ltu lux lva mlt nld pol prt rou svk svn swe]; %8为中国，43为美国。
share_eu_Ijp=Ijp(eu)/sum(Ijp(eu));
share_eu_Ijp1=Ijp(eu)/sum(Ijp);
Welfare_eu=[sum(Welfare(eu).*share_eu_Ijp);sum(ToTX(eu).*share_eu_Ijp);sum(ToTF(eu).*share_eu_Ijp);sum(VoTX(eu).*share_eu_Ijp);sum(VoTF(eu).*share_eu_Ijp);sum((W_all(eu)-1).*share_eu_Ijp)];
Welfare_eu1=[sum(Welfare(eu).*share_eu_Ijp1);sum(ToTX(eu).*share_eu_Ijp1);sum(ToTF(eu).*share_eu_Ijp1);sum(VoTX(eu).*share_eu_Ijp1);sum(VoTF(eu).*share_eu_Ijp1);sum((W_all(eu)-1).*share_eu_Ijp1)];
%表3
table3=[Welfare_world';Welfare_chn';Welfare_world_chn';Welfare_ame';Welfare_jpn';Welfare_kor';Welfare_mex';Welfare_ind';Welfare_eu'];

 %计算中美以及其他国家双边福利的变化，1030新增。计算VoTXca,VoTFca,ToTXca,ToTFca
ca= [8 22]; %8为中国，43为美国。 分别录入22（印度）、25日本、26韩国、30墨西哥 43美国
for n = ca
    for   s =  1:1:S
    VoTXca1(s,n) = nansum(votEijsX((ca-1)*S+s , n));%计算中国和美国分别来自于中国+美国的双边。其中ca的用法很好，每次都加在一起。
    VoTFca1(s,n) = nansum(votEijsF((ca-1)*S+s , n));%计算中国和美国分别来自于中国+美国的双边。其中ca的用法很好，每次都加在一起。
    
    ToTEjsFca1(s,n) = sum(EjisF_dc(s + S*(n-1),ca));
    ToTMjsFca1(s,n) = sum(EjisF_dc((ca-1)*S+s,n));
      
    ToTEjsXca1(s,n) = sum(EjisX_dc(s + S*(n-1),ca));
    ToTMjsXca1(s,n) = sum(EjisX_dc((ca-1)*S+s,n));
    end
end

VoTXca2 = sum(VoTXca1);
VoTFca2 = sum(VoTFca1);
VoTXca = VoTXca2(ca)'./Ij(ca);
VoTFca = VoTFca2(ca)'./Ij(ca);
VoTXcaROW= VoTX(ca)-VoTXca;
VoTFcaROW= VoTF(ca)-VoTFca;

ToTXca = (sum(ToTEjsXca1(:,ca))' - sum(ToTMjsXca1(:,ca))')./Ij(ca);
ToTFca = (sum(ToTEjsFca1(:,ca))' - sum(ToTMjsFca1(:,ca))')./Ij(ca);
ToTXcaROW = ToTX(ca)-ToTXca;
ToTFcaROW = ToTF(ca)-ToTFca;
welfare_chname=[ToTXca(1);ToTFca(1);VoTXca(1);VoTFca(1)];
%计算中欧盟福利的变化，1030新增。计算VoTXca,VoTFca,ToTXca,ToTFca
ce= [aut bel bgr chn cyp cze deu dnk esp est fin fra gbr grc hrv hun irl ita ltu lux lva mlt nld pol prt rou svk svn swe]; %8为中国，43为美国。
for n = ce
    for   s =  1:1:S
    VoTXce1(s,n) = nansum(votEijsX((ce-1)*S+s , n));%计算中国和美国分别来自于中国+美国的双边。其中ca的用法很好，每次都加在一起。
    VoTFce1(s,n) = nansum(votEijsF((ce-1)*S+s , n));%计算中国和美国分别来自于中国+美国的双边。其中ca的用法很好，每次都加在一起。
    
    ToTEjsFce1(s,n) = sum(EjisF_dc(s + S*(n-1),ce));
    ToTMjsFce1(s,n) = sum(EjisF_dc((ce-1)*S+s,n));
      
    ToTEjsXce1(s,n) = sum(EjisX_dc(s + S*(n-1),ce));
    ToTMjsXce1(s,n) = sum(EjisX_dc((ce-1)*S+s,n));
    end
end

VoTXce2 = sum(VoTXce1);
VoTFce2 = sum(VoTFce1);
VoTXce = VoTXce2(ce)'./Ij(ce);
VoTFce = VoTFce2(ce)'./Ij(ce);
VoTXceROW= VoTX(ce)-VoTXce;
VoTFceROW= VoTF(ce)-VoTFce;

ToTXce = (sum(ToTEjsXce1(:,ce))' - sum(ToTMjsXce1(:,ce))')./Ij(ce);
ToTFce = (sum(ToTEjsFce1(:,ce))' - sum(ToTMjsFce1(:,ce))')./Ij(ce);
ToTXceROW = ToTX(ce)-ToTXce;
ToTFceROW = ToTF(ce)-ToTFce;
welfare_chneu=[ToTXce(4);ToTFce(4);VoTXce(4);VoTFce(4)];
%计算行业对福利效应的贡献
ToTX_sharewelfar=(ToTEjsX(1:23,ca)-ToTMjsX(1:23,ca))./nansum(ToTEjsX(1:23,ca)-ToTMjsX(1:23,ca));
ToTF_sharewelfar=(ToTEjsF(1:23,ca)-ToTMjsF(1:23,ca))./nansum(ToTEjsF(1:23,ca)-ToTMjsF(1:23,ca));
VoTX_sharewelfar=VoTMjsX(1:23,ca)./nansum(VoTMjsX(1:23,ca));
VoTF_sharewelfar=VoTMjsF(1:23,ca)./nansum(VoTMjsF(1:23,ca));
table5=[ToTX_sharewelfar(:,1),ToTF_sharewelfar(:,1),VoTX_sharewelfar(:,1),VoTF_sharewelfar(:,1)];


% 贸易的变化，0604新增
%compute bilateral trade flows in the counterfactutal
%%注意，XXip与投入产出表的结构是一样的，行向量表示产出，即来源。列向量表示目的，即去处，因此其中一个元素表示为来源于行向量国家去往列向量国家的出口值。
XXip = Pip./tauXp.*repmat(gammas,N,1).*repmat(Youtput_test',J,1);
FFdip = PiFp./tauFp.*repmat(alphas,N,1).*repmat(Ijp',J,1);
%compute the export flows matrices in the base year and counterfactual eq.
FFdi_tiaozheng= zeros(J,N);%调整是为了去掉自己国家的内贸
FFdip_tiaozheng= zeros(J,N);
XXi_tiaozheng= zeros(J,J);
XXip_tiaozheng= zeros(J,J);
for j=1:N
    FFdi_tiaozheng(1+(j-1)*S:j*S,j)=FFdi(1+(j-1)*S:j*S,j);
    FFdip_tiaozheng(1+(j-1)*S:j*S,j)=FFdip(1+(j-1)*S:j*S,j);
    FFdi_exp(1+(j-1)*S:j*S,:)=FFdi(1+(j-1)*S:j*S,:)-FFdi_tiaozheng(1+(j-1)*S:j*S,:);
    FFdip_exp(1+(j-1)*S:j*S,:)=FFdip(1+(j-1)*S:j*S,:)-FFdip_tiaozheng(1+(j-1)*S:j*S,:);
    XXi_tiaozheng(1+(j-1)*S:j*S,1+(j-1)*S:j*S)=XXi(1+(j-1)*S:j*S,1+(j-1)*S:j*S);
    XXip_tiaozheng(1+(j-1)*S:j*S,1+(j-1)*S:j*S)=XXip(1+(j-1)*S:j*S,1+(j-1)*S:j*S);
    XXi_exp(1+(j-1)*S:j*S,:)=XXi(1+(j-1)*S:j*S,:)-XXi_tiaozheng(1+(j-1)*S:j*S,:);
    XXip_exp(1+(j-1)*S:j*S,:)=XXip(1+(j-1)*S:j*S,:)-XXip_tiaozheng(1+(j-1)*S:j*S,:);
end
%clear FFdi_tiaozheng FFdip_tiaozheng XXi_tiaozheng XXip_tiaozheng
%计算出口到全世界的变化
XXip_Ejn=sum(XXip_exp')';
FFdip_Ejn=sum(FFdip_exp')';
XXi_Ejn=sum(XXi_exp')';
FFdi_Ejn=sum(FFdi_exp')';
%出口到全世界扣除中国和美国，这个计算结果对中美是正确的。
gXXip_row_exp =XXip_Ejn-sum(XXip_exp(:,7*56+1:8*56)')'-sum(XXip_exp(:,42*S+1:43*S)')';
gXXi_row_exp =XXi_Ejn-sum(XXi_exp(:,7*56+1:8*56)')'-sum(XXi_exp(:,42*S+1:43*S)')';%只有中国和美国的可以用
gFFdip_row_exp =FFdip_Ejn-FFdip_exp(:,8)-FFdip_exp(:,43);
gFFdi_row_exp = FFdi_Ejn-FFdi_exp(:,8)-FFdi_exp(:,43);

gXXi_row_exp_chn =sum(gXXip_row_exp(7*S+1:8*S,1))./sum(gXXi_row_exp(7*S+1:8*S,1))-1;
gXXi_row_exp_ame =sum(gXXip_row_exp(42*S+1:43*S,1))./sum(gXXi_row_exp(42*S+1:43*S,1))-1;
gFFdi_row_exp_chn =sum(gFFdip_row_exp(7*S+1:8*S,1))./sum(gFFdi_row_exp(7*S+1:8*S,1))-1;
gFFdi_row_exp_ame =sum(gFFdip_row_exp(42*S+1:43*S,1))./sum(gFFdi_row_exp(42*S+1:43*S,1))-1;
gzong_row_exp_chn =(sum(gFFdip_row_exp(7*S+1:8*S,1))+sum(gXXip_row_exp(7*S+1:8*S,1)))./(sum(gFFdi_row_exp(7*S+1:8*S,1))+sum(gXXi_row_exp(7*S+1:8*S,1)))-1;
gzong_row_exp_ame =(sum(gFFdip_row_exp(42*S+1:43*S,1))+sum(gXXip_row_exp(42*S+1:43*S,1)))./(sum(gFFdi_row_exp(42*S+1:43*S,1))+sum(gXXi_row_exp(42*S+1:43*S,1)))-1;


Ejnp1=FFdip_Ejn+XXip_Ejn;
Ejnp = reshape(Ejnp1,S,N);
shareEjnp=Ejnp(1:23,:)./(ones(23,1)*sum(Ejnp(1:23,:))); %export shares by sectors in the counterfactual eq.
shareEjnp_all=shareEjnp;
HH_exp_allp=sum(shareEjnp.*shareEjnp);
HH_exp_allp=((HH_exp_allp-1/23)./(1-1/23))';  %Normalized Herfindahl index of export concentration in the counterfactual equilibrium

%compute total exports by country and sector in the base year
Ejn1=FFdi_Ejn+XXi_Ejn;
Ejn = reshape(Ejn1,S,N);
shareEjn=Ejn(1:23,:)./(ones(23,1)*sum(Ejn(1:23,:))); %export shares by sectors in the base year
shareEjn_all=shareEjn;
HH_exp_all=sum(shareEjn.*shareEjn);
HH_exp_all=((HH_exp_all-1/23)./(1-1/23))'; %Normalized Herfindahl index of export concentration in the base year

share_Ejn_bianhua=[shareEjn_all(:,8),shareEjnp_all(:,8)];
HH_exp_bianhua=[HH_exp_all(8),HH_exp_allp(8)];
table2=[share_Ejn_bianhua;HH_exp_bianhua];

%compute import growth
for i=1:N
    FFdi_imp(i,:)=sum(FFdi_exp(1+(i-1)*S:i*S,:));
    FFdip_imp(i,:)=sum(FFdip_exp(1+(i-1)*S:i*S,:));
    for j= 1:N
    XXi_imp(i,j)=sum(sum(XXi_exp(1+(i-1)*S:i*S,1+(j-1)*S:j*S),2));
    XXip_imp(i,j)=sum(sum(XXip_exp(1+(i-1)*S:i*S,1+(j-1)*S:j*S),2));
    end
end

gXXi_imp=100*(XXip_imp./XXi_imp-1); %import growth
gFFdi_imp=100*(FFdip_imp./FFdi_imp-1); %import growth
gzong_imp=100*((FFdip_imp+XXip_imp)./(FFdi_imp+XXi_imp)-1); %import growth
gXXi_row_imp =100*(((sum(XXip_imp)-XXip_imp(8,:)-XXip_imp(43,:))./(sum(XXi_imp)-XXi_imp(8,:)-XXi_imp(43,:)))-1);%只有中国和美国的可以用
gFFdi_row_imp =100*(((sum(FFdip_imp)-FFdip_imp(8,:)-FFdip_imp(43,:))./(sum(FFdi_imp)-FFdi_imp(8,:)-FFdi_imp(43,:)))-1);%只有中国和美国的可以用
gzong_row_imp=100*(((sum(XXip_imp)-XXip_imp(8,:)-XXip_imp(43,:))+(sum(FFdip_imp)-FFdip_imp(8,:)-FFdip_imp(43,:)))./((sum(XXi_imp)-XXi_imp(8,:)-XXi_imp(43,:))+(sum(FFdi_imp)-FFdi_imp(8,:)-FFdi_imp(43,:)))-1); %import growth
%计算中国总进口变化
gXXi_world_imp=(sum(XXip_imp)'./sum(XXi_imp)'-1); %import growth,注意，XXip_imp行表示出口国，列表示进口国
gFFdi_world_imp=(sum(FFdip_imp)'./sum(FFdi_imp)'-1); 
gzong_world_imp=(sum(FFdip_imp+XXip_imp)'./sum(FFdi_imp+XXi_imp)'-1); %import growth
chn_gzong_imp = [gXXi_world_imp(8);gFFdi_world_imp(8);gzong_world_imp(8)];
%计算中国总出口变化
gXXi_world_exp=(sum(XXip_imp')'./sum(XXi_imp')'-1); %import growth,注意，XXip_imp行表示出口国，列表示进口国
gFFdi_world_exp=(sum(FFdip_imp')'./sum(FFdi_imp')'-1); 
gzong_world_exp=(sum(FFdip_imp'+XXip_imp')'./sum(FFdi_imp'+XXi_imp')'-1); %import growth
chn_gzong_exp = [gXXi_world_exp(8);gFFdi_world_exp(8);gzong_world_exp(8)];
table1= [chn_gzong_imp';chn_gzong_exp'];
%计算中国出口价格变化
shareF_chukou=sum(FFdip_exp(7*56+1:7*56+23,:)')/sum(sum(FFdip_exp(7*56+1:7*56+23,:))');
pF_chukou=sum(shareF_chukou'.*PrFjHat(1:23,8))-1;
XX_chukou=zeros(56,56);
for i=1:N
    XX_chukou=XX_chukou+XXip_exp(1+(i-1)*56:i*56,7*56+1:8*56);    
end

shareXX_chukou=XX_chukou(1:23,1:23)/sum(sum(XX_chukou(1:23,1:23))');
PrsjHat(:,7*56+23)=ones(56,1);
pXX_chukou=sum(sum(shareXX_chukou.*PrsjHat(1:23,7*56+1:7*56+23))')-1;

%% 模型有效性检验 0830新增。
%计算中美双边中间品出口和最终品出口占总贸易的比重。(注意，计算货物贸易）
VAnp = Youtput_test-sum(XXip)';
for i=1:N
    z_FFdi_imp(i,:)=sum(FFdi_exp(1+(i-1)*S:(i-1)*S+23,:));
    z_FFdip_imp(i,:)=sum(FFdip_exp(1+(i-1)*S:(i-1)*S+23,:));
    VAnp_c(i,1)=sum(VAnp(1+(i-1)*S:i*56),1);
    for j= 1:N
    z_XXi_imp(i,j)=sum(sum(XXi_exp(1+(i-1)*S:23+(i-1)*S,1+(j-1)*S:23+(j-1)*S),2));
    z_XXip_imp(i,j)=sum(sum(XXip_exp(1+(i-1)*S:23+(i-1)*S,1+(j-1)*S:23+(j-1)*S),2));
    end
end
tottrade_bilat_sharechn= z_FFdi_imp(8,43)/sum(z_FFdi_imp(8,:)+z_XXi_imp(8,:),2);
tottrade_bilat_sharechn_zhongjianpin= z_XXi_imp(8,43)/sum(z_FFdi_imp(8,:)+z_XXi_imp(8,:),2);
tottrade_bilat_shareame= z_FFdi_imp(43,8)/sum(z_FFdi_imp(43,:)+z_XXi_imp(43,:),2);
tottrade_bilat_shareame_zhongjianpin= z_XXi_imp(43,8)/sum(z_FFdi_imp(43,:)+z_XXi_imp(43,:),2);

tottrade_bilat_sharechnp= z_FFdip_imp(8,43)/sum(z_FFdip_imp(8,:)+z_XXip_imp(8,:),2);
tottrade_bilat_sharechn_zhongjianpinp= z_XXip_imp(8,43)/sum(z_FFdip_imp(8,:)+z_XXip_imp(8,:),2);
tottrade_bilat_shareamep= z_FFdip_imp(43,8)/sum(z_FFdip_imp(43,:)+z_XXip_imp(43,:),2);
tottrade_bilat_shareame_zhongjianpinp= z_XXip_imp(43,8)/sum(z_FFdip_imp(43,:)+z_XXip_imp(43,:),2);
%计算中间品和最终品货物贸易占GDP的比重 44国。用中间品和最终品出口
%VAn为原数据中的GDP
XXi_gdp = sum(z_XXi_imp,2)./VAn;
FFdi_gdp = sum(z_FFdi_imp,2)./VAn;

XXip_gdp = sum(z_XXip_imp,2)./VAnp_c;
FFdip_gdp = sum(z_FFdip_imp,2)./VAnp_c;
