R 時間序列分析（一）

時間序列ts()介紹

這裡先介紹R中內建時間序列的物件格式 Time series。
我們來自己建立一組時間序列的資料，ts()建立，使用格式如下：
ts(data,frequency=12,start=c(year,month))
其中gm表示時間數據，frequency表示時間的單位(進位法)，例如季：4進位、月：12進位、小時：24、日：365…

這邊看一下官方文件
frequency: the number of observations per unit of time.
The value of argument frequency is used when the series is sampled an integral number of times in each unit time interval. For example, one could use a value of 7 for frequency when the data are sampled daily, and the natural time period is a week, or 12 when the data are sampled monthly and the natural time period is a year. Values of 4 and 12 are assumed in (e.g.) print methods to imply a quarterly and monthly series respectively.

start則表示時間序列開始的時間

例如：我們自訂一個模擬一個辦公室7月份一周的室溫資料，儀器紀錄時間是小時，並且來計算看看acf偏相關分析

temp_office <- ts(rnorm(24*7,mean=30,sd=3),frequency=1,start=c(1,1)) ;head(temp_office,20)

## Time Series:
## Start = 1 
## End = 20 
## Frequency = 1 
##  [1] 28.90558 28.41619 29.01357 30.80248 28.94694 25.43002 25.81873 34.62282
##  [9] 25.36138 24.15885 29.81578 30.61072 36.20670 27.06386 27.14283 27.53260
## [17] 31.88450 34.73242 32.01088 32.33917

acf(temp_office, plot=T)

acf()稍後就會介紹，這邊為了方便觀察資料隨著小時的自相關性，我直接將frequency設定為1，這樣acf的lag值較直覺（lag=1對應frquency一個週期）。
因為溫度是隨機產生的，因此看不到相關係數高的遲滯係數(lag)。

時間序列觀察

接下來正是進入本次練習，我們採用forecast套件包require(forecast)裡面的內建資料集wineind，該資料是從1980年01月到1994年08月，是葡萄酒商紀錄紅酒銷售的總量。另外在未來的文章也會用forecast套件包中的函式進行建模。

首先透過 str() 可以了解到這組資料是時間序列的資料。

str(wineind);wineind;summary(wineind)

##  Time-Series [1:176] from 1980 to 1995: 15136 16733 20016 17708 18019 ...

##        Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
## 1980 15136 16733 20016 17708 18019 19227 22893 23739 21133 22591 26786 29740
## 1981 15028 17977 20008 21354 19498 22125 25817 28779 20960 22254 27392 29945
## 1982 16933 17892 20533 23569 22417 22084 26580 27454 24081 23451 28991 31386
## 1983 16896 20045 23471 21747 25621 23859 25500 30998 24475 23145 29701 34365
## 1984 17556 22077 25702 22214 26886 23191 27831 35406 23195 25110 30009 36242
## 1985 18450 21845 26488 22394 28057 25451 24872 33424 24052 28449 33533 37351
## 1986 19969 21701 26249 24493 24603 26485 30723 34569 26689 26157 32064 38870
## 1987 21337 19419 23166 28286 24570 24001 33151 24878 26804 28967 33311 40226
## 1988 20504 23060 23562 27562 23940 24584 34303 25517 23494 29095 32903 34379
## 1989 16991 21109 23740 25552 21752 20294 29009 25500 24166 26960 31222 38641
## 1990 14672 17543 25453 32683 22449 22316 27595 25451 25421 25288 32568 35110
## 1991 16052 22146 21198 19543 22084 23816 29961 26773 26635 26972 30207 38687
## 1992 16974 21697 24179 23757 25013 24019 30345 24488 25156 25650 30923 37240
## 1993 17466 19463 24352 26805 25236 24735 29356 31234 22724 28496 32857 37198
## 1994 13652 22784 23565 26323 23779 27549 29660 23356

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   13652   22115   24669   25392   28461   40226

我們可以用start()及end()取出時間序列的開始和結束。

start(wineind);end(wineind);

## [1] 1980    1

## [1] 1994    8

end(wineind)[1] #年分

## [1] 1994

end(wineind)[2] #月份

## [1] 8

另外可以用window()來取出某一區間的資料。

window(wineind,1990,end(wineind)) #抓1990之後的資料

##        Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
## 1990 14672 17543 25453 32683 22449 22316 27595 25451 25421 25288 32568 35110
## 1991 16052 22146 21198 19543 22084 23816 29961 26773 26635 26972 30207 38687
## 1992 16974 21697 24179 23757 25013 24019 30345 24488 25156 25650 30923 37240
## 1993 17466 19463 24352 26805 25236 24735 29356 31234 22724 28496 32857 37198
## 1994 13652 22784 23565 26323 23779 27549 29660 23356

用time()取出時間欄位，這裡借用zoo函式庫幫我們轉換時間；另外可以用cycle()取出月份（精準地說應該是看你的時間序列週期單位）。

time(wineind) %>% zoo::as.yearmon() %>% head()

## [1] "一月 1980" "二月 1980" "三月 1980" "四月 1980" "五月 1980" "六月 1980"

#單單取出年分
x <- zoo::as.yearmon(time(wineind)) %>% format.Date("%Y") ;head(x)

## [1] "1980" "1980" "1980" "1980" "1980" "1980"

#取出月份
m <- cycle(wineind) ; head(m)

##      Jan Feb Mar Apr May Jun
## 1980   1   2   3   4   5   6

折線圖

接下來我們來觀察時序圖。

plot.ts(wineind) #用折線圖看銷售量時間序列

觀察可發現，我們的數據具有週期性變化，且隱約有長期增長的趨勢，並非一直維持在某一個常數上下波動，初步先判斷該序列爲非定態時間序列。
舉例來說，像股票大盤指數的價格會上下震盪，但是長期來看逐漸上升，這也是一種非定態時間序列。
而定態又可分為弱定態與嚴格定態等。
定態或非定態會影響我們對於一組數據的判斷，衡量到底這時間序列是否是可預測的。

ADF檢定

這時候我們可以先用一個常見的ADF檢定（Augmented Dickey-Fuller test）來驗證是否為定態。
引入tseries包中的adf.test()函數進行單位根檢驗

adf.test(wineind,k=0)

## Warning in adf.test(wineind, k = 0): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  wineind
## Dickey-Fuller = -11.325, Lag order = 0, p-value = 0.01
## alternative hypothesis: stationary

但結果顯示P值爲0.01<0.05，接受H1假設，檢測認爲該序列為定態(stationary)。
不過我們從觀察中認為數據有長期趨勢，故稍後還是可以嘗試當作為非定態數據進行差分處理。
這樣對我們之後要使用自相關與偏相關分析圖來判斷該用甚麼預測模型(AR,MA,ARIMA)時，可能會有所幫助。

推薦參考：
陳旭昇老師-單根與隨機趨勢
 Edward-資料的時間序列性質檢驗

自相關分析與偏相關分析

同一個時間序列在任意兩個不同時刻的取值之間的相關程度，自相關函式是描述隨機訊號x(t)在任意兩個不同時刻t1，t2的取值之間的相關程度。
自己本身的相關係數一定最高為1，因此lag=0時，固定會是1。

plot(wineind) #用折線圖看銷售量時間序列

acf(wineind,lag.max = 100, plot=T)

#lag.max最大遲滯階數，設定前後自相關計算的最大間隔
#plot預設為T
#除了相關係數外，還可以計算共變異數及偏自相關係數PCAF 
#透過設定 type = c("correlation", "covariance", "partial")

在此使用acf()查看wineind的自相關性(Autocorrelation)
從圖上來看，往後移動一個距離，相關係數約為0.2
直到lag=1時（因為我們時序是12進位的，因此lag=1代表12個月），可看到cor約等於0.8，可預期紅酒的銷售量有12個月的年度循環週期。
在此記錄下自相關顯著處，即Pearson相關係數約接近1(高度正相關)或-1(高度負相關)部分。

而偏自相關分析PCAF，可以直接使用pcaf()直接計算。
自相關分析會受到區間的取值影響，而偏相關的計算可以移除這一點，通常偏相關係數會小於對應的自相關係數。

在此若想知道相關分析原理，推薦參考：
知乎：自相關與偏相關係數
 AI Tech yuan：自相關與偏相關係數

我們一開始提到過，若判斷為非定態序列，可先進行差分處理，變成定態序列，這使用ndiffs()函數判斷差分次數。

ndiffs(wineind)

## [1] 1

建議出結果為一階差分，因此我們來試著對函式進行一階差分，並且來畫畫看做完差分後的acf和pacf圖。

winediff<-diff(wineind)
par(mfrow=c(1,2))
plot.ts(wineind,main="原始Wineind資料")
plot.ts(winediff,main="一階差分後的Wineind資料")

par(mfrow=c(1,2))
acf(wineind);acf(winediff)  #繪製自相關圖

par(mfrow=c(1,2))
pacf(wineind);pacf(winediff)  #繪製偏自相關圖

以上我們看了單一時間序列的自相關分析與偏相關分析。
從分析中我們預期wineind資料有年週期性的變化，之後可以納入預測模型的考量。接下來的文章再繼續介紹時間序列建模預測的相關處理。

周期製圖比較

我們在自相關分析中看到以12個月的相關係數蠻高的，另外當1、4、6、8亦有蠻高的相關係數。這跟直覺蠻符合的，銷售的淡旺季會有年度變化，如同業界常見年度同期比較；又或者跟前後個月有近似表現（環比），也可以是以季、半年的尺度做討論。

這裡先用完整的一個週期（12個月）將同期比較的圖繪出。

#按月份畫折線圖

y <- zoo::as.yearmon(time(wineind)) %>% format.Date("%Y") %>% as.numeric() #年分
m <- zoo::as.yearmon(time(wineind)) %>% format.Date("%m") %>% as.numeric() #月份
wine_df <- data.frame(sales=as.numeric(wineind),year=y,month=m) #轉換成dataframe
duration=max(y)-min(y)+1 #計算有多少年，方便用迴圈分年份繪圖

wine_df %>%
  ggplot( aes(x=month, y=sales, group=year, color=year)) +
    geom_line()+
    geom_point()+ 
    scale_x_continuous(breaks=wine_df$month, labels = wine_df$month)

#另外提供另外的繪圖方法，如果不想換成dataframe來畫的話，這裡參考"實戰R語言預測分析"書中的內容。
plot(1:12,ylim=range(wineind)+c(-100,100),col='white', xlab="月份", ylab="銷量")
for(i in 1:duration) {
points(1:12,wineind[(12*(i-1)+1):(12*i)])
lines(1:12,wineind[(12*(i-1)+1):(12*i)], lty=2 ,col=hcl.colors(duration)[i]) #用hcl.colors自動生成顏色
}

從圖上很明顯可以看出來，銷售數據有隨著月份的年週期性變化，而且隨著年份似乎有漸增的趨勢（非定態資料）
關於趨勢，我們之後再探討用移動平均Moving Average法來分析。
我們也好奇，是不是有跟著季產生的銷售變化，在自相關中4的相關係數蠻高，因此我們修改一下數據框，畫出來比較。

#按月份畫折線圖
q <- zoo::as.yearqtr(time(wineind)) %>% as.character()  #年與季
qm <- zoo::as.yearmon(time(wineind)) %>% format.Date("%m") %>%  as.numeric()%%4+1  #算出該月是在各季第幾個月
wine_df <- data.frame(sales=as.numeric(wineind),year=y,month=m,quarter=q,qutermonth=qm) #轉換成dataframe

wine_df %>%
  ggplot( aes(x=qutermonth, y=sales, group=quarter, color=quarter)) +
    geom_line()+
    geom_point()+ 
    scale_x_continuous(breaks=wine_df$qutermonth, labels = wine_df$qutermonth)

感覺每季的第一個月都是銷售峰值，圖說明資料也有季變化。

一維線性回歸

其他回歸方法在後續的文章介紹，這裡先一氣呵成的用簡單的線性回歸建模一下吧！
我們想用過去的資料來預測未來的銷量，在這裡的例子，就是1981年01月的銷量可以從1980的歷史銷量推算出。來整理一下資料

target_ind <-  wineind[13:length(wineind)] #要被預測的銷量，從第二年開始預測
#用過去1,4,6,8,12個月的值來預測當月份
dataset <- data.frame(target_ind,month=m[13:length(wineind)])
column_name <- c()
for (i in c(1,4,6,8,12)){
  dataset <-  cbind(dataset,col=wineind[(13-i):(length(wineind)-i)])
  column_name <- c(column_name,paste0('R',i))
}
dataset <- setNames(dataset,c('target_ind','month',column_name))

整完的dataset資料如下：

target_ind	month	R1	R4	R6	R8	R12
15028	1	29740	21133	22893	18019	15136
17977	2	15028	22591	23739	19227	16733
20008	3	17977	26786	21133	22893	20016
21354	4	20008	29740	22591	23739	17708
19498	5	21354	15028	26786	21133	18019
22125	6	19498	17977	29740	22591	19227
25817	7	22125	20008	15028	26786	22893
28779	8	25817	21354	17977	29740	23739
20960	9	28779	19498	20008	15028	21133
22254	10	20960	22125	21354	17977	22591
27392	11	22254	25817	19498	20008	26786
29945	12	27392	28779	22125	21354	29740
16933	1	29945	20960	25817	19498	15028
17892	2	16933	22254	28779	22125	17977
20533	3	17892	27392	20960	25817	20008
23569	4	20533	29945	22254	28779	21354
22417	5	23569	16933	27392	20960	19498
22084	6	22417	17892	29945	22254	22125
26580	7	22084	20533	16933	27392	25817
27454	8	26580	23569	17892	29945	28779
24081	9	27454	22417	20533	16933	20960
23451	10	24081	22084	23569	17892	22254
28991	11	23451	26580	22417	20533	27392
31386	12	28991	27454	22084	23569	29945
16896	1	31386	24081	26580	22417	16933
20045	2	16896	23451	27454	22084	17892
23471	3	20045	28991	24081	26580	20533
21747	4	23471	31386	23451	27454	23569
25621	5	21747	16896	28991	24081	22417
23859	6	25621	20045	31386	23451	22084
25500	7	23859	23471	16896	28991	26580
30998	8	25500	21747	20045	31386	27454
24475	9	30998	25621	23471	16896	24081
23145	10	24475	23859	21747	20045	23451
29701	11	23145	25500	25621	23471	28991
34365	12	29701	30998	23859	21747	31386
17556	1	34365	24475	25500	25621	16896
22077	2	17556	23145	30998	23859	20045
25702	3	22077	29701	24475	25500	23471
22214	4	25702	34365	23145	30998	21747
26886	5	22214	17556	29701	24475	25621
23191	6	26886	22077	34365	23145	23859
27831	7	23191	25702	17556	29701	25500
35406	8	27831	22214	22077	34365	30998
23195	9	35406	26886	25702	17556	24475
25110	10	23195	23191	22214	22077	23145
30009	11	25110	27831	26886	25702	29701
36242	12	30009	35406	23191	22214	34365
18450	1	36242	23195	27831	26886	17556
21845	2	18450	25110	35406	23191	22077
26488	3	21845	30009	23195	27831	25702
22394	4	26488	36242	25110	35406	22214
28057	5	22394	18450	30009	23195	26886
25451	6	28057	21845	36242	25110	23191
24872	7	25451	26488	18450	30009	27831
33424	8	24872	22394	21845	36242	35406
24052	9	33424	28057	26488	18450	23195
28449	10	24052	25451	22394	21845	25110
33533	11	28449	24872	28057	26488	30009
37351	12	33533	33424	25451	22394	36242
19969	1	37351	24052	24872	28057	18450
21701	2	19969	28449	33424	25451	21845
26249	3	21701	33533	24052	24872	26488
24493	4	26249	37351	28449	33424	22394
24603	5	24493	19969	33533	24052	28057
26485	6	24603	21701	37351	28449	25451
30723	7	26485	26249	19969	33533	24872
34569	8	30723	24493	21701	37351	33424
26689	9	34569	24603	26249	19969	24052
26157	10	26689	26485	24493	21701	28449
32064	11	26157	30723	24603	26249	33533
38870	12	32064	34569	26485	24493	37351
21337	1	38870	26689	30723	24603	19969
19419	2	21337	26157	34569	26485	21701
23166	3	19419	32064	26689	30723	26249
28286	4	23166	38870	26157	34569	24493
24570	5	28286	21337	32064	26689	24603
24001	6	24570	19419	38870	26157	26485
33151	7	24001	23166	21337	32064	30723
24878	8	33151	28286	19419	38870	34569
26804	9	24878	24570	23166	21337	26689
28967	10	26804	24001	28286	19419	26157
33311	11	28967	33151	24570	23166	32064
40226	12	33311	24878	24001	28286	38870
20504	1	40226	26804	33151	24570	21337
23060	2	20504	28967	24878	24001	19419
23562	3	23060	33311	26804	33151	23166
27562	4	23562	40226	28967	24878	28286
23940	5	27562	20504	33311	26804	24570
24584	6	23940	23060	40226	28967	24001
34303	7	24584	23562	20504	33311	33151
25517	8	34303	27562	23060	40226	24878
23494	9	25517	23940	23562	20504	26804
29095	10	23494	24584	27562	23060	28967
32903	11	29095	34303	23940	23562	33311
34379	12	32903	25517	24584	27562	40226
16991	1	34379	23494	34303	23940	20504
21109	2	16991	29095	25517	24584	23060
23740	3	21109	32903	23494	34303	23562
25552	4	23740	34379	29095	25517	27562
21752	5	25552	16991	32903	23494	23940
20294	6	21752	21109	34379	29095	24584
29009	7	20294	23740	16991	32903	34303
25500	8	29009	25552	21109	34379	25517
24166	9	25500	21752	23740	16991	23494
26960	10	24166	20294	25552	21109	29095
31222	11	26960	29009	21752	23740	32903
38641	12	31222	25500	20294	25552	34379
14672	1	38641	24166	29009	21752	16991
17543	2	14672	26960	25500	20294	21109
25453	3	17543	31222	24166	29009	23740
32683	4	25453	38641	26960	25500	25552
22449	5	32683	14672	31222	24166	21752
22316	6	22449	17543	38641	26960	20294
27595	7	22316	25453	14672	31222	29009
25451	8	27595	32683	17543	38641	25500
25421	9	25451	22449	25453	14672	24166
25288	10	25421	22316	32683	17543	26960
32568	11	25288	27595	22449	25453	31222
35110	12	32568	25451	22316	32683	38641
16052	1	35110	25421	27595	22449	14672
22146	2	16052	25288	25451	22316	17543
21198	3	22146	32568	25421	27595	25453
19543	4	21198	35110	25288	25451	32683
22084	5	19543	16052	32568	25421	22449
23816	6	22084	22146	35110	25288	22316
29961	7	23816	21198	16052	32568	27595
26773	8	29961	19543	22146	35110	25451
26635	9	26773	22084	21198	16052	25421
26972	10	26635	23816	19543	22146	25288
30207	11	26972	29961	22084	21198	32568
38687	12	30207	26773	23816	19543	35110
16974	1	38687	26635	29961	22084	16052
21697	2	16974	26972	26773	23816	22146
24179	3	21697	30207	26635	29961	21198
23757	4	24179	38687	26972	26773	19543
25013	5	23757	16974	30207	26635	22084
24019	6	25013	21697	38687	26972	23816
30345	7	24019	24179	16974	30207	29961
24488	8	30345	23757	21697	38687	26773
25156	9	24488	25013	24179	16974	26635
25650	10	25156	24019	23757	21697	26972
30923	11	25650	30345	25013	24179	30207
37240	12	30923	24488	24019	23757	38687
17466	1	37240	25156	30345	25013	16974
19463	2	17466	25650	24488	24019	21697
24352	3	19463	30923	25156	30345	24179
26805	4	24352	37240	25650	24488	23757
25236	5	26805	17466	30923	25156	25013
24735	6	25236	19463	37240	25650	24019
29356	7	24735	24352	17466	30923	30345
31234	8	29356	26805	19463	37240	24488
22724	9	31234	25236	24352	17466	25156
28496	10	22724	24735	26805	19463	25650
32857	11	28496	29356	25236	24352	30923
37198	12	32857	31234	24735	26805	37240
13652	1	37198	22724	29356	25236	17466
22784	2	13652	28496	31234	24735	19463
23565	3	22784	32857	22724	29356	24352
26323	4	23565	37198	28496	31234	26805
23779	5	26323	13652	32857	22724	25236
27549	6	23779	22784	37198	28496	24735
29660	7	27549	23565	13652	32857	29356
23356	8	29660	26323	22784	37198	31234

畫出散佈圖矩陣：

plot(dataset)

進階版本的散布矩陣，使用到car套件，自帶線性回歸線

scatterplotMatrix(dataset,regLine = list(method=lm, lty=1, lwd=2, col='red') )

當我們仔細看target_ind與R12的散佈圖時，發現有些明顯的離群值。這些離群值會影響建模的效果。

plot(dataset[,c('R12','target_ind')])
# 畫出回歸趨勢線
lm.model <- lm(target_ind~R12, dataset)    # 建立一個線性回歸
abline(lm.model,lwd=2,col='red')

庫克距離，去除異常值

我們在簡單線性回歸後看到有離群值，恐影響模型擬合，如果某個觀測值對應的期望值是異常值的話，這樣的值被認爲槓桿值很大。
這邊採用庫克距離來衡量異常值。

庫克距離：
當通過計算觀測值的槓桿值之後，發現具有較大槓桿值的那些觀測點，對模型穩定性有“潛在威脅”。
此時就可以計算庫克距離衡量一個觀測值對模型的影響大小，白話來說就是比較每個觀測值被移除前後，模型變化。
推薦閱讀：
庫克距離

lm.model <- lm(target_ind~R12, dataset) 
cook <- cooks.distance(lm.model) #計算庫克距離
head(cook) #看一下cook計算的結果

##            1            2            3            4            5            6 
## 7.400274e-03 1.037410e-04 1.060328e-03 7.571059e-03 3.470601e-05 3.286490e-03

plot(cook)
abline(h=0.10,lty=2,col='red')

dataset <- dataset[which(cook<0.10),] # remove 高槓桿的異常值

建立模型

我們將dataset分成訓練組與測試組，這裡就按照比例，前7成是訓練組，後3成是測試組。

#分離訓練與測試組
train_size = round(nrow(dataset)*0.7)
test_size = nrow(dataset)-train_size
train <- head(dataset,train_size)
test <- tail(dataset,test_size)

#建立模型
lm.fit <- lm(target_ind~.,data=train)
summary(lm.fit)

## 
## Call:
## lm(formula = target_ind ~ ., data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4825.9 -1542.3  -126.4  1495.8  6797.8 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.917e+03  2.123e+03   0.903 0.368605    
## month        4.170e+02  1.085e+02   3.844 0.000207 ***
## R1          -3.040e-04  3.814e-02  -0.008 0.993655    
## R4           6.858e-02  4.080e-02   1.681 0.095752 .  
## R6          -2.909e-02  4.303e-02  -0.676 0.500429    
## R8           1.298e-01  4.971e-02   2.612 0.010320 *  
## R12          6.705e-01  7.118e-02   9.419 1.14e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2074 on 106 degrees of freedom
## Multiple R-squared:  0.8576, Adjusted R-squared:  0.8495 
## F-statistic: 106.4 on 6 and 106 DF,  p-value: < 2.2e-16

模型結果可以看出R12與目標變數呈現明顯線性關係，這也是我們一開始自相關時知道的。
month呈現顯著，推測原因同上，month及R12此二變數的相依性高。
調整後R平方達到0.85，可以接受。
誤差項(Intercept)P值不顯著。
R1及R6的P值都很大，可以考慮剃除。

#試著剃除R1及R6在一次擬合
lm.fit <- lm(target_ind~.,data=train[,!(colnames(train) %in% c('R1','R6'))])
summary(lm.fit)

## 
## Call:
## lm(formula = target_ind ~ ., data = train[, !(colnames(train) %in% 
##     c("R1", "R6"))])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5071.5 -1409.8  -208.1  1497.1  6751.9 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 8.990e+02  1.360e+03   0.661  0.50991    
## month       4.441e+02  9.906e+01   4.483 1.84e-05 ***
## R4          7.455e-02  3.953e-02   1.886  0.06204 .  
## R8          1.371e-01  4.822e-02   2.842  0.00536 ** 
## R12         6.606e-01  6.922e-02   9.545 5.07e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2059 on 108 degrees of freedom
## Multiple R-squared:  0.8569, Adjusted R-squared:  0.8516 
## F-statistic: 161.7 on 4 and 108 DF,  p-value: < 2.2e-16

剔除後差異不大。 ***

多項式迴歸

接著我們嘗試使用非線性的方法改善模型，這裡採用多項式方法。因為R1和R6明顯不顯著，R12明顯顯著，這裡使用R4、R8兩變數做多項式，嘗試最高為5次。

lm.fit <- lm(target_ind~ month + R1 + R4 + I(R4^2) + I(R4^3) + I(R4^4) + I(R4^5) 
           + R6 + R8 + I(R8^2) + I(R8^3) + I(R8^4) + I(R8^5) + R12 ,data=train)
summary(lm.fit)

## 
## Call:
## lm(formula = target_ind ~ month + R1 + R4 + I(R4^2) + I(R4^3) + 
##     I(R4^4) + I(R4^5) + R6 + R8 + I(R8^2) + I(R8^3) + I(R8^4) + 
##     I(R8^5) + R12, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4408.1 -1240.3  -205.8  1498.0  5071.6 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.643e+05  3.639e+05   0.726    0.469    
## month        5.232e+02  1.146e+02   4.566 1.44e-05 ***
## R1           2.623e-02  4.058e-02   0.646    0.520    
## R4           5.429e+00  4.550e+01   0.119    0.905    
## I(R4^2)     -3.601e-04  3.583e-03  -0.100    0.920    
## I(R4^3)      1.122e-08  1.377e-07   0.081    0.935    
## I(R4^4)     -1.685e-13  2.588e-12  -0.065    0.948    
## I(R4^5)      1.035e-18  1.902e-17   0.054    0.957    
## R6          -9.190e-02  4.691e-02  -1.959    0.053 .  
## R8          -6.113e+01  5.192e+01  -1.177    0.242    
## I(R8^2)      4.969e-03  4.019e-03   1.236    0.219    
## I(R8^3)     -1.942e-07  1.522e-07  -1.276    0.205    
## I(R8^4)      3.679e-12  2.821e-12   1.304    0.195    
## I(R8^5)     -2.715e-17  2.050e-17  -1.324    0.188    
## R12          5.895e-01  7.686e-02   7.670 1.28e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2022 on 98 degrees of freedom
## Multiple R-squared:  0.8747, Adjusted R-squared:  0.8569 
## F-statistic: 48.89 on 14 and 98 DF,  p-value: < 2.2e-16

逐步回歸

由於涉及到的變數較多且彼此有相依性，我們用逐步回歸刪除影響較小的變數。

lm.fit <- step(lm.fit)

## Start:  AIC=1734.22
## target_ind ~ month + R1 + R4 + I(R4^2) + I(R4^3) + I(R4^4) + 
##     I(R4^5) + R6 + R8 + I(R8^2) + I(R8^3) + I(R8^4) + I(R8^5) + 
##     R12
## 
##           Df Sum of Sq       RSS    AIC
## - I(R4^5)  1     12101 400825065 1732.2
## - I(R4^4)  1     17340 400830305 1732.2
## - I(R4^3)  1     27148 400840113 1732.2
## - I(R4^2)  1     41308 400854272 1732.2
## - R4       1     58226 400871191 1732.2
## - R1       1   1708723 402521687 1732.7
## - R8       1   5669346 406482311 1733.8
## - I(R8^2)  1   6250455 407063419 1734.0
## - I(R8^3)  1   6659230 407472194 1734.1
## - I(R8^4)  1   6954305 407767269 1734.2
## <none>                 400812965 1734.2
## - I(R8^5)  1   7173149 407986114 1734.2
## - R6       1  15695998 416508963 1736.6
## - month    1  85278959 486091924 1754.0
## - R12      1 240586985 641399949 1785.3
## 
## Step:  AIC=1732.23
## target_ind ~ month + R1 + R4 + I(R4^2) + I(R4^3) + I(R4^4) + 
##     R6 + R8 + I(R8^2) + I(R8^3) + I(R8^4) + I(R8^5) + R12
## 
##           Df Sum of Sq       RSS    AIC
## - I(R4^4)  1    248012 401073077 1730.3
## - I(R4^3)  1    376496 401201561 1730.3
## - I(R4^2)  1    463875 401288941 1730.4
## - R4       1    498457 401323522 1730.4
## - R1       1   1786333 402611399 1730.7
## - R8       1   5766624 406591689 1731.8
## - I(R8^2)  1   6332455 407157521 1732.0
## - I(R8^3)  1   6722715 407547781 1732.1
## - I(R8^4)  1   6998947 407824013 1732.2
## <none>                 400825065 1732.2
## - I(R8^5)  1   7200642 408025708 1732.2
## - R6       1  15697552 416522617 1734.6
## - month    1  85566496 486391562 1752.1
## - R12      1 242716686 643541751 1783.7
## 
## Step:  AIC=1730.3
## target_ind ~ month + R1 + R4 + I(R4^2) + I(R4^3) + R6 + R8 + 
##     I(R8^2) + I(R8^3) + I(R8^4) + I(R8^5) + R12
## 
##           Df Sum of Sq       RSS    AIC
## - R4       1   1309323 402382400 1728.7
## - R1       1   1673950 402747027 1728.8
## - I(R4^2)  1   2255530 403328607 1728.9
## - I(R4^3)  1   3618906 404691983 1729.3
## - R8       1   5874109 406947186 1729.9
## - I(R8^2)  1   6432792 407505869 1730.1
## - I(R8^3)  1   6818164 407891240 1730.2
## - I(R8^4)  1   7092670 408165747 1730.3
## <none>                 401073077 1730.3
## - I(R8^5)  1   7295763 408368840 1730.3
## - R6       1  15458105 416531182 1732.6
## - month    1  85847576 486920653 1750.2
## - R12      1 243622053 644695129 1781.9
## 
## Step:  AIC=1728.66
## target_ind ~ month + R1 + I(R4^2) + I(R4^3) + R6 + R8 + I(R8^2) + 
##     I(R8^3) + I(R8^4) + I(R8^5) + R12
## 
##           Df Sum of Sq       RSS    AIC
## - R1       1   1843309 404225709 1727.2
## - R8       1   5426971 407809371 1728.2
## - I(R8^2)  1   5957139 408339539 1728.3
## - I(R8^3)  1   6328208 408710608 1728.4
## - I(R8^4)  1   6598847 408981247 1728.5
## - I(R8^5)  1   6805494 409187893 1728.6
## <none>                 402382400 1728.7
## - I(R4^2)  1  13555515 415937915 1730.4
## - R6       1  15189855 417572254 1730.8
## - I(R4^3)  1  17549001 419931400 1731.5
## - month    1  84792914 487175313 1748.3
## - R12      1 247577347 649959747 1780.8
## 
## Step:  AIC=1727.18
## target_ind ~ month + I(R4^2) + I(R4^3) + R6 + R8 + I(R8^2) + 
##     I(R8^3) + I(R8^4) + I(R8^5) + R12
## 
##           Df Sum of Sq       RSS    AIC
## - R8       1   5669333 409895041 1726.8
## - I(R8^2)  1   6159443 410385152 1726.9
## - I(R8^3)  1   6488205 410713913 1727.0
## - I(R8^4)  1   6715544 410941252 1727.0
## - I(R8^5)  1   6877913 411103621 1727.1
## <none>                 404225709 1727.2
## - I(R4^2)  1  12687002 416912711 1728.7
## - R6       1  13843149 418068858 1729.0
## - I(R4^3)  1  16600997 420826706 1729.7
## - month    1  86993437 491219146 1747.2
## - R12      1 251079708 655305416 1779.8
## 
## Step:  AIC=1726.75
## target_ind ~ month + I(R4^2) + I(R4^3) + R6 + I(R8^2) + I(R8^3) + 
##     I(R8^4) + I(R8^5) + R12
## 
##           Df Sum of Sq       RSS    AIC
## - I(R8^5)  1   2320810 412215851 1725.4
## - I(R8^4)  1   2679576 412574617 1725.5
## - I(R8^3)  1   3245144 413140186 1725.7
## - I(R8^2)  1   4167790 414062832 1725.9
## <none>                 409895041 1726.8
## - R6       1  12618388 422513429 1728.2
## - I(R4^2)  1  16344932 426239973 1729.2
## - I(R4^3)  1  20910198 430805240 1730.4
## - month    1  81705493 491600535 1745.3
## - R12      1 283965502 693860543 1784.2
## 
## Step:  AIC=1725.39
## target_ind ~ month + I(R4^2) + I(R4^3) + R6 + I(R8^2) + I(R8^3) + 
##     I(R8^4) + R12
## 
##           Df Sum of Sq       RSS    AIC
## - I(R8^4)  1   3559952 415775803 1724.4
## - I(R8^3)  1   5294134 417509985 1724.8
## <none>                 412215851 1725.4
## - I(R8^2)  1   7979387 420195238 1725.6
## - R6       1  13097591 425313442 1726.9
## - I(R4^2)  1  17607510 429823361 1728.1
## - I(R4^3)  1  22512166 434728017 1729.4
## - month    1  80326290 492542141 1743.5
## - R12      1 290296291 702512142 1783.6
## 
## Step:  AIC=1724.36
## target_ind ~ month + I(R4^2) + I(R4^3) + R6 + I(R8^2) + I(R8^3) + 
##     R12
## 
##           Df Sum of Sq       RSS    AIC
## <none>                 415775803 1724.4
## - R6       1  10308992 426084795 1725.1
## - I(R4^2)  1  16039899 431815702 1726.6
## - I(R4^3)  1  20793153 436568957 1727.9
## - I(R8^3)  1  21664667 437440470 1728.1
## - I(R8^2)  1  26622098 442397901 1729.4
## - month    1  79411931 495187734 1742.1
## - R12      1 296353662 712129465 1783.2

summary(lm.fit)

## 
## Call:
## lm(formula = target_ind ~ month + I(R4^2) + I(R4^3) + R6 + I(R8^2) + 
##     I(R8^3) + R12, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4984.6 -1318.5   -72.3  1510.7  5191.8 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.302e+03  2.531e+03   2.095   0.0386 *  
## month        4.907e+02  1.096e+02   4.478 1.92e-05 ***
## I(R4^2)     -1.169e-05  5.806e-06  -2.013   0.0467 *  
## I(R4^3)      3.077e-10  1.343e-10   2.292   0.0239 *  
## R6          -6.919e-02  4.288e-02  -1.614   0.1096    
## I(R8^2)      1.887e-05  7.276e-06   2.593   0.0109 *  
## I(R8^3)     -3.923e-10  1.677e-10  -2.339   0.0212 *  
## R12          6.271e-01  7.249e-02   8.651 6.40e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1990 on 105 degrees of freedom
## Multiple R-squared:  0.8701, Adjusted R-squared:  0.8614 
## F-statistic: 100.4 on 7 and 105 DF,  p-value: < 2.2e-16

我們估且接受這樣的結果。 lm.fit就是我們建立用於時間序列預測的線性迴歸模型。

test$pred <- predict(lm.fit,test)
#計算百分誤差率
test$diff <- abs(test$target_ind-test$pred)/test$target_ind

觀察預測結果

kable(test) %>% kable_styling() %>% scroll_box(height = "200px")

	target_ind	month	R1	R4	R6	R8	R12	pred	diff
115	27595	7	22316	25453	14672	31222	29009	29866.35	0.0823102
116	25451	8	27595	32683	17543	38641	25500	27798.10	0.0922204
117	25421	9	25451	22449	25453	14672	24166	23526.15	0.0745390
118	25288	10	25421	22316	32683	17543	26960	26142.57	0.0337935
119	32568	11	25288	27595	22449	25453	31222	32046.22	0.0160212
120	35110	12	32568	25451	22316	32683	38641	37836.65	0.0776601
121	16052	1	35110	25421	27595	22449	14672	15656.60	0.0246323
122	22146	2	16052	25288	25451	22316	17543	18062.03	0.1844114
123	21198	3	22146	32568	25421	27595	25453	25333.64	0.1950956
125	22084	5	19543	16052	32568	25421	22449	23587.93	0.0681006
126	23816	6	22084	22146	35110	25288	22316	23142.11	0.0282955
127	29961	7	23816	21198	16052	32568	27595	29068.61	0.0297850
128	26773	8	29961	19543	22146	35110	25451	27764.79	0.0370443
129	26635	9	26773	22084	21198	16052	25421	25046.73	0.0596311
130	26972	10	26635	23816	19543	22146	25288	27235.01	0.0097512
131	30207	11	26972	29961	22084	21198	32568	32121.74	0.0633874
132	38687	12	30207	26773	23816	19543	35110	33366.23	0.1375338
133	16974	1	38687	26635	29961	22084	16052	16285.53	0.0405605
134	21697	2	16974	26972	26773	23816	22146	21256.25	0.0203140
135	24179	3	21697	30207	26635	29961	21198	22426.49	0.0724808
136	23757	4	24179	38687	26972	26773	19543	23975.80	0.0092099
137	25013	5	23757	16974	30207	26635	22084	23622.76	0.0555806
138	24019	6	25013	21697	38687	26972	23816	24172.58	0.0063942
139	30345	7	24019	24179	16974	30207	29961	30269.35	0.0024928
140	24488	8	30345	23757	21697	38687	26773	27565.18	0.1256609
141	25156	9	24488	25013	24179	16974	26635	25769.57	0.0243906
142	25650	10	25156	24019	23757	21697	26972	27875.78	0.0867750
143	30923	11	25650	30345	25013	24179	30207	31233.41	0.0100383
144	37240	12	30923	24488	24019	23757	38687	36687.52	0.0148357
145	17466	1	37240	25156	30345	25013	16974	17504.66	0.0022133
146	19463	2	17466	25650	24488	24019	21697	21147.18	0.0865324
147	24352	3	19463	30923	25156	30345	24179	24530.17	0.0073165
148	26805	4	24352	37240	25650	24488	23757	25626.36	0.0439710
149	25236	5	26805	17466	30923	25156	25013	25069.21	0.0066092
150	24735	6	25236	19463	37240	25650	24019	24365.18	0.0149515
151	29356	7	24735	24352	17466	30923	30345	30510.43	0.0393253
152	31234	8	29356	26805	19463	37240	24488	26668.72	0.1461639
153	22724	9	31234	25236	24352	17466	25156	24977.09	0.0991502
154	28496	10	22724	24735	26805	19463	25650	26200.82	0.0805439
155	32857	11	28496	29356	25236	24352	30923	31581.82	0.0388100
156	37198	12	32857	31234	24735	26805	37240	36807.50	0.0104978
157	13652	1	37198	22724	29356	25236	17466	18000.20	0.3185030
158	22784	2	13652	28496	31234	24735	19463	19563.98	0.1413282
159	23565	3	22784	32857	22724	29356	24352	25105.21	0.0653603
160	26323	4	23565	37198	28496	31234	26805	28221.70	0.0721308
161	23779	5	26323	13652	32857	22724	25236	25050.73	0.0534812
162	27549	6	23779	22784	37198	28496	24735	24998.05	0.0925970
163	29660	7	27549	23565	13652	32857	29356	30189.14	0.0178402
164	23356	8	29660	26323	22784	37198	31234	30663.99	0.3128958

summary(test)

##    target_ind        month              R1              R4       
##  Min.   :13652   Min.   : 1.000   Min.   :13652   Min.   :13652  
##  1st Qu.:23565   1st Qu.: 4.000   1st Qu.:23565   1st Qu.:22724  
##  Median :25288   Median : 7.000   Median :25288   Median :25236  
##  Mean   :26041   Mean   : 6.592   Mean   :25986   Mean   :25694  
##  3rd Qu.:29356   3rd Qu.: 9.000   3rd Qu.:29356   3rd Qu.:28496  
##  Max.   :38687   Max.   :12.000   Max.   :38687   Max.   :38687  
##        R6              R8             R12             pred      
##  Min.   :13652   Min.   :14672   Min.   :14672   Min.   :15657  
##  1st Qu.:22316   1st Qu.:22449   1st Qu.:22449   1st Qu.:23623  
##  Median :25156   Median :25288   Median :25288   Median :25626  
##  Mean   :25533   Mean   :26212   Mean   :25745   Mean   :26235  
##  3rd Qu.:29356   3rd Qu.:30207   3rd Qu.:29009   3rd Qu.:29866  
##  Max.   :38687   Max.   :38687   Max.   :38687   Max.   :37837  
##       diff         
##  Min.   :0.002213  
##  1st Qu.:0.017840  
##  Median :0.053481  
##  Mean   :0.068024  
##  3rd Qu.:0.086532  
##  Max.   :0.318503

從結果可以看到最大百分誤差率為31.9%，平均誤差為7%，整體預測結果不錯。

推薦閱讀：
Edward - 迴歸模型介紹

本篇文章按照游皓麟老師的實戰R語言預測分析介紹進行實作，含金量相當高的好書。

時序思考題：內建資料集Airmiles

各位有興趣可以再用R內建的資料集 airmiles 思考看看，跟本文的範例有哪些特性相同或相異之處？
這個資料集是美國某航空公司的全年飛機里程資料。

airmiles;summary(airmiles)

## Time Series:
## Start = 1937 
## End = 1960 
## Frequency = 1 
##  [1]   412   480   683  1052  1385  1418  1634  2178  3362  5948  6109  5981
## [13]  6753  8003 10566 12528 14760 16769 19819 22362 25340 25343 29269 30514

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     412    1580    6431   10528   17532   30514

plot(airmiles) #先看一下折線圖

acf(airmiles) #自相關

pacf(airmiles) #偏相關