diary on ; % 使用当前 diary 日志文件名启用记录。
disp(['*********开始时间:',datestr(now,'yyyy-mm-dd HH:MM:SS'),'**********************************']);
*********开始时间:2025-07-08 14:27:51**********************************

%% 1.导入数据和参数

Input; %直接选中命令，右键点击后按选择执行即可。如果想查看命令内容，右键点击后选择"打开Input"即可，里面包含了数据和参数的导入程序。

%% 2.模拟矩方法求解相关参数T β α γ（fun_SMM为程序自编写，按说明文件选择对应命令右键点击选择运行即可）

tic 

[T,beta_new,alpha,gamma] = fun_SMM(par,data);

toc
历时 281.704590 秒。

%% 3.输出结果
% 对应表中θ ρ σ τ γ T κ α β 展示
fprintf('theta: %.2f\n' , theta);
theta: 4.50
fprintf('rho: %.2f\n'   , rho);
rho: 0.55
fprintf('sigma: %.2f\n' , sigma);
sigma: 4.00
fprintf('tau: %.2f\n'   , mean(mean(tau)));
tau: 1.76
fprintf('gamma: %.2f\n' , mean(mean(gamma)));
gamma: 2.00
fprintf('T: %.2f\n'     , mean(T));
T: 0.13
fprintf('kappa: %.2f\n' , K);
kappa: 0.10
fprintf('v: %.2f\n'     , v);
v: 0.56
fprintf('alpha: %.2f\n' , alpha);
alpha: 0.87
fprintf('beta: %.2f\n'  , beta_new);
beta: 2.22



Calibration       = [T;beta_new;alpha];
gamma_tau         = [gamma;tau];
[spsi, spsi_SOE, x_iln, x_SOE_iln ,Deficit,lambdaE,lambdaE_SOE,psi,psi_SOE] ...
                  = get_data_sim(s,tau,gamma,M,w,alpha,beta,MSOE,TSOE,K,T,X);
data.spsi         = spsi;
data.spsi_SOE     = spsi_SOE;
data.x_iln        = x_iln;
data.x_SOE_iln    = x_SOE_iln;
data.deficit      = Deficit;
data.X            = X;
data.T_X          = T_X;
data.lambdaE      = lambdaE;
data.lambdaE_SOE  = lambdaE_SOE;
data.psi          = psi;
data.psi_SOE      = psi_SOE;

save('Data','data');
save('Calibration','Calibration');
save('gamma_tau','gamma_tau');

%% 4.参数校准结果检验（运行后可以得到在stata生成对应图形的数据，运行方法同上）
T          = Calibration(1:m,1);
beta       = ones(m,1);
beta(s)    = Calibration(m+1,1);
alpha      = Calibration(m+2,1);
gamma      = gamma_tau(1:m,1:m);
tau        = gamma_tau(m+1:2*m,1:m);
[spsi, spsi_SOE, x_iln, x_SOE_iln ,Deficit,lambdaE,lambdaE_SOE] ...
           = get_data_sim(s,tau,gamma,M,w,alpha,beta,MSOE,TSOE,K,T,X);
MP         = [sum(x_iln+x_SOE_iln,3),sum(xiln,3)];
TR         = [permute(sum(x_iln+x_SOE_iln,1),[2,3,1]),permute(sum(xiln,1),[2,3,1])];
writematrix(MP,'MPCompare.txt');
writematrix(TR,'TRCompare.txt');


%% 5.矩条件的敏感性检验
h = 1e-10;
for i = 1:m
      for n = 1:m
        T_up    =  T;   
        T_up(i) = (1+h)*T(i);
        [spsi_or, spsi_SOE_or, x_iln_or, x_SOE_iln_or ,Deficit_or,lambdaE_or,lambdaE_SOE_or] = get_data_sim(s,tau,gamma,M,w,alpha,beta,MSOE,TSOE,K,T,X);
        [spsi_up, spsi_SOE_up, x_iln_up, x_SOE_iln_up ,Deficit_up,lambdaE_up,lambdaE_SOE_up] = get_data_sim(s,tau,gamma,M,w,alpha,beta,MSOE,TSOE,K,T_up,X);
        X_up  = permute(sum(sum(x_iln_up+x_SOE_iln_up,2),1),[1,3,2])';
        X_or  = permute(sum(sum(x_iln_or+x_SOE_iln_or,2),1),[1,3,2])';
        MP_up = sum(x_iln_up+x_SOE_iln_up,3);
        MPShare_up(i,n) = MP_up(i,n)./X_up(i);
        MP_or = sum(x_iln_or+x_SOE_iln_or,3);
        MPShare_or(i,n) = MP_or(i,n)./X_or(i);
        DeriveMP_T(i,n) = 1/h*(MPShare_up(i,n) - MPShare_or(i,n));
      end
end

[f1,xi1] = ksdensity(DeriveMP_T(:));

h = 1e-10;
for i = 1:m
      for n = 1:m
        T_up    =T;   
        T_up(i) = (1+h)*T(i);
        [spsi_or, spsi_SOE_or, x_iln_or, x_SOE_iln_or ,Deficit_or,lambdaE_or,lambdaE_SOE_or] = get_data_sim(s,tau,gamma,M,w,alpha,beta,MSOE,TSOE,K,T,X);;
        [spsi_up, spsi_SOE_up, x_iln_up, x_SOE_iln_up ,Deficit_up,lambdaE_up,lambdaE_SOE_up] = get_data_sim(s,tau,gamma,M,w,alpha,beta,MSOE,TSOE,K,T_up,X);;
        X_up  = permute(sum(sum(x_iln_up+x_SOE_iln_up,2),1),[1,3,2])';
        X_or  = permute(sum(sum(x_iln_or+x_SOE_iln_or,2),1),[1,3,2])';
        TR_up = permute(sum(x_iln_up+x_SOE_iln_up,1),[3,2,1])' ;
        TRShare_up(i,n)=TR_up(i,n)./X_up(i);
        TR_or = permute(sum(x_iln_or+x_SOE_iln_or,1),[3,2,1])';
        HGShare_or(i,n)=TR_or(i,n)./X_or(i);
        DeriveTR_T(i,n) = 1/h*(TRShare_up(i,n) - HGShare_or(i,n));
      end
end
[f2,xi2] = ksdensity(DeriveTR_T(:));

writematrix([f1',xi1',f2',xi2'],'sensitive.txt');
%% 
disp(['*********结束时间:', datestr(now,'yyyy-mm-dd HH:MM:SS'),'**********************************']);
*********结束时间:2025-07-08 14:33:56**********************************
diary off; %结束运行
