数字货币因子模型

Author: 小草, Created: 2022-09-27 16:10:28, Updated: 2023-09-15 20:58:26

[TOC]

img

因子模型框架

股市的多因子模型的研究报告可谓汗牛充栋,有着丰富的理论和实践。数字货币市场无论币种数量、总市值、交易额、衍生品市场等都足以满足进行因子研究,本文主要面对量化策略初学者,不会涉及复杂的数学原理和统计分析,将以币安永续期货市场为数据源,构建一个简单的因子研究的框架,方便对因子指标进行评价。

因子可以看作一个指标,可以写作表达式,因子不断变化,反映了未来的收益信息,通常因子代表了一种投资逻辑。

例如:收盘价close这个因子,背后的假设就是股价可以预测未来收益,股价越高未来收益越高(也可能越低),以此因子构建组合其实就是定期轮仓购买高价股的投资模式/策略。 通常而言那些能够持续产生超额收益的因子又常被称为Alpha。例如市值因子、动量因子等都被学术界和投资界验证过是曾经有效的因子。

无论是股市还是数字货币市场,都是一个复杂的系统,没有因子可以完全预测未来的收益,但仍然具有一定的可预测性。有效的alpha(投资模式)并随着更多资金的投入而逐渐失效。但这一过程将在市场上会产生其它的模式,从而诞生新的alpha。市值因子曾经在A股市场上是一个非常有效的策略 ,简单买入10个的最低市值的股票,每天调整一次,从2007年开始十年回测将获得超过400倍的收益,远远超越大盘。但2017年的白马股行情反映了就是小市值因子的失效,价值因子反而流行了起来。因此需要不断在对alpha的验证和使用之间权衡和尝试。

寻找的因子是建立策略的基础,通过组合多个不相关的有效因子可以构建比较好的策略。

import requests
from datetime import date,datetime
import time
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import requests, zipfile, io
%matplotlib inline

数据来源

币安USDT永续期货2022年初至今的小时K线数据,截止到目前,超过了150个币种。前面说过,因子模型是一种选币模型,是面向全部币种而不是某个币种。K线数据中包含高开低收的价格、成交量、成交笔数、主动买量等数据,这些数据当然不是所有因子的来源,比如美股指数、加息预期、盈利能力、链上数据、社交媒体热度等等。冷门的数据来源也可能发掘出有效的alpha,然而基础的量价数据也完全够用。

## 当前交易对
Info = requests.get('https://fapi.binance.com/fapi/v1/exchangeInfo')
symbols = [s['symbol'] for s in Info.json()['symbols']]
symbols = list(filter(lambda x: x[-4:] == 'USDT', [s.split('_')[0] for s in symbols]))
print(symbols)

Out:

['BTCUSDT', 'ETHUSDT', 'BCHUSDT', 'XRPUSDT', 'EOSUSDT', 'LTCUSDT', 'TRXUSDT', 'ETCUSDT', 'LINKUSDT',
'XLMUSDT', 'ADAUSDT', 'XMRUSDT', 'DASHUSDT', 'ZECUSDT', 'XTZUSDT', 'BNBUSDT', 'ATOMUSDT', 'ONTUSDT',
'IOTAUSDT', 'BATUSDT', 'VETUSDT', 'NEOUSDT', 'QTUMUSDT', 'IOSTUSDT', 'THETAUSDT', 'ALGOUSDT', 'ZILUSDT',
'KNCUSDT', 'ZRXUSDT', 'COMPUSDT', 'OMGUSDT', 'DOGEUSDT', 'SXPUSDT', 'KAVAUSDT', 'BANDUSDT', 'RLCUSDT',
'WAVESUSDT', 'MKRUSDT', 'SNXUSDT', 'DOTUSDT', 'DEFIUSDT', 'YFIUSDT', 'BALUSDT', 'CRVUSDT', 'TRBUSDT',
'RUNEUSDT', 'SUSHIUSDT', 'SRMUSDT', 'EGLDUSDT', 'SOLUSDT', 'ICXUSDT', 'STORJUSDT', 'BLZUSDT', 'UNIUSDT',
'AVAXUSDT', 'FTMUSDT', 'HNTUSDT', 'ENJUSDT', 'FLMUSDT', 'TOMOUSDT', 'RENUSDT', 'KSMUSDT', 'NEARUSDT',
'AAVEUSDT', 'FILUSDT', 'RSRUSDT', 'LRCUSDT', 'MATICUSDT', 'OCEANUSDT', 'CVCUSDT', 'BELUSDT', 'CTKUSDT',
'AXSUSDT', 'ALPHAUSDT', 'ZENUSDT', 'SKLUSDT', 'GRTUSDT', '1INCHUSDT', 'CHZUSDT', 'SANDUSDT', 'ANKRUSDT',
'BTSUSDT', 'LITUSDT', 'UNFIUSDT', 'REEFUSDT', 'RVNUSDT', 'SFPUSDT', 'XEMUSDT', 'BTCSTUSDT', 'COTIUSDT',
'CHRUSDT', 'MANAUSDT', 'ALICEUSDT', 'HBARUSDT', 'ONEUSDT', 'LINAUSDT', 'STMXUSDT', 'DENTUSDT', 'CELRUSDT',
'HOTUSDT', 'MTLUSDT', 'OGNUSDT', 'NKNUSDT', 'SCUSDT', 'DGBUSDT', '1000SHIBUSDT', 'ICPUSDT', 'BAKEUSDT',
'GTCUSDT', 'BTCDOMUSDT', 'TLMUSDT', 'IOTXUSDT', 'AUDIOUSDT', 'RAYUSDT', 'C98USDT', 'MASKUSDT', 'ATAUSDT',
'DYDXUSDT', '1000XECUSDT', 'GALAUSDT', 'CELOUSDT', 'ARUSDT', 'KLAYUSDT', 'ARPAUSDT', 'CTSIUSDT', 'LPTUSDT',
'ENSUSDT', 'PEOPLEUSDT', 'ANTUSDT', 'ROSEUSDT', 'DUSKUSDT', 'FLOWUSDT', 'IMXUSDT', 'API3USDT', 'GMTUSDT',
'APEUSDT', 'BNXUSDT', 'WOOUSDT', 'FTTUSDT', 'JASMYUSDT', 'DARUSDT', 'GALUSDT', 'OPUSDT', 'BTCUSDT',
'ETHUSDT', 'INJUSDT', 'STGUSDT', 'FOOTBALLUSDT', 'SPELLUSDT', '1000LUNCUSDT', 'LUNA2USDT', 'LDOUSDT',
'CVXUSDT']

print(len(symbols))

Out:

153

#获取任意周期K线的函数
def GetKlines(symbol='BTCUSDT',start='2020-8-10',end='2021-8-10',period='1h',base='fapi',v = 'v1'):
    Klines = []
    start_time = int(time.mktime(datetime.strptime(start, "%Y-%m-%d").timetuple()))*1000 + 8*60*60*1000
    end_time =  min(int(time.mktime(datetime.strptime(end, "%Y-%m-%d").timetuple()))*1000 + 8*60*60*1000,time.time()*1000)
    intervel_map = {'m':60*1000,'h':60*60*1000,'d':24*60*60*1000}
    while start_time < end_time:
        mid_time = start_time+1000*int(period[:-1])*intervel_map[period[-1]]
        url = 'https://'+base+'.binance.com/'+base+'/'+v+'/klines?symbol=%s&interval=%s&startTime=%s&endTime=%s&limit=1000'%(symbol,period,start_time,mid_time)
        res = requests.get(url)
        res_list = res.json()
        if type(res_list) == list and len(res_list) > 0:
            start_time = res_list[-1][0]+int(period[:-1])*intervel_map[period[-1]]
            Klines += res_list
        if type(res_list) == list and len(res_list) == 0:
            start_time = start_time+1000*int(period[:-1])*intervel_map[period[-1]]
        if mid_time >= end_time:
            break
    df = pd.DataFrame(Klines,columns=['time','open','high','low','close','amount','end_time','volume','count','buy_amount','buy_volume','null']).astype('float')
    df.index = pd.to_datetime(df.time,unit='ms')
    return df
start_date = '2022-1-1'
end_date = '2022-09-14'
period = '1h'
df_dict = {}
for symbol in symbols:
    df_s = GetKlines(symbol=symbol,start=start_date,end=end_date,period=period,base='fapi',v='v1')
    if not df_s.empty:
        df_dict[symbol] = df_s
symbols = list(df_dict.keys())
print(df_s.columns)

Out:

Index(['time', 'open', 'high', 'low', 'close', 'amount', 'end_time', 'volume',
       'count', 'buy_amount', 'buy_volume', 'null'],
      dtype='object')

初步从K线数据中提炼出我们感兴趣的数据:收盘价、开盘价、成交量、成交笔数、主动买入比例,以这些数据为基础,加工出需要的因子。

df_close = pd.DataFrame(index=pd.date_range(start=start_date, end=end_date, freq=period),columns=df_dict.keys())
df_open = pd.DataFrame(index=pd.date_range(start=start_date, end=end_date, freq=period),columns=df_dict.keys())
df_volume = pd.DataFrame(index=pd.date_range(start=start_date, end=end_date, freq=period),columns=df_dict.keys())
df_buy_ratio = pd.DataFrame(index=pd.date_range(start=start_date, end=end_date, freq=period),columns=df_dict.keys())
df_count = pd.DataFrame(index=pd.date_range(start=start_date, end=end_date, freq=period),columns=df_dict.keys())
for symbol in df_dict.keys():
    df_s = df_dict[symbol]
    df_close[symbol] = df_s.close
    df_open[symbol] = df_s.open
    df_volume[symbol] = df_s.volume
    df_count[symbol] = df_s['count']
    df_buy_ratio[symbol] = df_s.buy_amount/df_s.amount
df_close = df_close.dropna(how='all')
df_open = df_open.dropna(how='all')
df_volume = df_volume.dropna(how='all')
df_count = df_count.dropna(how='all')
df_buy_ratio = df_buy_ratio.dropna(how='all')

总览下市场指数的表现,可以说比较惨淡,年初至今下跌了60%。

df_norm = df_close/df_close.fillna(method='bfill').iloc[0] #归一化
df_norm.mean(axis=1).plot(figsize=(15,6),grid=True);
#最终指数收益图

img

因子有效性判定

  • 回归法 以下一期的收益率作为因变量,待测试的因子作为自变量,回归得到的系数也就是因子的收益率。构造回归方程之后,一般会参考系数t值的绝对值均值,系数t值绝对值序列大于2的占比,年化因子收益率,年化因子收益波动率,因子收益的夏普比等参数来看因子的有效性和波动性。可以一次回归多个因子,具体参考barra文档。

  • IC,IR等指标 所谓IC就是因子与下一期收益率的相关系数,现在一般也用RANK_IC,就是因子排名与下一期股票收益率的相关系数。IR一般来讲是IC序列的均值/IC序列的标准差。

  • 分层回归法 本文将用这种方法,就是根据待测试的因子排序,将币种分为N组来进行分组回测,使用固定的周期来进行调仓的操作。 如果情况理想,N 组币种的收益率会呈现较好的单调性,单调递增或递减,且每一组的收益差距较大。这样的因子体现为较好的区分度。假如第一组收益最高,最后一组收益最低,那么做多第一组组合并做空最后一组组合,最终得到的收益率,夏普比率的参考指标。

实际回测操作

根据因子从小到大排序把待选币种根据排序分为3组,每组币种大约占1/3,如果一个因子有效的话,每组分的越少往往收益率越高,但也意味着每个币种分配的资金相对较多,如果多空分别一倍杠杆,第一组和最后一组分别为10个币种,则一个比重占比10%,如果做空的某个币种上涨了2倍,则回撤20%;相应的如果分组数量为50,则回撤4%。分散币种可以降低黑天鹅的风险。做多第一组(因子值最小),做空第三组。如果因子越大收益越高,可以把多空反过来或者简单将因子变为负数或者倒数。

通常可以根据最终回测的收益率和夏普比例来粗略评估因子预测能力。此外还需要参考因子表达式是否简单、对分组的大小不敏感、对调仓间隔不敏感、对回测初始时间不敏感等。

关于调仓频率,股票市场往往是5天、10天以及一个月为周期,但对于数字货币市场,这样的周期无疑是太长了,并且实盘中行情是实时监控的,死守一个特定的周期再调仓没有必要,因此在实盘中我们是实时或者短时间周期调仓。

关于如何平仓,按照传统的方法,下次排序时不在组别中即可平仓。但是在实时调仓的情况下,一些币种可能正好处于分界处,会出现来回平仓的情况。因此本策略采用等待分组变化,需要开反方向仓位时再平仓,如第一组做多,当处于做多状态的币种被分为第三组后,此时再平仓做空。如果固定周期平仓,比如每天或每8小时,也可以采用不处于分组中即平仓的方式。可以多进行尝试。

#回测引擎
class Exchange:
    
    def __init__(self, trade_symbols, fee=0.0004, initial_balance=10000):
        self.initial_balance = initial_balance #初始的资产
        self.fee = fee
        self.trade_symbols = trade_symbols
        self.account = {'USDT':{'realised_profit':0, 'unrealised_profit':0, 'total':initial_balance, 'fee':0, 'leverage':0, 'hold':0}}
        for symbol in trade_symbols:
            self.account[symbol] = {'amount':0, 'hold_price':0, 'value':0, 'price':0, 'realised_profit':0,'unrealised_profit':0,'fee':0}
            
    def Trade(self, symbol, direction, price, amount):
        
        cover_amount = 0 if direction*self.account[symbol]['amount'] >=0 else min(abs(self.account[symbol]['amount']), amount)
        open_amount = amount - cover_amount
        self.account['USDT']['realised_profit'] -= price*amount*self.fee #扣除手续费
        self.account['USDT']['fee'] += price*amount*self.fee
        self.account[symbol]['fee'] += price*amount*self.fee

        if cover_amount > 0: #先平仓
            self.account['USDT']['realised_profit'] += -direction*(price - self.account[symbol]['hold_price'])*cover_amount  #利润
            self.account[symbol]['realised_profit'] += -direction*(price - self.account[symbol]['hold_price'])*cover_amount
            
            self.account[symbol]['amount'] -= -direction*cover_amount
            self.account[symbol]['hold_price'] = 0 if self.account[symbol]['amount'] == 0 else self.account[symbol]['hold_price']
            
        if open_amount > 0:
            total_cost = self.account[symbol]['hold_price']*direction*self.account[symbol]['amount'] + price*open_amount
            total_amount = direction*self.account[symbol]['amount']+open_amount
            
            self.account[symbol]['hold_price'] = total_cost/total_amount
            self.account[symbol]['amount'] += direction*open_amount
                    
    
    def Buy(self, symbol, price, amount):
        self.Trade(symbol, 1, price, amount)
        
    def Sell(self, symbol, price, amount):
        self.Trade(symbol, -1, price, amount)
        
    def Update(self, close_price): #对资产进行更新
        self.account['USDT']['unrealised_profit'] = 0
        self.account['USDT']['hold'] = 0
        for symbol in self.trade_symbols:
            if not np.isnan(close_price[symbol]):
                self.account[symbol]['unrealised_profit'] = (close_price[symbol] - self.account[symbol]['hold_price'])*self.account[symbol]['amount']
                self.account[symbol]['price'] = close_price[symbol]
                self.account[symbol]['value'] = abs(self.account[symbol]['amount'])*close_price[symbol]
                self.account['USDT']['hold'] += self.account[symbol]['value']
                self.account['USDT']['unrealised_profit'] += self.account[symbol]['unrealised_profit']
        self.account['USDT']['total'] = round(self.account['USDT']['realised_profit'] + self.initial_balance + self.account['USDT']['unrealised_profit'],6)
        self.account['USDT']['leverage'] = round(self.account['USDT']['hold']/self.account['USDT']['total'],3)

#测试因子的函数
def Test(factor, symbols, period=1, N=40, value=300):
    e = Exchange(symbols, fee=0.0002, initial_balance=10000)
    res_list = []
    index_list = []
    factor = factor.dropna(how='all')
    for idx, row in factor.iterrows():
        if idx.hour % period == 0:
            buy_symbols =  row.sort_values().dropna()[0:N].index
            sell_symbols = row.sort_values().dropna()[-N:].index
            prices = df_close.loc[idx,]
            index_list.append(idx)
            for symbol in symbols:
                if symbol in buy_symbols and e.account[symbol]['amount'] <= 0:
                    e.Buy(symbol,prices[symbol],value/prices[symbol]-e.account[symbol]['amount'])
                if symbol in sell_symbols and e.account[symbol]['amount'] >= 0:
                    e.Sell(symbol,prices[symbol], value/prices[symbol]+e.account[symbol]['amount'])
            e.Update(prices)
            res_list.append([e.account['USDT']['total'],e.account['USDT']['hold']])
    return pd.DataFrame(data=res_list, columns=['total','hold'],index = index_list)

简单的因子测试

成交量因子:简单做多成交量低的币种,做空成交量高的币种,表现非常好,这说明热门币更倾向于下跌。

成交价因子:做多价格低币种,做空价格高的币种,效果一般。

成交笔数因子:表现和成交量非常相似,。可以很明显的注意到成交量因子和成交笔数因子的相关性非常高,实际也是如此,它们的不同币种的平均相关性达到了0.97,这说明这两个因子非常相近,在合成多因子时,这个因素需要考虑在内。

3h动量因子:(df_close - df_close.shift(3))/df_close.shift(3)。即因子的3小时涨幅,回测结果可见3小时涨幅具有明显的回归特性,即上涨的在接下来更容易下跌。整体表现可以,但也有较长时间的回撤和振荡期。

24h动量因子:24h的调仓周期结果还不错,收益与3h动量相近,回撤更小。

成交额变动因子:df_volume.rolling(24).mean()/df_volume.rolling(96).mean(),即最近1天成交额与最近3天成交额的比值,每8h调仓一次。回测结果表现比较好,回撤也比较低,这说明成交量活跃的反而更倾向于下跌。

成交笔数变动因子:df_count.rolling(24).mean()/df_count.rolling(96).mean(),即最近1天成交笔数与最近3天成交笔数的比值,每8h调仓一次。回测结果表现比较好,回撤也比较低,这说明成交笔数增加活跃的反而更倾向于下跌。

单笔成交价值变动因子: -(df_volume.rolling(24).mean()/df_count.rolling(24).mean())/(df_volume.rolling(24).mean()/df_count.rolling(96).mean()) ,即最近1天成交价值与最近3天成交价值的比值,每8h调仓一次。这个因子也和成交量因子高度相关。

主动成交比例变动因子:df_buy_ratio.rolling(24).mean()/df_buy_ratio.rolling(96).mean(),即最近1天主动买入量与总成交量的比值与最近3天成交价值的比值,每8h调仓一次。这个因子表现尚可,和成交量因子相关性不大。

波动率因子:(df_close/df_open).rolling(24).std(), 做多波动性小的币种,有一定效果。

成交量和收盘价相关性因子:df_close.rolling(96).corr(df_volume), 最近4天收盘价有成交量的相关性因子,整体表现不错。

这里列出的只是一些基于量价的基础的因子,实际上因子公式组合可以非常复杂,可以没有明显的逻辑。可以参考著名的ALPHA101的因子构造方式:https://github.com/STHSF/alpha101

#成交量
factor_volume = df_volume
factor_volume_res = Test(factor_volume, symbols, period=4)
factor_volume_res.total.plot(figsize=(15,6),grid=True);

img

#成交价
factor_close = df_close
factor_close_res = Test(factor_close, symbols, period=8)
factor_close_res.total.plot(figsize=(15,6),grid=True);

img

#成交笔数
factor_count = df_count
factor_count_res = Test(factor_count, symbols, period=8)
factor_count_res.total.plot(figsize=(15,6),grid=True);

img

print(df_count.corrwith(df_volume).mean())

0.9671246744996017

#3小时动量因子
factor_1 =  (df_close - df_close.shift(3))/df_close.shift(3)
factor_1_res = Test(factor_1,symbols,period=1)
factor_1_res.total.plot(figsize=(15,6),grid=True);

img

#24小时动量因子
factor_2 =  (df_close - df_close.shift(24))/df_close.shift(24)
factor_2_res = Test(factor_2,symbols,period=24)
tamenxuanfactor_2_res.total.plot(figsize=(15,6),grid=True);

img

#成交量因子
factor_3 = df_volume.rolling(24).mean()/df_volume.rolling(96).mean()
factor_3_res = Test(factor_3, symbols, period=8)
factor_3_res.total.plot(figsize=(15,6),grid=True);

img

#成交笔数因子
factor_4 = df_count.rolling(24).mean()/df_count.rolling(96).mean()
factor_4_res = Test(factor_4, symbols, period=8)
factor_4_res.total.plot(figsize=(15,6),grid=True);

img

#因子相关性
print(factor_4.corrwith(factor_3).mean())

0.9707239580854841

#单笔成交价值因子
factor_5 = -(df_volume.rolling(24).mean()/df_count.rolling(24).mean())/(df_volume.rolling(24).mean()/df_count.rolling(96).mean())
factor_5_res = Test(factor_5, symbols, period=8)
factor_5_res.total.plot(figsize=(15,6),grid=True);

img

print(factor_4.corrwith(factor_5).mean())

0.861206620552479

#主动成交比例因子
factor_6 = df_buy_ratio.rolling(24).mean()/df_buy_ratio.rolling(96).mean()
factor_6_res = Test(factor_6, symbols, period=4)
factor_6_res.total.plot(figsize=(15,6),grid=True);

img

print(factor_3.corrwith(factor_6).mean())

0.1534572192503726

#波动率因子
factor_7 = (df_close/df_open).rolling(24).std()
factor_7_res = Test(factor_7, symbols, period=2)
factor_7_res.total.plot(figsize=(15,6),grid=True);

img

#成交量和收盘价相关性因子
factor_8 = df_close.rolling(96).corr(df_volume)
factor_8_res = Test(factor_8, symbols, period=4)
factor_8_res.total.plot(figsize=(15,6),grid=True);

img

多因子合成

不断挖掘出新的有效因子固然是构建策略过程中最为重要的一部分,但如果没有好的因子合成方法,优秀单个Alpha因子也不能发挥其最大作用。常见的多因子合成方法有:

等权法: 所有待合成因子等权重相加,得到新的合成后因子。

历史因子收益率加权法:所有待合成因子,按照最近一段时期内历史因子收益率的算术平均值作为权重进行相加,得到新的合成后因子。这种方法表现好的因子权重更高。

最大化IC_IR加权法:以历史一段时间的复合因子平均 IC值作为对复合因子下一期IC值的估计,以历史IC值的协方差矩阵作为对复合因下一期波动率的估计,根据IC_IR等于IC的期望值除以IC的标准差,可以得到最大化复合因子IC_IR的最优权重解。

主成分分析(PCA)法 :PCA 是数据降维的常用方法,因子之间的相关性可能比较高,使用降维后的主成分作为合成后的因子。

本文将手动参考因子有效性赋权。上面介绍的方法可参考:ae933a8c-5a94-4d92-8f33-d92b70c36119.pdf

在测试单因子的时候排序是固定的,但是多因子合成需要把完全不同的数据合并在一起,因此需要将所有的因子进行标准化处理,一般还需要去极值和缺失值处理。这里我们采用df_volume\factor_1\factor_7\factor_6\factor_8合成。

#标准化函数,去除缺失值和极值,并且进行标准化处理
def norm_factor(factor):
    factor = factor.dropna(how='all')
    factor_clip = factor.apply(lambda x:x.clip(x.quantile(0.2), x.quantile(0.8)),axis=1)
    factor_norm = factor_clip.add(-factor_clip.mean(axis=1),axis ='index').div(factor_clip.std(axis=1),axis ='index')
    return factor_norm


df_volume_norm = norm_factor(df_volume)
factor_1_norm = norm_factor(factor_1)
factor_6_norm = norm_factor(factor_6)
factor_7_norm = norm_factor(factor_7)
factor_8_norm = norm_factor(factor_8)
factor_total = 0.6*df_volume_norm + 0.4*factor_1_norm + 0.2*factor_6_norm + 0.3*factor_7_norm + 0.4*factor_8_norm
factor_total_res = Test(factor_total, symbols, period=8)
factor_total_res.total.plot(figsize=(15,6),grid=True);

img

总结

本文介绍了单因子的测试方法并测试了常见的单因子,初步介绍了多因子合成的方法,但多因子研究内容非常丰富,文中提到每一个点都能深入的展开,把对各种各样的策略研究转为对alpha因子的发掘是个可行的途径,使用因子的方法论,能极大的加快交易思想的验证,并且有很多可供参考的资料。

实盘地址:https://www.fmz.com/robot/486605


Related

More

Chanking 写得真好

去者伯仁 草神威武!!!最近刚好也在研究这个

轻轻的云 完了。。。。彻底看不懂。。。。。。

cjz140 草神威武!!!

jmxjqr0302 草神威武!!!

jmxjqr0302 草神威武!!!

jmxjqr0302 草神威武!!!

f_q 草神威武!!!

摩诃无量 草神威武!!!

Tututu001 草神威武!!!

xunfeng91 草神威武!!!