加载中…
个人资料
微博
  • 博客等级：
  • 博客积分：
  • 博客访问：
  • 关注人气：
  • 获赠金笔：0支
  • 赠出金笔：0支
  • 荣誉徽章：
正文 字体大小：大 中 小

MATLAB实现LMDI

(2017-07-11 17:08:13)
分类： IDA(IndexDecomposionAnalysis
分割线 加性函数========================================================
function out = LMDI_add(C,X,style)
%本函数计算LMDI的总效应Dtot，活动效应Dact，结构效应Dstr和强度效应Dint.
%C表示被研究的对象数据，如碳排，能量等，X表示因素数据，如经济总量，经济结构，消费强度等。
%在计算总效应时，X为空，即为LMDI(C,[],style),[]不可省略。
t=size(C,2);
disp(strcat('输出',style,':'));
if isempty(X) && (~isempty(C)||~isempty(style))
%% 判断并计算总效应
Dtot=sum(C(:,t))-sum(C(:,1));%t时期各部门消费总量与1时期消费总量之差；
out=Dtot;
elseif ~isempty(X) && (~isempty(C)||~isempty(style))
%% 判断并计算活动效应、结构效应和强度效应。
switch style
    case 'Dact'
        Q=X;
        Dact=sum(L(C(:,t),C(:,1)))*log(sum(Q(:,t))/sum(Q(:,1)));
        out=Dact;
    case 'Dstr'
        S=X;
        Dstr=sum(L(C(:,t),C(:,1)).*log(S(:,t)./S(:,1)));
        out=Dstr;
    case 'Dint'
        I=X;
        Dint=sum(L(C(:,t),C(:,1)).*log(I(:,t)./I(:,1)));
        out=Dint;
    otherwise
        disp('case选项无匹配！')
end
%可用DeltaEtot-(DeltaEact+DeltaEstr+DeltaEint)检验误差大小
else
%% 函数调用有误。
    disp('函数参数调用有误')
end
end

分割线 L函数====================================================
function L=L(x,y)
%% local function L(x,y)
L=(x-y)./(log(x)-log(y));
end
分割线 某案例 ====================================================
P=transpose([114333
117171
121121
123626
126743
128453
130756
132129
134091
135404
]);
E=transpose([3692.147499
5152.62969
6693.253419
8160.134644
10661.52595
13261.0888
19813.61982
27857.33094
39361.60549
47830.59612
]);
S=transpose([0.181506217 0.012589176 0.008599517 0.004024072 0.009887406 0.070008715 0.087148389 0.007725686 0.02828813 0.017184391 0.017280536 0.07777089 0.037348967 0.045282627 0.021382932 0.049475203 0.015764038 0.022784763 0.016115459 0.004021706 0.01405606 0.072086064 0.045106261 0.058841463 0.0757214
0.132693319 0.010600061 0.008909881 0.003362326 0.009337567 0.059350099 0.077601046 0.007109372 0.025767688 0.017206056 0.014733076 0.070935298 0.037034739 0.046357126 0.020854127 0.055335376 0.022514414 0.022486771 0.015128994 0.002943806 0.015957319 0.075996579 0.087477483 0.117756729 0.09108436
0.129937142 0.008505629 0.009604943 0.004749718 0.011997208 0.068459134 0.087278755 0.010105955 0.027764765 0.018741553 0.016894051 0.077200267 0.039739193 0.051169561 0.022321676 0.04904382 0.025875633 0.025451736 0.024046937 0.002286491 0.005776791 0.085610569 0.033666004 0.085980169 0.077792302
0.123483087 0.011146076 0.008163331 0.005978569 0.008880579 0.069016727 0.076892702 0.011214645 0.02211265 0.022172233 0.01550303 0.076120452 0.044071342 0.038905454 0.02493579 0.041165725 0.026589906 0.02782185 0.024500929 0.004154424 0.017820053 0.086995255 0.035153328 0.093820466 0.083381397
0.102690667 0.007857622 0.017012342 0.003801054 0.002715739 0.05688469 0.066352262 0.005817947 0.019051083 0.03727104 0.030841644 0.08381656 0.024364443 0.044757375 0.016304301 0.034459499 0.036945154 0.038211039 0.044995408 0.003754957 0.010345828 0.086029149 0.041043229 0.097306776 0.087370191
0.091180476 0.012796805 0.010411576 0.004634148 0.005074463 0.046200937 0.049884857 0.012599092 0.022493599 0.028210189 0.019412975 0.068827433 0.018519365 0.049030627 0.019135495 0.041467916 0.03077779 0.022721581 0.041404286 0.005389643 0.00922857 0.08975732 0.048931415 0.156128919 0.095780524
0.072685435 0.013425544 0.010539513 0.005957079 0.004190694 0.048298123 0.051666468 0.011022543 0.019995573 0.040507565 0.022879034 0.074322568 0.029321199 0.057911382 0.019421196 0.046129418 0.032543255 0.030346478 0.051925688 0.006930515 0.009142487 0.073593996 0.04626871 0.13145021 0.089525325
0.059708695 0.01177865 0.011644114 0.007509652 0.004703634 0.051034912 0.052841735 0.013425915 0.018236361 0.041244104 0.025736501 0.075712785 0.027848964 0.07461111 0.02162213 0.048221481 0.040273653 0.033162017 0.050302012 0.005959102 0.01288305 0.076596506 0.039604951 0.121777536 0.073560431
0.055332123 0.016096916 0.009317768 0.009107555 0.004302361 0.05382548 0.045316116 0.012010063 0.016597798 0.038099369 0.024065081 0.074434419 0.031979756 0.06551894 0.01955371 0.052897991 0.046857918 0.036634045 0.045153129 0.005697233 0.01090002 0.081691986 0.039326355 0.132298013 0.072866125
0.055831565 0.014053357 0.007657164 0.007793047 0.003961211 0.054918891 0.041384214 0.011706143 0.018326948 0.033414328 0.024982827 0.075563519 0.029098299 0.068750902 0.020121096 0.045833979 0.040369261 0.031221168 0.040459595 0.003416747 0.004794287 0.086544857 0.038689815 0.159776693 0.081330085
]);
I=transpose([0.000379353 0.002623701 0.001982877 0.02916343 0.001565364 0.000517993 0.000780485 0.003070864 0.001956593 0.001549375 0.001512427 0.002458407 0.002973471 0.026365363 0.010880333 0.006020744 0.005074065 0.006939539 0.003384405 0.00493289 0.00428945 0.005166309 0.001471247 0.000678048 0.001197431
0.000468547 0.002173363 0.00158526 0.027902304 0.001629578 0.000763024 0.000968214 0.002723186 0.00165192 0.001229514 0.00123948 0.001756611 0.001920778 0.017722863 0.010709775 0.00551776 0.003809734 0.005491104 0.002401703 0.003062186 0.003778012 0.004023615 0.000490012 0.001013654 0.001045104
0.000431433 0.001488977 0.00139167 0.01565566 0.001352803 0.000496815 0.000754947 0.00193013 0.001268816 0.001051461 0.001074414 0.001575647 0.001697056 0.014431596 0.009431051 0.004339645 0.003301481 0.00468262 0.002042677 0.002609298 0.00290135 0.003506416 0.000858833 0.000650513 0.000871457
0.000368577 0.001872341 0.00076698 0.018039486 0.001030817 0.000480802 0.000547846 0.001457843 0.00102856 0.001370127 0.000909031 0.001309445 0.001686681 0.016852505 0.008513956 0.005051211 0.004051239 0.005860906 0.002043182 0.002763639 0.001718852 0.00267027 0.000962771 0.000664239 0.000953794
0.00029126 0.001253352 0.000572445 0.015181181 0.000710208 0.000337963 0.000429228 0.001136124 0.000773273 0.000883828 0.000550202 0.000943647 0.001272176 0.01061956 0.006004318 0.003955138 0.002821997 0.004042586 0.001279635 0.002011208 0.001281889 0.002091652 0.000772701 0.000461463 0.000645222
0.000256218 0.000906227 0.000511072 0.005652832 0.000622662 0.000331462 0.000445391 0.000764027 0.000656217 0.000627343 0.000581281 0.000985858 0.001077852 0.009751129 0.005516943 0.003099063 0.002408313 0.003342466 0.001326573 0.001780248 0.000907735 0.002005339 0.000526918 0.000423542 0.000459428
0.000367424 0.00155723 0.000738462 0.006899065 0.001171824 0.000491837 0.00067753 0.0011706 0.001108871 0.001351143 0.000889532 0.001473532 0.001627104 0.01123031 0.006603611 0.004410654 0.003556095 0.004392116 0.002016517 0.002579893 0.001249476 0.002003175 0.000805707 0.00063844 0.000756442
0.000355156 0.00175256 0.001482259 0.007449722 0.001190209 0.000482034 0.000712161 0.001367609 0.001284183 0.001251222 0.001267258 0.001464975 0.001676854 0.013748651 0.006846637 0.005094754 0.003937223 0.005722795 0.00219584 0.002411667 0.001652853 0.003475933 0.000905411 0.000526519 0.00079022
0.000389894 0.001768284 0.0016 0.007698211 0.001581862 0.000539506 0.000758365 0.001579716 0.001432458 0.001395323 0.001317135 0.001587712 0.002002029 0.015220095 0.007364557 0.005587727 0.004124503 0.006044132 0.002399944 0.002672924 0.001851083 0.003328381 0.00108548 0.000513037 0.00077108
0.000322818 0.001643457 0.001316747 0.010105919 0.001672382 0.000419928 0.000594409 0.001288227 0.002210375 0.001222492 0.001114384 0.001456664 0.002124774 0.016460801 0.007687665 0.005201062 0.004288932 0.006366257 0.002444314 0.003099627 0.001907074 0.003988498 0.000963043 0.000495652 0.000725577
]);
C=transpose([29066.05187 13943.19377 7198.141702 49539.88972 6533.538127 15308.2929 28712.74232 10014.93706 23364.44429 11239.35207 11032.71365 80708.91302 46880.5829 503983.1417 98210.99032 125744.321 33765.57203 66746.08636 23023.73155 8374.566688 25451.64264 157210.7728 28013.87474 16842.0411 38275.39772
37536.27771 13908.79949 8527.493617 56640.7414 9186.664151 27340.63513 45361.5449 11688.46694 25698.84219 12772.14931 11025.08631 75229.2989 42947.27608 496020.2908 134840.832 184337.9558 51785.0443 74547.97444 21937.05579 5442.392501 36397.5653 184611.8497 25879.28451 72065.03637 57471.47934
45446.76621 10267.18083 10836.46579 60283.14623 13157.4371 27572.94988 53417.2531 15813.22827 28559.42221 15975.53653 14715.07117 98613.04371 54672.89059 598663.4633 170664.6675 172542.1328 69255.85588 96619.1074 39821.37175 4836.70915 13587.62178 243359.0526 23439.96064 45343.08931 54959.05895
45913.77953 21053.00146 6316.236618 108799.9222 9234.856982 33475.56093 42496.29571 16493.14371 22944.44143 30646.2471 14216.81336 100553.1637 74988.77438 661427.2631 214171.497 209767.6074 108670.523 164496.973 50500.62663 11582.41718 30899.73281 234346.1454 34142.58268 62867.88727 80228.88159
40416.11389 13307.84189 13159.55159 77974.5582 2606.254941 25978.10654 38484.61597 8931.799294 19906.55496 44512.57627 22929.94137 106876.5859 41883.90754 642265.3521 132284.5779 184167.9289 140882.6121 208733.1937 77803.24485 10204.81786 17920.86818 243152.3657 42854.50416 60676.88798 76175.53186
39795.51342 19754.30819 9064.038249 44623.06883 5382.277123 26086.0444 37847.19967 16397.25921 25143.73388 30146.33438 19222.1527 115584.4121 34002.33196 814415.0476 179829.7999 218910.248 126262.347 129368.6448 93562.07099 16344.21438 14269.77378 306606.3491 43919.21146 112642.6842 74958.07249
69189.55158 54163.90861 20163.88967 106475.3406 12722.53804 61542.74756 90690.64538 33428.40902 57443.43563 141795.6761 52725.99125 283730.3717 123601.1256 1684927.202 332263.7388 527116.2361 299819.4043 345308.7185 271274.8774 46322.57473 29594.95995 381932.6117 96580.66469 217423.5906 175447.3339
78053.92037 75981.18174 63528.43445 205919.5368 20606.0383 90548.81757 138513.8517 67583.9565 86199.12689 189947.6453 120047.2085 408260.2997 171886.5295 3775732.709 544895.7263 904276.8496 583644.8348 698532.7408 406559.1409 52897.59314 78377.36868 979981.6519 131987.4766 236003.5343 213958.619
113866.5805 150233.6111 78687.25298 370053.0988 35920.94228 153269.8956 181385.7915 100137.5297 125488.7399 280585.41 167297.7336 623760.7907 337923.1572 5263282.138 760060.5632 1560079.84 1020063.147 1168668.304 571954.4136 80375.38108 106494.1476 1435109.165 225308.6644 358239.9589 296549.8893
116727.9953 149580.7747 65299.13745 510058.9783 42904.29031 149359.7223 159315.3669 97665.99397 262357.4941 264555.0636 180307.4995 712867.48 400421.7792 7329369.95 1001805.364 1543891.979 1121339.723 1287272.803 640495.1323 68589.80954 59214.59739 2235568.247 241312.4515 512893.6827 382183.7941
]);
for i=2:10
pop(i-1)=LMDI_add(C(:,1:i),P(:,1:i),'Dact');
act(i-1)=LMDI_add(C(:,1:i),E(:,1:i),'Dint');
str(i-1)=LMDI_add(C(:,1:i),S(:,1:i),'Dstr');
int(i-1)=LMDI_add(C(:,1:i),I(:,1:i),'Dint');
end
其他函数：
分割线 乘数函数==============================================
function out = LMDI_mult(C,X,style)
%本函数计算LMDI的总效应Dtot，活动效应Dact，结构效应Dstr和强度效应Dint.
%C表示被研究的对象数据，如碳排，能量等，X表示因素数据，如经济总量，经济结构，消费强度等。
%在计算总效应时，X为空，即为LMDI(C,[],style),[]不可省略。
t=size(C,2);
disp(strcat('输出',style,':'))
if isempty(X) && (~isempty(C)||~isempty(style))
%% 判断并计算总效应
Dtot=sum(C(:,t))/sum(C(:,1));
out=Dtot;
elseif ~isempty(X) && (~isempty(C)||~isempty(style))
%% 判断并计算活动效应、结构效应和强度效应。
switch style
    case 'Dact'
        Q=X;
Dact=exp(sum(L(C(:,t),C(:,1)).*log(sum(Q(:,t))/sum(Q(:,1))))/L(sum(C(:,t)),sum(C(:,1))));
out=Dact;
    case 'Dstr'
        S=X;
Dstr=exp(sum(L(C(:,t),C(:,1)).*log(S(:,t)./S(:,1)))/L(sum(C(:,t)),sum(C(:,1))));
out=Dstr;
    case 'Dint'
        I=X;
Dint=exp(sum(L(C(:,t),C(:,1)).*log(I(:,t)./I(:,1)))/L(sum(C(:,t)),sum(C(:,1))));
out=Dint;
    otherwise
        disp('case选项无匹配！')
end
else
%% 函数调用有误。
    disp('函数参数调用有误')
end
end
分割线 MH函数==========================================
function out = MH(CH,XH,H,style)
%计算数据为三层，可计算两层之间的结构效应的LMDI方法。
%H为一行向量，表示数据层次状况，H=[12,4,12],表示C和X的Level1有三类，每类又可以分为12子类，4子类和12子类。
%% 将Level2数据分类加总为Level1.
%计数Level1的类别；.
L1=size(H,2);
%循环计算，得到Level1的数据矩阵
for i=1:L1
    if i=1
       C(i,:)=sum(CH(1:H(i),:));
        X(i,:)=sum(XH(1:H(i),:));
    elseif i~=1 && i<=L1
       C(i,:)=sum(CH(1+H(i-1):H(i),:));
        X(i,:)=sum(XH(1+H(i-1):H(i),:));
    else
        disp('L1有误！')
    end
end
%% 计算level2和level2之间的结构效应Dstr12和level2层次的强度效应DintL2
switch style
    case 'Dstr12'
t=size(C,2);
%内层sigma求和
for i=1:L1
   if i==1
      cont2=sum(L(CH(1:H(i),t)/X(i,t),CH(1:12,1)/X(i,1))/L(C(i,t)/X(i,t),C(i,1)/X(i,1))*log(XH(1:H(i),t)/X(i,t)/XH(1:12,1)/X(i,1)));
   elseif i~=1 && i<=L1
      r=sum(H(1:i-1)):sum(H(1:i))%表示第i个level1中所有level2数据的所在行
      const2=const2+...同i下对所有j的有关计算求和
      sum(...对第i个level1计算求和
      L(CH(r,t)./X(i,t),CH(r,1)./X(i,1))...L函数
      /L(C(i,t)/X(i,t),C(i,1)/X(i,1))...除以另一L函数
      .*log(XH(r,t)/X(i,t)/XH(r,1)/X(i,1))...乘以某log函数
      );
   else
       disp('Dstr12计算有误！');
   end
%外层Sigma求和
   const1=const1+...
           L(C(i,t),C(i,1))...
          /L(sum(C(:,t),sum(C(:,1))))...
          *const2;
end
    case 'DintL2'
        t=size(C,2);
        %内层sigma求和
        for i=1:L1
           if i==1
              cont2=sum(L(CH(1:H(i),t)/X(i,t),CH(1:12,1)/X(i,1))/L(C(i,t)/X(i,t),C(i,1)/X(i,1))*log(XH(1:H(i),t)/X(i,t)/XH(1:12,1)/X(i,1)));
           elseif i~=1 && i<=L1
              r=sum(H(1:i-1)):sum(H(1:i))%表示第i个level1中所有level2数据的所在行
              const2=const2+...同i下对所有j的有关计算求和
              sum(...对第i个level1计算求和
              L(CH(r,t)./X(i,t),CH(r,1)./X(i,1))...L函数
              /L(C(i,t)/X(i,t),C(i,1)/X(i,1))...除以另一L函数
              .*log(XH(r,t)/X(i,t)/XH(r,1)/X(i,1))...乘以某log函数
              );
           else
               disp('Dstr12计算有误！');
           end
        %外层Sigma求和
           const1=const1+...
                   L(C(i,t),C(i,1))...
                  /L(sum(C(:,t),sum(C(:,1))))...
                  *const2;
        end 
end
%% local function L(x,y)
function L=L(x,y)
%% local function L(x,y)
L=(x-y)./(log(x)-log(y));
end
end
分割线 PDM函数=============================================
%% Parametric Divisia method
% 数据格式：
% （1）E为矩阵，行为时间轴，列为部门轴，代表能源消耗。
% （2）Y为矩阵，行为时间轴，列为部门轴，代表产量。
% （3）S为矩阵，行为时间轴，列为部门轴，代表部门的产量份额。
% （4）I为矩阵，行为时间轴，列为部门轴，代表能源强度。
E=[];
Y=[];
S=[];
I=[];%输入数据；
t=size(E,2);%表示能源消耗的时间长度。
alpha=[];beta=[];tau=[];%设定参数；
%% 模型AWT-PDM的权重设置
% alpha=(sum(I(:,1))*(sum(Y(:,t))-sum(Y(:,1)))-sum(E(:,1))*log(sum(Y(:,t))/sum(Y(:,1))))/((sum(I(:,1))-sum(I(:,t)))*(sum(Y(:,1))-sum(Y(:,t)))-(sum(E(:,1))-sum(E(:,t)))*log(sum(Y(:,t))/sum(Y(:,1))));
% for i=1:size(E,1)
% beta(i)=(I(i,0)*sum(Y(:,1))*(S(i,t)-S(i,0))-E(i,0)*log(S(i,t)/S(i,1)))/((I(i,0)*sum(Y(:,1))-I(i,t)*sum(Y(:,t)))*(S(i,t)-S(i,t))-(E(i,0)-E(i,t))*log(S(i,t)/S(i,0)));
% end
% for i=1:size(E,1)
%     tau(i)=(Y(i,0)*(I(i,t)-I(i,0))-E(i,0)*log(I(i,t)/I(i,0)))/((Y(i,0)-Y(i,t))*(I(i,t)-I(i,0))-(E(i,0)-E(i,t)*log(I(i,t)/I(i,0))));
% end
%% Parametric Divisia method I（PDM I）：
DeltaEpdn = (sum(E(1))+alpha*(sum(E(t))-sum(E(1))))*log(sum(Y(:,t))/sum(Y(:,1)));
DeltaEstr=0;
for i=1:size(E,1)
DeltaEstr=DeltaEstr+(E(i,1)+beta(i)*(E(i,t)-E(i,1)))*log(S(i,t)/S(i,1));
end
DeltaEint=0;
for i=1:size(E,1)
DeltaEint=DeltaEint+(E(i,1)+tau(i)*(E(i,t)-E(i,1)))*log(I(i,t)/I(i,1));
end
DeltaEtot = DeltaEpdn+DeltaEstr+DeltaEint+D;%能量消耗可以分解为：产出效应、结构效应、强度效应和残差。
disp([DeltaEtot, DeltaEpdn,DeltaEstr,DeltaEint,D]);
%% Parametric Divisia method II (PDM II)
t=size(E,2);%表示能源消耗的时间长度。
DeltaEpdn = (sum(I(1))+alpha*(sum(I(t))-sum(I(1))))*(sum(Y(:,t))-sum(Y(:,1)));
DeltaEstr=0;
for i=1:size(E,1)
DeltaEstr=DeltaEstr+I(i,1)*sum(Y(:,1))+beta(i)*(I(i,t)*sum(Y(:,t))-I(i,1)*sum(Y(:,1)))*(S(i,t)-S(i,1));
end
DeltaEint=0;
for i=1:size(E,1)
DeltaEint=DeltaEint+(Y(i,1)+tau(i)*(Y(i,t)-Y(i,1)))*(I(i,t)-I(i,1));
end
DeltaEtot = DeltaEpdn+DeltaEstr+DeltaEint+D;%能量消耗可以分解为：产出效应、结构效应、强度效应和残差。
disp([DeltaEtot, DeltaEpdn,DeltaEstr,DeltaEint,D]);

喜欢

0

阅读┊ 收藏 ┊ 喜欢▼ ┊打印┊举报/Report
  

新浪BLOG意见反馈留言板　欢迎批评指正

新浪简介 | About Sina | 广告服务 | 联系我们 | 招聘信息 | 网站律师 | SINA English | 产品答疑

Copyright © 1996 - 2022 SINA Corporation, All Rights Reserved

新浪公司 版权所有