时间序列在资金计划预测中的应用

概述

公司历年预算与实际使用资金总是存在出入,做了预算但是没有被使用完的资金量高居不下,造成了不小的闲置成本,为了提高资金计划的准确性,减少资金预算与实际使用的差距,降低资金使用成本,并且在这个过程中探讨使用数据的手段对未来趋势进行预测的可行性,我们收集了最近几年的实际资金使用情况,试图通过对这些历史资金使用行为的分析,找出数据的规律,并利用这个规律对未来的资金计划做出预测,希望对后面的预算给出一些参考。

准备数据

我们收集了最近五年所有资金使用行为数据,这些数据包括了12万左右的合同付款记录,真实反映了资金计划的实际使用情况。

import os

work_dir = 'D:/ml/资金计划预测'
assert os.path.exists(work_dir), '工作目录不存在'

data_dir = os.path.join(work_dir, 'data')
assert os.path.exists(data_dir), '数据目录不存在'

output_dir = os.path.join(work_dir, 'output')
if not os.path.exists(output_dir):
    os.mkdir(output_dir)

os.chdir(work_dir)
import numpy as np
import pandas as pd

data_file_name = '付款列表.csv'
data_file_path = os.path.join(data_dir, data_file_name)

df = pd.read_csv(data_file_path, engine='python')

df.head()
资金日期 大区 金额
0 2013/8/29 22:57 大区1 9.270562e+04
1 2013/9/4 23:58 大区1 0.000000e+00
2 2013/8/29 22:24 大区1 0.000000e+00
3 2013/8/29 22:48 大区1 1.250790e+06
4 2013/9/2 19:34 大区1 3.347703e+04

数据预处理

  • 空值和无效值处理
  • 将数据转换为时间序列
  • 数据金额过大,为提高可读性,转换为以亿元为单位
  • 过滤数据,只考虑最近5年的数据
# 数据预处理
df.rename(columns={'大区':'area', '资金日期':'payment_date', '金额':'amount'}, inplace = True)
df.dropna(subset=['payment_date'], inplace=True)
df = df.loc[:, ['payment_date', 'amount', 'area']]
df_index = df.set_index('payment_date')
df_index.index = pd.DatetimeIndex(df_index.index)

# 转换为以亿元为单位
df_index['amount'] = df_index['amount'] / 100000000

# 只看最近5年的数据
df_index = df_index['2013':'2017']

分析数据

原始数据是以分钟维度的付款记录,这个时间粒度对我们来说太细了,不能说明什么问题,根据实际情况,我们分析和预测的过程以月的时间粒度来处理数据就足够了(当然我们也可以以天或周的维度来处理数据)。

我们把数据以月为维度进行汇总,得到的数据是什么样的,先来一个预览,看看整体曲线是个什么样子。

import matplotlib.pyplot as plt
import matplotlib.pylab as pylab

# 设置画图参数
params = {
        'axes.titlesize': '18',
        'axes.labelsize': '13',
        'xtick.labelsize': '13',
        'ytick.labelsize': '13',
        'lines.linewidth': '2',
        'legend.fontsize': '13',
        'figure.figsize': '12, 4',
#         'figure.facecolor': 'white',
        'figure.facecolor': 'snow',
        # 正常显示中文
        'font.sans-serif': 'SimHei',
        # 正常显示负号
        'axes.unicode_minus': False
    }
pylab.rcParams.update(params)
df_groupby = df_index.groupby(pd.TimeGrouper(freq='M')).sum()
plt.plot(df_groupby, label='付款记录')

plt.title('资金使用记录')
plt.xlabel('月份')
plt.ylabel('付款总额(亿元)')
plt.legend(loc='best')
plt.grid(True)
plt.show()

分析这个图,我们有几个结论

  • 数据以年为周期表现出周期性
  • 每年的12月至次年的1月为付款的高峰期,这符合春节前集中付款的实际情况,另外年中也有一个高峰期
  • 每年大概第一季度末到第二季初,第三季度末到第四季度初都是付款的淡季,体现在两个明显的波谷,而且前一个波谷比第二个波谷还要低
  • 总体来看每个周期的数据在纵轴上都整体往上迁移,这说明付款额每年都有增长,侧面反映出公司业务每年都保持一定的增长率,并且最近一年相对上一年的增长率有大幅提升

如果对上面总结的结论还找不到感觉的话,我们换个角度再来看看。按不同年份画图,看看不同年份的每个月的付款额走势是否具有相似性。

# 统计数据
def groupby_df(df, period='month'):
    assert len(df) != 0, '数据集不能为空'
    
    df.dropna(inplace=True)
    if(period == 'month'):
        # 按月查看数据
        df_groupby = df.groupby(pd.TimeGrouper(freq='M')).sum()
        df_groupby['period'] = df_groupby.index.month
    elif(period == 'week'):
        # 按周查看数据
        df_groupby = df.groupby(pd.TimeGrouper(freq='w')).sum()
        df_groupby['period'] = df_groupby.index.weekofyear
    elif(period == 'date'):
        # 按天查看数据
        df_groupby = df.groupby(pd.TimeGrouper(freq='D')).sum()
        df_groupby['period'] = df_groupby.index.dayofyear
    else:
        # 按月查看数据
        df_groupby = df.groupby(pd.TimeGrouper(freq='M')).sum()
        df_groupby['period'] = df_groupby.index.month

    return df_groupby
# 按大区分组数据
def df_area(area):
    if(area == '全国'):
        return df_index
    else:
        return df_index[(df_index['area'] == area)].copy()

area_names = ['全国', '大区1', '大区2', '大区3', '大区4', '大区5', '大区6']
area_titles = ('资金使用趋势(全国)', '资金使用趋势(大区1)', '资金使用趋势(大区2)', '资金使用趋势(大区3)', '资金使用趋势(大区4)',
          '资金使用趋势(大区5)', '资金使用趋势(大区6)')
df_areas = list(map(df_area, area_names))
def plot_dfs(dfs, titles, period='month'):
    assert len(dfs) != 0, '数据集不能为空'
    assert len(dfs) == len(titles), '数据集长度与标题集长度必须相等'

#     plt.figure(figsize=(12, 7 * len(dfs)))
    for i in range(0, len(dfs)):
        df = dfs[i]
        assert len(df) != 0, '数据集不能为空'
        df.dropna(inplace=True)

        df_groupby = groupby_df(df)

        # 按年汇总
        years = np.unique(df_groupby.index.year)

        plt.subplot(len(dfs), 1, i+1)

        for year in years:
            df_year = df_groupby[str(year)].set_index('period')
            plt.plot(df_year, label=year)

        fig = plt.gcf()
        fig.set_size_inches(12, 4 * len(dfs))
        fig.autofmt_xdate()
        
        plt.subplots_adjust(right=0.99, left=0.125, bottom=0.1, top=0.9)
        plt.title(titles[i])
        plt.xlabel('月份')
        plt.ylabel('付款总额(亿元)')
        plt.grid(True)
        plt.legend(loc='best')
        
    plt.show()
    
# plot_dfs([df_area('全国'), df_area('B16_HDDQ0_华东大区')], ['资金使用趋势(全国)', '资金使用趋势(华东)'])
plot_dfs([df_area('全国')], ['资金使用趋势(全国)'])

全国数据的规律是这个样子,那每个大区是不是这个样子呢?我们把各个大区的数据放在一起来做一个对比。

plot_dfs(df_areas, area_titles)

可以看出各个大区的付款风格就各有不同了,但总体规律还是保持不变。有部分月份缺数据,这是因为各个大区由于上系统的时间不同,较早的月份有些大区没有数据。

综上,我们可以得出一个基本结论,这是一个有规律的时间序列,具备可预测的前提。

但这个结论更多的是人为的判断,接下来我们通过技术手段来验证这个时间序列是否是稳定的。

验证时间序列的稳定性

时间序列可以进行分析和预测的前提,是时间序列必须是稳定的,或者它的差分是稳定的。

时间序列就是按时间顺序排列的,随时间变化的数据序列。生活中各行各业都有时间序列的场景,比如销售额、顾客数、访问量、股价、油价、GDP、气温等等。

那什么样的时间序列是稳定的?它需要具备如下的条件:

  • 稳定的移动平均值
  • 稳定的移动方差

下面开始验证我们要分析的时间序列是否是稳定的时间序列,以12个月为时间窗口计算移动平均值和移动方差,把移动平均值和移动方差都画到坐标系上,如果数据原值的移动平均值和移动方差随着时间剧烈波动,那么说明原始数据还不是稳定序列,这时候可以对原始数据做差分分析。

df_groupby = groupby_df(df_index)
ts = df_groupby['amount']
from statsmodels.tsa.stattools import adfuller
from pandas import Series

def rolling_ts(ts):
    ts.dropna(inplace=True)
    # 计算移动平均值
    rolling_mean = Series.rolling(ts, window=12, center=False).mean()
    rolling_std = Series.rolling(ts, window=12, center=False).std()

    # 画图
    plt.plot(ts, color='blue',label='原值')
    plt.plot(rolling_mean, color='red', label='移动平均值')
    plt.plot(rolling_std, color='black', label = '移动标准差')
    
    plt.grid(True)
    plt.legend(loc='best')
    plt.title('移动平均值和移动标准差')
    plt.show(block=False)
# 查看原始值的稳定性
rolling_ts(ts)

可以看到原始序列的移动平均值曲线和移动方差曲线的前半段并不稳定,后半段还算比较平稳,这是因为早期部分大区没有上系统收集数据,整体的全国曲线受到这部分数据缺失的干扰所致,不过这没有关系,这部分缺失数据的影响后面被整体平滑之后,影响会变得很小。

接下来对序列做一阶差分,再来观察稳定性。

# 查看一阶差分的稳定性
ts_diff = ts - ts.shift()
rolling_ts(ts_diff)

如上图,可见一阶差分之后的时间序列的移动平均值和移动方差随时间波动很小,基本保持为一个常数,可以认为是一个稳定的时间序列。

均值预测

确认了时间序列稳定性的前提,我们接下来对下一个周期的数据做预测。均值预测算法包括如下几个步骤:

  • 将数据以月为维度进行汇总,得到每个月的付款总额
  • 将月份付款数据以年为周期进行划分,每年12个月的数据
  • 以月份为维度,将每年各个月份数据分别汇总,计算平均值,得到每个月在所有年份的平均值
  • 以年为维度,汇总每年的所有付款总额
  • 计算每一年相比于上一年的增长比率,得到一个增长率列表,取增长率列表的平均值
  • 以月份平均值列表乘以平均增长率,得到每个月的预测值列表
# 预测算法
def means_predict(df):
    # 所有月份汇总
    month_sum = df.groupby('period').sum()
    
    # 按年汇总
    df['year'] = df.index.year
    df_groupby_year = df.groupby('year').sum()
    df_groupby_year = df_groupby_year.loc[:, ['amount']]
    df_groupby_year = np.reshape(df_groupby_year.values, (len(df_groupby_year)))

    # 计算所有年份的平均比率,作为年增长率
    increment_rate_sum = 0
    for i in range(1, len(df_groupby_year)):
        increment_rate_sum += df_groupby_year[i] / df_groupby_year[i-1]
    increment_rate = increment_rate_sum / len(df_groupby_year)
    
    # 计算所有月份均值
    month_means = month_sum * increment_rate / len(df_groupby_year)
    month_means = month_means.loc[:, ['amount']]
    month_means = np.reshape(month_means.values, (len(month_means)))

    return month_means

按照以上算法进行计算后,我们把结果画出来,来看看预测结果的曲线。

def plot_predicts(dfs, titles, period='month'):
    assert len(dfs) != 0, '数据集不能为空'
    assert len(dfs) == len(titles), '数据集长度与标题集长度必须相等'
    for i in range(0, len(dfs)):
        df = dfs[i]
        df.dropna(inplace=True)

        df_groupby = groupby_df(df, period)
        month_predict = means_predict(df_groupby)
        ts = df_groupby['amount']

        # 插入历史数据最后一个点,确保两个曲线能连接在一起
        df_groupby_values = np.reshape(ts.values, (len(ts)))
        month_predict = np.insert(month_predict, 0, df_groupby_values[-1])

        rng = pd.date_range(ts.index[-1], periods=len(month_predict), freq='M')
        month_predict = pd.Series(month_predict, index=rng)

        plt.subplot(len(dfs), 1, i+1)
        
        plt.plot(ts, label='原值')
        plt.plot(month_predict, label='预测值')

        ax=plt.gca()
        xticks_range = pd.date_range(ts.index[0], month_predict.index[-1], freq='q')
        ax.set_xticks(xticks_range.strftime('%Y-%m'))
        ax.set_xticklabels(xticks_range.strftime('%Y-%m'))

        fig = plt.gcf()
        fig.set_size_inches(12, 4 * len(dfs))
        fig.autofmt_xdate()
        
        plt.subplots_adjust(right=0.99, left=0.125, bottom=0.1, top=0.9)
        plt.title(titles[i])
        plt.xlabel('月份')
        plt.ylabel('付款总额(亿元)')
        plt.grid(True)
        plt.legend(loc='best')
        
    plt.show()
    
# plot_predicts([df_index.copy(), df_area('B16_HDDQ0_华东大区')], ['资金使用趋势(全国)', '资金使用趋势(华东)'])
plot_predicts([df_area('全国')], ['资金使用趋势(全国)'])

这个结果还不错,基本保留了原有时间序列的规律,在相同的季节表现出相似的波峰和波谷,并且整体相对去年保持了一定的增长率。
接下来我们分别对各个大区的数据做预测,把结果汇总在一起,看看算法有没有普适性。

plot_predicts(df_areas, area_titles)

各个大区的预测曲线基本上保留了历史上的曲线形状,区别在于各个大区的增长率不同,有部分大区像华南、华北、东北,增长率并不明显,各个周期基本持平,部分大区像华西、华中,甚至出现了负增长。为什么会有这个情况,我们认为有两个方面的原因。

  • 与各大区的实际业务情况有关

  • 部分大区上系统较迟导致早期数据缺失

误差分析

接下来估算预测值是否拟合得足够好,通常来讲,只有采用有监督学习的机器学习算法,才需要评估模型输出的预测值与事先确定的真值之间的拟合程度,而这里并没有建立模型,只是通过均值和增长率来预测(很多情况下这比复杂的模型更简单有效),所以并不存在评估模型拟合的问题。不过为了从技术上量化预测值的准确度,让我们心里有个数,我们还是把预测值和上一年的数据做个拟合评估,以防止看起来曲线很漂亮,实际与历史数据存在很大出入的情况。

R2确定系数用于计算R²(确定系数:coefficient of determination),它可以评估算法对历史数据的拟合程度,用来度量未来的样本是否可能通过算法被很好地预测。分值为1表示最好,它可以是0,可以是负数(因为算法可以很糟糕)。

from sklearn.metrics import r2_score

df_groupby = groupby_df(df_index)
month_predict = means_predict(df_groupby)

years = np.unique(df_groupby.index.year)
r2 = r2_score(month_predict, ts[str(years[-1])])

print('R2确定系数: %.2f' % r2)
R2确定系数: 0.85

R2系数为0.85,可见预测值相比于上一年的值没有偏离太远,事实上也不能太接近,因为我们希望预测值相比于上一年是要保持增长率的,R2系数达到一个很高的值反而不是我们希望的。

输出结果

合并历史数据和预测值,输出结果,保存为文件,方便从业务上人为评估数据准确性。

def output(dfs, titles, console=False):
    assert len(dfs) != 0, '数据集不能为空'
    assert len(dfs) == len(titles), '数据集长度与标题集长度必须相等'

    for i in range(0, len(dfs)):
        df = dfs[i]
        df.dropna(inplace=True)

        df_groupby = groupby_df(df)
        month_predict = means_predict(df_groupby)
        ts = df_groupby['amount']

        rng = pd.date_range(df_groupby.index[-1]+1, periods=len(month_predict), freq='M')
        predict_ts = pd.Series(month_predict, index=rng)

        # 预测结果输出为文件
        X_pred = pd.concat([df_groupby, predict_ts], axis=1)
        X_pred.fillna(0, inplace=True)
        X_pred.rename(columns={0:'预测值(亿元)', 'amount':'真实值(亿元)'}, inplace = True)
        X_pred = X_pred.drop(['period'], axis=1)

        # 计算月度结果
        X_pred_month = X_pred.drop(['year'], axis=1)
        X_pred_month.index = X_pred.to_period().index.strftime('%Y-%m')
        X_pred_month.index.names = ['月份']

        if(console):
            print(X_pred, end='\n\n')
            
        # 计算年度结果
        X_pred['year'] = X_pred.index.year
        df_groupby_year = X_pred.groupby('year').sum()
        df_groupby_year.index.names = ['年份']

        # 输出文件
        output_file_name = titles[i] + '.xlsx'
        output_path = os.path.join(output_dir, output_file_name)
        
        # 保存为excel文件,输出到两个sheet
        with pd.ExcelWriter(output_path) as writer:
            X_pred_month.to_excel(writer, sheet_name='月度')
            df_groupby_year.to_excel(writer, sheet_name='年度')
        print('预测结果保存为:' + output_path)
        
        if(console):
            print(df_groupby_year.head(10), end='\n\n')
output(df_areas, area_titles)
预测结果保存为:D:/ml/资金计划预测\output\资金使用趋势(全国).xlsx
预测结果保存为:D:/ml/资金计划预测\output\资金使用趋势(大区1).xlsx
预测结果保存为:D:/ml/资金计划预测\output\资金使用趋势(大区2).xlsx
预测结果保存为:D:/ml/资金计划预测\output\资金使用趋势(大区3).xlsx
预测结果保存为:D:/ml/资金计划预测\output\资金使用趋势(大区4).xlsx
预测结果保存为:D:/ml/资金计划预测\output\资金使用趋势(大区5).xlsx
预测结果保存为:D:/ml/资金计划预测\output\资金使用趋势(大区6).xlsx

总结

本文介绍了一个使用均值和增长率对未来数据规律进行预测的简单算法,实践中这个预测结果被确认是有效可行的。

事实证明有的时候对未来的预测并不需要复杂的数学建模、机器学习、深度学习等高深的算法和模型,实际上我们尝试过ARIMA模型和LSTM深度学习模型,分析了ARIMA模型的p,d,q参数,尝试了LSTM模型的多种epoch迭代次数、神经元个数、隐层层数的组合,最终也没有收到非常满意的效果。倒是花了一小时实现的简单算法,最终解决了问题并获得认可(其实这个简单算法还是有一些优化空间)。

ARIMA(p,d,q)模型,全称为自回归积分滑动平均模型(Autoregressive Integrated Moving Average Model),是由博克思(Box)和詹金斯(Jenkins)于20世纪70年代初提出的一种时间序列预测方法。ARIMA模型的基本思想是,将预测对象随时间推移而形成的数据序列视为一个随机序列,用一定的数学模型来近似描述这个序列。这个模型一旦被识别后就可以从时间序列的过去值及现在值来预测未来值。根据原序列是否平稳以及回归中所含部分的不同,包括移动平均过程(MA)、自回归过程(AR)、自回归移动平均过程(ARMA)和自回归滑动平均混合过程(ARIMA)。

LSTM模型,全称长短期记忆(Long Short Term Memory)网络,是一种特殊的RNN(Recurrent Neural Networks递归神经网络)模型,1997年首先被Sepp Hochreiter和Jurgen Schmidhuber提出。适用于处理和预测时间序列中间隔和延迟相对较长的重要事件。随着深度学习的流行,LSTM已经在各个领域有了多种应用,基于LSTM的系统可以完成语言翻译、控制机器人、图像分析、文档摘要、语音识别图像识别、手写识别、控制聊天机器人、预测疾病、点击率和股票、合成音乐等等任务。

值得一提的是在准备数据时花费了不少时间,在前期的分析中,反反复复剃除了不少无关和无效数据,之前有人说这个项目的数据没有规律,没有可预测的前提,这就是没有用心处理数据而得到的草率判断。在本项目中,数据清洗和处理至少花掉了一半时间,我们认为这才是重点所在,数据是分析和预测的源头,没有高质量的数据,后面的事情无从谈起。

下面是常被问及的问题。

  • 实际业务中经常有某一时间集中获取项目,集中付款的情况,造成某个时间点出现付款尖峰,这样的情况是否会影响数据预测的结果,能否兼容这样的情况?

从两个方面来考虑这个问题。

  1. 我们考虑的是大区层面的月度汇总数据,在这个区域和时间颗粒度上来看,某一个或几个时间点上集中获取项目,这些集中的付款数据在被连续几年同一月份的数据平均后,对总体的影响会变得比较小,预测值在时间点和数额上会有少量偏差,但幅度不会过大,这个问题随着我们积累的历史数据越来越多,它的影响会越来越小。如果上升到全国层面来考虑这个问题,这些数据被平均后,影响就更小了。所以,在大区以上的月度数据这样的颗粒度下,算法是可以适用的。实际上我们也可以看到有部分大区上系统比较晚,也缺了一部分数据,但算法对这个情况保持了兼容性,也获得了比较接近的结果。如果实在有很异常的情况,我们可以结合实际业务情况,把数据当作噪音,做一定的平滑处理。
  2. 在实际业务环境中,集中获取项目往往受当时形势和各种因素的影响,很难表现出规律性,对于没有规律的数据,就像股市一样,算法很难去分析和预测它,并把它考虑到模型中。如果我们能够从业务上找出集中付款的规律,那么是可以对这些规律加以处理的,如果真的有规律,往往它也已经体现在日常的数据中。
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 194,319评论 5 459
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 81,801评论 2 371
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 141,567评论 0 319
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 52,156评论 1 263
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 61,019评论 4 355
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 46,090评论 1 272
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 36,500评论 3 381
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 35,192评论 0 253
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 39,474评论 1 290
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 34,566评论 2 309
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 36,338评论 1 326
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 32,212评论 3 312
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 37,572评论 3 298
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 28,890评论 0 17
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,169评论 1 250
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 41,478评论 2 341
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 40,661评论 2 335

推荐阅读更多精彩内容