《python for data analysis》第十章，时间序列

《 python for data analysis 》一书的第十章例程， 主要介绍时间序列（time series）数据的处理。
label:
1. datetime object、timestamp object、period object
2. pandas的Series和DataFrame object的两种特殊索引：DatetimeIndex 和 PeriodIndex
3. 时区的表达与处理
4. imestamp object、period object的频率概念，及其频率转换
5. 两种频率转换：单个时间object——asfreq；以时间object为索引的时间序列——resample
6. 时间序列的移动窗口（rolling）

# -*- coding:utf-8 -*-
# 《python for data analysis》 第十章
# 时间序列
import pandas as pd
import numpy as np
import time
from datetime import datetime, timedelta
import matplotlib.pyplot as plt

start = time.time()
np.random.seed(10)

# 1、日期和时间数据类型及工具
# datetime以 year - month - date hour : minute : second 的格式存储时间
tnow = datetime.now()
print(tnow)
# 可调用datetime中的各个部分：年月日时分秒
print(tnow.year)
print(tnow.day)
# datetime.timedelta 可用于表示两个 datetime object之间的时间差， 以日，秒，毫秒的形式
print(tnow - timedelta(1, 10, 1000))  # 1天10秒100微秒之前
print('\n')
# 1.1、字符串和datetime的相互转换
# （1）str和strftime可以将datetime object转换成字符串
print(str(tnow))  # 按照年月日时分秒的格式全部转换
print(tnow.strftime('%Y-%m-%d'))  # 按照年月日格式导出
print('')
# （2）strptime可以将字符串转换成datetime object
print(datetime.strptime('1999/1/2', '%Y/%m/%d'))
print(datetime.strptime('1999/1/2', '%Y/%d/%m'))
print('')
# datetime.strptime方法进行时间格式解析最为精准，但是需要自定义格式。
# 第三方包dateutil中的parser.parse方法可自适应地去进行时间格式解析，但有些时候可能会出现问题
from dateutil.parser import parse

print('Jan 22, 1999, 11:23, pm')
print(parse('Jan 22, 1999, 11:23, pm'))
print('')
# 上述方法均是对单个字符串进行解析，pandas的to_datetime方法可以成组解析
timestring = ['2012/1/1', '2008/8/8']
df = pd.to_datetime(timestring)
print(df)
# pd.to-datetime还可以对缺失值进行处理，自动转成NaT(not a time)
timestring = timestring + [None]
df = pd.to_datetime(timestring)
print(df)
print('------------------------------------------------↑, section1\n\n')
# 2、时间序列基础
# pandas最基本的时间序列类型是以时间戳(timestamp)为索引(index)的Series
dates = [datetime.strptime('2000/1/' + str(i), '%Y/%m/%d') for i in range(1, 11)]
print(dates)
ts = pd.Series(np.random.randn(10), index=dates)
print(ts)
print(type(ts.index[1]))  # time series中的每个index均为timestamp object
print('')
# 2.1、索引、选取、子集构造（time series的）
# time series也是series的一种，所以原来的series的各种索引、切片方法等均适用
print(ts[2:3])
print('')
# 特别的，time series还可以通过日期（字符串）进行索引
print(ts['2000/01/05'])
# 对于时间跨度很大的时间序列，切片方式更加丰富
print('\n')
ts = pd.Series(np.random.randn(1000), index=pd.date_range('1/1/2000', periods=1000))
print(ts.describe())
# 通过年份切片
print(ts['2001'].describe())
# 通过年月切片
print(ts['2001/01'].describe())
# 时间段切片
print(ts['2002/05/01':'2002/05/06'])
print('\n')
# 上述索引、切片方法对dataframe同样适用
# 2.2、带有重复索引的时间序列
date = ['2001/02/01', '2001/02/01', '2001/02/02']
ts = pd.Series(range(3), index=date)
print(ts)
print(ts['2001/02/01'])  # 重复索引返回切片
print(ts['2001/02/02'])  # 非重复索引返回标量值
print('')
# 通过聚合可以消除重复索引
print(ts.groupby(level=0).count())
print('---------------------------------↑, section2\n\n')
# 3、日期的范围、频率以及移动
# resample的详细介绍见section 6
# 3.1、生成日期范围
# 用pd.date_range生成指定长度的DateTime Index
# （1）给定起止时间
index = pd.date_range('2000/1/1', '2000/2/1')
print(index)
print('')
# （2）给定起始时间（或终止时间）加时间长度
index = pd.date_range(start='2000/1/1', periods=10)  # 从2000/1/1开始，长度10天
print(index)
print('')
index = pd.date_range(end='2000/1/28', periods=10)  # 到2000/1/28为止，长度10天
print(index)
# 上述函数中为给出生成时间序列的间隔，为缺省值1天，通过freq关键字可以显式指定
print('')
print(pd.date_range('2000/1/1', '2000/12/28', freq='BM'))  # BM表示每个月的最后一个工作日
# freq的可选项：D（每天）、B（每工作日）、H（每小时）、T or min（每分钟）、S（每秒）……P314~315
# 3.2、频率与日期偏移量
# 上述提到的频率代码可以 以字符串形式 自由组合
print('')
print(pd.date_range('2000/1/10', '2000/1/15', freq='6H30T10S'))  # 以6小时30分钟10秒的间隔生成时间序列
# WOM日期（week of month）可以表示每月中的某些日子
print('')
print(pd.date_range('2000/1/10', '2001/1/10', freq='WOM-3FRI'))  # 每个月的第三个周五
# 3.3、移动（超前或滞后）数据
print('\n')
ts = pd.DataFrame(
    {'a': range(6),
     'b': np.random.randn(6)}
)
ts.index = pd.date_range('2000/10/1', '2000/10/6')
print(ts)
print(ts.shift(2))  # 数据往时间增大方向推移2天
print(ts.shift(-2))  # 往小了方向推移
# 上述函数中仅传入推移大小，则移动数据
# 若还传入频率，则移动time index
print(ts.shift(1, freq='M'))  # 数据不动，时间序列全部推移到本月最后一天
print(ts.shift(1, freq='3D'))  # 数据不动，时间序列全部推移3天
# shift最大的功用在于可以计算数据的百分比变化
print('')
print(ts / ts.shift(1) - 1)
# 时间偏移量还可以直接作用于timestamp和datetime object
# 先导入时间偏移量
from pandas.tseries.offsets import Day, BMonthBegin

print(tnow + 3 * Day())  # 向时间前进方向偏移3天
print(tnow + BMonthBegin())  # 偏移至下个月的第一个工作日
print('------------------------------------↑, section 3\n\n')

# 4、时区处理
# UTC，协调世界时间。时区以UTC的偏移量表示
# 4.1、本地化与转换
# 先按照之前的方法创建一个time series
ts = pd.Series(range(2))
ts.index = pd.date_range(tnow, periods=2, normalize=True)
print(ts)
# 上述创建过程中未指定时区，则默认时区未None
print('\ntime zone is %s' % ts.index.tz)
# 给时间序列加上时区数据，称为本地化。有两种方法。
# （1）tz_localize方法
ts_china = ts.tz_localize('Asia/Shanghai')
print(ts_china)  # time index 变成 本地时间+UTC偏移 的格式
print(ts_china.index.tz)  # 该时间序列的时区已经附上上海
# （2）在创建序列时直接显示指定时区
ts_china = pd.Series(range(2), index=pd.date_range(tnow, periods=2, normalize=True, tz='Asia/Shanghai'))
print(ts_china)
print(ts_china.index.tz)
# 时区可以通过tz_convert方法进行转换
print(ts_china.tz_convert('UTC'))  # 北京时间转换成UTC时间
print('\n')
# datetime、timestamp、DatetimeIndex这些objects均可以使用tz_localize、tz_convert这些方法
# 4.2、操作时区意识型TimeStamp Object
stamp = pd.Timestamp(str(tnow))
stamp_china = stamp.tz_localize('Asia/Shanghai')
print(stamp_china)
stamp_utc = stamp_china.tz_convert('utc')
print(stamp_utc)
# timestamp object 中有一个属性保存了utc时间戳值，即当前时间相对于UNIX纪元（1970年1月1日）的时间位移，以ns为单位
print(stamp_utc.value)
print(stamp_china.value)  # 这两个时间为stamp在不同时区的显示时间，故绝对位移相等，均为stamp相对UNIX纪元的时间位移
print('\n')
# 4.3、不同时区之间的运算
# 不同时区的时间序列的运算结果均以UTC标准显示
ts = pd.Series(range(2), index=pd.date_range(tnow, periods=2))
ts1 = ts.tz_localize('US/Eastern')
ts2 = ts.tz_localize('Asia/Shanghai')
print(ts1.index)
print(ts2.index)
print((ts1 + ts2).index)
print('----------------------------------↑, section4\n\n')

# 5、时期及其算术运算
# 一种新的object，时期（period）
# timestamp object表示某一时刻，相对的，period object是用来表示某一段周期的
# 创建一个period object
p = pd.Period(2000, freq='A-DEC')  # 以12月结尾的周期年
print(p)
# period object支持加减整数实现位移
print(p - 2)  # 2000年往前推2年
print(p + 1)  # 2000年的后面一年
# 同频率的period还支持加减
print(pd.Period(2005, freq='M') - pd.Period(2000, freq='M'))  # 5年共60个月
# period_range方法可以创建一组时期范围
plist = pd.period_range('2000Q1', '2002Q1', freq='Q')  # 从2000年1季度到2002年1季度，季度间隔
print(plist)
print('')
# 以datetime object为内容的index称为datetime index，同理，也有period index
ts = pd.Series(np.random.randn(len(plist)), index=plist)
print(ts)
print('')
# period index的构造还可以通过字符串转换完成
string = ['2000Q1', '2000Q2', '2000Q3']  # 2000年的前3季度
p = pd.PeriodIndex(string, freq='Q-DEC')
print(p)
print('\n')
# 5.1、时期的频率转换
# period 和 period index 均可通过asfreq进行频率转换
p = pd.Period(2000, freq='Q-DEC')
print(p)  # 2000年第一季度
print(p.asfreq('M'))  # 季转月，默认最后一个月
print(p.asfreq('M', how='start'))  # 用how关键字显式指定第一个月
print(p.asfreq('A'))  # 季转年
print('')
# 同理，period index或者包含period index的time series也可以这么操作
print(ts.asfreq('B', how='start'))  # ts时间序列中的period index（以季度为时期）转成以每季里面第一个工作日为时期的period index
print('')
# 5.2、按季度计算的时期频率
# Q-MAY 中的 MAY 表示该年末为五月，即6-8一个季度，9-11一个季度，12-2一个季度，3-5一个季度。 其余表示以此类推
# 单个 period object 和 一组period（如period index）均可通过运算来表示某个时刻，并通过to_timestamp方法变成timestamp object
a = pd.Period('2000Q1', freq='Q')
print(a)
# 通过period运算将2000年第一季度转换为2000年第一季度的倒数第三个工作日的上午9点30分的时间戳
tstamp = ((a.asfreq('B', 'e') - 2).asfreq('H', 's') + 9).asfreq('T', 's') + 30
print(tstamp)  # 还是一个period object，不过是分钟级的period
print(tstamp.to_timestamp())  # period 2 timestamp
print('')
# 5.3、将timestamp转换为period（及其反向过程）
# （1）to_period方法可以将timestamp转换为period
stamp = pd.date_range('2000/12/1', periods=2, freq='D')
print(stamp.to_period('M'))  # 指定period的频率为月
# （2）to_timestamp方法可以将period转换为timestamp
print(stamp.to_period('M').to_timestamp())
print('')
# 5.4、通过数组创建period index
# 很多时候时间数据是分成几个子部分（如年、月、日）存放在一张表格的某几列中的，PeriodIndex方法可以合并这样的列汇成一个period index
df = pd.DataFrame()
df['year'] = [2000] * 4 + [2001] * 4
df['month'] = range(1, 9)
df['day'] = range(22, 30)
print(df)
tindex = pd.PeriodIndex(year=df['year'], month=df['month'], day=df['day'], freq='D')
print(tindex)
print('---------------------------------------↑, section 5 \n\n')

# 6、重采样与频率转换
# 重采样是指从一个频率变换到另一种频率的过程，这里的频率并不是传统意义上的频率，而是指时间序列相关函数中freq关键字
# 重采样分3种：
# 1、升采样，频率变大，或者说period周期变小，如Q采样成D（季度->天）
# 2、降采样，频率变小，或者说period周期变大，如Q采样成A（季度->年）
# 3、同period周期的采样，如W-MON采样成W-WED（每周一->每周三）
# 重采样通过resample方法实现
# 6.1、降采样
# 先生成一个分钟序列数据
ts = pd.Series(range(15), index=pd.date_range('2010/10/1', periods=15, freq='T'))
print(ts)
print('')
# resample方法需要指定区间哪一边为闭区间（相应另一边为开区间），需要指定区间以左右哪个边界进行命名
# 在0.22.0版本（更高版本估计也类似）的pandas种，缺省情况比较复杂，视freq不同而不同，所以保险一点还是用关键字显示指定。
# 可通过closed关键字修改哪一边为闭区间，可通过label关键字指定以区间的哪个边界进行区间命名
print(ts.resample('10min', closed='left', label='left').count())
print('')
print(ts.resample('10min', closed='right', label='right').count())
# 为更清楚地显示区间，可使用loffset关键字进行偏移。其实也可以通过对整个time series进行shift实现相同的功能。
print('')
print(ts.resample('10min', closed='right', label='right', loffset='-1s').count())
# 特别地，降采样中有一种采样称为OHLC采样，特定用于金融数据，计算出每个区间的开盘价(open)，收盘价(close)，最高价(high)，最低价(low)
print('')
print(ts.resample('5min', closed='right', label='right').ohlc())
# 另一种实现降采样的方法是通过groupby方法实现，和resample各有适用情景
# 按照周几进行分组
ts = pd.Series(range(100), index=pd.date_range('2000/10/10', periods=100, freq='D'))
print(ts)
print('')
print(ts.groupby(lambda t: t.weekday).count())
# 6.2、升采样和插值
# 从大尺度采样到小尺度，自然而然会引入数值缺失的问题，故升采样需要插值处理
# 先创建一个周频率的时间序列
print('')
ts = pd.Series([1, 2], index=pd.date_range('2000/1/21', periods=2, freq='W-FRI'))  # 从2000/1/21开始的两个周五
print(ts)
# 将ts升采样到日频率
# 若不插值则会出现数值缺失，可通过前向插值(ffill)和后向插值(bfill)进行插值
print('')
print(ts.resample('D').ffill())
print(ts.resample('D').bfill())
# resample还可以实现既非升采样也非降采样
# 如将上述每周5的数据重采样到每周1
print('')
print(ts.resample('W-MON').ffill())
print(ts.resample('W-MON').bfill())
# 6.3、通过时期进行重采样
# 上述重采样都是对timestamp index的series进行的，重采样也可以对period index的series进行
# （1）period的降采样和timestamp一样，直接对其应用resample方法即可
# 先构造一个period index的时间序列
print('')
ts = pd.Series(range(10), index=pd.period_range('2010/1', periods=10, freq='M'))
print(ts)
# 然后在resample中传入一个更大周期的freq即可完成降采样
ts = ts.resample('Q').sum()
print(ts)
print('')
# （2）period的升采样，需要指定新区间中哪一端存放原来的值
# 比如，年到季度的升采样，原来的值放到第一个季度还是最后一个季度需要指定
# 通过关键字convention进行指定，'end'则放到最后一个季度，'start'则放到第一个季度。缺省为'start'
print(ts.resample('M', convention='end').ffill())
print('')
print(ts.resample('M').ffill())
print('--------------------------↑, section 6\n\n')

# 7、时间序列绘图
# pandas的时间序列数据可直接用plot()方法进行绘图，基于matplotlib包进行过处理的
# 导入数据，1990年至2010年的几只美股数据
stk = pd.read_csv('./data_set/stock_px.csv', parse_dates=True, index_col=0)
stk = stk[['AAPL', 'MSFT', 'SPX']]  # 从中取出3只股票
stk = stk.resample('B').ffill()  # 按工作日频率进行重采样，实现规则频率
print(stk.describe())
# 通过切片直接应用plot方法可完成时间序列的作图
fig, axes = plt.subplots(2, 2)  # figure1
stk['AAPL'].plot(ax=axes[0, 0])  # 品种切片
stk.ix['2005'].plot(ax=axes[0, 1])  # 时间切片
stk['AAPL'].ix['06/2006':'08/2008'].plot(ax=axes[1, 0])  # 双重切片
# 还可以对原数据重采样成季度数据，再作图
stk['AAPL'].resample('Q-DEC').ffill().plot(ax=axes[1, 1])  # 默认为每个季度最后一个工作日，若数据缺失则用前一天的数据填充
# plt.show() # 将这一行取消注释，使能作图功能
print('-------------------------------------↑, section 7\n\n')

# 8、移动窗口函数
# 移动窗口函数，用于在一个长序列中切割出一个子窗口进行相关量的统计
# 在金融数据中移动窗口应用较多，典型地，N日均线
# 对APPLE股价作250日移动平均
fig2, axes2 = plt.subplots(2, 2)  # figure2
stk['AAPL'].plot(ax=axes2[0, 0])
pd.Series.rolling(stk['AAPL'], 250).mean().plot(ax=axes2[0, 0])  # 高版本推荐写法，不同于书上例程
# 当rolling接受的数据很少时，将不返回移动平均值，通过min_periods关键字可以指定这个阈值
stk['AAPL'].plot(ax=axes2[0, 1])
pd.Series.rolling(stk['AAPL'], 250, min_periods=10).mean().plot(ax=axes2[0, 1])  # 最少有10个非NA值就返回移动平均
# 图（0，0）和图（0，1）的区别体现在图（0，1）更快出现250日均线
# 通过rolling还可以延申成 扩展窗口平均， 即窗口长度可变，相当于时间序列长度
expanding_mean = lambda ts: pd.Series.rolling(ts, len(ts), min_periods=1).mean()
expanding_mean(stk['AAPL']).plot(ax=axes2[1, 0])  # 全长度均线
stk['AAPL'].plot(ax=axes2[1, 0])
# 8.1、指数加权函数
# 移动窗口常搭配一个衰减因子使用，用于使近期的观测值有更大的权重，从而更快体现原数据的变化
# 通过ewm方法实现
ma = pd.Series.rolling(stk['AAPL'], 250, min_periods=50).mean()  # 均权的年均线
ewma = pd.Series.ewm(stk['AAPL'], span=250).mean()
ma.plot(ax=axes2[1, 1], style='--')
ewma.plot(ax=axes2[1, 1], style=':')  # 通过图可以看到带有衰减因子的均线更快出现拐点（反应更快）
# 8.2、二元移动窗口函数
# 某些统计变量用到两个数据，比如相关系数
spx = stk['SPX']
aapl = stk['AAPL']
# 两种方法计算股价的百分数变化
spx_pctc = spx / spx.shift(1) - 1
aapl_pctc = stk['AAPL'].pct_change()
fig3, axes3 = plt.subplots(2, 2)  # figure3
# 计算两者的移动窗口里的相关系数
corr = pd.Series.rolling(spx_pctc, window=125, min_periods=100).corr(aapl_pctc) # 6个月窗口期，移动相关系数
corr.plot(ax=axes3[0,0])
# 很多时候会以唯一数据作为标准，计算其余数据与标准数据的相关系数
# 此时，传入Series（标准数据）与DataFrame（其余数据）即可
corr = pd.DataFrame.rolling(stk[['AAPL', 'MSFT']].pct_change(), window=125, min_periods=60).corr(spx_pctc)
corr.plot(ax=axes3[0,1])
# 8.3、用户定义的移动窗口函数
ten_mean = lambda ts: np.mean(sorted(ts, reverse=True)[:10]) # 计算最大的前10个值得平均值
res = pd.DataFrame.rolling(stk[['AAPL', 'MSFT']], window=125, min_periods=50).apply(ten_mean)
print(res)
res.plot(ax=axes3[1,0])
plt.show()
# 其实，在高版本的pandas中，rolling的用法和groupby是比较接近的。
print('----------------------total time is %.5f s' % (time.time() - start))