【免責聲明：本文用于教學】

時間序列分析

基礎操作

數據輸入

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()

總結

以上是生活随笔為你收集整理的语言时间序列年月日_R语言系列 时间序列分析的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。