本文主要介紹BSTS模型原理以及CausalImpact對模型的代碼實現，旨在面對一些具有特定周期性特點的數據時，更精準科學地進行因果效應值的估計。
作者簡介
Yiwen，攜程數據分析師，專注用戶增長、因果推斷、數據科學等領域。
一、背景
如何科學地推斷某個產品策略對觀測指標產生的效應非常重要，這能夠幫助產品和運營更精準地得到該策略的價值，從而進行后續方向的迭代及調整。
在因果推斷框架下，效果評估的黃金準則一定是“AB實驗”，因為實驗的分流被認為是完全隨機且均勻的，在此基礎上對比實驗組與對照組的指標差異就可以體現某個干預帶來的增量值。但是很多場景下，我們較難進行嚴格的AB實驗，例如對于酒店的定價；現金獎勵的發放等等，不適宜向不同人群展現不同的內容。對于這些問題，我們會采取因果推斷的方法來進行策略的效果評估。
本文主要介紹BSTS模型原理以及CausalImpact對模型的代碼實現，旨在面對一些具有特定周期性特點的數據時，更精準科學地進行因果效應值的估計。下文將首先對模型原理進行簡要闡釋；隨后利用模擬數據展示代碼邏輯，最后在具體的業務場景中進行實踐。
二、現有方法及潛在問題
大部分運營和產品在評估一些全量上線的策略效果時，最常用的方式就是看上線前后的效果差異。但這種方法最大的問題在于其假設前提：假設上線的功能是唯一影響效果的變量（即沒有任何其他干預和混淆變量），但這個假設現實中往往很難實現。
于是我們嘗試使用更多因果推斷的方法，例如PSM（傾向分匹配法），在所有非實驗組的用戶群中，找到與實驗組用戶的特征非常相似的一群人，將他們的指標數據（例如下單率，訂單收益等等）與實驗組的用戶進行對比，從而體現出干預帶來的影響。但這個方法較為依賴選取的用戶特征與最后的匹配效果。
再例如SCM（合成控制方法），利用一些未受干預的地區合成一個“類似的虛擬地區”來與“上線策略的地區”進行整體的對比。但這也需要一個關鍵假設：可以找到長期變化趨勢高度同步的地區來進行對照，而這個條件往往也很難實現。
進而在傳統SCM的基礎上，我們企圖通過類似集成學習的方法，將多個未干預的對照組作為輸入值，再結合實驗組自身長期的時間序列波動情況，擬合出一個未受干預的虛擬對照組，從而將“對照組與實驗組高度同步”的強假設降為弱假設。本文介紹的BSTS模型就是用來刻畫某種“長期的時間序列波動”的數據模型，CausalImpact是用來針對這樣的數據進行因果效應值的估計。下文中我們將詳細介紹這兩個工具。
三、模型介紹
BSTS模型 （Bayesian Structured Time Series）稱為“貝葉斯結構化時間序列”，正如其名，它的主要特點體現在：
- 適用于有結構特征的時間序列數據
- 利用貝葉斯的思想來進行參數估計
結構化的時間序列數據在日常生活中不少見，尤其像攜程這樣的OTA行業，平臺的訂單情況其實是有一定時間規律的，例如周末和節假日是訂單高峰期；周中是訂單平峰期等。另一方面，貝葉斯的思想是指在得到樣本數據之前，即對要估計的參數有一些“先天認知”），隨后基于這樣的認知，結合樣本數據再得到后驗分布（如下方公式展示）
故BSTS模型主要即對結構化時序數據進行模型擬合及預測，在擬合的過程中使用到了貝葉斯的先驗思想。其好處是能夠給出預測值的置信區間，使得預測結果更科學可信。下文將對這幾種思想逐一進行介紹。
3.1 狀態空間模型
結構化的時間序列數據是指某一觀測數據的背后其實隱藏著隨時間變化而變化的不同狀態，其中觀測值與狀態值之間有對應關系；不同時刻的狀態之間也有轉換關系。我們一般用以下狀態空間模型來刻畫這兩種映射邏輯：
(1) 稱為觀測方程，反映觀測值與其背后隱藏狀態的關系；(2) 稱為狀態方程，反映隨時間推移各個狀態之間的轉換。；
都是不同變量之間的“關系映射矩陣”；
是獨立于其他變量且服從正態分布的噪聲。所謂數據的“結構化”，主要包括：
- Linear Local Trend（局部趨勢）：一定時間內的單調性（單調上升或下降）
- Seasonality（季節性因子）：固定長度的變化，類似于一年四季的溫度變化
- Cyclical（周期性）：類似季節性但波動時間不固定，波動頻率也不固定的變化
圖3-1：觀測數據及其結構化元素。第一張圖體現原數據的波動情況；第二張體現季節性因子的情況；第三張圖體現局部趨勢的情況。
如果希望在映射關系中加入協變量X，可以將(1)拓展為：
其中
表示協變量X與觀測數據之間的關系，如果協變量項表現很好（如有顯著影響）的話，那對應的local trend就會相對較弱。上述三個方程中的參數將在后文中展示估計方式。
3.2 貝葉斯及MCMC（馬爾可夫蒙特卡洛方法）
假設狀態方程(2)中各個時刻的狀態序列為
表示模型中所有的參數。我們現在希望對θ進行估計，核心步驟如下：
- 對θ設置先驗分布
以及初始狀態的分布 
- 構造馬爾科夫鏈，用MCMC方法得到
- 通過貝葉斯公式計算得到參數的后驗分布
下面對于各個步驟中用到的方法進行簡要說明：
1）貝葉斯估計：BSTS模型的一大特點就是在參數估計上使用了貝葉斯估計，即在估計之前先給出參數設置先驗分布，隨后再結合樣本數據給出參數的后驗分布。不同類型的參數一般有一些常用的先驗分布，例如均值一般使用正態分布，
，方差使用inverse-Gamma分布，
協方差矩陣可以使用IW分布等等。值得注意的是，先驗分布的設置一定程度上會影響后續MCMC收斂的情況以及后驗分布的準確性，因此并不能太過隨意地設置先驗分布，應盡可能多地根據實際數據推導出最合適的先驗分布，或是比較各先驗分布下后驗分布和似然函數的值來進行選擇。
2）MCMC方法：我們嘗試構造一條馬爾可夫鏈（一種特殊的序列，當前時刻的狀態值僅與前一時刻的狀態值有關，最終序列會收斂到某個穩定的分布），使得其最終收斂的穩態分布就是參數的后驗分布。這一過程我們可以通過Gibbs采樣實現：設置先驗分布之后，從初始狀態
出發，每次固定α采樣θ；再固定θ采樣α，逐漸一次次更新兩組參數，最終形成一條服從馬爾可夫性質的鏈路
，可以證明其穩態收斂的分布就是
，其中
代表所有的觀測數據。
3）預測值估計：得到
之后，我們從該分布中對（α，θ）進行采樣，再代入狀態空間方程(1)中對y進行預測，得到
，其中
表示時間點n之后y的預測值。
圖3-2：展示某結構化時序數據及其背后各個狀態轉換的過程。狀態α 包含 Local trend:
（局部趨勢); local level:
(局部趨勢的均值) 以及協變量x,
表示所有的觀測數據;
表示根據狀態模型得到的預測數據。
分別表示
的標準差這些參數均通過MCMC的方式得到估計。
四、模型應用與代碼實現
以上我們給出了BSTS模型及MCMC方法的簡要理論推導及結果輸出，核心目的就是對觀測值y做出預測。接下來我們將介紹如何在因果推斷場景中應用BSTS模型。
在對政策的效果評估上，我們核心想要的是觀測對象“反事實值”，例如“如果沒有這個廣告投放，用戶的瀏覽情況會怎樣？”相較于傳統的PSM或SCM方法，BSTS勝在其能夠對于時間序列數據進行效果評估；同時利用貝葉斯估計輸出反事實值y的預測，并給出預測值的置信區間，能一定程度上降低反事實值預測的波動性，提升效應評估的準確性與穩定性。
在實踐應用上，可以通過谷歌開源的CausalImpact包來實現BSTS模型，在Python/ target=_blank class=infotextkey>Python和R中均可調用，具體代碼實現詳見參考文獻[7][8]。
圖4-1：展示執行代碼時的三個步驟：訓練BSTS模型；反事實值預測；計算因果效應值，包括效應值的點估計及置信區間。
4.1 代碼實現
下面通過模擬數據展示代碼的具體命令
import tensorflow as tf
import tensorflow_probability as tfp
from causalimpact import CausalImpact
# 模型初始化 - 自定義時間序列數據：
def plot_time_series_components(ci):
component_dists = tfp.sts.decompose_by_component(ci.model, ci.observed_time_series, ci.model_samples)
num_components = len(component_dists)
mu, sig = ci.mu_sig if ci.mu_sig is not None else 0.0, 1.0
for i, (component, component_dist) in enumerate(component_dists.items()):
component_mean = component_dist.mean().numpy()
component_stddev = component_dist.stddev().numpy()
# 自定義觀測方程以及真實值y：
def plot_forecast_components(ci):
component_forecasts = tfp.sts.decompose_forecast_by_component(ci.model, ci.posterior_dist, ci.model_samples)
num_components = len(component_forecasts)
mu, sig = ci.mu_sig if ci.mu_sig is not None else 0.0, 1.0
for i, (component, component_dist) in enumerate(component_forecasts.items()):
component_mean = component_dist.mean().numpy()
component_stddev = component_dist.stddev().numpy()
# 生成模擬數據，包括一個實驗組數據（有干預）以及兩條對照組數據（無干預）
observed_stddev, observed_initial = (tf.convert_to_tensor(value=1, dtype=tf.float32),tf.convert_to_tensor(value=0., dtype=tf.float32))
level_scale_prior = tfd.LogNormal(loc=tf.math.log(0.05 * observed_stddev), scale=1, name='level_scale_prior') # 設置先驗分布
initial_state_prior = tfd.MultivariateNormalDiag(loc=observed_initial[..., tf.newaxis], scale_diag=(tf.abs(observed_initial) + observed_stddev)[..., tf.newaxis], name='initial_level_prior') # 設置先驗分布
ll_ssm = tfp.sts.LocalLevelStateSpaceModel(100, initial_state_prior=initial_state_prior, level_scale=level_scale_prior.sample()) #訓練時序模型
ll_ssm_sample = np.squeeze(ll_ssm.sample().numpy())
# 整合數據
x0 = 100 * np.random.rand(100) # 對照組1
x1 = 90 * np.random.rand(100) # 對照組2
y = 1.2 * x0 + 0.9 * x1 + ll_ssm_sample #生成真實值y
y[70:] += 10 #設置干預點
data = pd.DataFrame({'x0': x0, 'x1': x1, 'y': y}, columns=['y', 'x0', 'x1'])
圖4-2：展示模擬數據。虛線表示干預發生的時間點，藍線表示受到干預的觀測數據；黃線與綠線表示沒有受到干預的兩組對照數據。
# 調用模型：
pre_period = [0, 69] #設置干預前的時間窗口
post_period = [70, 99] #干預后的窗口
ci = CausalImpact(data, pre_period, post_period) #調用CausalImpact
# 對于causalImpact的使用我們核心需要填寫三個參數：觀測數據data、干預前的時間窗口、干預后的時間窗口。
# 輸出結果：
ci.plot()
ci.summary()圖4-3：展示CausalImpact輸出的結果圖，圖1表示真實值與模型擬合值的曲線；圖2表示每個時刻真實值與預測值的差異；圖3表示真實值與預測值的累計差值。表3-1：展示CausalImpact輸出的結果表格，量化效應值effect的估計及其置信區間，反映效應值是否具有顯著性。107.71表示干預之后實際值的平均；96.25表示干預之后預測值的平均，3.28表示估計的標準差，[89.77,102.64]表示反事實估計的置信區間。11.46表示實際值與預測值的差距，[5.07,17.94]表示差值的置信區間，由于差距的置信區間在0的右側，表示干預有顯著的提升作用。
圖4-3：展示CausalImpact輸出的結果圖，圖1表示真實值與模型擬合值的曲線；圖2表示每個時刻真實值與預測值的差異，橙色陰影部分表示置信區間；圖3表示真實值與預測值的累計差值。
表3-1：展示CausalImpact輸出的結果表格，量化效應值effect的估計及其置信區間，反映效應值是否具有顯著性。107.71表示干預之后實際值的平均；96.25表示干預之后預測值的平均，3.28表示估計的標準差，[89.77,102.64]表示反事實估計的置信區間。11.46表示實際值與預測值的差距，[5.07,17.94]表示差值的置信區間，由于差距的置信區間在0的右側，表示干預有顯著的提升作用。
4.2 模型校驗
對于模型擬合的結果，我們需要進行類似AB實驗的“AA校驗”。一般可以通過圖示的結果中的第二張圖，觀察干預之前真實值與預測值差值的置信區間是否包含0，如果包含0則說明通過檢驗，模型擬合效果不錯。上圖中，置信區間均含0，說明模型可用。
4.3 模型調整
- 過程參數：我們可以使用Tensorflow中的Decomposition來查看時序模型中各個結構元素，包括周期性/季節性等等。
seasonal_decompose(data)
圖4-4展示了生成數據背后的狀態元素。第一張圖反映原數據的走勢；第二張圖反映局部趨勢因子；第三張圖反映季節性因子?？梢钥闯鰯祿嬖诩竟澬越Y構且呈單調上升趨勢。
- 自定義參數：我們可以自定義參數的先驗分布；迭代次數；周期性的時間窗口長度等等。往往參數調整會對結果輸出有影響，例如正確的選取先驗分布會讓結果更準確；迭代次數更多能保證MCMC收斂更穩定（但也可能導致模型運行時間較長）等等。其中最重要的是對時間窗口長度的設置，需要正確地反映觀測數據的周期性。如果是年維度數據以星期為周期則設置neasnotallow=52；如果是天維度數據以小時為周期則設置neasnotallow=24等等。
CausalImpact(..., model.args = list(niter = 20000, nseasons = 24))
圖4-5展示CausalImpact包中各個參數含義及其默認值。
- 自定義時序模型：causalImpact的包中默認使用BSTS模型進行訓練，我們也可以改為其他的時序模型，但前提是需要對數據進行標準化。（如果使用默認的BSTS則不一定需要標準化）
from causalimpact.misc import standardize
normed_data, _ = standardize(data.astype(np.float32)) #標準化數據
obs_data = normed_data.iloc[:70, 0]
# 使用tfp中的其他模型來訓練時序數據
linear_level = tfp.sts.LocalLinearTrend(observed_time_series=obs_data)
linear_reg = tfp.sts.LinearRegression(design_matrix=normed_data.iloc[:, 1:].values.reshape(-1, normed_data.shape[1] -1))
model = tfp.sts.Sum([linear_level, linear_reg], observed_time_series=obs_data)
# 將自義定時序模型代入CausalImpact包中
ci = CausalImpact(data, pre_period, post_period, model=model)
五、業務場景實踐
用戶營銷是促進留存及轉化的重要方式，其中對用戶進行消息觸達是一大核心手段，尤其是在節假日的購票高峰期對用戶進行推送，方式包括站內push；微信生態中的小程序訂閱消息；公眾號或是企微環境等等，目的包括但不限于提醒用戶購票；宣傳品牌功能；發放優惠券吸引用戶轉化等等。在節假日之后，我們希望對這次的營銷觸達進行效果評估。
這是一個較為典型的不適合進行AB實驗的場景，首先因為節假日是流量高峰時期，如果嚴格預留50%用戶不觸達，可能會損失一批潛在的轉化用戶；如果改將對照組預留很少的人數，例如對照組：實驗組=1:9，那對于后續的轉化對比的科學性會產生影響。其次，節假日的推送策略往往非常精細化，總量達幾十條，我們較難保證對照組用戶的“純粹性”，用戶可能會被交叉觸達。
基于種種問題，我們較難通過傳統AB 的手段來評估推送帶來的轉化效果。因此我們考慮使用因果推斷的方式來解決。常規可選的方法和潛在問題如下：
- 如果使用PSM，需要在大盤中尋找與推送人群相似但是沒有被推送的用戶作為對照組。但一般節假日推送時都會有兜底策略，幾乎覆蓋了95%以上的平臺用戶，較難從中找到符合條件但未被推送的人群來進行對照。
- 如果使用SCM，我們較難找到合適的對照組來合成。如評估度假BU的推送效果時，我們不太可能用火車、機票、酒店等各個產線合成一個“虛擬度假BU”，因為本身各個產線的用戶需求就不同，使用這樣合成的虛擬對照組來對比度假訂單的轉化率是不夠科學的。
- DID的方式也同理，我們很難找到一個滿足平行趨勢假設且業務場景相似的對照組來進行推送前后的對比。
綜上所述，一些傳統的因果推斷方法縱使在技術上可行，在業務的解釋性上也有所欠缺。而且，以上三種方式都沒有考慮到用戶購票行為的“時間周期性”。因此即使合成了對照組也不一定能夠匹配到實驗組真正的結構特點，進而導致效應值計算有偏。于是我們考慮首先驗證用戶購票的數據周期性；在定位到周期規律之后嘗試使用BSTS模型結合CausalImpact來進行反事實值的預測。下文我們選擇2022年端午的火車票營銷推送場景進行實踐。
圖5-1展示端午期間對于用戶進行不同策略的推送觸達。
5.1 數據選取
- 我們以小時為周期窗口，通過簡單的圖像能夠看出大盤的火車票下單人數確實隨著時間推移呈現某種固定趨勢。
圖5-2展示選取端午周期內（正端午前后10天）每小時的火車票大盤支付人數
- 考慮到端午作為節假日本身就有的自然流量增長，支付人數的提升不能完全歸因于推送帶來的，因此訓練時序模型的時候，選取了19年-21年所有的端午數據（正端午前后10天）輸入BSTS模型進行訓練，得到端午這個窗口內的特有的結構狀態，隨后用這個結構化的模型來代入22年的端午數據，對2022年端午推送之后的轉化人數做出預測。
- 最后使用真實的轉化人數與預測人數作差體現本次營銷推送的效果。
5.2 R-代碼實現
# 選取19-22年每年的端午窗口，按照小時劃分，共960個數據點
y_hour=c(x1,x2,x3,x4)
x_time_hour=c(1:960)
data_hour <- cbind(y_hour, x_time_hour)
pre.period <- c(1, 808) # 2022年的推送發生在第808個時間點，故以此為干預節點。
post.period <- c(809, 960)
# nseasnotallow=24, 迭代次數2000，fit the model
impact_hour <- CausalImpact(data_hour, pre.period, post.period, model.args = list(niter = 20000, nseasons = 24))
summary(impact_hour)
plot(impact_hour)
圖5-3展示使用CausalImpact返回的結果圖。第一張圖表示真實支付人數與預測支付人數；第二張圖表示真實值與預測值的差值及置信區間；第三張圖是累計差值和置信區間。
圖像顯示模型能夠通過 AA校驗，模型有效。在干預點之后，實際值較預測值有所提升，但提升的置信區間含0，因此未達到顯著程度。體現出2022年端午營銷策略對于轉化人數有所促進作用，但是效果未達顯著。
六、方法優缺點
相較于傳統因果推斷方法，BSTS模型有2個主要優點：
- 能夠識別出數據背后的結構化特征，更好的做出預測；
- 利用了貝葉斯估計的思想，得到參數的后驗分布情況，計算效應值時能夠給出置信區間。但第（2）點對于BSTS模型是一把“雙刃劍”，如果先驗分布設置得不好，會影響MCMC的收斂速度和方向甚至最終的后驗分布情況。因此對于先驗分布的選取需謹慎。
七、方法拓展
本文介紹的結構化時序模型將數據的周期特點拆分成了趨勢項、季節項、周期項等等，每種元素挨個探究。更進一步，我們可以將時間序列按照周期性的長短來進行拆分，分為長周期項（使用大滑動窗口）、短周期項（使用小滑動窗口）、季節項等等。這樣的好處是防止一些小窗口內的周期情況被長周期的信息平滑掉，能夠更好的體現出數據在不同程度上的周期特點。具體的方程可以拆分成如下形式：
其中
表示不同時間點的狀態值；4個模塊分別代表長周期項/短周期項/季節項/序列相關性項（帶有協變量X）；每個結構模塊都有一個均值和一個標準差。
圖7-1展示了某時間序列背后的4個模塊：從上至下以此表示：原數據情況；季節性因子；短周期項；相關性項；長周期項。短周期來看數據的波動比較明顯；長周期來看數據波動較不明顯，因此這里需要考慮到短周期內的數據結構，避免被長周期的數據平滑。
接下來對以上4種模塊分別進行預測。針對長周期和季節性，由于他們在短時間內的變化不大，因此可以直接使用對應方程中的μ和σ來進行預測；針對短周期項和相關性項可以通過其他機器學習方式進行預測。得到各個模塊的預測結果之后，結合各模塊特征進行融合，得到整體的預測結果。參考文獻[4]中給出了更具體的預測方式和與傳統方式的對比結果。
圖7-2展示針對長短周期不同的預測方式：長周期項與季節項可以直接用μ表示預測；短周期及協變量相關項使用自定義的機器學習模型進行預測。
依照上述方法得到的時序預測模型后，我們再將其代入CausalImpact的代碼中，調整參數model為自定義的時序模型即可。
八、總結
本文介紹了用因果推斷的方式評估某一政策作用于時間序列數據帶來的效應，使用了BSTS的狀態空間模型來進行反事實值的預測，并通過CausalImpact的代碼進行實現。
不同于其他因果推斷方法的框架，本文中的方法對所有超參數進行貝葉斯估計，再對后驗分布進行MC采樣得到反事實的預測值，主要優勢是能夠根據最大程度考慮到所有對照組以及實驗組自身的結構特點，給出反事實預測值的估計以及置信區間，衡量效應值的顯著性。
同時，本文介紹的方法主要聚焦于結構化時序數據，利用BSTS模型識別觀測數據背后的狀態值以及各個狀態之間的轉化情況，進而在進行反事實預測時，盡可能消除由隱藏狀態帶來的影響。
需要注意的是，使用BSTS模型之前，需要驗證數據是否真的有周期性特點以及結構元素是怎樣的（是否是長短周期等等），再挑選合適的時序模型來訓練；同時對于參數的先驗分布設置也需謹慎，盡可能使得最終的效應估計值科學穩定。
