卷煙創(chuàng)新產(chǎn)品工商交易預測研究

韓偉民，楊應廣

卷煙創(chuàng)新產(chǎn)品工商交易預測研究

韓偉民1，楊應廣2*

1中煙商務物流有限責任公司，交易部，北京 100055；2北京中煙信息技術有限公司，第二業(yè)務服務中心，北京 100055

對卷煙創(chuàng)新產(chǎn)品的工商交易進行預測，有助于調控工商企業(yè)微觀層面市場策略，提升特色卷煙銷量。選定影響工商交易量指標，采用灰色理論計算這些指標的關聯(lián)矩陣，據(jù)此將指標分為5類，選取每類中與工商交易量關聯(lián)度最大的指標形成指標集，運用XGBoost建立分析模型，進行結果預測。以近十年卷煙創(chuàng)新產(chǎn)品的工商交易數(shù)據(jù)為訓練集和測試集，對XGBoost預測模型進行訓練和測試，并與時間序列模型、多元回歸模型和隨機森林模型的結果進行對比。實驗表明，4種模型的預測精度均在90%以上，以XGBoost為最佳，達到了96%的精度。因此，可以將XGBoost新增為卷煙銷量預測的一種參考算法。

卷煙創(chuàng)新產(chǎn)品；工業(yè)銷量；工商交易；銷量預測；XGBoost集成學習

卷煙銷量分為工業(yè)銷量、商業(yè)銷量和社會銷量，分別是指卷煙工業(yè)公司調撥給卷煙商業(yè)公司的數(shù)量（簡稱“工商交易量”）、商業(yè)公司批發(fā)給零售戶的數(shù)量、零售戶銷售給消費者的數(shù)量。關于商業(yè)銷量的預測，已經(jīng)存在大量的參考文獻[1-3]。一些研究認為，工業(yè)公司的生產(chǎn)和銷售具有很強的計劃經(jīng)濟特征，預測的必要性不大[4]，因此相關的研究報道[5]比較少。然而，工業(yè)公司作為一個獨立的市場主體，它需要同時面對不同商業(yè)公司對不同品規(guī)需求的波動變化，仍然具有一定的市場經(jīng)濟特征，因此有必要對工商交易量進行預測。從卷煙生產(chǎn)流通的全過程來看，工商交易環(huán)節(jié)具有重要的意義，一方面在于它可以提前感受未來市場的發(fā)展方向和競爭格局；另一方面在于它對市場的表征性，協(xié)議的增減調整與合同的執(zhí)行存廢均體現(xiàn)了各工商企業(yè)的生產(chǎn)經(jīng)營理念及市場布局。因此，深入分析卷煙交易行為，無論對行業(yè)的宏觀調控，還是對工商企業(yè)微觀層面的市場策略，均具參考價值。

卷煙銷量預測的傳統(tǒng)方法有時間序列分析法[6]、灰色系統(tǒng)理論[7]、多元線性回歸[8]等。隨著機器學習研究的興起，支持向量機[2]、神經(jīng)網(wǎng)絡[3]、隨機森林[9]等方法也被用于卷煙銷量預測。為了提高預測精度，一些文獻采用多個模型進行組合預測[10]。XGBoost[11]是近幾年才出現(xiàn)的新型集成機器學習算法，性能優(yōu)秀，已經(jīng)成為眾多文獻和人工智能大賽中的熱點工具[12]。據(jù)筆者調查，尚未發(fā)現(xiàn)使用XGBoost算法進行卷煙銷量預測的文獻。

本文以“細、短、中”為代表的卷煙創(chuàng)新產(chǎn)品交易行為作為研究對象，用XGBoost算法對其交易量進行預測分析，并與其他算法進行比較，結果顯示XGBoost算法有更好的預測效果。

1 相關理論

1.1 XGBoost模型

文獻[13]提出了弱學習和強學習理論，多個弱學習器可以組合成強學習器[14]。弱學習理論的代表性算法是Gradient Boosting[15]，其優(yōu)化目標是讓上次模型殘差在損失函數(shù)梯度方向上減少。文獻[11]在此基礎上提出了XGBoost算法（eXtreme Gradient Boosting），提升了模型的預測精度。

1.2 灰色關聯(lián)聚類

2 預測模型

2.1 樣本數(shù)據(jù)

以“細、短、中”創(chuàng)新產(chǎn)品的全國工商交易量作為研究對象，數(shù)據(jù)樣本為2011—2020年的月度數(shù)據(jù)，其中，2011—2019年的數(shù)據(jù)作為訓練集，2020年作為測試集，對它們在2021年的工商交易量進行預測。

“細、短、中”創(chuàng)新產(chǎn)品中，有些品規(guī)包含了爆珠煙，其銷量較為可觀，因此這里包含了爆珠煙的分析。由于純爆珠煙的體量較小，且2020年銷量增幅（約2%）較往年偏差較大，而XGBoost模型對時間序列的預測是基于歷史數(shù)據(jù)進行的，所以會對模型的預測結果造成較大干擾。因此，在本文中，關于各類型創(chuàng)新產(chǎn)品的數(shù)據(jù)統(tǒng)計，將只要包含該屬性的產(chǎn)品即計算在內，即細支煙包含純細支和細支爆珠，短支包括純短支和短支爆珠，中支煙包括純中支和中支爆珠，爆珠煙包含純爆珠煙、細支爆珠煙、短支爆珠煙和中支爆珠煙4個類型。

2.2 指標選取

1）指標集

卷煙交易的供求關系深受價格的影響，所以價格是首選指標。同時，商業(yè)存銷比和產(chǎn)品協(xié)議量都會直接影響商業(yè)公司的產(chǎn)品訂購策略。另外，商業(yè)公司也會受到終端市場的需求影響，而宏觀層面的某些因素會在一定程度上影響終端市場對卷煙產(chǎn)品的需求波動，包括人口因素、宏觀經(jīng)濟因素和居民消費水平因素等。結合統(tǒng)計數(shù)據(jù)進行處理，將影響因素歸納如下表：

表1 工商交易量的影響因素

Tab. 1 Influence factors of industrial sales volume

影響因素指標 商品因素工業(yè)銷量（Y）、協(xié)議量（X1）、協(xié)議剩余量（X2）、調撥價（X3）、條批價（X4）、毛利率（X5）、商業(yè)庫存（X6）、存銷比（X7） 人口因素人口總數(shù)（X8）、城鎮(zhèn)人口數(shù)（X9） 宏觀經(jīng)濟因素國內生產(chǎn)總值（X10）、人均國內生產(chǎn)總值（X11） 居民消費水平城鎮(zhèn)居民可支配收入（X12）、居民消費價格指數(shù)（X13）、商品零售價格指數(shù)（X14）、社會消費品零售總額（X15）

其中，協(xié)議量：工業(yè)公司與商業(yè)公司簽訂的在一定時期內銷售卷煙產(chǎn)品的數(shù)量；

協(xié)議剩余量：協(xié)議量與工業(yè)銷量的差值，即協(xié)議量－工業(yè)銷量；

調撥價：工業(yè)公司銷售給商業(yè)公司的卷煙條價；

條批價：商業(yè)公司批發(fā)給終端零售戶的卷煙條價；

毛利率：調撥價與條批價的差值比率，即（條批價－調撥價）/調撥價；

商業(yè)庫存：商業(yè)公司的庫存量；

存銷比：商業(yè)庫存與工業(yè)銷量的比值，即商業(yè)庫存/工業(yè)銷量；

人口總數(shù)和城鎮(zhèn)人口數(shù)：均為常住人口數(shù)量。

社會消費品零售總額：批發(fā)和零售業(yè)、住宿和餐飲業(yè)以及其他行業(yè)直接售給城鄉(xiāng)居民和社會集團的消費品零售額。

影響卷煙工商交易量的因素很多，需要篩選出主要因素，因此，需要對所采集的指標集進行聚類，降低指標維度，從而找出主要指標。

2）指標篩選

將上述選取的15個指標進行灰色關聯(lián)聚類，先計算其關聯(lián)矩陣，矩陣部分結果如下：

表2 指標關聯(lián)矩陣部分結果

Tab. 2 Partial results of index correlation matrix

X1X2X3X4X5 X110.810.620.570.59 X2--10.550.560.52 X3----10.790.86 X4------10.82 X5--------1

表3 指標集歸類結果

Tab. 3 Categorization results of index sets

類別指標 A類X7，X2，X6 B類X1 C類X3，X4，X5 D類X10，X11，X12 E類X8，X9，X13，X14，X15

分別計算各個指標與工業(yè)銷量的綜合關聯(lián)度，選取每一類別中與工業(yè)銷量灰色關聯(lián)度最大的指標，結果如下：

表4 指標與工業(yè)銷量Y綜合關聯(lián)度

Tab. 4 Comprehensive relevance between indexes and industrial sales volume Y

A類B類C類 X7X2X6X1X3X4X5 Y0.570.510.820.650.690.730.89 D類E類 X10X11X12X8X9X13X14X15 Y0.530.560.590.590.570.760.790.81

2.3 參數(shù)選取

為了使得預測模型的結果達到最優(yōu)，需要對XGBoost模型中的參數(shù)進行調優(yōu)，既可以防止模型過擬合，也可以提升模型預測效果，一般分為3類參數(shù)：通用參數(shù)、提升參數(shù)和學習任務參數(shù)。在本文中，通過采用遍歷的方式，對提升參數(shù)中的部分參數(shù)進行修正，分別為學習速率（learning_rate）、樹的最大深度（max_depth）、最小葉子節(jié)點（min_child_weight）、樹的棵數(shù)（n_estimators）和損失臨界值（gamma），從而提升模型預測效果。

對指定區(qū)間內的各個參數(shù)，借助工具Python按給定步長進行遍歷，分別對各類型的創(chuàng)新產(chǎn)品進行數(shù)據(jù)分析，得到調優(yōu)后的參數(shù)值，如下表所示：

表5 XGBoost模型參數(shù)

Tab. 5 Parameters of the XGBoost model

參數(shù)類別初始值細支煙調優(yōu)值短支煙調優(yōu)值中支煙調優(yōu)值爆珠煙調優(yōu)值 learning_rate提升參數(shù)0.10.050.10.050.08 max_depth提升參數(shù)35564 min_child_weight提升參數(shù)55576 n_estimators其他參數(shù)100100909080 gamma提升參數(shù)100.10.20.3

3 預測結果分析

根據(jù)所選取的指標集及調優(yōu)后的模型，對樣本進行訓練及預測。為了驗證本模型的更優(yōu)性，使用相同的樣本數(shù)據(jù)分別對時間序列模型（ARIMA）、多元回歸模型和隨機森林模型進行訓練和預測，具體結果如下：

3.1 調優(yōu)結果分析

根據(jù)調優(yōu)的XGBoost模型，各類型創(chuàng)新產(chǎn)品2020年四個季度的銷量預測值如表6所示?？梢?，大部分情況下，預測值比實際銷量偏低，各季度的誤差絕對值范圍控制在0.2%～4%之間，較為平穩(wěn)，因此該模型對測試集的預測效果較好。

表6 2020年各季度銷量預測結果表

Tab.6 Sales forecast results for each quarter in 2020 （箱，%）

第一季度第二季度第三季度第四季度 中支煙實際值412493326925396544273737 調優(yōu)模型預測值405141324255393546270831 誤差-1.78%-0.82%-0.76%-1.06% 細支煙實際值131104210098091081860652840 調優(yōu)模型預測值129936310086031060422644760 誤差-0.89%-0.12%-1.98%-1.24% 短支煙實際值18162713456614564790343 調優(yōu)模型預測值17869613353514436987226 誤差-1.61%-0.77%-0.88%-3.45% 爆珠煙實際值408014288109312088190411 調優(yōu)模型預測值408975284014304155193035 誤差0.24%-1.42%-2.54%1.38%

2020年各月的預測結果如圖1所示，默認參數(shù)的XGBoost模型的預測值較真實值普遍偏小，調參后的XGBoost預測模型的預測精度較高，修正誤差約為4%。

3.2 不同模型預測結果比較

將XGBoost模型與3個常用預測模型——ARIMA、多元線性回歸和隨機森林——進行對比，分別預測各類型創(chuàng)新產(chǎn)品2020年的銷量值，并進行分析比較。

1）ARIMA時間序列模型

適用于ARIMA模型分析預測的時間序列必須是平穩(wěn)非白噪聲序列，需要對原始數(shù)據(jù)進行差分去周期性等平穩(wěn)化處理；然后，計算樣本的自相關系數(shù)（）和偏自相關系數(shù)（），從而估計自相關階數(shù)和移動平均階數(shù)的值。

由下表7可知，各類型創(chuàng)新產(chǎn)品的ARIMA模型結果。其中，各參數(shù)值均小于0.05，調整2均大于0.9；以及各模型的統(tǒng)計值約為1.78，接近于2，在樣本容量足夠大時，可以認為不存在殘差自相關性。這一系列檢驗值表明了模型的擬合程度較好。

表7 ARIMA模型參數(shù)結果表

Tab.7 ARIMA model parameters

產(chǎn)品類型pdq 細支煙211 短支煙111 中支煙321 爆珠煙210

2）多元線性回歸模型

對選定的最終指標集關于銷量建立最優(yōu)的多元線性回歸方程，采用普通最小二乘法估計模型中的未知參數(shù)，根據(jù)判定系數(shù)（2）來判斷回歸方程對實際觀測值的擬合度效果，數(shù)值越大，擬合程度越好，趨勢線的可靠性越高。該預測模型的參數(shù)如下表所示：

表8 多元線性回歸模型變量系數(shù)表

Tab.8 Coefficients of multiple linear regression model variables

產(chǎn)品類型X6X1X5X15R2 細支煙-63.01620.108392.85620.09420.8621 短支煙-112.73120.0752258.85110.00210.8328 中支煙-27.98270.1184101.21010.10180.9155 爆珠煙-93.05610.0942216.01720.01120.8432

3）隨機森林模型

表9 隨機森林模型參數(shù)

Tab. 9 Parameters of the random forest model

產(chǎn)品類型mtryn_tree 細支煙4100 短支煙290 中支煙3100 爆珠煙380

4）預測結果比較

選用各產(chǎn)品的平均誤差（）和2作為該預測模型的評價指標，其中是預測值與真實值之間的誤差，其值越小預測效果越好；而2是趨勢線擬合程度的指標，其值越大預測效果越好。

在各對比模型使用最優(yōu)參數(shù)的情況下，對各類型產(chǎn)品2020年1—12月的銷量進行預測，結果如下表所示，調優(yōu)XGBoost模型的和2的值效果最好，尤其對于中支煙和細支煙，擬合程度達到0.97以上，預測誤差控制小于2%。因此，調優(yōu)XGBoost模型的預測結果相對較好，其次是隨機森林模型的預測效果。

表10 各預測模型的MAPE及判定系數(shù)

Tab. 10 MAPE and determination coefficients of various prediction models

模型系數(shù)細支煙短支煙中支煙爆珠煙 調優(yōu)XGBoostMAPE1.311%2.061%1.581%2.007% R20.98010.96770.97230.9693 ARIMAMAPE1.682%5.301%4.902%8.031% R20.97390.94210.95320.9122 多元回歸MAPE8.791%9.731%7.931%9.347% R20.86210.83280.91550.8432 隨機森林MAPE1.739%1.812%2.104%3.105% R20.97020.97010.96680.9599

以中支煙為例，各模型的誤差趨勢如下圖所示，除了在2月和8月這兩個時間點，其他月份里調優(yōu)XGBoost的誤差值均小于其他預測模型的結果，一直圍繞0上下波動。所以，從誤差絕對值的比較也可得到調優(yōu)XGBoost模型的預測效果優(yōu)于其他模型的結論。

圖2 中支煙各預測模型誤差值趨勢圖 Fig.2 Error trend of middle cigarettes by various prediction models

3.3 預測趨勢分析

對2021年各類型創(chuàng)新產(chǎn)品的交易進行趨勢分析，結果如圖3所示，各類型的創(chuàng)新產(chǎn)品將在2021年迎來迅猛發(fā)展，2021年創(chuàng)新產(chǎn)品全國工商間的交易總量有望突破800萬箱，預測其整體銷量增幅將不低于13個百分點。因此，在2021年，創(chuàng)新產(chǎn)品的市場需求仍然較為旺盛，市場發(fā)展空間相對較大。

根據(jù)近5年創(chuàng)新產(chǎn)品的工商交易量可知，短支煙的年銷量增幅一直呈現(xiàn)下降趨勢，目前已經(jīng)下降至10%以內了，年銷量增長趨勢逐漸趨于平穩(wěn)；雖然中支煙和細支煙的增幅也在逐漸下跌，但其年銷量增長趨勢仍較為明顯，根據(jù)預測結果，中支煙和細支煙在2021年繼續(xù)保持迅猛的發(fā)展勢頭，均能保持兩位數(shù)的增長速率，尤其中支煙更為明顯。

爆珠煙的整體銷量在逐年提升，增幅在逐年下降，不過，因為大部分爆珠煙是結合中支和細支進行銷售的，因此理論上其增量會從中支和細支的銷售增量中受益。預測結果顯示，2021年中支煙和細支煙會保持迅猛的增長勢頭，相應會帶來爆珠煙的銷售增量，由于爆珠煙的現(xiàn)有銷量較中支煙和細支煙差距較大，故推斷類似的增量會導致爆珠煙有更大的增幅。模型預測顯示，2021年爆珠煙的全年銷量增幅為19.88%，高于去年的實際值13.55%。2021年行業(yè)一季度經(jīng)濟運行數(shù)據(jù)顯示爆珠煙同比增長28.66%，超出了全年預測值19.88%約9個百分點，這一方面說明了模型的預測增長趨勢同真實情況是相符的，另一方面，根據(jù)表6所示近5年的第一季度增幅和全年增幅及其偏差來看，前者通常要比后者高出較多，因此，2021年全年預測值比第一季度實際值低9個百分點，一個重要原因就是第一季度因為包含元旦和春節(jié)，具有較大的特殊性。

4 結論

本文針對近十年卷煙創(chuàng)新產(chǎn)品的工商交易數(shù)據(jù)提出了一種基于XGBoost的預測方法，并對比分析了時間序列模型、多元線性回歸模型和隨機森林模型的預測效果。實驗表明，參數(shù)調優(yōu)后的XGBoost模型具有更好的預測效果，精度超過了96%。在后續(xù)研究中，可基于類似的思路和方法，推廣到具體卷煙品規(guī)的預測。

[1] 王詩豪，張曉妮，張云，等. 銅川市卷煙需求集成預測[J]. 中國煙草學報，2019, 25(6): 105-109.

WANG Shihao, ZHANG Xiaoni, ZHANG Yun, et al. Integrated forecast of cigarette demand in Tongchuan City[J]. Acta Tabacaria Sinica, 2019, 25(6): 105-109.

[2] 武牧，林慧蘋，李素科，等. 一種基于支持向量機的卷煙銷量預測方法[J]. 煙草科技，2016, 49(2): 87-91.

WU Mu, LIN Huiping, LI Suke, et al. An SVM-based method for predicting cigarette sales volume[J]. Tobacco Science & Technology, 2016, 49(2): 87-91.

[3] 齊志成. 基于 BP 神經(jīng)網(wǎng)絡模型的商洛市卷煙需求預測[J].湖南農(nóng)業(yè)科學，2017, (1): 86-89.

QI Zhicheng. Cigarette Demand Forecasting Based on BP Neural Network Model in Shangluo[J]. HUNAN AGRICULTURAL SCIENCES, 2017, (1):86-89.

[4] 沈秋云. 卷煙商業(yè)銷量預測方法的研究[D]. 上海: 上海交通大學, 2018.

SHEN Qiuyun. AN INVESTIGATION OF CIGARETTE SALES FORECAST[D]. Shanghai: Shanghai Jiao Tong University, 2018.

[5] 華勇. 基于工商協(xié)同的卷煙月度投放量預測模型[J]. 中國煙草學報, 2015, 21(4): 99-106.

HUA Yong. Monthly cigarette supply forecasting model based on coodination of producdtion and marketing[J]. Acta Tabacaria Sinica, 2015, 21(4): 99-106.

[6] 張笑通. 基于時間序列的卷煙需求預測模型的研究與應用[D]. 鄭州: 鄭州大學，2016.

ZHANG Xiaotong. Research and Application of Cigarette Demand Forecasting Model based Time Series[D]. Zhengzhou: Zhengzhou University, 2016.

[7] 袁妍，楊帆. 基于灰色關聯(lián)的卷煙銷售額預測及提升路徑研究[J]. 管理觀察，2015, (10): 129-131.

YUAN Yan, YANG Fan. [J]. Management Observer, 2015, (10): 129-131.

[8] 李瑩. 基于時間序列與多元線性回歸綜合模型的農(nóng)村卷煙銷量預測[D]. 昆明: 云南大學, 2015.

LI Ying. Based on time sequence and multiple linear regression model of rural cigarette sales forecasting[D].

[9] 鄒雯. 基于改進隨機森林的卷煙訂購量預測的研究與應用[D]. 南昌：南昌大學，2020.

ZOU Wen. Research and application of cigarette order forecasting based an improved random forests[D]. Nanchang: Nanchang University, 2020.

[10] 趙旻，張丹楓，曾中良，等. 基于組合模型的云南省卷煙需求預測與結果評價研究[J]. 中國煙草學報，2019, 25(1): 93-98.

ZHAO Min, ZHANG Danfeng, ZENG Zhongliang, et al. Prediction of cigarette demand in Yunnan province based on combination model and criteria of result evaluation[J]. Acta Tabacaria Sinica, 2019, 25(1): 93-98.

[11] Chen T，Tong H．Higgs boson discovery with boosted trees. In: NIPS Workshop on High-energy Physics and Machine Learning, 2015: 69–80.

[12] 毛開銀, 趙長名, 何嘉. 基于XGBoost的10m風速訂正研究[J]. 成都信息工程大學學報, 2020, 35(6): 604-609.

MAO Kaiyin, ZHAO Changming, HE Jia. A Research for 10 m Wind Speed Prediction based on XGBoost[J]. Journal Of Chengdu University Of Information Technology, 2020, 35(6): 604-609.

[13] L.G. Valiant, A Theory of the Learnable, Communications of the ACM, 1984, 27(11):1134-1142.

[14] ROBERT E. SCHAPIRE. The Strength of Weak Learnability[J]. Machine Learning, 1990, 5: 197-227.

[15] Jerome H. Friedman. Greedy function approximation: A gradient boosting machine[J]. Ann. Statist., 2001, 29(5): 1189-1232. DOI: 10.1214/aos/1013203451

[16] Djalel Benbouzid, Róbert Busa-Fekete, Norman Casagrande, et al. Multiboost: a multi-purpose boosting package. The Journal of Machine Learning Research, 2012, 13(1): 549- 553.

Forecast of industrial-commercial transaction of innovative cigarette products

HAN Weimin1, YANG Yingguang2*

1 Trading Department, China Tobacco E-commerce and Logistics Co. Ltd., Beijing 100055, China;2 Second business service center, Beijing China Tobacco Information Technology Co. Ltd., Beijing 100055, China

Forecasting industrial-commercial transactions of innovative cigarette products is beneficial for the regulation of the micro-level market strategies of industrial and commercial enterprises, thus increasing the sales of characteristic cigarette. The indexes influencing industrial sales were selected and then their correlation matrix was calculated by using the gray theory, based on which the indexes were divided into five categories. The indexes in each category with the largest correlation with industrial sales volume were selected to form an index set, and XGBoost was used to construct a model and make predictions. Taking the industrial sales data of innovative cigarette products in the past ten years as the training set and test set, the prediction model based on XGBoost was trained and tested in comparison with the time-series model, the multiple regression model and the random forest model. Experiments showed that the prediction accuracy of the four models were all above 90%, where XGBoost-based model was the best one, achieving an accuracy of 96%. Therefore, XGBoost can be adopted as a reference algorithm for cigarette sales forecast.

innovative cigarette products; industrial sales volume; industrial-commercial transaction; sales forecast; XGBoost ensemble learning

Corresponding author. Email：yangyg@ctitc.cn

韓偉民（1983—），碩士研究生，高級經(jīng)濟師，行業(yè)高級經(jīng)濟師評審專家?guī)斐蓡T，主要研究方向：電子商務與物流管理，Tel：010-63606847，Email：hanweimin@tobacco.gov.cn。

楊應廣（1974—），碩士研究生，行業(yè)物流分標委委員，主要研究方向：數(shù)據(jù)分析咨詢，Tel：010-63606134，Email：yangyg@ctitc.cn。

2021-05-10；

2021-12-28

韓偉民，楊應廣. 卷煙創(chuàng)新產(chǎn)品工商交易預測研究[J]. 中國煙草學報，2022，28（1）. HAN Weimin, YANG Yingguang. Forecast of industrial-commercial transaction of innovative cigarette products[J]. Acta Tabacaria Sinica, 2022, 28 (1). doi: 10.16472/j.chinatobacco. 2021.088