stata:データエンベロープ解析(DEA)チュートリアル
データエンベロープ解析(DEA)は、米国の著名なオペレーティング学者A.Charnes(チャーンズ)、W.W.Cooper(ライブラリプラチナ)、E.Rhodes(ローズ)が1978年に最初に提案し、相対効率評価概念に基づいて発展した非パラメータ検査方法である.この記事では、stataでDEA分析を行う方法やboostrap検査を行う方法について説明します.使用するコマンドは
DEAモデルの概要
技術効率の概念
データエンベロープ分析において、技術効率とは、1つの生産ユニット(DMU)の生産レベルが業界の技術レベルに達する程度を指す.技術効率は投入と産出の2つの角度から測定することができ、投入が既定の場合、技術効率は産出最大化の程度から測定する.生産が既定の場合、技術効率は投入の最小化の程度によって測定される.もちろん,TFPを計算する過程では,一般に既定のものが投入される.次に、技術効率の概念を理解するために、生産に投入された例を挙げます.
たんい
x x x(投入)
y y y(産出)
y/x y/x y/x
y/x y/x y/x(標準化)
A
2
1
0.5
0.625
B
3
2
0.667
0.533
C
4
3
0.75
0.938
D
5
4
0.8
1.00
E
5
2
0.4
0.5
この表において、y/x y/xは各製造単位の技術的効率の高低に反応し、y/x y/x(標準化)は各単位のy/x y/x y/xをその中の最大値で割ったものである.これは、この数値をよりよく比較するためです.複数の産出に関わる場合、各投入と産出に一定の重みを与え、それぞれ重み付けし、産出投入比を算出する.例えば、v=v 1 x 1+v 2 x 2+..+v n x n v = v_1x_1+v_2x_2+...+v_nx_n v=v1x1+v2x2+...+vnxn u = u 1 y 1 + u 2 y 2 + . . . + u n y n u=u_1y_1+u_2y_2+...+u_ny_n u=u1y1+u2y2+...+un ynがu/vu/vu/vの投入比を生成するデータエンベロープ解析は,データ自体によって重みを得る方法を議論し,各DMUの技術効率を計算することである.
ラジアル距離モデル
このコマンドの半径方向効率の測定方法は、Debreu-Farrell(Debreu 1951;Farrell 1957)メソッドを使用します.k k個のD M U DMUがあると仮定する.D M U K DMU_についてK DMUK,N N N種投入あり,x k=(x k 1,.,x k N)∈R N x_と記す.k=(x_{k 1},…,x_{kN})in R^N xk=(xk 1,…,xkN)∈RN,MM種産出あり,y k=(x k 1,..,x k M)∈R y_k =(x_{k1},...,x_{kM})\in R^M yk=(xk1,...,xkM)∈RM.そして、技術条件T T T T T T T Tの下でy y yが投入x x xから産出されると仮定し、数学表現はT={(x,y):y a r e p r o d u c i b l e b y x}T={(x,y):yareproduciblebyx}T={(x,y):y are producible by x}では、科学技術T T T T T Tの下で、生産可能セットは、P(x)={y:(x,y)∈T}P(x)={y:(x,y)in T}P(x)={y:(x,y)∈T}投入の需要セットは、P(y)={x:(x,y)∈T}P(y)={x:(x,y)in T}P(y)={x:(x,y)in T}P(y)={x:(x,y)∈T}生産可能セットを例に、技術効率は、ある所与のデータ点と生産可能セットとの境界の距離として表される.DEAモデルでこの技術効率を測定すると、k k個のD M U DMUに対して、D M U DMU毎にN N種の投入があり、M M M種が生成したデータセットにとって.Debreu–Farrell(Debreu 1951;Farrell 1957)の産出指向の推定方法は、各データ点k(k=1,2,3...K)k(k=1,2,3...K)k(k=1,2,3...K)F k o(y k,x k,y,xθ s . t . ∑ k = 1 K z k y k m ≥ y k m θ m , m = 1 , . . . , M ∑ k = 1 K z k x k n ≤ x k n θ n , n = 1 , . . . , N z k ≥ 0 F_k^o(y_k,x_k,y,x|CRS)=max\theta\\s.t.\sum_{k=1}^Kz_ky_{km}\geq y_{km}\theta_m,m=1,...,M\\\sum_{k=1}^Kz_kx_{kn}\leq x_{kn}\theta_n,n=1,...,N\\z_k\geq 0 Fko(yk,xk,y,x∣CRS)=maxθs.t.k=1∑Kzkykm≥ykmθm,m=1,...,Mk=1∑Kzkxkn≤xknθn,n=1,...,Nzk≧0 yがKである× M K\times M K×Mの生成行列,x x xはKである× N K\times N K×Nの投入行列.P(x)P(x)P(x)は最小の包囲面(smallest convex free−disposal hull)であると推定した.上記の線形計画は,規模報酬不変(CRS)の技術効率を解く.他の規模報酬に関する仮定では、z k z_を変更するだけです.k zkの制約、例えば規模報酬可変(VRS)、Σk=1 Kz k=1sum_の設定{k=1}^Kz_k=1Σk=1 K zk=1でよい.
ひラジアルこうりつモデル
この命令の非半径方向効率測定法はRussell法である.この場合、生成をガイドとする非径方向メトリックは、R M k o(y K,x K,y,x∣C R S)=m a x{M−1Σm=1 Mと定義されるθ m : ( θ 1 y k 1 , . . . , θ m y k m ∈ P ( x ) , θ m ≥ 0 , m = 1 , . . . M ) } RM_k^o(y_K,x_K,y,x|CRS)=max\{M^{-1}\sum_{m=1}^M\theta_m: (\theta_1y_{k1},...,\theta_my_{km}\in P(x),\theta_m\geq0,m=1,...M)\} RMko(yK,xK,y,x∣CRS)=max{M−1m=1∑Mθm:(θ1yk1,...,θmykm∈P(x),θm≥0,m=1,...M)}その線形計画方程式は、RMk o(yK,xK,y,x∣C R S)=M−1 m a xΣm=1 Mθ m s . t . ∑ k = 1 K z k y k m ≥ y k m θ m , m = 1 , . . . , M ∑ k = 1 K z k x k n ≤ x k n θ n , n = 1 , . . . , N z k ≥ 0 RM_k^o(y_K,x_K,y,x|CRS)=M^{-1}max\sum_{m=1}^M\theta_m\\s.t.\sum_{k=1}^Kz_ky_{km}\geq y_{km}\theta_m,m=1,...,M\\\sum_{k=1}^Kz_kx_{kn}\leq x_{kn}\theta_n,n=1,...,N\\z_k\geq 0 RMko(yK,xK,y,x∣CRS)=M−1maxm=1∑Mθms.t.k=1∑Kzkykm≥ykmθm,m=1,...,Mk=1∑Kzkxkn≤xknθn,n=1,...,Nzk≥0
ラジアルモデルに対してBoostrapを用いて仮定検査を行う
主にこの命令規模収益の2つの仮定検査を紹介し、まず、T e s t≠1:H 0:T i s g l o b a l y C R S H 1:T i s V R S Test eq 1:H_0: T\is\globally\CRS\\H_1:TisVRS Test=1:H 0:T is globally CRSH 1:T is VRS H 0 H_と仮定すると0 H 0が拒否されると、T e s t≠2:H 0’:T i s g l o b a l y N I R S H 1:T i s V R S Test eq 2:H_0': T\is\globally\NIRS\\H_1:TisVRS Test=2:H 0’:T is globally NIRSH 1:T is VRSは、この技術条件下で規模報酬が変わらないかどうかを先に確認し、NIRS(Non-Decreasing Returns to scale)であるかどうかを見なければならない.
stataコマンドの実装
tenonradialコマンド
teradialbcコマンド
teradialbcコマンドの関数サポート
stataは、
nptestindコマンド
nptestrtsコマンド
インスタンスの適用
データはCharnes,Cooper,and Rhodes(1981)に由来し,新しいコマンドの機能を説明するために人為的に新しい変数
3つの規模収益状況における技術効率を,半径方向モデルと非半径方向モデルをそれぞれ用いて測定した.次に
P値の結果は以下の通りである.
Malmquist指数
Malmquist指数の計算と分解式はここではリストされません.プログラムで例を挙げます.
PS:文章で使われている学習資料は『データエンベロープ分析方法とMaxDeaソフトウェア』と文章「Nonparametric frontier analysis using Stata」で、この資料は私の公衆番号「肖夕木の自習室」でdeaに返信して取得することができる.文の中のデータファイルは自分でstataの上でダウンロードすることができて、もちろん、私のdeaの学習資料の中にも含まれています.
tenonradial、teradialbcなどです.なお、このコマンドの演算速度とマトリクスの最大処理量はdeaコマンドよりも優れているが、いくつかの制限がある.DEAモデルの概要
技術効率の概念
データエンベロープ分析において、技術効率とは、1つの生産ユニット(DMU)の生産レベルが業界の技術レベルに達する程度を指す.技術効率は投入と産出の2つの角度から測定することができ、投入が既定の場合、技術効率は産出最大化の程度から測定する.生産が既定の場合、技術効率は投入の最小化の程度によって測定される.もちろん,TFPを計算する過程では,一般に既定のものが投入される.次に、技術効率の概念を理解するために、生産に投入された例を挙げます.
たんい
x x x(投入)
y y y(産出)
y/x y/x y/x
y/x y/x y/x(標準化)
A
2
1
0.5
0.625
B
3
2
0.667
0.533
C
4
3
0.75
0.938
D
5
4
0.8
1.00
E
5
2
0.4
0.5
この表において、y/x y/xは各製造単位の技術的効率の高低に反応し、y/x y/x(標準化)は各単位のy/x y/x y/xをその中の最大値で割ったものである.これは、この数値をよりよく比較するためです.複数の産出に関わる場合、各投入と産出に一定の重みを与え、それぞれ重み付けし、産出投入比を算出する.例えば、v=v 1 x 1+v 2 x 2+..+v n x n v = v_1x_1+v_2x_2+...+v_nx_n v=v1x1+v2x2+...+vnxn u = u 1 y 1 + u 2 y 2 + . . . + u n y n u=u_1y_1+u_2y_2+...+u_ny_n u=u1y1+u2y2+...+un ynがu/vu/vu/vの投入比を生成するデータエンベロープ解析は,データ自体によって重みを得る方法を議論し,各DMUの技術効率を計算することである.
ラジアル距離モデル
このコマンドの半径方向効率の測定方法は、Debreu-Farrell(Debreu 1951;Farrell 1957)メソッドを使用します.k k個のD M U DMUがあると仮定する.D M U K DMU_についてK DMUK,N N N種投入あり,x k=(x k 1,.,x k N)∈R N x_と記す.k=(x_{k 1},…,x_{kN})in R^N xk=(xk 1,…,xkN)∈RN,MM種産出あり,y k=(x k 1,..,x k M)∈R y_k =(x_{k1},...,x_{kM})\in R^M yk=(xk1,...,xkM)∈RM.そして、技術条件T T T T T T T Tの下でy y yが投入x x xから産出されると仮定し、数学表現はT={(x,y):y a r e p r o d u c i b l e b y x}T={(x,y):yareproduciblebyx}T={(x,y):y are producible by x}では、科学技術T T T T T Tの下で、生産可能セットは、P(x)={y:(x,y)∈T}P(x)={y:(x,y)in T}P(x)={y:(x,y)∈T}投入の需要セットは、P(y)={x:(x,y)∈T}P(y)={x:(x,y)in T}P(y)={x:(x,y)in T}P(y)={x:(x,y)∈T}生産可能セットを例に、技術効率は、ある所与のデータ点と生産可能セットとの境界の距離として表される.DEAモデルでこの技術効率を測定すると、k k個のD M U DMUに対して、D M U DMU毎にN N種の投入があり、M M M種が生成したデータセットにとって.Debreu–Farrell(Debreu 1951;Farrell 1957)の産出指向の推定方法は、各データ点k(k=1,2,3...K)k(k=1,2,3...K)k(k=1,2,3...K)F k o(y k,x k,y,xθ s . t . ∑ k = 1 K z k y k m ≥ y k m θ m , m = 1 , . . . , M ∑ k = 1 K z k x k n ≤ x k n θ n , n = 1 , . . . , N z k ≥ 0 F_k^o(y_k,x_k,y,x|CRS)=max\theta\\s.t.\sum_{k=1}^Kz_ky_{km}\geq y_{km}\theta_m,m=1,...,M\\\sum_{k=1}^Kz_kx_{kn}\leq x_{kn}\theta_n,n=1,...,N\\z_k\geq 0 Fko(yk,xk,y,x∣CRS)=maxθs.t.k=1∑Kzkykm≥ykmθm,m=1,...,Mk=1∑Kzkxkn≤xknθn,n=1,...,Nzk≧0 yがKである× M K\times M K×Mの生成行列,x x xはKである× N K\times N K×Nの投入行列.P(x)P(x)P(x)は最小の包囲面(smallest convex free−disposal hull)であると推定した.上記の線形計画は,規模報酬不変(CRS)の技術効率を解く.他の規模報酬に関する仮定では、z k z_を変更するだけです.k zkの制約、例えば規模報酬可変(VRS)、Σk=1 Kz k=1sum_の設定{k=1}^Kz_k=1Σk=1 K zk=1でよい.
ひラジアルこうりつモデル
この命令の非半径方向効率測定法はRussell法である.この場合、生成をガイドとする非径方向メトリックは、R M k o(y K,x K,y,x∣C R S)=m a x{M−1Σm=1 Mと定義されるθ m : ( θ 1 y k 1 , . . . , θ m y k m ∈ P ( x ) , θ m ≥ 0 , m = 1 , . . . M ) } RM_k^o(y_K,x_K,y,x|CRS)=max\{M^{-1}\sum_{m=1}^M\theta_m: (\theta_1y_{k1},...,\theta_my_{km}\in P(x),\theta_m\geq0,m=1,...M)\} RMko(yK,xK,y,x∣CRS)=max{M−1m=1∑Mθm:(θ1yk1,...,θmykm∈P(x),θm≥0,m=1,...M)}その線形計画方程式は、RMk o(yK,xK,y,x∣C R S)=M−1 m a xΣm=1 Mθ m s . t . ∑ k = 1 K z k y k m ≥ y k m θ m , m = 1 , . . . , M ∑ k = 1 K z k x k n ≤ x k n θ n , n = 1 , . . . , N z k ≥ 0 RM_k^o(y_K,x_K,y,x|CRS)=M^{-1}max\sum_{m=1}^M\theta_m\\s.t.\sum_{k=1}^Kz_ky_{km}\geq y_{km}\theta_m,m=1,...,M\\\sum_{k=1}^Kz_kx_{kn}\leq x_{kn}\theta_n,n=1,...,N\\z_k\geq 0 RMko(yK,xK,y,x∣CRS)=M−1maxm=1∑Mθms.t.k=1∑Kzkykm≥ykmθm,m=1,...,Mk=1∑Kzkxkn≤xknθn,n=1,...,Nzk≥0
ラジアルモデルに対してBoostrapを用いて仮定検査を行う
主にこの命令規模収益の2つの仮定検査を紹介し、まず、T e s t≠1:H 0:T i s g l o b a l y C R S H 1:T i s V R S Test eq 1:H_0: T\is\globally\CRS\\H_1:TisVRS Test=1:H 0:T is globally CRSH 1:T is VRS H 0 H_と仮定すると0 H 0が拒否されると、T e s t≠2:H 0’:T i s g l o b a l y N I R S H 1:T i s V R S Test eq 2:H_0': T\is\globally\NIRS\\H_1:TisVRS Test=2:H 0’:T is globally NIRSH 1:T is VRSは、この技術条件下で規模報酬が変わらないかどうかを先に確認し、NIRS(Non-Decreasing Returns to scale)であるかどうかを見なければならない.
stataコマンドの実装
tenonradialコマンド
tenonradial非ラジアルモデルRMを用いて技術効率を推定し、構文詳細:tenonradial outputs = inputs [(ref outputs = ref inputs)] [if] [in] [,rts(rtsassumption) base(basetype) reference(varname) tename(newvar) noprint]そのうちoutputは産出変数inputは投入変数ref outputsは産出変数の個数ref inputsは投入変数の個数rts(rtsassumption)は規模収益仮定を指定し、CRS、VRS、NIRSの3種類のbase(basetype)は最適化の方向を設定し、すなわち産出base(output)に向かい、非径方向測定の技術効率を含む投入base(input)reference(varname)向け技術参照セットtename(newvar)生成newvarを設定します.noprint推定詳細、データ記述、および参照セットをキャンセルします.teradialbcコマンド
teradialbcコマンドラジアルモデルを使用して技術効率を推定します.構文の詳細:teradialbc outputs = inputs [(ref outputs = ref inputs)] [if] [in] [,rts(rtsassumption)base(basetype) reference(varname) subsampling kappa(#) smoothed heterogeneous reps(#) level(#) tename(newvar) tebc(newvar) biasboot(newvar) varboot(newvar) biassqva(newvar)telower(newvar) teupper(newvar) noprint nodots]その違いは主にBoostrap部分にあります.詳細はhelp teradialbcを参照してください.teradialbcコマンドの関数サポート
stataは、
kdens bw()とmm quantile()をダウンロードする必要があります.nptestindコマンド
nptestind独立性検査を行うnptestrtsコマンド
nptestrts規模収益検査を行うインスタンスの適用
データはCharnes,Cooper,and Rhodes(1981)に由来し,新しいコマンドの機能を説明するために人為的に新しい変数
drefを生成する.set seed 717117
use ccr81
generate dref = x5 != 10
teradial y1 y2 y3 = x1 x2 x3 x4 x5, rts(crs) base(output) reference(dref) tename(TErdCRSo)
teradial y1 y2 y3 = x1 x2 x3 x4 x5, rts(nirs) base(output) reference(dref) tename(TErdNRSo) noprint
teradial y1 y2 y3 = x1 x2 x3 x4 x5, rts(vrs) base(output) reference(dref) tename(TErdVRSo) noprint
tenonradial y1 y2 y3 = x1 x2 x3 x4 x5, rts(crs) base(output) reference(dref) tename(TEnrCRSo) noprint
tenonradial y1 y2 y3 = x1 x2 x3 x4 x5, rts(nirs) base(output) reference(dref) tename(TEnrNRSo) noprint
tenonradial y1 y2 y3 = x1 x2 x3 x4 x5, rts(vrs) base(output) reference(dref) tename(TEnrVRSo) noprint
list TErdCRSo TErdNRSo TErdVRSo TEnrCRSo TEnrNRSo TEnrVRSo in 1/7
3つの規模収益状況における技術効率を,半径方向モデルと非半径方向モデルをそれぞれ用いて測定した.次に
nptestindを用いてboostrap検査を行うmatrix testsindpv = J(2, 3, .)
matrix colnames testsindpv = CRS NiRS VRS
matrix rownames testsindpv = output-based input-based
nptestind y1 y2 y3 = x1 x2 x3 x4 x5, rts(crs) base(output) reps(999) alpha(0.05)
matrix testsindpv[1,1] = e(pvalue)
nptestind y1 y2 y3 = x1 x2 x3 x4 x5, rts(nirs) base(output) reps(999) alpha(0.05) noprint
matrix testsindpv[1,2] = e(pvalue)
nptestind y1 y2 y3 = x1 x2 x3 x4 x5, rts(vrs) base(output) reps(999) alpha(0.05) noprint
matrix testsindpv[1,3] = e(pvalue)
nptestind y1 y2 y3 = x1 x2 x3 x4 x5, rts(crs) base(input) reps(999) alpha(0.05) noprint
matrix testsindpv[2,1] = e(pvalue)
nptestind y1 y2 y3 = x1 x2 x3 x4 x5, rts(nirs) base(input) reps(999) alpha(0.05) noprint
matrix testsindpv[2,2] = e(pvalue)
nptestind y1 y2 y3 = x1 x2 x3 x4 x5, rts(vrs) base(input) reps(999) alpha(0.05) noprint
matrix testsindpv[2,3] = e(pvalue)
matrix list testsindpv
P値の結果は以下の通りである.
CRS NiRS VRS
output-based .06906907 .25625626 .04804805
input-based .03703704 .001001 .23323323
Malmquist指数
Malmquist指数の計算と分解式はここではリストされません.プログラムで例を挙げます.
use pwt56, clear
reshape wide y k l, i(nu country) j(year)
teradial y1965 = k1965 l1965 (y1965 = k1965 l1965), rts(crs) base(output) tename(F11) noprint
teradial y1990 = k1990 l1990 (y1965 = k1965 l1965), rts(crs) base(output) tename(F21) noprint
teradial y1965 = k1965 l1965 (y1990 = k1990 l1990), rts(crs) base(output) tename(F12) noprint
teradial y1990 = k1990 l1990 (y1990 = k1990 l1990), rts(crs) base(output) tename(F22) noprint
generate mpi = sqrt(F12 / F22 * F22 / F21)
generate effch = F11 / F22
generate techch = mpi / effch
PS:文章で使われている学習資料は『データエンベロープ分析方法とMaxDeaソフトウェア』と文章「Nonparametric frontier analysis using Stata」で、この資料は私の公衆番号「肖夕木の自習室」でdeaに返信して取得することができる.文の中のデータファイルは自分でstataの上でダウンロードすることができて、もちろん、私のdeaの学習資料の中にも含まれています.