title	author	date	output
2021项目评估时间序列预测	张剑	2021/4/23	html_document

```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` ## 1. Plot the Time series graphics > Graphs enable many features of the data to be visualised, including patterns, unusual observations, changes over time, and relationships between variables. ### tsibble objects > A time series can be thought of as a list of numbers, along with some information about what times those numbers were recorded. 假设我们有近些年的年度数据： | year | Obs | |:----:|:---:| | 2012 | 123 | | 2013 | 39 | | 2014 | 78 | | 2015 | 52 | | 2016 | 110 | We turn this into a tsibble object using the tsibble() function: ```{r message=FALSE, warning=FALSE} library(fpp3) y % mutate(Month_new=yearmonth(Month)) %>% as_tsibble(index = Month_new) %>% select(Month_new,obs) z_ts ``` > 首先，将“月”列从文本转换为带有yearmonth()的月时间对象。然后，我们通过使用as _ tsibble()标识索引变量，将数据帧转换为tsible。注意第一行加“[1M]”，表示这是月度数据，每月观察一次。 ```{r} library(tsibbledata) olympic_running ``` > olympic_running 数据集是一个 tsibble object，包含了312行，4列。“[4Y]”告知我们这些观测的间隔是每四年一次。 可以使用distinct()函数显示类别变量的特征 ```{r} olympic_running %>% distinct(Sex) ``` ```{r} help(PBS) PBS %>% filter(ATC2 == "A10") %>% select(Month, Concession, Type, Cost) ``` ```{r} PBS %>% filter(ATC2 == "A10") %>% select(Month, Concession, Type, Cost) %>% summarise(TotalC = sum(Cost)) %>% mutate(Cost = TotalC / 1e6) -> a10 ``` 读取csv,该数据集按州，性别，法律地位和土著地位提供了有关澳大利亚监狱人口规模的信息。 ```{r message=FALSE, warning=FALSE} library(readr) prison % mutate(Quarter = yearquarter(Date)) %>% # 使用yearquarter函数将原始的date转换为时间序列格式 select(-Date) %>% as_tsibble(key = c(State, Gender, Legal, Indigenous), index = Quarter) prison ``` Indigenous [64] 代表这个tisbble包含了64个独立时间。 ## 时间趋势图 ansett dataset： Passenger numbers on Ansett airline flights autoplot函数可以非常方便的绘制时间序列图，放入参数数据，需要关注的变量。 ```{r} help(ansett) melsyd_economy % filter(Airports == "MEL-SYD", Class == "Economy") %>% mutate(Passengers = Passengers/1000) melsyd_economy autoplot(melsyd_economy, Passengers) + labs(title = "Ansett airlines economy class", subtitle = "Melbourne-Sydney", y = "Passengers ('000)") ``` #### 时间图立即显示出一些有趣的特征： - 1989年有一段时期，没有乘客载运，这是由于劳资纠纷造成的。 - 在1992年有一段时间，载客量有所减少 - 1991年下半年，客运量大大增加。 - 每年年初，负载都有一些大的下降。 这些是由于假期的影响。 - 该系列的水平存在长期波动，该波动在1987年增加，在1989年减少，并在1990年和1991年再次增加。 #### 上面做的a10数据的时间图 ```{r} autoplot(a10, Cost) + labs(y = "$ (millions)", title = "Australian antidiabetic drug sales") ``` - 在这里，有一个明显且不断增长的趋势。 - 还有一个强烈的季节性模式，随着时间的推移，其大小也会增加。 - 每年年初会出现一个大幅度的上升。 ### 时间序列的特征 #### 接下来对趋势trend，季节seasonal，周期Cyclic进行界定。 #### 趋势Trend > A trend exists when there is a long-term increase or decrease in the data. It does not have to be linear. Sometimes we will refer to a trend as “changing direction,” when it might go from an increasing trend to a decreasing trend. #### 季节seasonal > A seasonal pattern occurs when a time series is affected by seasonal factors such as the time of the year or the day of the week. Seasonality is always of a fixed and known period. #### 周期Cyclic > A cycle occurs when the data exhibit rises and falls that are not of a fixed frequency. These fluctuations are usually due to economic conditions, and are often related to the “business cycle.” The duration of these fluctuations is usually at least 2 years. > 许多人将周期性行为与季节性行为混为一谈，但实际上却大不相同。 如果波动的频率不是固定的，则它们是周期性的。 如果频率是不变的并且与日历的某些方面相关联，则该模式是季节性的。 通常，周期的平均长度比季节性样式的长度长，周期的大小往往比季节样式的大小可变。 ### 季节图 使用feasts包下面的gg_season函数可以非常轻松的画出优美的季节图 ```{r} a10 %>% gg_season(Cost, labels = "both") + labs(y = "$ (millions)", title = "Seasonal plot: Antidiabetic drug sales") + expand_limits(x = ymd(c("1972-12-28", "1973-12-04"))) ``` 这些数据与之前显示的完全相同，但是现在每个季节的数据已经重叠。 季节性图可以更清楚地看到基本的季节性模式，并且在确定模式发生变化的年份时特别有用。 在这种情况下，很明显，每年一月份的销售额都有很大的增长。 实际上，这些可能是12月下旬的销售，因为客户在日历年末之前有库存，但是直到一两周后才向政府注册销售。 该图还显示，2008年3月的销售数量异常少（大多数其他年份显示2月至3月之间有所增长）。 2008年6月的销售量很少，可能是由于在收集数据时对销售的计数不完整。 ### 不同时间段的季节图 使用gg_season函数是，如果数据具有多个季节模式，则可以使用period参数选择所需的季节图。比如我们可以尝试日，周，月，季等。 我们使用vic_elec数据库来展示多个季节模式，澳大利亚维多利亚州每半小时的用电需求. - Demand - Temperature - Holiday 数据频率30m一次，第一张图用了period day的参数，展示的是日频率数据，实际上看不清楚。 可以尝试使用周或者月 ```{r} vic_elec vic_elec %>% gg_season(Demand, period = "day") + theme(legend.position = "none") + labs(y="MW", title="Electricity demand: Victoria") ``` ```{r} vic_elec %>% gg_season(Demand, period = "week") + theme(legend.position = "none") + labs(y="MW", title="Electricity demand: Victoria") ``` ```{r} vic_elec %>% gg_season(Demand, period = "year") + labs(y="MW", title="Electricity demand: Victoria") ``` ### 季节性子系列图 强调季节模式的另一种方法是将每个季节的数据一起收集在单独的迷你时间图中。这样有可能会更清楚。 ```{r} a10 %>% gg_subseries(Cost) + labs( y = "$ (millions)", title = "Australian antidiabetic drug sales" ) ``` 蓝色水平线表示每个月的平均值。这种形式的图可以清楚地看到潜在的季节模式，也显示了季节性随时间的变化。它在识别特定季节内的变化时特别有用。这个例子中，并不是特别具有启发性；但是在某些情况下，这是观察季节变化最有用的方法。 #### 使用tourism数据库展示季节性子图的威力 tourism {tsibble} 数据集包含从1998年Q1到2016年Q5澳大利亚的季度夜间旅行。Pourpose，筛选出holiday的，分类汇总一下每个state，每个时间频率上的的trips人数 ```{r} tourism %>% filter(Purpose == 'Holiday') %>% group_by(State) %>% summarise(Trips = sum(Trips)) -> holiday holiday holiday %>% autoplot(Trips) + labs(y = "Overnight trips ('000)", title = "Australian domestic holidays") ``` ```{r} holiday %>% gg_subseries(Trips) + labs(y = "Overnight trips ('000)", title = "Australian domestic holidays") ``` 这一数字表明，西澳大利亚旅游业近年来已显著跃升，而维多利亚州旅游业在Q1和第四季度有所增长，但不是在年中。 ### Scatterplots散点图 #### 2014年用电情况 ```{r} vic_elec %>% filter(year(Time) == 2014) %>% autoplot(Demand) + labs(y = "GW", title = "Half-hourly electricity demand: Victoria") ``` #### 2014年温度情况 ```{r} vic_elec %>% filter(year(Time) == 2014) %>% autoplot(Temperature) + labs( y = "Degrees Celsius", title = "Half-hourly temperatures: Melbourne, Australia" ) ``` #### 将气温和用电量放在要给坐标系下 ```{r} vic_elec %>% filter(year(Time) == 2014) %>% ggplot(aes(x = Temperature, y = Demand)) + geom_point() + labs(x = "Temperature (degrees Celsius)", y = "Electricity demand (GW)") ``` > 这个散点图有助于我们把变量之间的关系形象化。很明显，当温度由于空调的影响而很高时，就会出现高需求。但也有一个加热效应，在那里对极低温度的需求增加。 ### Correlation 相关系数 $$ r=\frac{\sum\left(x_{t}-\bar{x}\right)\left(y_{t}-\bar{y}\right)}{\sqrt{\sum\left(x_{t}-\bar{x}\right)^{2}} \sqrt{\sum\left(y_{t}-\bar{y}\right)^{2}}} $$ ### 矩阵散点图 ```{r message=FALSE, warning=FALSE} holiday %>% pivot_wider(values_from = Trips,names_from = State) %>% GGally::ggpairs(columns = 2:9) ``` #### 散点图矩阵的价值在于其使得能够快速查看所有变量对之间的关系。 ### Lag plots 滞后图 使用 aus_production aus_production {tsibbledata} Quarterly estimates of selected indicators of manufacturing production in Australia. Beer: Beer production in megalitres. 使用2000年后的数据 使用gg_lag函数，非常容易的绘制出之后图 ```{r} recent_production % filter(year(Quarter) >= 2000) recent_production %>% gg_lag(Beer, geom = "point") + labs(x = "lag(Beer, k)") ``` 这里的颜色表示垂直轴上变量的季度。这种关系在滞后4和8处呈强正相关，反映了数据的强烈季节性。 ### Autocorrelation 自相关 $$ r_{k}=\frac{\sum_{t=k+1}^{T}\left(y_{t}-\bar{y}\right)\left(y_{t-k}-\bar{y}\right)}{\sum_{t=1}^{T}\left(y_{t}-\bar{y}\right)^{2}} $$ 其中 T 是时间序列的长度 ```{r} recent_production %>% ACF(Beer, lag_max = 9) recent_production %>% ACF(Beer) %>% autoplot() + labs(title="Australian beer production") ``` #### 通过上面的数据和图，可以发现： - r4高于其他滞后。这是由于数据中的季节性模式:波峰往往相隔四个季度，波谷往往相隔四个季度 - r2比其他滞后指标更负，因为低谷期往往比峰值期晚四分之二。 - 蓝色虚线表示相关性是否与零显著不同 #### 当数据具有趋势时，小滞后的自相关往往较大且为正，因为时间上邻近的观察值也在价值上邻近。因此，趋势时间序列的ACF往往具有正值，该正值随滞后增加而缓慢减小。 #### 当数据是季节性的时，季节性滞后(以季节性周期的倍数)的自相关将大于其他滞后。 #### 当数据同时具有趋势和季节性时，就会看到这些影响的组合。 ```{r} a10 %>% ACF(Cost, lag_max = 48) %>% autoplot() + labs(title="Australian antidiabetic drug sales") ``` ### White noise 白噪声 没有自相关的时间序列被称为白噪音 ```{r} set.seed(30) y % autoplot(wn) + labs(title = "White noise", y = "") ``` ```{r} y %>% ACF(wn) %>% autoplot() + labs(title = "White noise") ``` > 对于白噪声序列，我们期望每个自相关接近于零。当然，它们不会完全等于零，因为存在一些随机变化。 ### 练习 #### 1. 调入pp3包，用help函数观察gafa_stock, PBS, vic_elec和 pelt数据库。 - 使用 autoplot()函数去这些数据集中的某些感兴趣的变量的时间图 - 这些你感兴趣的变量的时间间隔是多少？ #### 2. 使用filter()查找gafa_stock中四只股票中每只股票的峰值收盘价对应的日期。 #### 3. 下载 tute1.csv 文件，[文件下载地址](https://gitee.com/jefeerzhang/forcast/raw/master/data/tute1.csv)，在excel中打开，观察tute1数据的结构。使用readr包中的read_csv读入数据，并观察数据结构。 - 将数据改为时间序列结构 tsibble结构。 ```R mytimeseries % mutate(Quarter = yearmonth(Quarter)) %>% as_tsibble(index = Quarter) ``` - 构建三个系列中每个系列的时间序列图 ```R mytimeseries %>% pivot_longer(-Quarter) %>% ggplot(aes(x = Quarter, y = value, colour = name)) + geom_line() + facet_grid(name ~ ., scales = "free_y") ``` 测试一下，如果不加最后一行 facet_grid(name ~ ., scales = "free_y")，会发生什么？ #### 4. - 从[这里下载](https://gitee.com/jefeerzhang/forcast/raw/master/data/tourism.xlsx)tourism.xlsx数据。 - 读入数据，并将其转换为tsibble格式的数据。 #### 5. - 创建以下四个时间序列的时间图:来自aus_production的Bricks、来自pelt的Lynx、来自gafa_stock的Close、来自vic _ elec的Demand。 - 使用help函数了解每个数据库的结构和各类指标 - 对于最后一个图，修改轴标签和标题。 #### 6. 使用以下图形函数: autoplot(), gg_season(), gg_subseries(), gg_lag(), ACF() 并浏览下列时间序列中的特征: “Total Private”来自us_employment，Bricks来自aus_production，Hare来自pelt,以及us_gasoline。 - 你能看出任何季节性、周期性和趋势吗？ - 每个序列，你发现了什么？ - 对于每个系列的季节性趋势有什么看法？ - 是否发现了一些不同寻常的年份？