付舉望， 孔新兵

(南京審計大學 統(tǒng)計與數(shù)據(jù)科學學院，南京 211815)

0 引 言

統(tǒng)計學中，相關關系能夠通過相關系數(shù)進行度量，而因果關系卻很難有一個明確的度量指標，由此衍生出的因果推斷(Causal Inference)問題成為統(tǒng)計學者們關注的焦點。Rubin[1]在反事實理論基礎上構(gòu)建了潛在結(jié)果模型(Rubin Causal Model)，其核心是比較同一個研究對象在接受干預和不接受干預時的結(jié)果差異，即該干預的因果效應(Causal Effects)。在該反事實框架下，因果推斷問題轉(zhuǎn)變?yōu)榉词聦嵵档墓烙媶栴}。目前，關于受數(shù)據(jù)限制較小，能夠服務于高維數(shù)據(jù)，且估計效果更好的因果推斷方法，還有待進一步研究。

Rubin關注單個協(xié)變量情形下平均處理效應的估計問題，通過處理組與控制組的分離與再匹配估計反事實值； Heckman等[2-3]使用雙重差分法(Difference-in-Difference) 估計反事實值，并將其應用于社會公共政策評估；Abadie 等[4]提出合成控制法(Synthetic Control Method)，將控制組個體加權(quán)，合成一個與處理組相似的虛擬組，通過比較干預前后的處理組和虛擬組的變量變化，得出平均處理效應, 改進了雙重差分法的內(nèi)生性問題；Zheng 等[5]又在此基礎上，使用機器學習二次規(guī)劃法來確定控制組的權(quán)重并重構(gòu)虛擬組，以預測反事實值；Doudchenko 等[6]比較了雙重差分法、合成控制法、約束回歸法、最優(yōu)子集選擇法對于估計參數(shù)的不同限制條件，并使用彈性網(wǎng)法(Elastic Net)設置懲罰項，以此篩選控制變量構(gòu)造模型，獲得反事實值預測。這些方法都是基于觀測到的面板數(shù)據(jù)特點進行模型假設，導致其受限于所得觀測數(shù)據(jù)，在面對具有不同特征的數(shù)據(jù)時穩(wěn)健性較弱。

因子分析作為常見的宏觀經(jīng)濟變量預測方法，面對高維數(shù)據(jù)表現(xiàn)優(yōu)異且能應用于不同特征的數(shù)據(jù)。因此在反事實框架下，學者們通過非隨機觀測數(shù)據(jù)，在因子分析的視角下進行因果推斷，提出了新的反事實值估計方法。Xiong 等[7]提出帶有缺失值情況下的潛在因子模型，其通過行列調(diào)整后的協(xié)方差矩陣估計獲得公共因子與因子載荷，由此得到反事實值估計；Bai 等[8]提出TW算法(Tall-Wide algorithm)，其假設面板數(shù)據(jù)具有強因子結(jié)構(gòu)，并將數(shù)據(jù)劃分為bal、tall、wide、miss 4塊(block)，分別利用tall-block 估計公共因子、wide-block估計旋轉(zhuǎn)后的因子載荷，從而得到平均處理效應的估計。這些方法都是在因子分析前對面板數(shù)據(jù)進行調(diào)整，將缺失值所在行、列或是周圍一整塊數(shù)據(jù)丟棄，導致已有觀測數(shù)據(jù)信息無法完全利用。

本文從優(yōu)化的角度提出L2因子分析方法，獲得平均處理效應估計，且避免了信息丟失問題。與Xiong 等不同，該優(yōu)化方法的顯著優(yōu)勢在于勿需調(diào)整面板數(shù)據(jù)的行和列，也勿需為了構(gòu)造滿足奇異值分解的子矩陣而整行或整列地丟棄數(shù)據(jù)；而與Bai等的區(qū)別在于：該優(yōu)化方法將除缺失值以外的所有數(shù)據(jù)作為整體進行因子分析，而非將整個面板數(shù)據(jù)劃分后根據(jù)特定數(shù)據(jù)塊分別進行潛在因子與因子載荷估計，勿需舍棄未使用數(shù)據(jù)塊的信息。該方法僅需舍棄面板中的缺失值，即可通過優(yōu)化一步得到潛在數(shù)據(jù)生成效應，即潛在結(jié)果，從而提高估計效率。另外，與上述文獻關心平均處理效應不同，本文還通過引入L1損失函數(shù)并優(yōu)化L1風險，得到中位數(shù)處理效應的估計。

1 處理效應模型

本文以政策評估為例，簡述反事實框架與隨機對照實驗。假設某個城市(i)在某個時刻(t)被政策介入，在該地區(qū)對感興趣的某項指標yI,t進行研究?？紤]擁有關于該城市以及其他未被政策介入的城市的面板數(shù)據(jù)Y∈RN×T，其中在政策實施前的數(shù)據(jù)一共有a年，實施后的數(shù)據(jù)一共有b年,共計T年(T=a+b)。

借鑒已有因果模型變量設置方法，記城市i=1,2,…,n1為被政策介入的城市,作為實驗組(treat unit)；城市j=1,2,…,n2為未被政策介入的城市,作為控制組(control unit)，共計N個城市(N=n1+n2)。時間t=T1,T2,…,Ta,Ta+1,…,Ta+b,其中t=T1,T2,…,Ta為介入前，t=Ta+1,Ta+2,…,Ta+b為介入后，將面板數(shù)據(jù)YN×T劃分為4個部分，具體變量設置如表1所示。

表1 變量設置

(1)

由式(1)可以看出：只需要多次進行實驗估計反事實值，將其與真實觀測值的差求期望即可得到政策的因果效應。Rubin(1973)提出了一般化的估計方法框架：

(2)

即實驗組的值可以表示為控制組值的線性組合。由此一來，就可以通過對處理前實驗組和控制組進行回歸，得到參數(shù)ω=(ω1,ω2,…,ωn2)和μ，再使用估計出的參數(shù)和處理后控制組的值估計處理后實驗組的反事實預測值。這樣一來，問題就由求反事實估計值變?yōu)榍蠡貧w估計參數(shù)θ=(μ,ω)。一個簡單的想法是利用簡單最小二乘法最小化以下目標函數(shù)Q0，來得到參數(shù)ω和μ的估計:

但當θ=(μ,ω)的維數(shù)較高時，最小二乘法將導致較大的預測誤差，使得政策效應的估計變得非常不穩(wěn)定。因此，在這基礎上還需要一些假設條件對參數(shù)θ=(μ,ω)加以限制。而目前的各種方法，均是基于已有面板數(shù)據(jù)設立模型假設，在該假設下對參數(shù)θ=(μ,ω)施加條件，以獲得更好的估計效果。

2 基于因子分析的反事實值估計

2.1 因子模型及其估計

近年來，因子模型的理論和應用已經(jīng)得到了很大的完善和發(fā)展，動態(tài)因子模型常被應用于宏觀經(jīng)濟政策的評估、經(jīng)濟指數(shù)的構(gòu)建和經(jīng)濟指標的預測。統(tǒng)計學中，學者們關注因子個數(shù)的確定研究以及公共因子與因子載荷的估計方法研究。Bai[9]提出確定靜態(tài)因子個數(shù)的信息準則，該準則在保證因子模型的擬合優(yōu)度前提下，通過面板數(shù)據(jù)結(jié)構(gòu)特征得到因子個數(shù)的無偏估計；Fan等[10]討論具有條件稀疏結(jié)構(gòu)的高維協(xié)方差估計問題，通過引入POET(Principal Orthogonal compl Ement Thresholding)方法來探索這種近似因子結(jié)構(gòu)；Kong[11]提出一種局部主成分分析方法，用高頻數(shù)據(jù)確定具有時變因子載荷的連續(xù)時間因子模型的公共因子數(shù),該模型采用離散時間因子模型在收縮塊上進行局部近似。

(3)

通常，一個數(shù)據(jù)集總是由若干隨機變量的若干觀測組成。因子分析的目標就是將原始數(shù)據(jù)集進行降維，將這些觀測投射到一個低維因子空間中。這樣的投射有無數(shù)種，主成分方法希望找到這樣一種投射，可以使得數(shù)據(jù)在低維空間的投影擁有最大的方差。而L2因子分析問題則可以表述為目標函數(shù)Q1的最小化問題：

其中：Λ=(λ1,…,λn1,λn1+1,…,λN)T,F=(f1,…,fa,fa+1,…,fT)。對于目標函數(shù)Q1的最小化問題，常用方法為對協(xié)方差矩陣YTY進行奇異值分解(SVD)，尋找最大的r個特征值，再用其對應的r個特征向量構(gòu)成的矩陣做低維投影降維。

2.2 L2范數(shù)因子分析

在對帶有大量缺失值的面板數(shù)據(jù)進行因子分析時，基于最小二乘法的經(jīng)典方法會遇到穩(wěn)健性問題，難以得到滿意的公共因子與因子載荷。對于缺失值處理，Xiong 等提出帶有缺失值情況下的潛在因子模型，在協(xié)方差矩陣估計時，對所有個體重加權(quán)，刪除缺失值所在的行和列，從而將面板數(shù)據(jù)調(diào)整為滿足奇異值分解的子矩陣，再根據(jù)所得矩陣的特征向量進行潛在因子與因子載荷的估計； Bai 等將面板數(shù)據(jù)劃分為4塊(block)，分別利用tall-block 估計公共因子、wide-block 估計旋轉(zhuǎn)后的因子載荷；而本文則從優(yōu)化的角度直接舍棄目標函數(shù)中對應缺失的累加項，此時的目標函數(shù)Q2改寫為

s.t.ΛTΛ/N=IN,FFT/T為對角陣。

使用該方法的優(yōu)點是勿需對面板數(shù)據(jù)的行列進行拼湊和刪減操作，且使用了所有未缺失數(shù)據(jù)的信息。

2.3 L1范數(shù)因子分析

L2因子分析的缺陷是對離群值十分敏感，為了解決該問題，本文使用更具有穩(wěn)健性的L1因子分析進行研究。

將目標函數(shù)Q2更換為使用L1范數(shù),即可得目標函數(shù)Q3：

s.t.ΛTΛ/N=IN,FFT/T為對角陣。

2.4 交替凸優(yōu)化算法

目標函數(shù)Q1、Q2與Q3的最小化問題實質(zhì)為無約束最優(yōu)化問題。其中Q1為均值插補后的因子分析,可使用奇異值分解進行直接計算；Q2與Q3則需要使用優(yōu)化算法解決。

一種經(jīng)典的方法是使用交替凸優(yōu)化算法求解該問題。該算法的思想是:當面臨一個兩維變量的優(yōu)化問題，而該問題不是凸優(yōu)化問題因此無法求其最優(yōu)解時，可以采用交替迭代的方法，每一步將其中一維未知變量的值看作是常數(shù)(使用該變量上次迭代的取值)，來求解另一維未知量。由目標函數(shù)Q2與Q3給出式(4)、式(5)：

(4)

(5)

改寫目標函數(shù)Q2與Q3，有

Q2(Λ,F)=‖Y-Λ(m-1)F‖L2=

Q3(Λ,F)=‖Y-Λ(m-1)F‖L1=

(6)

(7)

同理將式(5)轉(zhuǎn)化為下面的N個獨立的子優(yōu)化問題:

(8)

(9)

其中:λs是Λ的第s行。式(6)—式(9)都可以采用線性規(guī)劃問題求解。在帶有缺失值的情況下，直接舍棄目標函數(shù)中的對應累加項，在上述算法中對應的做法就是直接刪除一個約束條件。

需要注意的是，交替凸優(yōu)化算法只能保證在每一步求得當前最優(yōu)解，并不能保證最后得到全局最優(yōu)解。具體算法步驟如下(以L2因子分析為例)：

Step1 初始化參數(shù)：給出Λ、Σ的初始值Λ(0)、Σ(0)。

Step2 交替凸優(yōu)化：對于迭代次數(shù)m(m=1,2,…,M)，有

將算法中L2范數(shù)改為L1范數(shù)即為交替凸優(yōu)化求解L1因子分析。關于初始值的選取，Σ為對角矩陣，使用單位矩陣作為初始值Σ(0)；而Λ可以采用簡單隨機數(shù)進行初始化，這里為了加快收斂速度，使用均值插補后通過奇異值分解(SVD)算法得到的因子載荷矩陣作為初始值Λ(0)。

2.5 預測效果評價

由于反事實值永遠無法得到真實的觀測值，所以在隨機對照實驗中無法獲得預測誤差，從而無法比較各個方法的預測效果。常用的解決方法是Abadie提出的安慰劑檢驗法，其在因果研究中用于檢驗反事實估計量是否具有穩(wěn)健性。安慰劑(placebo)源于醫(yī)學上的隨機實驗。研究者為檢驗藥物是否有效，將實驗人群隨機分為服用真藥的實驗組與服用安慰劑的控制組，通過不讓實驗者知道自己服用的是真藥還是安慰劑，來避免主觀心理作用的影響。以此為基礎的安慰劑檢驗，其核心思想在于從控制組中選取偽實驗組，并用相同的方法估計虛擬的“反事實值”，這樣能同時獲得真實的觀測值與估計的預測值，即可對預測結(jié)果進行評價。由此得到的安慰劑效應(即真實值與虛擬反事實值之差)越趨于0，說明其與該政策的因果效應差距越大，也就說明估計方法越穩(wěn)健。在實證研究中，安慰劑檢驗受到了廣泛使用。Abadie通過假定控制組城市受到政策影響估計安慰劑效應，并作圖比較，說明合成控制法的穩(wěn)健性;Athey等[13]則分別假設虛擬實驗組與虛擬政策實施時間，對雙重差分、帶懲罰項的橫向遞減法(horizontal regression)、矩陣補全方法(Matrix Completion Methods)等多種方法進行了比較。

本文具體考慮以下兩種情況：

(1) 隨機選取控制組城市j*，假設其在t=Ta+1,Ta+2，…,Ta+b時間段中受到政策介入影響，成為偽實驗組(pseudo-treat unit)。其余控制組城市(j(.))仍然為控制組，政策介入時間不變。具體變量設置由表2所示。

表2 偽實驗組假設

(2) 隨機選取時刻Ta+c與Ta+d(其中Ta+c滿足時間順序Ta+1,Ta+2，…,Ta+c,…,Ta+d;1

表3 偽時間假設

在以上兩種情形中，同樣地將“介入后”的“實驗組”數(shù)據(jù)視作缺失值，用上節(jié)的方法進行因子分析，并將得到的L1因子和L2因子以及對應因子載荷建立預測模型，預測“反事實值”。在安慰劑檢驗中，能同時獲得真實的觀測值與估計的預測值，對預測結(jié)果進行評價。下面選取3種指標：FMSE、FMAE和FMPAE，并由這3種指標比較各方法的預測精度。計算方法如式(10)所示(eit為實驗組城市i在t時刻的預測誤差)：

(10)

3 實證分析

3.1 數(shù)據(jù)選取與變量設置

本文選取因果推斷研究中部分學者使用的關于加利福尼亞州限制吸煙政策的數(shù)據(jù)(Abadie等,2010; Doudchenko等,2016; Athey等,2021)。使用該數(shù)據(jù)的優(yōu)勢在于：可以通過本文使用方法得出的結(jié)果與以往已有估計方法得到的結(jié)果進行比較，對因果效應是否存在進行基本驗證。在該數(shù)據(jù)中，加利福尼亞州于1988年被限制吸煙政策介入，作為實驗組；選取另外38個未被任何吸煙管控政策介入的州作為控制組。同時選取了1970—2000年共計31a間各州的煙草銷量數(shù)據(jù)，并設定1988年為政策介入時刻(Ta)，該政策于1989年(Ta+1)起對煙草銷量產(chǎn)生因果效應。

3.2 模型設置

實證研究中，面板數(shù)據(jù)Y為39行、31列的矩陣。在具體因子分析模型中，N=39，n1=1，n2=38;T=31，a=16，b=15。采用Bai等提出的信息準則確定靜態(tài)因子個數(shù)。在該信息準則下，因子個數(shù)需要最小化：

其中，V(Λ,F)為因子殘差平方和，G(N,T)為懲罰函數(shù)，其使得在N,T→∞時，G(N,T)→0，且min(N,T)G(N,T)→∞。參考Bai&Ng(2002)建議，本文在實證研究中選擇式(11)作為懲罰函數(shù)：

(11)

之后對目標函數(shù)Q2與Q3利用交替凸優(yōu)化算法進行迭代，分別得到對應的潛在因子Ft與因子載荷Λs；再通過Ft與Λs進行估計，獲得加州在沒有被政策介入情況下的反事實預測值；最后由式(1)得到加州限制吸煙政策在不同方法下計算出的政策效應。

3.3 預測結(jié)果分析

實證研究中，除了本文因子分析相關的3種方法外，使用因果推斷中常見的雙重差分法(Difference-in-Difference，DID)以及Doudchenko(2016)使用的彈性網(wǎng)絡法(Elastic Networks)作為比較，由式(4)計算所得的政策效應(表4)。

表4 政策效應估計

表4給出了用不同方法得到的限制吸煙政策對于煙草價格的因果效應。可以看出政策效應均為正，且本文所使用的3種方法效應更顯著。具體每年的預測值如圖1所示。

圖1 各方法每年預測值

其中，線y為加州煙草銷售額的實際變化。由圖1可以看出，雙重差分法(線yhatDID)全時間段估計效果均較差，不能說明政策效應是否顯著。經(jīng)過均值插補后的簡單因子分析(線yhat)雖然在1989年政策介入后的估計與實際值有顯著差異，能夠說明政策的因果效應，但在政策介入前的時間段(1970—1988)與實際值相差較大。而L1因子分析(線yhatL1)、L2因子分析(線yhatL2)以及彈性網(wǎng)絡方法(線yhatEN)均能在政策介入前與實際值保持一致，且在介入后與實際值表現(xiàn)出顯著差距。

3.4 預測誤差分析

偽實驗組假設中，隨機抽取5個原控制組州作為偽實驗組，其余33個州仍為控制組。政策介入時刻、介入前后總時間不變，由式(10)計算所得誤差如表5所示。

表5 偽實驗組假設預測誤差

由表5可以看出：通過均值插補后的簡單因子分析得出的結(jié)果誤差較大，而L1、L2因子分析的效果均不錯，且L1范數(shù)較L2范數(shù)具有更小的平均預測誤差。具體到每個偽實驗組的平均預測誤差，L1、L2因子分析所得的箱線圖如圖2所示。

圖2 偽實驗組假設預測誤差

由圖2可以看出：L1因子分析離群值點更少，最大值也更小，總體而言預測誤差也小于L2因子分析，可以認為該情況下，L1范數(shù)預測效果略好于L2范數(shù)。

偽時間假設情況中，將1999年作為偽政策偽介入時刻(Ta+c)，該政策于2000年(Ta+c+1)起對煙草銷量產(chǎn)生因果效應?？刂平M與實驗組不變，選取1989—2019年共計31a間各州的煙草銷量數(shù)據(jù)，由式(10)計算所得誤差如表6所示。

表6 偽時間假設預測誤差

由表6可以看出：通過均值插補后的簡單因子分析得出的結(jié)果誤差較大，而L1、L2因子分析的效果均不錯，且L1范數(shù)與偽實驗組假設一樣具有更小的預測誤差。具體到實驗組每年的預測誤差，L1、L2因子分析所得的箱線圖如圖3所示。

圖3 偽時間假設預測誤差

由圖3可以看出:L1因子分析離群值點情況與L2因子分析相差無幾，而總體預測誤差明顯小于L2因子分析，可以認為該情況下L1范數(shù)預測效果同樣好于L2范數(shù)。綜上所述，經(jīng)均值插補后的簡單因子分析表現(xiàn)較差，各個誤差顯著大于另外兩種方法。而L1、L2因子分析的預測效果在不同情況下表現(xiàn)優(yōu)秀，且L1因子分析在兩種假設下效果均略好于L2因子分析。通過比較，可以認為L2因子分析在進行因果推斷時，雖然具有一定穩(wěn)健性，但較之于L1因子分析仍有所欠缺。

4 總 結(jié)

首先簡要介紹了因果推斷的提出與發(fā)展歷程。通過已有文獻，概述了因果推斷的相關模型，確定其最終目的是反事實值估計，并介紹了已有的估計反事實值的方法與理論模型；之后引入因子模型的基礎理論，包括基本概念和模型參數(shù)的估計方法:正交因子模型主要通過主成分分析法來估計，而主成分估計得出的因子得分可以很好地應用于宏觀經(jīng)濟變量預測；由此，結(jié)合因果推斷與正交因子模型，將因果推斷反事實值估計轉(zhuǎn)變?yōu)閹в腥笔е档臐撛谝蜃幽Ｐ凸烙嫛?/p>

基于以上理論依據(jù)，舍棄面板數(shù)據(jù)中的缺失值，通過優(yōu)化一步得到潛在結(jié)果與平均處理效應；并采用L1因子分析代替L2因子分析來估計模型，做出穩(wěn)健性上的改進，獲得中位數(shù)處理效應；通過對L1因子分析的問題進行表述，介紹了一種交替凸優(yōu)化算法并給出其實現(xiàn)步驟。

最后，為了驗證L1因子分析能否代替L2因子分析作為因子模型的估計，基于加利福尼亞州限制煙草政策案例作了實證研究，將兩種不同的主成分估計應用在反事實值的預測上。實證研究的結(jié)果表明：因子模型的L1估計量同樣適用于宏觀經(jīng)濟變量預測。后又通過偽實驗組與偽介入的假設，分析比較了L1、L2因子分析的預測效果，結(jié)果表明L1因子分析較之其他方法具有更穩(wěn)健的預測效果。