【免责声明：本文用于教学】

时间序列分析

基础操作

数据输入

d <- c(10,15,10,10,12,10,7,7,10,14,8,17,14,18,3,9,11,10,6,12,14,10,25,29,33,33,12,19,16,19,19,12,34,15,36,29,26,21,17,19,13,20,24,12,6,14,6,12,9,11,17,12,8,14,14,12,5,8,10,3,16,8,8,7,12,6,10,8,10,5)
##convert it to a time series
is.ts(d)

## [1] FALSE

d <- ts(d, start=min(201401))
d

## Time Series:
## Start = 201401
## End = 201470
## Frequency = 1
## [1] 10 15 10 10 12 10 7 7 10 14 8 17 14 18 3 9 11 10 6 12 14 10 25 29 33
## [26] 33 12 19 16 19 19 12 34 15 36 29 26 21 17 19 13 20 24 12 6 14 6 12 9 11
## [51] 17 12 8 14 14 12 5 8 10 3 16 8 8 7 12 6 10 8 10 5

plot.ts(d,ylab="y", xlab="day")

ACF和PACF

par(mfrow=c(1,2))
##autocovariance function ACF 自相关
acf(d)
acf(d,plot=F)

##
## Autocorrelations of series 'd', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 0.506 0.539 0.374 0.291 0.258 0.148 0.270 0.186 0.178 0.258
## 11 12 13 14 15 16 17 18
## 0.207 0.226 0.138 -0.027 -0.053 -0.112 -0.139 -0.155

##partial autocovariance function PACF 偏自相关
pacf(d)

pacf(d,plot=F)

##
## Partial autocorrelations of series 'd', by lag
##
## 1 2 3 4 5 6 7 8 9 10 11
## 0.506 0.380 0.020 -0.055 0.046 -0.062 0.207 0.032 -0.069 0.163 0.037
## 12 13 14 15 16 17 18
## -0.023 -0.055 -0.318 -0.101 0.087 -0.091 -0.060

检验

单位根检验 : 非平稳 vs : 平稳

library(fUnitRoots)

## Loading required package: timeDate

## Loading required package: timeSeries

## Loading required package: fBasics

unitrootTest(d)

##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 1
## STATISTIC:
## DF: -1.2139
## P VALUE:
## t: 0.2041
## n: 0.4322
##
## Description:
## Mon Dec 7 00:41:14 2020 by user:

纯随机性 :白噪声 vs : 非白噪声

Box.test(d,type="Ljung-Box")

##
## Box-Ljung test
##
## data: d
## X-squared = 18.7, df = 1, p-value = 1.53e-05

Box.test(d,type="Ljung-Box",lag=6)

##
## Box-Ljung test
##
## data: d
## X-squared = 64.016, df = 6, p-value = 6.85e-12

Box.test(d,type="Ljung-Box",lag=12)

##
## Box-Ljung test
##
## data: d
## X-squared = 88.975, df = 12, p-value = 7.794e-14

做差分

d1 <- diff(d)
plot.ts(d1,ylab="diff", xlab="day")

par(mfrow=c(1,2))
##autocovariance function ACF 自相关
acf(d1)
acf(d1,plot=F)

##
## Autocorrelations of series 'd1', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 -0.529 0.195 -0.080 -0.059 0.092 -0.256 0.216 -0.075 -0.070 0.101
## 11 12 13 14 15 16 17 18
## -0.048 0.104 0.075 -0.142 0.045 -0.032 -0.026 -0.022

##partial autocovariance function PACF 偏自相关
pacf(d1)

pacf(d1,plot=F)

##
## Partial autocorrelations of series 'd1', by lag
##
## 1 2 3 4 5 6 7 8 9 10 11
## -0.529 -0.119 -0.039 -0.142 -0.011 -0.283 -0.095 0.015 -0.175 -0.091 -0.015
## 12 13 14 15 16 17 18
## 0.024 0.272 0.038 -0.121 0.043 0.002 -0.057

##单位根检验 H0: 非平稳
library(fUnitRoots)
unitrootTest(d1)

##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 1
## STATISTIC:
## DF: -7.9309
## P VALUE:
## t: < 2.2e-16
## n: 0.0461
##
## Description:
## Mon Dec 7 00:41:14 2020 by user:

##纯随机性 H0:白噪声
Box.test(d1,type="Ljung-Box")

##
## Box-Ljung test
##
## data: d1
## X-squared = 20.195, df = 1, p-value = 6.993e-06

Box.test(d1,type="Ljung-Box",lag=6)

##
## Box-Ljung test
##
## data: d1
## X-squared = 29.458, df = 6, p-value = 4.982e-05

Box.test(d1,type="Ljung-Box",lag=12)

##
## Box-Ljung test
##
## data: d1
## X-squared = 35.943, df = 12, p-value = 0.0003308

ts函数

R中自带的ts()函数可以直接生成时序数据，调用公式如下：

ts(data = NA, start = , end = , frequency = , ...)
frequency：取值包括1,4,12,分别代表年度,季度,月度数据

##生成季度数据,且以1959年第二季度为起始点(月度数据同理)
ts(1:10, frequency = 4, start = c(1959, 2))

## Qtr1 Qtr2 Qtr3 Qtr4
## 1959 1 2 3
## 1960 4 5 6 7
## 1961 8 9 10

# 同时生成多条时间序列
data <- ts(matrix(rnorm(180), 60, 3), start = c(2012, 1), frequency = 12,names = c('a','b','c'))
head(data)

## a b c
## [1,] 1.38880495 -1.7780026 -0.8546904
## [2,] -0.03141699 0.9623493 1.2816379
## [3,] 0.69086949 0.5687210 -2.0918558
## [4,] 0.86851666 -0.9937985 2.2763149
## [5,] -0.58213129 -0.2361856 0.2417787
## [6,] 0.69791200 -1.3831174 0.6326396

plot(data)

ts.data <- ts(data,start=c(1996,10),freq=12)
plot(ts.data)

# 使用window函数选取部分数据
data_subset <- window(ts.data,start=c(1998,1),end=c(1998,6))
data_subset

## a b c
## Jan 1998 0.09813527 0.2172911 0.10958001
## Feb 1998 1.78146726 -0.7813561 0.39675998
## Mar 1998 -1.07601484 0.2697066 -0.87731481
## Apr 1998 0.53832758 1.5250992 -0.09629001
## May 1998 0.71610457 0.8583480 -0.29970496
## Jun 1998 -0.33019984 0.2186358 -1.36307364

时间序列分解

简单移动平均

R中的forecast包的ma()函数、zoo包的rollmean()函数和TTR包的SMA()函数均可实现简单移动平均处理，这里我们选用TTR包的SMA函数

# nhtemp数据集为康涅狄格州纽黑文地区从1912年至1971年每年的平均气温时序数据
data(nhtemp)
ylim <- c(45,55)
library(TTR)

##
## Attaching package: 'TTR'

## The following object is masked from 'package:fBasics':
##
## volatility

par(mfrow=c(2,2))
plot(nhtemp, main="Raw time series")
plot(SMA(nhtemp,3), main="SMA(k=3)",ylim=ylim)
plot(SMA(nhtemp,5), main="SMA(k=5)",ylim=ylim)
plot(SMA(nhtemp,7), main="SMA(k=7)",ylim=ylim)

stl函数

stl(ts, s.window=, t.window=)

s.window：控制季节效应变化的速度,设定为"periodic"时可使得季节效应在各年间都一样
t.window：控制趋势项变化的速度,较小的值意味着更快的变化速度

data(nottem)
head(nottem,24)

## [1] 40.6 40.8 44.4 46.7 54.1 58.5 57.7 56.4 54.3 50.5 42.9 39.8 44.2 39.8 45.1
## [16] 47.0 54.1 58.7 66.3 59.9 57.0 54.2 39.7 42.8

stl(nottem, s.window = "periodic", t.window = 50)

## Call:
## stl(x = nottem, s.window = "periodic", t.window = 50)
##
## Components
## seasonal trend remainder
## Jan 1920 -9.3542441 49.94158 0.01266149
## Feb 1920 -9.8589768 49.89992 0.75905251
## Mar 1920 -6.8533560 49.85827 1.39509007
## Apr 1920 -2.7528761 49.81661 -0.36373152
## May 1920 3.5211867 49.77495 0.80386403
## Jun 1920 8.9998830 49.73329 -0.23317393
## Jul 1920 12.8601452 49.68456 -4.84470197
## Aug 1920 11.4833902 49.63582 -4.71921287
## Sep 1920 7.4481630 49.58709 -2.73525153
## Oct 1920 0.4621026 49.53835 0.49954308
## Nov 1920 -6.4543306 49.48962 -0.13528963
## Dec 1920 -9.5010869 49.43627 -0.13518287
## Jan 1921 -9.3542441 49.38292 4.17132468
## Feb 1921 -9.8589768 49.32957 0.32940779
## Mar 1921 -6.8533560 49.27622 2.67713744
## Apr 1921 -2.7528761 49.22287 0.53000795
## May 1921 3.5211867 49.17359 1.40522344
## Jun 1921 8.9998830 49.12431 0.57580543
## Jul 1921 12.8601452 49.07503 4.36482150
## Aug 1921 11.4833902 49.02576 -0.60914528
## Sep 1921 7.4481630 48.97648 0.57536017
## Oct 1921 0.4621026 48.94145 4.79644894
## Nov 1921 -6.4543306 48.90642 -2.75208962
## Dec 1921 -9.5010869 48.87139 3.42969504
## Jan 1922 -9.3542441 48.83636 -1.98211951
## Feb 1922 -9.8589768 48.80134 -0.24235850
## Mar 1922 -6.8533560 48.72863 -2.37527800
## Apr 1922 -2.7528761 48.65593 -3.80305665
## May 1922 3.5211867 48.58323 3.59558183
## Jun 1922 8.9998830 48.51053 0.28958681
## Jul 1922 12.8601452 48.43783 -4.49797412
## Aug 1922 11.4833902 48.34906 -5.53244838
## Sep 1922 7.4481630 48.26029 -1.40845040
## Oct 1922 0.4621026 48.17152 -1.53361914
## Nov 1922 -6.4543306 48.08275 0.17158479
## Dec 1922 -9.5010869 47.99397 3.20711193
## Jan 1923 -9.3542441 47.96499 3.18925666
## Feb 1923 -9.8589768 47.93600 2.02297694
## Mar 1923 -6.8533560 47.90701 1.84634377
## Apr 1923 -2.7528761 47.87802 0.67485145
## May 1923 3.5211867 47.84904 -2.17022374
## Jun 1923 8.9998830 47.88889 -4.18877107
## Jul 1923 12.8601452 47.92874 3.41111570
## Aug 1923 11.4833902 47.96859 0.14801960
## Sep 1923 7.4481630 48.00844 -1.05660426
## Oct 1923 0.4621026 48.04829 0.68960516
## Nov 1923 -6.4543306 48.10019 -5.34585958
## Dec 1923 -9.5010869 48.15209 -1.05100111
## Jan 1924 -9.3542441 48.20399 0.45025814
## Feb 1924 -9.8589768 48.25588 -0.89690704
## Mar 1924 -6.8533560 48.30778 -3.15442568
## Apr 1924 -2.7528761 48.33679 -0.08391365
## May 1924 3.5211867 48.36580 1.31301551
## Jun 1924 8.9998830 48.39481 0.30531116
## Jul 1924 12.8601452 48.42381 -0.48395910
## Aug 1924 11.4833902 48.45282 -1.73621222
## Sep 1924 7.4481630 48.49354 0.45829771
## Oct 1924 0.4621026 48.53426 0.80364092
## Nov 1924 -6.4543306 48.57497 2.27935680
## Dec 1924 -9.5010869 48.61569 4.48539590
## Jan 1925 -9.3542441 48.65641 0.69783579
## Feb 1925 -9.8589768 48.70660 1.65237963
## Mar 1925 -6.8533560 48.75679 -1.10342997
## Apr 1925 -2.7528761 48.80697 -0.95409873
## May 1925 3.5211867 48.85716 1.42164964
## Jun 1925 8.9998830 48.90735 1.49276451
## Jul 1925 12.8601452 48.93118 1.70866992
## Aug 1925 11.4833902 48.95502 0.56159246
## Sep 1925 7.4481630 48.97885 -3.42701276
## Oct 1925 0.4621026 49.00268 0.53521530
## Nov 1925 -6.4543306 49.02651 -4.47218397
## Dec 1925 -9.5010869 49.00862 -3.20752918
## Jan 1926 -9.3542441 48.99072 -0.43647360
## Feb 1926 -9.8589768 48.97282 4.28615754
## Mar 1926 -6.8533560 48.95492 1.29843523
## Apr 1926 -2.7528761 48.93702 2.71585376
## May 1926 3.5211867 48.91284 -1.83402983
## Jun 1926 8.9998830 48.88866 -1.08854692
## Jul 1926 12.8601452 48.86448 0.77537007
## Aug 1926 11.4833902 48.84031 1.67630420
## Sep 1926 7.4481630 48.81613 1.23571057
## Oct 1926 0.4621026 48.80908 -2.57117867
## Nov 1926 -6.4543306 48.80203 -0.74769523
## Dec 1926 -9.5010869 48.79498 0.50611142
## Jan 1927 -9.3542441 48.78793 -0.03368114
## Feb 1927 -9.8589768 48.78087 -0.42189814
## Mar 1927 -6.8533560 48.75952 3.39383915
## Apr 1927 -2.7528761 48.73816 1.11471730
## May 1927 3.5211867 48.71680 -0.53798743
## Jun 1927 8.9998830 48.69544 -2.69532566
## Jul 1927 12.8601452 48.67408 -1.13422980
## Aug 1927 11.4833902 48.62411 0.39249606
## Sep 1927 7.4481630 48.57414 -1.32230584
## Oct 1927 0.4621026 48.52417 1.31372554
## Nov 1927 -6.4543306 48.47420 0.28012960
## Dec 1927 -9.5010869 48.42423 -3.72314313
## Jan 1928 -9.3542441 48.40188 1.75236119
## Feb 1928 -9.8589768 48.37954 2.57944107
## Mar 1928 -6.8533560 48.35719 1.29616749
## Apr 1928 -2.7528761 48.33484 1.71803477
## May 1928 3.5211867 48.31249 -0.93368082
## Jun 1928 8.9998830 48.32363 -0.92351403
## Jul 1928 12.8601452 48.33477 1.00508685
## Aug 1928 11.4833902 48.34590 0.67070487
## Sep 1928 7.4481630 48.35704 -0.40520487
## Oct 1928 0.4621026 48.36818 1.36971866
## Nov 1928 -6.4543306 48.39305 1.06128043
## Dec 1928 -9.5010869 48.41792 -1.61683458
## Jan 1929 -9.3542441 48.44279 -4.28854880
## Feb 1929 -9.8589768 48.46766 -7.30868747
## Mar 1929 -6.8533560 48.49254 -0.63917959
## Apr 1929 -2.7528761 48.50596 -1.85308478
## May 1929 3.5211867 48.51939 1.05942716
## Jun 1929 8.9998830 48.53281 -0.63269441
## Jul 1929 12.8601452 48.54624 1.09361812
## Aug 1929 11.4833902 48.55966 0.25694778
## Sep 1929 7.4481630 48.56926 3.78257706
## Oct 1929 0.4621026 48.57886 0.15903962
## Nov 1929 -6.4543306 48.58846 0.76587484
## Dec 1929 -9.5010869 48.59805 2.80303329
## Jan 1930 -9.3542441 48.60765 2.34659252
## Feb 1930 -9.8589768 48.62978 -1.67080301
## Mar 1930 -6.8533560 48.65191 -0.59855199
## Apr 1930 -2.7528761 48.67404 0.97883988
## May 1930 3.5211867 48.69616 -1.01735113
## Jun 1930 8.9998830 48.71829 2.68182437
## Jul 1930 12.8601452 48.73106 -1.49120982
## Aug 1930 11.4833902 48.74384 1.37277313
## Sep 1930 7.4481630 48.75661 0.79522832
## Oct 1930 0.4621026 48.76938 1.66851679
## Nov 1930 -6.4543306 48.78215 0.67217793
## Dec 1930 -9.5010869 48.76693 -0.46583922
## Jan 1931 -9.3542441 48.75170 -2.29745557
## Feb 1931 -9.8589768 48.73647 -0.47749637
## Mar 1931 -6.8533560 48.72125 -3.46789062
## Apr 1931 -2.7528761 48.70602 0.54685598
## May 1931 3.5211867 48.70229 1.27652278
## Jun 1931 8.9998830 48.69856 0.70155608
## Jul 1931 12.8601452 48.69483 -0.95497654
## Aug 1931 11.4833902 48.69110 -1.97449201
## Sep 1931 7.4481630 48.68737 -2.33553525
## Oct 1931 0.4621026 48.73109 -2.59319717
## Nov 1931 -6.4543306 48.77482 3.17951358
## Dec 1931 -9.5010869 48.81854 1.28254755
## Jan 1932 -9.3542441 48.86226 2.89198231
## Feb 1932 -9.8589768 48.90598 -0.64700737
## Mar 1932 -6.8533560 48.96480 -1.81144612
## Apr 1932 -2.7528761 49.02362 -1.67074402
## May 1932 3.5211867 49.08244 -1.70362479
## Jun 1932 8.9998830 49.14126 -1.14113906
## Jul 1932 12.8601452 49.20007 0.03978075
## Aug 1932 11.4833902 49.25096 2.76565272
## Sep 1932 7.4481630 49.30184 -0.45000308
## Oct 1932 0.4621026 49.35272 -2.51482559
## Nov 1932 -6.4543306 49.40361 0.65072456
## Dec 1932 -9.5010869 49.45449 1.84659794
## Jan 1933 -9.3542441 49.50632 -3.95207345
## Feb 1933 -9.8589768 49.55815 -0.39916928
## Mar 1933 -6.8533560 49.60997 1.74338143
## Apr 1933 -2.7528761 49.66180 1.79107300
## May 1933 3.5211867 49.71363 0.96518169
## Jun 1933 8.9998830 49.77232 2.02779414
## Jul 1933 12.8601452 49.83101 2.80884066
## Aug 1933 11.4833902 49.88971 3.52690433
## Sep 1933 7.4481630 49.94840 2.70344024
## Oct 1933 0.4621026 50.00709 -0.26919058
## Nov 1933 -6.4543306 50.04874 -1.49440788
## Dec 1933 -9.5010869 50.09039 -4.78930196
## Jan 1934 -9.3542441 50.13204 -1.37779526
## Feb 1934 -9.8589768 50.17369 -2.11471300
## Mar 1934 -6.8533560 50.21534 -2.96198419
## Apr 1934 -2.7528761 50.19648 -0.54360561
## May 1934 3.5211867 50.17762 -0.29880990
## Jun 1934 8.9998830 50.15876 0.44135231
## Jul 1934 12.8601452 50.13991 3.49994860
## Aug 1934 11.4833902 50.12105 -1.20443797
## Sep 1934 7.4481630 50.06654 1.68529384
## Oct 1934 0.4621026 50.01204 0.72585893
## Nov 1934 -6.4543306 49.95753 -0.70320331
## Dec 1934 -9.5010869 49.90303 5.39805767
## Jan 1935 -9.3542441 49.84852 -0.49428057
## Feb 1935 -9.8589768 49.81024 2.64873548
## Mar 1935 -6.8533560 49.77196 0.58139808
## Apr 1935 -2.7528761 49.73367 0.11920153
## May 1935 3.5211867 49.69539 -3.21657789
## Jun 1935 8.9998830 49.65711 1.84300919
## Jul 1935 12.8601452 49.63032 2.10953856
## Aug 1935 11.4833902 49.60352 2.91308507
## Sep 1935 7.4481630 49.57673 -0.22489618
## Oct 1935 0.4621026 49.54994 -1.41204416
## Nov 1935 -6.4543306 49.52315 1.13118054
## Dec 1935 -9.5010869 49.47717 -3.57608455
## Jan 1936 -9.3542441 49.43119 -2.77694885
## Feb 1936 -9.8589768 49.38521 -4.52623760
## Mar 1936 -6.8533560 49.33924 1.51412020
## Apr 1936 -2.7528761 49.29326 -2.64038114
## May 1936 3.5211867 49.26631 -0.08749256
## Jun 1936 8.9998830 49.23935 0.36076252
## Jul 1936 12.8601452 49.21240 -2.07254831
## Aug 1936 11.4833902 49.18545 0.43115800
## Sep 1936 7.4481630 49.15850 1.49333654
## Oct 1936 0.4621026 49.16944 -0.03154470
## Nov 1936 -6.4543306 49.18038 -1.12605326
## Dec 1936 -9.5010869 49.19133 1.60976139
## Jan 1937 -9.3542441 49.20227 0.95197682
## Feb 1937 -9.8589768 49.21321 1.64576782
## Mar 1937 -6.8533560 49.25402 -4.00066275
## Apr 1937 -2.7528761 49.29483 0.85804753
## May 1937 3.5211867 49.33564 1.24317495
## Jun 1937 8.9998830 49.37645 0.22366885
## Jul 1937 12.8601452 49.41726 -0.87740315
## Aug 1937 11.4833902 49.45101 0.86560245
## Sep 1937 7.4481630 49.48476 -0.63291971
## Oct 1937 0.4621026 49.51851 0.91939141
## Nov 1937 -6.4543306 49.55226 -1.69792480
## Dec 1937 -9.5010869 49.58600 -2.98491779
## Jan 1938 -9.3542441 49.58897 1.86527399
## Feb 1938 -9.8589768 49.59194 1.46704133
## Mar 1938 -6.8533560 49.59490 4.55845521
## Apr 1938 -2.7528761 49.59787 -0.24499005
## May 1938 3.5211867 49.60083 -0.72201819
## Jun 1938 8.9998830 49.60194 0.39818022
## Jul 1938 12.8601452 49.60304 -2.86318728
## Aug 1938 11.4833902 49.60415 -0.68753765
## Sep 1938 7.4481630 49.60525 -0.05341578
## Oct 1938 0.4621026 49.60636 0.63153937
## Nov 1938 -6.4543306 49.60619 4.64814451
## Dec 1938 -9.5010869 49.60601 -0.90492712
## Jan 1939 -9.3542441 49.60584 -0.85159798
## Feb 1939 -9.8589768 49.60567 1.15330673
## Mar 1939 -6.8533560 49.60550 -0.35214202
## Apr 1939 -2.7528761 49.60300 0.94987333
## May 1939 3.5211867 49.60051 -0.72169420
## Jun 1939 8.9998830 49.59801 -0.59789523
## Jul 1939 12.8601452 49.59552 -1.75566217
## Aug 1939 11.4833902 49.59302 0.72358803
## Sep 1939 7.4481630 49.58770 1.16413231
## Oct 1939 0.4621026 49.58239 -3.34449013
## Nov 1939 -6.4543306 49.57707 3.47726010
## Dec 1939 -9.5010869 49.57175 -2.27066645

plot(stl(nottem, s.window = "periodic", t.window = 50))

上图的四个组成部分分别是时序图、季节效应图、趋势图以及随机波动项，由于这是一个关于平均气温的时间序列，所以季节成分还是很稳定的；另外还需要注意的是，每个图中的y轴尺度都是不同的，因此我们通过图中右侧的灰色长条来指示量级，即每个长条代表的量级一样。

• 使用monthplot()函数查看时序的月度图，该图将不同年份的相同月份分类汇总在一起并合成到同一幅图上，图中的各个横线为各个子序列的均值，从这个图中也可以看到整体的趋势走向与季节变化情况:

par(mfrow=c(1,1))
monthplot(nottem, xlab="", ylab="")

decompose函数

• stl只做相加模型，需考虑相乘模型时，需要进行对数变换 log(Yt) = log(Trend * Seasonal * Irregular)= log(Trend) + log(Seasonal) + log(Irregular)转化为相加模型。
• R中自带的decompose()函数对相加与相乘模型都可以直接进行季节分解
• decompose(x, type = c(“additive”, “multiplicative”), filter = NULL)

# r中自带的AirPassengers数据集就是一个典型的相乘模型
data(AirPassengers)
data <- decompose(AirPassengers,type='multiplicative')
# decompose函数返回的是一个列表,包含季节性成分，趋势成分和不规则成分等信息
data

## $x
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1949 112 118 132 129 121 135 148 148 136 119 104 118
## 1950 115 126 141 135 125 149 170 170 158 133 114 140
## 1951 145 150 178 163 172 178 199 199 184 162 146 166
## 1952 171 180 193 181 183 218 230 242 209 191 172 194
## 1953 196 196 236 235 229 243 264 272 237 211 180 201
## 1954 204 188 235 227 234 264 302 293 259 229 203 229
## 1955 242 233 267 269 270 315 364 347 312 274 237 278
## 1956 284 277 317 313 318 374 413 405 355 306 271 306
## 1957 315 301 356 348 355 422 465 467 404 347 305 336
## 1958 340 318 362 348 363 435 491 505 404 359 310 337
## 1959 360 342 406 396 420 472 548 559 463 407 362 405
## 1960 417 391 419 461 472 535 622 606 508 461 390 432
##
## $seasonal
## Jan Feb Mar Apr May Jun Jul
## 1949 0.9102304 0.8836253 1.0073663 0.9759060 0.9813780 1.1127758 1.2265555
## 1950 0.9102304 0.8836253 1.0073663 0.9759060 0.9813780 1.1127758 1.2265555
## 1951 0.9102304 0.8836253 1.0073663 0.9759060 0.9813780 1.1127758 1.2265555
## 1952 0.9102304 0.8836253 1.0073663 0.9759060 0.9813780 1.1127758 1.2265555
## 1953 0.9102304 0.8836253 1.0073663 0.9759060 0.9813780 1.1127758 1.2265555
## 1954 0.9102304 0.8836253 1.0073663 0.9759060 0.9813780 1.1127758 1.2265555
## 1955 0.9102304 0.8836253 1.0073663 0.9759060 0.9813780 1.1127758 1.2265555
## 1956 0.9102304 0.8836253 1.0073663 0.9759060 0.9813780 1.1127758 1.2265555
## 1957 0.9102304 0.8836253 1.0073663 0.9759060 0.9813780 1.1127758 1.2265555
## 1958 0.9102304 0.8836253 1.0073663 0.9759060 0.9813780 1.1127758 1.2265555
## 1959 0.9102304 0.8836253 1.0073663 0.9759060 0.9813780 1.1127758 1.2265555
## 1960 0.9102304 0.8836253 1.0073663 0.9759060 0.9813780 1.1127758 1.2265555
## Aug Sep Oct Nov Dec
## 1949 1.2199110 1.0604919 0.9217572 0.8011781 0.8988244
## 1950 1.2199110 1.0604919 0.9217572 0.8011781 0.8988244
## 1951 1.2199110 1.0604919 0.9217572 0.8011781 0.8988244
## 1952 1.2199110 1.0604919 0.9217572 0.8011781 0.8988244
## 1953 1.2199110 1.0604919 0.9217572 0.8011781 0.8988244
## 1954 1.2199110 1.0604919 0.9217572 0.8011781 0.8988244
## 1955 1.2199110 1.0604919 0.9217572 0.8011781 0.8988244
## 1956 1.2199110 1.0604919 0.9217572 0.8011781 0.8988244
## 1957 1.2199110 1.0604919 0.9217572 0.8011781 0.8988244
## 1958 1.2199110 1.0604919 0.9217572 0.8011781 0.8988244
## 1959 1.2199110 1.0604919 0.9217572 0.8011781 0.8988244
## 1960 1.2199110 1.0604919 0.9217572 0.8011781 0.8988244
##
## $trend
## Jan Feb Mar Apr May Jun Jul Aug
## 1949 NA NA NA NA NA NA 126.7917 127.2500
## 1950 131.2500 133.0833 134.9167 136.4167 137.4167 138.7500 140.9167 143.1667
## 1951 157.1250 159.5417 161.8333 164.1250 166.6667 169.0833 171.2500 173.5833
## 1952 183.1250 186.2083 189.0417 191.2917 193.5833 195.8333 198.0417 199.7500
## 1953 215.8333 218.5000 220.9167 222.9167 224.0833 224.7083 225.3333 225.3333
## 1954 228.0000 230.4583 232.2500 233.9167 235.6250 237.7500 240.5000 243.9583
## 1955 261.8333 266.6667 271.1250 275.2083 278.5000 281.9583 285.7500 289.3333
## 1956 309.9583 314.4167 318.6250 321.7500 324.5000 327.0833 329.5417 331.8333
## 1957 348.2500 353.0000 357.6250 361.3750 364.5000 367.1667 369.4583 371.2083
## 1958 375.2500 377.9167 379.5000 380.0000 380.7083 380.9583 381.8333 383.6667
## 1959 402.5417 407.1667 411.8750 416.3333 420.5000 425.5000 430.7083 435.1250
## 1960 456.3333 461.3750 465.2083 469.3333 472.7500 475.0417 NA NA
## Sep Oct Nov Dec
## 1949 127.9583 128.5833 129.0000 129.7500
## 1950 145.7083 148.4167 151.5417 154.7083
## 1951 175.4583 176.8333 178.0417 180.1667
## 1952 202.2083 206.2500 210.4167 213.3750
## 1953 224.9583 224.5833 224.4583 225.5417
## 1954 247.1667 250.2500 253.5000 257.1250
## 1955 293.2500 297.1667 301.0000 305.4583
## 1956 334.4583 337.5417 340.5417 344.0833
## 1957 372.1667 372.4167 372.7500 373.6250
## 1958 386.5000 390.3333 394.7083 398.6250
## 1959 437.7083 440.9583 445.8333 450.6250
## 1960 NA NA NA NA
##
## $random
## Jan Feb Mar Apr May Jun Jul
## 1949 NA NA NA NA NA NA 0.9516643
## 1950 0.9626030 1.0714668 1.0374474 1.0140476 0.9269030 0.9650406 0.9835566
## 1951 1.0138446 1.0640180 1.0918541 1.0176651 1.0515825 0.9460444 0.9474041
## 1952 1.0258814 1.0939696 1.0134734 0.9695596 0.9632673 1.0003735 0.9468562
## 1953 0.9976684 1.0151646 1.0604644 1.0802327 1.0413329 0.9718056 0.9551933
## 1954 0.9829785 0.9232032 1.0044417 0.9943899 1.0119479 0.9978740 1.0237753
## 1955 1.0154046 0.9888241 0.9775844 1.0015732 0.9878755 1.0039635 1.0385512
## 1956 1.0066157 0.9970250 0.9876248 0.9968224 0.9985644 1.0275560 1.0217685
## 1957 0.9937293 0.9649918 0.9881769 0.9867637 0.9924177 1.0328601 1.0261250
## 1958 0.9954212 0.9522762 0.9469115 0.9383993 0.9715785 1.0261340 1.0483841
## 1959 0.9825176 0.9505736 0.9785278 0.9746440 1.0177637 0.9968613 1.0373136
## 1960 1.0039279 0.9590794 0.8940857 1.0064948 1.0173588 1.0120790 NA
## Aug Sep Oct Nov Dec
## 1949 0.9534014 1.0022198 1.0040278 1.0062701 1.0118119
## 1950 0.9733720 1.0225047 0.9721928 0.9389527 1.0067914
## 1951 0.9397599 0.9888637 0.9938809 1.0235337 1.0250824
## 1952 0.9931171 0.9746302 1.0046687 1.0202797 1.0115407
## 1953 0.9894989 0.9934337 1.0192680 1.0009392 0.9915039
## 1954 0.9845184 0.9881036 0.9927613 0.9995143 0.9908692
## 1955 0.9831117 1.0032501 1.0003084 0.9827720 1.0125535
## 1956 1.0004765 1.0008730 0.9835071 0.9932761 0.9894251
## 1957 1.0312668 1.0236147 1.0108432 1.0212995 1.0005263
## 1958 1.0789695 0.9856540 0.9977971 0.9802940 0.9405687
## 1959 1.0531001 0.9974447 1.0013371 1.0134608 0.9999192
## 1960 NA NA NA NA NA
##
## $figure
## [1] 0.9102304 0.8836253 1.0073663 0.9759060 0.9813780 1.1127758 1.2265555
## [8] 1.2199110 1.0604919 0.9217572 0.8011781 0.8988244
##
## $type
## [1] "multiplicative"
##
## attr(,"class")
## [1] "decomposed.ts"

plot(data)

在分解季节成分的基础上，如果有需要的话，我们可以对时间序列进行季节因素调整，将这一部分信息从原始数据中去除。

data <- decompose(nottem,type='additive')
data2 <- nottem-data$seasonal
par(mfrow=c(2,1))
plot(nottem)
plot(data2)

上边的图形为原始数据，下边的图形则为去除掉季节成分后的修正数据，此时时序中仅包含趋势成分与随机波动成分。

forecast包等其他操作

• 我们想使用指数平滑法进行预测的话，必须设定时序数据中的预测误差是不相关的，且服从零均值、方差不变的正态分布。
• 自回归移动平均模型（ARIMA）：时间序列的预测值表示为由最近的真实值与最近的预测误差所组成的线性函数

lag函数用于实现多阶滞后

nhtemp

## Time Series:
## Start = 1912
## End = 1971
## Frequency = 1
## [1] 49.9 52.3 49.4 51.1 49.4 47.9 49.8 50.9 49.3 51.9 50.8 49.6 49.3 50.6 48.4
## [16] 50.7 50.9 50.6 51.5 52.8 51.8 51.1 49.8 50.2 50.4 51.6 51.8 50.9 48.8 51.7
## [31] 51.0 50.6 51.7 51.5 52.1 51.3 51.0 54.0 51.4 52.7 53.1 54.6 52.0 52.0 50.9
## [46] 52.6 50.2 52.6 51.6 51.9 50.5 50.9 51.7 51.4 51.7 50.8 51.9 51.8 51.9 53.0

# 滞后2阶
lag(nhtemp,2)

## Time Series:
## Start = 1910
## End = 1969
## Frequency = 1
## [1] 49.9 52.3 49.4 51.1 49.4 47.9 49.8 50.9 49.3 51.9 50.8 49.6 49.3 50.6 48.4
## [16] 50.7 50.9 50.6 51.5 52.8 51.8 51.1 49.8 50.2 50.4 51.6 51.8 50.9 48.8 51.7
## [31] 51.0 50.6 51.7 51.5 52.1 51.3 51.0 54.0 51.4 52.7 53.1 54.6 52.0 52.0 50.9
## [46] 52.6 50.2 52.6 51.6 51.9 50.5 50.9 51.7 51.4 51.7 50.8 51.9 51.8 51.9 53.0

ACF和PACF

自相关度量的是时间序列中各个预测值之间的相关性，为一系列观测值与时期之前的观测值之间的相关性；这样，就是一阶滞后序列和0阶滞后序列间的相关性，是二阶滞后序列和0阶滞后序列之间的相关性，以此类推。根据这些相关性, , , 所绘制的图即为自相关函数图（ACF图），它可用于为模型选择合适的参数，并评估最终模型的拟合效果。

需要明确的是，此时我们所得到的自相关系数，实际上并不是与之间单纯的相关关系。因为同时还会受到中间k-1个随机变量, , 的影响，而这k-1个随机变量又都和具有相关关系，所以自相关系数里实际掺杂了其他变量对与的影响。

在R中，自相关与偏自相关可使用forecast包中的函数来实现：

library(forecast)

## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo

par(mfrow=c(1,2))
a <- Acf(nhtemp)
a

##
## Autocorrelations of series 'nhtemp', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 0.315 0.375 0.264 0.241 0.106 0.292 0.161 0.259 0.068 0.165
## 11 12 13 14 15 16 17
## 0.077 0.107 -0.094 0.051 0.011 0.126 0.073

b <- Pacf(nhtemp)

b

##
## Partial autocorrelations of series 'nhtemp', by lag
##
## 1 2 3 4 5 6 7 8 9 10 11
## 0.315 0.307 0.105 0.065 -0.077 0.215 0.032 0.109 -0.139 0.026 0.020
## 12 13 14 15 16 17
## -0.001 -0.210 -0.016 0.101 0.154 0.038

从上述图形中我们可以看到，自相关图中的滞后一阶、二阶、三阶和六阶均超出了置信边界；偏相关图除了滞后一阶与二阶外都在置信区间内。这些信息将会为我们后面的模型参数选取提供指导。

检验

• 时间序列的平稳性一般可以通过图形来直接判断，观察是否有明显的趋势向和不规则波动项，
• 也可使用单位根（ADF）统计检验来验证平稳性，原假设为该序列存在单位根，如果经过原假设不能被拒绝的话，那么我们就可以认为该序列为非平稳时间序列。
• H0:非平稳； H1:平稳

R中的ADF检验可使用tseries包中的adf.test()函数来实现：

library(tseries)
adf.test(nhtemp)

##
## Augmented Dickey-Fuller Test
##
## data: nhtemp
## Dickey-Fuller = -3.2773, Lag order = 3, p-value = 0.08376
## alternative hypothesis: stationary

# 备择假设为'stationary',即平稳时序,零假设也就自然为非平稳时序
# 已知P值为0.08>0.05,

差分

使用R中自带的diff()函数与forecast包中的ndiffs()函数均可以进行差分，diff函数会返回差分后的数据，ndiffs函数可以帮助我们最优的d值。

ndiffs(nhtemp)

## [1] 1

#ndiffs(test=c("kpss","adf","pp"))
diff(nhtemp,1)

## Time Series:
## Start = 1913
## End = 1971
## Frequency = 1
## [1] 2.4 -2.9 1.7 -1.7 -1.5 1.9 1.1 -1.6 2.6 -1.1 -1.2 -0.3 1.3 -2.2 2.3
## [16] 0.2 -0.3 0.9 1.3 -1.0 -0.7 -1.3 0.4 0.2 1.2 0.2 -0.9 -2.1 2.9 -0.7
## [31] -0.4 1.1 -0.2 0.6 -0.8 -0.3 3.0 -2.6 1.3 0.4 1.5 -2.6 0.0 -1.1 1.7
## [46] -2.4 2.4 -1.0 0.3 -1.4 0.4 0.8 -0.3 0.3 -0.9 1.1 -0.1 0.1 1.1

#par(mfrow=c(1,2))
plot(nhtemp)

plot(diff(nhtemp,1))

arima函数

平稳性是模型中非常重要的一个假设，只有在满足这一假设并结合前面我们所学习的相关概念，我们才能进一步拟合出具有自回归项（AR）、移动平均项（MA）的模型，并在此基础上进行差分，得到最终的ARIMA模型。

ARMA的全称是自回归移动平均模型，进一步的，ARIMA(p,d,q)又称为差分自回归移动平均模型

一个完整的模型的应用由以下步骤构成： - 1. 时间序列的平稳性检验（若不平稳，则进行差分） - 2. 寻找最为合适的参数与进行模型构建 - 3. 从统计假设和预测准确性等角度评估模型 - 4. 将得到的模型应用于预测

同样，R当中实现模型需要加载forecast包，使用其中的arima函数：

library(forecast)
data(nhtemp)
library(tseries)
dnhtemp <- diff(nhtemp)
adf.test(dnhtemp)

## Warning in adf.test(dnhtemp): p-value smaller than printed p-value

##
## Augmented Dickey-Fuller Test
##
## data: dnhtemp
## Dickey-Fuller = -4.6366, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationary

par(mfrow=c(1,2))
acf(dnhtemp) ##q: q阶后减小到零
Pacf(dnhtemp) ##p: p阶后减小到零

除了这一种方法外，forecast包中的auto.arima()函数可以帮助我们自行选取最优参数值，这比观察ACF图自行确定要便捷的多。

model_fit <- auto.arima(nhtemp)
model_fit

## Series: nhtemp
## ARIMA(0,1,1)
##
## Coefficients:
## ma1
## -0.7983
## s.e. 0.0956
##
## sigma^2 estimated as 1.313: log likelihood=-91.76
## AIC=187.52 AICc=187.73 BIC=191.67

模型选定的参数值为(0,1,1)

model_fit1 <- auto.arima(dnhtemp)
model_fit1

## Series: dnhtemp
## ARIMA(0,0,1) with zero mean
##
## Coefficients:
## ma1
## -0.7983
## s.e. 0.0956
##
## sigma^2 estimated as 1.313: log likelihood=-91.76
## AIC=187.52 AICc=187.73 BIC=191.67

模型的构建、评估与预测

model_fit <- arima(nhtemp,order=c(0,1,1))
model_fit

##
## Call:
## arima(x = nhtemp, order = c(0, 1, 1))
##
## Coefficients:
## ma1
## -0.7983
## s.e. 0.0956
##
## sigma^2 estimated as 1.291: log likelihood = -91.76, aic = 187.52

accuracy(model_fit)

## ME RMSE MAE MPE MAPE MASE
## Training set 0.1285716 1.126722 0.8952077 0.2080785 1.749865 0.7513123
## ACF1
## Training set -0.008258617

一般来说，一个模型如果合适，那模型的残差应该满足均值为0的正态分布，并且对于任意 的滞后阶数，残差自相关系数都应该为零。换句话说，模型的残差应该满足独立正态分布（即残差间没有关联）

这一点我们可以使用QQ图进行验证:

qqnorm(model_fit$residuals)
qqline(model_fit$residuals)

Box.test函数可以检验残差的自相关系数是否都为零

Box.test(model_fit$residuals,type='Ljung-Box')

##
## Box-Ljung test
##
## data: model_fit$residuals
## X-squared = 0.0043004, df = 1, p-value = 0.9477

# P值为0.95,说明模型的残差没有通过显著性检验(即可以认为残差的自相关系数为零),ARIMA模型能较好地拟合本数据

以下值分别为预测值及80%,95%的置信区间值

forecast(model_fit,5)

## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 1972 51.90008 50.44396 53.35619 49.67314 54.12701
## 1973 51.90008 50.41463 53.38552 49.62828 54.17187
## 1974 51.90008 50.38586 53.41429 49.58429 54.21586
## 1975 51.90008 50.35764 53.44252 49.54112 54.25904
## 1976 51.90008 50.32991 53.47024 49.49872 54.30143

plot(forecast(model_fit,5))

时间序列分析作图

ggplot

library(ggplot2)
# 查询数据集类型
class(Nile)

## [1] "ts"

# 对ts类型的数据进行df类型转换
a <- data.frame(Time=c(time(Nile)),Nile=c(Nile))
p <- ggplot(a,aes(x=Time,y=Nile))
p + geom_line(colour = 'green') + xlab('The Time Series of Date') + ylab('The Time Series of Nile')

ggfortify包

ggfortify包对于Time Series图表的绘制将使用更为简洁的代码，并且更容易实现不同需求下的Time Series图表的绘制。

##单变量
library(ggfortify)

## Registered S3 methods overwritten by 'ggfortify':
## method from
## autoplot.Arima forecast
## autoplot.acf forecast
## autoplot.ar forecast
## autoplot.bats forecast
## autoplot.decomposed.ts forecast
## autoplot.ets forecast
## autoplot.forecast forecast
## autoplot.stl forecast
## autoplot.ts forecast
## fitted.ar forecast
## fortify.ts forecast
## residuals.ar forecast

autoplot(Nile,ts.colour = 'green') + xlab("Date of Nile") + ylab("Number of Nile") + ggtitle("The Time Series of Nile")

# ts.linetype可以改变线条形状
##多变量
data <- ts(matrix(rnorm(180), 60, 3), start = c(2012, 1), frequency = 12,names = c('a','b','c'))
ts.data <- ts(data,start=c(1996,10),freq=12)
autoplot(ts.data,ts.colour = 'green') + ggtitle("The Time Series")

## Warning: Ignoring unknown parameters: ts.colour

autoplot(ts.data,facets=F) + ggtitle("The Time Series")

使用移动平均进行平滑处理

library(forecast)
a <- autoplot(ma(Nile,3))
b <- autoplot(ma(Nile,7))
c <- autoplot(ma(Nile,15))
d <- autoplot(Nile)
library(gridExtra)
grid.arrange(d,a,b,c,ncol=2)

## Warning: Removed 2 row(s) containing missing values (geom_path).

## Warning: Removed 6 row(s) containing missing values (geom_path).

## Warning: Removed 14 row(s) containing missing values (geom_path).

##从图像来看，随着k的增大，图像变得越来越平滑。

分解时间序列

# stl:Decompose a time series into seasonal, trend and irregular components
autoplot(stl(AirPassengers, s.window = 'periodic'), ts.colour = 'green')

Plotting with changepoint package

The changepoint package provides a simple approach for identifying shifts in mean and/or variance in a time series.

library(changepoint)

## Loading required package: zoo

##
## Attaching package: 'zoo'

## The following object is masked from 'package:timeSeries':
##
## time<-

## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric

## Successfully loaded changepoint package version 2.2.2
## NOTE: Predefined penalty values changed in version 2.2. Previous penalty values with a postfix 1 i.e. SIC1 are now without i.e. SIC and previous penalties without a postfix i.e. SIC are now with a postfix 0 i.e. SIC0. See NEWS and help files for further details.

autoplot(cpt.meanvar(Nile),colour='green') + ggtitle("The Time Series of Nile")

## Warning: `filter_()` is deprecated as of dplyr 0.7.0.
## Please use `filter()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

Plotting with strucchange package

strucchange is an R package for detecting jumps in data.

library(strucchange)

## Loading required package: sandwich

autoplot(breakpoints(Nile ~ 1),colour = 'green')

Time Series的指数预测模型

• 单指数模型（simple exponential model）：拟合的是只有常数水平项和时间点处随机项的时间序列，并不考虑时间序列的趋势项和季节效应。
• 双指数模型（double exponential model）：拟合的是有水平项和趋势项的时序，也叫Holt指数平滑。
• 三指数模型（triple exponential model)：拟合的是有水平项、趋势项以及 季节效应的时序。也叫Holt-Winters指数平滑（Holt-Winters exponential smoothing）

使用forecast包ets函数绘制单指数平滑观测图（指定model=‘ANN’)

AirPassengers_fit <- ets(AirPassengers,model='ANN')
plot(forecast(AirPassengers_fit,1))

对nottem数据集进行H-W指数平滑预测（model=‘AAA’)

• alpha：控制水平项的指数型下降
• beta：控制斜率的指数型下降
• gamma：光滑参数控制季节项的指数下降
• 以上三个参数值越大,意味着越近的观测值的权重越大
• 以上三个参数取值均在[0,1]之间

fit <- ets(nottem,model='AAA')
fit

## ETS(A,A,A)
##
## Call:
## ets(y = nottem, model = "AAA")
##
## Smoothing parameters:
## alpha = 0.0894
## beta = 1e-04
## gamma = 1e-04
##
## Initial states:
## l = 48.9625
## b = 0.0019
## s = -9.6666 -6.4508 0.7628 7.4973 11.4723 12.9467
## 8.9557 3.4648 -2.8831 -6.8772 -9.7282 -9.4937
##
## sigma: 2.3405
##
## AIC AICc BIC
## 1740.964 1743.721 1800.135

# 对接下来五个月的数据进行预测
nottem_forecast <- forecast(ets(nottem),5)
# 后四列数据分别为80%和95%区间的上下界
nottem_forecast

## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 1940 40.17228 37.19869 43.14587 35.62456 44.72000
## Feb 1940 39.76945 36.79358 42.74532 35.21825 44.32065
## Mar 1940 42.71854 39.74039 45.69669 38.16385 47.27323
## Apr 1940 46.77634 43.79591 49.75676 42.21817 51.33450
## May 1940 52.94793 49.96523 55.93063 48.38628 57.50957

plot(nottem_forecast)

autoplot函数

autoplot函数也可以理解其他的时间序列类型，支持的R包包括以下： - xts包：xts - zoo包：zooreg - tseries包：irts - timeSeries包：timeSeries

library(zoo)
autoplot(as.zooreg(AirPassengers),colour = 'green')

library(xts)
autoplot(as.xts(AirPassengers),colour = 'green') # 等价于as.zooreg()