基于Python与Stata的股利模型搭建

近些年，股利研究是会计研究的焦点之一，2006年的股权分置改革也对上市公司股利发放产生了重大影响，本文在文献《股权分置改革与上市公司股利政策》（链接如下）的基础上研究了2016年-2020年的数据，并用代码与统计软件实现出来。

​​​​​​​股权分置改革与上市公司股利政策_基于迎合理论的证据_支晓强 - 道客巴巴

（一）逻辑梳理

文章的逻辑较为简单：

1. 股利支付意愿(PTP) = 实际发放股利概率 - 预测发放股利概率

用上市公司的历史指标数据作为自变量，是否支付股利作为因变量；逻辑回归后对后续年份的股利支付概率做出预测。如果后续年份实际发放的股利 > 预测发放的股利，则有“股利支付意愿(PTP) > 0”，即为“股利超预期发放”；若后续年份实际发放的股利 < 预测发放的股利，则有“股利支付意愿(PTP) < 0”，即为“股利发放不及预期”。

 

 2. 股利溢价(PDND) = 发股利公司的Avg(ln(市值/净资产)) - 不发股利公司的Avg(ln(市值/净资产))

容易理解，若PDND > 0，说明发股利的公司估值更高，不发股利的公司估值更低；若PDND < 0，说明不发股利的公司估值更高，发股利的公司估值更低。 

 

 3. 迎合理论

为了探寻PTP和PDND在股改前后的相关性变化，文章引入了迎合理论：

（1）对PTP和上年PDND做线性回归，之后对回归参数做显著性检验，来判定二者的相关性；

（2）“基于时序和截面混合数据的检验”，说白了就是把 “PDND（主要）” 和各种乱七八糟的参数（如：LNA, EA... 还有时间Timetrend）放进自变量池，统统做回归看显著性。

 

（二）数据获取与清洗

要获取的数据主要就是文章给出的8个参数，以及计算PDND时要用的市值和账面价值：

首先，文章对个股样本提出了要求：

 我们从CSMAR中获取所有的上市公司代码，CSMAR可以完美的帮我们完成(2)(3)条件的剔除：

 条件(1)没有现成的筛选工具，这里手写了一个api接上东财的ajax：

def examine_bh_stock(code):  # 检验是否存在B股或H股

    query_url = f'http://emweb.securities.eastmoney.com/PC_HSF10/CompanySurvey/CompanySurveyAjax?code={code}'
    # print(query_url)

    ua = 'Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/99.0.4844.51 Safari/537.36'
    referer = f'http://emweb.securities.eastmoney.com/PC_HSF10/CompanySurvey/Index?type=web&code={code}'
    host = 'emweb.securities.eastmoney.com'
    cookie = ''

    headers = {'User-Agent': ua, 'Host': host, 'Referer': referer, }

    response = requests.get(query_url, headers=headers).text
    # print(response)

    if re.findall(r'"bgdm":"--"', response) and re.findall(r'"hgdm":"--"', response):
        return False  # 无B或H股
    else:
        return True  # 有B或H股

note：api建议做东财，目前运行下来几乎没有访问限制；某顺有JS混淆加密，破解的话要copy一份JS做exec，很费时间，不推荐。

条件(4)(5) 的筛选 和 数据指标获取 使用了python的tushare包，这个包功能强大，二级市场上涉及到的指标应有尽有（唯一的美中不足的可能就是2015年之前的部分财务数据还没有收录）。其中托宾Q用市净率代替，H10用前十大流通股东持股比例计算：

def get_financial_data(code, yr):
    pro = ts.pro_api(token='')

    trade_date = str(yr) + '1231'
    if yr == 2017:
        trade_date_adjusted = '20171229'
    elif yr == 2016 or yr == 2011:
        trade_date_adjusted = str(yr) + '1230'
    elif yr == 2018:
        trade_date_adjusted = '20181228'
    else:
        trade_date_adjusted = trade_date

    df0 = pro.fina_indicator(ts_code=code, period=trade_date)
    df1 = pro.bak_basic(ts_code=code, trade_date=trade_date_adjusted)
    daily_df = pro.daily_basic(ts_code=code, trade_date=trade_date_adjusted)
    holders_df = pro.top10_holders(ts_code=code, period=trade_date)
    income_df = pro.income(ts_code=code, period=trade_date)
    dividend_df = pro.dividend(ts_code=code, fields='ts_code,div_proc,cash_div_tax,stk_div,record_date,ex_date')
    # print(dividend_df)
    try:
        div_date = np.array(dividend_df['record_date'])
        div_date = div_date[div_date != None]
        div_arr = np.unique(np.array([ele for ele in div_date if ele.startswith(str(yr + 1))]))
        if div_arr:
            dividend_data = dividend_df[dividend_df['record_date'] == div_arr[0]]
            dividend_tup = (dividend_data['cash_div_tax'].values[0], dividend_data['stk_div'].values[0])
        else:
            dividend_tup = (0, 0)

        stock_price = daily_df['close'].values[0]
        pb = daily_df['pb'].values[0]

        bps = df0['bps'].values[0]
        total_share = daily_df['total_share'].values[0]
        M = (daily_df['circ_mv'].values[0] + (total_share - daily_df['float_share'].values[0]) * bps) * 10000
        B = bps * total_share * 10000

        ni = income_df['n_income'].values[0]

        return df0, df1, stock_price, holders_df, bps, pb, dividend_tup, M, B, ni

    except Exception:
        return ['x'] * 10

需要注意的是，部分年份的最终交易日不是12月31日，需要调整（或者也可以用其他包里的方法自动获取）

（三）数据处理

首先计算PTP，如文章所说，“采用逐年回归的方式，估计股利支付决策的公司特征模型”。即将往年的数据集视作训练集，logit回归出参数后对当年进行预测，与实际比较作差得到PTP：

def sigmoid(z):  # 定义sigmoid函数
    z_ravel = z.ravel()
    length = len(z_ravel)
    y = []
    for index in range(length):
        if z_ravel[index] >= 0:
            y.append(1.0 / (1 + np.exp(-z_ravel[index])))
        else:
            y.append(np.exp(z_ravel[index]) / (np.exp(z_ravel[index]) + 1))
    return np.array(y).reshape(z.shape)


def initialize_with_zeros(dim):
    w = np.zeros((dim, 1))
    b = 0
    return w, b


def propagate(w, b, X, Y):
    # 获取样本数m：
    m = X.shape[1]

    A = sigmoid(np.dot(w.T, X) + b)    # 调用sigmoid函数
    cost = -(np.sum(Y * np.log(A) + (1-Y) * np.log(1-A))) / m

    dZ = A-Y
    dw = (np.dot(X, dZ.T)) / m
    db = (np.sum(dZ)) / m

    grads = {"dw": dw,
             "db": db}

    return grads, cost


def optimize(w, b, X, Y, num_iterations, learning_rate, print_cost=False):
    # costs数组用于存放每若干次迭代后的cost便于后续看趋势
    costs = []
    # 迭代
    for i in range(num_iterations):
        # propagate计算每次迭代后的cost和梯度
        grads, cost = propagate(w, b, X, Y)
        dw = grads["dw"]
        db = grads["db"]

        # 用得到的梯度更新参数
        w = w - learning_rate*dw
        b = b - learning_rate*db

        if i % 100 == 0:
            costs.append(cost)

        # 实时打印
        if print_cost and i % 100 == 0:
            print(f"Cost after iteration {i}: {cost}")

    params = {"w": w,
              "b": b}
    grads = {"dw": dw,
             "db": db}
    return params, grads, costs


def to_verify(X_train, Y_train, W, b):
    m = X_train.shape[1]
    Y_prediction = np.zeros((1, m))

    A = sigmoid(np.dot(W.T, X_train) + b)
    for i in range(m):
        if A[0, i] > 0.5:
            Y_prediction[0, i] = 1
        else:
            Y_prediction[0, i] = 0

    err_arr = Y_prediction[0] - Y_train[0]
    accurate_rate = np.sum(err_arr == 0) / len(err_arr)

    return Y_prediction, accurate_rate

def training_package(yr_list, indicator):
    train_set_cash_div, train_set_stk_div = extract_training_sets(yr_list)
    # print(train_set_cash_div)
    if indicator == 'stk':
        train_set = train_set_stk_div
    else:
        train_set = train_set_cash_div

    X_train, Y_train = [], []
    for i in range(9):
        temp_arr = []
        for j in train_set:
            temp_arr.append(j[i])
        if i != 8:
            X_train.append(temp_arr)
        else:
            Y_train.append(temp_arr)

    X_train = np.array(X_train)
    Y_train = np.array(Y_train)

    # Logit回归
    # print(X_train, Y_train)
    print('Logit数据集开始训练...')
    dim = X_train.shape[0]
    # dim_Y = Y_train.shape[0]
    W, b = initialize_with_zeros(dim)

    num_iterations = 2000
    learning_rate = 0.1
    print_cost = False
    params, grads, costs = optimize(W, b, X_train, Y_train, num_iterations, learning_rate, print_cost)
    W = params['w']
    b = params['b']
    print(W, b)

    # 检验样本准确性
    self_prediction_arr, accurate_rate = to_verify(X_train, Y_train, W, b)
    print(self_prediction_arr, '\r\n', f'准确率：{accurate_rate * 100}%')

    return W, b, self_prediction_arr, accurate_rate

当然，也可以使用统计软件如stata；stata的回测准确性更高：

做完回归后直接matplot绘图，或许是训练集太少的缘故，2017-2020数据的正相关性是惊人的不显著：（甚至感觉是负相关？离大谱...）

 哈哈，这么烂的数据，检验也不用做啦，直接宣告失败吧 :-D

算了还是做一下

 现金股利支付意愿(cash_PTP)与上一年的投资者现金股利偏好(cash_PDND t-1)回归结果：

 

 股票股利支付意愿(stk_PTP)与上一年的投资者股票股利偏好(stk_PDND t-1)回归结果：

 

意料之中，结果令所有人失望。现金股利的回归系数不显著，股票股利的回归系数显著了，但是是负相关。

至于文章后面的时间截面Tobit检验，这里由于只有4年的样本，多重共线性严重，就不做了：

 （四）结论

很遗憾的，论文给出了4个主结论到这里只剩1个存活了下来，它就是：

 看来我国的投资者真的很厌恶现金股利...

当然，受制于训练集，本文做出的研究本身就存在较大的偏误（尤其是PTP的估计），从而可能导致结论与原文的截然不同；本人小白一枚，如果做的不好还请各位读者海涵，虚心接受建议和批评，谢谢！