1 課程目標

有些資料可以按照時間先後排列，我們想描述資料點之間的關聯，判斷趨勢的型態，也想預測未來的資料點，或者想找出其他的資料跟時間數列資料的關係，這都是時間數列分析的目標。

2 數列分析的基本原理

  • 依照資料發生時間排序的資料，稱為時間數列。
  • 時間數列中的觀察值大小，決定於長期趨勢、循環變動、季節變動、不規則變動。觀察值之間有關連，時間隔得越久，關聯將會越小。
  • 長期趨勢指的是長時間的趨勢，忽視短時間或者循環的變動。比方說，長期以來識字率一直上升，即使有某幾年可能下降或者沒有變動。
  • 循環變動指的是相對長期的資料變動，可能需要好幾年才會出現一次循環，例如經濟景氣循環。
  • 季節變動指的是每一年都會出現的變動，例如溫度隨著季節升降，每一年都有類似的循環。
  • 不規則變動也稱為雜訊(noice)，可能因為一些突發事件或者沒有控制到的變數而出現的變化。
  • 我們用\(T_{t}\)代表長期趨勢，\(C_{t}\)代表循環變動， \(S_{t}\)代表季節變動，\(N_{t}\)代表不規則變動，資料由4種因素組成。

    • \[\begin{equation} Y_{t}=T_{t}\times C_{t}\times S_{t}\times I_{t} \tag{2.1} \end{equation}\]

    方程式(2.1)顯示，時間序列的相乘模型。我們也可以假設資料呈現相加模型如方程式(2.2)：

    \[\begin{equation} Y_{t}=T_{t}+C_{t}+S_{t}+I_{t} \tag{2.2} \end{equation}\]

  • 相加模型假設長期趨勢、季節變動等四個變數之間不相關、互相獨立，也就是季節的變動是固定的，例如啤酒的銷售在夏季最好，每一年都是如此，例如圖 2.1。如果隨著時間，季節變動會增加，例如每個月搭飛機的乘客人數，有可能因為每年的各種狀況，出現每個月的乘客人數上升或者下降，例如圖 2.2。相乘模型比較適用這種資料。 
    • library(fpp)
plot(as.ts(fpp::ausbeer))
    澳洲啤酒的每一季銷售量

    Figure 2.1: 澳洲啤酒的每一季銷售量 

    plot(as.ts(fpp::ausair))
    澳洲註冊的航空公司的每年載客人數

    Figure 2.2: 澳洲註冊的航空公司的每年載客人數

    3 分析步驟

  • 時間數列分析的目標之一是分離長期趨勢以及短期變動的作用，也就是拆開\(T\cdot C\)以及 \(S\cdot I\)。從方程式(2.1)可得出方程式(3.1)：

    • \[\begin{equation} S_{t}\times I_{t}=\frac{Y_{t}}{ T_{t}\times C_{t}} \tag{3.1} \end{equation}\]

  • 按照以下步驟分析時間數列：
    1. 計算\(K\)期的移動總和(moving totals)
    2. 計算\(K\)期的移動平均(moving average)
    3. 計算中央移動平均數
    4. 根據方程式(3.1)，計算\(S_{t}\times I_{t}\)

    3.1 手動計算

  • 我們使用教科書的資料如表 3.1： 
    • d0<-matrix(c(288.84,293.74,287.43,277.29,
               311.88,309.4,283.59,282.4,
               292.09,293.53,281.95,275.26,
          278.45,279.57,269.94, 259.9),ncol=4,nrow=4, byrow = T)
colnames(d0)<-c('Spring','Summer','Autumn','Winter')
rownames(d0)<-c('104','105','106','107')
kbl(d0, caption="正新橡膠營業收入") %>%
kable_styling('striped',font_size = 20)
    Table 3.1: 正新橡膠營業收入
    Spring Summer Autumn Winter
    104 288.8 293.7 287.4 277.3
    105 311.9 309.4 283.6 282.4
    106 292.1 293.5 281.9 275.3
    107 278.4 279.6 269.9 259.9
  • 為了分析方便，我們修改一下資料格式如表 3.2： 
    • df<-data.frame(spring=c(288.84,293.74,287.43,277.29),
               summer=c(311.88,309.4,283.59,282.47),
               autumn=c(292.09,293.53,281.95,275.26),
               winter=c(278.45,279.57,269.94, 259.9))
kbl(df, caption="正新橡膠資料") %>%
kable_styling('striped',font_size = 20)
    Table 3.2: 正新橡膠資料
    spring summer autumn winter
    288.8 311.9 292.1 278.4
    293.7 309.4 293.5 279.6
    287.4 283.6 281.9 269.9
    277.3 282.5 275.3 259.9
    1. 轉換為向量資料，再用ts指令轉換為時間數列如表 3.3。請注意轉換的過程。 
    df<-as.matrix(df)
t.df<-t(df)
dat<-as.vector(t.df)
dat_ts<-ts(dat, frequency = 4, start=c(104,1))
kbl(dat_ts, align='l', caption = "轉換成向量的正新橡膠資料") %>%
kable_styling('striped',font_size = 20, position='left')
    Table 3.3: 轉換成向量的正新橡膠資料
    x
    288.8
    311.9
    292.1
    278.4
    293.7
    309.4
    293.5
    279.6
    287.4
    283.6
    281.9
    269.9
    277.3
    282.5
    275.3
    259.9
  • 計算四季的移動平均，第一個四期平均為104年的1到4季的平均，第二個則是去掉104年的第1季之後，104年的2到3季加上105年的第1季計算得到的平均。依此類推，得到的平均數如表 3.4： 
    • ma <- function(x, n = 4){stats::filter(x, rep(1 / n, n), sides = 2)}
dat.ma<-ma(dat)
kbl(dat.ma, align='l', caption = "四季移動平均資料") %>%
kable_styling('striped',font_size = 20, position='left')
    Table 3.4: 四季移動平均資料
    x
    NA
    292.8
    294.0
    293.4
    293.8
    294.1
    292.5
    286.0
    283.1
    280.7
    278.2
    277.9
    276.2
    273.7
    NA
    NA
  • 根據表 3.4的資料，可以畫四季的移動平均數列如圖 3.1 
    • plot.ts(dat.ma)
    四季移動平均線

    Figure 3.1: 四季移動平均線

  • 第二個步驟是根據表 3.4的資料，計算中央移動平均。第一個資料值為前兩個四季平均數的平均數，接下來去掉第一個四季平均數，往下加一個四季平均數，計算出第二個資料值，依此類推，得到結果如表 3.5： 
    • macenter <- function(x, n = 2){stats::filter(x, rep(1 / n, n), sides = 2)}
dat.center<-macenter(dat.ma)
kbl(dat.center, align='l', caption = "平均移動資料")%>%
kable_styling('striped',font_size = 20, position='left')
    Table 3.5: 平均移動資料
    x
    NA
    293.4
    293.7
    293.6
    293.9
    293.3
    289.3
    284.6
    281.9
    279.5
    278.1
    277.1
    275.0
    NA
    NA
    NA
  • 畫中央移動平均趨勢如圖 3.2： 
    • plot.ts(dat.center)
    中央移動平均線

    Figure 3.2: 中央移動平均線

    1. 第三個步驟是根據表 3.5的資料，計算季節不規則成分。第1個資料值為第3季原始資料/第1個中央移動平均值，第2個資料值為第4季原始資料/第2個中央移動平均值，依此類推，得到結果如表 3.6： 
    options(digits = 3)
#原始資料
df<-data.frame(spring=c(288.84,293.74,287.43,277.29),
               summer=c(311.88,309.4,283.59,282.47),
               autumn=c(292.09,293.53,281.95,275.26),
               winter=c(278.45,279.57,269.94, 259.9))
df<-as.matrix(df)
t.df<-t(df)
revenue<-as.vector(t.df)

#中央平均
ma <- function(x, n = 4){stats::filter(x, rep(1 / n, n), sides = 2)}
dat.ma<-ma(revenue)
macenter <- function(x, n = 2){stats::filter(x, rep(1 / n, n), sides = 2)}
dat.center<-macenter(dat.ma)
MA=as.numeric(format(dat.center[2:13], digits=2))

#季節不規則
season=paste(c(rep("104",2),rep("105",4),rep("106",4),rep("107",2)), 
             c(3,4,rep(c(1,2,3,4),2), 1,2),
             sep=".")
dat.ir <- data.table::data.table(season,revenue=revenue[3:14], MA, 
            irregular=as.numeric(format(revenue[3:14]/MA, digits=2)))

kbl(dat.ir, align='l', caption = "季節不規則成分")%>%
kable_styling('striped',font_size = 20, position='left')
    Table 3.6: 季節不規則成分
    season revenue MA irregular
    104.3 292 293 1.00
    104.4 278 294 0.95
    105.1 294 294 1.00
    105.2 309 294 1.05
    105.3 294 293 1.00
    105.4 280 289 0.97
    106.1 287 285 1.01
    106.2 284 282 1.01
    106.3 282 279 1.01
    106.4 270 278 0.97
    107.1 277 277 1.00
    107.2 282 275 1.03
    1. 第四個步驟是根據表 3.6的季節不規則成分計算季節因子。計算方式為加總每一年的特定季節的不規則成分，然後求平均數。我們先整理表 3.6得到表 3.7： 
    #季節不規則
year=c(rep("104",2),rep("105",4),rep("106",4),rep("107",2))
season=c("autumn","winter",rep(c("spring","summer","autumn","winter"),2), "spring","summer")
dat.irt=data.frame(year,season, x=dat.ir$irregular)
kbl(dat.irt, align='l', caption = "個別年度季節不規則成分") %>%
kable_styling('striped',font_size = 20, position='left')
    Table 3.7: 個別年度季節不規則成分
    year season x
    104 autumn 1.00
    104 winter 0.95
    105 spring 1.00
    105 summer 1.05
    105 autumn 1.00
    105 winter 0.97
    106 spring 1.01
    106 summer 1.01
    106 autumn 1.01
    106 winter 0.97
    107 spring 1.00
    107 summer 1.03
    1. 根據表 3.7，加總每一年的特定季節的不規則成分，然後求平均數。得到結果如表 3.8： 
    library(dplyr)
dat.irt$season<-factor(dat.irt$season, levels =  c("spring","summer","autumn","winter"))
ir.table<-dat.irt %>% group_by(season)%>%
          summarize(n=n(), sum=sum(x), avg=sum/n)

kbl(ir.table, align='l', caption = "季節因子")%>%
kable_styling('striped',font_size = 20, position='left')
    Table 3.8: 季節因子
    season n sum avg
    spring 3 3.01 1.003
    summer 3 3.09 1.030
    autumn 3 3.01 1.003
    winter 3 2.89 0.963
    1. 根據表 3.8，可計算季節指數如 \[\begin{equation} \textrm{season index}=\frac{\textrm{season factor}_{i}}{\sum \textrm{season factor}}\times 4 \end{equation}\] 得到結果如表 3.9，夏季的營業收入會比較高，冬季比較低。 
    options(digits = 4)
S=sum(ir.table$avg)
seasonix<-4*ir.table$avg/S
Table<-rbind(ir.table$avg, seasonix)
colnames(Table)<-c("spring","summer","autumn","winter")
rownames(Table)<-c("季節因子","季節指數")
kbl(Table, align='l', caption = "季節因子與指數")%>%
kable_styling('striped',font_size = 20, position='left')
    Table 3.9: 季節因子與指數
    spring summer autumn winter
    季節因子 1.003 1.03 1.003 0.9633
    季節指數 1.003 1.03 1.003 0.9633
    1. 根據表 3.9，可計算消除季節性因素的營業收入\(Y^{'}_{t}\)，按照方程式(3.2)得到結果如表 3.10：

    \[\begin{equation} Y^{'}_{t}=Y_{t}/S_{t}=T_{t}\cdot C_{t}\cdot I_{t} \tag{3.2} \end{equation}\]

    df<-data.frame(spring=c(288.84,293.74,287.43,277.29),
               summer=c(311.88,309.4,283.59,282.47),
               autumn=c(292.09,293.53,281.95,275.26),
               winter=c(278.45,279.57,269.94, 259.9))
df<-as.matrix(df)
t.df<-t(df)
revenue<-as.vector(t.df)
yprime<-revenue/seasonix
DF<-data.frame(revenue, seasonix, yprime)
kbl(DF, align='l', caption = "消除季節因素的營業收入")%>%
kable_styling('striped',font_size = 20, position='left')
    Table 3.10: 消除季節因素的營業收入
    revenue seasonix yprime
    288.8 1.0033 287.9
    311.9 1.0300 302.8
    292.1 1.0033 291.1
    278.4 0.9633 289.0
    293.7 1.0033 292.8
    309.4 1.0300 300.4
    293.5 1.0033 292.6
    279.6 0.9633 290.2
    287.4 1.0033 286.5
    283.6 1.0300 275.3
    281.9 1.0033 281.0
    269.9 0.9633 280.2
    277.3 1.0033 276.4
    282.5 1.0300 274.2
    275.3 1.0033 274.3
    259.9 0.9633 269.8
    1. 根據表 3.10，可計算消除季節性因素後的營業收入長期趨勢，也就是用季節(\(t\))當做自變項解釋\(Y^{'}_{t}\)：

    \[\begin{equation} Y^{'}_{t}=Y_{t}/S_{t}=\alpha+\beta\cdot t+\epsilon \tag{3.3} \end{equation}\]

    t<-c(1:16)
DFm<-data.frame(yprime, t)
m1 <- lm(yprime ~ t, data=DFm)
DFols<-summary(m1)
stargazer::stargazer(m1, type = "html", 
          title = "消除季節因素的營業收入的迴歸模型", 
          header = FALSE, style='apsr',
          single.row = TRUE)
    (#tab:)消除季節因素的營業收入的迴歸模型
    yprime
    t -1.737*** (0.291)
    Constant 300.000*** (2.812)
    N 16
    R2 0.718
    Adjusted R2 0.698
    Residual Std. Error 5.363 (df = 14)
    F Statistic 35.660*** (df = 1; 14)
    p < .1; p < .05; p < .01
  • 根據表 ，該公司的營業收入\(Y^{'}_{t}=300-1.737\cdot t\)。-1.737表示在過去16個月，每個月少了1.737億元。
  • 根據\(Y^{'}_{t}=300-1.737\cdot t\)，代入\(t\)可得到每一季的趨勢值。趨勢值乘以季節指數，即為預測值。以108年為例，\(t\)代入17,18,19,20，計算108年四季的預測值如下： 
    • t_108<-c(17:20)
yprimenew=300-1.737*t_108
pred=seasonix[1:4]*yprimenew
season=c("spring","summer","autumn","winter")
predt<-data.frame(season, y_prime=yprime[13:16], season_t=seasonix[1:4], predict=pred)
 kbl(predt, align='l', caption = "108年的營業收入的預測值") %>%
 kable_styling('striped',font_size = 20, position='left' )
    Table 3.11: 108年的營業收入的預測值
    season y_prime season_t predict
    spring 276.4 1.0033 271.4
    summer 274.2 1.0300 276.8
    autumn 274.3 1.0033 267.9
    winter 269.8 0.9633 255.5

    4 迴歸分析方法

  • 以上的分析基本上假設每一段時期（年、季、月、週）對於依變數\(y\)的期望值有固定的作用()。我們稱之為deterministic （確定性）方法。例如每多一歲，BMI上升0.1；年資每增加一年，薪水增加3%。
  • 但是\(y\)有可能受到\(y_{t-p}\)的影響或者\(e_{t-q}\)的影響，其中，\(\{t,q\}\in1,2,\ldots\)這個方法稱為dynamic stochastic processes(Lags)，是最有彈性的方法。
  • 第三個方法稱為seasonality (cycles)，\(y\)是固定循環\(c\)的加法或者乘法，也就是\(y=f(y_{t-c})\)。

    • 4.1 確定性模型

  • 確定性模型指的是時間對依變數有固定的作用，可以寫成：

    • \[y_{t}=\beta_{0}+t\beta_{1}+\epsilon_{t},\quad \epsilon_{t}\sim N(0,\sigma^2)\]

    • 例如我們建立一個24個時間點的資料，：
    set.seed(02138)
t=c(1:24)
Y=0.5+t*2.4
plot(t, Y, col='#FF11AA', pch=20, ylim=c(0,40))
abline(lm(Y~t), lty=1, col='#11CC11')
axis(1, c(1:24))
    24個時間點的趨勢圖1

    Figure 4.1: 24個時間點的趨勢圖1

    • 加上誤差項，也就是雜訊，出現上下的波動如圖4.2：
    set.seed(02138)
t=c(1:24)
Y=0.5 + t*2.4 + rnorm(24, 2, 5) #time series
Yseries<-ts(Y) #time series
plot.ts(Yseries, col='#DD1111', pch=20, 
              ylim=c(0,60), lwd=1.5)
Y1=0.5 + t*2.4 #no noisy
points(t, Y1,  col='#201BCC', lwd=1.5, type='l')

axis(1, c(1:24))
legend(1, 62, lty=c(1,1,1),c('time series',
              'trend'), bty='n', col=c('#DD1111','#201BCC'))
    24個時間點的趨勢圖2

    Figure 4.2: 24個時間點的趨勢圖2

    • 再加上線性迴歸，線性迴歸的線跟趨勢線很接近，但是還是有所區別：
    set.seed(02138)
t=c(1: 24)
Y=0.5 + t*2.4 + rnorm(24, 2, 5) #time series
Yseries<-ts(Y) #time series
plot.ts(Yseries, col='#DD1111', pch=20, 
              ylim=c(0,60), lwd=1.5)
Y1=0.5 + t*2.4 #no noisy
points(t, Y1,  col='#201BCC', lwd=1.5, type='l')
abline(lm(Y ~ t), lty=2, col='#312C00', lwd=1.5)
axis(1, c(1:24))
legend(1, 62, lty=c(1,1,2),c('time series',
              'trend','linear estimate'), bty='n', col=c('#DD1111','#201BCC', '#312C00'))
    24個時間點的趨勢圖3

    Figure 4.3: 24個時間點的趨勢圖3

    • 在時間數列的分析中，不同起始點與終點可能得到不同的\(\hat{\beta}\)，所以挑選起始點跟終點很重要。圖4.4顯示，假如我們
    set.seed(02138)
t=c(1:24)
Y=0.5 + t*2.4 + rnorm(24, 2, 5) #time series
Yseries<-ts(Y) #time series
plot.ts(Yseries, col='#33EE11', pch=21, type='p',
              ylim=c(0,60), lwd=1.5)
points(t[c(4:10)], Yseries[c(4:10)], col='#DD1111', pch=16, type='p')

Y1=0.5 + t*2.4 #no noisy
points(t, Y1,  col='#201BCC', lwd=1.5, type='l')
abline(lm(Y[c(4:10)] ~ t[c(4:10)]), lty=2, col='#312C00', lwd=1.5)
axis(1, c(1:24))
legend(1, 62, lty=c(1,1,2),c('time series',
              'trend','linear estimate'), bty='n', col=c('#DD1111','#201BCC', '#312C00'))
    24個時間點的趨勢圖3

    Figure 4.4: 24個時間點的趨勢圖3

    • 在時間數列的分析中，不同起始點與終點可能得到不同的\(\hat{\beta}\)，所以挑選起始點跟終點很重要。時間數列分析要注意selection bias.
  • 當觀察的時期(period)越來越多，例如超過100，\(\hat{\beta}\)會非常接近\(\beta\)，離散程度也會變小。
  • 因此，假設時間對於依變數有固定的作用，前提是需要有許多連續的時間點，而且不能有selection bias。

    • 4.2 實例

  • 我們用一筆42位英國國王死亡的年齡的資料說明。首先讀入資料，然後用42位國王為X軸，Y軸為年齡，圖4.5顯示，死亡年齡似乎越來越高，表示國王越來越長壽。 
    • kings <- scan("http://robjhyndman.com/tsdldata/misc/kings.dat", skip=3)
t<-c(1:42)
plot(t, kings, type='l', col='#DD1111', lwd=1.5)
abline(lm(kings~t), col='#200BCC', lty=2, lwd=1.5)
    英國國王的死亡年齡趨勢圖1

    Figure 4.5: 英國國王的死亡年齡趨勢圖1

  • 我們用\(\tt{ts}\)指令，把資料變成時間數列，圖4.6顯示，死亡年齡有升有降，並不是像是線性模型所顯示的越來越高的趨勢，表示並不是每一個國王都越來越長壽。 
    • kings <- scan("http://robjhyndman.com/tsdldata/misc/kings.dat",skip=3)
kingsseries <-ts(kings)
plot.ts(kingsseries, type='l', col='#DD1111', lwd=1.5)
abline(lm(kings ~ t), col='#0AC200', add=T,
       lty=2, lwd=1.5)
    英國國王的死亡年齡趨勢圖2

    Figure 4.6: 英國國王的死亡年齡趨勢圖2

    4.3 Serial correlation

  • 如果把時間數列資料去掉時間的影響，而且誤差項是隨機分配，理論上應該是一直線，也就是\(y-t\hat{\beta}_{1}=\hat{\beta}_{0}+\hat{\epsilon}_{t}\)，如同圖4.7。 
    • set.seed(02138)
t=c(1:24)
e=rnorm(24, 2, 5)
Y=0.5 + t*2.4 + e #time series
Yseries<-ts(Y) #time series
plot.ts(Yseries, col='#22AACC', type='l',
              ylim=c(0,60), lwd=1.5, ylab=bquote(y[i]-t*beta))
points(t, Y - 2.4 * t - e, col='#DD1111', type='l')

axis(1, c(1:24))
    24個時間點的趨勢圖4

    Figure 4.7: 24個時間點的趨勢圖4

  • 如果我們去除時間的影響之後，還是發現時間數列隨著時間變動，那就是誤差之間可能有相關。也就是\(\epsilon_{t}\)跟\(\epsilon_{t-1}\)或者\(\epsilon_{t-2}\)相關。

    • \[\epsilon_{t}=\rho\epsilon_{t-1}+u_{t},\quad -1<\rho<1\]

    • \(\rho\)稱為first-order autocorrelation coefficient。這個時間數列稱為first-order autoregressive process(AR1)。當\(\rho>0\)，\(\epsilon\)固定為正或者負，當\(\rho<0\)，\(\epsilon\)有時候為正有時候是負。
  • 除了殘差互相影響外，\(y_{t-1}\)可能影響\(y_{t}\)。\(\epsilon_{t-1}\)也可能影響\(y_{t}\)。E\((y_{t-1})\)可能影響\(y_{t}\)。
  • \(y_{t-1}\)可能影響\(y_{t}\)稱為自我迴歸(autoregressive process)，可寫成：

    • \[y_{t}=y_{t-1}\phi_{1}+\epsilon_{t}\]

    • 因為\(y_{t-1}\)是\(y_{t-2}\)的函數，一直延續下去，所以：

    \[\begin{align*} y_{t}&=y_{t-1}\phi_{1} +\epsilon_{t}\\ &=(y_{t-2}\phi_{1}+\epsilon_{t-1})\phi_{1}+\epsilon_{t}\\ &= y_{t-2}\phi_{1}^2+\epsilon_{t-1}\phi_{1}+\epsilon_{t}\\ &=(y_{t-3}\phi_{1}+\epsilon_{t-2})\phi_{1}^2+\epsilon_{t-1}\phi_{1}+\epsilon_{t}\\ &\ldots \rm{substitue} \hspace{.2em}\rm{through}\hspace{.2em}y_{t-k}\\ &=y_{t-k}\phi_{1}^k+\sum\limits_{j=0}^{k-1}\epsilon_{t-j}\phi^{j}\\ &\ldots \rm{substitue} \hspace{.2em}\rm{through}\hspace{.2em}y_{t-\infty}\\ &=\sum\limits_{j=0}^{\infty}\epsilon_{t-j}\phi^{j} \end{align*}\]

  • 假如\(y_{t}\)是AR(1)，那麼\(y_{t}\)包括從時間數列一開始的誤差項的影響，也就是\(\epsilon_{t-1}\ldots\epsilon_{t-\infty}\)
    • 我們可以用\(\tt{arima.sim()}\)函式模擬自我相關的資料如圖 4.8：
    set.seed(02138)
y <- arima.sim(list(order=c(1,0,0),
    ar = 0.5,
    ma = NULL), n=200)
y=ts(y)
plot.ts(y, col='#CC22AA')
    ARIMA時間數列圖1

    Figure 4.8: ARIMA時間數列圖1

    • Christopher Adolph建立了\(\tt{makeTS()}\)函式。圖 4.9顯示\(\tt{makeTS()}\)模擬的兩筆時間數列資料的趨勢：
    adolph='http://faculty.washington.edu/cadolph/panEssex/makeTS.R'
source(adolph)
set.seed(02138)
Y <- makeTS(n=200,
             ar= 0, #phi=0
             ma= NULL,
            trend=0,
            seasons=c(0,0,0,0,0,0,0,0,0,0,0,0),
            adjust=c('level'),
              varY = 0, #Error Term variance
initY = 0,
initE = 0,
burnin = 10)
Y1 <- makeTS(n=200,
             ar= 0.5, #phi=0
             ma= NULL,
            trend=0,
            seasons=c(0,0,0,0,0,0,0,0,0,0,0,0),
            adjust=c('level'),
              varY = 1, #Error Term variance
initY = 0,
initE = 0,
burnin = 10)
plot.ts(Y1, col='#CC22AA')
t<-c(1:200)
lines(t, Y, type='l', col='#123CDA', lty=2)
    ARIMA時間數列圖2

    Figure 4.9: ARIMA時間數列圖2

    4.4 ACF, PACF

  • 上述的\(y_{t}\)與\(y_{t-1}\)之間的自我相關稱為Autocorrelation Function(ACF)，也就是自我相關函數，ACF可以寫成：

    • \[\begin{align*} \rm{ACF}_{j}&=\rm{cor}(y_{t},y_{t-j})=\frac{cov(y_{t},y_{t-j})}{\sqrt{var(y)}\sqrt{var(y)}} \\ &=\frac{cov(y_{t},y_{t-j})}{var(y)} \end{align*}\]

    • 用R可以畫圖：
    acf(Y1)

  • 如果我們可以控制從\(y_{t-1}\)到\(y_{t-j+1}\)，求出\(y_{t}\)以及\(y_{t-j}\)的相關，我們稱為淨自我相關，PACF。 
    • pacf(Y1)

  • 對於一個時間數列，我們需要嘗試不同的\(\phi_{1}\)，從0到0.95都可能，才可以確定到底是deterministic還是stochastic。
  • 自我相關可以是隔好幾個時間差，也就是AR(\(p\))：

    • \[\begin{align*} y_{t}&=y_{t-1}\phi_{1}+y_{t-2}\phi_{2}+\ldots+y_{t-p}\phi_{p}+\epsilon_{t} \end{align*}\]

  • 如果是AR(1)模型，我們可以用最小平方法估計\(Y_{t}\)跟\(Y_{t-1}\)之間的關係。
    • 例如，前一期的通貨膨脹變動率可以預測當期的通貨膨脹變動率。
    • 如果是AR(p)的迴歸模型可寫成：

    \[Y_{t}=\beta_{0}+\beta_{1}Y_{t-1}+\beta_{2}Y_{t-2}+\ldots+\beta_{p}Y_{t-p}+\epsilon_{t}\]

    • 我們可以用\(F\)檢定係數全部是0的虛無假設。

    4.5 MA模型

  • 如果數列受到不同期的誤差項影響，稱為移動平均(moving average, MA)模型。例如目前的經濟受到去年疫情的意外影響，或者今年的身體健康狀況受到去年一場車禍的影響等等，也就是誤差項的意外(shock)。MA模型可寫成：

    • \[Y_{t}=\epsilon_{t-1}\phi_{1}+\epsilon_{t-2}\phi_{2}+\ldots+\epsilon_{t-q}\phi_{q}+\epsilon_{t}\]

  • 在AR模型，\(y_{t}\)與\(y_{t-k}\)相關，但是在MA模型，因為意外影響在一期之後消失，所以\(y_{t}\)與\(y_{t-k}\)不相關。
    • 模擬MA(1)的數列，其中\(\phi_{1}=0.57\)。
    set.seed(02138)
y <- arima.sim (list(order=c(0,0,1),
                ar = NULL,
                ma = c(0.57)),
                n= 1000)    
plot(y, col='#DD22AA')

    • 畫ACF圖。
    acf(y)

    • 畫PACF圖。
    pacf(y)

    4.6 週期模型

  • 週期模型假設有類似特徵的時間點每隔一些時間點就會出現，例如\(t=1\)跟\(t=13\)隔了12個時間點有類似的影響，\(t=2\)跟\(t=14\)。但是我們不能排除\(t=1\)跟\(t=2\)之間也有相關。
    1. 相加(additive)模型：每一個循環同一個時間點增加\(k_{k}\)單位，例如\(k_{1}=15\)，\(k_{2}=-5\)，\(\bar{k}=6\)等等。\(k_{1}-\bar{k}=9\)，表示1月比其他月份平均來講多了9單位。
    • 假如我們知道真實的\(k\)，我們可以對每一個觀察值減去\(k\)。我們也可以用最小平方法，以某一個月作為參考組，估計\(k-1\)參數。
    • 如果我們確定時間數列有明顯的循環，而且沒有確定性模型的趨勢，數值之間有固定的離散程度，那麼資料適用於相加模型。
    1. 倍數(multiplicative)模型：每一個循環同一個時間點增加\(p\)倍單位，0<p<1。結合AR與MA模型的ARMA模型有類似的模式，但是循環的時間點比較多，例如ARMA(12,0)代表\(y_{t}\)受到\(y_{t-12}\)的影響，\(\phi_{1}=\phi_{2}=\ldots=\phi_{11}=0\)，假設\(\phi_{12}=0.9\)，我們模擬5年也就是60個月的時間數列。
    adolph='http://faculty.washington.edu/cadolph/panEssex/makeTS.R'
source(adolph)
set.seed(02138)
Y <- makeTS(n= 120,
             ar= c(rep(0,11),0.85), #phi=0.85
             ma= NULL,
            trend=0.001,
            seasons=c(0,0,0,0,0,
                      0,0,0,0,0,
                      0,0),
            adjust=c('level'),
            x = NULL,
            beta = NULL,
              varY = 1, #Error Term variance
initY = 0,
initE = 0,
burnin = 10)
plot.ts(Y, col='#CC22AA')
    週期性時間數列圖1

    Figure 4.10: 週期性時間數列圖1

    • 我們把120個月的資料轉成10年也就是10條線，圖4.11顯示，10條線出現循環的現象。
    Y.ts<-ts(Y, start=1, frequency = 12)

t<-c(1: 12)
plot(t, Y.ts[1:12], type='n')

dt <- matrix(Y.ts, ncol=12, nrow=length(Y.ts)/12,
             byrow=TRUE)

for (i in 1: 10){
  lines(t, dt[i,], type='l', lty=i, lwd=1.5,
    col=RColorBrewer::brewer.pal(10, 'BuPu')[i])
}
    週期性時間數列圖2

    Figure 4.11: 週期性時間數列圖2

  • 確定時間數列之後，我們可以就不同類型的數列進行多變量分析，找出其他變數與時間數列變數的關係。
  • 如果假設時間數列資料的誤差項跟前一期的誤差項相關，可以用一階自我迴歸模型(first-order autoregressive model)。
  • 自我迴歸可寫成：
    • \[\begin{equation} y_{t}=\beta_{0}+\beta_{1}y_{t-1}+\epsilon_{t}\quad \epsilon_{t}\sim (0,\sigma^2) \#(eq:auto) \end{equation}\]
  • \(y_{t}\)的期望值：\(E[y_{t}]=\mu\quad if\quad |\beta_{1}|<1\)
  • 我們還可以使用滯後回歸(Lagged Regression)，簡單來說就是將當期變數與上X期的資料去做相關，舉例來說我們想創建一個落後一期的迴歸模型，使用最開始的lm()線性迴歸即可。

    • 5 使用forecast套件分析

  • 接下來參考網站分析資料。首先畫完整的資料趨勢如圖5.1： 
  • plot.ts(dat)
    正新橡膠營業收入趨勢

    Figure 5.1: 正新橡膠營業收入趨勢

  • 畫四季的移動平均趨勢如圖5.2，得到的結果如同圖3.1： 
    • library(TTR)
datsma3 <-TTR::SMA(dat, n=4)
plot.ts(datsma3)
    TTR產出的四季移動平均線

    Figure 5.2: TTR產出的四季移動平均線 

    library(forecast)
trend_dat<-forecast::ma(dat_ts, order=4, centre=T)
kbl(trend_dat, align='l', caption = "四季移動平均值")%>%
  kable_styling('striped',font_size = 20, position='left')
    Table 5.1: 四季移動平均值
    x
    NA
    NA
    293.4
    293.7
    293.6
    293.9
    293.3
    289.3
    284.6
    281.9
    279.5
    278.1
    277.1
    275.0
    NA
    NA
  • 再來用forecast套件畫圖 5.3，比對圖 3.1。 
    • plot(as.ts(trend_dat))
    forecast產出的移動平均線與長期趨勢

    Figure 5.3: forecast產出的移動平均線與長期趨勢 

    plot(as.ts(dat_ts))
lines(trend_dat, col='tomato')
    forecast產出的移動平均線與長期趨勢

    Figure 5.4: forecast產出的移動平均線與長期趨勢

  • decompose指令可以求出移動平均(trend)、季節因子(figure)、季節指數(seasonal)等等。 
    • xt<-decompose(dat_ts, 'multiplicative')
kbl(data.frame(xt$trend), caption='decompose功能')%>%
          kable_styling('striped', font_size=18, position='center', full_width = F)
    Table 5.2: decompose功能
    xt.trend
    NA
    NA
    293.4
    293.7
    293.6
    293.9
    293.3
    289.3
    284.6
    281.9
    279.5
    278.1
    277.1
    275.0
    NA
    NA 
    kbl(data.frame(xt$figure), caption='figure')%>%
            kable_styling('striped', font_size=18, position='center', full_width = F)
    Table 5.2: figure
    xt.figure
    1.0048
    1.0297
    1.0028
    0.9628 
    kbl(data.frame(xt$seasonal), caption='Seasonal')%>%
        kable_styling('striped',font_size=18, position='center', full_width = F)
    Table 5.2: Seasonal
    xt.seasonal
    1.0048
    1.0297
    1.0028
    0.9628
    1.0048
    1.0297
    1.0028
    0.9628
    1.0048
    1.0297
    1.0028
    0.9628
    1.0048
    1.0297
    1.0028
    0.9628
  • 我們也可以直接畫圖decompose指令所得到的可以求出移動平均(trend)、季節因子(figure)、季節指數(seasonal)等等。 
    • plot(xt)

    6 美國的GDP資料分析

  • 除了forecast, TTR之外，還可以使用\(\tt{zoo}\), \(\tt{quantmod}\)等套件分析時間序列資料。
  • 我們參考Christoph Hanck的說明，使用美國的GDP資料，示範用\(\tt{quantmod}\), \(\tt{zoo}\)分析時間序列資料。
  • \(\tt{zoo}\)可以根據資料對應的時間或其他順序排序向量，例如： 
    • x <- zoo(rnorm(5), as.Date("2021-02-01")+seq(1,10, by=2) )
kbl(x) %>% kable_styling('striped',font_size = 18)
    x
    2021-02-02 1.2864
    2021-02-04 0.1515
    2021-02-06 0.1596
    2021-02-08 -0.1517
    2021-02-10 -0.3866
  • \(\tt{quantmod}\)特別適用於財務分析、股票分析，例如建立超前一期或是落後一期的資料、建立移動平均線的圖形等等。 
    • library(AER);library(quantmod)
data('NYSESW')
NYSESW <- xts(Delt(NYSESW))
plot(as.zoo(NYSESW),
     col = "steelblue",
     lwd = 2,
     ylab = "Percent per Day",
     xlab = "Date",
     main = "New York Stock Exchange Composite Index",
     cex.main = 1)

  • 首先，我們讀取資料，然後把資料轉換成年度與季的格式，並且更換每一欄資料的名稱。 
    • #load libraries
library(readxl)
#specify path
usquar<-here::here('data','us_Macro_Quarterly.xlsx')
# load US macroeconomic data
USMacroSWQ <- read_xlsx(usquar,
              sheet = 1,
             col_types = c("text", rep("numeric", 9)))

# format date column
USMacroSWQ$...1 <- zoo::as.yearqtr(USMacroSWQ$...1, format = "%Y:0%q")

# adjust column names
colnames(USMacroSWQ) <- c("Date", "GDPC96", "JAPAN_IP", "PCECTPI",  "GS10", "GS1", "TB3MS", "UNRATE", "EXUSUK", "CPIAUCSL")
  • 接著，我們選擇1960到2013的資料，並且計算GDP成長的數列。因為經濟成長資料取自然對數之後，會接近線性，而且標準差會接近常數，所以在此取自然對數。
  • \(Y_{t}\)與\(Y_{t-1}\)各取自然對數之後，兩者之間的差異可寫成\(\Delta\rm{log}(Y_{t})=\rm{log}(Y_{t})-\rm{log}(Y_{t-1})\)，計算方式為\(\rm{log}(Y/\rm{lag}(Y))\)
  • \(100\Delta\rm{log}(Y_{t})\)等於\(Y_{t}\)與\(Y_{t-1}\)以百分比表示的變動。 
    • # GDP series as xts object
GDP <- xts(USMacroSWQ$GDPC96, USMacroSWQ$Date)["1960::2013"]

# GDP growth series as xts object
GDPGrowth <- xts(400 * log(GDP/stats::lag(GDP)))
    1. 畫圖表示GDP的成長曲線： 
    plot(log(as.zoo(GDP)),
     col = "steelblue",
     lwd = 2,
     ylab = "Logarithm",
     xlab = "Date",
     main = "U.S. Quarterly Real GDP")
    GDP成長曲線

    Figure 6.1: GDP成長曲線

  • 圖 6.1 顯示，從1960年代到2010年代，美國的國民生產毛額越來越高。 
    • plot(as.zoo(GDPGrowth),
     col = "steelblue",
     lwd = 2,
     ylab = "Logarithm",
     xlab = "Date",
     main = "U.S. Real GDP Growth Rates")
    GDP成長率之成長曲線

    Figure 6.2: GDP成長率之成長曲線

  • 圖 6.2 顯示，取自然對數後的每一季的國民生產毛額的成長率有高有低。我們可以自定函數，產生如表6.1。注意年度的成長率要乘以400，因為一年有四季。 
    • # compute logarithms, annual growth rates and 1st lag of growth rates
quants <- function(series) {
  s <- series
  return(
    data.frame("Level" = s,
               "Logarithm" = log(s),
               "AnnualGrowthRate" = 400 * log(s / stats::lag(s)),
               "1stLagAnnualGrowthRate" = stats::lag(400 * log(s / stats::lag(s))))
    )
}
# obtain a data.frame with level, logarithm, annual growth rate and its 1st lag of GDP
Qd<-quants(GDP["2011-07::2013-01"])
kbl(Qd, caption='GDP自然對數') %>%
  kable_styling(font_size = 20, position='left' )
    Table 6.1: GDP自然對數
    Level Logarithm AnnualGrowthRate X1stLagAnnualGrowthRate
    2011 Q3 15062 9.620 NA NA
    2011 Q4 15242 9.632 4.7518 NA
    2012 Q1 15382 9.641 3.6422 4.7518
    2012 Q2 15428 9.644 1.1972 3.6422
    2012 Q3 15534 9.651 2.7470 1.1972
    2012 Q4 15540 9.651 0.1453 2.7470
    2013 Q1 15584 9.654 1.1392 0.1453

    7 結論

  • 時間數列資料可拆解為長期趨勢、季節變動、隨機變動等等。
  • 時間數列資料可分為相加模型以及相乘模型。前者假設季節變動為固定。
  • 時間數列也可以分為deterministic, stochastic, seasonlity等模型，其中stochastic 模型又可分為AR, MA等模型。
  • decompose指令可以求出移動平均(trend)、季節因子(figure)、季節指數(seasonal)等等。

    • 8 作業

    1. 請畫圖分析fpp套件中的cafe這筆資料的移動平均以及原始資料，並且求出季節因子，請問哪一季花在餐廳、外帶的費用最高？

    2. 在R的資料庫中，lynx這筆資料記錄西元1821至1934年的加拿大西北部馬更些(Mackenzie)流域附近每年的山貓數量。請嘗試畫出加拿大山貓(lynx)這筆資料的時間數列圖、自我相關函數圖，以及淨自我相關函數圖。

    3. 請練習現代統計學第二版第16.4題。（提示：2019年的資料如表8.1；趨勢值乘以季節指數等於預測值。）

    d2019<-c(41,42,43,44)
kbl(d2019, caption='某地區2019年資料') %>%
  kable_styling(position='center',full_width = F, font_size = 18)
    Table 8.1: 某地區2019年資料
    x
    41
    42
    43
    44

    4. 請練習現代統計學第二版第16.12題。（提示：先轉換表格資料為向量，然後設定12個時間點的自變數\(t\)，以及設定3個季節的虛擬變數，再以最小平方法估計模型。）

    5. 請練習現代統計學第二版第16.14題。（提示：比照前一題整理資料，然後以最小平方法估計模型，最後以Year=11預測2020年的平均單價。）

    9 更新日期

    最後更新日期： 06/17/2021