三葉草之家: 6月 2021

2021年6月11日星期五

R：使用ARIMA對股價做預測分析

時間序列分析是根據系統觀測得到的時間序列數據，通過曲線擬合和參數估計來建立數學模型的理論和方法。時間序列分析常用在企業經營管理、市場潛量預測、農作物病蟲災害預報、環境污染控制、生態平衡、天文學和海洋學等方面。

今天來針對台灣的護國神山股價的走勢作預測分析，雖然購買股票聽起來很有誘惑力，但我們更應該進行詳細的深入分析，避免以投機為基礎的股票購買。

1.載入相關套件
library(PerformanceAnalytics)
library(quantmod)
library(stringr)
library(forecast)

2.取得TSMC自 2020/1/1 開始的股價 (yahoo的台股股價名稱要用代號.TW)
getSymbols("2330.TW", src="yahoo", from="2020-01-01")
TSMC.TW <- get("2330.TW")
將第四欄位(收盤價)存在變數TSMC_close_price
TSMC_close_price <- TSMC.TW[, 4]

3.對TSMC收盤價做 auto arima，存於變數 fit.a
fit.a <- auto.arima(TSMC_close_price, seasonal=F)
可以先看一下時間序列分析圖
tsdisplay(residuals(fit.a), lag.max=40, main="(3,1,4 Model Residuals)")

最上圖是時間序列的趨勢圖，判斷趨勢是否平穩，
下圖左ACF是畫自相關圖，右邊PACF是畫偏相關圖，用來判斷拖尾和截尾。

4.對 fit.a 做10天的預測(forecast)，結果存於 fc1
fc1 <- forecast(fit.a, h=7)

畫出預測趨勢圖
plot(fc1)

fc1的預測低值/高值可以輸入fc1$lower 或 fc1$upper 查看:

ACF和PACF概念：
ACF（自相關系數）描述的是現在股價與過去某時期之間的所有價格是否相關。PACF（偏自相關系數）描述的是現在的價格單純地和過去某一個價格的相關關系。兩者均為取值範圍 -1~1，絕對值越靠近1，說明相關關系越明顯；絕對值越靠近0，說明兩者之間線性相關很弱。藍色虛線表示的是誤差範圍，當值沒有超過藍色虛線的時候，基本可以認為他們不相關。如：

par(mfrow = c(1, 2))
acf(TSMC.TW$`2330.TW.Close`, main="Series tsmc")
pacf(TSMC.TW$`2330.TW.Close`, main="Series tsmc")
auto.arima(TSMC_close_price, seasonal = F)

對於 MA 模型，我們將使用 ACF 圖來識別 (q) 階，並且 PACF 將呈指數衰減。
如果我們查看 PACF 圖，我們可以注意到它僅在第一次滯後時有一個顯著的尖峰，
這意味著所有高階自相關都可以通過第一次滯後自相關有效地解釋。

如果要尋找最優模型：

先列出 AUTO ARIMA模型：
auto.arima(TSMC_close_price, trace=T)

雖然購買股票聽起來很有誘惑力，但我們更應該進行詳細的深入分析，避免以投機為基礎的股票購買。

原文網址：https://kknews.cc/invest/apy6rlv.html

data_d1=diff(TSMC_close_price, 1) #進行一階差分
tsdisplay(data_d1) #進行一階差分得到acf和pacf
data_d2=diff(TSMC_close_price, 2) #進行二階差分
tsdisplay(data_d2) #得到acf和pacf

由上圖可以看出二階差分的PACF有明顯 spike，由於我們使用 AUTO ARIMA 函數為我們提供了更好的數據集方法，因此我們就可以不用再花時間去找模型參數。

我們可以利用這個方式來檢驗過去的資料，例如在2021年5月份台股經歷一次大跌，TSMC股價也是無法倖免，在5/13最低價$547，我們可以觀測TSMC 2020/1/1到2021/5/2股價，然後預測2021/5/3~2021/6/11 的股價(30天)，驗證是否落在真實的股價區間。

將觀測區間股價轉成 ts 時間序列(time series)
tsmc.ts <- ts(TSMC_close_price)
時間序列長度就是觀測區間的天數(348天)，存於test.end.ts
(test.end.ts <- length(TSMC_close_price))
將天數往前30天(第318天)
(test.start.ts <- end.ts - 30)
測試資料從318~348天
stock_test<- window(tsmc.ts, start=test.start.ts, end=test.end.ts)
訓練模型從第一天(筆)資料至317天
stock_train <- window(tsmc.ts, start=1, end=(test.start.ts - 1))
對 stock_train 做auto arima
(fit <- auto.arima(stock_train))
然後forecast 30天的資料
fc <- forecast(fit, h=30)
使用自動擬合模型並檢驗殘差
checkresiduals(fit)
繪圖
plot(fc)
將預測期間的股價用紅色線段繪出
lines(stock_test, col="red")

利用ARIMA(0,1,1) 模型預測30個值(藍線)，與之前我們預留的30個觀測值(紅線)進行對比，如下圖：

自動擬合模型並檢驗殘差(如上圖)

我們可以看出在5月份股價雖然跌至547，仍然在我們做的預測模型區間內，我們可以看forecast的upper與lower值，最低預測在539：

我們發現雖然模型從統計學意義上來說是擬合的。預測表明TSMC股價應該呈現上漲趨勢，但是由於突發事件的影響，預測值與觀測值的趨勢並沒有很好地重合，但是區間是吻合的。

2021年6月9日星期三

R：使用ARIMA對COVID-19做時間序列分析

ARIMA模型（Autoregressive Integrated Moving Average model），差分整合移動平均自迴歸模型，又稱整合移動平均自回歸模型，利用這個模型可以對時間序列資料做預測分析，像是股票價格，營業銷售預測，今天就來預測 COVID-19 的疫情。

首先，前面幾篇已經談到如何將約翰·霍普金斯大學的疫情統計抓資料下來，程式碼:

library(tidyverse)
library(ggplot2)
library(cowplot)
library(patchwork)
library(forecast)
path <- "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/"
confirm_file <- "time_series_covid19_confirmed_global.csv"
deaths_file <- "time_series_covid19_deaths_global.csv"
recovered_file <- "time_series_covid19_recovered_global.csv"

confirmed_df <- read_csv(paste0(path, confirm_file))
deaths_df <- read_csv(paste0(path, deaths_file))
recovered_df <- read_csv(paste0(path, recovered_file))

range <- 60
sd <- length(confirmed_df) - range; sd
ed <- length((confirmed_df)); ed
dates <- colnames(confirmed_df[, sd:ed]); dates
dates <- as.Date(dates,format = "%m/%d/%y"); dates

#date.list <- seq(from=as.Date(dates[1]),
#                 to=as.Date(dates[range+1]), by=1); date.list
getCountrydata <- function(Country,
                           dates = dates,
                           confirmed_df = confirmed_df,
                           deaths_df = deaths_df,
                           recovered_df = recovered_df,
                           sd = sd, ed = ed) {
if (Country == "all") {
    cases <- confirmed_df %>%
      #select(-(1:400)) %>%
      select(sd:ed) %>%
      colSums()
    death <- deaths_df %>%
      #select(-(1:400)) %>%
      select(sd:ed) %>%
      colSums()
    recovered <- recovered_df %>%
      #select(-(1:400)) %>%
      select(sd:ed) %>%
      colSums()
}
else {
    Country <- enquo(Country)
    cases <- confirmed_df %>%
      filter(`Country/Region` == !! Country) %>%
      #select(-(1:400)) %>%
      select(sd:ed) %>%
      colSums()
    death <- deaths_df %>%
      filter(`Country/Region` == !! Country) %>%
      #select(-(1:400)) %>%
      select(sd:ed) %>%
      colSums()
    recovered <- recovered_df %>%
      filter(`Country/Region` == !! Country) %>%
      #select(-(1:400)) %>%
      select(sd:ed) %>%
      colSums()
}
res.df <- tibble(dates,
                   cases = cases,
                   death = death,
                   recovery = recovered,
                   mortality_rate = death/cases,
                   recovery_rate = recovery/cases)
return(res.df)
}
world.df <- getCountrydata(Country = "all",
                           dates = dates,
                           confirmed_df = confirmed_df,
                           deaths_df = deaths_df,
                           recovered_df = recovered_df, sd, ed)
#Taiwan
taiwan.df <- getCountrydata(Country = "Taiwan*",
                            dates = dates,
                            confirmed_df = confirmed_df,
                            deaths_df = deaths_df,
                            recovered_df = recovered_df, sd, ed)
getdailycase <- function(comfirmed_df, new.df, country, sd, ed, dates) {
sd <- sd - 1
ed <- sd
dates <- colnames(confirmed_df[, sd:ed])
df <- getCountrydata(Country = country,
                       dates = dates,
                       confirmed_df = confirmed_df,
                       deaths_df = deaths_df,
                       recovered_df = recovered_df, sd, ed)
first_data <- df$cases[[1]] ; first_data
new.df['daily'] <- NA;
for(i in 1:nrow(new.df)) {
    if(i == 1)
      new.df$daily[i] <- new.df$cases[i]- first_data
    else
      new.df$daily[i] <- new.df$cases[i]- new.df$cases[i - 1]
}
new.df
}
taiwan.df <- getdailycase(confirmed_df, taiwan.df, "Taiwan*", sd, ed, dates)

我們取最近60天的資料來做分析，資料名稱為 taiwan.df，輸入tail(taiwan.df)可以看到最後6筆(天)的資料：

接著要對這資料做一些處理，
1.載入zoo 套件包，以便將資料轉成可視的時間序列
library(zoo)
2.將taiwan.df的 death (死亡組數欄位)依日期做成zoo.ts 時間序列
zoo.ts <- zoo(as.numeric(taiwan.df$death), dates)
3.自動ARIMA模型分析
arc <- auto.arima(zoo.ts)
summary(arc)

4.對分析後的數據 arc 預測4天後的資料存於 fc1
fc1 <- forecast(arc, h=4)

5.畫圖：
plot(fc1, xlab="Date", ylab="Comfirmed Cases")

我們可以看到預測往後四天死亡總人數最高最低的預測數字。

因為這圖是對時間序列來繪圖，因此要對這圖表的時間軸轉成日期的格式。

先把圖的x軸座標隱藏
plot(fc1, xaxt="n", xlab="Date", ylab="Death")
用函數 time() 將資料zoo.ts的日期取出存於tt
tt <- time(zoo.ts)
將這日期依間隔10天做成一個向量
ds.lst <- seq(1, length(tt), by=10)
日期格式變數為 YYYY-MM 的格式
fmt <- "%Y-%m"
依時間格式將日期轉成字串標籤
labs <- format(tt[ds.lst], fmt)
畫出X軸的格線
axis(side=1, at = tt[ds.lst], labels=F)
將日期用text函數補上去，角度45度，字體縮小0.8
text(x=tt[ds.lst], y=par()$usr[3]-15, labels=labs,
srt=45, adj=1, cex=0.8, xpd=TRUE)

下面是對每日確診數做7日的預測:

library(zoo)
zoo.ts <- zoo(as.numeric(taiwan.df$cases), dates)
arc <- auto.arima(zoo.ts)
summary(arc)
fc1 <- forecast(arc, h=7)
plot(fc1, xaxt="n", xlab="Date", ylab="Comfirmed Cases")
tt <- time(zoo.ts)
ds.lst <- seq(1, length(tt), by=10)
fmt <- "%Y-%m"
labs <- format(tt[ds.lst], fmt)
axis(side=1, at = tt[ds.lst], labels=F)
text(x=tt[ds.lst], y=par()$usr[3]-500, labels=labs,
srt=45, adj=1, cex=0.8, xpd=TRUE)

輸出結果：