MATLAB实现LMDI分解方法（含零值处理）

影响因素分解分析常使用对数平均迪氏指数法 (Logarithmic Mean Index Method, LMDI)，其理论基础和Python实现请参考 LMDI 理论推导详解【从理论到Python-MATLAB实现（理论）】 和 LMDI 理论推导详解【从理论到Python-MATLAB实现（编程实现）】 ，0值部分的处理参考 LMDI 理论推导详解【从理论到Python-MATLAB实现（理论）】，本文仅讨论LMDI的MATLAB实现。

MATLAB参考代码如下：

% X 为多列（含1列）的自变量数据，按照时间升序排列，第一行数据为第一年，第二行数据为第二年，...
X = xlsread(Xfilepath);
% Y 为因变量数据，按照时间升序排列，第一行数据为第一年，第二行数据为第二年，...
Y = xlsread(Yfilepath);

dertX = [];
for i = 2 : length(Y)
    curdertX = [];
    for j = 1 : size(X, 2)
        curdert = (Y(i) - Y(1)) / (log(Y(i) / Y(1))) * log(X(i, j) / X(1, j));
        
        % 0 value
        if Y(i) * Y(1) * X(i, j) * X(1, j) == 0
            curdert = 0;
        end
        if Y(i) > 0 && Y(1) == 0 && X(i, j) > 0 && X(1, j) == 0
            curdert = Y(i);
        end
        if Y(i) == 0 && Y(1) > 0 && X(i, j) == 0 && X(1, j) > 0
            curdert = -Y(1);
        end
        
        curdertX = [curdertX, curdert];
    end
    dertX = [dertX; curdertX];
end

结合实例说明，选取2000~2014年间15年的数据，因变量是Y，自变量为X（共7个因素），年份为Years，参考代码如下：

%%
Y = [
    45.095; 49.556; 55.905; 64.706; 65.834; 74.213; 82.912; 90.793; ...
    98.330; 108.624; 130.310; 210.232; 212.156; 245.704; 256.856
    ];

X = [
    2.375 	0.551 	0.301 	0.522 	0.302 	0.096 	0.027
    2.259 	0.576 	0.308 	0.496 	0.311 	0.109 	0.028
    2.146 	0.553 	0.313 	0.513 	0.322 	0.122 	0.029
    2.054 	0.640 	0.310 	0.437 	0.332 	0.135 	0.031
    1.990 	0.682 	0.327 	0.388 	0.343 	0.148 	0.033
    1.944 	0.557 	0.399 	0.443 	0.354 	0.161 	0.034
    1.898 	0.516 	0.399 	0.480 	0.364 	0.174 	0.038
    1.848 	0.457 	0.408 	0.525 	0.382 	0.186 	0.040
    1.809 	0.432 	0.426 	0.533 	0.402 	0.187 	0.044
    1.781 	0.421 	0.436 	0.534 	0.424 	0.197 	0.050
    1.760 	0.407 	0.441 	0.544 	0.444 	0.214 	0.054
    1.745 	0.415 	0.413 	0.558 	0.462 	0.198 	0.058
    1.726 	0.401 	0.434 	0.554 	0.477 	0.181 	0.061
    1.697 	0.373 	0.483 	0.542 	0.492 	0.181 	0.069
    1.652 	0.370 	0.470 	0.554 	0.509 	0.160 	0.074
    ];

Years = (2000 : 2014)';

dertX = [];
for i = 2 : length(Y)
    curdertX = [];
    for j = 1 : size(X, 2)
        curdert = (Y(i) - Y(1)) / (log(Y(i) / Y(1))) * log(X(i, j) / X(1, j));
        
        % 0 value
        if Y(i) * Y(1) * X(i, j) * X(1, j) == 0
            curdert = 0;
        end
        if Y(i) > 0 && Y(1) == 0 && X(i, j) > 0 && X(1, j) == 0
            curdert = Y(i);
        end
        if Y(i) == 0 && Y(1) > 0 && X(i, j) == 0 && X(1, j) > 0
            curdert = -Y(1);
        end
        
        curdertX = [curdertX, curdert];
    end
    dertX = [dertX; curdertX];
end

outresult = [Years(2 : length(Years)), dertX];

输出结果为：

>> outresult

outresult =
  
    2001	-2.368      2.098       1.087 	-2.416 	1.389 	6.006 	1.720 
    2002	-5.101      0.182       1.967 	-0.875 	3.226 	12.057 	3.595 
    2003	-7.887      8.132       1.600 	-9.653 	5.144 	18.516 	7.503 
    2004	-9.694      11.691      4.541 	-16.261 6.978 	23.726 	10.999 
    2005	-11.705 	0.633       16.474 	-9.592 	9.286 	30.222 	13.474 
    2006	-13.922 	-4.075      17.502 	-5.209 	11.595 	36.929 	21.221 
    2007	-16.383 	-12.215 	19.862 	0.374 	15.345 	43.190 	25.666 
    2008	-18.590 	-16.615 	23.719 	1.424 	19.532 	45.532 	33.349 
    2009	-20.799 	-19.446 	26.776 	1.642 	24.520 	51.948 	44.528 
    2010	-24.066 	-24.326 	30.671 	3.315 	30.949 	64.375 	55.663 
    2011	-33.065 	-30.407 	33.934 	7.154 	45.605 	77.655 	82.020 
    2012	-34.435 	-34.282 	39.478 	6.419 	49.312 	68.413 	87.928 
    2013	-39.774 	-46.167 	55.958 	4.449 	57.750 	75.038 	111.024 
    2014	-44.186 	-48.473 	54.241 	7.242 	63.540 	62.177 	122.721 

对于部分人说结果不一样：

拜托看看前面乘以了一个系数。