系統設置

#載入packages
library(tseries)
library(dplyr)
library(ggplot2)
library(forecast)

對原始資料進行趨勢分析

資料來源:行政院環境資源資料庫

由時間序列圖觀察可發現，我們的數據具有週期性變化，且有下降的趨勢，並非一直維持在某一個常數上下波動，初步先判斷該序列爲非平穩時間序列。

也可以另外從ACF圖觀察，當Lag=12、24、36、48時可以預期具有相隔12期的季節性循環且遞減緩慢，因此我們判斷應該對資料進行差分，轉換為平穩型序列。

TSP<- read.csv("C:/Users/a1205/Desktop/TSP.csv")
tsTSP<- ts(TSP$總懸浮微粒,frequency = 12,start=c(2001,1))
plot(tsTSP)

acf(TSP$總懸浮微粒, lag.max = 50,plot = TRUE)

轉換為平穩型序列

對原始資料進行相隔12期差分一次，所得結果如下表，可由觀測值的趨勢圖及ACF圖看出資料已轉換為平穩型。

data_d1=diff(tsTSP,12,1)
par(mfrow=c(1,2))
plot(tsTSP);plot(data_d1)  

adf.test(data_d1)
## Warning in adf.test(data_d1): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  data_d1
## Dickey-Fuller = -5.3715, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

建立與比較模型

模式鑑定方式:

  • AR模式 : ACF遞減，PACF切斷
  • MA模式 : ACF切斷，PACF遞減
  • ARMA模式 : ACF遞減，PACF遞減

我們觀察到ACF表的第12期顯著、PACF第12、20、24、34、36期顯著，由此判斷ACF切斷、PACF逐漸遞減。 因此，我們考慮先建立兩個模型，分別為季節性ARIMA過程(0,1,1)[12]與變化型的季節性ARIMA過程(20,0,0)(0,1,1)[12]，再將兩個模型進行比較何者模式較佳。

tsdisplay(data_d1,xlab="year",ylab="TSP")

我們觀察兩模式，發現模型2的AIC較小，且殘差自相關檢定的部分，可以發現模型二p值為不顯著，故殘差為隨機跳動。

結論:選擇ARIMA(20,0,0)(0,1,1)[12]。

fit1 <- arima(TSP$總懸浮微粒, order = c(0, 0, 0), seasonal = list(order = c(0, 1, 1), period=12))
fit2 <- arima(TSP$總懸浮微粒, order = c(20, 0, 0), seasonal = list(order = c(0, 1, 1), period=12))
fit1
fit2
Box.test(residuals(fit1),type="Box-Pierce")
Box.test(residuals(fit2),type="Box-Pierce")
## 
## Call:
## arima(x = TSP$總懸浮微粒, order = c(0, 0, 0), seasonal = list(order = c(0, 1, 
##     1), period = 12))
## 
## Coefficients:
##          sma1
##       -0.3564
## s.e.   0.0555
## 
## sigma^2 estimated as 142.3:  log likelihood = -795.96,  aic = 1595.93
## 
## Call:
## arima(x = TSP$總懸浮微粒, order = c(20, 0, 0), seasonal = list(order = c(0, 
##     1, 1), period = 12))
## 
## Coefficients:
##          ar1     ar2     ar3     ar4     ar5     ar6     ar7      ar8     ar9
##       0.2169  0.0896  0.0648  0.0584  0.0028  0.0250  0.0478  -0.0523  0.0435
## s.e.  0.0702  0.0709  0.0725  0.0732  0.0717  0.0737  0.0714   0.0705  0.0713
##         ar10    ar11     ar12    ar13    ar14     ar15     ar16    ar17    ar18
##       0.0211  0.1816  -0.0061  0.0813  0.0979  -0.0178  -0.0159  0.0398  0.1297
## s.e.  0.0703  0.0714   0.0822  0.0711  0.0719   0.0741   0.0755  0.0733  0.0749
##         ar19    ar20     sma1
##       0.0922  -0.105  -0.9053
## s.e.  0.0730   0.072   0.0857
## 
## sigma^2 estimated as 81.76:  log likelihood = -749.47,  aic = 1542.93
## 
##  Box-Pierce test
## 
## data:  residuals(fit1)
## X-squared = 3.5113, df = 1, p-value = 0.06095
## 
## 
##  Box-Pierce test
## 
## data:  residuals(fit2)
## X-squared = 0.0040926, df = 1, p-value = 0.949