clear;clc
cd '/Users/yuwenfang/Desktop/2023-01295_拟刊用稿/2023-01295+数据/2023-01295+运行数据'
addpath('/Users/yuwenfang/Desktop/2023-01295_拟刊用稿/2023-01295+数据/2023-01295+程序代码')
addpath('/Users/yuwenfang/Desktop/2023-01295_拟刊用稿/2023-01295+数据/2023-01295+日志文件')
load('Data_use.mat');
%% 中国经济内循环贸易壁垒分析 %% 
%% 中国经济内循环贸易壁垒的时间特征（图1）
figure(1)
clf
subplot(1,2,1) % 地区层面
year =1:1:7;
SDTF=plot(year,DTB_CI(1,:),'-k','LineWidth',1);
SDTF.Marker = 'o';
SDTF.MarkerFaceColor = [0 0 0];
SDTF.MarkerSize = 5;
hold on
LDTF=plot(year,DTB_CI(2,:),'-k','LineWidth',1);
LDTF.Marker = 'diamond';
LDTF.MarkerFaceColor = [0 0 0];
LDTF.MarkerSize = 5;
hold on
CDTF=plot(year,DTB_CI(3,:),'-k','LineWidth',1);
CDTF.Marker = 'square';
CDTF.MarkerFaceColor = [0 0 0];
CDTF.MarkerSize = 5;
hold off
axis([0.5,7.5,1,3.5])
legend('全国','沿海','内陆','FontSize',10,'FontName','Songti SC')
legend('NumColumns',3,'FontName','Songti SC')
ylabel('内循环贸易壁垒','Fontsize',8,'FontName','Songti SC')
set(gca,'XTick',0:1:7.5)
set(gca,'YTick',1:0.5:3.5)
set(gca,'xticklabel',{'0','1987','1992','1997','2002','2007','2012','2017'},'FontSize',10,'FontName','Times New Roman');
set(gca,'LineWidth',1,'box','off') 

subplot(1,2,2) % 部门层面
SDTF=plot(year,DTB_LC(1,:),'-k','LineWidth',1);
SDTF.Marker = 'o';
SDTF.MarkerFaceColor = [0 0 0];
SDTF.MarkerSize = 5;
hold on
LDTF=plot(year,DTB_LC(2,:),'-k','LineWidth',1);
LDTF.Marker = 'diamond';
LDTF.MarkerFaceColor = [0 0 0];
LDTF.MarkerSize = 5;
hold on
CDTF=plot(year,DTB_LC(3,:),'-k','LineWidth',1);
CDTF.Marker = 'square';
CDTF.MarkerFaceColor = [0 0 0];
CDTF.MarkerSize = 5;
hold off
axis([0.5,7.5,1,3.5])
legend('工业','劳动','资本','FontSize',10,'FontName','Songti SC')
legend('NumColumns',3,'FontName','Songti SC')
ylabel('内循环贸易壁垒','Fontsize',8,'FontName','Songti SC')
set(gca,'XTick',0:1:7.5)
set(gca,'YTick',1:0.5:3.5)
set(gca,'xticklabel',{'0','1987','1992','1997','2002','2007','2012','2017'},'FontSize',10,'FontName','Times New Roman');
set(gca,'LineWidth',1,'box','off') 
set(gcf,'position',[0,0,1000,260],'color','white'); 
fn_1 = '/Users/yuwenfang/Desktop/2023-01295_拟刊用稿/2023-01295+数据/2023-01295+日志文件/figure_1.tif';
saveas(gcf,fn_1,'tif');

%% 中国经济内循环贸易壁垒的空间特征（图2）
figure(2)
clf
M_0 = 0.75:-0.02:0.2;
M_1 = [1 1 1];
map = M_0'*M_1;
subplot(2,3,1)
NA_DTB = [Reg_dtbp(:) Reg_dtbp_Adj(:)];
[~,density,X,Y]=kde2d(NA_DTB);                      
contour(X,Y,density,75)
axis([-1 8 -1 8])
hold on 
x=-1:1:8;
y=x;
plot(x,y,'--k');
hold off
xlabel('t年本省内循环壁垒','FontSize',12,'FontName','Songti SC');
ylabel('t年邻省内循环壁垒','FontSize',12,'FontName','Songti SC');
title('（a）密度等高线（全国）','FontSize',12,'FontName','Songti SC')

subplot(2,3,2)
Coa_DTB = [Coa_dtb(:) Coa_dtb_adj(:)]; 
[~,density,X,Y]=kde2d(Coa_DTB); 
contour(X,Y,density,75)
axis([-1 8 -1 8])
hold on 
x=-1:1:8;
y=x;
plot(x,y,'--k');
hold off
xlabel('t年本省内循环壁垒','FontSize',12,'FontName','Songti SC');
ylabel('t年邻省内循环壁垒','FontSize',12,'FontName','Songti SC');
title('（b）密度等高线（沿海）','FontSize',12,'FontName','Songti SC')

subplot(2,3,3)
IN_DTB = [In_dtb(:) In_dtb_adj(:)];
[~,density,X,Y]=kde2d(IN_DTB); 
contour(X,Y,density,75)
axis([-1 8 -1 8])
hold on 
x=-1:1:8;
y=x;
plot(x,y,'--k');
hold off
xlabel('t年本省内循环壁垒','FontSize',12,'FontName','Songti SC');
ylabel('t年邻省内循环壁垒','FontSize',12,'FontName','Songti SC');
title('（c）密度等高线（内陆）','FontSize',12,'FontName','Songti SC')

subplot(2,3,4)
NAS_DTB = [Indu_dtbs(:) Indu_dtbs_Adj(:)]; 
[~,density,X,Y]=kde2d(NAS_DTB);       
contour(X,Y,density,75)
axis([-1 9 -1 9])
hold on 
x=-1:1:9;
y=x;
plot(x,y,'--k');
hold off
xlabel('t年本省内循环壁垒','FontSize',12,'FontName','Songti SC');
ylabel('t年邻省内循环壁垒','FontSize',12,'FontName','Songti SC');
title('（d）密度等高线（工业）','FontSize',12,'FontName','Songti SC')

subplot(2,3,5)
Labor_DTB = [Labor_dtbs(:) Labor_dtbs_Adj(:)];
[~,density,X,Y]=kde2d(Labor_DTB); 
contour(X,Y,density,75)
axis([-1 9 -1 9])
hold on 
x=-1:1:9;
y=x;
plot(x,y,'--k');
hold off
xlabel('t年本省内循环壁垒','FontSize',12,'FontName','Songti SC');
ylabel('t年邻省内循环壁垒','FontSize',12,'FontName','Songti SC');
title('（e）密度等高线（劳动）','FontSize',12,'FontName','Songti SC')

subplot(2,3,6)
IN_DTB = [Cap_dtbs(:) Cap_dtbs_Adj(:)]; 
[~,density,X,Y]=kde2d(IN_DTB); 
contour(X,Y,density,75)
axis([-1 9 -1 9])
hold on 
x=-1:1:9;
y=x;
plot(x,y,'--k');
hold off
xlabel('t年本省内循环壁垒','FontSize',12,'FontName','Songti SC');
ylabel('t年邻省内循环壁垒','FontSize',12,'FontName','Songti SC');
title('（f）密度等高线（资本）','FontSize',12,'FontName','Songti SC')

colormap(map)
set(gcf,'position',[0,0,1200,550],'color','white'); 
fn_2 = '/Users/yuwenfang/Desktop/2023-01295_拟刊用稿/2023-01295+数据/2023-01295+日志文件/figure_2.tif';
saveas(gcf,fn_2,'tif');

%% 中国经济内循环贸易壁垒的稳态分析（表1）
% 地区层面
WT_PIC = [machain(Reg_dtbp) machain(Coa_dtb) machain(In_dtb)];
WT_PIC = round(WT_PIC,2);
IntraColname = {'低';'中';'高'};
Varname_region ={'全国_低','全国_中','全国_高','沿海_低','沿海_中','沿海_高','内陆_低','内陆_中','内陆_高'}; 
Table_1_1 = table(WT_PIC(:,1),WT_PIC(:,2),WT_PIC(:,3),WT_PIC(:,4),WT_PIC(:,5), ...
    WT_PIC(:,6),WT_PIC(:,7),WT_PIC(:,8),WT_PIC(:,9),'RowNames',IntraColname,'VariableNames',Varname_region);
tname_1_1 = '/Users/yuwenfang/Desktop/2023-01295_拟刊用稿/2023-01295+数据/2023-01295+日志文件/table_1_1.csv';
writetable(Table_1_1, tname_1_1,'WriteRowNames',true);
% 部门层面
Sector_dtb = DTBS;
[La_dtb,Ca_dtb]=Sec_Di(Sector_dtb);
Dsec = [La_dtb;Ca_dtb];
WT_SIC = [machain(Dsec) machain(La_dtb) machain(Ca_dtb)];
WT_SIC = round(WT_SIC,2);
IntraColname = {'低';'中';'高'};
Varname_sector ={'工业_低','工业_中','工业_高','劳动_低','劳动_中','劳动_高','资本_低','资本_中','资本_高'}; 
Table_1_2 = table(WT_SIC(:,1),WT_SIC(:,2),WT_SIC(:,3),WT_SIC(:,4),WT_SIC(:,5), ...
    WT_SIC(:,6),WT_SIC(:,7),WT_SIC(:,8),WT_SIC(:,9),'RowNames',IntraColname,'VariableNames',Varname_sector);
tname_1_2 = '/Users/yuwenfang/Desktop/2023-01295_拟刊用稿/2023-01295+数据/2023-01295+日志文件/table_1_2.csv';
writetable(Table_1_2, tname_1_2,'WriteRowNames',true);

%% 中国经济外循环贸易壁垒分析 %% 
%% 中国经济外循环贸易壁垒的时间特征（图3）
figure(3)
clf
subplot(1,2,1) % 地区层面
year =1:1:5;
SDTF=plot(year,FTB_CI(1,:),'-k','LineWidth',1);
SDTF.Marker = 'o';
SDTF.MarkerFaceColor = [0 0 0];
SDTF.MarkerSize = 5;
hold on
LDTF=plot(year,FTB_CI(2,:),'-k','LineWidth',1);
LDTF.Marker = 'diamond';
LDTF.MarkerFaceColor = [0 0 0];
LDTF.MarkerSize = 5;
hold on
CDTF=plot(year,FTB_CI(3,:),'-k','LineWidth',1);
CDTF.Marker = 'square';
CDTF.MarkerFaceColor = [0 0 0];
CDTF.MarkerSize = 5;
hold off
axis([0.5,5.5,1,5])
legend('全国','沿海','内陆','FontSize',10,'FontName','Songti SC')
legend('NumColumns',3,'FontName','Songti SC')
ylabel('外循环贸易壁垒','Fontsize',8,'FontName','Songti SC')
set(gca,'XTick',0:1:5.5)
set(gca,'YTick',1:0.5:5)
set(gca,'xticklabel',{'0','1997','2002','2007','2012','2017'},'FontSize',10,'FontName','Times New Roman');
set(gca,'LineWidth',1,'box','off') 

subplot(1,2,2) % 部门层面
SDTF=plot(year,FTB_LC(1,:),'-k','LineWidth',1);
SDTF.Marker = 'o';
SDTF.MarkerFaceColor = [0 0 0];
SDTF.MarkerSize = 5;
hold on
LDTF=plot(year,FTB_LC(2,:),'-k','LineWidth',1);
LDTF.Marker = 'diamond';
LDTF.MarkerFaceColor = [0 0 0];
LDTF.MarkerSize = 5;
hold on
CDTF=plot(year,FTB_LC(3,:),'-k','LineWidth',1);
CDTF.Marker = 'square';
CDTF.MarkerFaceColor = [0 0 0];
CDTF.MarkerSize = 5;
hold off
axis([0.5,5.5,1,4])
legend('工业','劳动','资本','FontSize',10,'FontName','Songti SC')
legend('NumColumns',3,'FontName','Songti SC')
ylabel('外循环贸易壁垒','Fontsize',8,'FontName','Songti SC')
set(gca,'XTick',0:1:5.5)
set(gca,'YTick',1:0.5:4)
set(gca,'xticklabel',{'0','1997','2002','2007','2012','2017'},'FontSize',10,'FontName','Times New Roman');
set(gca,'LineWidth',1,'box','off') 
set(gcf,'position',[0,0,1000,260],'color','white'); 
fn_3 = '/Users/yuwenfang/Desktop/2023-01295_拟刊用稿/2023-01295+数据/2023-01295+日志文件/figure_3.tif';
saveas(gcf,fn_3,'tif');


%% 中国经济外循环贸易壁垒的空间特征（图4）
figure(4)
clf
M_0 = 0.75:-0.02:0.2;
M_1 = [1 1 1];
map = M_0'*M_1;
subplot(2,3,1)
NA_FTB = [Reg_ftbp(:) Reg_ftbp_Adj(:)]; 
[~,density,X,Y]=kde2d(NA_FTB); 
contour(X,Y,density,75)
axis([-1 8 -1 8])
hold on 
x=-1:1:8;
y=x;
plot(x,y,'--k');
hold off
xlabel('t年本省外循环壁垒','FontSize',12,'FontName','Songti SC');
ylabel('t年邻省外循环壁垒','FontSize',12,'FontName','Songti SC');
title('（a）密度等高线（全国）','FontSize',12,'FontName','Songti SC')

subplot(2,3,2)
Coa_FTB = [Coa_ftb(:) Coa_ftb_adj(:)]; 
[~,density,X,Y]=kde2d(Coa_FTB); 
contour(X,Y,density,75)
axis([-1 8 -1 8])
hold on 
x=-1:1:8;
y=x;
plot(x,y,'--k');
hold off
xlabel('t年本省外循环壁垒','FontSize',12,'FontName','Songti SC');
ylabel('t年邻省外循环壁垒','FontSize',12,'FontName','Songti SC');
title('（b）密度等高线（沿海）','FontSize',12,'FontName','Songti SC')

subplot(2,3,3)
IN_FTB = [In_ftb(:) In_ftb_adj(:)]; 
[~,density,X,Y]=kde2d(IN_FTB); 
contour(X,Y,density,75)
axis([-1 8 -1 8])
hold on 
x=-1:1:8;
y=x;
plot(x,y,'--k');
hold off
xlabel('t年本省外循环壁垒','FontSize',12,'FontName','Songti SC');
ylabel('t年邻省外循环壁垒','FontSize',12,'FontName','Songti SC');
title('（c）密度等高线（内陆）','FontSize',12,'FontName','Songti SC')

subplot(2,3,4)
NAS_FTB = [Indu_ftbs(:) Indu_ftbs_Adj(:)]; 
[~,density,X,Y]=kde2d(NAS_FTB);  
contour(X,Y,density,75)
axis([-1 15 -1 15])
hold on 
x=-1:1:15;
y=x;
plot(x,y,'--k');
hold off
xlabel('t年本省外循环壁垒','FontSize',12,'FontName','Songti SC');
ylabel('t年邻省外循环壁垒','FontSize',12,'FontName','Songti SC');
title('（d）密度等高线（工业）','FontSize',12,'FontName','Songti SC')

subplot(2,3,5)
Labor_FTB = [Labor_ftbs(:) Labor_ftbs_Adj(:)];
[~,density,X,Y]=kde2d(Labor_FTB); 
contour(X,Y,density,75)
axis([-1 15 -1 15])
hold on 
x=-1:1:15;
y=x;
plot(x,y,'--k');
hold off
xlabel('t年本省外循环壁垒','FontSize',12,'FontName','Songti SC');
ylabel('t年邻省外循环壁垒','FontSize',12,'FontName','Songti SC');
title('（e）密度等高线（劳动）','FontSize',12,'FontName','Songti SC')

subplot(2,3,6)
IN_FTB = [Cap_ftbs(:) Cap_ftbs_Adj(:)]; 
[~,density,X,Y]=kde2d(IN_FTB); 
contour(X,Y,density,75)
axis([-1 15 -1 15])
hold on 
x=-1:1:15;
y=x;
plot(x,y,'--k');
hold off
xlabel('t年本省外循环壁垒','FontSize',12,'FontName','Songti SC');
ylabel('t年邻省外循环壁垒','FontSize',12,'FontName','Songti SC');
title('（f）密度等高线（资本）','FontSize',12,'FontName','Songti SC')
colormap(map)
set(gcf,'position',[0,0,1200,550],'color','white'); 
fn_4 = '/Users/yuwenfang/Desktop/2023-01295_拟刊用稿/2023-01295+数据/2023-01295+日志文件/figure_4.tif';
saveas(gcf,fn_4,'tif');

%% 中国经济外循环贸易壁垒的稳态分析（表2）
% 地区层面
WT_pftb = machain(Reg_ftbp);
WT_coa_ftb = machain(Coa_ftb);
WT_inl_ftb = machain(In_ftb);
WT_PFC = [WT_pftb WT_coa_ftb WT_inl_ftb];
WT_PFC =round(WT_PFC,2);
IntraColname = {'低';'中';'高'};
Varname_region ={'全国_低','全国_中','全国_高','沿海_低','沿海_中','沿海_高','内陆_低','内陆_中','内陆_高'}; 
Table_2_1 = table(WT_PFC(:,1),WT_PFC(:,2),WT_PFC(:,3),WT_PFC(:,4),WT_PFC(:,5), ...
    WT_PFC(:,6),WT_PFC(:,7),WT_PFC(:,8),WT_PFC(:,9),'RowNames',IntraColname,'VariableNames',Varname_region);
tname_2_1 = '/Users/yuwenfang/Desktop/2023-01295_拟刊用稿/2023-01295+数据/2023-01295+日志文件/table_2_1.csv';
writetable(Table_2_1, tname_2_1,'WriteRowNames',true);
% 部门层面
Sector_ftb = FTBS;
Sector_ftb(1:2,:)=[];
[La_ftb,Ca_ftb]=Sec_Di(Sector_ftb);
Fsec = [La_ftb;Ca_ftb];
WT_SFC = [machain(Fsec) machain(La_ftb) machain(Ca_ftb)];
WT_SFC = round(WT_SFC,2);
IntraColname = {'低';'中';'高'};
Varname_sector ={'工业_低','工业_中','工业_高','劳动_低','劳动_中','劳动_高','资本_低','资本_中','资本_高'}; 
Table_2_2 = table(WT_SFC(:,1),WT_SFC(:,2),WT_SFC(:,3),WT_SFC(:,4),WT_SFC(:,5), ...
    WT_SFC(:,6),WT_SFC(:,7),WT_SFC(:,8),WT_SFC(:,9),'RowNames',IntraColname,'VariableNames',Varname_sector);
tname_2_2 = '/Users/yuwenfang/Desktop/2023-01295_拟刊用稿/2023-01295+数据/2023-01295+日志文件/table_2_2.csv';
writetable(Table_2_2, tname_2_2,'WriteRowNames',true);

%% 中国内循环贸易壁垒的分解 %% (图5) 
RNGE_REGION = NGETB_Region./(GETB_Region+NGETB_Region)*100; % 各地区非地理壁垒占比
RNGE_SECTOR = NGETB_Sector./(GETB_Sector+NGETB_Sector)*100; % 各部门非地理壁垒占比
figure(5)
clf
subplot(1,2,1)
years = 1987:5:2017;
yyaxis left; 
Ngefig_reigon=bar(years, NGETB_Region'); 
ylabel('非地理壁垒','Fontsize',12,'FontName','Songti SC'); 
ylim([1 3]);
set(gca,'ycolor','k');
set(gca,'YTick',1:0.5:3,'Fontsize',12,'FontName','Times New Roman');
colorList = [0.4 0.4 0.4;0.6 0.6 0.6; 0.9 0.9 0.9];
for i=1:1:3
  Ngefig_reigon(i).FaceColor = colorList(i,:)  ; 
end
yyaxis right;
NA = plot(years,RNGE_REGION(1,:),'-k','Linewidth',1,'Markersize',6);
NA.Marker = 'o';
NA.MarkerFaceColor = [0 0 0];
NA.MarkerSize = 5;
hold on;
Coa = plot(years,RNGE_REGION(2,:),'-k','Linewidth',1,'Markersize',6);
Coa.Marker = 'diamond';
Coa.MarkerFaceColor = [0 0 0];
Coa.MarkerSize = 5;
hold on;
Inl = plot(years,RNGE_REGION(3,:),'-k','Linewidth',1,'Markersize',6);
Inl.Marker = 'square';
Inl.MarkerFaceColor = [0 0 0];
Inl.MarkerSize = 5;
hold off;
ylim([10 70]);
set(gca,'ycolor','k');
set(gca,'YTick',10:10:70,'Fontsize',12,'FontName','Times New Roman');
ylabel('非地理壁垒占比','Fontsize',12,'FontName','Songti SC'); 
ylim([10 70]);
legend('全国','沿海','内陆','全国占比','沿海占比','内陆占比','FontName','Songti SC','NumColumns',6,'Fontsize',10,'Location','northeast')
xticks(years); 
set(gca, 'xcolor','k');

subplot(1,2,2)
years = 1987:5:2017;
yyaxis left; 
Ngefig_reigon=bar(years, NGETB_Sector'); 
ylabel('非地理壁垒','Fontsize',12,'FontName','Songti SC'); 
ylim([1 3]);
set(gca,'ycolor','k');
set(gca,'YTick',1:0.5:3,'Fontsize',12,'FontName','Times New Roman');
colorList = [0.4 0.4 0.4;0.6 0.6 0.6; 0.9 0.9 0.9];
for i=1:1:3
  Ngefig_reigon(i).FaceColor = colorList(i,:)  ; 
end

yyaxis right;
NA = plot(years,RNGE_SECTOR(1,:),'-k','Linewidth',1,'Markersize',6);
NA.Marker = 'o';
NA.MarkerFaceColor = [0 0 0];
NA.MarkerSize = 5;
hold on;
Coa = plot(years,RNGE_SECTOR(2,:),'-k','Linewidth',1,'Markersize',6);
Coa.Marker = 'diamond';
Coa.MarkerFaceColor = [0 0 0];
Coa.MarkerSize = 5;
hold on;
Inl = plot(years,RNGE_SECTOR(3,:),'-k','Linewidth',1,'Markersize',6);
Inl.Marker = 'square';
Inl.MarkerFaceColor = [0 0 0];
Inl.MarkerSize = 5;
hold off;
ylim([10 70]);
set(gca,'ycolor','k');
set(gca,'YTick',10:10:70,'Fontsize',12,'FontName','Times New Roman');
ylabel('非地理壁垒占比','Fontsize',12,'FontName','Songti SC');
ylim([10 70]);
legend('工业','劳动','资本','工业占比','劳动占比','资本占比','FontName','Songti SC','Fontsize',10,'NumColumns',6,'Location','northeast')
xticks(years); 
set(gca, 'xcolor','k');
set(gcf,'position',[0,0,1400,400],'color','white'); 
fn_5 = '/Users/yuwenfang/Desktop/2023-01295_拟刊用稿/2023-01295+数据/2023-01295+日志文件/figure_5.tif';
saveas(gcf,fn_5,'tif');