• Posts tagged "时间序列"

Blog Archives

用timetk进行时间序列可视化

R的极客理想系列文章,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,互联网….都在使用R语言。

要成为有理想的极客,我们不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让我们一起动起来吧,开始R的极客理想。

关于作者

  • 张丹,数据分析师/程序员/Quant: R,Java,Nodejs
  • blog: http://fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-timetk/

前言

时间序列分析,是统计分析一个重要的分支,我们日常生活的很多数据,都是时间序列数据,至少包含一个时间项和数据值。如果能很好地观察时间序列数据,对于我们理解数据能起到非常重要的作用。

timetk就是一个非常方便地观察时间序列数据的包,从多个角度让我们快速理解数据。

目录

  1. timetk包是介绍
  2. timetk的安装
  3. timetk的核心函数
  4. 可视化:基本时间序列图
  5. 可视化:时间序列分组
  6. 可视化:时间序列分箱
  7. 可视化:时间序列回归
  8. 可视化:时间序列异常检测
  9. 可视化:时间序列季节性分解
  10. 可视化:时间序列重采样

1. timetk包是介绍

timetk是R语言的一个时间序列处理工具包,可轻松实现时间序列数据的可视化、整理和特征工程,用于预测和机器学习建模。整合并扩展了包括 dplyr、stats、xts、forecast、slider、padr、recipes 和 rsample 等包中的时间序列功能。Timetk 是一个功能强大的包,属于 modeltime 生态系统的一部分,用于时间序列分析和预测。

官方网站:https://business-science.github.io/timetk/

timetk 包让用户在 R 语言中更便捷地处理时间序列对象。该工具包提供了用于检查和处理基于时间的索引、扩展用于数据挖掘和机器学习的时间特征,以及在不同时间序列类之间进行转换的功能。

主要优势如下:

  • 索引提取:从任意时间序列对象中提取时间序列索引。
  • 时间序列理解:基于时间序列索引生成特征摘要和汇总信息。
  • 构建未来时间序列:根据已有索引创建未来时间序列。
  • 基于时间的tibble类型与主流时间序列数据类型(xts、zoo、zooreg、ts)之间的相互转换:简化转换过程,并在转换为规则时间序列(如 ts)时最大程度保留基于时间的数据信息。

2. timetk的安装

timetk包的安装,安装timetk包非常简单,2条命令就可以完成,安装和加载。


# 安装
> install.packages("timetk")

# 加载
> library(timetk)

3. timetk的核心函数

使用timetk包,可以针对不同R语言时间序列类型做适配,如tk_tbl()适配tibble类型, tk_zoo()适配zoo类型,tk_xts()适配xts类型等。

timetk包中最核心的函数就是画图函数,支持11种可视化的画图效果。

  • plot_time_series 一个或多个时间序列的交互式绘图
  • plot_time_series_boxplot 交互式时间序列箱线图
  • plot_time_series_regression 可视化时间序列线性回归公式
  • plot_anomalies 可视化一个或多个时间序列的异常值
  • plot_anomalies_cleaned 可视化一个或多个时间序列的异常值(清理后)
  • plot_anomalies_decomp 可视化一个或多个时间序列的异常值(分解后)
  • plot_anomaly_diagnostics 可视化一个或多个时间序列的异常值诊断
  • plot_seasonal_diagnostics 可视化一个或多个时间序列的多种季节性特征
  • plot_stl_diagnostics 可视化一个或多个时间序列的 STL 分解特征
  • plot_acf_diagnostics 可视化一个或多个时间序列的 ACF、PACF 和 CCF
  • plot_time_series_cv_plan 可视化时间序列重采样方案

4. 可视化:基本时间序列图

以timetk包自带的数据集taylor_30_min为例,taylor_30_min数据集。为半小时电力需求,2000年6月5日星期一至2000年8月27日星期日期间英格兰和威尔士的半小时电力需求。

taylor_30_min数据集为 tibble类型, 4,032 行 × 2 列:

  • date:日期时间变量,间隔为 30 分钟
  • value:电力需求,单位为兆瓦

查看数据情况


# 查看数据前6条
> head(taylor_30_min)
# A tibble: 6 × 2
  date                value
                
1 2000-06-05 00:00:00 22262
2 2000-06-05 00:30:00 21756
3 2000-06-05 01:00:00 22247
4 2000-06-05 01:30:00 22759
5 2000-06-05 02:00:00 22549
6 2000-06-05 02:30:00 22313

# 查看数据2列统计概览
> summary(taylor_30_min)
      date                         value      
 Min.   :2000-06-05 00:00:00   Min.   :18640  
 1st Qu.:2000-06-25 23:52:30   1st Qu.:24272  
 Median :2000-07-16 23:45:00   Median :29485  
 Mean   :2000-07-16 23:45:00   Mean   :29617  
 3rd Qu.:2000-08-06 23:37:30   3rd Qu.:35132  
 Max.   :2000-08-27 23:30:00   Max.   :38777  

基于数据画图,使用基本plot()函数画图。


> plot(taylor_30_min,type="l")

使用timetk画图,并进行数据拟合


# 加载包
> library(dplyr)
> library(ggplot2)
> library(lubridate)
> library(magrittr)

# 生成基于ggplot2的静态图
> interactive <- FALSE # 画图拟合 > taylor_30_min %>% 
+   plot_time_series(date, value, 
+                    .interactive = interactive,
+                    .plotly_slider = TRUE)
Ignoring unknown labels:
• colour : "Legend"

按月度进行数据拆分拟合


# 画图拟合
> taylor_30_min %>%
+   plot_time_series(date, value, 
+                    .color_var = month(date, label = TRUE),
+                    .interactive = interactive,  
+                    .title = "Taylor's MegaWatt Data",
+                    .x_lab = "Date (30-min intervals)",
+                    .y_lab = "Energy Demand (MW)",
+                    .color_lab = "Month") +
+   scale_y_continuous(labels = scales::label_comma())

使用可交互的画图输出, interactive=TRUE。可生成基于plotly的动态交互的结果。


> taylor_30_min %>% 
+     plot_time_series(date, value, 
+                      .interactive = TRUE, 
+                      .plotly_slider = TRUE)
Ignoring unknown labels:
• colour : "Legend"

5. 可视化:时间序列分组

1. 分组展示

我们再换一组数据集m4_daily,进行演示测试。m4_daily数据集来源:M4 竞赛(2018 年 1 月 – 5 月),完整竞赛总序列数:100,000 条,本样本:4 条日度时间序列,行数:9,743,列:id、date、value。

  • id 因子型 4 条序列各自的唯一标识符
  • date 日期型 每日时间戳
  • value 数值型 对应时间戳的数值

由于这些数据来自 M4 竞赛,每条序列很可能来自不同的领域(例如微观、工业、宏观、金融、人口统计等),并可能表现出:趋势,季节性(周度、年度),异方差性,零值或间歇性数值,不同长度(日度序列的历史数据长度往往各不相同)。


# 查看数据集
> head(m4_daily)
# A tibble: 6 × 3
  id    date       value
        
1 D10   2014-07-03 2076.
2 D10   2014-07-04 2073.
3 D10   2014-07-05 2049.
4 D10   2014-07-06 2049.
5 D10   2014-07-07 2006.
6 D10   2014-07-08 2018.

# 统计概览
> summary(m4_daily)
    id            date                value      
 D10 : 674   Min.   :1978-06-23   Min.   : 1735  
 D160:4197   1st Qu.:2003-01-12   1st Qu.: 5719  
 D410: 676   Median :2006-05-14   Median : 8253  
 D500:4196   Mean   :2005-02-14   Mean   : 8280  
             3rd Qu.:2009-09-13   3rd Qu.:11219  
             Max.   :2016-05-06   Max.   :19433  

分组数据的可视化只需在将数据传入 plot_time_series() 函数之前使用 group_by() 对数据集进行分组即可。关键点如下:

  • 可以通过两种方式添加分组:使用 group_by(),或通过 … 参数添加分组。
  • 分组随后会转换为分面图。
  • .facet_ncol = 2 生成一个两列的分面图。
  • .facet_scales = “free” 允许每个图的 x 轴和 y 轴独立于其他图进行缩放。

> m4_daily %>%
+   group_by(id) %>%
+   plot_time_series(date, value,
+                    .facet_ncol = 2, .facet_scales = "free",
+                    .interactive = interactive)
Ignoring unknown labels:
• colour : "Legend"

2. 对数变换与子分组

让我们切换到一个包含多个分组的小时级数据集。我们可以展示以下内容。其意图是在分面图中展示分组,同时突出数据中的每周窗口(子分组),并对数值进行 log() 变换。

  • .value = log(value):应用对数变换
  • .color_var = week(date):将日期列转换为 lubridate::week() 的周数,颜色将应用于每个周数。

换成m4_hourly,M4 竞赛始于 2018 年 1 月 1 日,于 2018 年 5 月 31 日结束。该竞赛包含 100,000 个时间序列数据集。本数据集包含来自该竞赛的 4 个小时度时间序列样本。


> m4_hourly %>%
+   group_by(id) %>%
+   plot_time_series(date, log(value),             # Apply a Log Transformation
+                    .color_var = week(date),      # Color applied to Week transformation
+                    # Facet formatting
+                    .facet_ncol = 2,
+                    .facet_scales = "free",
+                    .interactive = interactive)

6. 可视化:时间序列分箱

plot_time_series_boxplot() 函数可用于绘制箱线图。箱线图使用聚合方式,这是由 .period 参数定义的分箱参数。

换成m4_monthly,M4 竞赛始于 2018 年 1 月 1 日,于 2018 年 5 月 31 日结束。该竞赛包含 100,000 个时间序列数据集。本数据集包含来自该竞赛的 4 个月度时间序列样本。


> m4_monthly %>%
+   group_by(id) %>%
+   plot_time_series_boxplot(
+     date, value,
+     .period      = "1 year",
+     .facet_ncol  = 2,
+     .interactive = FALSE)
Ignoring unknown labels:
• colour : "Legend"

7. 可视化:时间序列回归

使用时间序列回归图函数 plot_time_series_regression() 可用于快速评估与时间序列相关的关键特征。该函数内部会将一个公式传递给 stats::lm() 函数。通过设置 show_summary = TRUE,可以输出线性回归的摘要信息。

换成m4_monthly数据画图


> m4_monthly %>%
+   group_by(id) %>%
+   plot_time_series_regression(
+     .date_var     = date,
+     .formula      = log(value) ~ as.numeric(date) + month(date, label = TRUE),
+     .facet_ncol   = 2,
+     .interactive  = FALSE,
+     .show_summary = FALSE
+   )

8. 可视化:时间序列异常检测

使用walmart_sales_weekly数据集, 可视化一个或多个时间序列的异常值。

walmart_sales_weekly数据集,沃尔玛招聘——门店销售额预测”竞赛使用了零售数据,涉及门店与店内各个部门的组合。该竞赛于2014年2月20日开始,于2014年5月5日结束。竞赛数据包含位于不同地区的45家零售门店。数据集包含多种外部特征,包括节假日信息、气温、燃油价格和降价信息。本数据集包含门店ID为1的7个部门的样本(共7个时间序列)。

查看数据walmart_sales_weekly


# 查看数据
> head(walmart_sales_weekly)
# A tibble: 6 × 17
  id    Store  Dept Date       Weekly_Sales IsHoliday Type    Size Temperature Fuel_Price MarkDown1 MarkDown2 MarkDown3 MarkDown4
                                                          
1 1_1       1     1 2010-02-05       24924. FALSE     A     151315        42.3       2.57        NA        NA        NA        NA
2 1_1       1     1 2010-02-12       46039. TRUE      A     151315        38.5       2.55        NA        NA        NA        NA
3 1_1       1     1 2010-02-19       41596. FALSE     A     151315        39.9       2.51        NA        NA        NA        NA
4 1_1       1     1 2010-02-26       19404. FALSE     A     151315        46.6       2.56        NA        NA        NA        NA
5 1_1       1     1 2010-03-05       21828. FALSE     A     151315        46.5       2.62        NA        NA        NA        NA
6 1_1       1     1 2010-03-12       21043. FALSE     A     151315        57.8       2.67        NA        NA        NA        NA
# ℹ 3 more variables: MarkDown5 , CPI , Unemployment 

# 统计概览
> summary(walmart_sales_weekly)
       id          Store        Dept            Date             Weekly_Sales    IsHoliday           Type          
 1_1    :143   Min.   :1   Min.   : 1.00   Min.   :2010-02-05   Min.   :  6166   Mode :logical   Length:1001       
 1_3    :143   1st Qu.:1   1st Qu.: 3.00   1st Qu.:2010-10-08   1st Qu.: 28257   FALSE:931       Class :character  
 1_8    :143   Median :1   Median :13.00   Median :2011-06-17   Median : 39886   TRUE :70        Mode  :character  
 1_13   :143   Mean   :1   Mean   :35.86   Mean   :2011-06-17   Mean   : 54646                                     
 1_38   :143   3rd Qu.:1   3rd Qu.:93.00   3rd Qu.:2012-02-24   3rd Qu.: 77944                                     
 1_93   :143   Max.   :1   Max.   :95.00   Max.   :2012-10-26   Max.   :148798                                     
 (Other):143                                                                                                       
      Size         Temperature      Fuel_Price      MarkDown1         MarkDown2          MarkDown3          MarkDown4      
 Min.   :151315   Min.   :35.40   Min.   :2.514   Min.   :  410.3   Min.   :    0.50   Min.   :    0.25   Min.   :    8.0  
 1st Qu.:151315   1st Qu.:57.79   1st Qu.:2.759   1st Qu.: 4039.4   1st Qu.:   40.48   1st Qu.:    6.00   1st Qu.:  577.1  
 Median :151315   Median :69.64   Median :3.290   Median : 6154.1   Median :  144.87   Median :   25.96   Median : 1822.5  
 Mean   :151315   Mean   :68.31   Mean   :3.220   Mean   : 8090.8   Mean   : 2941.32   Mean   : 1225.40   Mean   : 3746.1  
 3rd Qu.:151315   3rd Qu.:80.49   3rd Qu.:3.594   3rd Qu.:10122.0   3rd Qu.: 1569.00   3rd Qu.:  101.64   3rd Qu.: 3750.6  
 Max.   :151315   Max.   :91.65   Max.   :3.907   Max.   :34577.1   Max.   :46011.38   Max.   :55805.51   Max.   :32403.9  
                                                  NA's   :644       NA's   :707        NA's   :651        NA's   :644      
   MarkDown5            CPI         Unemployment  
 Min.   :  554.9   Min.   :210.3   Min.   :6.573  
 1st Qu.: 3127.9   1st Qu.:211.5   1st Qu.:7.348  
 Median : 4325.2   Median :215.5   Median :7.787  
 Mean   : 5018.7   Mean   :216.0   Mean   :7.610  
 3rd Qu.: 6222.2   3rd Qu.:220.6   3rd Qu.:7.838  
 Max.   :20475.3   Max.   :223.4   Max.   :8.106  
 NA's   :644                          

> walmart_sales_weekly %>%
+   filter(id %in% c("1_1", "1_3")) %>%
+   group_by(id) %>%
+   anomalize(Date, Weekly_Sales) %>%
+   plot_anomalies(Date, .facet_ncol = 2, .ribbon_alpha = 0.25, .interactive = FALSE)
frequency = 13 observations per 1 quarter
trend = 52 observations per 1 year
frequency = 13 observations per 1 quarter
trend = 52 observations per 1 year

可视化一个或多个时间序列的异常值,自动进行数据清理后的结果。


> walmart_sales_weekly %>%
+   filter(id %in% c("1_1", "1_3")) %>%
+   group_by(id) %>%
+   anomalize(Date, Weekly_Sales, .message = FALSE) %>%
+   plot_anomalies_cleaned(Date, .facet_ncol = 2, .interactive = FALSE)

把数据按时间序列的多维度进行分解。


> walmart_sales_weekly %>%
+   filter(id %in% c("1_1", "1_3")) %>%
+   group_by(id) %>%
+   anomalize(Date, Weekly_Sales, .message = FALSE) %>%
+   plot_anomalies_decomp(Date, .interactive = FALSE)

把数据按时间序列进行异常诊断。


> walmart_sales_weekly %>%
+   group_by(id) %>%
+   plot_anomaly_diagnostics(Date, Weekly_Sales,
+                            .message = FALSE,
+                            .facet_ncol = 3,
+                            .ribbon_alpha = 0.25,
+                            .interactive = FALSE)

9. 可视化:时间序列季节性分解

季节分析,可视化一个或多个时间序列的多种季节性特征,分别按日期,星期,周,月的4个维度,生成可视化看板。


> taylor_30_min %>%
+   plot_seasonal_diagnostics(date, value, .interactive = FALSE)

STL分解,按STL分解特征拆分,基于局部加权回归的季节性趋势分解。STL分解是一种将时间序列分解为三个组成部分的统计方法:

  • 趋势分量 (Trend) 数据的长期变化方向,反映平滑的总体走势
  • 季节分量 (Seasonal) 周期性重复的模式,如每周、每月或每年的规律性波动
  • 余项分量 (Remainder) 剔除趋势和季节后剩余的随机波动或噪声

数学表达式为:
观测值 = 趋势 + 季节 + 余项


> m4_hourly %>%
+   filter(id == "H10") %>%
+   plot_stl_diagnostics(
+     date, value,
+     # Set features to return, desired frequency and trend
+     .feature_set = c("observed", "season", "trend", "remainder"),
+     .frequency   = "24 hours",
+     .trend       = "1 week",
+     .interactive = FALSE)
frequency = 24 observations per 24 hours
trend = 168 observations per 1 week

可视化ACF, PACF, and CCFs检测


> m4_hourly %>%
+   group_by(id) %>%
+   plot_acf_diagnostics(
+     date, value,               # ACF & PACF
+     .lags = "7 days",          # 7-Days of hourly lags
+     .interactive = FALSE
+   )

10. 可视化:时间序列重采样

使用FANG数据集,FANG是包含科技股(“FB”、“AMZN”、“NFLX”和“GOOG”)每日历史股价的数据集,时间跨度从2013年初到2016年底。FANG包含 4,032 行和 8 个变量的“tibble”(“整洁”数据框):symbol股票代码符号,开高低收等信息。


> head(FANG)
# A tibble: 6 × 8
  symbol date        open  high   low close    volume adjusted
                     
1 FB     2013-01-02  27.4  28.2  27.4  28    69846400     28  
2 FB     2013-01-03  27.9  28.5  27.6  27.8  63140600     27.8
3 FB     2013-01-04  28.0  28.9  27.8  28.8  72715400     28.8
4 FB     2013-01-07  28.7  29.8  28.6  29.4  83781800     29.4
5 FB     2013-01-08  29.5  29.6  28.9  29.1  45871300     29.1
6 FB     2013-01-09  29.7  30.6  29.5  30.6 104787700     30.6

# 数据概览
> summary(FANG)
    symbol               date                 open              high              low              close       
 Length:4032        Min.   :2013-01-02   Min.   :  22.99   Min.   :  23.09   Min.   :  22.67   Min.   :  22.9  
 Class :character   1st Qu.:2014-01-01   1st Qu.: 106.31   1st Qu.: 107.75   1st Qu.: 104.50   1st Qu.: 106.2  
 Mode  :character   Median :2015-01-01   Median : 335.67   Median : 340.68   Median : 331.20   Median : 335.1  
                    Mean   :2015-01-01   Mean   : 382.70   Mean   : 386.38   Mean   : 378.67   Mean   : 382.6  
                    3rd Qu.:2016-01-01   3rd Qu.: 581.47   3rd Qu.: 585.86   3rd Qu.: 576.52   3rd Qu.: 582.1  
                    Max.   :2016-12-30   Max.   :1226.80   Max.   :1228.88   Max.   :1218.60   Max.   :1220.2  
     volume             adjusted     
 Min.   :     7900   Min.   : 13.14  
 1st Qu.:  2699875   1st Qu.: 74.82  
 Median :  7285600   Median :190.75  
 Mean   : 16489749   Mean   :297.12  
 3rd Qu.: 21942350   3rd Qu.:531.37  
 Max.   :365457900   Max.   :844.36  

数据处理,重采样


> FB_tbl <- FANG %>%
+   filter(symbol == "FB") %>%
+   select(symbol, date, adjusted)

# 重采样参数配置
> resample_spec <- time_series_cv(
+   FB_tbl,
+   initial = "1 year",
+   assess  = "6 weeks",
+   skip    = "3 months",
+   lag     = "1 month",
+   cumulative  = FALSE,
+   slice_limit = 6
+ )
Using date_var: date

# 重采样输出
> resample_spec %>% tk_time_series_cv_plan()
# A tibble: 1,812 × 5
   .id    .key     symbol date       adjusted
                   
 1 Slice1 training FB     2015-11-19     106.
 2 Slice1 training FB     2015-11-20     107.
 3 Slice1 training FB     2015-11-23     107.
 4 Slice1 training FB     2015-11-24     106.
 5 Slice1 training FB     2015-11-25     105.
 6 Slice1 training FB     2015-11-27     105.
 7 Slice1 training FB     2015-11-30     104.
 8 Slice1 training FB     2015-12-01     107.
 9 Slice1 training FB     2015-12-02     106.
10 Slice1 training FB     2015-12-03     104.
# ℹ 1,802 more rows
# ℹ Use `print(n = ...)` to see more rows

# 可视化
> resample_spec %>%
+   tk_time_series_cv_plan() %>%
+   plot_time_series_cv_plan(
+     date, adjusted, # date variable and value variable
+     # Additional arguments passed to plot_time_series(),
+     .facet_ncol = 2,
+     .line_alpha = 0.5,
+     .interactive = FALSE
+   )

timetk体验起来是非常棒的,短信几行代码,就把如何观察时间序列数据看得明明白白。

转载请注明出处:
http://blog.fens.me/r-timetk/

plot.xts时间序列可视化

R的极客理想系列文章,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,互联网….都在使用R语言。

要成为有理想的极客,我们不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让我们一起动起来吧,开始R的极客理想。

关于作者:

  • 张丹(Conan), 程序员Java,R,PHP,Javascript
  • weibo:@Conan_Z
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-xts-xtsextra/

xtsExtra-r

前言

r-bloggers的一篇博文,让我有动力继续发现xts的强大。xts扩展了zoo的基础数据结构,并提供了更丰富的功能函数。xtsExtra补充库,从可视化的角度出发,提供了一个简单而效果非凡的作图函数plot.xts。

本文将用plot.xts来演示,xts对象的时间序列可视化!

目录

  1. xtsExtra介绍
  2. xtsExtra安装
  3. plot.xts的使用

1. xtsExtra介绍

xtsExtra是xts包的功能补充包,该软件包在Google Summer of Code 2012被开发,最终将合并到xts包。xtsExtra提供的主要功能就是plot.xts。

注:我发现xts::plot.xts的函数,与xtsExtra::plot.xts还是有差别的。

关于xts包的介绍,请参考文章:可扩展的时间序列xts

下面我们安装xtsExtra包。

2. xtsExtra安装

由于xtsExtra没有发布到CRAN,我们要从R-Forge下载。


~ R

> install.packages("xtsExtra", repos="http://R-Forge.R-project.org")
Warning in install.packages :
  package ‘xtsExtra’ is not available (for R version 3.0.1)
trying URL 'http://R-Forge.R-project.org/bin/windows/contrib/3.0/xtsExtra_0.0-1.zip'
Content type 'application/zip' length 242682 bytes (236 Kb)
opened URL
downloaded 236 Kb

package ‘xtsExtra’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
	C:\Users\Administrator\AppData\Local\Temp\Rtmp04stLd\downloaded_packages

加载xtsExtra


> library(xtsExtra)
载入需要的程辑包:zoo

载入程辑包:‘zoo’

下列对象被屏蔽了from ‘package:base’:

    as.Date, as.Date.numeric

载入需要的程辑包:xts

载入程辑包:‘xtsExtra’

下列对象被屏蔽了from ‘package:xts’:

    plot.xts

Warning messages:
1: 程辑包‘zoo’是用R版本3.0.2 来建造的 
2: 程辑包‘xts’是用R版本3.0.2 来建造的 

plot.xts函数被用来,覆盖xts::plot.xts函数。

3. plot.xts的使用

  • 1). plot.xts的参数列表
  • 2). 简单的时间序列
  • 3). K线图
  • 4). panel配置
  • 5). screens配置
  • 6). events配置
  • 7). 双时间序列
  • 9). barplot

1). plot.xts的参数列表


> names(formals(plot.xts))
 [1] "x"              "y"              "screens"        "layout.screens" "..."           
 [6] "yax.loc"        "auto.grid"      "major.ticks"    "minor.ticks"    "major.format"  
[11] "bar.col.up"     "bar.col.dn"     "candle.col"     "xy.labels"      "xy.lines"      
[16] "ylim"           "panel"          "auto.legend"    "legend.names"   "legend.loc"    
[21] "legend.pars"    "events"         "blocks"         "nc"             "nr"       

2). 简单的时间序列


> data(sample_matrix)
> sample_xts <- as.xts(sample_matrix)
> plot(sample_xts[,1]) 
> class(sample_xts[,1])
[1] "xts" "zoo"

plot.xts-basic

3). K线图

红白色


> plot(sample_xts[1:30, ], type = "candles")

plot.xts-red

自定义颜色


> plot(sample_xts[1:30, ], type = "candles", bar.col.up = "blue", bar.col.dn = "violet", candle.col = "green4")

plot.xts-col

4). panel配置
基本面板


> plot(sample_xts[,1:2]) 

plot.xts-panel

多行面板


> plot(sample_xts[,rep(1:4, each = 3)]) 

plot.xts-mpanel

自由组合面板


> plot(sample_xts[,1:4], layout.screens = matrix(c(1,1,1,1,2,3,4,4),ncol = 2, byrow = TRUE))

plot.xts-cpanel

5). screens配置

双屏幕显示,每屏幕2条线


> plot(sample_xts, screens = 1:2) 

plot.xts-screen

双屏幕显示,指定曲线出现的屏幕和颜色


> plot(sample_xts, screens = c(1,2,1,2), col = c(1,3,2,2))

plot.xts-screen-4

双屏幕显示,指定不同的坐标系


> plot(10^sample_xts, screens = 1:2, log= c("","y"))

plot.xts-screen-axis

双屏幕显示,指定不同的输出图形


> plot(sample_xts[1:75,1:2] - 50.5, type = c("l","h"), lwd = c(1,2))

plot.xts-screen-chart

多屏幕,分组显示


> plot(sample_xts[,c(1:4, 3:4)], layout = matrix(c(1,1,1,1,2,2,3,4,5,6), ncol = 2, byrow = TRUE), yax.loc = "left")

plot.xts-screen-mg

6). events配置

基本事件分割线


> plot(sample_xts[,1], events = list(time = c("2007-03-15","2007-05-01"), label = "bad days"), blocks = list(start.time = c("2007-03-05", "2007-04-15"), end.time = c("2007-03-20","2007-05-30"), col = c("lightblue1", "lightgreen")))

plot.xts-event

7). 双时间序列

双坐标视图


> plot(sample_xts[,1],sample_xts[,2])

plot.xts-scatter

双坐标梯度视图


> cr <- colorRampPalette(c("#00FF00","#FF0000"))
> plot(sample_xts[,1],sample_xts[,2], xy.labels = FALSE, xy.lines = TRUE, col = cr(NROW(sample_xts)), type = "l")

plot.xts-gradient

8). xts类型转换作图
ts类型作图


> tser <- ts(cumsum(rnorm(50, 0.05, 0.15)), start = 2007, frequency = 12)
> class(tser)
[1] "ts"
> plot(tser)

plot.xts-ts

以xts类型作图


> plot.xts(tser)

plot.xts-ts-xts

9). barplot


> x <- xts(matrix(abs(rnorm(72)), ncol = 6), Sys.Date() + 1:12)
> colnames(x) <- LETTERS[1:6]
> barplot(x)

plot.xts-barplot

我们看到xtsExtra::plot.xts提供了强大的作图功能,很容易做出可视的时间序列!

转载请注明出处:
http://blog.fens.me/r-xts-xtsextra/

打赏作者

可扩展的时间序列xts

R的极客理想系列文章,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,互联网….都在使用R语言。

要成为有理想的极客,我们不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让我们一起动起来吧,开始R的极客理想。

关于作者:

  • 张丹(Conan), 程序员Java,R,PHP,Javascript
  • weibo:@Conan_Z
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-xts/

r-xts

前言

本文是继R语言zoo时间序列基础库的扩展实现。看上去简单的时间序列,内藏复杂的规律。zoo作为时间序列的基础库,是面向通用的设计,可以用来定义股票数据,也可以分析天气数据。但由于业务行为的不同,我们需要更多的辅助函数,来帮助我们更高效的完成任务。

xts扩展了zoo,提供更多的数据处理和数据变换的函数。

目录

  1. xts介绍
  2. xts安装
  3. xts数据结构
  4. xts的API介绍
  5. xts使用

1. xts介绍

xts是对时间序列数据(zoo)的一种扩展实现,目标是为了统一时间序列的操作接口。实际上,xts类型继承了zoo类型,丰富了时间序列数据处理的函数,API定义更贴近使用者,更实用,更简单!

xts项目地址:http://r-forge.r-project.org/projects/xts/

2. xts安装

系统环境

  • Win7 64bit
  • R: 3.0.1 x86_64-w64-mingw32/x64 b4bit

xts安装


> install.packages("xts")
also installing the dependency ‘zoo’

trying URL 'http://mirror.bjtu.edu.cn/cran/bin/windows/contrib/3.0/zoo_1.7-10.zip'
Content type 'application/zip' length 875046 bytes (854 Kb)
opened URL
downloaded 854 Kb

trying URL 'http://mirror.bjtu.edu.cn/cran/bin/windows/contrib/3.0/xts_0.9-7.zip'
Content type 'application/zip' length 661664 bytes (646 Kb)
opened URL
downloaded 646 Kb

package ‘zoo’ successfully unpacked and MD5 sums checked
package ‘xts’ successfully unpacked and MD5 sums checked

3. xts数据结构

xts-structure

xts扩展zoo的基础结构,由3部分组合。

  • 索引部分:时间类型向量
  • 数据部分:以矩阵为基础类型,支持可以与矩阵相互转换的任何类型
  • 属性部分:附件信息,包括时区,索引时间类型的格式等

4. xts的API介绍

xts基础

  • xts: 定义xts数据类型,继承zoo类型
  • coredata.xts: 对xts部分数据赋值
  • xtsAttributes: xts对象属性赋值
  • [.xts: 用[]语法,取数据子集
  • dimnames.xts: xts维度名赋值
  • sample_matrix: 测试数据集,包括180条xts对象的记录,matrix类型
  • xtsAPI: C语言API接口

类型转换

  • as.xts: 转换对象到xts(zoo)类型
  • as.xts.methods: 转换对象到xts函数
  • plot.xts: 为plot函数,提供xts的接口作图
  • .parseISO8601: 把字符串(ISO8601格式)输出为,POSIXct类型的,包括开始时间和结束时间的list对象
  • firstof: 创建一个开始时间,POSIXct类型
  • lastof: 创建一个结束时间,POSIXct类型
  • indexClass: 取索引类型
  • .indexDate: 取索引的
  • .indexday: 索引的日值
  • .indexyday: 索引的年(日)值
  • .indexmday: 索引的月(日)值
  • .indexwday: 索引的周(日)值
  • .indexweek: 索引的周值
  • .indexmon: 索引的月值
  • .indexyear: 索引的年值
  • .indexhour: 索引的时值
  • .indexmin: 索引的分值
  • .indexsec: 索引的秒值

数据处理

  • align.time: 以下一个时间对齐数据,秒,分钟,小时
  • endpoints: 按时间单元提取索引数据
  • merge.xts: 合并多个xts对象,重写zoo::merge.zoo函数
  • rbind.xts: 数据按行合并,为rbind函数,提供xts的接口
  • split.xts: 数据分隔,为split函数,提供xts的接口
  • na.locf.xts: 替换NA值,重写zoo:na.locf函数

数据统计

  • apply.daily: 按日分割数据,执行函数
  • apply.weekly: 按周分割数据,执行函数
  • apply.monthly: 按月分割数据,执行函数
  • apply.quarterly: 按季分割数据,执行函数
  • apply.yearly: 按年分割数据,执行函数
  • to.period: 按期间分割数据
  • period.apply: 按期间执行自定义函数
  • period.max: 按期间计算最大值
  • period.min: 按期间计算最小值
  • period.prod: 按期间计算指数
  • period.sum: 按期间求和
  • nseconds: 计算数据集,包括多少秒
  • nminutes: 计算数据集,包括多少分
  • nhours: 计算数据集,包括多少时
  • ndays: 计算数据集,包括多少日
  • nweeks: 计算数据集,包括多少周
  • nmonths: 计算数据集,包括多少月
  • nquarters: 计算数据集,包括多少季
  • nyears: 计算数据集,包括多少年
  • periodicity: 查看时间序列的期间

辅助工具

  • first: 从开始到结束,设置条件取子集
  • last: 从结束到开始,设置条件取子集
  • timeBased: 判断是否是时间类型
  • timeBasedSeq: 创建时间的序列
  • diff.xts: 计算步长和差分
  • isOrdered: 检查向量是否是顺序的
  • make.index.unique: 强制时间唯一,增加毫秒随机数
  • axTicksByTime: 计算X轴刻度标记位置按时间描述
  • indexTZ: 查询xts对象的时区

5. xts使用

  • 1). xts类型基本操作
  • 2). xts的作图
  • 3). xts类型转换
  • 4). xts数据处理
  • 5). xts数据统计计算
  • 6). xts时间序列工具使用

1). xts类型基本操作

测试数据集sample_matrix


> library(xts)
> data(sample_matrix)
> head(sample_matrix)
               Open     High      Low    Close
2007-01-02 50.03978 50.11778 49.95041 50.11778
2007-01-03 50.23050 50.42188 50.23050 50.39767
2007-01-04 50.42096 50.42096 50.26414 50.33236
2007-01-05 50.37347 50.37347 50.22103 50.33459
2007-01-06 50.24433 50.24433 50.11121 50.18112
2007-01-07 50.13211 50.21561 49.99185 49.99185

定义xts类型对象


> sample.xts <- as.xts(sample_matrix, descr='my new xts object')
> class(sample.xts)
[1] "xts" "zoo"

> str(sample.xts)
An ‘xts’ object on 2007-01-02/2007-06-30 containing:
  Data: num [1:180, 1:4] 50 50.2 50.4 50.4 50.2 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : chr [1:4] "Open" "High" "Low" "Close"
  Indexed by objects of class: [POSIXct,POSIXt] TZ: 
  xts Attributes:  
List of 1
 $ descr: chr "my new xts object"

> head(sample.xts)
               Open     High      Low    Close
2007-01-02 50.03978 50.11778 49.95041 50.11778
2007-01-03 50.23050 50.42188 50.23050 50.39767
2007-01-04 50.42096 50.42096 50.26414 50.33236
2007-01-05 50.37347 50.37347 50.22103 50.33459
2007-01-06 50.24433 50.24433 50.11121 50.18112
2007-01-07 50.13211 50.21561 49.99185 49.99185

> attr(sample.xts,'descr')
[1] "my new xts object"

xts数据查询


> head(sample.xts['2007'])
               Open     High      Low    Close
2007-01-02 50.03978 50.11778 49.95041 50.11778
2007-01-03 50.23050 50.42188 50.23050 50.39767
2007-01-04 50.42096 50.42096 50.26414 50.33236
2007-01-05 50.37347 50.37347 50.22103 50.33459
2007-01-06 50.24433 50.24433 50.11121 50.18112
2007-01-07 50.13211 50.21561 49.99185 49.99185

> head(sample.xts['2007-03/'])
               Open     High      Low    Close
2007-03-01 50.81620 50.81620 50.56451 50.57075
2007-03-02 50.60980 50.72061 50.50808 50.61559
2007-03-03 50.73241 50.73241 50.40929 50.41033
2007-03-04 50.39273 50.40881 50.24922 50.32636
2007-03-05 50.26501 50.34050 50.26501 50.29567
2007-03-06 50.27464 50.32019 50.16380 50.16380

> head(sample.xts['2007-03-06/2007'])
               Open     High      Low    Close
2007-03-06 50.27464 50.32019 50.16380 50.16380
2007-03-07 50.14458 50.20278 49.91381 49.91381
2007-03-08 49.93149 50.00364 49.84893 49.91839
2007-03-09 49.92377 49.92377 49.74242 49.80712
2007-03-10 49.79370 49.88984 49.70385 49.88698
2007-03-11 49.83062 49.88295 49.76031 49.78806

> sample.xts['2007-01-03']
              Open     High     Low    Close
2007-01-03 50.2305 50.42188 50.2305 50.39767

2). 操作xts的作图

曲线图


> data(sample_matrix)
> plot(sample_matrix)

> plot(as.xts(sample_matrix))
Warning message:
In plot.xts(as.xts(sample_matrix)) :
  only the univariate series will be plotted

plot-line

K线图


> plot(as.xts(sample_matrix), type='candles')

plot-candles

3). xts类型转换

分别创建首尾时间:firstof, lastof


> firstof(2000)
[1] "2000-01-01 CST"

> firstof(2005,01,01)
[1] "2005-01-01 CST"

> lastof(2007)
[1] "2007-12-31 23:59:59.99998 CST"

> lastof(2007,10)
[1] "2007-10-31 23:59:59.99998 CST"

创建首尾时间


> .parseISO8601('2000')
$first.time
[1] "2000-01-01 CST"

$last.time
[1] "2000-12-31 23:59:59.99998 CST"

> .parseISO8601('2000-05/2001-02')
$first.time
[1] "2000-05-01 CST"

$last.time
[1] "2001-02-28 23:59:59.99998 CST"

> .parseISO8601('2000-01/02')
$first.time
[1] "2000-01-01 CST"

$last.time
[1] "2000-02-29 23:59:59.99998 CST"

> .parseISO8601('T08:30/T15:00')
$first.time
[1] "1970-01-01 08:30:00 CST"

$last.time
[1] "1970-12-31 15:00:59.99999 CST"

取索引类型


> x <- timeBasedSeq('2010-01-01/2010-01-02 12:00')
> x <- xts(1:length(x), x)

> head(x)
                    [,1]
2010-01-01 00:00:00    1
2010-01-01 00:01:00    2
2010-01-01 00:02:00    3
2010-01-01 00:03:00    4
2010-01-01 00:04:00    5
2010-01-01 00:05:00    6

> indexClass(x)
[1] "POSIXt"  "POSIXct"

索引时间格式化


> indexFormat(x) <- "%Y-%b-%d %H:%M:%OS3"
> head(x)
                          [,1]
2010-一月-01 00:00:00.000    1
2010-一月-01 00:01:00.000    2
2010-一月-01 00:02:00.000    3
2010-一月-01 00:03:00.000    4
2010-一月-01 00:04:00.000    5
2010-一月-01 00:05:00.000    6

取索引时间


> .indexhour(head(x))
[1] 0 0 0 0 0 0

> .indexmin(head(x))
[1] 0 1 2 3 4 5

4). xts数据处理
数据对齐


> x <- Sys.time() + 1:30

#整10秒对齐
> align.time(x, 10)
 [1] "2013-11-18 15:42:30 CST" "2013-11-18 15:42:30 CST"
 [3] "2013-11-18 15:42:30 CST" "2013-11-18 15:42:40 CST"
 [5] "2013-11-18 15:42:40 CST" "2013-11-18 15:42:40 CST"
 [7] "2013-11-18 15:42:40 CST" "2013-11-18 15:42:40 CST"
 [9] "2013-11-18 15:42:40 CST" "2013-11-18 15:42:40 CST"
[11] "2013-11-18 15:42:40 CST" "2013-11-18 15:42:40 CST"
[13] "2013-11-18 15:42:40 CST" "2013-11-18 15:42:50 CST"
[15] "2013-11-18 15:42:50 CST" "2013-11-18 15:42:50 CST"
[17] "2013-11-18 15:42:50 CST" "2013-11-18 15:42:50 CST"
[19] "2013-11-18 15:42:50 CST" "2013-11-18 15:42:50 CST"
[21] "2013-11-18 15:42:50 CST" "2013-11-18 15:42:50 CST"
[23] "2013-11-18 15:42:50 CST" "2013-11-18 15:43:00 CST"
[25] "2013-11-18 15:43:00 CST" "2013-11-18 15:43:00 CST"
[27] "2013-11-18 15:43:00 CST" "2013-11-18 15:43:00 CST"
[29] "2013-11-18 15:43:00 CST" "2013-11-18 15:43:00 CST"

#整60秒对齐
> align.time(x, 60)
 [1] "2013-11-18 15:43:00 CST" "2013-11-18 15:43:00 CST"
 [3] "2013-11-18 15:43:00 CST" "2013-11-18 15:43:00 CST"
 [5] "2013-11-18 15:43:00 CST" "2013-11-18 15:43:00 CST"
 [7] "2013-11-18 15:43:00 CST" "2013-11-18 15:43:00 CST"
 [9] "2013-11-18 15:43:00 CST" "2013-11-18 15:43:00 CST"
[11] "2013-11-18 15:43:00 CST" "2013-11-18 15:43:00 CST"
[13] "2013-11-18 15:43:00 CST" "2013-11-18 15:43:00 CST"
[15] "2013-11-18 15:43:00 CST" "2013-11-18 15:43:00 CST"
[17] "2013-11-18 15:43:00 CST" "2013-11-18 15:43:00 CST"
[19] "2013-11-18 15:43:00 CST" "2013-11-18 15:43:00 CST"
[21] "2013-11-18 15:43:00 CST" "2013-11-18 15:43:00 CST"
[23] "2013-11-18 15:43:00 CST" "2013-11-18 15:43:00 CST"
[25] "2013-11-18 15:43:00 CST" "2013-11-18 15:43:00 CST"
[27] "2013-11-18 15:43:00 CST" "2013-11-18 15:43:00 CST"
[29] "2013-11-18 15:43:00 CST" "2013-11-18 15:43:00 CST"

按时间分割数据,并计算


> xts.ts <- xts(rnorm(231),as.Date(13514:13744,origin="1970-01-01"))
> apply.monthly(xts.ts,mean)
                  [,1]
2007-01-31  0.17699984
2007-02-28  0.30734220
2007-03-31 -0.08757189
2007-04-30  0.18734688
2007-05-31  0.04496954
2007-06-30  0.06884836
2007-07-31  0.25081814
2007-08-19 -0.28845938

> apply.monthly(xts.ts,function(x) var(x))
                [,1]
2007-01-31 0.9533217
2007-02-28 0.9158947
2007-03-31 1.2821450
2007-04-30 1.2805976
2007-05-31 0.9725438
2007-06-30 1.5228904
2007-07-31 0.8737030
2007-08-19 0.8490521

> apply.quarterly(xts.ts,mean)
                 [,1]
2007-03-31 0.12642053
2007-06-30 0.09977926
2007-08-19 0.04589268

> apply.yearly(xts.ts,mean)
                 [,1]
2007-08-19 0.09849522

按期间分隔:to.period


> data(sample_matrix)
> to.period(sample_matrix)
           sample_matrix.Open sample_matrix.High sample_matrix.Low sample_matrix.Close
2007-01-31           50.03978           50.77336          49.76308            50.22578
2007-02-28           50.22448           51.32342          50.19101            50.77091
2007-03-31           50.81620           50.81620          48.23648            48.97490
2007-04-30           48.94407           50.33781          48.80962            49.33974
2007-05-31           49.34572           49.69097          47.51796            47.73780
2007-06-30           47.74432           47.94127          47.09144            47.76719
> class(to.period(sample_matrix))
[1] "matrix"

> samplexts <- as.xts(sample_matrix)
> to.period(samplexts)
           samplexts.Open samplexts.High samplexts.Low samplexts.Close
2007-01-31       50.03978       50.77336      49.76308        50.22578
2007-02-28       50.22448       51.32342      50.19101        50.77091
2007-03-31       50.81620       50.81620      48.23648        48.97490
2007-04-30       48.94407       50.33781      48.80962        49.33974
2007-05-31       49.34572       49.69097      47.51796        47.73780
2007-06-30       47.74432       47.94127      47.09144        47.76719
> class(to.period(samplexts))
[1] "xts" "zoo"

按期间分割索引数据


> data(sample_matrix)

> endpoints(sample_matrix)
[1]   0  30  58  89 119 150 180

> endpoints(sample_matrix, 'days',k=7)
 [1]   0   6  13  20  27  34  41  48  55  62  69  76  83  90  97 104 111 118 125
[20] 132 139 146 153 160 167 174 180

> endpoints(sample_matrix, 'weeks')
 [1]   0   7  14  21  28  35  42  49  56  63  70  77  84  91  98 105 112 119 126
[20] 133 140 147 154 161 168 175 180

> endpoints(sample_matrix, 'months')
[1]   0  30  58  89 119 150 180

数据合并:按列合并


> (x <- xts(4:10, Sys.Date()+4:10))
           [,1]
2013-11-22    4
2013-11-23    5
2013-11-24    6
2013-11-25    7
2013-11-26    8
2013-11-27    9
2013-11-28   10

> (y <- xts(1:6, Sys.Date()+1:6))
           [,1]
2013-11-19    1
2013-11-20    2
2013-11-21    3
2013-11-22    4
2013-11-23    5
2013-11-24    6

> merge(x,y)
            x  y
2013-11-19 NA  1
2013-11-20 NA  2
2013-11-21 NA  3
2013-11-22  4  4
2013-11-23  5  5
2013-11-24  6  6
2013-11-25  7 NA
2013-11-26  8 NA
2013-11-27  9 NA
2013-11-28 10 NA

#取索引将领合并
> merge(x,y, join='inner')
           x y
2013-11-22 4 4
2013-11-23 5 5
2013-11-24 6 6

#以左侧为基础合并
> merge(x,y, join='left')
            x  y
2013-11-22  4  4
2013-11-23  5  5
2013-11-24  6  6
2013-11-25  7 NA
2013-11-26  8 NA
2013-11-27  9 NA
2013-11-28 10 NA

数据合并:按行合并


> x <- xts(1:3, Sys.Date()+1:3)

> rbind(x,x)
           [,1]
2013-11-19    1
2013-11-19    1
2013-11-20    2
2013-11-20    2
2013-11-21    3
2013-11-21    3

数据切片:按行切片


> data(sample_matrix)
> x <- as.xts(sample_matrix)

按月切片
> split(x)[[1]]
               Open     High      Low    Close
2007-01-02 50.03978 50.11778 49.95041 50.11778
2007-01-03 50.23050 50.42188 50.23050 50.39767
2007-01-04 50.42096 50.42096 50.26414 50.33236
2007-01-05 50.37347 50.37347 50.22103 50.33459
2007-01-06 50.24433 50.24433 50.11121 50.18112
2007-01-07 50.13211 50.21561 49.99185 49.99185
2007-01-08 50.03555 50.10363 49.96971 49.98806
2007-01-09 49.99489 49.99489 49.80454 49.91333
2007-01-10 49.91228 50.13053 49.91228 49.97246
2007-01-11 49.88529 50.23910 49.88529 50.23910
2007-01-12 50.21258 50.35980 50.17176 50.28519
2007-01-13 50.32385 50.48000 50.32385 50.41286
2007-01-14 50.46359 50.62395 50.46359 50.60145
2007-01-15 50.61724 50.68583 50.47359 50.48912
2007-01-16 50.62024 50.73731 50.56627 50.67835
2007-01-17 50.74150 50.77336 50.44932 50.48644
2007-01-18 50.48051 50.60712 50.40269 50.57632
2007-01-19 50.41381 50.55627 50.41278 50.41278
2007-01-20 50.35323 50.35323 50.02142 50.02142
2007-01-21 50.16188 50.42090 50.16044 50.42090
2007-01-22 50.36008 50.43875 50.21129 50.21129
2007-01-23 50.03966 50.16961 50.03670 50.16961
2007-01-24 50.10953 50.26942 50.06387 50.23145
2007-01-25 50.20738 50.28268 50.12913 50.24334
2007-01-26 50.16008 50.16008 49.94052 50.07024
2007-01-27 50.06041 50.09777 49.97267 50.01091
2007-01-28 49.96586 50.00217 49.87468 49.88096
2007-01-29 49.85624 49.93038 49.76308 49.91875
2007-01-30 49.85477 50.02180 49.77242 50.02180
2007-01-31 50.07049 50.22578 50.07049 50.22578

按周切片
> split(x, f="weeks")[[1]]
               Open     High      Low    Close
2007-01-02 50.03978 50.11778 49.95041 50.11778
2007-01-03 50.23050 50.42188 50.23050 50.39767
2007-01-04 50.42096 50.42096 50.26414 50.33236
2007-01-05 50.37347 50.37347 50.22103 50.33459
2007-01-06 50.24433 50.24433 50.11121 50.18112
2007-01-07 50.13211 50.21561 49.99185 49.99185
2007-01-08 50.03555 50.10363 49.96971 49.98806
> split(x, f="weeks")[[2]]
               Open     High      Low    Close
2007-01-09 49.99489 49.99489 49.80454 49.91333
2007-01-10 49.91228 50.13053 49.91228 49.97246
2007-01-11 49.88529 50.23910 49.88529 50.23910
2007-01-12 50.21258 50.35980 50.17176 50.28519
2007-01-13 50.32385 50.48000 50.32385 50.41286
2007-01-14 50.46359 50.62395 50.46359 50.60145
2007-01-15 50.61724 50.68583 50.47359 50.48912

NA值处理


> x <- xts(1:10, Sys.Date()+1:10)
> x[c(1,2,5,9,10)] <- NA
> x
           [,1]
2013-11-19   NA
2013-11-20   NA
2013-11-21    3
2013-11-22    4
2013-11-23   NA
2013-11-24    6
2013-11-25    7
2013-11-26    8
2013-11-27   NA
2013-11-28   NA

#取前一个
> na.locf(x)
           [,1]
2013-11-19   NA
2013-11-20   NA
2013-11-21    3
2013-11-22    4
2013-11-23    4
2013-11-24    6
2013-11-25    7
2013-11-26    8
2013-11-27    8
2013-11-28    8

#取后一个
> na.locf(x, fromLast=TRUE)
           [,1]
2013-11-19    3
2013-11-20    3
2013-11-21    3
2013-11-22    4
2013-11-23    6
2013-11-24    6
2013-11-25    7
2013-11-26    8
2013-11-27   NA
2013-11-28   NA

5). xts数据统计计算

取开始时间,结束时间


> xts.ts <- xts(rnorm(231),as.Date(13514:13744,origin="1970-01-01"))

> start(xts.ts)
[1] "2007-01-01"
> end(xts.ts)
[1] "2007-08-19"

> periodicity(xts.ts)
Daily periodicity from 2007-01-01 to 2007-08-19 

计算时间区间


> data(sample_matrix)
> ndays(sample_matrix)
[1] 180
> nweeks(sample_matrix)
[1] 26
> nmonths(sample_matrix)
[1] 6
> nquarters(sample_matrix)
[1] 2
> nyears(sample_matrix)
[1] 1

按期间计算统计指标


> zoo.data <- zoo(rnorm(31)+10,as.Date(13514:13744,origin="1970-01-01"))

#按周获得期间
> ep <- endpoints(zoo.data,'weeks')
> ep
 [1]   0   7  14  21  28  35  42  49  56  63  70  77  84  91  98 105 112 119
[19] 126 133 140 147 154 161 168 175 182 189 196 203 210 217 224 231

#计算周的均值
> period.apply(zoo.data, INDEX=ep, FUN=function(x) mean(x))
2007-01-07 2007-01-14 2007-01-21 2007-01-28 2007-02-04 2007-02-11 2007-02-18 
 10.200488   9.649387  10.304151   9.864847  10.382943   9.660175   9.857894 
2007-02-25 2007-03-04 2007-03-11 2007-03-18 2007-03-25 2007-04-01 2007-04-08 
 10.495037   9.569531  10.292899   9.651616  10.089103   9.961048  10.304860 
2007-04-15 2007-04-22 2007-04-29 2007-05-06 2007-05-13 2007-05-20 2007-05-27 
  9.658432   9.887531  10.608082   9.747787  10.052955   9.625730  10.430030 
2007-06-03 2007-06-10 2007-06-17 2007-06-24 2007-07-01 2007-07-08 2007-07-15 
  9.814703  10.224869   9.509881  10.187905  10.229310  10.261725   9.855776 
2007-07-22 2007-07-29 2007-08-05 2007-08-12 2007-08-19 
  9.445072  10.482020   9.844531  10.200488   9.649387 

#计算周的最大值
> head(period.max(zoo.data, INDEX=ep))
               [,1]
2007-01-07 12.05912
2007-01-14 10.79286
2007-01-21 11.60658
2007-01-28 11.63455
2007-02-04 12.05912
2007-02-11 10.67887

#计算周的最小值
> head(period.min(zoo.data, INDEX=ep))
               [,1]
2007-01-07 8.874509
2007-01-14 8.534655
2007-01-21 9.069773
2007-01-28 8.461555
2007-02-04 9.421085
2007-02-11 8.534655

#计算周的一个指数值
> head(period.prod(zoo.data, INDEX=ep))
               [,1]
2007-01-07 11140398
2007-01-14  7582350
2007-01-21 11930334
2007-01-28  8658933
2007-02-04 12702505
2007-02-11  7702767

6). xts时间序列工具使用

检查时间类型


> timeBased(Sys.time())
[1] TRUE
> timeBased(Sys.Date())
[1] TRUE
> timeBased(200701)
[1] FALSE

创建时间序列


#按年
> timeBasedSeq('1999/2008')
 [1] "1999-01-01" "2000-01-01" "2001-01-01" "2002-01-01" "2003-01-01"
 [6] "2004-01-01" "2005-01-01" "2006-01-01" "2007-01-01" "2008-01-01"

#按月
> head(timeBasedSeq('199901/2008'))
[1] "十二月 1998" "一月 1999"   "二月 1999"   "三月 1999"   "四月 1999"  
[6] "五月 1999" 

#按日
> head(timeBasedSeq('199901/2008/d'),40)
 [1] "十二月 1998" "一月 1999"   "一月 1999"   "一月 1999"   "一月 1999"  
 [6] "一月 1999"   "一月 1999"   "一月 1999"   "一月 1999"   "一月 1999"  
[11] "一月 1999"   "一月 1999"   "一月 1999"   "一月 1999"   "一月 1999"  
[16] "一月 1999"   "一月 1999"   "一月 1999"   "一月 1999"   "一月 1999"  
[21] "一月 1999"   "一月 1999"   "一月 1999"   "一月 1999"   "一月 1999"  
[26] "一月 1999"   "一月 1999"   "一月 1999"   "一月 1999"   "一月 1999"  
[31] "一月 1999"   "一月 1999"   "二月 1999"   "二月 1999"   "二月 1999"  
[36] "二月 1999"   "二月 1999"   "二月 1999"   "二月 1999"   "二月 1999" 

#按数量创建,100分钟的数据集
> timeBasedSeq('20080101 0830',length=100)
$from
[1] "2008-01-01 08:30:00 CST"
$to
[1] NA
$by
[1] "mins"
$length.out
[1] 100

按索引取数据first, last


> x <- xts(1:100, Sys.Date()+1:100)

> head(x)
           [,1]
2013-11-19    1
2013-11-20    2
2013-11-21    3
2013-11-22    4
2013-11-23    5
2013-11-24    6

> first(x, 10)
           [,1]
2013-11-19    1
2013-11-20    2
2013-11-21    3
2013-11-22    4
2013-11-23    5
2013-11-24    6
2013-11-25    7
2013-11-26    8
2013-11-27    9
2013-11-28   10

> first(x, '1 day')
           [,1]
2013-11-19    1

> last(x, '1 weeks')
           [,1]
2014-02-24   98
2014-02-25   99
2014-02-26  100

计算步长和差分


> x <- xts(1:5, Sys.Date()+1:5)
#正向
> lag(x)
           [,1]
2013-11-19   NA
2013-11-20    1
2013-11-21    2
2013-11-22    3
2013-11-23    4

#反向
> lag(x, k=-1, na.pad=FALSE) 
           [,1]
2013-11-19    2
2013-11-20    3
2013-11-21    4
2013-11-22    5

#1阶差分
> diff(x)
           [,1]
2013-11-19   NA
2013-11-20    1
2013-11-21    1
2013-11-22    1
2013-11-23    1

#2阶差分
> diff(x, lag=2)
           [,1]
2013-11-19   NA
2013-11-20   NA
2013-11-21    2
2013-11-22    2
2013-11-23    2

检查向量是否排序好的


> isOrdered(1:10, increasing=TRUE)
[1] TRUE

> isOrdered(1:10, increasing=FALSE)
[1] FALSE

> isOrdered(c(1,1:10), increasing=TRUE)
[1] FALSE

> isOrdered(c(1,1:10), increasing=TRUE, strictly=FALSE)
[1] TRUE

强制唯一索引


> x <- xts(1:5, as.POSIXct("2011-01-21") + c(1,1,1,2,3)/1e3)
> x
                        [,1]
2011-01-21 00:00:00.000    1
2011-01-21 00:00:00.000    2
2011-01-21 00:00:00.000    3
2011-01-21 00:00:00.002    4
2011-01-21 00:00:00.003    5

> make.index.unique(x)
                           [,1]
2011-01-21 00:00:00.000999    1
2011-01-21 00:00:00.001000    2
2011-01-21 00:00:00.001001    3
2011-01-21 00:00:00.002000    4
2011-01-21 00:00:00.003000    5

查询xts对象时区


> x <- xts(1:10, Sys.Date()+1:10)

> indexTZ(x)
[1] "UTC"
> tzone(x)
[1] "UTC"

> str(x)
An ‘xts’ object on 2013-11-19/2013-11-28 containing:
  Data: int [1:10, 1] 1 2 3 4 5 6 7 8 9 10
  Indexed by objects of class: [Date] TZ: UTC
  xts Attributes:  
 NULL

xts给了zoo类型时间序列更多的API支持,这样我们就有了更方便的工具,可以做各种的时间序列的转换和变形了。

转载请注明出处:
http://blog.fens.me/r-xts/

打赏作者

R语言时间序列基础库zoo

R的极客理想系列文章,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,互联网….都在使用R语言。

要成为有理想的极客,我们不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让我们一起动起来吧,开始R的极客理想。

关于作者:

  • 张丹(Conan), 程序员Java,R,PHP,Javascript
  • weibo:@Conan_Z
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-zoo/

zoo-r

前言

时间序列分析是一种动态数据处理的统计方法,通过对时间序列数据的分析,我们可以感觉到世界正改变着什么!R语言作为统计分析的利器,对时间序列处理有着强大的支持。在R语言中,单独为时间序列数据定义了一种数据类型zoo,zoo是时间序列的基础,也是股票分析的基础。

本文将介绍zoo库在R语言中的结构和使用。

目录

  1. zoo介绍
  2. zoo安装
  3. zoo的API介绍
  4. zoo使用

1. zoo介绍

zoo是一个R语言类库,zoo类库中定义了一个名为zoo的S3类型对象,用于描述规则的和不规则的有序的时间序列数据。zoo对象是一个独立的对象,包括索引、日期、时间,只依赖于基础的R环境,zooreg对象继承了zoo对象,只能用于规则的的时间序列数据。

R语言的其他程序包,都是以zoo, zooreg为时间序列数据的基础!

zoo的项目主页:http://zoo.r-forge.r-project.org/

2. zoo安装

系统环境

  • Win7 64bit
  • R: 3.0.1 x86_64-w64-mingw32/x64 b4bit

zoo安装


~ R
> install.packages("zoo")
> library(zoo)

3. zoo的API介绍

基础对象

  • zoo: 有序的时间序列对象
  • zooreg: 规则的的时间序列对象,继承zoo对象

类型转换

  • as.zoo: 把一个对象转型为zoo类型
  • plot.zoo: 为plot函数,提供zoo的接口
  • xyplot.zoo: 为lattice的xyplot函数,提供zoo的接口
  • ggplot2.zoo: 为ggplot2包,提供zoo的接口

数据操作

  • coredata: 获得和修改zoo的数据部分
  • index: 获得和修改zoo的索引部分
  • window.zoo: 按时间过滤数据
  • merge.zoo: 合并多个zoo对象
  • read.zoo: 从文件读写zoo序列
  • aggregate.zoo: 计算zoo数据
  • rollapply: 对zoo数据的滚动处理
  • rollmean: 对zoo数据的滚动,计算均值

NA值处理

  • na.fill: NA值的填充
  • na.locf: 替换NA值
  • na.aggregate: 计算统计值替换NA值
  • na.approx: 计算插值替换NA值
  • na.StructTS: 计算seasonal Kalman filter替换NA值
  • na.trim: 过滤有NA的记录

辅助工具

  • is.regular: 检查是否是规则的序列
  • lag.zoo: 计算步长和分差
  • MATCH: 取交集
  • ORDER: 值排序,输出索引

显示控制

  • yearqtr: 以年季度显示时间
  • yearmon: 以年月显示时间
  • xblocks: 作图沿x轴分隔图型
  • make.par.list: 用于给plot.zoo 和 xyplot.zoo 数据格式转换
  • 4. zoo使用

    • 1). zoo函数
    • 2). zooreg函数
    • 3). zoo的类型转换
    • 4). ggplot2画时间序列
    • 5). 数据操作
    • 6). 数据滚动处理
    • 7). NA值处理
    • 8). 数据显示格式
    • 9). 按时间分隔做衅
    • 10). 从文件读入zoo序列

    1). zoo函数

    zoo对象包括两部分组成,数据部分、索引部分。

    函数定义:

    zoo(x = NULL, order.by = index(x), frequency = NULL)

    参数列表:

    • x: 数据部分,允许向量,矩阵,因子
    • order.by: 索引部分,唯一字段,用于排序
    • frequency: 每个时间单元显示的数量

    构建一个zoo对象,以时间为索引

    
    > x.Date <- as.Date("2003-02-01") + c(1, 3, 7, 9, 14) - 1
    > x.Date
    [1] "2003-02-01" "2003-02-03" "2003-02-07" "2003-02-09" "2003-02-14"
    > class(x.Date)
    [1] "Date"
    
    > x <- zoo(rnorm(5), x.Date)
    > x
    2003-02-01 2003-02-03 2003-02-07 2003-02-09 2003-02-14 
    0.01964254 0.03122887 0.64721059 1.47397924 1.29109889 
    > class(x)
    [1] "zoo"
    
    > plot(x)
    

    zoo1

    以数学为索引的,多组时间序列

    
    > y <- zoo(matrix(1:12, 4, 3),0:30)
    > y        
    0  1 5  9
    1  2 6 10
    2  3 7 11
    3  4 8 12
    4  1 5  9
    5  2 6 10
    6  3 7 11
    7  4 8 12
    8  1 5  9
    9  2 6 10
    10 3 7 11
    
    > plot(y)
    

    zoo2

    2). zooreg函数

    函数定义:

    zooreg(data, start = 1, end = numeric(), frequency = 1,
    deltat = 1, ts.eps = getOption("ts.eps"), order.by = NULL)

    参数列表:

    • data: 数据部分,允许向量,矩阵,因子
    • start: 时间部分,开始时间
    • end: 时间部分,结束时间
    • frequency: 每个时间单元显示的数量
    • deltat: 连续观测之间的采样周期的几分之一,不能与frequency同时出现,例如1/2
    • ts.eps: 时间序列间隔,在时间间隔大于ts.eps时认为是相等的。通过getOption(“ts.eps”)设置,默认是1e-05
    • order.by: 索引部分,唯一字段,用于排序, 继承zoo的order.by

    构建一个zooreg对象

    
    > zooreg(1:10, frequency = 4, start = c(1959, 2))
    1959(2) 1959(3) 1959(4) 1960(1) 1960(2) 1960(3) 1960(4) 1961(1) 1961(2) 
          1       2       3       4       5       6       7       8       9 
    1961(3) 
         10 
    
    > as.zoo(ts(1:10, frequency = 4, start = c(1959, 2)))
    1959(2) 1959(3) 1959(4) 1960(1) 1960(2) 1960(3) 1960(4) 1961(1) 1961(2) 
          1       2       3       4       5       6       7       8       9 
    1961(3) 
         10
    
    > zr<-zooreg(rnorm(10), frequency = 4, start = c(1959, 2))
    > plot(zr) 
    

    zoo-zr1

    3). zoo的类型转换

    转型到zoo类型

    
    > as.zoo(rnorm(5))
             1          2          3          4          5 
    -0.4892119  0.5740950  0.7128003  0.6282868  1.0289573 
    > as.zoo(ts(rnorm(5), start = 1981, freq = 12))
       1981(1)    1981(2)    1981(3)    1981(4)    1981(5) 
     2.3198504  0.5934895 -1.9375893 -1.9888237  1.0944444 
    

    从zoo类型转型到其他类型

    
    > x <- as.zoo(ts(rnorm(5), start = 1981, freq = 12))
    > x
       1981(1)    1981(2)    1981(3)    1981(4)    1981(5) 
     1.8822996  1.6436364  0.1260436 -2.0360960 -0.1387474 
    
    > as.matrix(x)
                     x
    1981(1)  1.8822996
    1981(2)  1.6436364
    1981(3)  0.1260436
    1981(4) -2.0360960
    1981(5) -0.1387474
    
    > as.vector(x)
    [1]  1.8822996  1.6436364  0.1260436 -2.0360960 -0.1387474
    
    > as.data.frame(x)
                     x
    1981(1)  1.8822996
    1981(2)  1.6436364
    1981(3)  0.1260436
    1981(4) -2.0360960
    1981(5) -0.1387474
    
    > as.list(x)
    [[1]]
       1981(1)    1981(2)    1981(3)    1981(4)    1981(5) 
     1.8822996  1.6436364  0.1260436 -2.0360960 -0.1387474 
    

    4). ggplot2画时间序列

    ggplot2::fortify函数,通过zoo::ggplot2.zoo函数,转换成ggplot2可识别的类型。

    
    library(ggplot2)
    library(scales)
    
    x.Date <- as.Date(paste(2003, 02, c(1, 3, 7, 9, 14), sep = "-"))
    x <- zoo(rnorm(5), x.Date)
    xlow <- x - runif(5)
    xhigh <- x + runif(5)
    z <- cbind(x, xlow, xhigh)
    
    g<-ggplot(aes(x = Index, y = Value), data = fortify(x, melt = TRUE))
    g<-g+geom_line()
    g<-g+geom_line(aes(x = Index, y = xlow), colour = "red", data = fortify(xlow))
    g<-g+geom_ribbon(aes(x = Index, y = x, ymin = xlow, ymax = xhigh), data = fortify(x), fill = "darkgray") 
    g<-g+geom_line()
    g<-g+xlab("Index") + ylab("x")
    g
    
    > z
                         x        xlow       xhigh
    2003-02-01 -0.36006612 -0.88751958 0.006247816
    2003-02-03  1.35216617  0.97892538 2.076360524
    2003-02-07  0.61920828  0.23746410 1.156569424
    2003-02-09  0.27516116  0.09978789 0.777878867
    2003-02-14  0.02510778 -0.80107410 0.541592929
    

    zoo-ggplot2

    5). 数据操作

    修改zoo的数据部分coredata

    
    > x.date <- as.Date(paste(2003, rep(1:4, 4:1), seq(1,20,2), sep = "-"))
    > x <- zoo(matrix(rnorm(20), ncol = 2), x.date)
    > coredata(x)
                 [,1]        [,2]
     [1,] -1.04571765  0.92606273
     [2,] -0.89621126  0.03693769
     [3,]  1.26938716 -1.06620017
     [4,]  0.59384095 -0.23845635
     [5,]  0.77563432  1.49522344
     [6,]  1.55737038  1.17215855
     [7,] -0.36540180 -1.45770721
     [8,]  0.81655645  0.09505623
     [9,] -0.06063478  0.84766496
    [10,] -0.50137832 -1.62436453
    
    > coredata(x) <- matrix(1:20, ncol = 2)
    > x        
    2003-01-01  1 11
    2003-01-03  2 12
    2003-01-05  3 13
    2003-01-07  4 14
    2003-02-09  5 15
    2003-02-11  6 16
    2003-02-13  7 17
    2003-03-15  8 18
    2003-03-17  9 19
    2003-04-19 10 20
    

    修改zoo的索引部分index

    
    > x.date <- as.Date(paste(2003, rep(1:4, 4:1), seq(1,20,2), sep = "-"))
    > x <- zoo(matrix(rnorm(20), ncol = 2), x.date)
    
    > index(x)
     [1] "2003-01-01" "2003-01-03" "2003-01-05" "2003-01-07" "2003-02-09"
     [6] "2003-02-11" "2003-02-13" "2003-03-15" "2003-03-17" "2003-04-19"
    
    > index(x) <- 1:nrow(x)
    > index(x)
     [1]  1  2  3  4  5  6  7  8  9 10
    

    按时间过滤数据window.zoo

    
    > x.date <- as.Date(paste(2003, rep(1:4, 4:1), seq(1,20,2), sep = "-"))
    > x <- zoo(matrix(rnorm(20), ncol = 2), x.date)
    
    > window(x, start = as.Date("2003-02-01"), end = as.Date("2003-03-01"))
    2003-02-09  0.7021167 -0.3073809
    2003-02-11  2.5071111  0.6210542
    2003-02-13 -1.8900271  0.1819022
    
    > window(x, index = x.date[1:6], start = as.Date("2003-02-01"))
    2003-02-09 0.7021167 -0.3073809
    2003-02-11 2.5071111  0.6210542
    
    > window(x, index = x.date[c(4, 8, 10)])
    2003-01-07  1.4623515 -1.198597
    2003-03-15 -0.5898128  1.318401
    2003-04-19 -0.4209979 -1.648222
    

    合并多个zoo对象merge.zoo

    
    > y1 <- zoo(matrix(1:10, ncol = 2), 1:5)
    > y2 <- zoo(matrix(rnorm(10), ncol = 2), 3:7)
    
    > merge(y1, y2, all = FALSE)
      y1.1 y1.2       y2.1       y2.2
    3    3    8  0.9514985  1.7238941
    4    4    9 -1.1131230 -0.2061446
    5    5   10  0.6169665 -1.3141951
    
    > merge(y1, y2, all = FALSE, suffixes = c("a", "b"))
      a.1 a.2        b.1        b.2
    3   3   8  0.9514985  1.7238941
    4   4   9 -1.1131230 -0.2061446
    5   5  10  0.6169665 -1.3141951
    
    > merge(y1, y2, all = TRUE)
      y1.1 y1.2       y2.1       y2.2
    1    1    6         NA         NA
    2    2    7         NA         NA
    3    3    8  0.9514985  1.7238941
    4    4    9 -1.1131230 -0.2061446
    5    5   10  0.6169665 -1.3141951
    6   NA   NA  0.5134937  0.0634741
    7   NA   NA  0.3694591 -0.2319775
    
    > merge(y1, y2, all = TRUE, fill = 0)
      y1.1 y1.2       y2.1       y2.2
    1    1    6  0.0000000  0.0000000
    2    2    7  0.0000000  0.0000000
    3    3    8  0.9514985  1.7238941
    4    4    9 -1.1131230 -0.2061446
    5    5   10  0.6169665 -1.3141951
    6    0    0  0.5134937  0.0634741
    7    0    0  0.3694591 -0.2319775
    

    计算zoo数据aggregate.zoo

    
    > x.date <- as.Date(paste(2004, rep(1:4, 4:1), seq(1,20,2), sep = "-"))
    > x <- zoo(rnorm(12), x.date); x
     2004-01-01  2004-01-03  2004-01-05  2004-01-07  2004-02-09  2004-02-11 
     0.67392868  1.95642526 -0.26904101 -1.24455152 -0.39570292  0.09739665 
     2004-02-13  2004-03-15  2004-03-17  2004-04-19 
    -0.23838695 -0.41182796 -1.57721805 -0.79727610 
     
    > x.date2 <- as.Date(paste(2004, rep(1:4, 4:1), 1, sep = "-")); x.date2
     [1] "2004-01-01" "2004-01-01" "2004-01-01" "2004-01-01" "2004-02-01"
     [6] "2004-02-01" "2004-02-01" "2004-03-01" "2004-03-01" "2004-04-01"
    
    > x2 <- aggregate(x, x.date2, mean); x2
    2004-01-01 2004-02-01 2004-03-01 2004-04-01 
     0.2791904 -0.1788977 -0.9945230 -0.7972761 
    

    6). 数据滚动处理

    对zoo数据的滚动处理rollapply

    
    > z <- zoo(11:15, as.Date(31:35))
    > rollapply(z, 2, mean)
    1970-02-01 1970-02-02 1970-02-03 1970-02-04 
          11.5       12.5       13.5       14.5 
    

    等价操作:rollapply , aggregate

    
    > z2 <- zoo(rnorm(6))
    > rollapply(z2, 3, mean, by = 3) # means of nonoverlapping groups of 3
             2          5 
    -0.3065197  0.6350963 
    > aggregate(z2, c(3,3,3,6,6,6), mean) # same
             3          6 
    -0.3065197  0.6350963 
    

    等价操作:rollapply, rollmean

    
    > rollapply(z2, 3, mean) # uses rollmean which is optimized for mean
             2          3          4          5 
    -0.3065197 -0.7035811 -0.1672344  0.6350963 
    > rollmean(z2, 3) # same
             2          3          4          5 
    -0.3065197 -0.7035811 -0.1672344  0.6350963 
    

    7). NA值处理
    NA填充na.fill

    
    > z <- zoo(c(NA, 2, NA, 3, 4, 5, 9, NA))
    > z
     1  2  3  4  5  6  7  8 
    NA  2 NA  3  4  5  9 NA 
    
    > na.fill(z, "extend")
      1   2   3   4   5   6   7   8 
    2.0 2.0 2.5 3.0 4.0 5.0 9.0 9.0 
    
    > na.fill(z, c("extend", NA))
     1  2  3  4  5  6  7  8 
     2  2 NA  3  4  5  9  9 
    
    > na.fill(z, -(1:3))
     1  2  3  4  5  6  7  8 
    -1  2 -2  3  4  5  9 -3 
    

    NA替换na.locf

    
    > z <- zoo(c(NA, 2, NA, 3, 4, 5, 9, NA, 11));z
     1  2  3  4  5  6  7  8  9 
    NA  2 NA  3  4  5  9 NA 11 
    
    > na.locf(z)
     2  3  4  5  6  7  8  9 
     2  2  3  4  5  9  9 11 
    
    > na.locf(z, fromLast = TRUE)
     1  2  3  4  5  6  7  8  9 
     2  2  3  3  4  5  9 11 11 
    

    统计值替换NA值na.aggregate

    
    > z <- zoo(c(1, NA, 3:9),
    +          c(as.Date("2010-01-01") + 0:2,
    +            as.Date("2010-02-01") + 0:2,
    +            as.Date("2011-01-01") + 0:2))
    > z
    2010-01-01 2010-01-02 2010-01-03 2010-02-01 2010-02-02 2010-02-03 2011-01-01 
             1         NA          3          4          5          6          7 
    2011-01-02 2011-01-03 
             8          9 
    
    > na.aggregate(z)
    2010-01-01 2010-01-02 2010-01-03 2010-02-01 2010-02-02 2010-02-03 2011-01-01 
         1.000      5.375      3.000      4.000      5.000      6.000      7.000 
    2011-01-02 2011-01-03 
         8.000      9.000 
    
    > na.aggregate(z, as.yearmon)
    2010-01-01 2010-01-02 2010-01-03 2010-02-01 2010-02-02 2010-02-03 2011-01-01 
             1          2          3          4          5          6          7 
    2011-01-02 2011-01-03 
             8          9 
    
    > na.aggregate(z, months)
    2010-01-01 2010-01-02 2010-01-03 2010-02-01 2010-02-02 2010-02-03 2011-01-01 
           1.0        5.6        3.0        4.0        5.0        6.0        7.0 
    2011-01-02 2011-01-03 
           8.0        9.0 
    
    > na.aggregate(z, format, "%Y")
    2010-01-01 2010-01-02 2010-01-03 2010-02-01 2010-02-02 2010-02-03 2011-01-01 
           1.0        3.8        3.0        4.0        5.0        6.0        7.0 
    2011-01-02 2011-01-03 
           8.0        9.0 
    

    计算插值替换NA值

    
    > z <- zoo(c(2, NA, 1, 4, 5, 2), c(1, 3, 4, 6, 7, 8));z
     1  3  4  6  7  8 
     2 NA  1  4  5  2 
    
    > na.approx(z)
           1        3        4        6        7        8 
    2.000000 1.333333 1.000000 4.000000 5.000000 2.000000 
    
    > na.approx(z, 1:6)
      1   3   4   6   7   8 
    2.0 1.5 1.0 4.0 5.0 2.0 
    

    计算seasonal Kalman filter替换NA值

    
    z <- zooreg(rep(10 * seq(4), each = 4) + rep(c(3, 1, 2, 4), times = 4),
                start = as.yearqtr(2000), freq = 4)
    z[10] <- NA
    zout <- na.StructTS(z);zout
    plot(cbind(z, zout), screen = 1, col = 1:2, type = c("l", "p"), pch = 20)
    

    zoo-zout

    过滤有NA的行

    
    > xx <- zoo(matrix(c(1, 4, 6, NA, NA, 7), 3), c(2, 4, 6));xx
    2 1 NA
    4 4 NA
    6 6  7
    
    > na.trim(xx)
    6 6 7
    

    8). 数据显示格式
    以年+季度格式输出

    
    > x <- as.yearqtr(2000 + seq(0, 7)/4)
    > x
    [1] "2000 Q1" "2000 Q2" "2000 Q3" "2000 Q4" "2001 Q1" "2001 Q2" "2001 Q3"
    [8] "2001 Q4"
    
    > format(x, "%Y Quarter %q")
    [1] "2000 Quarter 1" "2000 Quarter 2" "2000 Quarter 3" "2000 Quarter 4"
    [5] "2001 Quarter 1" "2001 Quarter 2" "2001 Quarter 3" "2001 Quarter 4"
    
    > as.yearqtr("2001 Q2")
    [1] "2001 Q2"
    
    > as.yearqtr("2001 q2") 
    [1] "2001 Q2"
    
    > as.yearqtr("2001-2")
    [1] "2001 Q2"
    

    以年+月份格式输出

    
    > x <- as.yearmon(2000 + seq(0, 23)/12)
    > x
     [1] "一月 2000"   "二月 2000"   "三月 2000"   "四月 2000"   "五月 2000"  
     [6] "六月 2000"   "七月 2000"   "八月 2000"   "九月 2000"   "十月 2000"  
    [11] "十一月 2000" "十二月 2000" "一月 2001"   "二月 2001"   "三月 2001"  
    [16] "四月 2001"   "五月 2001"   "六月 2001"   "七月 2001"   "八月 2001"  
    [21] "九月 2001"   "十月 2001"   "十一月 2001" "十二月 2001"
    
    > as.yearmon("mar07", "%b%y")
    [1] NA
    
    > as.yearmon("2007-03-01")
    [1] "三月 2007"
    
    > as.yearmon("2007-12")
    [1] "十二月 2007"
    

    9). 按时间分隔线

    使用xblock函数, 以不同的颜色划分3个区间(-Inf,15),[15,30],(30,Inf)

    
    set.seed(0)
    flow <- ts(filter(rlnorm(200, mean = 1), 0.8, method = "r"))
    rgb <- hcl(c(0, 0, 260), c = c(100, 0, 100), l = c(50, 90, 50), alpha = 0.3)
    plot(flow)
    xblocks(flow > 30, col = rgb[1]) ## high values red
    xblocks(flow < 15, col = rgb[3]) ## low value blue
    xblocks(flow >= 15 & flow <= 30, col = rgb[2]) ## the rest gray
    

    zoo-xblocks

    10). 从文件读入zoo序列
    创建文件:read.csv

    
    ~ vi read.csv
    
    2003-01-01,1.0073644,0.05579711
    2003-01-03,-0.2731580,0.06797239
    2003-01-05,-1.3096795,-0.20196174
    2003-01-07,0.2225738,-1.15801525
    2003-02-09,1.1134332,-0.59274327
    2003-02-11,0.8373944,0.76606538
    2003-02-13,0.3145168,0.03892812
    2003-03-15,0.2222181,0.01464681
    2003-03-17,-0.8436154,-0.18631697
    2003-04-19,0.4438053,1.40059083
    

    读文件并生成zoo序列

    
    > r <- read.zoo(file="read.csv",sep = ",", format = "%Y-%m-%d")
    
    > r
                       V2          V3
    2003-01-01  1.0073644  0.05579711
    2003-01-03 -0.2731580  0.06797239
    2003-01-05 -1.3096795 -0.20196174
    2003-01-07  0.2225738 -1.15801525
    2003-02-09  1.1134332 -0.59274327
    2003-02-11  0.8373944  0.76606538
    2003-02-13  0.3145168  0.03892812
    2003-03-15  0.2222181  0.01464681
    2003-03-17 -0.8436154 -0.18631697
    2003-04-19  0.4438053  1.40059083
    
    > class(r)
    [1] "zoo"
    

    我们已经完全掌握了zoo库及zoo对象的使用,接下来就可以放手去用R处理时间序列了!

    转载请注明出处:
    http://blog.fens.me/r-zoo/

    打赏作者