#======================================================
#	10 時系列分析
#======================================================
#======================================================
#	10.1 ARIMAモデル
#======================================================
#パッケージのロード
library(forecast)

#-----------<定常性、自己共分散及び自己相関>-----------
#-----------<ホワイトノイズ>-----------

set.seed(456)
w <- rnorm(200)    #ホワイトノイズ（正規乱数の場合）
x <- cumsum(w)     #ランダム・ウォーク
wd <- w + 0.2      #ホワイトノイズ＋ドリフト
xd <- cumsum(wd)   #ドリフト付きランダム・ウォーク

#-----------<時系列プロット>-----------

#図10.1
op <- par(mfrow=c(2, 1), mar=c(4.1, 4.1, 3.1, 1.1))
plot.ts(w, ylim=c(-10, 10), main="ホワイトノイズ")
plot.ts(xd, ylim=c(-10, 40), 
        main="ランダム・ウォークとドリフト付きランダム・ウォーク")
lines(x)
abline(0, 0.2, lty=2)
par(op)

#-----------<ARモデル>-----------
#-----------<AR(1)>-----------

set.seed(456)
ar1.1 <- arima.sim(n=1000, model=list(ar=c(0.9)))
ar1.2 <- arima.sim(n=1000, model=list(ar=c(-0.9)))

#-----------<ACFプロット>-----------

#図10.2
op <- par(mfrow=c(2, 1), mar=c(4.1, 4.1, 3.1, 1.1))
acf(ar1.1, lag.max=30, main=expression(AR(1)~~~phi[1]==+0.9))
acf(ar1.2, lag.max=30, main=expression(AR(1)~~~phi[1]==-0.9))
par(op)

#-----------<特性方程式>-----------

(phi.z <- polyroot(c(1, -0.9, -0.5))) #多項式の根
Mod(phi.z)                            #根の絶対値
arima.sim(n=1000, model=list(ar=c(0.9, 0.5)))

#-----------<根が実数の場合と複素数の場合>-----------

set.seed(456)
(phi.re <- polyroot(c(1, -0.6, -0.3))) #実数根
ar2.1 <- arima.sim(n=1000, model=list(ar=c(0.6, 0.3)))
(phi.im <- polyroot(c(1, -0.6, 0.3)))  #複素根
ar2.2 <- arima.sim(n=1000, model=list(ar=c(0.6, -0.3)))

#-----------<ACFプロット>-----------

#図10.3
op <- par(mfrow=c(2, 1), mar=c(4.1, 4.1, 3.1, 1.1))
acf(ar2.1, lag.max=30, 
       main=expression(AR(2)~~~phi[1]==+0.6~~phi[2]==+0.3))
acf(ar2.2, lag.max=30, 
       main=expression(AR(2)~~~phi[1]==+0.6~~phi[2]==-0.3))
par(op)

#-----------<PACFプロット>-----------

#図10.4
op <- par(mfrow=c(2, 1), mar=c(4.1, 4.1, 3.1, 1.1))
pacf(ar1.1, lag.max=30, main=expression(AR(1)~~~phi[1]==+0.9))
pacf(ar2.1, lag.max=30, 
       main=expression(AR(2)~~~phi[1]==+0.6~~phi[2]==+0.3))
par(op)

#-----------<ar.mleの当てはめ>-----------

ar2.ar <- ar.mle(ar2.2) #当てはめ
ar2.ar                  #推定結果
ar2.ar$aic              #AIC

#-----------<MAモデル>-----------
#-----------<MA(2)>-----------

set.seed(456)
ma2 <- arima.sim(n=1000, model=list(ma=c(0.6, -0.3)))

#-----------<ACFプロットとPACFプロット>-----------

#図10.5
op <- par(mfrow=c(2, 1), mar=c(4.1, 4.1, 3.1, 1.1))
acf(ma2, lag.max=30, 
       main=expression(MA(2)~~~theta[1]==+0.6~~theta[2]==-0.3))
pacf(ma2, lag.max=30, 
       main=expression(MA(2)~~~theta[1]==+0.6~~theta[2]==-0.3))
par(op)

#-----------<ARIMAモデルとSARIMAモデル>-----------
#-----------<差分>-----------

#図10.6
op <- par(mfrow=c(2, 1), mar=c(4.1, 4.1, 3.1, 1.1))
plot(LakeHuron, main="原系列")
plot(diff(LakeHuron), main="差分系列")
par(op)

#-----------<ARIMAモデルとauto.arima>-----------

Huron.0 <- window(LakeHuron, end=1962) #最後の10観測値を除く
Huron.arima <- auto.arima(Huron.0, seasonal=F, allowdrift=F) #モデル選択
summary(Huron.arima)           #要約

#図10.7
tsdiag(Huron.arima)            #モデル診断プロット

Huron.fore <- forecast(Huron.arima, 10) #10期先までの予測
Huron.fore

#図10.8
windows(10, 7)
plot(Huron.fore)               #予測値のプロット
lines(LakeHuron, lty=2)        #実績値

#-----------<SARIMAモデルとauto.arima>-----------

#図10.9
op <- par(mfrow=c(2, 1), mar=c(4.1, 4.1, 3.1, 1.1))
plot(diff(AirPassengers), main="差分系列")
plot(diff(log(AirPassengers)), main="対数変換")
par(op)

Air.0 <- window(AirPassengers, end=c(1958, 12)) #最後の24観測値を除く
Air.sarima <- auto.arima(log(Air.0)) #モデル選択
summary(Air.sarima)           #要約

#Box and Jenkins
Air.sarima2 <- Arima(log(Air.0), c(0, 1, 1),
        seasonal = list(order = c(0, 1, 1), period = 12))
summary(Air.sarima2)

#図10.10
tsdiag(Air.sarima2)            #モデル診断プロット

Air.fore <- forecast(Air.sarima2, 24) #24期先までの予測

#図10.11
windows(10, 7)
plot(Air.fore)               #予測値のプロット
lines(log(AirPassengers), lty=2)        #実績値

#-----------<外生変数を伴うモデル>-----------

library(mlbench)#パッケージのロード
data(Ozone)     #データのロード
summary(Ozone)  #要約
?Ozone          #変数の確認

#時系列データへの変換と補間

library(zoo)                      #パッケージのロード
start <- "1976-01-01"             #開始日
x.Date <- as.Date(start) + 0:365  #１年間の日付文字列 
オゾン <- na.spline(zoo(Ozone[ ,4], x.Date))
気圧高度 <- na.spline(zoo(Ozone[ ,5], x.Date))
風速 <- na.spline(zoo(Ozone[ ,6], x.Date))
湿度 <- na.spline(zoo(Ozone[ ,7], x.Date))
気温 <- na.spline(zoo(Ozone[ ,8], x.Date))

#データの概観

#図10.12
op <- par(mfrow=c(2, 2), mar=c(4.1, 4.1, 3.1, 1.1))
plot(気圧高度, オゾン, pch=16, col="brown", lwd=2)
plot(風速, オゾン, pch=16, col="darkgreen", lwd=2)
plot(湿度, オゾン, pch=16, col="steelblue", lwd=2)
plot(気温, オゾン, pch=16, col="tomato", lty=2, lwd=2)
par(op)

#auto.arimaによるモデル選択

X <- data.matrix(cbind(気圧高度, 風速, 風速^2, 湿度, 気温))
n <- 366             #データ数
idx <- 1:(n-31)      #最後の31日分の実績を除くためのインデックス

ozone.arimax <- auto.arima(オゾン[idx], xreg=X[idx, ])
ozone.arimax

#t値を計算するための関数

ts.t_value <- function(fit) {
    coef <- fit$coef
    se <- sqrt(diag(fit$var.coef))
    t_value <- coef/se
    t_value
}

#t値
ts.t_value(ozone.arimax)

#風速を除く

ozone.arimax2 <- auto.arima(オゾン[idx], xreg=X[idx, -(2:3)])
ozone.arimax2

#図10.13
tsdiag(ozone.arimax)            #モデル診断プロット

#対数変換

ozone.arimax3 <- auto.arima(log(オゾン[idx]), xreg=X[idx, -(2:3)])
ozone.arimax3
ts.t_value(ozone.arimax3)

#図10.14
tsdiag(ozone.arimax3)            #モデル診断プロット

#予測

ozone.fore <- forecast(ozone.arimax3, 31, 
        xreg=X[-idx, -(2:3)])      #最後の31日分の予測
#図10.15
windows(10, 7)
plot(ozone.fore)                   #予測値のプロット
lines(log(オゾン[-idx]), lty=2)    #実績値

#信頼区間

実績値 <- オゾン[-idx]
点推定 <- exp(ozone.fore$mean)       #点推定
下限 <- exp(ozone.fore$lower)        #下限
上限 <- exp(ozone.fore$upper)[ ,2:1] #上限

#予測一覧

tbl <- cbind(上限, cbind(実績値, 点推定), 下限)
dimnames(tbl)[[2]] <- c("95％上限", "80％上限", "実績値",
        "予測値", "80％下限", "95％下限")
round(tbl, 1)

#======================================================
#	10.2 GARCHモデル
#======================================================
#パッケージのロード
library(quantmod)    #ファイナンス時系列を処理するため
library(ghyp)        #一般化双曲型分布のため
library(rugarch)     #GARCHモデルのため

#-----------<quantmod>-----------

#インターネット環境が使用できる場合

AAPL <- getSymbols("AAPL")   #YahooからApple株の日次価格をダウンロードする。
UNRATE <- getSymbols("UNRATE", src="FRED") #FREDから失業率をダウンロードする。

#インターネット環境が使用できる別の環境における、データ保存

saveSymbols("AAPL", file.path = "~/r/適切なフォルダへのパス")

#仕方がないので、こちらを使用する。

load("C:/r/users/book/data/AAPL.RData")  #データのロード
head(AAPL)

#プロット
#図10.16
windows(10, 7)
chartSeries(AAPL, theme = "white")               #この図の表示は割愛
chartSeries(AAPL$AAPL.Adjusted, theme = "white") #補正系列

#対数収益

AAPL.rtn <- diff(log(AAPL$AAPL.Adjusted)) *100 #対数収益を計算
AAPL.rtn <- as.zoo(AAPL.rtn)[-1]

#対数収益のプロット

#推移
#図10.17
plot(AAPL.rtn, col="chocolate", xlab="年", ylab="対数収益率")
abline(h=0)

#分布
#図10.18
dens <- density(AAPL.rtn)
plot(dens, main="Apple株の日次対数収益の分布")
rug(AAPL.rtn, side=1, col="springgreen")
abline(v=0, lty=2, col="blue")
lines(dens$x, dnorm(dens$x, mean=mean(AAPL.rtn), sd=sd(AAPL.rtn)), 
	lty=3, lwd=3, col="red")
legend("topleft", legend=c("原系列", "正規分布による近似"), 
	 bty="n",
        lty=c(1, 3), lwd=c(1, 3), col=c("black", "red"))

#一般化双曲型分布の当てはめ

fit.ghd <- fit.ghypuv(AAPL.rtn, symmetric = FALSE,
                     control = list(maxit = 1000), silent=T)
ghd.dens <- dghyp(dens$x, fit.ghd)

#図10.19
plot(dens, main="Apple株収益への一般化双曲型分布の当てはめ")
rug(AAPL.rtn, side=1, col="springgreen")
lines(dens$x, ghd.dens, lty=3, col="darkgreen",lwd=3)
legend("topleft", legend=c("原系列", "一般化双曲型分布による近似"),
       bty="n", 
	lty=1:3, col=c("black", "darkgreen"), lwd=c(1,3))

#-----------<ARCHモデル>-----------
#-----------<ACFプロットとPACFプロット>-----------
#図10.20
op <- par(mfrow=c(2, 1), mar=c(4.1, 4.1, 3.1, 1.1))
acf(AAPL.rtn, lag.max=30, na.action=na.pass)
pacf(AAPL.rtn, lag.max=30, na.action=na.pass)
par(op)

#-----------<ARIMA>-----------

AAPL.arima <- Arima(AAPL.rtn, order=c(3, 0, 2), 
        optim.control=list(maxit=500))
AAPL.arima
ts.t_value(AAPL.arima)

#図10.21
tsdiag(AAPL.arima) #診断プロット

#-----------<標準化残差の２乗>-----------

res <- resid(AAPL.arima, standerdize=T) #標準化残差

#図10.22
op <- par(mfrow=c(2, 1), mar=c(4.1, 4.1, 3.1, 1.1))
acf(res^2, lag.max=30, na.action=na.pass)
pacf(res^2, lag.max=30, na.action=na.pass)
par(op)

#-----------<ARCH(10)>-----------

spec <- ugarchspec(variance.model=list(model="sGARCH", 
    garchOrder=c(10, 0)),
    mean.model=list(armaOrder=c(2, 2)), 
    distribution.model="norm")
AAPL.arch10 <- ugarchfit(AAPL.rtn, spec=spec)
AAPL.arch10

#-----------<GARCHモデル>-----------
#-----------<GARCH(1,1)-norm>-----------

spec.norm <- ugarchspec(variance.model=list(model="sGARCH", 
    garchOrder=c(1, 1)),
    mean.model=list(armaOrder=c(2, 2)), 
    distribution.model="norm")
AAPL.garch.norm <- ugarchfit(AAPL.rtn, spec=spec.norm)
AAPL.garch.norm

#プロット
#図10.23
op <- par(mfrow=c(2, 2), mar=c(4.1, 4.1, 3.1, 1.1))
plot(AAPL.garch.norm, which=10)
plot(AAPL.garch.norm, which=11)
plot(AAPL.garch.norm, which=8)
plot(AAPL.garch.norm, which=9)
par(op)

#-----------<GARCH(1,1)-ghyp>-----------

spec.ghyp <- ugarchspec(variance.model=list(model="sGARCH", 
    garchOrder=c(1, 1)),
    mean.model=list(armaOrder=c(2, 2)), 
    distribution.model="ghyp")
AAPL.garch.ghyp <- ugarchfit(AAPL.rtn, spec=spec.ghyp)
AAPL.garch.ghyp

#プロット
#図10.24
op <- par(mfrow=c(2, 2), mar=c(4.1, 4.1, 3.1, 1.1))
plot(AAPL.garch.ghyp, which=10)
plot(AAPL.garch.ghyp, which=11)
plot(AAPL.garch.ghyp, which=8)
plot(AAPL.garch.ghyp, which=9)
par(op)

#-----------<ローリング予測>-----------

AAPL.roll <- ugarchroll(spec.ghyp, AAPL.rtn, refit.window="moving")
AAPL.roll
report(AAPL.roll, type="VaR", VaR.alpha=0.05, conf.level=0.95)

#プロット
#図10.25
op <- par(mfrow=c(2, 1), mar=c(4.1, 4.1, 3.1, 1.1))
plot(AAPL.roll, which=4, VaR.alpha=0.05) #５％水準
plot(AAPL.roll, which=4, VaR.alpha=0.01) #１％水準
par(op)