Untitled

1 時序分析步驟

• 觀察原始資料

x <- lynx
logx <- log(x)
plot(x,type = 'o')

qqnorm(x,pch = 16,main = "x")
qqline(x)

qqnorm(logx,pch = 16,main = "log x")
qqline(logx)

par(mfrow = c(1,2),oma = c(0,0,0,0))
acf(x);pacf(x)

• acf上下震盪遞減至0，代表有週期性

• 沒有趨勢，有週期，約10年左右為一個循環

plot(log(x))

par(mfrow = c(1,2),oma = c(0,0,0,0))
acf(log(x));pacf(log(x))

fit_ar <- arima(logx,order = c(11,0,0))
fit_ar

## 
## Call:
## arima(x = logx, order = c(11, 0, 0))
## 
## Coefficients:
##          ar1      ar2     ar3      ar4     ar5      ar6     ar7      ar8
##       1.1676  -0.5446  0.2661  -0.3093  0.1540  -0.1463  0.0569  -0.0294
## s.e.  0.0877   0.1394  0.1474   0.1489  0.1509   0.1498  0.1511   0.1483
##          ar9    ar10     ar11  intercept
##       0.1346  0.2021  -0.3394     6.6678
## s.e.  0.1455  0.1383   0.0885     0.1086
## 
## sigma^2 estimated as 0.1915:  log likelihood = -70.07,  aic = 166.13

plot(x,type = 'p',main = "fitted value");lines(fitted(fit_ar) %>% exp)

Box.test(fit_ar$residuals,type = "Ljung")

## 
##  Box-Ljung test
## 
## data:  fit_ar$residuals
## X-squared = 0.088022, df = 1, p-value = 0.7667

par(mfrow = c(2,2),oma = c(0,0,0,0))
acf(fit_ar$residuals)
pacf(fit_ar$residuals)

par(mfrow = c(2,2),oma = c(0,0,0,0))
qqnorm(as.vector(fit_ar$residuals),datax = T,pch = 16,main = "",xlab = "Residual")
qqline(as.vector(fit_ar$residuals),datax = T)

plot(as.vector(fitted(fit_ar)),as.vector(fit_ar$residuals),xlab = "Fitted value",ylab = "Residual")
abline(h = 0)

hist(as.vector(fit_ar$residuals),xlab = "Residual",main = "")

plot(as.vector(fit_ar$residuals),type = "o",ylab = "Residual")
abline(h = 0)

2 預測

f_ar <- forecast(fit_ar,h = 20) %>% as.array()


plot(x,type = "o",xlim = c(1820,1955),main = "forecasts")
lines(1935:1954,f_ar$lower[,2] %>% exp,lty = 2)
lines(1935:1954,f_ar$upper[,2] %>% exp,lty = 2)
lines(1935:1954,f_ar$mean %>% exp,col = 2)

3 用差分配ar

logx_diff <- diff(logx,lag = 11)
fit_ar <- arima(logx_diff,order = c(11,0,0))

plot(logx_diff %>% exp,type = 'p',main = "fitted value");lines(fitted(fit_ar) %>% exp)

# Box.test(fit_ar$residuals,type = "Ljung")

par(mfrow = c(2,2),oma = c(0,0,0,0))
acf(fit_ar$residuals)
pacf(fit_ar$residuals)

par(mfrow = c(2,2),oma = c(0,0,0,0))
qqnorm(as.vector(fit_ar$residuals),datax = T,pch = 16,main = "",xlab = "Residual")
qqline(as.vector(fit_ar$residuals),datax = T)

plot(as.vector(fitted(fit_ar)),as.vector(fit_ar$residuals),xlab = "Fitted value",ylab = "Residual")
abline(h = 0)

hist(as.vector(fit_ar$residuals),xlab = "Residual",main = "")

plot(as.vector(fit_ar$residuals),type = "o",ylab = "Residual")
abline(h = 0)

4 決定 p & q 值

5 尚未理解的問題

• acf 及 pacf 所包含的資訊
• Box.test 的 fitdf 代表甚麼
• 差分的使用
• ar(p) ma(q) , p、q 更嚴謹的求法
• arima 如何從 acf 及 pacf 看出 p、d、q